Turning catalytically active pores into active pumps

We develop a semi-analytical model of self-diffusioosmotic transport in active pores, which includes advective transport and the inverse chemical reaction which consumes solute. In previous work (Phys. Rev. Lett. 129, 188003, 2022), we have demonstrated the existence of a spontaneous symmetry breaking in fore-aft symmetric pores that enables them to function as a micropump. We now show that this pumping transition is controlled by three timescales. Two timescales characterize advective and diffusive transport. The third timescale corresponds to how long a solute molecule resides in the pore before being consumed. Introducing asymmetry to the pore (either via the shape or the catalytic coating) reveals a second type of advection-enabled transitions. In asymmetric pores, the flow rate exhibits discontinuous jumps and hysteresis loops upon tuning the parameters that control the asymmetry. This work demonstrates the interconnected roles of shape and catalytic patterning in the dynamics of active pores, and shows how to design a pump for optimum performance.


I. INTRODUCTION
Over the last decades, novel methods to manipulate fluid flows on the micro-and nano-scale have enabled a technological leap [1][2][3][4][5] , with microfluidic devices being employed in fields such as chemical synthesis 6 , and biomedical engineering 7 .The current need for increased sustainability has spawned a new field of chemistry concerned with reducing the environmental footprint of the chemical industry 8,9 .One way is to reduce the size of the chemical reactors themselves.It has been shown that chemical reactors on the scale of dozens to hundreds of micrometers can produce compounds in a way that is safer, faster, and more efficient when compared to traditional batch reactors [10][11][12] .Microfluidic devices can also be used to produce and assemble nanoparticles 13,14 , or for 3D fabrication 15 via e.g.inkjet printing 16 .Biomedical researchers often use microfluidic devices to e.g.develop novel drug delivery systems 17 , or to mimic biological systems such as the vascular system 18,19 .Indeed, microfluidic chips may provide very quick and reliable diagnostics while requiring a minimal amount of biological sample from a patient 20 .
A common need in all of these applications is the need to pump fluid in a controlled fashion.Thus, micropumps are commonly found in many microfluidic devices 21 .A common challenge is the decreased effectiveness of traditional hydraulic pumping arising from the increased surface-to-volume ratio 22 .An alternative is to use surface-driven flows, such as electroosmosis.This phenomenon occurs when a fluid in contact with a charged surface is pumped via an application of an external electric field 23 .The elevated voltage needed to operate an electroosmotic pump can however be an issue 24 .An alternative form of pumping that bypasses the need for an external field entirely is diffusioosmosis 25 .When a solution is placed in contact with a solid phase such as the walls of a channel, solute-solid and solvent-solid interactions play a role.As the solvent is chemically distinct from the solute, these interactions will be of different strength and range.Thus, if the concentration of solute is spatially-varying, the overall interaction between solute and wall will also be spatially-varying, leading to local pressure imbalances near the solid.The end result is a fluid flow parallel to the solute concentration gradient 23 .One way to enforce an inhomogenous concentration of solute is by coating the solid with a catalyst, and to choose a chemically reactive solution.The solution is decomposed in the presence of the catalyst, generating a solute.Covering only part of the solid in catalyst ensures that the solute concentration is inhomogeneous.This technique has been used to fabricate self-motile Janus particles, which are half-coated in catalyst [26][27][28][29][30] .It has also been used to fabricate active, self-diffusioosmotic channels with self-pumping walls 31 .Such diffusioosmotic channels can act as micropumps (i.e.exhibit a non-zero net flow rate) by breaking fore-aft symmetry either via their catalytic coating 32 , or by the shape of their corrugation 81 .Recently, it has been shown that breaking the channel's fore-aft symmetry by design is not necessary, as a spontaneous symmetry breaking occurs when the advection of solute plays a large enough role 34 .It is worth noting that the advection-enabled spontaneous symmetry breaking occurs for colloids as well, leading to diffusiophoretic isotropic colloids [35][36][37] .The interplay of geometrical and chemical inhomogeneities are crucial for the performance of these active pumps, as it has been shown that, while flat channels with homogeneous chemical properties may exhibit local convection currents, they do not amount to a net flow 34,38,39 .
Advection-enabled pumping in fore-aft symmetric channels is expected to occur when the timescale for advective transport matches the one of diffusive transport 34,35 .Using typical experimental values for platinum-coated channels holding a hydrogen peroxide solution, this matching is expected to occur for channels on a scale of a hundred to a thousand micrometers 34 .To our knowledge however, there is no work regarding the effect of advection for a general channel (symmetric or otherwise), despite many microfluidic devices operating on precisely that scale.An understanding of transport in catalytically active materials is itself a subject of interest, as these materials play a large role in the chemical industry [40][41][42][43][44] , and are expected to be vital for the decarbonization of society 45,46 .Transport on the scales ranging from micrometers to hundreds of micrometers is of great importance to the functioning of catalysts [47][48][49] , and it is at these scales where diffusioosmosis plays a role 28,31 .
One often desires micropumps with flow rates that are stable over time.Thus, the most relevant properties of the active channel are those of its steady states.However, a confined system such as the active channel cannot reach a useful steady state unless some mechanism stops the product concentration from rising.Thus, implementations of the active pore as a micropump must include a sink of solute, regardless of its physical nature.This sink is often imposed in theoretical work by imposing a constant solute concentration in sections of the channel wall.This is unlikely to be the case in experimental setups.One way to ensure the existence of a steady state is to take into account the often-neglected inverse reaction that consumes solute.This reaction is always present 50 , even if the rate may be significantly smaller than the rate of the forward reaction.Including microscopic reversibility has already been shown to alter the dynamics of active colloids 51 .
With this paper, we seek to shed light on the dynamics of active channels/pores by deriving a model of diffusioosmotic transport for corrugated pores.The chemical activity which drives the flows is spatially-varying, and the inverse chemical reaction is taken into account.Both advective and diffusive transport play a role, and our model captures the spontaneous symmetry breaking in fore-aft symmetric pores.We show that when the inverse chemical reaction is taken into account, the control parameter that controls the transition from non-pumping to pumping in symmetric channels is a combination of three timescales: the advective timescale, the diffusive timescale, and the time a solute molecule resides in the pore before being consumed.We further find a second type of advection-enabled transition.This transition occurs only for fore-aft asymmetric channels, which may exhibit bistability, a discontinuous jump in the flow rate, and hysteresis loops.
The remainder of this paper is structured as follows.Our model of the catalytically active pore is presented in Section II.In Section III, the governing equations are systematically reduced to an effective 1D model by employing the Fick-Jacobs and the lubrication approximations.The resulting model is employed in the study of a sinusoidally-shaped pore in Section IV, which is split between fore-aft symmetric pores (Section IV A) and asymmetric pores (Section IV B).We relax the condition of sinusoidal shape and examine pores with a ratcheted shape in Section V. Finally, in Section VI, a summary of the results is given, and a connection to the experimental setup of platinum-coated pores containing a hydrogen peroxide solution is made.

