Abdessamad Nait-Ouhra

Alexander Farutin

Hamid Ez-Zahraouy

Abdelilah Benyoussef

Chaouqi Misbah

Rheology of a confined vesicle suspension

Rheology of a confined suspension of vesicles (model for blood) is analyzed numerically in two dimensions as a function of the viscosity contrast λ (defined as the ratio between the inner and the outer viscosities ηin and ηout). We consider both the dilute, semi dilute and more concentrated regimes up to about 35% of vesicle area fraction. This paper is a follow up of our recent work [A. Nait-Ouhra et al., Phys. Rev. Fluids 3, 123601 (2018).] where we have found that depending on λ the vesicle can either be centered or off-centered in a channel. Here we will analyze in particular how the existence of the new attractor (off-centered solution) affects rheology. In the dilute regime, it is found that the effective viscosity corresponding to the off-centered solution may have a lower value than that of the centered solution. Moreover the effective viscosity shows a discontinuity as a function of λ. This discontinuity is traced back to the existence of a saddle node bifurcation of the lateral position. When the flow strength is increased the viscosity behavior becomes continuous, owing to the evolution of the transition of the off-centered solution from a saddle-node to a pitchfork bifurcation. We analyze then rheology as a function of the suspension concentration ϕ for different values of λ. For large enough λ the viscosity may be lower or larger than that corresponding to small λ, depending on initial configuration. There is a critical concentration beyond which the initial configuration is irrelevant, but still the behavior of the suspension viscosity with ϕ for low and large λ are quite different both qualitatively and quantitatively. This is traced back to special spatial organization of the confined suspensions, which depends on λ.

I. INTRODUCTION

