Mass transfer enhancement and surface functionalization in digital microfluidics using AC electrowetting: the smaller, the better

This paper addresses modelling and computa-tion of stirring and chemical transport at drop scale due to shape oscillations of small amplitude induced by coplanar AC electrowetting on dielectrics (EWOD). An axisymmetric flow inside a drop


Introduction
Substantial benefits of microfluidics microsystems were revealed during the last decades. Low-cost, portable or high-throughput microdevices were developed offering the possibilities to handle small amounts of biological samples and to perform chemical reactions with small time response for biosensing and environmental monitoring (Pollack et al. 2011;Delattre et al. 2012). Due to the small length scale of liquid samples, biosensing performances are mostly limited by a small Peclet number or equivalently by the diffusive transport of the reagents (antigenes, antibodies, enzymes, cells, functionnalized beads . . . ) to the sensing surfaces. Because chemical or biological processes are commonly diffusion limited, stirring stands as a worthwhile mechanism for speeding up biosensing. At microscale, the Reynolds number is small and the stirring based on turbulent flows seems, if not impossible, at least very difficult to achieve. One first reason for this is the recourse to internal obstacles when the goal is to promote passive mixing. Alternatively, if the goal is to promote active mixing, a second reason is the driving source of the flow, which commonly relies on body forces or pressure gradients. Current means of realizing active mixing or active stirring at microscale were increasingly developed during the last decades (Hessel et al. 2005;Capretto et al. 2011) making use of external sources of energy such as ultrasounds, acoustic waves or surface acoustic waves, electrokinetics, electrohydrodynamics, electrowetting on dielectrics (EWOD), piezoelectricity, electrothermal effects, magnetohydrodynamics, micro valves, steady or unsteady pumping . . . . Despite the fact that with miniaturization, the surface to volume ratio becomes increasingly high, only few of these external sources (electrokinetics, EWOD, Taylor-Melcher electrohydrodynamics, see for instance Chang and Yang 2007;Ramos 2011) exploit surface phenomena (electric double layer, electric surface density) in order to generate internal flows at microfluidics scale. In microfluidics applications, steady microstreaming is one of the rare phenomena based on inertial effects, which does not require a high value of the Reynolds number. When it is generated by means of acoustic waves (Sadhal 2012) or surface acoustic waves (Li et al. 2009), it can be used to pump a liquid, to trap biochemical reagents within vortices, enhancing thus chemical reactions (Lutz and Chen 2003), to favour stirring or mixing in microchambers (Liu and Yang 2002;Wiklund et al. 2012) with sometimes the objective to enhance biosensing performance (Ducloux et al. 2010).
• This paper focuses on coplanar oscillating electrowetting on dielectrics (coplanar AC EWOD) widely used in the development of modern digital microsystems (Pollack et al. 2011). In contrast to the common EWOD configuration with a needle electrode, there is no need for taking into account the capillary effects due to a moving meniscus along the needle. This paper addresses the ability of AC EWOD to enhance mixing and to modify surface ageing of a sessile micro-drop in the coplanar configuration. By using the generic term surface ageing, it is suggested that liquid and solid surfaces, eventually functionnalized, can adsorb target reagents initially solubilized in the drop. Heterogeneous biosensing is thus thought of as a possible example of surface ageing. By focusing on surface ageing of both the drop surface and the supporting solid surface, we investigate the impact of flow patterns generated by the oscillatory drop deformation at different mode numbers, n = {2, 4, 6}. Mass transfers and surface ageing are computed by a finite element method, benchmarked with one-dimensional analytical calculations .
Base configurations are investigated: a first one with only the liquid drop surface behaving as an active surface, and a second one with only the wetted solid surface behaving as an active surface. Then, the competition between both active surfaces is finally investigated. Axisymmetry is assumed in this paper with no significant loss of generality because we know that whatever the geometry of the electrode pair, oscillating electrowetting remains radially isotropic provided that: • the radius of the wetting area remains much larger than the electrode gap while being smaller than the outer boundary of the electrode pair (Malk et al. 2012), • the electric Bond number, defined as the ratio between the electric stress at the electrode gap and the capillary stress, remains small enough, Bo e = εV 2 γ g ≪ 1, in order to avoid drop pinch-off at the electrode gap (Davoust et al. 2013) (g, V , ε and γ denote, respectively, the electrode gap, the voltage, the absolute dielectric permittivity of the drop and the surface tension), • the oscillation amplitude remains much smaller than the drop radius (Hong et al. 2012). Note that this last assumption is hopefully consistent with the following perturbative theory based on a large Strouhal number. • These conditions are fulfilled in standard electrowetting applications with millimetre drops in air deformed by shape oscillations of about 10 µm, surface tensions between 10 and 100 mN, millimetre electrodes, an electrode gap of about 10 µm and an actuation voltage ranging from 10 to 100 V RMS (see, e.g. Davoust et al. 2013 for more details on possible sources of symmetry breakup).