II. THE MODEL
We consider a three-dimensional, axially symmetric pore with its axis along the z direction.The cross-section is thus circular with a radius which varies with z.The ends of the pore are located at z = ±L.The pore walls are located at x 2 + y 2 = R 2 (z), with the average radius R 0 defined as Throughout this work, we will consider R(z) to be a 2Lperiodic function in the z direction.See Fig. 1 for a sketch of the pore geometry.Inside this pore is a fluid that undergoes a chemical reaction due to the presence of catalytic material, which coats the pore walls inhomogeneously.We model this fluid as a mixture of three chemical species, the solvent, the reactant, and the product.In the vicinity of catalytic material, reactant is turned into product.The rate of production depends on factors such as the concentration of the catalytic material, and (in principle) of the reactant.However, in order to keep the model as simple as possible, in the following we specialize to the reaction-limited case.Indeed, this is the case for typical experimental setups 52 .Accordingly, we assume the reactant concentration to be homogeneous in space.The concentration of solvent is then determined by the concentration of product, so as to keep the overall mixture incompressible.
In order to attain a steady state, the reactant must be replenished, and the product must be disposed of.In cases where the fluid is unbound, these processes occur due to diffusion from (for reactant) and into (for product) a chemical reservoir that is placed far away from the catalytic material 53 .In the case where the reaction takes place in a confined space, transfer of reactant (product) into (out of) the reactor has been obtained by exploiting the porosity of the confining walls 54 .Regardless of its physical nature, a mechanism that replenishes reactant and disposes of product is needed to obtain a steady state, and thus steady pumping.We model such a mechanism by accounting for an inverse chemical reaction that converts the product back into reactant.This inverse reaction occurs with a rate χρ(r,t), where ρ is the number density of the product, and χ is a sink constant of dimension s −1 .We further assume that the concentration of reactant is much larger than that of product, and thus does not appreciably change in time.As a consequence, the reaction that generates product occurs with a rate ξ (r) with dimension m −3 s −1 .The spatial dependence of the reaction rate originates for a possibly spatially-varying concentration of catalytic material.Under the assumptions of our model, the influence of the reactant is entirely captured by specifying the function ξ (r).As such, we will no longer mention the reactant concentration.To match common terminology, we will henceforth refer to the product as "solute".The mixture is modelled as a Newtonian, incompressible fluid.Due to the small scales of the pore, the Reynolds number is much smaller than one, and so we describe the fluid flow via the Stokes equation together with the condition for incompressibility where v(r,t) is the velocity profile, η is the dynamic viscosity, and P(r,t) is the pressure.The solute (and therefore, the FIG. 1. Cartoon of the longitudinal cross-section of the three-dimensional, axially-symmetric, active pore.The pore length is 2L and its variable radius is R(z), with an average radius of R 0 .The deviation from this average is δ R(z).The catalysis of the chemical solute occurs near the pore wall with a possibly spatially inhomogeneous rate ξ (r).The deeper the color red, the higher the value of ξ (r).The reverse chemical reaction removes solute in the bulk with a spatially homogeneous rate χ (see inset).The pore walls interact with the solution via the effective potential U wall .A flow field v emerges due to an effective diffusioosmotic slip velocity (blue arrows), which may lead to the onset of a net non-zero flow rate Q.
mixture) interacts with the wall with an effective potential that accounts for the soft interaction with the pore walls, U wall (r), as well as for the confining of the solute inside the pore.Thus, an inhomogeneous concentration of solute leads to an inhomogeneous body force −ρ∇U wall acting on the mixture.At steady state, a pressure gradient develops along the wall, with a corresponding flow that travels up/down the gradient of solute concentration if the potential is repulsive/attractive.This is the phenomenon of selfdiffusioosmosis 23 .Self-diffusioosmosis is operational within a thin layer, of thickness comparable to the potential range r λ , near to the channel walls.Accordingly, for what concerns the velocity profile at positions r ≪ R(z) − r λ far away from the walls, we can approximate the contribution of self-diffusioosmosis as an effective slip velocity located at the channel walls 23 .We will see that only the longitudinal component of this slip velocity will be needed, and this component is given by where r w is the set of points defining the channel wall, and β = 1/(k B T ) is the inverse thermal energy.Furthermore, ∇ || is the derivative along the surface evaluated at the wall, and e z is the unit vector in the z direction.L is the phoretic mobility 23 and is given by where r ′ is the distance to the wall.This expression for L assumes a potential U wall that depends only on the distance to the wall, and a radius of curvature of the wall that is much larger than the potential range 23 .Additionally, the solute concentration is assumed to be in thermal equilibrium within this thin layer where U wall ̸ = 0 23 .From Eq. ( 6), one may see that L may be positive or negative, depending on whether the interaction potential between solute molecules and the pore wall is attractive or repulsive, respectively.Finally, periodic boundary conditions are used in the pore inlet (z = −L) and outlet (z = L).We focus on periodic solutions for two reasons.On the one hand they allow us to keep the model simple.
On the other hand, they allow us to characterize the intrinsic dynamics of the system i.e,.far away from additional physical boundaries along the longitudinal direction.These boundary conditions allow us to recover the spontaneous symmetry breaking observed in the simulations of Ref. [ 34].From Eqs. (2) and ( 5), we see that the steady state flows are entirely determined by the concentration of solute.This concentration evolves in time according to ρ(r,t) = −∇ • j(r,t) + ξ (r) − χρ(r,t) , where j(r,t) is the flux of solute, ξ (r) is a source term, and χρ(r,t) is a sink term.As discussed before, solute production occurs at the pore wall (x 2 + y 2 = R 2 (z)), and so ξ (r) is represented with a Dirac delta where ξ S (z) has dimensions m −2 s −1 , and encodes for the inhomogeneous catalytic coating along the pore's axis of symmetry.It is further assumed to be 2L-periodic.For what concerns the flux of solute, it contains three contributions: one coming from diffusion, one coming from advection due to fluid flow, and one due to the interaction with the pore walls.Thus, the flux j(r,t) is given by j(r,t) = −D∇ρ(r,t) − β Dρ(r,t)∇W (r) + v(r,t)ρ(r,t), (9)  where D is the diffusion coefficient.Owing to the confinement given by the potential (Eq. ( 4)), there is no flux through the pore walls.At the pore inlet (z = −L) and outlet (z = L), the concentration of solute obeys periodic boundary conditions.The governing equations presented above are not analytically solvable.In order to get insight into the physical mechanisms responsible for the onset of the diverse scenarios that we have observed numerically 34 , we aim at a reduced model that provides a deeper understanding of the mechanisms involved in active pumping.It will also enable a much faster and simpler computation of the estimated pumping rates.We proceed by reducing Eqs. ( 2) and ( 7) to effective 1D equations in a procedure analogous to the Fick-Jacobs approach for diffusion in corrugated pores [55][56][57] .A graphical sketch of the following derivation and its assumptions can be found in Appendix C.