Rheology of confined suspensions, where the suspended entity size becomes of the order of the available gap, has received so far much less attention than their unconfined counterparts. There are many examples where the confinement effects come to the fore. For example, in microfluidics and nanofluidics [START_REF] Kirby | Micro-and nanoscale fluid mechanics: transport in microfluidic devices[END_REF][START_REF] Andrew | Modification of flow through silica micro-capillaries via polymer brushes[END_REF][START_REF] Wang | Covalent modified hydrophilic polymer brushes onto poly (dimethylsiloxane) microchannel surface for electrophoresis separation of amino acids[END_REF][START_REF] Constable | Characterization of polymer brushes in capillaries[END_REF][START_REF] El-Ali | Cells on chips[END_REF][START_REF] Huh | Microfluidics for flow cytometric analysis of cells and particles[END_REF][START_REF] Stein | Pressure-driven transport of confined dna polymers in fluidic channels[END_REF][START_REF] Howard A Stone | Engineering flows in small devices: microfluidics toward a lab-on-a-chip[END_REF][START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Jäggi | Microfluidic depletion of red blood cells from whole blood in high-aspect-ratio microchannels[END_REF][START_REF] Michael | Fluid dynamics of dissolved polymer molecules in confined geometries[END_REF] the transported entities (cells, macromolecules, such as DNA, and so on) have often sizes which can be comparable to the typical size of the transport devices. We therefore expect the interaction of the suspended entities with walls to play decisive roles. For example, walls may create depletion (zones), due either to entropic or hydrodynamic effects [START_REF] Jäggi | Microfluidic depletion of red blood cells from whole blood in high-aspect-ratio microchannels[END_REF][START_REF] Starkey | The laminar flow of streams of suspended particles[END_REF][START_REF] Segre | Radial particle displacements in poiseuille flow of suspensions[END_REF][START_REF] Howard | A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure[END_REF][START_REF] Ma | Theory of shear-induced migration in dilute polymer solutions near solid boundaries[END_REF][START_REF] Secomb | Blood flow in the microcirculation[END_REF][START_REF] Pranay | Depletion layer formation in suspensions of elastic capsules in newtonian and viscoelastic fluids[END_REF] that should affect transport properties. A prominent example is blood flow where red blood cells (RBCs) interaction with the vessel walls leads to a cell-free layer close to the wall. This cell-free layer plays a pivotal role in oxygen transport, since it leads to a drastic reduction of the effective viscosity with decreasing vessel diameter, allowing thus for a more efficient transport in small vessels of the vascular network where oxygen (and other elements, such as ATP [START_REF] Zhang | Atp release by red blood cells under flow: Model and simulations[END_REF]) is delivered to tissues and organs. The decrease of blood viscosity with vessel diameter is known as the Fåhraeus-Lindqvist effect [START_REF] Fåhraeus | The viscosity of the blood in narrow capillary tubes[END_REF]. A basic understanding of the behavior of confined suspensions requires, on the one hand, the identification of the key parameters that play a relevant role, and on the other hand linking the microscale (e.g. the spatiotemporal organization) to the macroscale. The present study is directed along this objective. We will see that the viscosity contrast λ (defined as the ratio between the inner viscosity η in and the outer one η out ) plays an essential role in the spatial organization, which in turn affects rheology. The rheology of confined suspensions has been addressed quite sparsely to date. Regarding suspensions of rigid spheres, Peyla et al. [START_REF] Peyla | New confinement effects on the viscosity of suspensions[END_REF][START_REF] Davit | Intriguing viscosity effects in confined suspensions: A numerical study[END_REF][START_REF] Ashok S Sangani | Roles of particle-wall and particle-particle interactions in highly confined suspensions of spherical particles being sheared at low reynolds numbers[END_REF] have shown that for a certain range of ratio between the gap size and the particle size, the interactions between particles leads to a decrease of the total dissipation. 2D and 3D numerical simulations of confined RBCs [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF][START_REF] Shen | Blood crystal: emergent order of red blood cells under wall-confined shear flow[END_REF] have been performed using vesicles and capsules as a model. It has been shown that confinement leads to a subtle ordering with a strong impact on rheology. Similar ordering [START_REF] Fornari | Rheology of confined non-brownian suspensions[END_REF] has been reported on confined suspensions of rigid paricles in the presence of inertia.

Our study is motivated by blood flow. However, since our objective is to analyze the basic fundamental features, we focus here on a simplified model: a 2D suspension of vesicles (a vesicle is a closed phospholipid membrane [START_REF] Lipowsky | Structure and dynamics of membranes: I. from cells to vesicles/II. generic and specific interactions[END_REF]). It has been shown that 2D vesicles capture many features known for 3D models of RBCs. Indeed, both models show several common shapes and dynamics, namely: tank-treading and tumbling [START_REF] Thierry Biben | Three-dimensional vesicles under shear flow: Numerical study of dynamics and phase diagram[END_REF][START_REF] Ghigliotti | Rheology of a dilute two-dimensional suspension of vesicles[END_REF], slipper and parachute shapes [START_REF] Farutin | Analytical and numerical study of three main migration laws for vesicles under flow[END_REF][START_REF] Kaoui | Why do red blood cells have asymmetric shapes even in a symmetric flow?[END_REF], flow alignment [START_REF] Thierry Biben | Three-dimensional vesicles under shear flow: Numerical study of dynamics and phase diagram[END_REF][START_REF] Nait-Ouhra | Lateral vesicle migration in a bounded shear flow: Viscosity contrast leads to off-centered solutions[END_REF] and very recently polylobes shapes [START_REF] Mauer | Flow-induced transitions of red blood cell shapes under shear[END_REF][START_REF] Zhang | Simulation of Blood Microcirculation and Its Coupling to Biochemical Signaling[END_REF] etc. In addition, both models show the same rheological behavior [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF]. These observations provide confidence that the adopted 2D system may retain the essential effects of a more realistic model, and at the same time offers low enough numerical cost for a more systematic investigation.

An important factor in our study is the role of the viscosity contrast λ. It is traditionally reported that under shear flow (in 2D) two basic dynamics prevail: tank-treading (TT, obtained for low λ) and tumbilng (TB, obtained for large λ). In the presence of confining walls, a vesicle under shear flow (obtained by two countertranslating plates) is known to migrate towards the centerline, be it in the TT or TB regime (in the TB regime the migration velocity is smaller than in the TT one) [START_REF] Nait-Ouhra | Lateral vesicle migration in a bounded shear flow: Viscosity contrast leads to off-centered solutions[END_REF][START_REF] Coupier | Noninertial lateral migration of vesicles in bounded poiseuille flow[END_REF]. The migration due to a wall has been the subject of many studies [START_REF] Coupier | Noninertial lateral migration of vesicles in bounded poiseuille flow[END_REF][START_REF] Olla | The lift on a tank-treading ellipsoidal cell in a shear flow[END_REF][START_REF] Cantat | Lift force and dynamical unbinding of adhering vesicles under shear flow[END_REF][START_REF] Seifert | Hydrodynamic lift on bound vesicles[END_REF][START_REF] Sukumaran | Influence of shear flow on vesicles near a wall: a numerical study[END_REF][START_REF] Olla | The behavior of closed inextensible membranes in linear and quadratic shear flows[END_REF][START_REF] Olla | The role of tank-treading motions in the transverse migration of a spheroidal vesicle in a shear flow[END_REF][START_REF] Olla | Simplified model for red cell dynamics in small blood vessels[END_REF][START_REF] Kaoui | Lateral migration of a twodimensional vesicle in unbounded poiseuille flow[END_REF].

Recently [START_REF] Nait-Ouhra | Lateral vesicle migration in a bounded shear flow: Viscosity contrast leads to off-centered solutions[END_REF], we have discovered that a vesicle under a confined linear shear flow in the Stokes limit has other final solutions than settling in the centerline. More precisely, for a small enough flow strength and beyond a certain λ the central position coexists with an off-centered one. The off-centered solution arises as a saddle-node bifurcation. When the flow strength increases, the saddle node bifurcation transforms into a pitchfork bifurcation. Here we study how the emergence of the new attractors affects rheology. We shall see that the new solutions significantly impact the rheology, provided that the suspension concentration is small enough. When the concentration exceeds a certain value, we shall see that the coexistence of attractors in the dilute regime is overridden by hydrodynamic interactions. However, the viscosity contrast still plays an important role by changing qualitatively and quantitatively the rheological properties.

II. MODEL AND METHOD

Here we describe briefly the membrane model and the simulation method (for further information see [START_REF] Nait-Ouhra | Lateral vesicle migration in a bounded shear flow: Viscosity contrast leads to off-centered solutions[END_REF][START_REF] Thiébaud | Rheology of a vesicle suspension with finite concentration: A numerical study[END_REF]). We consider a set of 2D vesicles inside a channel, under a planar shear flow (v 0 = γye x ). The vesicle membrane influences the ambient fluid via bending and tension forces, and the total membrane force is given by:

f mem→flu = κ d 2 c ds 2 + 1 2 c 3 n -ζcn + dζ ds t (1)
which includes the local effective tension term ζ in order to fulfill the local arclength conservation constraint (ζ is a Lagrange multiplier). Here c is the local curvature, s is the curvilinear coordinate along the membrane, n and t are the unit normal and tangent vectors, and κ is the membrane bending rigidity. The enclosed area A is conserved automatically thanks to the fluid incompressibility. However, very slight variations of area and perimeter are observed because of numerical errors (see Fig. 1

in [44]

). A derivation of Eq.( 1) can be found in [START_REF] Kaoui | Lateral migration of a twodimensional vesicle in unbounded poiseuille flow[END_REF], obtained from the functional derivative of the Helfrich [START_REF] Zhong-Can | Bending energy of vesicle membranes: General expressions for the first, second, and third variation of the shape energy and applications to spheres and cylinders[END_REF] bending energy E = κ 2 mem c 2 ds + mem ζds. We consider the limit of zero Reynolds number (the Stokes limit). The full problem is described by four dimensionless numbers, which are:

λ = η in η out , C κ = η out γR 3 0 κ ≡ γt c , C n = 2R 0 W , τ = 4πA L 2 (2)
λ is viscosity contrast , where η in , η out denote the inner and the outer viscosities. C κ is the capillary number measuring the flow strength over the bending energy of the membrane. In other words, C κ controls how the shapes of vesicles deform in response to an applied external flow, where t c is the typical time needed for the vesicle to recover its equilibrium shape after cessation of flow. γ is the applied shear rate. C n is the degree of confinement, where W is the channel width. τ is the reduced area, where L is the vesicle perimeter, the effective vesicle radius is defined as R 0 = A/π. Using the Green's function techniques, the velocity v(r 0 ) of a point on the membrane can be written in the following dimensionless expression:

v(r 0 ) = 2 1 + λ v 0 (r 0 ) + 1 2πC κ (1 + λ) mem ds f mem→flu (r).G 2w (r, r 0 ) + 1 -λ 2π(1 + λ) mem ds v(r).T 2w (r, r 0 ).n(r) (3)
where we have used the following scales: R 0 for the length, U = γR 0 for the velocity and κ/R 3 0 for the force per unit area. G 2w (r, r 0 ) and T 2w (r, r 0 ) are the Green's functions (second-and third-order tensors) satisfying the no-slip boundary condition at the bounding walls (see [START_REF] Thiébaud | Rheology of a vesicle suspension with finite concentration: A numerical study[END_REF] for more details). v 0 is the imposed velocity.

The dimensionless force f mem→flu is given formally by the same expression in Eq.( 1) in which membrane bending rigidity κ is set to unity (we keep the same notation for simplicity).

The still unknown Lagrange multiplier ζ is obtained by imposing a divergence-free velocity along the membrane, expressing the membrane incompressibility. The numerical method of imposing membrane incompressibility follows closely that presented in [START_REF] Ghigliotti | Rheology of a dilute two-dimensional suspension of vesicles[END_REF].

The membrane displacement in time is obtained by updating the discretization points after each time iteration, using a Euler scheme, r 0 (t + dt) = v(t)dt + r 0 (t).

III. RESULTS

A. Rheology and its relation to attractors of vesicle positions in the channel

In a recent work [START_REF] Nait-Ouhra | Lateral vesicle migration in a bounded shear flow: Viscosity contrast leads to off-centered solutions[END_REF] we have seen that the vesicle final lateral position in the channel (denoted as h f ) depends on λ. We have found that beyond a critical λ (and for a small enough C κ ) there is emergence of a non centered solution coexisting with the centered one. In this section, we consider the very dilute regime (where the additivity effect due to each vesicle is legitimate) in order to study the impact of the lateral position on the rheology. The effective viscosity is given by:

η ef f = η out (1 + [η]ϕ) (4)
Where [η] is the normalized viscosity (is called the intrinsic viscosity when ϕ → 0), represents the contribution of the vesicles to the viscosity. It is given by

[η] = 1 η out S γ mem ds yf flu→mem,x +η out (λ -1) mem ds (n x v y + n y v x ) ( 5
)
where S is the total area of the vesicles. This expression is the extension of the Batchelor result for suspensions of rigid particles to vesicles [START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF]. For a small enough λ the vesicle always reaches the centerline regardless of initial position [START_REF] Nait-Ouhra | Lateral vesicle migration in a bounded shear flow: Viscosity contrast leads to off-centered solutions[END_REF]. Once the final position is reached the normalized viscosity is extracted from [START_REF] El-Ali | Cells on chips[END_REF]. Figure 1 shows the result. Here parameters are: τ = 0.7, C κ = 1 and C n = 0.4. Below a critical value of λ (λ c 16, a value which depends on the other parameters) the vesicle migrates towards the center of the channel. In this case (centered vesicle), the suspension shows a classical behavior [START_REF] Ghigliotti | Rheology of a dilute two-dimensional suspension of vesicles[END_REF][START_REF] Thiébaud | Rheology of a vesicle suspension with finite concentration: A numerical study[END_REF][START_REF] Danker | Rheology of a dilute suspension of vesicles[END_REF][START_REF] Zhao | The dynamics of a non-dilute vesicle suspension in a simple shear flow[END_REF][START_REF] Rahimian | Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method[END_REF][START_REF] Vitkova | Micro-macro link in rheology of erythrocyte and vesicle suspensions[END_REF][START_REF] Danker | Dynamics and rheology of a dilute suspension of vesicles: Higher-order theory[END_REF] : in the tank-treading motion (TT) [η] decreases with λ until a minimal value. Then [η] increases in the tumbling regime (TB) (Fig. 1a). Above λ c , if the initial position is close enough to the wall (the distance from the lower wall in the simulation is 0.72R 0 ), the vesicle settles in an off-center position [START_REF] Nait-Ouhra | Lateral vesicle migration in a bounded shear flow: Viscosity contrast leads to off-centered solutions[END_REF] as shown in Fig. 1b. This solution coexists with the centered solution. The vesicle aligns with the flow and we have referred to this solution as flow alignment (FA), whereas in the coexisting centered regime the vesicle shows TB. The two corresponding branches are labeled FA and TB in Fig. 1a. In other words, the suspension presents two possible values of viscosity due to the existence of two different attractors (see Fig. 2 in [44] for an averaged viscosity over all possible final configurations). Whether the system selects one solution or the other depends on initial conditions. The viscosity in the FA regime is significantly smaller than that in the TB regime. This is quite intuitive since a TB vesicle has a larger cross section in the channel than a FA vesicle, opposing thus higher resistance to the flow. When λ increases further the viscosity of the FA solution shows an increase of the suspension viscosity. This is due to the fact that the FA vesicle is closer and closer to the wall, and thus the cell free layer becomes smaller and smaller. This is reminiscent of the Fåhraeus-Lindqvist effect [START_REF] Fåhraeus | The viscosity of the blood in narrow capillary tubes[END_REF].

B. Effect of degree of confinement

In this section we investigate the effect of the degree of confinement on the dynamics and rheology of a very dilute suspension of vesicles. We take a single vesicle initially placed close to the lower wall. We set τ = 0.7, C κ = 1 and vary λ (key parameter in this study) for several C n . The results are depicted in Fig. 2 which shows the same qualitative behavior as seen in the previous section. The same dynamics of the vesicle for C n = 0.4 (previous section) persists for all values of confinement investigated here, that is, TT, TB and FA motion are observed. For each value of C n (Fig. 2) there is a different critical viscosity contrast λ c , below which the vesicle migrates to the centerline and [η] shows a non-monotonic behavior as already discussed in the previous section. Beyond λ c the vesicle settles at an off-center position, exhibiting a sudden decrease of the normalized viscosity [η] of the suspension (Fig. 2), which then increases with increasing λ. In conclusion, the same trend is observed pointing to the robustness of the phenomenon.

C. Effect of capillary number

In this section we consider the effect of capillary number C κ on the dynamics and rheology of the vesicle by fixing the other parameters, C n = 0.4, τ = 0.7. Figure 3 illustrates the effect of C κ on the dynamics and rheology. The main difference with previous results is that upon increasing C κ the gray line representing the off-centered position as a function of λ exhibits weaker and weaker jumps at the bifurcation point, until the transition becomes continuous. The saddle-node bifurcation transforms into a pitchfork bifurcation. This transition also affects the behavior of the normalized viscosity [η] which shows weaker jumps until showing a continuous behavior. A point worth of mention is that the TB regime completely disappears in favor of a continuous transition from a TT solution to a FA one. We also note that for a given λ in the off-centered case (FA regime in Fig. 3), the distance between the vesicle and the wall increases as a function of capillary number C κ implying a decrease of the normalized viscosity with C κ (shear-thinning behavior).

D. Effect of reduced area

We carried out a systematic study on the effect of reduced area τ (for a circle τ = 1, otherwise τ < 1) on the equilbrium lateral position of a vesicle and its effect on the rheology of the suspension. The results are reported in Fig. 4. The other parameters are fixed to C n = 0.4 and C κ = 1. When τ is small enough, we find the same trends as before. However, as τ increases (approaching 1), the situation changes significantly. For example, for τ = 0.9 while the vesicle migrates out of center (for λ > 40), the normalized viscosity still exhibits a jump. Contrary to the other values of τ, the viscosity is higher for the off-centered solution at the transition point. This is due to two effects: (i) the higher value of τ lowers the lift due to wall so that the fluid left between the vesicle and the wall is smaller, leading to higher dissipation. (ii) The vesicle has a shape close to a circle causing a higher cross-section against the flow. These two effects together yield higher viscosity than a tumbling vesicle at the center. Even if the vesicle tumbles at the center, due to its quasi-circular shape, its elongation is weaker now causing less and less resistance against the flow.

IV. SUSPENSIONS

We have seen above that depending on λ a single vesicle may migrate to the center or towards the wall and this significantly alters the rheology of the suspension. The question naturally arises of whether this off-centered position may survive for higher concentrations. We have investigated this question by varying concentration and by choosing various initial vesicle positions with the aim to see whether or not the final configuration is sensitive to initial conditions. Periodic boundaries are used along x direction, the box length is L x = 40.5R 0 and the enclosed area of a single vesicle is A = πR 2 0 . The other parameters are kept to λ = 20, C n = 0.4 and C κ = 1 for which a single vesicle would select an off-centered position. We first discuss the dilute and the semi-dilute regimes. We compare the results to those obtained in Ref. [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF] for another set of parameters (λ = 1, C n = 0.4 and C κ = 1) where only the centered solution exists; i.e. for this viscosity contrast λ = 1, there exists no second attractor (off-centered) for a single vesicle.

The results are presented in Fig. 5. We see that for λ = 20, and below a certain concentration (ϕ 12.4%), the normalized viscosity [η] may have several values for the same concentration (ϕ) depending on the final configuration (which depends on the initial configuration) of the suspension. This means that for dilute and semi-dilute suspensions, several stable spatial configurations are possible which give rise to different values of normalized viscosity [η] or shows in the FA regime a jump with higher values than those of the TB regime for τ = 0.9. This contrasts with the behavior for other value of τ, see text. Note that the TB branch continues for larger λ but here we only show the results until the FA branches appears.

effective viscosity η ef f (Fig. 6). It is interesting to note that, below a certain concentration (ϕ 6.2%), the suspension of relatively rigid vesicles (λ = 20) may have a much lower value of [η] than that of a suspension of soft vesicles (λ = 1 Fig. 5). This lower value for more rigid particles is associated with the existence of stable off-centered positions, as described above. The highly non-monotonic behavior of [η] for λ = 1, described in [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF], as a function of ϕ is due to the existence of ordered solutions with a single (for low enough ϕ) and double (for larger ϕ) file configurations. For a large enough λ for which two attractors exist for a single vesicle, the situation is quite different. We have fixed λ = 20. We find, depending on initial conditions (as described in the caption of Fig. 5) that for a small enough ϕ (< 6.2%) there exist two solutions (i) centered single file (as found for λ = 1) undergoing tumbling giving rise to the viscosity values shown as points (a , b and c ) in Fig. 5 and (ii) off-centered single file solution exhibiting flow alignment having the viscosity shown as points (a, b and c) in Fig. 5. In other words, for low enough ϕ the single ordered file seems to exists whatever the value of λ is. For λ = 20 the existence of FA single file exists only for a low enough concentration (points a, b and c in Fig. 5). In this low concentration regime the normalized viscosity shows a plateau owing to the weak hydrodynamic interaction.

However, when the concentration increases slightly beyond point c in Fig. 5 the hydrodynamic interaction among vesicles starts to play a role, but still the existence of two attractors continues to affect rheology up to about ϕ = 12%. The situation becomes more complex than for a lower concentration. Indeed, the single file FA solution undergoes a transition towards the centered TB file (point d and inset d in Fig. 5). From point c (FA file) to point d (TB centered solution, see inset d) the viscosity increases as expected, since TB vesicles in the centerline have a higher cross section. When the concentration is increased further (from point d to e) the single file TB solution still exists, but [η] decreases. The reason is that the insertion of a new vesicles fills the gap between already existing vesicles and reduces the strength of recirculation zones between vesicles. This reasoning is similar to that evoked for λ = 1 (single file TT solution) [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF]. For solutions d and e the initial condition corresponds to a centered set of solution. It seems thus that the central attractor is robust despite the presence of hydrodynamic interaction. This is, however, not the case for the the second (off-centered) attractor. We have prepared half of vesicles with the off-centered position (at the attractor position for a single vesicle) and the other half at the opposite position (which is also an attractor as well) with the hope of obtaining a double file solution (as with λ = 1 [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF]). All simulations with these initial conditions led to a disordered final solution (represented by inset f in Fig. 5). The same outcome is obtained with random initial position. It seems thus that the second attractor has become unstable in the presence of hydrodynamic interactions.

A further increase of ϕ (for λ = 20) beyond point e, destroys the single file TB solution in favor of an apparently disordered solution (see inset f ), leading to higher and higher viscosity. This is drastically different from the case of λ = 1, where the system organizes into a double file that survives up to about ϕ = 35%. At that concentration, the suspension for λ = 20 shows a significantly higher η ef f , which is about 1.45 times that corresponding to λ = 1.

In summary, the suspension viscosity shows a rich behavior depending on the viscosity contrast λ. For small enough ϕ, the more rigid particles can either show lower or higher normalized viscosity depending on initial conditions. This points to a non trivial fact that data handling in experiments inherently contain statistical variances not only due to noise and natural imperfections and variety of samples, but also due to the very nature of multi-stability of solutions reported on here. This sheds a new light for future interpretations of rheology of soft suspensions. Beyond a certain concentration, the softer particles exhibit order (e.g. double file), whereas more rigid particles show disorder, resulting in a much higher effective viscosity. Note finally, that the appearance of order (soft enough particles) and disorder (for more rigid particles) is not a 2D property, but also a 3D property, shown both theoretically for model systems (drops, capsules, vesicles) and experimentally for red blood cells [START_REF] Shen | Blood crystal: emergent order of red blood cells under wall-confined shear flow[END_REF]. The present study can thus serve as a guide for a systematic analysis in 3D.

V. CONCLUSION

A major central point of this study is that the existence of two attractors leads to a rheology which depends on initial configurations. This happens in the dilute and the semi-dilute regimes, up to about a concentration of 12%. Beyond that concentration, the attractors loose their identities. Still, the viscosity contrast affects significantly the value of the effective viscosity as shown in Fig. 6. Several diseases (such as malaria) are accompanied by hardening of the RBC cytoplasm, a fact which can be mimicked by a higher internal viscosity (and thus higher viscosity contrast). This leads to a higher effective viscosity in microcirculation compromising the efficiency of oxygen transport. This study was dedicated to 2D and to a linear shear flow. An interesting extension should now be made for 3D and for a pipe flow which is a more relevant scenario for blood flow. Ordering of RBCs (and its impact on rheology), similar to the reported result in 2D has already been identified in 3D [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF][START_REF] Shen | Blood crystal: emergent order of red blood cells under wall-confined shear flow[END_REF]. This fact provides some confidence to the present study for its relevance to a more general situations. . This is traced back to the existence of off-centered positions for λ = 20 as discussed in the text. In the more concentrated regime, the suspension of relatively soft vesicles (λ = 1) has a much lower η ef f than that of a suspension of relatively rigid vesicles. This is traced back to the existence of ordered solutions with a double file configuration for λ = 1 and disordered solutions for λ = 20. See the text for more details.

FIG. 1 .

 1 FIG. 1. (Color online)The evolution of the normalized viscosity [η] and equilibrium lateral position h f (a and b resp.) as a function of the bifurcation parameter λ. TT: Tank-treading, TB: Tumbling, FA: Flow alignment. Here τ = 0.7, Cκ = 1 and Cn = 0.4. The dashed gray lines correspond to the results where the vesicle initially started at the centerline or at its vicinity and ultimately gets centered. The centerline is at y/W = 0 on the vertical axis in (Fig.1b) and the lower wall is at y/W = -1/2.

9 FIG. 2 .

 92 FIG. 2. (Color online)The evolution of [η] and equilibrium lateral position h f as a function of λ for several degrees of confinement Cn. Here we set Cκ = 1 and τ = 0.7. The vesicle is initially set at a small enough distance (0.72R0) from the lower wall. TT, TB and FA are shown. The channel center is at zero and the walls are at y/W = ±1/2 on the right vertical axis. Note that the TB branch continues for larger λ but here we only show the results until the FA branches appears.

FIG. 3 .

 3 FIG. 3. (Color online)The evolution of [η] and equilibrium lateral position h f according the bifurcation parameter λ for several capillary numbers Cκ. The vesicle is initially set at a small enough distance (0.72R0) from the lower wall. TT, TB and FA are shown. Beyond a certain value of Cκ (c and d panels) , the TB motion is suppressed and the evolution of [η] as a function of λ does no longer exhibit any jump. The channel center is at zero and the walls are at y/W = ±1/2 on the right y vertical axis. Note that the TB branch (when it exists) continues for larger λ but here we only show the results until the FA branches appears. For the continuous curves, TB regime does not exist.

9 FIG. 4 .

 94 FIG. 4. (Color online)The evolution of [η] and equilibrium lateral position h f as a function of the bifurcation parameter λ for several reduced areas τ. Here we set Cκ = 1 and Cn = 0.4. The vesicle is initially set at a small enough distance (0.72R0) from the lower wall. TT, TB and FA are shown. The channel center is at zero and the walls are at y/W = ±1/2 on the right y vertical axis. The normalized viscosity[η] shows in the FA regime a jump with higher values than those of the TB regime for τ = 0.9. This contrasts with the behavior for other value of τ, see text. Note that the TB branch continues for larger λ but here we only show the results until the FA branches appears.

FIG. 5 .

 5 FIG. 5. (Color online) The evolution of [η] as a function of concentration ϕ for typical values of λ. Here we set Cκ = 1 and Cn = 0.4. The initial position of the vesicles has impact on the rheology up to ϕ 12.4%. Points a, b and c on the curve correspond to FA configuration (a typical configuration is shown in inset c). Points a , b and c on the curve correspond to another type of configuration (centered TB; a typical configuration is shown in inset c ) obtained from different initial conditions as compared to a, b and c configurations. Point b on the curve corresponds to another configuration which is a mixed state, in which a vesicle undergoes an off-centered TB, whereas the second one shows FA (see inset b ). Points d and e on the curve correspond to centered TB configuration (see inset d and e resp.). Points d , e and f on the curve correspond to disordered configuration (a typical configuration is shown in inset f ). The blue curve (λ = 1) is taken from Ref.[23]. Configurations a, b and c are obtained from initial conditions where the vesicle are equidistant and set close to lower wall (0.72R0 distance from the wall). Configuration d and e are obtained with initial set of equidistant vesicles at the center. Configuration f is obtained with random initial positions. Configurations a , b and c are obtained from initial conditions where the vesicle are equidistant and set at the center. It is noteworthy that we may also obtain the points b and c with initial random positions with a much longer time span. Configurations d and e are obtained with initial random positions. Configuration b is obtained with a vesicle at the center and the other close to the wall (0.72R0 distance from the wall).

FIG. 6 .

 6 FIG. 6. (Color online) The evolution of the effective viscosity η ef f as a function of concentration ϕ for two values of λ. The filled circles (red curve) and the stars correspond to λ = 20, a value which leads to different spatial configurations of the suspension which in turn give rise to different values of η ef f (or [η]) as explained in the caption of Fig.5. The stars are the counterparts of the points a -b and b in Fig.5. In the dilute regime the suspension of relatively rigid vesicles (λ = 20) may have a lower value of η ef f than that of a suspension of soft vesicles (λ = 1). This is traced back to the existence of off-centered positions for λ = 20 as discussed in the text. In the more concentrated regime, the suspension of relatively soft vesicles (λ = 1) has a much lower η ef f than that of a suspension of relatively rigid vesicles. This is traced back to the existence of ordered solutions with a double file configuration for λ = 1 and disordered solutions for λ = 20. See the text for more details.

VI. ACKNOWLEDGEMENTS

We thank CNES (Centre National d'Etudes Spatiales), and the French-German university programme "Living Fluids" (Grant CFDA-Q1-14) for a financial support. A. N-O thanks CNRST Morocco for awarding a scholarship of excellence (Grant Ref. k 1/045). The simulations were performed on the Cactus cluster of the CIMENT infrastructure, which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA).


Wei-Fan Hu

Te-Sheng Lin

Salima Rafai

Chaouqi Misbah

Chaotic swimming of phoretic self-propelled particles

Swimming of a rigid phoretic particle in an isotropic fluid is studied numerically as a function of the dimensionless solute emission rate (or Péclet number Pe). The particle sets into motion at a critical Pe. Whereas the particle trajectory is straight at small enough Pe, it is found that it looses its stability at a critical Pe in favor of a meandering motion. When Pe is increased further the particle meanders at short scale but its trajectory wraps into a circle at larger scale. Increasing even further Pe causes the swimmer to escape momentarily the circular trajectory in favor of chaotic motion lasting for a certain time, before regaining a circular trajectory, and so on. The chaotic bursts become more and more frequent as Pe increases, until the trajectory becomes fully chaotic, via intermittency scenario. The statistics of the trajectory is found to be of run and tumble-like nature at short enough time, and of diffusive nature at long time without any source of noise.

Introduction.-The field of microswimmers covers today an extremely wide panel of systems going from prokaryotic and eukaryotic microorganisms (such as bacteria, algae, leukocytes, and so on) to artificial microswimmers [1,2]. A prototypical example of the later is Janus-like particle, named after the two-faced Roman god, their motion originates from the asymmetry of their surface properties. Other popular designs of artificial swimmers consists of active droplets [3][4][5][6][7], involving a chemical reaction, or a constant local source/sink of the solute. In these circumstances a spontaneous symmetry breaking of the solute distribution field on the surface (or in the bulk) may emerge in the form of a solute polarity giving rise to hydrodynamical flow (of Marangoni type) leading to self-propulsion.

Several studies have been devoted to the elucidation of the mechanisms of such Marangoni driven propulsion [8,9]. In particular, the onset of self propulsion has been described in terms of a coupling between a chemical reaction that produces a surface tension gradient at the surface of the droplet and a mechanism of advection that allows to sustain this gradient against diffusion. The description of the surface activity of the droplet involves in general advection-diffusion-reaction contributions for each species involved in the system (for a recent review, see [10]). Combined theoretical and numerical efforts have allowed elucidating the basic elements at the origin of the occurence of spontaneous locomotion [3,4,6,9]. The outcome is that beyond a critical value of a dimensionless solute emission compared to diffusion (the Péclet number Pe) the droplet sets into a spontaneous unidirectional motion. However, not only straight propulsion but also meandering trajectories have been observed in different experiments [11][12][13][14].

So far, the explanation of complex motions, such as circular and helical trajectories, relies on the complexity of the swimming particle (e.g. chiral shape) [12] or that of the suspending medium, such as viscoelasticty [14]. These various rich behaviors raise naturally the question of whether or not this complexity may already be hidden in purely isotropic media due to the intrinsic nonlinearities of the problem. This is the main objective of this Letter. We show that circular tra-jectories appear in the full model without any restriction on the model symmetry. Moreover, a transition to chaos via intermittency scenario is observed and quantified. The chaotic behavior can lead to a run and tumble-like dynamics of purely deterministic origin.

The model.-We consider a two-dimensional circular shaped particle of radius a immersed in a fluid of dynamic viscosity η. The fluid is assumed to be incompressible and obeys the Stokes equations. Surrounding solute particles are emitted or adsorbed on the particle with isotropic emission rate A . The solute interacts with the particle through a short-range potential with a characteristic length λ that is much smaller than a. In the so-called sharp interface limit this amounts to a slip tangential velocity [6] on the particle surface (see below). The concentration of the solute, denoted by c, diffuse with diffusivity D and is advected by the fluid flow.

As in [6], the length, fluid velocity, the concentration, and the pressure are scaled by the characteristic values a, |A M |/D, a|A |/D, and η|A M |/aD respectively. Here M = ±k B T λ 2 /η is the signed mobility (see [15]) defined with the Boltzmann constant k B and the temperature T . The model equations in dimensionless form are described in the following: In a co-moving frame attached to the particle center the tangential slip velocity on the particle surface is given, in a polar coordinate system, by u(1, θ ) = M∇ s c, where u(r, θ ) is the fluid velocity, ∇ s = (∂ /∂ θ )τ (where τ is the unit tangent vector) is the surface gradient operator and M = M /|M | = ±1 is the dimensionless mobility. In the far field (r → ∞), the velocity field converges to the translational phoretic velocity -U. The solute concentration c(r, θ ,t) is governed by the following advection-diffusion equation and the Stokes equation

∂ c ∂t + u • ∇c = 1 Pe ∆c, ∆u -∇p = 0, ∇ • u = 0, (1)
where the Péclet number is defined by

Pe = |A M |a/D 2 .
The associated boundary condition of surface activity reads

∂ c ∂ r (1, θ ,t) = -A, where A = A /|A | = ±1
is the dimensionless emission rate (A > 0: emission, A < 0: adsorption). One should note that, in most previous studies that considered the three-dimensional cases, the fluid domain is unbounded and the concentration is assumed to attenuate in the far-field limit (c → 0 as r → ∞). However, such an assumption is not appropriate in 2D. Indeed, it is known that, in 2D the solute transport equation (1) for Pe = 0 with the far-field attenuation condition does not support steady-state solutions due to the logarithm divergence [16]. In order to cure this problem, we assume that the size of the fluid domain is a finite number R and the attenuation condition is thus replaced by the outer Dirichlet boundary condition c(R, θ ,t) = 0. A preliminary 3D study recovers the occurence of complex motion, presented below.

The force-free constraint on the solid particle (using the reciprocal theorem) to relate the phoretic velocity U(t) to the slip velocity (and thus to concentration from the above slip condition) [17]:

U(t) = - 1 2π 2π 0 M∇ s c(1, θ ,t) dθ . (2)
Since the concentration c is periodic in θ -direction, it can be expressed by Fourier series as

c(r, θ ,t) = ∞ ∑ k=-∞ ĉk (r,t)e ikθ ,
where ĉk are the Fourier modes. Substituting this representation into Eq. ( 2), one finds (in Cartesian components):

U(t) = M [-Re( ĉ1 (1,t)), Im( ĉ1 (1,t))] , (3)
where Re and Im denote real and imaginary parts, respectively. The velocity field can be expressed as u = (

1 r ∂ ψ ∂ θ , -∂ ψ ∂ r ) (in polar components)
, where ψ is the stream function and has the following analytical form [18,19]

ψ(r, θ ,t) = ∞ ∑ k=-∞ 1 -r 2 2r |k| ikM ĉk (1,t)e ikθ . ( 4
)
where we have used the boundary condition u(1, θ ) = M∇ s c to express the series coefficients in terms of c. Thus, the velocity u (as well as are U) is given in terms of c, and when injected into (1) yields a closed nonlinear equation for c. Solution ( 4) is valid for infinite system, but for a finite size R = 200 the residual errors are of order 0.005 (see [START_REF](i) the scaling of the radius of the circular trajectory as a function of Pe and (ii) the evolution of the spectrum of different harmonics as a function of Pe[END_REF]). In the following we will see how the full model, not restricted to symmetric solutions, leads to the emergence of rich dynamics, turning into chaos at large Péclet numbers.

Stationary solution and linear stability analysis.-A stationary solution where there is no net flow and zero phoretic velocity U(t) ≡ 0 exists at all Péclet numbers with the solute concentration c 0 (r) = A ln(R/r). We perform a linear stability analysis of this solution by introducing an infinitesimal perturbation φ (r, θ ,t) = ∑ ∞ k=-∞ φk (r,t)e ikθ with φ c 0 . Neglecting higher order terms, the following relation for the Fourier mode φk is obtained:

∂ φk ∂t = -AMk 2 1 -r 2 2r |k|+2 φk (1,t) + 1 Pe ∂ 2 ∂ r 2 + 1 r ∂ ∂ r - k 2 r 2 φk . (5)
Since r > 1, the coefficient of the first term on the right hand side of the equation is positive/negative when A and M have the same/opposite sign, respectively. That is, the first term may be stabilizing or destabilizing, while the second term represents diffusion that is always stabilizing. Whether or not a small perturbation of the concentration is able to trigger the spontaneous symmetry-breaking swimming motion depends on the competition between these two mechanisms. In the case when A and M have opposite signs any fluctuation acting on the stationary solution should die out and hence for any Péclet number there is only one observable solution that is stationary. Henceforth, we will concentrate on the nontrivial case where A and M have the same sign. For definitness we set

A = M = 1.
We first perform the linear stability of the stationary solution. The eigenvalue problem reads σ k φk = L k φk , where σ k is the eigenvalue (we look for perturbations in the form e σ k t ) and L k is the linear operator on the right hand side of Eq. ( 5). The eigenvalue problem is solved numerically using a Chebyshevspectral method. In all simulations, we set the size of the outer domain R = 200. It is found that, for k = 1, one of the eigenvalues becomes positive at Pe ≈ 0.466 so that the stationary solution loses its stability at this Pe value. This marks the transition between a nonmotile to motile particle. The mode k = 1 corresponds to polarity of the concentration field.

Spontaneous symmetry-breaking autophoretic motion, Pe > Pe 1 .-Once the instability threshold is reached (Pe 1 ∼ 2/ ln(R); see the full expression in [START_REF](i) the scaling of the radius of the circular trajectory as a function of Pe and (ii) the evolution of the spectrum of different harmonics as a function of Pe[END_REF]) nonlinear terms are necessary in order to fix the amplitude of the swimming speed. This task is handled numerically. We employ the second-order Runge-Kutta time-marching scheme while spatial derivatives in r-direction are performed either by second-order finite difference method or Chebyshev spectral differentiation, derivatives in θ -direction are achieved through Fourier spectral differentiation. In the present time-stepping method, the nonlinear advection term is treated in an explicit manner while the diffusion term is solved implicitly (see [START_REF] Lai | [END_REF] for more detail about the fast solver of the present scheme).

A symmetry-breaking solution occurs when Pe ≥ Pe 1 0.466. We find that in the range 0.466 Pe 4.65 the particle swims in a given direction (fixed by initial conditions) with constant velocity. The concentration distribution at Pe = 3, exhibiting a comet-like pattern. In this case, there exists a nonzero concentration gradient on the particle surface and the particle propels itself in a straight direction with the constant phoretic velocity U ≈ 0.148.

Figure 1 shows the phoretic velocity as a function of the Péclet number. The particle sets into a directed motion in the form of a supecritical bifurcation. Note that a non-monotonic variation of the phoretic velocity is observed in the range 0.466 Pe 4.65 and the particle attains the highest swimming speed at the optimal value of Pe ≈ 2.1. This result is in qualitative agreement with that found in [6] for threedimensional axisymmetric phoretic particles.

Meandering motion, Pe > Pe 2 .-The directed swimming solution looses its stability (a secondary instability) at Pe 2 ≈ P e 4.65 in favor of a meandering motion: the particle follows a meandering path as shown in Fig. 3 (a). In this regime the swimming velocity is no longer a constant but periodically oscillates in time. As an example, we show in Fig. 2 the time-periodic concentration distribution at Pe = 6. These solutions are shown as the red solid line in Fig. 1. For such time-periodic solution the swimming speed is defined by

U = T -1 T 0 U(t)
dt, which is measured over one time period, T . This transverse instability is the first initial stage for a more complex dynamics, as discussed below. The meandering period scales as aD/(AM), which is a compromise between diffusion, flow mobility along the surface and the emission rate. An interesting feature has emerged: the meandering swimmer exhibits a higher velocity than the directed one. Although the directed motion is unstable, we can still follow this branch by solving the steady-state problem (a steady directed motion always exist, though it becomes unstable) by imposing to the particle to set at the same rectilinear direction. By this manner, we could determine its velocity. The dashed line shown in Fig. 1 indicates the unstable equilibrium solution corresponding to a directional movement. The meandering solution swims faster than the directed one, presumably due to subtle nonlinear effects.For Pe Pe 2 the period of meandering diverges as (Pe -Pe 2 ) 1/2 , which is symptomatic of a supercritical bifurcation (see [START_REF](i) the scaling of the radius of the circular trajectory as a function of Pe and (ii) the evolution of the spectrum of different harmonics as a function of Pe[END_REF]) . Circular motion, Pe > Pe 3 .-By increasing the Péclet number above Pe 3 12.5145 we find that the above meandering solution becomes unstable in a somewhat remarkable way (in fact, a period doubling bifurcation occurs at Pe 3 ). Figure 2 shows the bead with an asymmetric concentration field leading to trajectory deviation. As seen in Fig. 3(b) for the case of Pe = 12.86291, although the local movement is again zigzaglike (or meander) at short scale, the velocity direction continuously varies with time, and the global dynamics of swimming pattern converges to a circle for the long time behavior (radius is insensitive to system size [START_REF](i) the scaling of the radius of the circular trajectory as a function of Pe and (ii) the evolution of the spectrum of different harmonics as a function of Pe[END_REF]). The particle short scale meandering is maintained, whereas a longer scale dynamics drives the particle trajectory to wrap around a circular trajectory. Note that here the studied system is fully isotropic in a marked contrast with the observed complex trajectories (like meander, circles, and so on) for chiral particles [12]. Here, the present complex motion emerges from the inherent nonlinear advection-diffusion coupling. This coupling leads to an axial symmetry-breaking bifurcation of the comet (Fig. 1), so that in the bead frame the comet drift sideways along the bead surface leading to a circular trajectory. The radius of the circle diverges at Pe 3 (and so does the period of revolution) like (Pe -Pe 3 ) -α , with α 0.2, acting as a macroscopic length scale (a coherence-like length) much larger than the swimmer size. Further development would be needed in order to explain the numerical exponent (α ∼ 0.2; see Fig. S1 in [START_REF](i) the scaling of the radius of the circular trajectory as a function of Pe and (ii) the evolution of the spectrum of different harmonics as a function of Pe[END_REF]).

Chaotic swimming motion, Pe > Pe 4 12.86293.-Increasing the Péclet number even further reveals another instability: the particle quits the circular trajectory and enters into an apparently chaotic regime, where it exhibits few runs and tumbles during some time, before it regains a circular trajectory (Fig. 3(c)). This type of chaotic bursts appears in a non-regular manner, as also shown in Fig. 4 which represents the swimming speed of the particle close to the transition towards chaos. This nonregularity is due to the subcritical nature of the bifurcation towards chaos. Indeed, during a certain time the motion is regular, and the degree of natural intrinsic fluctuations can cause the system to jump from the regular solution to the apparently chaotic one. As the Péclet number is increased further and further the chaotic bursts become more and more frequent, until a critical value Pe 4 13 where the whole trajectory becomes erratic in time. In Fig. 1 we shows the magnitude of the averaged velocity as a function of Pe. This is symptomatic of the intermittency scenario [START_REF] Bergé | Order within chaos: towards a deterministic approach to turbulence[END_REF]. Its particularity is that the limit cycle (in our case the meandering-circular regular solution) looses its stability at a given critical value of Pe in a subcritical fashion (the analogue of a first order transition for systems at equilibrium). That is, the regular solution may coexist with new solutions and the particle can jump back and forth from one solution to the other, when higher modes (k > 1) approaches the instability conditions. This manifests itself in temporary chaotic bursts until coexisting solutions (the analogue of metastability in equilibrium systems) is fully lost and the regime becomes fully chaotic.

In the chaotic regime, the particle exhibits an apparently random motion (of purely deterministic origin), see Fig. 3

(d).

To quantify this specific random walk, we measure the mean square displacement MSD(τ) =< r(t + τ) -r(t) 2 >, where r(t) is the location of the particle at time t and < • > denotes the average along the entire trajectory. Figure 5 reports the mean square displacement (msd) corresponding to the trajectory at Pe = 13. At short times the particle produces a persistent swimming motion and msd is then quadratic in time. At longer times, a de-correlation process due to chaotic turns takes place with msd proportional to t, typical of a classical random walk. Actually, it is not obvious that a chaotic motion is equivalent (at long time) to normal diffusion. There are several chaotic maps yielding anomalous diffusion [START_REF] Geisel | [END_REF]. This is attributed to the topology of the Poincaré map, in which the trajectory may spend long time periods in some regions of phase space. Despite that chaos via intermittency favors long period stay in specific regions of phase space [START_REF] Bergé | Order within chaos: towards a deterministic approach to turbulence[END_REF], a diffusionlike behavior still persists here.

Discussion.-We have identified the emergence of complex dynamics going from straight swimming to chaotic trajectories for a rigid particle in a purely isotropic fluid medium. Chaos is also found in our preliminary 3D study [START_REF](i) the scaling of the radius of the circular trajectory as a function of Pe and (ii) the evolution of the spectrum of different harmonics as a function of Pe[END_REF]). This highlights the fact that the complexity can arise from a minimal version. The chaotic regime can be quantified as runand-tumble-like events. The first unstable mode is k = 1, that provides directed swimming. As Pe increases higher harmonics becomes unstable, first k = 2 and k = 3, and so on. The activity of these three modes is sufficient to generate chaos. Interestingly, Pe k (k > 1) are independent of system size R.

The time over which the diffusion-like behavior is reached (as well as the duration of the circulation trajectory) are much larger than the typical time for the growth of the linear modes ∼ 1/σ 1 . This feature may be attributed to the very nature of chaos via intermittency. Indeed, in this regime it is known that trajectory in phase space alternates intermittently between slow regular motion close to the marginally stable fixed point (the regular solution), and chaotic bursts. As a consequence, the correlation decay may exhibit long range behavior. From physical point of view, the longtime behavior may be related to the persistence of circular orbits (albeit transients) having a radius much larger than the particle radius Non trivial motions have been observed so far with complex shapes or in complex media [12,14]. The main message of this Letter is the demonstration that a complex motion is possible in an isotopic, nonviscoelastic phoretic system. Experiments on non colloidal systems (autophoretic drops) [7] have provided some evidence of complex trajectories without diging further into their properties. The present results may serve as a guide for future systematic analyses.

15 FIG. 1 .FIG. 2 .

 1512 FIG. 1. The bifurcation diagram of the phoretic velocity U as a function of Pe. The empty circle indicates the bifurcation point Pe 1 = 0.466 and the filled-square indicates the bifurcation point Pe 2 = 4.65. The dashed line corresponds to a steady-state but unstable straight motion and the red line to non straight motions (see below).

FIG. 3 .

 3 FIG. 3. Different trajectories: meander, circular motion with short scale meander, intermittency and chaos.

15 P e = 13 timeFIG. 4 .FIG. 5 .

 151345 FIG.4. The swimming speed U(t) from regular meander to chaos. Chaotic bursts are more and more frequent as Pe increases.

We thank CNES (Centre National d'Etudes Spatiales) (C.M.) and the French-German university program "Living Fluids" (grant CFDA-Q1-14) (C.M. and S.R.). W.-F. Hu acknowledges support by Ministry of Science and Technology, Taiwan under research grant MOST-107-2115-M-005-004-MY2. T.-S. Lin acknowledges support from MOST-107-2115-M-009-008-MY2.


L Khomenkova

K Makasheva

P Petrik

C Bonafos

N Ratel-Ramond

O Melnichuk

L Melnichuk

I Balberg

F Gourbilleau

N Korsunska

Thermally Stimulated Evolution of Optical and Structural Properties of Germanium-Doped Alumina Films

Ge-doped Al 2 O 3 films were deposited on long, non-rotated quarts substrate by magnetron co-sputtering of Ge and Al 2 O 3 space-apart targets. This approach allowed to vary Ge content along the film length in the range 5-85 at.%. Simultaneous annealing treatment of such long samples permitted to investigate the evolution of optical and structural properties of the films with different Ge content. It was observed that all as-deposited films and those annealed at T A =550 °C during 50 min in nitrogen atmosphere were amorphous whatever Ge content. The formation of Ge nanocrystals was observed after treatment at T A =600 °C during 1-15 min. The higher Ge content, the shorter annealing time was required to form Ge nanocrystals. This annealing stimulated also an appearance of photoluminescence in visible-near infrared spectral range. Three specific bands centered at about 580-600 nm, 620-630 nm and 700-730 nm were detected. Their relative contribution was depended on Ge content determining the shape and peak position of luminescence spectrum. Two former components were ascribed to the carrier recombination via the defects in Ge oxide shell covered Ge nanocrystals, whereas the latter one was assumed to originate from carrier recombination in amorphous Ge clusters or via defects at Ge-cluster/Al 2 O 3 interface. This statement was supported by the heat treatment of the films at T A =700-800 °C caused significant decrease of Ge content due to Ge out-diffusion from film volume via formation of volatile GeO.

10.

Introduction

In recent decades, significant attention of research groups has been focused on silicon or germanium nanostructured materials. Since the first report on nanostructured Si-based materials (1,2), optical and structural properties of free standing silicon nanocrystallites (Si-ncs) or those embedded in dielectric hosts were widely studied (3)(4)(5)(6)(7)(8)(9)(10). A lot of efforts were concentrated on the development of different technological approaches for fabrication of Si-ncs, monitoring their size distribution and surface chemistry. An optical gain and electroluminescent emission were achieved offering the perspective to create all-in-one-chip Si-based optoelectronic devices. Significant knowledge on the light emission processes of Si nanostructured materials was obtained. The mechanisms of Sincs formation in different hosts were proposed.

It is interesting that size-dependent optical properties of germanium nanocrystals (Ge-ncs) were reported by Hayashi et al. much earlier than those of Si-ncs (11). Although the Ge-ncs were also studied by different research groups (see, for instance, [12][13][14][15][16][17][18][19][20][21][22], the corresponding publications were not so numerous unlike those on the Si-ncs embedded in different hosts. Recently, the attention of the researchers has been particularly renewed towards the Ge-ncs since they are environmentally friendlier alternatives to classical compound semiconductor nanocrystals being nontoxic, biocompatible, electrochemically stable and compatible with current microelectronics. The investigation of either freestanding Ge-ncs or those embedded in dielectric hosts was reported (23)(24)(25)(26)(27)(28).

As compared to silicon, germanium is especially appealing due to the larger dielectric constant, higher electron and hole mobility, smaller effective masses of electrons and holes. A narrower bulk band gap (0.66 eV for Ge against 1.1 eV for Si) and a larger absorption coefficient offered a possibility for tuning of light emission in wider spectra range. Bohr exciton radius of germanium, r B Ge =24.3 nm, is significantly larger than that of silicon (r B Si =4.5 nm) (12). This means that Ge-ncs will impart stronger, more easily identified quantum confinement effects than Si-ncs with the same sizes. Moreover, quantum confinement effects will emerge for larger Ge-ncs than that for Si-ncs (25). Due to this property, the Ge-ncs are considered to be ideal nodes for memory devices (17,21,23,26). Being crystallized at lower temperature and offering a negative conduction band offset with respect to the Si substrate, the Ge-ncs provide a larger memory window than devices based on Si-ncs with improved retention time (16,21).

The Ge-ncs embedded in SiO 2 host were mainly studied. However, Al 2 O 3 is an interesting material not only for microelectronic application, but also for optical communication due to its relatively higher refractive index (1.73 at 1.95 eV) in comparison with that of SiO 2 (1.46 at 1.95 eV) offering better light confinement. Although Al 2 O 3 waveguides doped with Si or Ge and rare-earth ions were developed for optical communications (28,31), such materials are still less addressed. The present paper reports the fabrication of Ge-doped Al 2 O 3 films with different Ge content by magnetron co-sputtering and the effect of post-deposition processing on the structural and optical properties of these materials.

Experimental details

Ge-rich Al 2 O 3 films were deposited on long non-rotated quarts substrates by radiofrequency magnetron sputtering of two spatially separated targets, pure Ge and pure Al 2 O 3 . The diameter of each target was 2 inches. The length of quartz substrate was 15 cm. The distance between the targets and the substrate was 13 cm. The sputtering was performed in an argon atmosphere at a pressure of 1.1 Pa for 250 minutes. The power densities applied to Al 2 O 3 and Ge targets were 0.98 and 0.49 W/cm 2 , respectively. The substrate temperature was 15 °C. Using this deposition approach, the films with variation of Ge content along the film length were obtained. Besides, pure Ge and pure Al 2 O 3 films were also produced by the sputtering of corresponding target only at the same conditions. More details about this approach can be found in (27). To form Ge-ncs, the heat treatment of the films was performed in conventional horizontal furnace at T A =500-800 °C during 1-50 min in nitrogen atmosphere.

Optical and structural properties of the films were studied by spectroscopic ellipsometry, Raman scattering, X-ray diffraction (XRD) and photoluminescence (PL) methods. Spectroscopic ellipsometry measurement was performed by means of a Jobin-Yvon ellipsometer (UVISEL, HORIBA Ltd., Kyoto, Japan), where the incident light was scanned in the range of 1.5 to 4.5 eV under an incident angle of 66.3°. The fitting of the experimental data was performed using DeltaPsi2 software (HORIBA Ltd., Kyoto, Japan). The refractive indexes are reported being taken at 1.95 eV. Raman scattering and photoluminescence spectra were recorded with a double grating monochromator equipped with CCD camera. Solid-state laser with the 532-nm wavelength was used as excitation source. The power density of this excitation light was chosen to prevent any heating of the films during spectrum recording. XRD data were collected in a 2θ range of 20°-80° with Cu K α radiation (0.1514 nm) of a Philips X'PERT HPD Pro device. An asymmetric grazing geometry (ω~0.5°) was chosen to increase the volume of material interacting with the X-ray beam. The data were compared with standard cards of Powder Diffraction File Database (#4-0545 for cubic Ge). All measurements were performed at room temperature. The combination of these non-destructive methods allows information about microstructural and optical properties of the films to be obtained.

Results and discussion

Spectroscopic ellipsometry

Being a fast, sensitive, and non-destructive method, spectroscopic ellipsometry is used for the determination of the thickness of thin film and its optical parameters (such as refractive index, n, and absorption coefficient, α). The fitting of experimental data obtained for Ge-Al 2 O 3 films as well as pure Ge and pure Al 2 O 3 counterparts was performed using the approach described in (4,7). For pure Ge and pure Al 2 O 3 films the refractive indexes (taken at 1.95 eV) were estimated as n Ge =4.21±0.01 and n Al2O3 =1.75±0.01, respectively. The fitting of the spectra of Ge x (Al 2 O 3 ) 1-x film deposited on long substrate revealed the variation of the refractive index (taken at 1.95 eV) from n=1.81±0.01 (Al 2 O 3 -rich side) up to n=3.85±0.01 (Ge-rich side) along the film length. The thickness of the films was varied from 450±5 nm (Al 2 O 3 -rich side) up 1300±10 nm (Ge-rich side). This thickness variation was caused by the higher sputtering yield of Ge target in comparison with that of Al 2 O 3 one.

Earlier, we have shown that for composite Si-rich Al 2 O 3 and Si-rich HfO 2 thin films grown by magnetron sputtering, the use of the effective medium approximation permits to determine the relationship between the refractive index n and film stoichiometry x (4,7,30,31). The same procedure was applied in present study for Ge-rich Al 2 O 3 films. Since the experimental spectra were featureless, this allowed assuming that the films are amorphous and can be considered as an effective medium which dielectric constant is determined by the relative contribution of Ge and Al 2 O 3 components as

2 3 2 3 Al O Ge Ge Al O (1 ) 0 2 2 eff eff eff eff x x ε -ε ε -ε + - = ε + ε ε + ε [1]
where ε eff is the dielectric permittivity of composite material, ε Ge and ε Al2O3 are the dielectric permittivity of Ge and Al 2 O 3 , respectively, and x is the fraction of the Ge component. Thus, our Ge-rich films can be considered as Ge x (Al 2 O 3 )

Raman scattering spectra

As-deposited Ge x (Al 2 O 3 ) 1-x films were found to be amorphous. For the films with x>0.6, the peak position of Ge TO phonon was found at ω TO ≈270-274 cm -1 (Fig. 1,a). This spectral position is similar to that observed in amorphous germanium. When Ge content decreases to x=0.1, the broadening of Raman scattering spectrum and shift of the ω TO down to 268 cm -1 were observed. This shift is the evidence that the structure of the films was unrelaxed. This can be caused by the presence of point defects as well as by the tensile stress appeared due to mismatching between the lattice parameters of Ge (a=0.566 nm) and Al 2 O 3 (a=0.792 nm for cubic; a=0.475 nm, c=1.299 nm for hexagonal alumina (32)). It should be noted that Raman scattering spectra for the films with x<0.1 were not recorded due to high transparency of the films. Significant contribution of amorphous Ge phase in the films with x≥0.2 (Fig. 1,a) testifies to the formation of Ge amorphous inclusions during deposition process similar to the effect reported for Si-rich Al 2 O 3 films (4,7,10,27). The heat treatment of the films at T A =500 °C did not affect film structure. Whatever annealing time, all these films kept amorphous nature. The T A rise up to 550°C led to the increase of the intensity of Raman signal and its narrowing accompanied by the shift of peak position towards higher wavenumbers (ω TO ≈272-275 cm -1 ) that is the evidence of the stress relaxation. This transformation of Raman spectra was more pronounced for the films with x<0.6 (Fig. 1,b). The longer annealing time enhanced this effect. It should be noted that these annealed films became more transparent that can be caused by the formation of Ge clusters as well as appearance of GeO 2 due to phase separation process. It is worth to point out that the formation of volatile GeO can be ruled out for this annealing temperature.

Annealing of the films at T A =600 °C resulted in further narrowing of Raman signal and its shift towards ω TO ≈275 cm -1 for x<0.6 (Fig. 1,c). However, this final peak position showed that whatever annealing time, the Ge inclusions are rather amorphous than nanocrystalline. This result demonstrates that in the contrast to the ability of aluminum to enhance Ge crystallization (33), in our films with x<0.6 and annealed at T A =600 °C amorphous Ge phase is stable. This result was similar to those reported in (34) for Ge-Al 2 O 3 /Al 2 O 3 superlattices. At the same time the shift of peak position to higher wavenumbers than that of amorphous Ge is the evidence of full relaxation of tensile stresses and, even, an appearance of compressive ones that can be caused by the transformation of Al 2 O 3 as well as by underlying quarts substrate.

In the films with x>0.6, the Raman peak showed its narrowing and shift towards ω TO ≈295 cm -1 (Fig. 1,c) that is the evidence of the crystallization of Ge clusters. It should be noted that the higher Ge content, the more pronounced Ge crystallization takes place and the shorter time is required for Ge-ncs formation. Further increase of T A up to 700-800 °C caused not only Ge-ncs formation, but also the decrease of Raman signal intensity (not shown here). This phenomenon was explained by the Ge out-diffusion from the films via formation of volatile GeO. The results described above were supported by XRD data.

XRD data

Figure 2,a shows the XRD patterns of as-deposited films with different Ge content. All films demonstrate abroad scattering peak indicating its amorphous nature. Similar XRD patterns were recorded for the films with x<0.5 annealed at T A =550 and 600 °C (Fig. 2,b) At the same time, for the films with x>0.5, the appearance of the features at 2θ=27.23° and 50.0° as well as their enhancement after film annealing at 600 °C was observed testifying to the formation of Ge inclusions via phase separation. This result is in agreement with Raman scattering data. As one can see from Fig. 2, the higher Ge content, the shorter time is required for the formation of Ge crystallites.

Three peaks marked as (111), ( 220) and (311) in Fig. 2b,c are the signature of Gencs formation that are randomly oriented. However, corresponding peak positions are observed at 2θ=27.23, 45.36 and 53.70° for x=0.58 (Fig. 2,b) and 2θ=27.32, 45.49 and 53.80° for x=0.85 (Fig. 2,c) being shifted in comparison with those reported for Ge standard (2θ=27.28, 45.30 and 53.68°). This indicates in the presence of residual compressive stress in the films with x=0.85. This finding is in agreement with the data that obtained from Raman scattering method (Fig. 1,c). The estimation of coherent domain sizes, d, using Scherrer formula revealed that in the films with x=0.58 annealed at T A =600 °C during 15 min d=3.8-4.0 nm. At the same time, for the film with x=0.85 annealed at T A =600 °C during 1 min d=4.9-5.0 nm. Longer annealing time (5 min) resulted in the d increase up to 8.0-8.2 nm. The analysis of the intensities of XRD peaks showed that the Ge-ncs in these films may have elongated shape with large size in (111) direction.

The comparison of XRD and Raman scattering data revealed the discrepancy between Ge-ncs sizes and corresponding Raman peak positions. This is an additional confirmation of the presence of residual stress in Ge-doped Al 2 O 3 films that brings some difficulty in the estimation of Ge-ncs sizes from Raman scattering spectra. Similar effects have been reported for Ge-doped HfO 2 and Si-doped Al 2 O 3 films by different groups (7,8,18,27). To determine real sizes of Si or Ge crystallites embedded in these hosts, the investigation of such films with transmission electron microscopy was performed. These studies for our Ge-doped Al 2 O 3 films are in progress.

Photoluminescence

Photoluminescence (PL) study of as-deposited films revealed a weak PL emission only for the samples with low Ge content (x<0.2). It is peak position was observed at 560-570 nm and this band was found to be similar to the PL band reported in (4,7,35) for the F 2+ centers in Al 2 O 3 host. The heat treatment at T A =500-550 °C did not changed the PL spectra of the films (whatever annealing time and Ge content). Further T A increase up to 600 °C caused an appearance of three PL bands peaked at 580-600 nm, 620-630 nm and 700-730 nm (Fig. 3). Thus, it is clear that several recombination channels contribute to PL emission. The shape and the intensity of PL spectrum are determined by their relative contribution.

Usually, for the films with embedded Si or Ge nanocrystals, the shift of PL peak position with variation of nanocrystal sizes is expected because of quantum confinement effect. Thus, the increase of nanocrystal sizes has to result in the shift of PL peak position towards longer wavelengths (1)(2)(3)(4)7,18). The increase of Ge content from x=0.05 up to x=0.30 caused the non-monotonic variation of the intensity of the PL component at 700-730 nm and its quenching with further x rise (Fig. 3,a). The intensities of other two PL bands changed also nonmonotonically with Ge content whereas for the samples with x>0.3 only these bands were detected being significantly overlapped (Fig. 3,b). Similar PL emission centered at about 600 nm was reported in (36) for large Ge-ncs. It was shown that this PL band is due to Ge oxide shell covered Ge-ncs. Thus, one can assume that in our films with x≥0.58, the PL emission at 580-620 nm (Fig. 3b) is due to Ge oxide shell covered Ge-ncs. This statement can be supported by the fact that the variation of the Ge-ncs sizes from 3.8 nm (x=0.58) up to 8.2 nm (x=0.85) results in the variation of PL intensity only without any significant shift of PL peak position (Fig. 3,b). The higher PL intensity for x=0.58 can be caused by the larger amount of corresponding emitting centers located in Ge oxide shell.

The PL band centered around 700-730 nm was found to be similar to PL emission reported for Ge-ncs embedded in Al 2 O 3 host in (18). In that work, this PL band was ascribed to exciton recombination in Ge-ncs with the sizes d=6.6-7.1 nm. However, simple confinement model did not consider that the carriers confined in the nanocrystals are under a finite potential (18). At the same time, the band at 700-730 nm (1.69-1.77 eV) is significantly shifted in comparison with Ge bandgap (0.66 eV). This fact was explained by the dependence of confinement potential on surrounding matrix caused by the conduction and valence band offsets between the Ge-ncs and Al 2 O 3 matrix. Besides, the modification of Ge-ncs bandgap induced by the polarization interface charge was also taken into account since it may also play an important role due to the difference in dielectric constant of the host matrix and the Ge nanocrystals.

As it was shown by XRD and Raman scattering study described above, our samples did not demonstrate the formation of Ge-ncs in the films with x<0.3. In this regard, one can assume that the PL band at 700-730 nm can be caused by the recombination of carriers in amorphous Ge clusters being affected by surrounded Al 2 O 3 host. At the same time, this PL band can originate also from the defects located in the vicinity of Ge-cluster/Al 2 O 3 interface. Since PL band centered at 700-730 nm was not observed in pure Al 2 O 3 films, one can assume that the defect responsible for such emission is constructed from Ge, Al and O atoms. Based on the obtained results one can concluded that the complex nature of PL spectrum is determined by the contribution of different recombination channels. However, additional study of PL emission in infrared spectral range as well as temperature behavior of PL spectrum can help to discriminate the contribution of carrier recombination via the Ge-ncs and different defects in light emission process. These experiments are in progress.

Conclusions

The evolution of optical and structural properties of Ge-doped Al 2 O 3 films with different Ge content was investigated. As-deposited films and those annealed at T A =550 °C were found to be amorphous whatever Ge content. The formation of Ge nanocrystals was observed after treatment at T A =600 °C during 1-15 min. The higher Ge content, the shorter annealing time was required to form Ge nanocrystals. This annealing stimulated also an appearance of photoluminescence in visible-near infrared spectral range. Three PL bands peaked at about 590-600 nm, 620-630 nm and 700-730 nm were detected. Their relative contribution was depended on Ge content determining the shape and peak position of luminescence spectrum. Two former components were ascribed to the carrier recombination via Ge oxide shell covered Ge-ncs. The band peaked at 700-730 nm was assumed to originate from carrier recombination in amorphous Ge clusters or via defects at Ge-cluster/Al 2 O 3 interface. The approach for the rough express estimation of the Ge excess content in the films based on the application of the spectroscopic ellipsometry technique was proposed.

Figure 1 .

 1 Figure 1. Raman scattering spectra of as-deposited samples (a) and those annealed at T A =550 °C for 50 min (b) and 600 °C for 15 min (c). The Ge content is mentioned in the figures. The spectra are vertically shifted for clarity. The vertical dashed lines show the peak position of Raman signal for amorphous (a-Ge) and bulk (c-Ge) germanium.

Figure 2 .

 2 Figure 2. XRD patterns as-deposited (a) and annealed (b,c) films. Annealing parameters and Ge content in the films are mentioned in the figures. The patterns are vertically shifted for clarity.

Figure 3 .

 3 Figure 3. PL spectra of the films with different Ge content annealed at 600 °C for 15 min in nitrogen atmosphere. The Ge content is mentioned in the figures. Curves 4 in both figures correspond to the same sample. Excitation wavelength is 532 nm.

  1-x medium. Taking into account the n values obtained by the fitting of experimental ellipsometry spectra, the x values of Ge x (Al 2 O 3 ) 1-x films were found to vary from x=0.05 (at Al 2 O 3 -rich side) up to x=0.85 (at Ge-rich side). It should be noted that the films with х>0.7 can be considered rather as a Ge film enriched with Al 2 O 3 , than Al 2 O 3 host doped with Ge.

Acknowledgments

The authors would like to thank for financial support the National Academy of Sciences of Ukraine and the Ministry of Education and Sciences of Ukraine (the project 89452), the National Research, Development and Innovation Office of Hungary (the project OTKA K131515 "Low-dimensional nanomaterials for the optical sensing of organic molecules on liquid and gas interfaces").