Large Strouhal number approximation
Streaming is defined as a net steady flow of a liquid, induced by the oscillatory motion of its wavy boundary surface. It is made possible due to (inertial) nonlinear terms in the Navier-Stokes equation (Stuart 1966;Riley 1965;Batchelor 1967). Considering the literature devoted to oscillating drops actuated by coplanar AC electrowetting Mugele et al. 2011;Malk et al. 2011Malk et al. , 2012Davoust et al. 2013), the system under consideration here is an incompressible sessile drop deformed by a spherical capillary wave, oscillating radially at the eigen frequency ω n during the capillary resonance of the oscillatory mode n (Fig. 1). In this case, the normal momentum at the surface is balanced by the normal viscous stress and the surface tension. The tangential momentum along the surface is balanced by viscous shear and possibly the damping Marangoni force due to stretching of the aged drop surface. We consider here that these momentum balances at the surface are taken into account from the outset: the focus is put on drop oscillation-induced streaming, and especially its weak coupling with chemical transport in a time-averaged liquid half-sphere of radius R (Fig. 1). At the resonant mode, n, we therefore introduce the timescale ω n −1 , always considered as very small with respect to all other timescales in the following, the magnitude of the normal velocity at the drop surface, U osc,n , and the time-averaged drop radius, R, as characteristic scales to write the Strouhal and Reynolds numbers: and where η and ρ are the Newtonian viscosity and the density of the drop, respectively.
For a two-dimensional axisymmetric flow in spherical coordinates (R, θ , ϕ) with no dependence of the velocity on the azimuthal angle θ (ϕ is referred to as the colatitude, see, e.g. Fig. 1), it is worthy to introduce the stream function, Φ, such that the 2-D velocity writes as u = ∇ × (e θ �(r,ϕ) R cos(ϕ) ) in the cross-section (R, ϕ) and continuity is automatically satisfied. 1 By introducing the only nonzero component of vorticity, � = e θ .(∇ × u), the vorticity equation writes in terms of the stream function: S = ω n R U osc,n , Re = ρRU osc,n η , 1 However, for experimental work using coplanar electrowetting (work referred to in this paper), it is worthwhile to note that the gap g defining the geometry between two half-circular electrodes is required to be much smaller than the drop diameter (g/2R < 1 %) in order for axisymmetry to be guaranteed (Malk 2011). with the following dimensionless numbers: and with the differential operator E 2 , When ξ ≪ 1, or equivalently S ≫ 1, it is possible to account for the nonlinear terms in Eq. (1) by using a regular perturbation method to find a solution (see Schlichting and Gersten (1979) Injecting the latter series into the vorticity Eq. (1), a system of equations is obtained, which was solved by Riley (1967) in cartesian geometry. For ξ = 0, the vorticity at leading order Φ 0 is found to be the harmonic solution of the associated linear equation. The first-order term, Φ 1 , contains both an unsteady component, Φ u 1 (t) and a steady component, Φ s 1 ; the latter represents the net microstreaming flow (Stokes drift) of interest in this paper: This solution exhibits a first component of unsteady vorticity confined to a Stokes boundary layer, of thickness δ s = η ρω, whose time-averaged value is again zero. At the outer boundary of the Stokes layer, in the inner of the drop, the second component is responsible for a net steady streaming flow (also referred to as drift flow in the literature), which is characterized by the steady tangential velocity at the outer boundary of the Stokes layer, (a) (b) Fig. 1 Spherical and 2-D axisymmetric cylindrical coordinates used in this paper, (R, θ, ϕ) and (r, z) respectively, with r = R sin(ϕ) and √ z 2 + r 2 = R. Note that ϕ is referred to as either the colatitude or the polar angle, depending on the coordinates considered in this paper.
a Side view drop shape deformation for mode n = 4. Inset sketch of the Stokes layer (thickness: δ S ) and the Stuart layer (thickness: δ Stu ) standing along the oscillating drop surface. b Top view geometry of the coplanar electrode pair topped with the centred drop This last expression can be obtained by making use of complex notations (see Riley 1965;Batchelor 1967 for instance). Here, ϕ is the colatitude, U osc,n is the (small) amplitude of the spherical capillary wave, and ψ n is the phase lag. The first term on the right hand side of (2) corresponds to the case of a drop deformed by a standing capillary wave, while the second term corresponds to the case of a traveling capillary wave. For R s ≫ 1, Stuart (1966) shows that the latter steady streaming flow persists beyond the Stokes layer over the typical depth δ Stu = O( η ρω n ξ). For R s ≪ 1, it is induced in the whole drop volume.

Drop (micro)streaming
With electrowetting in the oscillatory regime, coplanar EWOD sustains a single spherical standing wave along the drop surface with capillary resonance properly achieved at the apex, provided that the wave frequency is small enough (n 6 Davoust et al. 2013). Hence, only the first term of the right hand side of Eq. (2) is relevant here. In the resonance state of the eigenmode n, the deformed shape of a drop oscillating under such a standing spherical wave at the frequency, ω n , is given by Kang et al. (2008) in terms of the colatitude ϕ and the time-averaged drop radius R: Here, a n is the amplitude of the standing wave, ψ n is the associated phase lag, and P n is the nth-degree Legendre polynomial. The velocity due to such an oscillatory motion is calculated using the flow potential properly expressed in the inner of the drop (Prosperetti 1974): Hence, the time-dependent surface velocity of the associated standing wave writes: with the following amplitude: Hence, according to Eq. (2), for a drop only deformed by a standing wave, the steady streaming flow at the outer of the Stokes layer writes: (2) u s,n = − 3 8Rω n dU 2 osc,n dϕ + 2U 2 osc,n dψ n dϕ .
To calculate the internal drop convective flows, the latter expression can be applied as a boundary condition along the time-averaged drop surface (Fig. 2), provided that the normal surface deformation and the thickness of the Stokes layer remain much smaller than the drop radius. Typically, the standing nature of the drop oscillations is a valid assumption if the viscosity of the drop and the frequency remain small enough (ω n ∼ 10-100 Hz for a water drop), so that there is no significant damping of the mechanical energy of the surface wave. Fortunately, such conditions are not restrictive since they are currently encountered in applications of oscillating electrowetting. Note also that a typical scale for the streaming velocity writes: U s,n = 3(a n ) 2 ω n 8Rn 2 .

Scaling analysis
It is now worthwhile checking how coplanar EWOD experiments comply with the latter assumptions of both a large Strouhal number and a small thickness of the Stokes layer. The aim is to make sure that it is consistent to impose a streaming velocity along the drop surface. To this end, we focus on the most energetic mode (n = 2), for which the two latter assumptions could be difficult to fulfil. The classical situation of a 1.5 µL drop is addressed, oscillating with a mode n = 2 in the resonance state, the Fig. 2 Surface velocity, u s,n (ϕ), along the time-averaged drop surface for modes n = 2, 4 and 6 (drop apex: ϕ = 0, triple contact line: ϕ = π/2). For negative values, the streaming velocity is directed towards the drop apex along the drop surface eigenfrequency and the amplitude to radius ratio are typically, ω 2 = 1500 rad/s and a 2 /R = 0.1, respectively (R ∼ 1 mm, a 2 ∼ 100 µm at most, see Davoust et al. 2013  which exhibit an agreement between experimental imaging of tracer particles and numerical simulations based on the present streaming model. The velocity of 15.8 mm/s is issued from theory and experimentally verified for a common oscillation regime.

Flow pattern as calculated for the most energetic modes
The calculation of steady streaming flow is performed for a laminar regime in a half-spherical (time-averaged) drop with the Stokes layer disregarded. Whichever the mode number 2 (n = 2, 4 or 6), it is checked experimentally and numerically that the thickness of the Stokes layer, though larger than the size of the biological molecules considered, is always much smaller than the thickness of the Stuart layer. The boundary conditions for the flow velocity u are as follows: • a symmetry condition along the vertical axis (r = 0), ∂u ∂n = 0, with n, the normal vector, • a slip condition 3 along the base axis (z = 0), ∇u.t = 0, with t, the tangential vector, together with an impermeability condition, u.n = 0, • and a moving wall condition, u.t = u s,n , along the (timeaveraged) liquid surface (R = √ r 2 + z 2 ) where u s,n is the streaming velocity (see Eq. 5) prescribed along the drop surface (Fig. 2).
• The steady momentum so imposed from this last boundary condition along the drop surface is finally transmitted by viscosity in the inner of the drop.
By switching the latter moving wall boundary condition along the drop surface from one resonant mode to another one ( Fig. 2), flow patterns are accordingly found to be a series of n superimposed toroidal vortices. In all cases, in the meridian cross-section investigated here (r = 0 or ϕ = [0 . . . π/2]), the upper toroidal vortex near the drop apex is always found counter-clockwise and successive vortices alternate between clockwise and counter-clockwise (see Figs. 3,4).
In Fig. 4a-d, calculations are performed for coplanar electrowetting experiments in typical conditions with water-based or PBS solutions at room temperature (viscosity: η = 1 mPa s, density: ρ = 1000 kg/m 3 ).

Chemical transport at drop scale
In contrast to a previous paper by the present authors , devoted to the coupling between evaporation processes, diffusion and adsorption of surfactant molecules in a drop under static conditions, steady streaming is now introduced in this paper due to AC electrowetting. This time, chemical gradients are no more isotropic since they are strongly modified by internal flows: streaming current transports molecules and thus creates a non-isotropic distribution of the chemical concentration in the bulk as well as along the surface of the drop. Several effects might influence the efficiency of molecular transport, namely the sorption kinetics and the total amount of adsorbed molecules. First of all, the streaming current might accelerate the dispersion of molecules throughout the droplet, stimulating stirring, to outrun diffusion. Secondly, stirring might help to reinforce the surface concentration in one or several locations (rings) corresponding to stagnation points of streaming-induced internal flows. This might lead to flow focussing of adsorbed molecular complexes. Since this could be expected to cause surface concentration gradients, the effect of molecular 2-D diffusion along the drop surface may no longer be neglected.

Chemical transport within the drop volume
In order to prepare the presentation of chemical surface transport, and to estimate whether streaming is able to outrun diffusion or not, this section gives some statements on how diffusive and convective transport act in the drop volume. The considered system is still a halfsphere with the radius R, but now use is made of 2-D axisymmetric cylindrical coordinates (r, z) normalized by the drop radius R (Fig. 1). The diffusion flux, j diff , is commonly given by the non-dimensionalized Fick's first law: where the ratio of timescales and the magnitude of finiteness  are defined, respectively, as η Td = T /τ diff , and ε 0 = C 0 R/Γ e . Here, C 0 is the initial surfactant volume concentration as well as the concentration scale for defining the normalized bulk concentration, C. Use is made of the arbitrary timescale T selected among the possible timescales governing surface ageing kinetics (discussed later). The diffusion timescale is defined as τ diff = R 2 /D, with D the bulk diffusivity. The scale for the surface concentration, Γ e , as calculated from Henry's law at thermodynamical equilibrium, Γ e = C 0 k a /k d , is used to normalize the surface concentration, Γ , with k a and k d the adsorption and desorption constants, respectively. The unsteady balance between the divergence of the diffusion flux and the divergence of the streaming-induced molecular transport writes in non-dimensional form, where the ratio of timescales, η T s = T /τ stream , is defined from the dynamic timescale, τ stream = R/U s .
The non-dimensional initial condition writes, A first boundary condition writes, along the vertical axis of symmetry, Depending on whether the surface under consideration is subjected to biological ageing or not, two other boundary conditions can be written either at the wetted area (z = 0), or at the drop surface ( √ z 2 + r 2 = 1), The terms j ad, 0 and j ad, 1 are referred to as the ad/desorption flux (see following section for further details on j ad,1 for instance), and the Thiele number, Th, is defined as Th 2 = τ diff /τ ads = R 2 k d /D, with τ ads = 1/k d , the timescale of sorption processes, and k d , the desorption coefficient in the Henry's law. Because use is made of Henry's law, it is found that the time scales for adsorption and desorption processes τ ads and τ des are same.

Chemical transport along the drop surface
The equation for surfactant transport along a drop surface as presented in Theisen and Davoust (2012) does not reflect the whole truth in this case: steady streaming as well as surface diffusion must be taken into account by considering the mass balance and the non-dimensional surface concentration Γ . As shown in Stone (1990) for instance, a general transport equation along a moving deformed interface writes in its non-dimensional form, with u σ s , the 2-D (in-plane) surface velocity along the deformed drop surface, and j ad,1 , the nondimensional sorption flux at the drop surface: . Note that the surface velocity u σ s must be distinguished from the streaming velocity, u s,n t, along the time-averaged drop surface.
Here, η T a is a dimensionless ratio defined by the timescale of sorption processes, τ ads . The dimensionless ratio, η T ,stream and η T ,diff s are defined from the timescales of steady streaming and surface diffusion, τ stream = R/U s and τ diff,s = R 2 /D s , respectively.
The surface del operator ∇ s is the del operator projected onto the tangent plane at any point of the deformed surface. To this end, use is made of the tensor operator, (11 − n ⊗ n) applied either to the volume velocity, u, or to the del operator, ∇, all of them being valued at the deformed surface: with n the normal vector at the deformed surface, and ⊗ the tensor product.
It must be kept in mind that in this paper, the Strouhal approximation is done with a negligible thickness of the Stokes layer. For sake of consistency with streaming theory, for which the (thin) Stokes layer is governed by a zero net flow, surfactant convection along the time-averaged surface is finally not taken into account (no retraction nor elongation along the time-averaged drop surface). Hence, the second and third (convective) terms are disregarded in Eq. (6): u σ s = 0. The fourth term in Eq. (6), Γ (∇ · n)(u · n), represents inflation/deflation of the surface, and it can be again disregarded in this paper since the drop geometry is time-averaged. 4 Nevertheless, steady streaming is expected to have a significant impact upon surface transport because it can modify the ad/desorption source term, j ad,1 , in Eq. (6): drop convection causes gradients of chemical concentration along the drop surface via the ad/desorption term described by Henry's law (Γ ∼ C). This is also why it is necessary to consider surface diffusion. Hence, the transport equation (6) includes a term for surface diffusion (Picard and Davoust 2006), η T ,diff s ∇ s · (∇ s Γ ), equivalent to the 2-D projection of Fick's second law. From a physical point of view, this means that any streaming-induced surface gradient of the surface concentration can be damped if the surface diffusion term is large enough.
Finally, the useful form of Eq. (6) considered in this paper writes accordingly: Validity tests of its numerical implementation are provided in Appendix 2.
4 Coupling chemical transport phenomena: the role of all relevant dimensionless numbers for the most energetic mode, n = 2 Among remaining physical phenomena to be coupled, there is steady drop streaming, governed by the outer boundary condition (Eq. 5) for the streaming velocity, and bulk diffusion as governed by Fick's law. 5 Furthermore, bulk and surface transport mechanisms are coupled by Henry's law for ad-and desorption processes. The adsorbed surfactant molecules can also be transported under surface diffusion along the drop surface from high-to low-concentrated spots. Streaming, volume diffusion, sorption processes as well as surface diffusion imply their own timescales, τ stream = R/U s , τ diff = R 2 /D, τ ads = 1/k d and τ diff,s = R 2 /D s , respectively, all of them being related to each other by means of the following non-dimensional numbers: • The Thiele number, • Clearly, only three of them could be sufficient in describing the whole system since for instance, Th 2 = Pe. Da. But our goal is to present explicitly the contribution of each kinetics from non-dimensional numbers well referenced in the literature.
The analysis of streaming efficiency is only relevant for diffusion-limited systems, because otherwise the concentration gradient at the vicinity of the drop surface would be too weak for streaming-induced convection to have a significant effect on surface ageing. This means that the Thiele number is assumed to be large in the following analysis (Th = 10 3 ).
Finally, as regard the transient chemical regime considered in this paper, the initial conditions for volume and surface concentrations write, respectively: 5 Note that the Stokes layer is strongly stirred, and therefore, its chemical impedance can be disregarded.

Peclet number
The Peclet number, Pe = τ diff /τ stream , characterizes a system regarding the efficiency of (streaming-induced) convection versus the efficiency of diffusional transport. Systems with a high Peclet number are generally governed by convection, whereas a low Peclet number is preponderant in systems in which diffusion acts on a smaller timescale than convection. Figure 5 displays the normalized variance, σ 2 C /C 0 , of the volume concentration as a function of time, t/τ diff , when surface diffusion is neglected: Pe Stu → ∞. The variance bears information on the mixing of solute molecules (Liu and Haller 2004;Rothstein et al. 1999) and is defined as follows: where N vol (t) := V C(r, z, t)dV corresponds to the amount of solute molecules in the volume, 6 and where C (t) := N vol (t)/V is the mean concentration. The graph states that for a high Peclet number, the streaming flow generates a high variance of volume concentration, which falls quickly to zero at times about 0.01 to 0.05 τ diff : EWOD-induced mixing is efficient, and the system tends to quickly homogenize molecular distribution. Intermediate systems with a Peclet number close to Pe 1 are hardly distinguishable from purely diffusional regime, Pe ≤ 1, for C(r, t = 0) = 1 and Γ (t = 0) = 0. which σ 2 C does not depend on the Peclet number. It is entirely governed by diffusion, and convection does no longer play a significant role: the relaxation time of σ 2 C is close to the diffusion timescale, τ diff .
Practically speaking, convection and diffusion act as follows in the two regimes: • For high Peclet numbers (Pe = 10 2 -10 3 , τ diff ≫ τ stream ), adsorption of surfactants to the drop surface will create a depleted boundary layer at the vicinity of the drop surface (Fig. 6). This boundary layer can be convected by the streaming current towards the apex and the triple contact line, respectively. In this case, the streaminginduced counter-rotating vortex pair is responsible for a beautiful mushroom-like pattern, as demonstrated in Fig. 6d. The diffusion boundary layer remains of small thickness, as diffusion is outrun by convection. • For low Peclet numbers, adsorption also takes place, but this time the boundary layer is not convected. The socreated concentration gradient generates diffusive flux from the drop centre to the surface and thus sustains the adsorption flux. This extends the boundary layer more and more, until the drop volume is emptied from solute molecules (if ǫ 0 small enough).
• For intermediate Peclet numbers, convection and diffusion exhibit approximately the same timescales and the boundary layer expands to the centre of the droplet, while being convected slowly at the same time.
Regarding the curve for a very high Peclet number, Pe = 10 4 (τ diff ≫ τ stream ), a shoulder arises around t/τ diff ∼ 10 −3 , right after the local maximum at t/τ diff ∼ 10 −4 . The emergence of this shoulder is induced by the counter-rotating (toroidal) vortex pair, which refills transiently-after one first rotation of the vortex pair (t ∼ τ stream )-the chemical boundary layer near the drop surface. This is illustrated by an outwards jet along the axis ϕ = π/4 ( Fig. 6b-d). Figure 7 gives a first overview of the surface ageing process in function of the Peclet number though it needs to be taken with care, as the Damköhler number, Da, is not yet taken into account. Nevertheless, Fig. 6 Concentration plots during transient regime, at t = 0.5 τ stream and t = 1.5 τ stream , for mode n = 2, Th 2 = 10 3 and Pe Stu → ∞, and two Peclet numbers: for a high Peclet number, Pe = 10 4 , the boundary layer is convected right away and has little time to develop. For an intermediate Peclet number, Pe = 10 2 , the boundary layer develops further. The concentration is coded in grey scale: white stands for C = 1, black stands for C = 0, separated by fifteen equidistant levels on the grey scale. a Pe = 10 4 , time t = 0.5 τ stream , b Pe = 10 4 , time t = 1.5 τ stream , c Pe = 10 2 , time t = 0.5 τ stream , d Pe = 10 2 , time t = 1.5 τ stream (a) (b) (c) (d) Fig. 7 indicates that the time required to reach thermodynamical equilibrium is sped up in systems with a high Peclet number. Further analysis needs to be made in order to confirm these preliminary results.

Damköhler number
The Damköhler number, Da = τ stream /τ ads , characterizes a system regarding the efficiency of streaming-induced convection compared to adsorption. In order to facilitate the analysis, we focus on the case Pe = 10 4 (and Th ≫ 1). The surface Peclet number exhibits a high value, Pe Stu → ∞, and hence, surface diffusion is again neglected here. In Fig. 8, the non-dimensional variance of the volume concentration, σ 2 C /C 0 , is displayed for several values of the Damköhler number, Da = 0.01-100. The peak of variance can be observed and is approximately situated at t/τ stream ≈ 1, indicating that in this case, mixing kinetics is essentially dictated by streaming. Since the system is qualified by Pe = 10 4 , essentially the same mixing behaviour as in Fig. 5 (black circles) can be observed.
In the case of a diffusion-limited system (Th ≫ 1: adsorption kinetics is quasi-instantaneous), streaminginduced convection is able to provide a faster surface ageing provided Pe remains large enough (τ ads ≪ τ diff and τ stream ≪ τ diff ). Now, if Da is small enough (τ stream < τ ads ), surface ageing enhancement by streaming is secured because diffusion transport is by-passed and the chemical flux required by sorption processes is sustained by convection: thermodynamic equilibrium is free to take place since it is only governed by the timescale of streaming-induced stirring. This is also found when representing the total amount of adsorbed molecules as a function of the nondimensional time with T = τ ads (Fig. 9). Time for complete adsorption of the major amount of surfactant molecules is clearly shortened when reducing Da.

Surface Peclet number
The surface Peclet number, Pe Stu = τ diff,s /τ stream = RU s /D s , characterizes the balance between streaming-induced convection and surface diffusion; this is a special focus of the present paper. For a high surface Peclet number, surface diffusion is slow compared to streaming-induced effects, and hence, the distribution of adsorbed surfactant molecules is expected to be strongly influenced by the steady streaminginduced flow pattern. A low surface Peclet number means a fast surface diffusion kinetics, such that surface gradients of surface concentration, ∇ s Γ , due to streaming should be efficiently damped. It is also worthy to note that whichever the value of Pe Stu , the timescale for achieving thermodynamical equilibrium is independent of Pe Stu , a fact which is due to the application of Henry's law (data not shown). As shown in the following, some of these expectations are clearly demonstrated in Fig. 10 for a diffusion-limited system, Th 2 = 10 3 , with a large Peclet number, Pe = 10 4 , and a rather small Damköhler number, Da = 0.1: τ stream ≪ τ ads ≪ τ diff (surface ageing enhanced by streaming).
It is also worthy to note that in order to observe a dependence on Pe Stu , the focus can be placed on the modified surface concentration of surfactant molecules, Γ µ . In digital microfluidics, the modified surface concentration, Γ µ , originally introduced by Theisen and Davoust (2012), permits to take into account the finite size of the drop and therefore the possibility to change surface ageing kinetics when the amount of adsorbed molecules becomes of same order as the amount of molecules within the drop volume C /C 0 . The system is characterized by a diffusion-limited regime (Th 2 ≫ 1), and a strong mixing (Pe = 10 4 ). Surface diffusion is disregarded (Pe Stu → ∞). Time is scaled by T = τ stream (ε 0 ∼ 1). When ε 0 assumes small to moderate values (liquid drops at microscale), the surface concentration at thermodynamical equilibrium must be corrected by considering the mass balance of molecules in its dimensionless form. Based on these considerations, Theisen and Davoust (2012) introduce a corrected equilibrium surface concentration at drop scale, reflecting the high surface to volume ratio in digital microfluidics. Accordingly, the normalized surface concentration, Γ µ , in the following figures is consistently written as: The modified surface concentration of surfactant molecules, Γ µ , is represented in Fig. 10 for different surface Peclet numbers, Pe Stu = 10 −2 to Pe Stu = 10 5 , as a function of the polar angle, ϕ. Surface concentration profiles are displayed every time step of 20 τ stream . Surprisingly, three regimes are revealed: • For a low surface Peclet number, Pe Stu 1 (τ stream < τ diff,s ), surface diffusion is efficient and the distribution of molecules is consistently found nearly horizontal. All concentration gradients due to steady streaming along the surface are quickly damped under surface diffusion, such that the dimensionless surface concentration, Γ µ , is found homogeneous. • For a surface Peclet number of the order Pe Stu ≈ 10 2 to Pe Stu ≈ 10 3 , a nonuniformity is initially imposed by streaming flow, but beyond a delay of several τ stream , it is further damped with time. The system recovers the latter steady uniform distribution obtained with a low surface Peclet number. • For a large enough surface Peclet number, Pe Stu 10 4 (τ diff,s ≫ τ stream ), surface diffusion is not efficient enough to dissipate surface gradients imposed by the streaming flow. The nonuniformity of the surface concentration profile grows and generates densified spots, which are not dissipated after thermodynamical equilibrium (dN s /dt = 0) is reached, due to the slow process of surface diffusion. This stands as an interesting finding strengthening the conclusion that thermodynamical equilibrium is independent of surface diffusion effects. Again, for highly loaded surfaces, this could be questionable since the Henry's adsorption law used in this paper does not take into account steric effects. This last regime with a very large surface Peclet number illustrates clearly the mechanism of flow focussing. As a typical application, it could be applied to biological molecules: for instance, bovine serum albumin (BSA) exhibits surface Peclet numbers of approximately Pe Stu ≈ 10 4 -10 5 .
In the analysis based on the volume Peclet number, Pe, it is shown how steady streaming causes the transport of a volume concentration gradient along the streamlines. For a high volume Peclet number, Pe, the adsorption flux is mainly supplied by the streaming-induced recirculating flow (not by the diffusive flux). As already mentioned when analysing the variance of the volume concentration for the mode n = 2, the counter-rotating vortex pair is loaded with surfactant molecules. On its way from the drop centre to the surface, along the direction ϕ = π/4, this converging flow delivers molecules in the vicinity of the stagnation point ( √ z 2 + r 2 = R, ϕ = π/4) where it also splits into two outwards flows along the drop surface: the first surface flow is directed towards the apex and the second one is directed towards the triple contact line. With such a particular flow pattern, it is therefore consistent to notice in the early time a bump centred at ϕ = π/4 on the surface concentration profiles (see Fig. 10 when Pe Stu = 10 3 -10 5 ).

What about adsorption to the bottom? "Real life" hybridization assays
Many biochips and microdevices for molecular recognition at drop scale with or without EWOD technology [surface plasmon resonance (Malic et al. 2009), field effect transistor (Choi et al. 2012), cantilever detection (Leichle et al. 2006) or microresonators (Royal et al. 2013), QCM detection (Fuchiwaki et al. 2014)] rely on adsorption of ligand molecules to a biochemical receptor grafted to a solid substrate. Stirring is particularly useful to enhance biochemical transfers and gain a quicker response of the transduction mechanism. By stirring, faster ligand-receptor binding may be achieved, especially if the heterogenous reaction is diffusion limited (Th ≫ 1). First, adsorption is enabled at the solid substrate, but no more at the drop surface. The system is again characterized by the following dimensionless numbers: Th 2 = 10 3 , Pe = 10 4 , Da = 0.1 (τ diff ≫ τ ads > τ stream ) and Pe Stu = 10 2 (τ diff,s ≫ τ stream ). In Fig. 11a, the curves for the amount of adsorbed molecules are shown for the first three resonant modes, n = 2, 4 and 6. In the case of a solid substrate, stirring efficiency of the different resonant modes is strongly different from what happens with a liquid surface: modes 2 and 4 exhibit the fastest kinetics, and surface ageing in the case of mode 6 requires more time to reach thermodynamical equilibrium. The reason for this lies in the internal flow pattern and more specifically in the properties of the Stuart layer, as suggested by the analysis of the streamlines (Fig. 4). The clockwise vortex near the solid surface contributes mainly to mixing and to solid surface ageing enhancement, before the concentration gradient is propagated to the rest of the vortices. Whereas for mode n = 2, the concentration gradient along the solid substrate is convected radially inwards up to the middle of the drop, since the clockwise vortex is of large extent. For mode n = 6, it is concentrated near the contact line due to the reducing of the thickness of the Stuart layer.
Regarding ageing of the drop surface, it is worthy to note that all vortices contribute equally to stirring and molecular renewal, and a streaming-induced concentration gradient is therefore sustained all along the drop surface. The influence of higher-order modes (n = 4, 6) on drop surface ageing and drop mixing is presented in details in Appendix 1.
Furthermore, it is interesting to know whether steady streaming can have an effect on the bulk/surface repartition of the molecular content when, this time, adsorption is permitted both at the drop surface and at the solid substrate. Amounts of adsorbed molecules are shown in Fig. 11b for equal adand desorption coefficients, k a and k d , on both interfaces. The equilibrium (steady) value of the amount of adsorbed molecules on the respective interfaces is different, since the available developed surface for the drop surface (top) is greater than for the liquid/solid interface (bottom). A second point, which is not completely intuitive, is that the timescale related to ageing of the solid surface is not strongly altered by streaming for a specific resonant mode, while it is clearly reduced as far as ageing of the liquid surface is concerned. This is especially true for high-order modes. The fact that the thickness of the Stuart layer is reducing with a growing resonant mode contributes to increase the proximity of the streaming-induced mixing near the drop surface, with subsequent surface ageing enhancement for modes n = 4 and 6.

Conclusion
In a companion paper , a foregoing analysis was performed on unsteady diffusion and sorption phenomena at drop scale with finite size effects included. In this paper, this analysis is extended to the case of a drop oscillating under coplanar AC electrowetting with the arising of a steady streaming-induced internal flow provided at a small enough frequency ( f ∼ 100-500 Hz for a 1 µL drop). Drop oscillations are considered to be standing and characterized by a small amplitude, so that the large Strouhal approximation can be fairly developed.
The approximation of a very thin Stokes layer standing along the oscillating drop surface allows us to apply a moving wall as a steady boundary condition at the time-averaged drop surface for flow calculations in spherical coordinates. The way the drop streaming flows enhance chemical transfers within (and between) the bulk and the surface of the drop is analysed, based on unsteady convective-diffusive transport equations written for the drop volume as well as for the drop surface, coupled each other by serial coupling between diffusive flux and sorption flux at the drop surface. This mathematical model is numerically solved by a finite element method in 2-D axisymmetric cylindrical coordinates. A dedicated weak form is developed especially for including surface diffusion effects when solving the transport equation along the drop surface. The numerical implementation of the model is successfully compared to the foregoing analytical and numerical results established in Theisen and Davoust (2012), when the drop does not experience oscillating electrowetting and subsequent streaming flows. The impact of the drop oscillating modes is analysed using an approach based on four dimensionless numbers: the Peclet number, the Thiele number, the surface Peclet and Damköhler numbers.
Several conclusions can be drawn. Surface ageing enhancement from oscillating EWOD is expected to be particularly reliable for diffusion-limited systems (Th ≫ 1), especially for smaller and smaller drops (R → 0) because the Peclet number remains unchanged (potentially large) while the Damköler number reduces according to Da ∼ R 2 . This is essentially due to steady streaming at large Strouhal number which behaves as a surface generator of momentum particularly suitable for stirring in large S/V microsystems: the streaming velocity u s,n scales as 1/R until the amplitude of the drop oscillations is much smaller than the drop radius, R.
As demonstrated from the variance of the volume concentration, σ 2 C , which gives a hint on the mixing efficiency, higher mode numbers privilege stirring in the vicinity of the drop surface and seem to be more efficient for homogenization of the volume concentration. Conversely, for adsorption of molecules solely on the supporting solid substrate, small mode numbers seem to be preferential. An interpretation based on the Stuart layer is proposed to interpret these findings. Drop streaming is shown to particularly enhance surface ageing at one or several points (rings) along the drop surface, corresponding to stagnation points of drop internal flows. This flow focussing phenomenon could be used with substantial benefit in order to functionalize active spots more efficiently for a fast sensitive detection.
For future work, as a way to include steric effects, it could be worthwhile implementing the Langmuir law rather than the Henry's law used in this paper. Another aspect for future investigations could be the implementation of a stirring protocol at fixed dimensionless numbers, providing more efficient stirring by exploiting for instance a sequence of different resonant modes. Finally, it is believed that the dependence of surface ageing kinetics on the finiteness of the system is an interesting point, which remains to be investigated: a methodology close to the one developed in Theisen and Davoust (2012), based on the dimensionless number ε 0 , can be recommended.

Appendix 1: Drop surface ageing beyond two vortices: higher-order modes
To investigate the dependence of drop surface ageing and drop mixing on the oscillating mode, two successive modes, n = {4, 6}, are investigated with the same maximum Fig. 12 Variance of the volume concentration, σ 2 C /C 0 , as a measurement parameter of the mixing efficiency due to steady streaming. The system is characterized by: Th 2 = 10 3 , Pe = 10 4 , Da = 0.1 and Pe Stu = 10 2 . Time is non-dimensionalized by T = τ stream streaming velocity along the drop surface as the one for the previous mode, n = 2. As suggested by experiments (Davoust et al. 2013), modes of higher order, n 8, no longer follow the streaming theory of this paper essentially based on the presence of a standing capillary wave. The system under consideration here is characterized by the following dimensionless numbers: Th 2 = 10 3 , Pe = 10 4 , Da = 0.1 (τ diff ≫ τ ads > τ stream ) and Pe Stu = 10 2 (τ diff,s ≫ τ stream ).
It is found that the time required to reach thermodynamical equilibrium is reduced with higher-order modes, n = {4, 6}, especially the sixth mode which warrants the fastest kinetics (Fig. 13a). The time required to make the spatial distribution of adsorbed molecules uniform is also lowered in comparison with mode n = 2 as demonstrated in Fig. 13b-d. An indication of how mixing kinetics is modified may also be given by the variance of the volume concentration (Fig. 12) whose profile for the mode n = 6 exhibits the most efficient mixing ability with the smallest peak.
As suggested by the surface velocity profiles (Fig. 2) and the streamlines (Fig. 4b-d), (n + 1) stagnation points (apex and contact line included) can be pointed out along the drop surface in the meridian plane under consideration (ϕ = 0 . . . π/2), of which n/2 are expected to generate a local growth of the surface concentration due to internal flow converging towards the surface. As seen in Fig. 2, the streaming velocity is also damped along the drop surface from the apex to the contact line when n = 4, especially when n = 6. This could explain why only one single bump is clearly observed in Fig. 13b-d, as it is associated with the more intense outward flow among the n/2 stagnation points made evident from a: • bump located at ϕ = π/4 if n = 2 (Fig. 13b), • bump located at ϕ ∼ π/8 if n = 4 (Fig. 13c), the bump expected at ϕ ∼ 3π/8 is questionable, • bump located at ϕ ∼ 3π/32 if n = 6 (Fig. 13d), no bump is clearly observed at ϕ ∼ π/4 and a last one expected at ϕ ∼ 5π/12 seems to arise.
Appendix 2: Implementation of the model without streaming (Pe → 0) and benchmarking procedure

Validity test for sorption phenomena
In this appendix, use if made of the one-dimensional model presented in Theisen and Davoust (2012) as a benchmark for the present model when drop streaming and surface diffusion can be disregarded (Pe << 1). This time, the 2-D axial symmetry is preserved as well as the polar rotational symmetry which is otherwise broken by the streaming current.
The strong form of the surfactant transport equation here used for benchmarking is written using the Henry's law for sorption phenomena (surface diffusion is switched off in the transport equation Eq. (7): with j ad,1 the ad(de)sorption flux at(from) the interface and η Ta , a dimensionless number defined from the timescale of adsorption/desorption equilibrium, τ ads = 1/k d , where k d is the desorption coefficient. The 2-D results of the simulation done after implementing the corresponding weak form, S w(r, z) ∂Γ ∂t − (η Ta (C| r=R − Γ )) dS = 0 in Comsol Multiphysics, are compared with analytical and numerical solutions of the one-dimensional model in Fig. 14: Panel (a) shows the comparison for an adsorption-limited regime with a quasi homogeneous bulk concentration. The analytical solution and the 1-D or 2-D numerical solutions agree perfectly. Panel (b) shows the comparison for a diffusionlimited regime with quasi-instantaneous adsorption/desorption equilibrium. The analytical long-and short-time approximations, referred to as Eqs. (2.23)-(2.24) in Theisen and , are perfectly matched by the present 2-D numerical solutions. This good agreement also suggests a successful implementation of the present 2-D axisymmetric model.

Validity test for mass balance
The mass balance calculated via (N − N 0 )/N 0 , using the total amount of molecules, N(t) = V CdV + S Γ dS, and the initial amount of molecules, N 0 , is constantly found to be smaller than 1 %. This suggests again a successful implementation of the diffusion-adsorption model using ∂Γ ∂t = j ad,1 = η T a (C| r=R − Γ ),

Validity test for surface diffusion
It may now be considered to implement a last physical mechanism, i.e. surface diffusion which is an important ingredient before implementing the effects of streaming current. The reason for this lies in the fact that streaming is expected to generate a concentration gradient along the drop surface, which may complicate the numerical solution. By integrating surface diffusion, all concentration gradients are smoothed, such that the convergence of the calculation is favoured. This means a rougher time stepping and thus less calculation time.
In order to implement and test surface diffusion, ad/desorption processes from the volume to the surface is now switched off by suppressing the term j ads,1 in the transport equation (7). The initial condition for the modified surface concentration is set to Γ µ (ϕ, t = 0) = 2ϕ/π, which corresponds to a linear distribution of the surface concentration with Γ µ (ϕ = 0) = 0 and Γ µ (ϕ = π/2) = 1 at t=0. The result of the simulation is shown in Fig. 15: the concentration gradient is damped during the course of time and a uniform steady value is found to be, Γ µ (ϕ, t → ∞) = 0.64. This is to be compared to the theoretical value found by calculating the total amount of surfactant molecules at the surface, N s : The equilibrium value of the modified surface concentration as analytically calculated writes therefore: N s = π/2 0 Γ µ (ϕ, t = 0)2π R 2 sin ϕdϕ = 4R 2 .
The relative error between the numerical and theoretical values is lower than 1 × 10 −5 , which suggests again a successful implementation of the weak formulation for surface diffusion.