III. REDUCTION TO EFFECTIVE 1D EQUATIONS
Equation ( 7) is a 3D partial differential equation whose analytical solution is not easy to find.Moreover, here we are interested in the transport properties of the pore along its axis of symmetry.Therefore we aim at projecting Eq. ( 7) onto this axis (the z direction).The procedure that we aim at is based on a length scale separation between the (long) longitudinal length scale and the (short) transverse length scale.As such, the theory we now derive is valid for narrow pores (R 0 ≪ L).Such an approach (called Fick-Jacobs approximation 55,56,[59][60][61][62][63] ) has been widely used to model the transport of different systems, spanning from colloidal particles [64][65][66][67][68] to polymers 69,70 and including electrolytes [71][72][73][74] and active systems 34,[75][76][77] .

A. Equation for solute concentration
In this section, we derive an effective 1D equation for the solute concentration.This equation governs the time evolution of the solute concentration integrated along the pore cross-section.As such, the dynamics in the transverse degrees of freedom is integrated over and does not need to be explicitly solved for.To do so, we will exploit the length scale separation between the potential range r λ , the average pore radius R 0 , and the pore length L.
As discussed, the range r λ of the interaction potential between the solute and the wall is typically much smaller than the pore dimensions R 0 and L, and so we focus on the regime in which the advective and diffusive contributions vρ and D∇ρ in Eq. ( 9) dominate the term β Dρ(r,t)∇W (r) in Eq. ( 9).Accordingly, Eq. ( 9) reduces to Due to the axial symmetry of the pore, it is convenient to switch to cylindrical coordinates.Let r = (x, y, 0) be a vector running from the axis of symmetry (x, y) = (0, 0) to a generic point inside the pore.This vector is perpendicular to the axis of symmetry, and its magnitude r is the distance from the aforementioned point to the axis of symmetry.We also define φ as the angle formed by r with the xz-plane.Integrating Eq. ( 7) over the cross-section and exploiting the axial symmetry of the system leads to The above step exploits the fact that, due to Eq. ( 4), there is no solute anywhere outside the pore, i.e., and and therefore, As discussed, we will now focus on the case of narrow pores such that which introduces a separation of length scales between the radius and the length of the channel.The result is also a separation between the relaxation timescale in the radial and in the longitudinal direction, with the relaxation of ρ occurring much faster in the former.Thus, we employ the Fick-Jacobs approximation [55][56][57] , which consists in assuming no net transport in the radial direction, i.e.
As a result of this fast relaxation, we may examine the dynamics in the radial direction separately from the longitudinal direction.We begin by estimating the magnitude of the advective transport via the incompressibility equation (Eq.( 3)), i.e., where in the last step we have assumed that the axial symmetry of the pore leads to an axial symmetry of the velocity profile, i.e., ∂ φ v = 0. We now estimate how large the velocities in the radial direction are, when compared to those in the longitudinal direction.We introduce dimensionless variables where v * z and v * r are the order of magnitude of the longitudinal and transverse components of the velocity, respectively.Thus, hatted variables are dimensionless and take values of order unity.Using these definitions, we may write Eq. ( 17) as where K is a quantity of order unity given by Due to the form of Eq. ( 22), K must be a constant.We now introduce the radial Péclet number Pe r by quantifying the relative magnitude of diffusive and advective timescales in the radial direction: and as per Eq. ( 22) In a sufficiently narrow pore (R 0 ≪ L), one has Pe r ≪ 1.This means that diffusive transport dominates the advection due to the fluid flow in the radial direction.
For pores that are narrow enough, diffusion leads to a solute concentration that depends weakly on the radial direction, despite solute being produced only at the pore wall.Furthermore, the flow field is entirely determined by the solute concentration at the wall alone (via the slip boundary condition of Eq. ( 5)).We make use of this property, and neglect the variation of solute concentration in the radial direction: Accordingly, inserting Eq. ( 26) into Eq.( 11) and using Eq. ( 10) leads to with We note that Q is the volumetric fluid flow which, due to incompressibility (see Eq. ( 3)), does not vary along the z direction.For later use, we define the integrated source ξ (z) and integrated solute concentration p(z,t) which, when plugged-in to Eq. ( 27) yields

