Directed motion of spheres induced by unbiased driving forces in viscous fluids beyond the Stokes' law regime

PHYSICAL REVIEW E 97, 032219 (2018) Directed motion of spheres induced by unbiased driving forces in viscous fluids beyond the Stokes’ law regime Jesús Casado-Pascual* Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, 41080 Sevilla, Spain (Received 17 September 2017; revised manuscript received 27 February 2018; published 30 March 2018) The emergence of directed motion is investigated in a system consisting of a sphere immersed in a viscous fluid and subjected to time-periodic forces of zero average. The directed motion arises from the combined action of a nonlinear drag force and the applied driving forces, in the absence of any periodic substrate potential. Necessary conditions for the existence of such directed motion are obtained and an analytical expression for the average terminal velocity is derived within the adiabatic approximation. Special attention is paid to the case of two mutually perpendicular forces with sinusoidal time dependence, one with twice the period of the other. It is shown that, although neither of these two forces induces directed motion when acting separately, when added together, the resultant force generates directed motion along the direction of the force with the shortest period. The dependence of the average terminal velocity on the system parameters is analyzed numerically and compared with that obtained using the adiabatic approximation. Among other results, it is found that, for appropriate parameter values, the direction of the average terminal velocity can be reversed by varying the forcing strength. Furthermore, certain aspects of the observed phenomenology are explained by means of symmetry arguments. DOI: 10.1103/PhysRevE.97.032219 I. INTRODUCTION The generation of directed motion of particles in nonlinear systems subjected to deterministic and/or stochastic unbiased driving forces has been an active research topic over the past decades [1–3]. This phenomenon—commonly known as ratchet effect—has been studied extensively in fields as diverse as biophysics [4], nanotechnology [5], granular media [6], and spatially extended nonlinear systems [7,8]. On the theoreticalside,muchworkhasbeenfocusedonunderstanding the mechanisms by which nonlinearity and symmetry breaking cooperate to favor motion in one direction [2,3]. One class of ratchet models that has received considerable attention in the literature is the so-called rocking ratchet [9]. In the rocking ratchet, particles moving in a periodic substrate potential also experience a periodic or quasiperiodic timedependent force of zero average. Although most of the works also consider the presence of random forces, noise is not a crucial element [10,11]. In this class of models, the mechanism behind the generation of directed motion is basically harmonic mixing [3,8,12]. For this mechanism to be effective, certain spatiotemporal symmetries [13], supersymmetries [14], and hidden symmetries [15] must be broken. Alternative ways of breaking these symmetries have also been investigated in a two-state Brownian motor, realized with Brownian particles alternating between two phase-shifted, symmetric potentials [16]. It should be noted that, in the traditional rocking ratchet, the friction forces, if any, are assumed to be linear in the velocities of the particles. Consequently, the appearance of harmonic mixing is solely due to the substrate potential, which is the only source of nonlinearity in the system. *[email protected] However, there are situations in which friction cannot be assumed to be linear in the velocities of the particles. This is the case, for instance, of the drag force—a type of friction that acts on bodies moving in viscous fluids—at Reynolds numbers beyond the Stokes range [17]. Nonlinear friction forces also appear in other fields such as in the active Brownian motion [18], which describes the motion of self-propelled organisms, in the relativistic Brownian motion [19], as well as in the Brownian motion with dry friction [20]. In these situations, friction constitutes an additional source of nonlinearity that may give rise to harmonic mixing, thus playing a major role in the ratchet effect. In particular—and contrary to what intuition might suggest—nonlinear friction may contribute positively to the emergence of directed motion under appropriate conditions. In the case of random forces, this possibility has been analyzedinRef.[21],wherethefrictionconsideredisnonlinear and anisotropic. The primary aim of the present work is to investigate how directed motion emerges from the combined action of nonlinear friction and zero-mean oscillating forces. To this end, we focus our analysis on a simple but realistic model of a physical system in which friction is the only nonlinear element present. More specifically, we examine the motion of a sphere immersed in a viscous fluid and subjected to a time-periodic force of zero average, covering a wide range of Reynolds numbers. We pay special attention to the case of two mutually perpendicular forces with sinusoidal time dependence, one with twice the period of the other. In particular, we show that, although neither of these two forces induces directed motion when acting separately, the resultant force obtained by adding them together causes a net motion of the sphere along the direction of the force with the shortest period [22]. The outline of the remainder of this paper is as follows. In Sec. II, we introduce the system under consideration and 2470-0045/2018/97(3)/032219(8) 032219-1 ©2018 American Physical Society JESÚS CASADO-PASCUAL PHYSICAL REVIEW E 97, 032219 (2018) define the quantities of interest, namely, the time-dependent terminal velocity and the average terminal velocity. We also provide necessary conditions for the appearance of directed motion. The convergence of the solutions of the equation of motion to the time-dependent terminal velocity is analyzed in the Appendix. In Sec. III, we derive an expression for the average terminal velocity within the adiabatic approximation. In Sec. IV, the theoretical results of the previous sections are illustrated by numerical simulations. Finally, in Sec. V,we present conclusions for the main findings of our work. II. PROBLEM FORMULATION We consider the motion of a sphere of mass mand radius r, immersed in a steady fluid of density ρfand viscosity η, and subjected to a time-periodic force F(t)ofperiodTand zero time average, i.e., with T 0dt F(t)/T =0. Since the objective is to study how directed motion emerges from the combined actionofthedragforceandF(t),forexpositionalclaritywewill assume that these are the only forces acting on the sphere. This implies that the system is in a microgravity environment and that the size of the sphere is large enough to neglect the effects of the Brownian forces. The time evolution of the velocity of the sphere relative to the fluid, v(t), is governed by the equation of motion m˙ v(t)=Fd(t)+F(t),(1) where the overdot indicates derivative with respect to time and Fd(t) is the hydrodynamic drag force exerted on the sphere at time t. Hereinafter, we will also assume that the hydrodynamic drag force can be expressed in terms of the steady drag coefficient Cd(Re)[17]as Fd(t)=−πrη 4Re(t)Cd[Re(t)]v(t),(2) where Re(t)=2rρf|v(t)|/η is the Reynolds number of the sphere at time t(throughout this paper, a pair of vertical bars indicatesthemagnitudeoftheenclosedvector).Inotherwords, we will assume that the hydrodynamic drag force at any given instant is what it would be if the sphere were moving uniformly with its instantaneous velocity v(t). This assumption is strictly valid in the low-frequency limit defined by the condition r2ρfω/η 1, where ω=2π/T is the fundamental angular frequency of the driving force (see, e.g., §24 of Ref. [17]). In this limit, the velocity varies so slowly in time that the flow can be regarded as steady at any given instant. In practice, the low-frequency limit can be achieved by sufficiently decreasing the driving frequency and/or the size of the sphere (but keeping it larger than the Brownian size). In addition, it is more easily achieved in fluids with high kinematic viscosities η/ρf. By way of example, in the case of a sphere of radius r= 10−4m moving in water, ethanol, or air at normal temperature and pressure (i.e., 20◦C and 1 atm), the low-frequency limit is valid for ω1.0×102rad/s, ω1.5×102rad/s, and ω1.5×103rad/s, respectively [23]. The use of Eqs. (1) and (2) implies that the effects of the virtual mass force and the Basset history force are negligible [24,25]. In addition, if the density of the sphere, ρs=3m/(4πr3), is less than or of the order of ρf, the inertial term m˙ v(t) is also negligible, since it is proportional to the virtual mass force −2πr3ρf˙ v(t)/3 with proportionality constant −2ρs/ρf. In fact, it can be shown (see Sec. III) that, in dimensionless units, the inertial term is proportional to the dimensionless parameter r2ρsω/η and, consequently, it is negligible in the low-frequency limit r2ρfω/η 1ifρs/ρf is less than or of the order of unity. Thus, for inertial effects to be non-negligible, it is necessary to assume that ρf/ρs1 (heavy-particle limit). According to Eq. (2), the drag force is a nonlinear function of v(t). More specifically, the ratio |Fd(t)|/|v(t)|increases with increasing |v(t)|if, as will be assumed henceforth, |v(t)|< 105η/(rρf). This is so given that ReCd(Re) is an increasing function of Reuntil the onset of the “drag crisis,” which occurs for Re≈2×105[25]. It is only in the limit of vanishing ReynoldsnumberthatReCd(Re)tendsto24andEq.(2)reduces to the linear expression Fd,St(t)=−6πrηv(t) (Stokes’ law) [17]. In practice, the departure from Stokes’ law is already quite significant for Reynolds numbers of the order of unity or evensmaller(see,forinstance,Fig.3.9inRef.[24]).Therefore, for the nonlinear nature of the drag force to become apparent, it is sufficient that |v(t)|η/(2rρf); this is precisely the regime of interest in this work. For example, in the case considered above, the nonlinearity becomes significant when |v(t)| 5.0×10−3m/s(water),|v(t)|7.6×10−3m/s (ethanol), and |v(t)|7.6×10−2m/s (air)—velocities that can be reasonably achieved in the laboratory. Given an initial condition for the velocity at some initial instant t0, the value of the velocity at any later time t>t 0 can be calculated by integrating the equation of motion (1). In the Appendix it is shown that, as happens in the linear case, as the time interval t−t0increases, the solutions of Eq. (1) become independent of the initial conditions and converge exponentially to a single time-dependent terminal velocity, which will be denoted by V(t). Moreover, it is also shown that the relaxation time to reach this terminal velocity is less than, or of the same order as, the characteristic timescale τ=m/(6πrη). The time-dependent terminal velocity V(t) is uniquely determined by the equation of motion (1). Consequently, if Eq. (1) is invariant under some transformation, so will be V(t). For example, using Eq. (2) and the fact that F(t)is periodic, it is easy to verify that Eq. (1) is invariant under the transformation v(t)→ v(t)=v(t+T), in the sense that it remains unchanged if v(t) is replaced by v(t)=v(t+T). It then follows that V(t) must also be invariant under the same transformation, i.e., V(t)=V(t)=V(t+T). Therefore, the time-dependent terminal velocity is periodic in time with the same period Tas F(t). We are interested in studying the average terminal velocity V=1 TT 0 dt V(t)(3) and, more specifically, the conditions for Vto be nonzero. Firstly, it should be pointed out that the existence of nonzero average terminal velocities is an unequivocal signal that Stokes’ law is no longer applicable. Indeed, if Stokes’ law wereapplicable, then the time-dependent terminal velocity 032219-2 DIRECTED MOTION OF SPHERES INDUCED BY … PHYSICAL REVIEW E 97, 032219 (2018) would satisfy the linear differential equation ˙ V(t)=−V(t) τ+F(t) m.(4) By integrating the above equation from 0 to T, and taking into account the periodicity of V(t) and that F(t) has zero time average, it would necessarily follow that V=0. In addition to the violation of Stokes’ law, there is another necessary condition for the existence of nonzero average terminal velocities. Let us assume that the force F(t) fulfills the time-shift symmetry F(t)=−F(t+T/2).(5) In this case, it can be easily verified that the equation of motion(1)isinvariantunderthetransformationv(t)→ v(t)= −v(t+T/2), so V(t+T/2) =−V(t). This last result, together with Eq. (3), implies that V=0. Consequently, nonzero average terminal velocities are only possible if the time-shift symmetry (5) is broken. As a particularly illustrative example, let us consider the force F(t)=F1(t)e1+F2(t)e2, where F1(t) and F2(t)are periodic functions with periods T1=Tand T2=T/2, respectively, and e1and e2are two mutually perpendicular unit vectors [22]. The period of F(t) is, thus, equal to T. Let us assume that both F1(t) and F2(t) satisfy the aforementioned time-shift symmetry, i.e., Fj(t+Tj/2) =−Fj(t)forj=1,2. Consequently, if each component of F(t) were considered separately, the resulting average terminal velocity would be zero. By contrast, if both components are added together, a nonzero value of Vis possible given that F(t+T/2) = −F1(t)e1+F2(t)e2=−F(t). Moreover, in this case, the average terminal velocity is necessarily parallel to e2. Indeed, for the considered F(t), the equation of motion (1)isclearly invariant under the transformation v(t)→ v(t)=v1(t)e1+ v2(t)e2−v3(t)e3, where e3is a unit vector perpendicular to both e1and e2. As a consequence, V3(t)=−V3(t)=0. Since Eq. (1) is also invariant under the transformation v(t)→ v(t)=−v1(t+T/2)e1+v2(t+T/2)e2+v3(t+T/2)e3,it then follows that V(t+T/2) =−V1(t)e1+V2(t)e2.Using this last expression in Eq. (3), we obtain that, if V= 0, then it is necessarily parallel to e2[22]. In addition, we have also proved that the period of V2(t)isT/2. III. ADIABATIC LIMIT An explicit expression for the average terminal velocity can be obtained in the adiabatic limit ωτ 1. Notice that, in the heavy-particle limit, the adiabatic limit is more restrictive than the low-frequency limit mentioned in the previous section, since r2ρfω/η =9ρfωτ/(2ρs)ωτ if ρf/ρs1. To study the adiabatic limit, we first introduce the dimensionless quantities θ=ωt,ν(θ)=2rρfv(θ/ω)/η, and f(θ)=F(θ/ω)/F0, with F0being a typical value of |F(t)|such that |f(θ)|⩽1for all θ. The value of the Reynolds number at the dimensionless time θis thus given by |ν(θ)|. In dimensionless variables, the equation of motion (1) becomes ωτ dν(θ) dθ =− 1 24 Cd[|ν(θ)|]|ν(θ)|ν(θ)+f0f(θ),(6) where f0=ρfF0/(3πη2) is a dimensionless parameter characterizing the strength of the driving force. The adiabatic limit of ν(θ), denoted hereafter as νad(θ), can be found by solving the equation Cd[|νad(θ)|]|νad(θ)|νad(θ)=24f0f(θ),(7) which is obtained by taking the limit ωτ →0inEq.(6). Once Eq. (7) is solved and νad(θ) is known, we can immediately determine the adiabatic limit of the average terminal velocity, Vad, using the expression Vad =η 4πrρf2π 0 dθ νad(θ),(8) which is just Eq. (3) rewritten in terms of dimensionless quantities. In order to solve Eq. (7) and obtain νad(θ), an explicit expression for the steady drag coefficient Cd(Re) is required. For the present purposes, we will use the semiempirical expression Cd(Re)=24 Re1+√Re δ02 ,(9) with δ0=9.06 (see Ref. [26] for a heuristic derivation). This expression is in remarkable agreement with the experiments for Re5×103[26]. Consequently, in the case under consideration, its use is justified provided that |νad(θ)|5×103 for all θ. After inserting Eq. (9) into Eq. (7), one obtains 1+√|νad(θ)| δ02 νad(θ)=f0f(θ).(10) It follows from the above equation that the vector νad(θ) points in the same direction as the unit vector f(θ)/|f(θ)|. Furthermore,italsofollowsthatthemagnitudeofνad(θ)isasolutionofthe algebraic equation [1 +|νad(θ)|1/2/δ0]2|νad(θ)|= f0|f(θ)|. Thus, taking into account that |νad(θ)|must be real and non-negative, it is straightforward to show that the only physically meaningful solution of Eq. (10)is νad(θ)=δ2 0 4⎡ ⎣−1+1+4√f0|f(θ)| δ0⎤ ⎦ 2 f(θ) |f(θ)|.(11) This last expression also allows us to estimate the maximum value of f0consistent with the condition |νad(θ)|5×103. Indeed, since |f(θ)|⩽1 for all θ, it is then clear that |νad(θ)|⩽ δ2 0[−1+(1 +4f1/2 0/δ0)1/2]2/4. From this inequality it can be readily verified that the condition |νad(θ)|5×103is automatically fulfilled for all θif f03.87 ×105. For larger f0values the use of Eq. (9) is no longer justified, and an alternative expression for the steady drag coefficient must be used (see, e.g., Ref. [24] for a discussion of possible choices). Finally, according to Eq. (8), the adiabatic limit of the average terminal velocity, Vad, can be calculated by evaluating the integral Vad =δ2 0η 16πrρf2π 0 dθ⎡ ⎣−1+1+4√f0|f(θ)| δ0⎤ ⎦ 2 f(θ) |f(θ)|. (12) 032219-3 JESÚS CASADO-PASCUAL PHYSICAL REVIEW E 97, 032219 (2018) As shown in the following section, even though 2π 0dθ f(θ)= 0,thevalueofVad obtained from the above expression is in general nonzero. IV. RESULTS In this section, we illustrate our results using the particular case of a dimensionless biharmonic force of the form f(θ)=ζcos (θ)e1+(1 −ζ) cos (2θ+ϕ)e2,(13) where ζis a parameter that allows us to simultaneously vary the amplitudes of the two harmonic components, and ϕis the phase difference between them. The parameter ζtakes values between 0 and 1, with ζ=0 and ζ=1 corresponding, respectively, to monochromatic forces along the directions of the unit vectors e2and e1. In order to determine the time-dependent terminal velocity, we have numerically integrated the differential equation obtained by substituting Eqs. (9) and (13) into Eq. (6), with the initial condition ν(θ0)=0. The parameter θ0has been chosen to be negative and much larger in magnitude than the dimensionless relaxation time ωτ, so as to ensure that the asymptotic time-periodic regime has been reached for θ⩾ 0. Once the time-dependent terminal velocity is known, the averageterminalvelocity can be easily computed by evaluating numerically the integral appearing in Eq. (3). In Fig. 1, the method described above has been used to determine the dependence of the dimensionless terminal velocity 2rρfV/η on the dimensionless time θfor three values of ωτ, namely, ωτ =0.1 (dotted lines), ωτ =1.1 (dashed lines), and ωτ =2.1 (dot-dashed lines). Only the first two components of 2rρfV/η have been plotted, since the third one is identically zero. The results obtained by using the adiabatic expression in Eq. (11) are indicated with solid lines. The values of the remaining parameters are f0=100, ζ=0.5, and ϕ=π.We have chosen a relatively large value of f0in order to highlight the effect of nonlinearity. However, this does not necessarily mean that the magnitude of the driving force is also large. For -20 -10 0 10 20 0 -20 -10 0 10 20 θ π/2π3π/22π 2rρfV2/η2rρfV1/η Adiabatic limit ωτ =0.1 ωτ =1.1 ωτ =2.1 FIG. 1. Dependence of the components of the dimensionless time-dependent terminal velocity 2rρfV/η on the dimensionless time θfor ωτ =0.1, 1.1, and 2.1. The results obtained by using the adiabatic expression in Eq. (11) are depicted with solid lines. The remaining parameter values are f0=100, ζ=0.5, and ϕ=π. example, if the fluid is air at normal temperature and pressure, it can be easily verified from the definition of f0that the value f0=100 corresponds to a value of F0of approximately 2.6×10−7N. In the case of a solid iron sphere of radius r=10−4m, this is about 80.4% of its weight on Earth. The results in Fig. 1confirm that the first component of 2rρfV/η (bottom panel) reverses sign every half period and that the second component (top panel) has half the period of the first, as was shown at the end of Sec. II. Furthermore, it is observed that the agreement between the analytical results obtained from Eq. (11) and the numerical results is quite good for the lowest dimensionless frequency, ωτ =0.1, but becomes progressively worse as ωτ increases, as is to be expected in an adiabatic approximation. More specifically, with increasing the value of ωτ, the curves shift to the right and the amplitudes of the oscillations decrease. At this point, it is convenient to recall that, as pointed out in Sec. III, the adiabatic limit is more restrictive than the low-frequency limit if ρf/ρs1; as a result, in the heavyparticle limit, the low-frequency range extends beyond the adiabaticregime.Toillustratethisfact,letusconsideragainthe example of a solid iron sphere of radius r=10−4m immersed in air at normal temperature and pressure. In this case, it is easy to verify that ρf/ρs≈1.5×10−4and r2ρfω/η = 9ρfωτ/(2ρs)≈6.9×10−4ωτ. From this last expression, it is clear that the three values of ωτ considered in Fig. 1are in the low-frequency range. However, as can be seen in Fig. 1,the adiabatic approximation is no longer valid for ωτ =1.1 and ωτ =2.1. Figure 2depicts the dependence of the second component of 2rρfV/η on the relative phase ϕfor the same values of ωτ, f0, and ζas in Fig. 1. In this and the following figures, the first component of 2rρfV/η is not shown as it is identically zero (see the proof in Sec. II). We have restricted the values of ϕto the interval [0,2π] since the average terminal velocity is 2πperiodic in ϕ. This periodicity follows from the invariance of the dynamics under the transformation ϕ→ ϕ+2π[see Eqs. (6) and (13)]. Again, the adiabatic approximation (solid 0 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 π/2π3π/22π 2rρfV2/η ϕ Adiabatic limit ωτ =0.1 ωτ =1.1 ωτ =2.1 FIG. 2. Dependence of the second component of the dimensionless average terminal velocity 2rρfV/η on the relative phase ϕfor ωτ =0.1, 1.1, and 2.1. The result obtained in the adiabatic limit by using Eq. (12) is depicted with a solid line. The remaining parameter values are f0=100 and ζ=0.5. 032219-4 DIRECTED MOTION OF SPHERES INDUCED BY … PHYSICAL REVIEW E 97, 032219 (2018) line) provides a very satisfactory description of the numerical results for the lowest dimensionless frequency (dotted line). A close examination of Fig. 2reveals that the function V2(ϕ), as well as its adiabatic limit Vad,2(ϕ), satisfies the phase-shift symmetry V2(ϕ+π)=− V2(ϕ).(14) To elucidate the origin of this symmetry, observe that replacing ϕby ϕ+πin Eq. (13) is equivalent to changing the sign of the second component of f(θ). Using this fact, it is easy to verify that the equation of motion (1) is invariant under the transformationv(t,ϕ)→ v(t,ϕ)=v1(t,ϕ +π)e1−v2(t,ϕ + π)e2+v3(t,ϕ +π)e3, where the dependence of the velocity on the phase difference ϕhas been explicitly indicated. Since the time-dependent terminal velocity is uniquely determined by Eq. (1), it is clear that V(t,ϕ)=V1(t,ϕ +π)e1−V2(t,ϕ + π)e2. This, together with Eq. (3), leads to Eq. (14). Note that, according to this reasoning, Eq. (14) is valid independently of the value of ωτ and, in particular, in the limit ωτ →0 [i.e., when V2(ϕ) is replaced by Vad,2(ϕ)]. The results in Fig. 2also reveal that the adiabatic limit of V2(ϕ) satisfies the symmetry relation Vad,2(ϕ)=Vad,2(−ϕ),(15) whereas outside the adiabatic regime, generally V2(ϕ)= V2(−ϕ).Tounderstandwhythisisso,observethatthefunction f(θ)inEq.(13) is unchanged if the sign of both θand ϕis switched. From this, it readily follows that Eq. (7) is invariant under the time-and-phase-reversal transformation νad(θ,ϕ)→ ν ad(θ,ϕ)=νad(−θ, −ϕ) and, consequently, that νad(θ,ϕ)= νad(−θ, −ϕ). After an appropriate change of variables, this lastexpression,togetherwith Eq. (8),leads to Eq.(15). Outside the adiabatic regime, however, this argument fails because the invariance of Eq. (6) under time-and-phase reversal is broken by the presence of the first-order time derivative. Thesymmetryrelation(15) is but a special case of the more general relation Vad,2(nπ/2+ϕ)=(−1)nVad,2(nπ/2−ϕ),(16) valid for any integer n.ToproveEq.(16), it suffices to note that, according to Eq. (14), the right-hand side of the equality Vad,2(nπ/2+ϕ)=Vad,2(−nπ/2−ϕ) is equal to (−1)nVad,2(nπ/2−ϕ). By setting ϕ=0inEq.(16), we conclude that Vad,2(nπ/2)=0 (17) if nis odd. The properties (16) and (17) are visible in Fig. 2. As a final comment on Fig. 2, it is worth mentioning that the curves shown can be very well fitted by an expression of the form A1cos(ϕ+χ1)+A3cos(3ϕ+χ3), where A1,χ1,A3, and χ3are fitting parameters which depend on the value of ωτ; thefittingcurvesarenot showninthefigure,astheyarevisually indistinguishable from the original ones. The functional form of this fitting function is a direct consequence of the system symmetries, being independent of the details of the dynamics (see, e.g., Refs. [27,28]). In particular, in the adiabatic limit, it is only necessary to calculate two fitting parameters since, according to Eq. (15), χ1and χ3can be chosen to be zero. In Fig. 3, the dependence of 2rρfV2/η on the parameter ζ is shown for the same values of ωτ,f0, and ϕas in Fig. 1. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 ζ 2rρfV2/η Adiabatic limit ωτ =0.1 ωτ =1.1 ωτ =2.1 FIG. 3. Dependence of the second component of the dimensionless average terminal velocity 2rρfV/η on the parameter ζfor ωτ =0.1, 1.1, and 2.1. The result obtained in the adiabatic limit by using Eq. (12) is depicted with a solid line. The remaining parameter values are f0=100 and ϕ=π. The numerical results obtained for ωτ =0.1 (dotted line) are indistinguishable from those provided by the adiabatic expression (12) (solid line). Notice that, independently of the value of ωτ, the average terminal velocity vanishes for ζ=0 and ζ=1. This is so because the dimensionless force in Eq. (13) satisfies the time-shift symmetries f(θ+π/2) =−f(θ), if ζ= 0, and f(θ+π)=−f(θ), if ζ=1 (see discussion at the end of Sec. II). The curves in Fig. 3also reveal that, for fixed values of the other parameters, there exists an optimal value of ζwhich maximizes the second component of the average terminal velocity. Furthermore, as ωτ increases, the maximum velocity decreases and its location shifts toward lower values of ζ. It should be noted here that, in the lowest order, the general formalism developed in Refs. [27,28] leads to the approximate expression V2(ζ)≈Cζ2(1 −ζ), where Cis independent of ζ. This expression vanishes at ζ=0 and ζ=1, and displays a maximum at ζ=2/3, thus qualitatively resembling the behavior seen in Fig. 3. However, it is unable to account for the dependence of the location of the maximum velocity on ωτ. This deficiency is not surprising, given that the above approximation is expected to be accurate only for small values of f0and, in Fig. 3, we have taken f0=100. In Fig. 4, we plot the dimensionless average terminal velocity 2rρfV2/η as a function of the dimensionless driving strength f0for the same values of ωτ,ϕ, and ζas in Fig. 1.We have limited the range of f0to values well below 3.87 ×105, so as to ensure the applicability of Eq. (9) (see Sec. III). In particular, the large panel shows the results obtained for values of f0varying from 0 to 104, while the inset zooms in on the range 0 ⩽f0⩽400. Again the numerical results obtained for ωτ =0.1 (dotted line) are indistinguishable from those provided by the adiabatic expression (12) (solid line). A glance at the large panel might tempt one to conclude that, for the parameter values considered, an increase in ωτ causes an increase in 2rρfV2/η. However, this conclusion is refuted by the data shown in the inset. 032219-5 JESÚS CASADO-PASCUAL PHYSICAL REVIEW E 97, 032219 (2018) 0 2000 4000 6000 8000 10000 0 5 10 15 20 25 30 0 100 200 300 400 0 0.5 1 1.5 2 f0 f0 2rρfV2/η 2rρfV2/η Adiabatic limit ωτ =0.1 ωτ =1.1 ωτ =2.1 FIG. 4. Dependence of the second component of the dimensionless average terminal velocity 2rρfV/η on the parameter f0for ωτ =0.1, 1.1, and 2.1. The result obtained in the adiabatic limit by using Eq. (12) is depicted with solid lines. The remaining parameter values are ζ=0.5andϕ=π. The inset shows a zoomed-in view of the curves in the range 0 ⩽f0⩽400. In the case shown in Fig. 4the average terminal velocity is an increasing function of f0.FromEq.(14) it is evident that, if we had used ϕ=0 instead of ϕ=π, we would have observed that the average terminal velocity decreases monotonically with f0. The question then arises as to whether there are parameter values for which the average terminal velocity exhibits a nonmonotonic dependence on f0.The answer to this question is affirmative, as can be seen in Fig. 5 for ϕ=π/2 and ωτ =1.1 and ωτ =2.1. Notice that the nonmonotonic behavior is accompanied by the appearance of a current reversal as a function of the dimensionless driving strength f0. As the value of the dimensionless frequency ωτ decreases, the position of the current reversal shifts toward f0=0, becoming indistinguishable from zero for ωτ =0.1. This type of current reversal is quite common in underdamped rocking ratchets (see, e.g., Refs. [10,29]). The above results clearly show that, when nonlinear friction is present, current 0 2000 4000 6000 8000 10000 -5 -4 -3 -2 -1 0 0 100 200 300 400 -0.2 0 0.2 0.4 f0 f0 2rρfV2/η 2rρfV2/η Adiabatic limit ωτ =0.1 ωτ =1.1 ωτ =2.1 FIG. 5. The same as in Fig. 4but now ϕ=π/2. reversals can be observed even in the absence of any periodic substrate potential. V. CONCLUSIONS In this paper, a theoretical study of the motion of a sphere immersed in a viscous fluid and subjected to a time-periodic force of zero average has been presented. Our focus has been on situations in which Stokes’ law is not applicable, so the drag force depends nonlinearly on the velocity of the sphere relative to the fluid. Let us summarize the main results of this work. (i) It has been shown that, when the time-shift symmetry (5) is broken, the combined action of the zero-mean oscillating force and the nonlinear drag force is able to induce a directed motion of the sphere. Unlike in the traditional rocking ratchet, in this case the directed motion emerges in the absence of any periodic substrate potential. (ii) Explicit expressions for the terminal velocity and the average terminal velocity have been derived within the adiabatic approximation. A comparison between the predictions of these expressionsandtheresultsobtainedbynumericallysolving the equation of motion has been carried out. As expected, it has been found that the lower the frequency of the driving force, the more accurate the adiabatic approximation becomes. (iii) By way of example, the case of two mutually perpendicular forces with sinusoidal time dependence has been considered. Although neither of these two forces induces directed motion when acting separately, it has been shown that the resultant force obtained by adding them together causes a net motion of the sphere along the direction of the force with the shortest period [22]. (iv) A detailed analysis of the dependence of the average terminal velocity on the system parameters has been made and some aspects of the observed phenomenology, such as the suppressionoftransportforparticular valuesofthe parameters, have been rationalized using symmetry arguments. (v) A remarkable finding of this analysis is that, for appropriate parameter values, the average terminal velocity exhibits a nonmonotonic behavior as a function of the forcing strength, resulting in the appearance of current reversal. This kind of behavior resembles that observed in underdamped rocking ratchets [10,29]. A natural extension of this work would be to empirically verify the theoretical results reported here. It is hoped that the present paper will provide the stimulus to do experimental research in this area. ACKNOWLEDGMENTS I thank the Junta de Andalucía for funding support. I also want to thank M. L. Olivera-Atencio, N. R. Quintero, R. Alvarez-Nodarse, J. A. Cuesta, and A. T. Pérez for their critical reading of the manuscript and fruitful discussions. APPENDIX: CONVERGENCE OF THE SOLUTIONS OF EQ. (1) TO THE TIME-DEPENDENT TERMINAL VELOCITY Let v(t) and v(t) be two solutions of the equation of motion (1) corresponding to the initial conditions v(t0)=v0 032219-6 DIRECTED MOTION OF SPHERES INDUCED BY … PHYSICAL REVIEW E 97, 032219 (2018) and v(t0)=v 0= v0. By using Eqs. (1) and (2), it is not hard to verify that d dt |v(t)−v(t)|2=−(t)|v(t)−v(t)|2−(t),(A1) where (t)=πrη 4m{Re(t)Cd[Re(t)] +R e(t)Cd[R e(t)]}(A2) and (t)=3πη3 2mrρ2 f{[Re(t)]2−[R e(t)]2} ×{Re(t)Cd[Re(t)] −R e(t)Cd[R e(t)]},(A3) with Re(t)=2rρf|v(t)|/η and R e(t)=2rρf|v(t)|/η. As mentioned in Sec. II, in the range of Reynolds numbers considered in this work, ReCd(Re) is an increasing function of Re. Thus, from Eq. (A3) it follows that (t)⩾0. In addition, since limRe→0ReCd(Re)=24, it is clear from Eq. (A2) that (t)⩾12πrη/m. Using these results, Eq. (A1) leads to the inequality d dt |v(t)−v(t)|2⩽−2 τ|v(t)−v(t)|2,(A4) with τ=m/(6πrη). It then follows from Gronwall’s inequality (see, e.g., Ref. [30]) that |v(t)−v(t)|2⩽|v0−v 0|2e−2(t−t0)/τ .(A5) Therefore, as the time interval t−t0increases, the solutions of Eq. (1) become independent of the initial conditions and converge exponentially to a single time-dependent terminal velocity. Furthermore, according to Eq. (A5), the relaxation time to reach this terminal velocity is less than, or of the same order as, the characteristic timescale τ=m/(6πrη). [1] R. D. Astumian, Thermodynamics and kinetics of a Brownian motor, Science 276,917 (1997); R. D. Astumian and P. Hänggi, Brownian motors, Phys. Today 55(11), 33 (2002); P. Hänggi, F. Marchesoni,andF.Nori, Brownianmotors,Ann.Phys.(Leipzig) 14,51 (2005); D. Cubero and F. Renzoni, Brownian ratchets: From statistical physics to bio and nano-motors (Cambridge University Press, Cambridge, 2016). [2] P. Reimann, Brownian motors: Noisy transport far from equilibrium, Phys. Rep. 361,57 (2002). [3] P. Hänggi and F. Marchesoni, Artificial Brownian motors: Controlling transport on the nanoscale, Rev. Mod. Phys. 81,387 (2009). [4] L. Mahadevan and P. Matsudaira, Motility powered by supramolecular springs and ratchets, Science 288,95 (2000); R. Di Leonardo, L. Angelani, D. Dell’Arciprete, G. Ruocco, V. Iebba,S.Schippa,M.P.Conte, F. Mecarini,F. DeAngelis, andE. DiFabrizio,Bacterialratchetmotors,Proc.Natl. Acad.Sci. USA 107,9541 (2010); D. Chowdhury, Stochastic mechano-chemical kinetics of molecular motors: A multidisciplinary enterprise from a physicist’s perspective, Phys. Rep. 529,1(2013). [5] A. van Oudenaarden and S. G. Boxer, Brownian ratchets: Molecular separations in lipid bilayers supported on patterned arrays, Science 285,1046 (1999); H. Linke, T. E. Humphrey, A. Lofgren, A. O. Sushkov, R. Newbury, R. P. Taylor, and P. Omling, Experimental tunneling ratchets, ibid. 286,2314 (1999); V. Serreli, C. F. Lee, E. R. Kay, and D. A. Leigh, A molecular information ratchet, Nature (London) 445,523 (2007); C. Drexler, S. A. Tarasenko, P. Olbrich, J. Karch, M. Hirmer, F. Müller, M. Gmitra, J. Fabian, R. Yakimova, S. Lara-Avila, S. Kubatkin, M. Wang, R. Vajtai, P. M. Ajayan, J. Kono, and S. D. Ganichev, Magnetic quantum ratchet effect in graphene, Nat. Nanotechnol. 8,104 (2013). [6] D. C. Rapaport, The wonderful world of granular ratchets, Comput. Phys. Commun. 147,141 (2002); G. Costantini, U. Marini Bettolo Marconi, and A. Puglisi, Granular Brownian ratchet model, Phys.Rev.E75,061124 (2007); V. Berardi, J. Lydon, P. G. Kevrekidis, C. Daraio, and R. Carretero-González, Directed ratchet transport in granular chains, ibid. 88,052202 (2013). [7] F.Marchesoni,ThermalRatchetsin1+1Dimensions,Phys.Rev. Lett. 77,2364 (1996); M. Salerno and N. R. Quintero, Soliton ratchets, Phys. Rev. E 65,025602 (2002); L. Morales-Molina, F. G. Mertens, and A. Sánchez, Soliton ratchets out of pointlike inhomogeneities, Eur.Phys.J.B37,79 (2004); Ratchet behaviorinnonlinearKlein-Gordonsystemswithpointlikeinhomogeneities, Phys. Rev. E 72,016612 (2005); B. Sánchez-Rey, J. Casado-Pascual, and N. R. Quintero, Kink ratchet induced by a time-dependent symmetric field potential, ibid. 94,012221 (2016). [8] M. Salerno and Y. Zolotaryuk, Soliton ratchetlike dynamics by acforces withharmonic mixing,Phys.Rev.E65,056603(2002). [9] R. Bartussek, P. Hänggi, and J. G. Kissner, Periodically rocked thermal ratchets, Europhys. Lett. 28,459 (1994). [10] P. Jung, J. G. Kissner, and P. Hänggi, Regular and Chaotic Transport in Asymmetric Periodic Potentials: Inertia Ratchets, Phys.Rev.Lett.76,3436 (1996);J. L.Mateos,ChaoticTransport and Current Reversal in Deterministic Ratchets, ibid. 84,258 (2000). [11] A. Sarmiento and H. Larralde, Deterministic transport in ratchets, Phys. Rev. E 59,4878 (1999); D. Cubero, J. CasadoPascual, A. Alvarez, M Morillo, and P. Hänggi, Overdamped deterministic ratchets driven by multifrequency forces, Acta Phys. Pol. B 37, 1467 (2006); J. Casado-Pascual, D. Cubero, and F. Renzoni, Universal asymptotic behavior in nonlinear systems driven by a two-frequency forcing, Phys. Rev. E 88, 062919 (2013); D. Cubero, J. Casado-Pascual, and F. Renzoni, Irrationality and Quasiperiodicity in Driven Nonlinear Systems, Phys. Rev. Lett. 112,174102 (2014). [12] F. Marchesoni, Harmonic mixing signal: Doubly dithered ring laser gyroscope, Phys. Lett. A 119,221 (1986). [13] S. Flach, O. Yevtushenko, and Y. Zolotaryuk, Directed Current due to Broken Time-Space Symmetry, Phys. Rev. Lett. 84, 2358 (2000); S. Denisov, S. Flach, A. A. Ovchinnikov, O. Yevtushenko, and Y. Zolotaryuk, Broken space-time symmetries and mechanisms of rectification of ac fields by nonlinear (non)adiabatic response, Phys.Rev.E66,041104 (2002). [14] P. Reimann, Supersymmetric Ratchets, Phys.Rev.Lett.86,4992 (2001). 032219-7 JESÚS CASADO-PASCUAL PHYSICAL REVIEW E 97, 032219 (2018) [15] D. Cubero and F. Renzoni, Hidden Symmetries, Instabilities, and Current Suppression in Brownian Ratchets, Phys.Rev.Lett. 116,010602 (2016). [16] H. Hagman, M. Zelan, and C. M. Dion, Breaking the symmetry of a Brownian motor with symmetric potentials, J. Phys. A: Math. Theor. 44,155002 (2011). [17] L. D. Landau and E. M. Lifshitz, Course of theoretical physics. Fluid mechanics, 2nd ed. (Pergamon, Oxford, 2003). [18] F. Schweitzer, W. Ebeling, and B. Tilch, Complex Motion of Brownian Particles with Energy Depots, Phys. Rev. Lett. 80, 5044 (1998); A. Fiasconaro, W. Ebeling, and E. GudowskaNowak, Active Brownian motion models and applications to ratchets, Eur. Phys. J. B 65,403 (2008); P. Romanczuk and L. Schimansky-Geier, Brownian Motion with Active Fluctuations, Phys. Rev. Lett. 106,230601 (2011). [19] J.Dunkeland P. Hänggi,Theoryofrelativistic Brownianmotion: The (1+1)-dimensional case, Phys. Rev. E 71,016124 (2005); Theory of relativistic Brownian motion: The (1+3)-dimensional case, 72,036106 (2005); Relativistic Brownian motion, Phys. Rep. 471,1(2009). [20] P. G.de Gennes,Brownianmotionwithdryfriction,J.Stat. Phys. 119,953 (2005); A. Pototsky and F. Marchesoni, Periodically driven Brownian motion with dry friction and ultrarelativistic Langevin equations, Phys.Rev.E87,032132 (2013). [21] Y. L. Bolotin, A. V. Tur, and V. V. Yanovsky, Nonlinear friction as a mechanism of directed motion, Tech. Phys. 47,803 (2002). [22] S. Denisov, Y. Zolotaryuk, S. Flach, and O. Yevtushenko, Vortex and Translational Currents due to Broken Time-Space Symmetries, Phys.Rev.Lett.100,224102 (2008); V. Lebedev andF. Renzoni,Two-dimensionalrocking ratchetfor coldatoms, Phys.Rev.A80,023422 (2009). [23] The viscosity and density data used in the present paper were obtained from Wolfram|Alpha. [24] R. Clift, J. R. Grace, and M. E. Weber, Bubbles, drops, and particles (Academic, New York, 1978). [25] C. T. Crowe, T. R. Troutt, and J. N. Chung, Particle interaction with vortices, in Fluid vortices (fluid mechanics and its applications), edited by Sheldon I. Green (Springer-Science+Business MediaB. V.,Dordrecht, 1995),Vol. 30,Chap. XIX,pp. 829–861. [26] Farid F. Abraham, Functional dependence of drag coefficient of a sphere on Reynolds number, Phys. Fluids 13,2194 (1970). [27] J. A. Cuesta, N. R. Quintero, and R. Alvarez-Nodarse, TimeShift Invariance Determines the Functional Shape of the Current in Dissipative Rocking Ratchets, Phys.Rev.X3,041014 (2013). [28] J. Casado-Pascual, J. A. Cuesta, N. R. Quintero, and R. AlvarezNodarse, General approach for dealing with dynamical systems with spatiotemporal periodicities, Phys. Rev. E 91,022905 (2015). [29] D. Cubero, V. Lebedev, and F. Renzoni, Current reversals in a rocking ratchet: Dynamical versus symmetry-breaking mechanisms, Phys. Rev. E 82,041116 (2010); A. Wickenbrock, D. Cubero, N. A. Abdul Wahab, P. Phoonthong, and F. Renzoni, Current reversals in a rocking ratchet: The frequency domain, ibid. 84,021127 (2011). [30] N. A. A. Wahab, Infinite-dimensional dynamical systems in mechanics and physics (Springer-Verlag, New York, 1988). 032219-8