B. Equation for flow rate
In this section, we handle the Stokes equation to obtain the flow rate Q, as required for the solution of Eq. ( 31).To do so, we once again make use of the length separation between the average pore radius R 0 and the pore length L. Such a separation leads to the velocity in the longitudinal direction varying strongly in the radial direction, as compared with the longitudinal direction.Further requesting that the pore radius varies smoothly ((∂ z R(z)) 2 ≪ 1) enables Q to be written as a functional of the solute concentration.
As discussed, to solve Eq. ( 31), one has to express Q(t) in terms of p(z,t).As mentioned in section III A, the pore is considered to be narrow (L ≫ R 0 ) and axially symmetric, as is the velocity profile.Within such a regime, we exploit the lubrication approximation 78 .To do so, we write Eq. ( 2) for both longitudinal and transverse coordinates and introduce the dimensionless pressure P(r, ẑ,t) = P(r, z,t) P * , where P * is the order of magnitude of the pressure.With this definition and those of Eqs. ( 18) -( 21) , Eq. ( 32) can be written as and since the pore is narrow (R 0 ≪ L), the left hand side of Eq. ( 35) may be approximated as We now write Eq. ( 33) as where we have used Eq.(22).By comparing Eqs. ( 35) and (37), we may write and thus the variation of the pressure along the transverse direction is negligible As a result, we may write Eq. ( 36) as where we have restored the dimensional quantities.Integrating Eq. ( 40) twice along the radial direction leads to where we have used the diffusioosmotic slip boundary conditions v z (r = R, z,t) = v 0 (z,t).Note that ∂ z P(z,t) contains, in addition to a possible external pressure drop, also contributions stemming from fluid incompressibility that will ensure that Q does not depend on z.The volumetric fluid flow is then As discussed, we focus on cases where P(z) fulfills periodic boundary conditions.As such, the integral of ∂ z P(z,t) over the pore length has to vanish, i.e.
This allows one to determine Q(t): Finally, we need to write the slip velocity v 0 (z,t) as a function of p(z,t).The vector perpendicular to the pore wall is: where e r is the unit vector pointing in the radial direction.With this expression Eq. ( 5) turns into (46) We recall that Eq. ( 10) is particularly valid in the regime of weakly varying pore radii.In this regime one has (∂ z R) 2 ≪ 1, and Eq. ( 46) can be approximated as where we have used ∂ r ρ = 0 coming from Eq. ( 26).By using Eqs.( 26) and ( 47), Eq. ( 44) becomes

C. The case of weakly-corrugated pores
The 1D integro-differential equation obtained by plugging Eq. ( 48) in Eq. ( 31) remains challenging to solve analytically.Accordingly, we will expand both equations to linear order in the corrugation height.This truncation is a good approximation for weakly-corrugated pores, and allows for an analytical solution for the steady state.We thus define where where in the limit of low corrugation We also expand p(z,t) From Eqs.( 31) and ( 48) we obtain Using the fact that both R(z) and p(z) are periodic, the leading term in the right-hand-side of Eq. ( 53) is From Eq. ( 54), we find the contribution p 0 (z) as the solution of (56) We will henceforth neglect contributions to p(z,t) of higher order in δ R. As the values p l for l > 0 will no longer be considered, we simplify the notation by dropping the subscript and renaming p 0 (z,t) as P(z,t) P(z,t) ≡ p 0 (z,t). ( We now perform a Fourier expansion in space of P(z,t) and ξ S (z), such that where we define the Fourier coefficients {P j } and {ξ j } associated to the symmetric part of the expansion, and the coefficients { P j } and { ξ j } associated to the antisymmetric part of the expansion.We have also defined k j = (π/L) j.Plugging Eqs. ( 58) and (59) in Eq. ( 56) leads to Ṗ j (t) = a j P j (t) + Qb j P j (t) + ξ j , where We insert Eq. ( 58) into Eq.( 55) to obtain where which depends only on the geometrical properties of the pore.
we obtain In steady-state, Eqs. ( 60) and ( 61) result in which, if a 2 j + b 2 j Q 2 ̸ = 0, is solved by We note that in the non-pumping case Q = 0, even/odd contributions to the solute distribution depend solely on the even/odd contributions of the source function ξ (z).
IV. SINUSOIDAL PORE Equation ( 64) together with Eqs. ( 60) and ( 61) results in a polynomial equation of (in principle) infinite degree for the flow rate Q.While for a general shape R(z) and source ξ (z) finding the roots of such an equation analytically is not possible, the problem becomes much simpler in the case of a sinusoidal pore such that Without loss of generality, we may place the center of our frame of reference where the pore is narrowest (bottleneck), meaning R1 = 0 and R 1 < 0. The larger the value of |R 1 |, the larger the height of the corrugation.Note that our assumption of weak corrugation limits our approach to values of With the definition of the shape R(z) of Eq. ( 74), Eq. ( 64) now yields The fact that the shape R(z) only contains one Fourier mode leads to the conclusion that only the very same mode of the source function contributes to the flow rate.The reverse is also true: should the source function ξ (z) only contain one Fourier mode, only the same Fourier mode of the shape R(z) contributes to Q. Thus, without loss of generality, we may write the source function as where ξ ′ > 0 is the amplitude of the variation of the source, and z 0 is the shift between the maximum of the source and the minimum of R(z).It is useful to write Eq. ( 76) in the form of where A. In-phase pore shape and source (z 0 = 0) We first investigate the case z 0 = 0, where the pore is foreaft symmetric, i.e.R(z) = R(−z) and ξ (z) = ξ (−z).In this case, Eq. ( 75) reduces to where Π(Q) is the polynomial on the right-hand-side of Eq. ( 80).First, we note that both Q and −Q are solutions, as required from the symmetry of the pore.Furthermore, the non-pumping state Q = 0 is always a solution, albeit not always a stable one.A linear stability analysis reveals that Q = 0 is only stable if it is the only real solution (see Appendix D).
There is only a non-zero solution if Π(Q) is such that as seen graphically in Fig. 2 (a).This condition yields which, combined with Eqs. ( 62) and ( 63) results in a minimum corrugation height R c such that pumping only occurs when with the critical value It is important to note the special case of a flat pore (R 1 = 0), for which no pumping is possible, as seen in Fig. 2 (b).The solutions Q ̸ = 0 can be obtained from Eq. ( 80) yielding It must be noted that there can only be pumping This means that we can only obtain pumping in two situations: If the potential is repulsive 79 (L < 0), then the source has to be strongest at the bottleneck (ξ 1 > 0).This result is consistent with previous works on active pores 34 .If the potential is attractive (L > 0), the exact opposite has to occur.Pumping then occurs only when the source is strongest where the pore is widest (ξ 1 < 0).This is analogous to the case of the isotropic diffusiophoretic colloids 35,36 .From here on out, for simplicity and clarity of the results, we will assume L < 0, ξ 1 > 0, but the generalization is straightforward.To understand Eq. (85), we define two timescales.One timescale corresponds to the average time between creation and destruction of a solute molecule, which we call the average lifetime The second timescale τ D characterizes diffusion.More precisely, it is the time that it takes for a solute molecule to diffuse for a length L: Substituting Eqs. ( 87) and (88) in Eq. (85) yields 1.5 1.0 0.5 0.0 0.5 1.0 1.5 0.00 0.02 0.04 0.06 0.08 0.10 The smaller either τ χ or τ D are, the smaller the flow rate will be.Indeed, diffusion acts to inhibit the spontaneous symmetry breaking.The smaller the value of τ D , the stronger the effect of diffusion.Alongside that process, the timescale τ χ limits the length over which a solute molecule may be advected: if a solute molecule does not reside for long enough to reach the bottleneck, the spontaneous symmetry breaking cannot occur either.Thus, we define a new timescale τ a which combines these two effects: Plugging this quantity into Eq.(89) yields We have as of yet not characterized the timescale for advection τ Q , defined as the time that a solute molecule takes to be advected over a length L. Together with τ D and τ χ , τ Q characterizes all processes/terms in Eq. ( 7).Therefore, it is natural to identify τ Q from the terms in Eq. (91) that are not yet accounted for.The shape of Eq. ( 91) suggests the definition such that one may write where we have defined the flow rate scale as which roughly corresponds to the flow rate needed to empty half of the pore volume in the time τ a .Indeed, τ Q is the only term that depends on the phoretic mobility L , and therefore on the flow velocities.Equation (93) describes a master curve for Q that is displayed in Fig. 2 (c).It can be seen from Eq. ( 93) that the appropriate adimensional number that controls the pumping transition is the ratio τ a /τ Q .In the limit of weak inverse chemical reaction τ χ ≫ τ D , the timescale τ a is well approximated by τ D and lim where Pe = τ D /τ Q is the Péclet number.In the limit of negligible inverse chemical reactions, the pumping condition yields a critical Péclet number of π 2 for all sinusoidal pores (with z 0 = 0).
To further clarify the influence of the pore geometry on the pumping rate, let us examine the role of the pore length L independently of all other parameters.Note that both τ Q and τ a depend on L, but with a different functional form.The timescale τ χ is independent of L, the timescale τ D goes with L 2 , and the timescale τ Q goes with L. Thus, the ratio τ a /τ Q is non-monotonic with L, as seen in Fig. 3 (a).This results in two critical values of L such that pumping can only occur if L min < L < L max , where where and which corresponds to the average length a solute molecule diffuses in its lifetime (multiplied by π).Furthermore, the ratio τ a /τ Q (and thus Q/Q a ) attains a maximum when L = L χ , independently of any other pore parameters, as can be seen in Fig. 3 (b).This diffusive length scale is thus the optimum pore length, no matter the amplitude of the corrugation, the average pore length, or the magnitude of the phoretic mobility.
B. Phase-shift between shape and source (z 0 ̸ = 0) We now break the fore-aft symmetry of the system by introducing an offset z 0 ̸ = 0 between the shape R(z) and the source function ξ (z).Note that both are still sinusoidals, but the bottleneck is no longer located where the source is most intense.We obtain the possible steady-state values of the pumping rate Q from Eq. ( 64) as  which has only one real solution if Note that if z 0 = 0, and ξ = 0, then Eq. ( 101) reduces to Eq. ( 83).If z 0 ̸ = 0, it is still possible to obtain three real (only two of them stable) solutions of Q, meaning pumping in both directions is still possible even when the fore-aft symmetry is broken.The fact that ξ1 ̸ = 0 means that Eq. ( 100) cannot be reduced to a quadratic equation, as could Eq. ( 80).As finding the steady state values of Q now requires finding the roots of a 0.0 0.4 0.8 1.2 1.6 2.0 (a) cubic polynomial (and the degree increases when considering a non-sinusoidal pore), from here on out, we rely on numerical tools to obtain Q.The stability of the steady states is obtained from a linear stability analysis (see Appendix D).In Fig. 4 (a), we report the case z 0 = 0 again, as a reference.Slightly shifting the source function (z 0 /L ≈ 10 −2 ), leads to an imperfect bifurcation 80 i.e., a deformation and detachment of the bifurcation branches.A stable branch emerges such that Q > 0 for all values of |R 1 |/R 0 .This flow is directed from the bottleneck towards the maximum of the source.A second stable branch with flow going in the opposite direction (Q < 0) is still present, but the magnitude of its flow rates, as well as the optimum value of |R 1 |/R 0 for which it appears, decreases as compared to the symmetric z 0 = 0 case.Note that the two stable branches no longer meet and that the nonpumping state Q = 0 is no longer a solution.Increasing the offset between source and shape even further (z 0 /L = 10 −1 ) causes this second stable branch to be pushed to larger values of |R 1 |/R 0 that are no longer within the regime of validity of this model.
It needs to be noted that the values of Q are antisymmetric with z 0 , as required from the symmetry of the problem.This is clearly seen in Fig. 4 (b).This figure shows also that there can be a range of z 0 /L for which two stable branches exist, if the corrugation height is large enough.These two branches also do not meet.Indeed, this implies a discontinuous transition and a sudden inversion of flow if, in an experimental realization, one slowly changes the value of z 0 such as to cross from one stable branch to the other.Figure 4 (c) shows how such a procedure leads to a hysteresis loop, with a range of flow rates not accessible to the active pore.
Conversely, if z 0 /L = ±1, there is no pumping as this corresponds to the z 0 = 0 case discussed before, but with ξ 1 < 0 i.e., when the minimum value of the source is located in the bottleneck.To characterize the pumping performance of the channel, we now focus on the highest value of |Q| from the three possible solutions, which we call max(|Q|).The dependence of max(|Q|) on |R 1 |/R 0 and z 0 /L is plotted on Fig. 4  (d).While increasing |R 1 |/R 0 generally increases the flow rate, there is an optimum value of z 0 /L beyond which increasing the shift z 0 between source and shape decreases the flow rate.Finally, we can see how the maximum values of |Q| are maximized for different offsets z 0 /L depending on the value of |R 1 |/R 0 in Fig. 4 (d).
For the symmetric (z 0 = 0) system, the pumping transition is captured entirely by the ratio τ a /τ Q .For the present asymmetric (z 0 ̸ = 0) case, while there is no pumping transition as before, there is still a transition between one and three real solutions of Q.This transition corresponds then to the point at which the pore may pump in both ways rather than just one, and is determined by Eq. ( 101).Plugging in the expressions for the timescales τ a and τ Q shows that they are not enough to capture this new transition.Indeed, one has to introduce a second advective timescale τQ , associated to the magnitude of the anti-symmetric part of the source ξ1 rather than the symmetric part of the source ξ 1 .With this new timescale, one identifies the condition for which the pore may pump in both ways as which is verified in Fig. 5.To shed light on the meaning of Eq. ( 103), let us expand ξ 1 and ξ1 to linear order in z 0 /L.Expanding Eqs. ( 78) and ( 79), one obtains and Eq.(103) becomes In the strong pumping regime where τ Q ≫ τ a , we obtain and where Q sym is the flow rate obtained when z 0 = 0 as per Eq.(93).Thus, Eq. ( 107) can be written as with the right-hand-side being roughly the distance a solute molecule is advected in a time τ A .Thus, pumping in both directions requires the flow rate to be high enough to carry solute from the maximum of the source function to the other side of the bottleneck (a distance z 0 ) before the solute is destroyed or diffuses away.In the case of a symmetric channel where z 0 = 0, Eq. ( 108) reduces to the condition |Q sym | > 0, which is equivalent to Eq. (83).

V. FORE-AFT ASYMMETRIC SHAPE
To further clarify the role of the pore shape on the pumping dynamics, we now study a pore with more complex geometry, described by a function R(z) that is non-sinusoidal and asymmetric.We now define the radius as which has its minimum value at z = L r , and exhibits a form similar to a smoothened sawtooth (see Fig. 6 (a)).The asymmetry is controlled by the length L r .Indeed, if L r = 0, the function R(z) becomes sinusoidal, and thus fore-aft symmetric.In order to assess the role of shape asymmetry, we need to take into account all the possible Fourier modes.Thus, we require a source function that contains all Fourier modes as well.To such effect, we introduce a step function source where Ω sym is the region of points at a distance smaller than L s from the bottleneck, when accounting for periodic boundary conditions.A diagram of the source function can be found in Fig. 6 (a).Both the source and the shape defined above exhibit non-zero Fourier components for all wavenumbers.This means the sum in Eq. ( 64) generally contains an infinite number of terms.We thus truncate the sum to include only the first fifteen terms, and we solve for Q numerically.Including further terms does not significantly change the obtained solutions.We begin by investigating the effect of the extent of chemical activity in Fig. 6 (b).If L s /L = 1, the entire pore produces solute, and no pumping is obtained for any value of L r .This is because a homogeneous source promotes a constant solute density, which results in zero slip velocity, and thus, no flow.This result holds up to linear order in δ R, but may not hold to quadratic order or above.Indeed, Michelin et al. studied homogeneously-coated ratcheted channels, and found a flow rate that grows non-linearly with the corrugation height 81 .
In the symmetric case (L r = 0), pumping occurs only in a range of L s /L centered around 1/2.Indeed, the flow rate exhibits two bifurcations and is symmetric with respect to L s /L = 1/2 (see Fig. 6 (b)).When asymmetry in introduced (L r /L ̸ = 0), the symmetry with respect to L s /L = 1/2 is lost, and the two stable branches are no longer symmetric with respect to Q = 0.The pitchfork bifurcations disappear and are replaced by imperfect bifurcations 80 , as for the sinusoidal pore with z 0 ̸ = 0. Again there is a discontinuous jump with changing L s as the system jumps from the stable branch with Q < 0 to the stable branch with Q > 0. This discontinuity now arises from R(z) being itself asymmetric, rather than from a shift between two symmetric functions R(z) and ξ (z), as was the case in the sinusoidal pore.For large enough asymmetries, however, there is only one solution, with the maximum flow rate depending on the values of L r /L.This can be seen easily in Fig. 6 (c), where the flow rate grows with L r /L until a maximum is reached, after which the flow rate decreases as L r /L → 1 (note that our assumption of weakly-varying R(z) limits our approach to values of L r /L lower than about 0.9).Indeed, pumping is promoted when high gradients in solute density and in the radius coexist in the same position, as seen in Eq. ( 55).If the drop in source strength (and thus, the drop in solute density) is located in a steep region, one observes higher values of Q, hence why an increasing value of L s /L leads to smaller values of the optimal value of L r /L.
In the pores explored until now, the discontinuities in the flow rate manifested only for a prohibitively small degree of asymmetry.We now show that the discontinuities are a general phenomenon in asymmetric phoretic pores, and can be obtained for strongly asymmetric pores as well.To this effect, we consider solute production in only one of the halves of the pore as shown in Fig. 7 (a).
Here, Ω asym is the upper half of Ω sym (accounting for periodic boundary conditions).Note that L s may now be negative (indicating that there is production only in a region in the weaklysloped half of the pore z < L r ).It can be seen in Fig. 7 (b) that the symmetric pore (L r /L = 0) exhibits only one branch for Q, and it is antisymmetric.The maximum value of |Q| is found at L s /L ≈ ±0.8.The antisymmetry is broken when pore asymmetry is introduced (L r /L > 0).The maximum value of Q is located at increasingly lower values of L s /L the more asymmetric the pore is (higher L r /L).All curves show no pumping (Q = 0) when L s = 0.This is because the pore is chemically inert in this case.Remarkably, asymmetry induces a minimum at Q = 0, with positive flow rates for both positive and negative L s .This behavior acquires added complexity for sufficiently asymmetric pores (L r /L = 0.9), where two stable branches can be found in a region of L s .A second maximum arises for L s < 0 and the flow rate there is even higher than the maximum for L s > 0. Therefore, sufficiently asymmetric pores pump optimally when chemical activity is present in the weakly-sloped half (z < L r ).This is in contrast with symmetric / weakly-asymmetric pores which pump optimally when chemical activity is present in the strongly-sloped half (z > L r ). Figure 7 (c) shows how for L s /L = −0.1, it is still possible to obtain a stable non-pumping state (Q = 0), even when the pore shape is asymmetric.Furthermore, it shows how, for the same value of L s /L, the extent of the asymmetry dictates the direction of pumping.Two stable branches appear for L s /L = −0.5, and a hysteresis loop is possible, now by varying the asymmetry L r /L, rather than the active length L s /L.

VI. CONCLUSIONS
We have derived a model to describe self-diffusioosmotic transport of a chemically-reactive solution in an active pore.We have extended the classical self-diffusioosmotic/phoretic models to include the inverse chemical reaction, which con-sumes solute.By assuming a narrow, weakly-corrugated pore, we have reduced the 3D governing equations to an effective 1D system of advection-diffusion equations (Fick-Jacobs approximation [55][56][57] ).The pore walls are coated in a catalytic material which triggers the production of solute.As the catalytic coating is spatially-varying, thus so is the composition of the mixture.By tuning the extent of catalytic coating, as well as the pore's geometrical properties, one may harness diffusioosmosis to pump fluid across the pore.We have shown that this is true even in the case of a fore-aft symmetric pore, as a spontaneous symmetry breaking occurs when advection plays a large enough role.The model is analytically solv-able for a sinusoidal pore, revealing that the appropriate control parameter is the ratio τ a /τ Q of two timescales: the advective timescale τ Q , which determines how quickly a solute molecule is advected over the pore; and a second timescale τ a , which combines both the diffusive timescale, as well as the timescale associated to the consumption of solute.In the limit of negligible inverse chemical reaction, the control parameter τ a /τ Q reduces to the Péclet number.In Appendix A, we show that this control parameter can be used to approximately predict the onset of pumping for non-sinusoidal pores as well.Pores that break fore-aft symmetry, either due to their shape or due to their catalytic coating, pump as well, even if advection plays no role.Furthermore, asymmetric pores can still display bistability, where pumping in both directions is permitted.The bistability is such that the flow rate exhibits discontinuous jumps and hysteresis loops upon tuning the parameters that control the asymmetry.For sinusoidal pores asymmetrically coated in solute, it was found that the ratio τ a /τ Q is not sufficient to describe this second transition between one and two possible directions.Instead, the appropriate control parameter requires the introduction of a fourth timescale, also associated to advective transport.Typically, experimental realizations of diffusioomosis make use of a hydrogen peroxide solution which is decomposed upon contact with platinum.We note that, while the inverse chemical reaction (formation of hydrogen peroxide from water and oxygen) is often not a decisive factor for experimental conditions, oxygen must nonetheless be removed from the pore to achieve a steady state.One possibility is to construct part of the pore's walls out of a porous gel, which allows oxygen to diffuse through 54 .In such a case, the sink constant should be understood as an effective parameter.To obtain analytical results, we approximated the solute concentration as independent of the transverse direction.This approximation is better at smaller scales.Indeed, for an average radius of 1µm, the solute dynamics is expected to be reaction-limited 52 , promoting a solute concentration that is homogeneous in the transverse direction.At larger scales, the Fick-Jacobs theory may still be used, given the appropriate choice of the source function, which should then be taken as an effective parameter.Such possibility arises from the flow field being entirely determined by the solute concentration near the wall.So long as that value is correct, it does not matter what values the solute distribution takes in the remaining pore volume.Flows inside such pores show characteristic velocities of v ≤ 10µm/s 82 .We thus expect the flux associated with pumping active pores to be of magnitude of 10 −2 pL/µm 2 /s.Furthermore, symmetric active pores are expected to exhibit a pumping transition when the Péclet number ≈ 1, and thus when L ≈ 10 2 −10 3 µm 34 .We expect the discontinuous jumps in the flow rates of asymmetric pores to manifest at similar scales.An example of a catalytic material which exhibits pores of such sizes is that of metal foam catalysts 48,49 .We thus predict that advection-enabled instabilities will play a role in the functioning of such catalysts.To experimentally verify the discontinuous jumps and hysteresis in the flow rate, we thus propose experiments with platinum-coated pores containing hydrogen peroxide.These pores should showcase a length in the hundreds of microme-ters.We suggest the fabrication of pores with slight asymmetry in the pore shape, as in Fig. 6 (a).To access both branches of Fig. (6) (c), one may initially apply a pressure drop across the pore, to favor relaxation to either positive or negative values of Q (depending on the sign of the pressure drop).Measuring the values of Q for pores with increasing L r /L, one may identify the end of a stable branch as the value of L r /L beyond which the pore only pumps in one direction, regardless of the pressure drop applied at the start of the experiment.steady state, rather than sustained oscillations in the flow rate.
As explained in Ref. [ 34], said oscillations arise when the fluid relaxation time is finite.In the Fick-Jacobs theory, this is not the case, as we assume Stokes flow (Eq.( 2)).Therefore, the oscillatory regime is beyond the scope of the current theory.Finally, it must be said that the active pore in the simulations is wide (R 0 /L = 1.7), and thus out of the regime in which the approximations that lead to the Fick-Jacobs theory are valid.Nonetheless, we have shown that the theory may still show semi-quantitative agreement beyond this regime.where P 0 , P j and P j are the coefficients in Eqs. ( 72) and (73).We can now write (D15) Note that the matrix M is not diagonal in general.In simple cases, it is possible to obtain the eigenvalues of M analytically, and thus to determine the stability of each steady state.In the most general case, we determine stability by solving the dynamical system of equation (D14) using forward Euler time integration.Starting from a random set of amplitudes at time t = 0, we determine stability by whether the point {q j }(t) diverges from the origin given a certain time.
Consider the case where B j is zero or both Γ j and Θ j are zero for a specific value of j.This ensures that the offdiagonal element is zero.This is the case if the source function ξ (x) and the shape R(x) are both of the form a 1 + a 2 cos(k r x) + a 3 sin(k r x).Then B j , Γ j , Θ j = 0, for all values of j ̸ = r.Under these conditions, the wavenumbers are decoupled, and we obtain a matrix equation for q j : where The real part of M ′ 11 and M ′ 22 is negative and therefore these modes are stable.The only mode that can possibly be unstable is that corresponding to k r , the one that describes the source and shape.We obtain for j = r, where for relaxation to steady state with no ringing (Im(λ ) = 0).

2 0
FIG. 4. (a) The possible steady state flow rates as function of |R 1 |/R c for different values of z 0 /L.Open symbols correspond to unstable steady states and filled-in symbols correspond to stable steady states.Parameters are (L /β η)(τ χ /L 5 ) = −10 −6 , ξ 0 Lτ χ = ξ 1 Lτ χ = 5 × 10 4 , R 0 /L = 0.1, Dτ χ /L 2 = 10 −3 .(b) The possible steady state flow rates as function of z 0 /L for different values of |R 1 |/R c .Parameters as in (a).(c) Zoom-in of the data in panel (b).Arrows mark the discontinuities that give rise to a hysteresis loop.(d) The maximum flow rate one can obtain as a function of z 0 /L and |R 1 |/R 0 .The dotted line marks the value of z 0 /L that maximizes Q/Q a for each value of |R 1 |/R 0 .Parameters as in panel (a).

FIG. 10 .
FIG. 10.Diagram of the derivation of the model, including necessary assumptions.Arrows indicate which assumptions are required to proceed in the derivation.The acronym "P.B. C." stands for "periodic boundary conditions".

TABLE I .
FIG. 9. (a) Classification of the asymptotic dynamics seen from the Lattice Boltzmann simulations of Ref. [ 34] into non-pumping states (×), steady pumping states (•), and oscillating states (•).The line is the Fick-Jacobs theory's prediction for the onset of pumping for α f it = 17.1.(b) Absolute value of state flow rate |Q| as function of Pe ′ as obtained from the Fick-Jacobs theory for α f it as per Table I (×) and Lattice Boltzmann simulations (•, •).Lines are guides to the eye connecting simulation data.Colors indicate value of L Values of α f it for data in Fig. 9 (b).