Programmable rheotaxis of magnetic rollers in time-varying fields

Abstract

Magnetic microrobots capable of navigating complex fluid environments typically rely on real-time feedback to adjust external fields for propulsion and guidance. As an alternative, we explore the use of field-programmable rheotaxis, in which time-periodic magnetic fields drive directional migration of ferromagnetic particles in simple shear flows. Using a deterministic model that couples magnetic torques to hydrodynamic interactions near a surface, we show that the frequency, magnitude, and waveform of the applied field can encode diverse rheotactic behaviors—including downstream, upstream, and cross-stream migration relative to the flow. We analyze the mechanisms underlying these responses for canonical fields and use this understanding to design complex waveforms that optimize migration speed and direction. Our results reveal a tradeoff between performance and robustness: high-performance designs enable upstream motion but are sensitive to system parameters, whereas robust designs operate in the linear response regime with more modest performance gains. These findings establish a general strategy for programming flow-guided navigation in magnetic colloids and suggest routes toward self-guided microrobots that respond predictably to fluid environments without external feedback.

---

1 Introduction

Time-varying magnetic fields offer a powerful means of actuating microrobots at the scale of living cells.¹ Such fields can drive a variety of functions, including propulsion via swimming or surface rolling,^2–5 reconfigurable shape changes and cargo transport,^6–8 and mechanical deformation of surrounding materials.^9–11 Because magnetic fields penetrate biological media with minimal attenuation and negligible physiological impact,¹² they provide a biocompatible energy source for medical applications such as targeted delivery, minimally invasive surgery, and local diagnostics.^13–15 However, such fields act uniformly over large spatial domains, thereby limiting their ability to direct the independent behaviors of multiple particles in heterogeneous surroundings.

One strategy to overcome these limitations is to develop microrobots that guide their own motion in response to local environmental cues.^16–18 While feedback control systems can dynamically adjust magnetic fields to direct individual particles along prescribed trajectories,^19,20 these approaches require continuous tracking and are difficult to scale in systems with many particles or limited visual access. Autonomous behaviors inspired by biological taxis—such as artificial chemotaxis,²¹ phototaxis,²² and rheotaxis²³—offer a compelling alternative. In such systems, the self-propelled motion of an active particle is directed by local gradients without the need for global information or external actuation.²⁴ Achieving similar self-guided behavior in magnetically driven particles presents difficulties, as they typically move along field-specified directions and thus lack the autonomy characteristic of self-propelled particles.

Nevertheless, recent work suggests that time-varying fields with appropriate symmetries can encode the self-guided navigation of magnetic particles in response to local gradients.¹ In particular, toggled rotating fields direct the migration of ferromagnetic spheres up topographic gradients on solid substrates, thereby enabling multiple particles to navigate complex landscapes in a common field.²⁵ As explained by predictive models, the time-averaged response derives from interactions between field-driven particle rotation and asymmetries in the particle environment. By tuning the waveform of the driving field, these models predict other responses—including rapid migration up, down, or perpendicular to a local incline.²⁶ These results motivate efforts to extend field programmable migration strategies to other types of environmental gradients, including fluid shear.

Rheotaxis describes the directed migration of micro-organisms^27,28 or active particles^23,29–33 in response to external shear flows. Such flows create a locally anisotropic environment defined by three mutually orthogonal directions: the flow direction, the velocity gradient direction, and the vorticity direction. Upstream migration against the flow direction—termed positive rheotaxis—is often a primary focus; however, other responses are possible, such as cross-stream migration in the vorticity direction.³¹ While some micro-organisms exhibit rheotaxis in bulk flows,³⁴ it is more commonly observed near solid boundaries, where hydrodynamic interactions break symmetry and bias propulsion. For example, spermatozoa and flagellated bacteria migrate upstream due to the “weathervane effect,” in which shear flow reorients the swimmer about a pivot near the wall.^29,35 Analogous behavior is observed in phoretic colloids, where competition between external shear and self-generated interfacial flows leads to a stable upstream orientation near a boundary.^23,29,33

Here, we investigate the possibility of field-programmable rheotaxis in a ferromagnetic sphere positioned above a planar wall in simple shear flow (Fig. 1). We construct a deterministic model of the particle's rotational and translational dynamics by coupling magnetic torques to low-Reynolds-number hydrodynamics in the sphere-wall geometry. By tuning the frequency, magnitude, and waveform of the time-periodic driving field, we identify a range of rheotactic behaviors—including upstream, downstream, and cross-stream migration—relative to the flow direction. For each driving field, we quantify the rheotactic velocity as a function of two key dimensionless parameters: the ratio of the field frequency to the particle's magnetic relaxation rate, and the ratio of shear-induced and magnetic torques. We show that oscillating fields directed normal to the surface enhance downstream migration, while in-plane rotating fields induce cross-stream migration perpendicular to the flow direction. Using evolutionary optimization algorithms, we identify driving protocols that achieve upstream migration and evaluate their sensitivity to system parameters. We discuss the feasibility of realizing these behaviors experimentally, as well as strategies for enhancing their magnitude through anisotropic particle shapes. These results demonstrate a general strategy for designing time-periodic fields that enable self-guided microrobots with programmable responses to asymmetric fluid environments.

Fig. 1 (a) A spherical particle of radius a with magnetic moment m is immersed in a viscous fluid above a plane wall. The particle is subject to both an external shear flow v(x) with shear rate [small gamma, Greek, dot above] and a time-varying magnetic field B(t). (b) Depending on the waveform of the magnetic field, the particle can exhibit (i) enhanced migration in the direction of flow, (ii) lateral migration in the ± vorticity directions, and (iii) upstream migration against the direction of flow. Such field-driven rheotaxis is invariant to changes in the direction of flow.

2 Model

We consider a spherical particle with radius a and magnetic moment m immersed in viscous fluid at a height h above a solid plane at z = 0. The sphere is subject to a time-varying magnetic field B(t) and to an external shear flow with velocity v(x) = [small gamma, Greek, dot above]ze_y where [small gamma, Greek, dot above] is the constant shear rate (Fig. 1). In a uniform magnetic field, the particle experiences a torque, T_m = m × B(t), but no magnetic force, F_m = 0. Other forces and torques—for example, those due to gravity and thermal fluctuations—are assumed negligible by comparison. We further neglect the influence of magnetic and hydrodynamic interactions among neighboring particles, focusing instead on the dynamics of an individual particle within a sufficiently dilute suspension. The validity of these assumptions under realistic experimental conditions is discussed further in the SI.

In the absence of inertial effects (i.e., at low Reynolds number), the linear and angular velocity of the particle depend linearly the external forces and torques as

where a, b, etc. are components of the hydrodynamic mobility matrix.³⁶ Here, F_s and T_s denote the hydrodynamic force and torque on a stationary sphere above a plane wall in a uniform shear flow. In the limit of small surface separations (h → a), these quantities have been computed previously to be F_s = 1.701 × 6πηa²[small gamma, Greek, dot above]e_y and T_s = −0.6293 × 6πηa³[small gamma, Greek, dot above]e_x where η is the fluid viscosity.³⁷ The components of the mobility matrix have the formwhere δ is the identity tensor, and ε is the Levi-Civita tensor. The mobility coefficients Y_a, X_a, etc. depend on the surface separation scaled by the particle radius ξ = (h − a)/a as described by semi-analytical solutions.^38–42 We consider a constant surface separation of ξ = 0.01 for which the relevant mobility coefficients are Y_a = 0.2961, Y_b = 0.03654, X_c = 0.6339, Y_c = 0.3392.²⁵ Together with the kinematics of rigid body motion, eqn (1) fully describes the dynamics of the sphere.

2.1 Non-dimensionlization

To facilitate numerical solution and analytical analysis of the particle dynamics, it is convenient to non-dimensionalize the problem using the following scales:

Here, B₀ is the characteristic magnitude of the applied field, and ω₀ is a magnetic relaxation rate for particle rotation about an axis parallel to the surface. With these units, the shear-induced force and torque are expressed as

where Γ = [small gamma, Greek, dot above]/ω₀ is the ratio of the shear rate and the magnetic relaxation rate. The dimensionless coefficients α and β are approximated as α = 0.577 and β = 0.213 for the chosen surface separation of ξ = 0.01 and depend weakly on this variable when ξ ≪ 1. The dimensionless dynamics of eqn (1) is further parameterized by the mobility ratios μ = Y_a/Y_c, κ = Y_b/Y_c, and λ = X_c/Y_c, which also depend weakly on surface separation with estimates κ = 0.108, μ = 0.873, and λ = 1.87 for ξ = 0.01 (see ref. 25 and 26 for plots of these ratios versus ξ). Below, we use the same notation to refer to dimensionless quantities, which corresponds to setting the characteristic scales to unity in the dynamical eqn (1) (e.g., a → 1).

2.2 Euler angle dynamics

The particle orientation is parameterized by the Euler angles u = [ϕ,θ,ψ]^T (using the most common 313 rotation sequence⁴³); these angles evolve in time as

[small phi (variant), Greek, dot above] = −cot θ(B_x cos ψ + B_y sin ψ) − (β + κα)Γ csc θ sin ψ (7)

[small theta, Greek, dot above] = −cos θ(B_y cos ψ − B_x sin ψ) − B_z sin θ − (β + κα)Γ cos ψ (8)

where B_x(t), B_y(t), and B_z(t) are the time-dependent components of the magnetic field (in the lab frame). The in-plane components of the particle velocity are given by

U_x = κ(B_x cos θ − B_z sin ψ sin θ) (10)

ΔU_y = U_y − U^nm_y = κ(B_y cos θ + B_z cos ψ sin θ) (11)

where U^nm_y = (μα + κβ)Γ is the linear velocity of a non-magnetic, force- and torque-free particle in the shear field. These equations are integrated numerically for a specified field B(t) to determine the particle's transient position x_p(t) and orientation u(t).

2.3 Experimental relevance

For comparison to experiments on magnetic Janus spheres,^44–46 we consider a = 2 μm particles with magnetic moment m = 3 × 10⁻¹⁴ A m² subject to fields of strength B₀ = 5 mT, typical of air-cooled electromagnets. Assuming a surface separation of ξ = 0.01, the magnetic relaxation rate is ω₀ = 340 rad s⁻¹ in water (η = 1 mPa s). More generally, the magnetic moment m of an engineered particle is typically chosen to prohibit irreversible dipole–dipole interactions between contacting particles, which requires that the field-induced torque mB₀ exceed the dipole–dipole energy μ₀m²/16πa³ such that m < 16πa³B₀/μ₀ where μ₀ is the vacuum permeability.¹ Using this upper bound for the magnetic moment, the magnetic relaxation rate becomes independent of particle size: ω₀ < 8Y_cB₀²/3μ₀η ≈ 18 000 rad s⁻¹.

For further context, we consider the flow conditions encountered by engineered particles within the human vasculature. Such particles should be comparable in size to red blood cells (a ≈ 2 μm) to allow access to capillary vessels while minimizing the effects of Brownian motion. As they move through the vasculature, particles encounter shear rates ranging from [small gamma, Greek, dot above] = 20 s⁻¹ in veins to 2000 s⁻¹ in small arteries.⁴⁷ Using the above estimate for the relaxation rate of a magnetic Janus sphere, the dimensionless shear rate varies from Γ = 0.06 to 6 across these conditions. Depending on the targeted shear conditions, the parameter Γ can be further tuned by varying the magnetic moment m through addition of more or less magnetic material during particle fabrication.

3 Results & discussion

To understand how magnetic actuation directs rheotactic motion, we begin by analyzing the response of a ferromagnetic sphere to a set of simple, time-periodic fields. These include static fields, which induce weak upstream migration relative to a non-magnetic particle; oscillating fields normal to the wall, which enhance downstream migration via field-driven rotation; and in-plane rotating fields, which generate cross-stream migration perpendicular to the flow direction (Table 1). These canonical cases reveal how magnetic actuation couples with shear-induced torques to produce time-averaged particle migration. Building on this foundation, we design and interpret more complex, optimized fields that modulate both the speed and direction of rheotactic motion. Using evolutionary algorithms, we identify waveforms that amplify cross-stream transport or drive upstream migration against the imposed flow.

Table 1 Summary of simple driving fields B(t) and the time-averaged rheotactic response 〈ΔU〉 relative to a non-magnetic sphere for flow in the e_y direction

Field	B(t) =	Response
Static	ez	Upstream, −ey
Oscillating	cos ωtez	Downstream, +ey
Rotating	cos ωtex + sin ωtey	Cross-stream, ±ex, ±ey

---

3.1 Static field

We first consider the particle dynamics in a static field of unit magnitude directed normal to the surface, B = e_z. This configuration serves as a baseline to understand how magnetic torques compete with shear-induced rotation to influence particle migration. Under this field, the particle's magnetic moment relaxes to the yz plane, where it's described by the Euler angle θ representing the inclination from the vertical z-direction such that m = −sin θe_y + cos θe_z. The angular dynamics of eqn (8) simplifies to

[small theta, Greek, dot above] = Ω^nm_x − sin θ (12)

where Ω^nm_x = −(β + κα)Γ is the angular velocity of an otherwise identical non-magnetic particle in the same shear flow. The translational dynamics of eqn (11) simplify to

ΔU_y = κ sin θ (13)

with no cross-stream motion in the x-direction perpendicular to the flow (U_x = 0).

Fig. 2 shows the time-averaged velocity difference 〈ΔU_y〉 between a magnetic and non-magnetic particle as a function of shear rate Γ. For Γ < Γ_c = (β + κα)⁻¹, the magnetic torque is sufficient to suppress shear-induced rotation, locking the particle orientation at a fixed angle. In this regime, the magnetic particle translates more slowly than a non-magnetic particle by ΔU_y = κΩ^nm_x = −κ(β + κα)Γ, which increases linearly with shear rate. Above the critical shear rate, shear-induced torques overcome magnetic alignment, leading to time-periodic particle rotation along the surface. This “slipping” behavior leads to a reduction in the time-averaged velocity difference, which follows the relation . At the transition point Γ = Γ_c, the velocity difference reaches a maximum value of −κ.

Fig. 2 Static field. Time-averaged velocity difference 〈ΔU_y〉 between a magnetic and non-magnetic particle as a function of shear rate Γ, for a static magnetic field B = e_z directed normal to the surface. The critical shear rate is Γ_c = (β + κα)⁻¹ at which the velocity difference is −κ.

In dimensional units, this maximum velocity difference is ΔU_y = −κaω₀ = −Y_bmB₀/6πηa², which corresponds to −73 μm s⁻¹ for the magnetic Janus sphere described above. For comparison, a non-magnetic particle under the same flow conditions moves at U_y = 1.3 mm s⁻¹, indicating that torque-driven motions are relatively slow compared to particle convection with the flow. For an engineered particle with moment m ∝ a³B₀, this response scales linearly with particle size a and quadratically with the applied field strength B₀. The particle-wall surface separation ξ influences the strength of rotation-translation coupling as described by the mobility coefficient Y_b. The maximum velocity difference is achieved at a critical shear rate of [small gamma, Greek, dot above]_c = ω₀(β + κα)⁻¹, which corresponds to 1200 s⁻¹ for the magnetic Janus sphere and can be further enhanced by increasing the magnetic moment m. Below, we focus on field-driven particle dynamics at smaller shear rates Γ < Γ_c, where the characteristic magnetic torque is stronger than the shear-induced torque.

3.2 Oscillating field

We next consider particle dynamics under a time-periodic field with constant direction and oscillating magnitude

B(t) = cos ωte_z (14)

As in the static case, the particle's magnetic moment relaxes to and evolves within the yz plane as described by the Euler angle θ. The angular dynamics of eqn (8) simplify to

[small theta, Greek, dot above] = Ω^nm_x − cos ωt sin θ (15)

while the translational dynamics reduce to

ΔU_y = κ cos ωt sin θ (16)

There is no cross-stream motion in the x-direction (U_x = 0). Using the estimated mobility ratios, we integrate these equations numerically to determine the time-averaged velocity difference 〈ΔU_y〉. Further analytical analysis of the particle dynamics in the oscillating field is provided in the SI.

Fig. 3a shows the time-averaged velocity difference versus the driving frequency ω for different values of the shear rate Γ. These calculations identify two critical frequencies ω₁ and ω₂, which mark qualitative changes in the particle's rotational dynamics. Between these frequencies (ω₁ < ω < ω₂), the particle completes one full rotation per cycle of the oscillating field, resulting in synchronized dynamics with 〈[small theta, Greek, dot above]〉 = −ω. Averaging eqn (15) and (16) and substituting this result gives the average translational velocity difference

〈ΔU_y〉 = κ(ω + Ω^nm_x) for ω₁ < ω < ω₂ (17)

shown as dashed lines in Fig. 3a. In this synchronized regime (shaded in Fig. 3b), downstream migration is significantly enhanced, reaching a maximum as the frequency approaches ω₂.

Fig. 3 Oscillating field. (a) Time-averaged velocity difference 〈ΔU_y〉 between a magnetic and non-magnetic particle as a function of field frequency ω for different shear rates Γ. Markers denote results from numerical integration; the dashed line shows the approximation of eqn (17), valid for frequencies ω₁ < ω < ω₂. (b) The frequencies ω₁ and ω₂ increase with shear rate and bound the shaded region, where particle rotation is synced to the oscillating field.

At higher frequencies (ω > ω₂), synchronized rotation is lost when the combination of the magnetic and shear torques become insufficient to rotate the particle at speeds required by the oscillating field. Because the shear flow enhances particle rotation, the step-out frequency ω₂ increases with increasing shear rate Γ (Fig. 3b). Above this critical value, the field moves through multiple oscillation cycles for every full rotation of the particle. In the limit as ω → ∞, the rapidly oscillating field allows for shear-induced particle rotation, and the average particle velocity approaches that of a non-magnetic particle (〈ΔU_y〉 → 0; see SI for details).

At lower frequencies (ω < ω₁), synchronized rotation is lost when shear flow alone rotates the particle through multiple revolutions during that portion of the oscillation cycle when the field strength is low enough to allow shear-induced rotation. In the low shear regime (Γ ≪ 1), the critical frequency ω₁ scales quadratically with the shear rate (Fig. 3b), since shear-induced rotation scales like Γ as does the fraction of the oscillation cycle during which such rotation is allowed (see SI). In this regime, synchronized rotation in the oscillating field can enable large enhancements in the particle's angular velocity relative to that of a non-magnetic particle. For Γ = 10⁻³, the average angular velocity increases 500-fold from 〈[small theta, Greek, dot above]〉 = Ω^nm_x = −2.7 × 10⁻⁴ for a non-magnetic particle to 〈[small theta, Greek, dot above]〉 = −ω for a magnetic particle driven by an oscillating field at the optimal frequency ω = ω₂(Γ) = 0.13.

These results demonstrate that oscillating fields can significantly increase the speed of particle migration along the flow direction. In the absence of magnetic actuation, the particle travels at the baseline velocity U^nm_y = (μα + κβ)Γ. In the absence of shear, the oscillating field alone produces no net migration without inertial effects (neglected here), which can enable symmetry-breaking instabilities.⁴⁸ Only through the combination of oscillatory torque and steady shear does the particle achieve directed migration—enhancing downstream motion by ΔU_y ≈ 0.042 at Γ = 1. The relative magnitude of such enhancements (ΔU_y/U^nm_y) is greatest at low shear rates, where the speed of field-driven particle translation can exceed that of shear-driven motion by more than an order-of-magnitude. Importantly, while the field powers particle motion, the direction of migration is dictated by the flow, a hallmark of rheotaxis.

3.3 Rotating field

Other magnetic field protocols can direct cross-stream rheotaxis, whereby particles migrate perpendicular to the flow direction. The simplest such protocol is a rotating field with unit magnitude and angular velocity ω = ωe_z

B(t) = cos ωte_x + sin ωte_y (18)

In the absence of the external shear (Γ = 0), the particle rotates about the z-axis but does not translate. For rotation frequences below the step-out threshold ω_c = λ, the particle rotates in phase with the applied field, Ω = ωe_z. Above this critical frequency, the magnetic torque cannot overcome viscous resistance, and the average angular velocity decreases as .

In the presence of shear, the flow tilts the axis of rotation, resulting in net particle motion along the vorticity direction (±e_x), depending on the sign of ω. Fig. 4 shows the time-averaged x-velocity as a function of the rotation frequency ω and the shear rate Γ, computed by numerical integration. At low frequencies and shear rates, the oscillation-averaged velocity is well approximated by the following expression derived using perturbation analysis (see SI).²⁶

Fig. 4 Rotating field. Time-averaged particle velocity in the x-direction, perpendicular to flow, as a function of (a) frequency ω and (b) shear rate Γ. Markers show numerical integration results; solid lines show the approximation of eqn (19), valid at low ω and Γ. All quantities are made dimensionless using the characteristic scales of eqn (5).

For positive rotation frequencies and shear rates, the particle translates in the e_x direction perpendicular to the direction of flow. Reversing the rotation direction drives the particle to translate in the opposite −e_x direction.

The rotating field also modulates particle motion in the flow direction (i.e., e_y). At low frequencies and shear rates, the oscillation-averaged velocity is approximated by

The first term (order ω⁰Γ¹) describes the suppression of shear-driven translation due to field-induced resistance to rotation. This slow down is half of that created by a static field, as the particle's shear-induced rotation is prohibited only during a fraction of the rotation period. The second term (order ω²Γ¹) describes a small frequency-dependent enhancement to the downstream velocity due to the rotating field. The perturbative predictions are in quantitative agreement with numerical simulations (Fig. S2).

The rotating field breaks the mirror symmetry of the particle-wall-field system thereby enabling cross-stream migration in the x-direction perpendicular to both the flow direction (y-direction) and the surface normal (z-direction). By contrast, the static and oscillating fields detailed above are symmetric to reflection about the x-direction. The time-averaged particle velocity shares this mirror symmetry thereby prohibiting cross-stream migration for those fields. Both the rotating magnetic field and the linear shear field exhibit a different symmetry that combines 180° rotation about the x-direction followed by reflection about that direction. In the absence of symmetry-breaking instabilities, this rotation-reflection symmetry is also incompatible with cross-stream migration in the x-direction. The presence of the solid wall serves to break the rotation-reflection symmetry of the external fields to allow for the cross-stream migration described by eqn (19). For large surface separations, the mobility coefficient κ approaches zero thereby eliminating the hydrodynamic coupling between magnetic torque and particle translation.

3.4 Designed fields

The simple fields above demonstrate the feasibility of field-driven rheotaxis directed upstream, down-stream, and cross-stream relative to the motion of a non-magnetic sphere. To enhance these rheotactic responses, it is necessary to design the driving field through search and optimization within a space of candidate fields.

Design space. To explore programmable rheotaxis, we define a design space that describes the time-varying magnetic field and its coupling to particle motion under shear. Candidate driving fields are assumed to be time-periodic with frequency ω, satisfying B(ωt + 2π) = B(ωt). Each waveform B(ωt) is constructed to exhibit m-fold rotational symmetry about the surface normal e_z, such that R_z(ν_m)B^±(ωt) = B^±(ωt ∓ ν_m), where the matrix R_z(ν_m) describes a coordinate rotation by an angle ν_m = 2π/m about the z-axis.²⁶ This symmetry prohibits net migration in the absence of flow, as rotation of the field by ν_m is equivalent to a phase shift of ±ν_m. We distinguished between “right-handed” fields B⁺(ωt) and “left-handed” fields B⁻(ωt) based on the sign of this phase shift.

As a representative case, we consider a family of right-handed, time-periodic fields with m = 3 fold rotational symmetry, defined by

where T = ωt is the oscillation phase, and the nine coefficients a₀, a₁, b₁, a₂, b₂, a₃, b₃, a₄, b₄ are tunable parameters. More generally, the design space for the driving field is characterized by the symmetry order m, the oscillation frequency ω, and a corresponding set of Fourier coefficients that define the waveform.

The design space also includes the shear rate Γ, which sets the strength of the shear-induced torque relative to the characteristic magnetic torque. While simulations fix the flow direction along e_y, this direction is assumed to be unknown in practice. Accordingly, we represent the orientation and phase of the driving field relative to the flow using an in-plane rotation angle ν and a phase offset σ, such that B(ωt) = R_z(ν)B_ref(ωt + σ). Finally, we standardize the initial conditions by orienting the particle's magnetic moment along the e_z direction in a static field prior to application of the time-periodic waveform at t = 0. Given these inputs, we integrate the particle dynamics numerically in time to determine the transient orientation u(t) and position x_p(t). The key figure of merit is the oscillation-averaged particle velocity 〈U〉, which characterizes the rheotactic response in the imposed shear flow.

Direct optimization. We begin by designing driving fields that promote upstream rheotaxis—particle migration against the direction of shear flow. To this end, we maximize the oscillation-averaged velocity component −〈U_y〉, which quantifies steady motion in the upstream direction. For each design instance, we fix the oscillation frequency ω, shear rate Γ, and field orientation ν, while varying the waveform coefficients a₀, a₁, b₁, etc. that define the time-periodic field. The average velocity 〈U〉 is computed by numerically integrating the particle dynamics over several cycles. We optimize the waveform using a differential evolution algorithm,⁴⁹ which searches the design space for fields that maximize upstream drift (see SI for details).

Fig. 5 shows the optimized driving field and the resulting particle trajectory for a representative parameters set: m = 3, ω = 0.2, Γ = 0.02, and ν = 0. In the absence of shear, the field drives periodic particle rotation without net translation. The applied shear field breaks the symmetry of the particle environment, resulting in net upstream migration—against the flow—at an average speed of −0.0089 (equivalent to 0.28 particle radii per cycle). In addition to upstream drift, the particle exhibits cross-stream migration along the x-direction. This lateral motion can be eliminated if desired by alternating between the optimized right-handed field and its mirror image reflected about the e_x direction.²⁵

Fig. 5 Optimized upstream rheotaxis. (a) Designed right-handed field with m = 3-fold rotational symmetry that maximizes −〈U_y〉 for ω = 0.2 and Γ = 0.02. (b) Simulated particle trajectory in the xy plane showing three oscillation cycles. The shear flow is directed in the positive y-direction; the inset shows the driving field from above. (c) Oscillation-averaged particle velocity 〈U_y〉 as a function of orientation angle ν between the driving field and flow direction (green markers). Shaded regions indicate bistability: particle dynamics approach one of two stable attractors (black markers), depending on initial orientation or phase.

While effective in achieving the desired rheotactic response, the designed field is sensitive to its orientation ν relative to the flow. As shown in Fig. 5c, variations in ν can reverse the direction of migration, with 〈U_y〉 switching from negative (upstream) to positive (downstream). For certain orientations, the dynamics become bistable: the migration direction depends not only on ν but also on the particle's initial orientation or the phase of the applied waveform. While the orientation-averaged drift remains biased in the upstream direction, this behavior does not constitute true “self-guided” rheotaxis, as it depends on the relative orientation of the field and the flow.

Furthermore, the migration velocity is sensitive to both the oscillation frequency ω and the shear rate Γ, highlighting the complexity of the design landscape. Deviations from the parameter values used during optimization can substantially reduce or eliminate the desired upstream migration. These sensitivities arise from the nonlinear coupling of magnetic actuation to shear-induced torques, which gives rise to nonlinear rheotactic responses. As discussed below, more robust “self-guided” behaviors can be achieved by designing fields that operate in the linear response regime, where the migration velocity scales linearly with the applied shear.

Robust optimization for small ω & Γ. The sensitivities described above can be mitigated by restricting the optimization to the linear response regime—specifically, to small oscillation frequencies ω ≪ ω* and shear rates Γ ≪ Γ*, where ω* and Γ* are field-dependent thresholds to be determined. In this regime, the oscillation-averaged particle velocity 〈U〉 can be expanded as a Taylor series in both parameters

〈U〉 = (〈U₀₁〉 + 〈U₁₁〉ω + 〈U₂₁〉ω²)Γ + … (22)

Here, the vector coefficients 〈U_ij〉 quantify contributions of order ωⁱΓ^j. Notably, the expansion excludes field-driven motion in the absence of shear (Γ = 0), since such contributions vanish for fields with m-fold rotational symmetry. Higher-order terms in ω and Γ are neglected under the assumption that both parameters remain sufficiently small. Importantly, the leading-order coefficients 〈U_ij〉 can be computed efficiently for any candidate waveform B(ωt) using a multiple scale perturbation scheme detailed in the SI.^26,50 By operating within this linear response regime—where 〈U〉 ∝ Γ—the particle's migration velocity becomes effectively independent of the field orientation angle ν, enabling robust and orientation-agnostic rheotactic behavior.

Using the perturbative expansion in eqn (22), we optimize specific components of the migration velocity—such as cross-stream or upstream motion—by varying the oscillation frequency ω, shear rate Γ, and field waveform B(ωt), parameterized by the Fourier coefficients a₀, a₁, b₁, etc. To ensure the validity of the truncated expansion, both ω and Γ are constrained to lie within bounds ω < ω* and Γ < Γ*. As detailed in the SI, these limits are chosen such that higher order terms are successively smaller than the leading order contributions by a specified factor 0 < f < 1. For instance, the time-averaged coefficients in the x direction satisfy the condition, |〈U^x_ij〉| < f^i+j−2|〈U^x₁₁〉|, for positive integers i and j with 0 ≤ i + j ≤ 4. This heuristic ensures that the particle's magnetic moment can continuously follow the trajectory of the driving field, thereby preserving the assumptions underpinning the linear response approximation. Within these bounds, the frequency ω and shear rate Γ are optimized via linear programming, while the waveform coefficients are optimized using differential evolution.

Fig. 6 shows the result of this procedure for maximizing the rate of cross-stream migration, 〈U^x₁₁〉ωΓ. For the optimized field (Fig. 6a), the linear rheotactic response 〈U^x₁₁〉ωΓ, with 〈U^x₁₁〉 = 0.0318, remains valid up to frequencies ω* = 0.122 and shear rates Γ* = 2.29. Under these conditions, the orientation of the driving field relative to the flow direction does has no effect on the oscillation-averaged particle migration (Fig. 6b). Fig. 6c confirms the validity of the perturbative approximation: the numerically computed velocity 〈U_x〉 closely matches the predicted linear trend for small ω and Γ. The optimized cross-stream velocity is ∼220% larger than that of the simpler rotating field analyzed above (Fig. 6c asterisk * vs. circle ○). This modest improvement illustrates a key trade-off: while restricting the dynamics to the linear response regime yields robust and orientation-independent behavior, it limits the performance gains that could otherwise be achieved using more complex time-varying fields operating in the nonlinear regime.

Fig. 6 Optimized cross-stream rheotaxis. (a) Designed right-handed field with m = 3-fold rotational symmetry that maximizes 〈U^x₁₁〉ωΓ. (b) Simulated particle position x_p(t) in the cross-stream direction for the driving field in (a) rotated by different angles ν about the z-axis. The time-averaged migration (black dashed line) is independent of the field's orientation relative to the flow. (c) Time-averaged particle velocity 〈U_x〉 as a function of frequency ω for different shear rates Γ. Markers denote denote the result of numerical integration; solid lines show the perturbation approximation. The asterisk (*) denotes optimum U^x₁₁ω*Γ*; the open circle (○) denotes the analogous quantity for the (sub-optimum) rotating field.

Attempts to enhance upstream rheotaxis in the linear response regime were less successful. At zeroth order in frequency, fields designed to maximize −〈U^y₀₁〉Γ were inferior to the simple static field, which enables upstream migration with a maximum velocity, ΔU_y = −κ, relative to that of a non-magnetic particle (Fig. 2). At second order, fields designed to maximize −〈U^y₂₁〉ω²Γ produce only negligible upstream contributions. These results suggest that other strategies are needed—for example, non-spherical particles—to enable field-driven upstream migration that is invariant to flow orientation.

3.5 Experimental relevance

For context, we compare our results for cross-stream rheotaxis to that of an active particle reported by Upsal, Sánchez, and co-workers.³¹ In their system, the authors consider the migration of platinum-coated Janus spheres (a = 1 μm) across a solid wall under the combined influence of hydrogen peroxide fuel (propulsion speed U_p = 6 μm s⁻¹) and external shear flow (translation speed U* = 14 μm s⁻¹). These fields act to orient the particle's Janus director perpendicular to the flow direction, resulting in particle migration at an angle of approximately arctan (U*/U_p) ≈ 20° from the downstream direction. By contrast, a time-periodic magnetic field optimized for cross-stream rheotaxis (Fig. 6a) drives particle migration at much smaller angles of approximately 〈U_x〉/〈U_y〉 ≈ 0.4° relative to the flow.

The slow rates of cross-stream rheotaxis for the magnetically driven sphere compared to that of an active particle highlight fundamental differences in their respective mechanisms. For the field-driven particle, there is no propulsion velocity in the absence of shear, and the velocity components—both perpendicular and parallel to the flow—scale linearly with shear rate. In this linear response regime, the migration angle for magnetic rheotaxis is independent of the shear rate. By contrast, for the active particle, the shear field orients the propulsion velocity U_p perpendicular to the flow direction; however, its magnitude is largely shear-independent—determined instead by chemically fueled phoretic propulsion. As a result, the migration angle can be increased by reducing the magnitude of the share rate (assuming other contributions due to Brownian motion and gravity remain negligible).

For the prototypical magnetic Janus sphere^44–46 detailed above, cross-stream rheotaxis under optimal conditions drives particle migration at speeds 〈U_x〉 = 6 μm s⁻¹ for shear rates of [small gamma, Greek, dot above] = 800 s⁻¹, typical of the human vasculature.⁴⁷ The migration rate is similar to the propulsion speed of an active particle but small compared to the convective velocity U_y = 800 μm s⁻¹ in the applied shear flow—consistent with the small migration angles noted above.

In principle, such rheotactic responses could be applied to microrobot navigation or particle separations; however, further enhancements are needed to motivate experimental demonstrations. For instance, consider a dilute mixture of magnetic and non-magnetic spheres traveling along the floor of a spiral-shaped microfluidic channel. Application of a time-periodic field optimized for cross-stream rheotaxis (Fig. 6a) would drive the selective migration of the magnetic particles across the full channel width W after traveling a distance of ∼140W along the flow direction. Future enhancements in field-driven rheotaxis could avoid the need for such long channels and enable downstream, upstream, and cross-stream responses at speeds exceeding the flow velocity.

The possibility of such enhancements is motivated by our present results for the oscillating field (Fig. 3) for downstream migration and the designed field for upstream migration (Fig. 5). These fields operate outside of the linear response regime to enable field-driven particle migration at speeds exceeding the flow velocity. Future work should explore the design of such fields in a manner that preserves the rotational-invariance of the rheotactic response as demonstrated by the oscillating field but not by the designed field.

Further enhancements in field-programmable rheotaxis could be achieved by expanding the design space to include anisotropic particle shapes. In particular, ellipsoidal particles introduce new couplings in the hydrodynamic mobility matrix that link torque-driven rotation to particle translation above a solid wall. Unlike spheres, prolate ellipsoids experience orientation-dependent torques in shear flows that tend to align their long axis with the flow direction, thereby enabling new modes of field-responsive behavior.³⁶ In practice, magnetic colloids have been prepared by two common strategies: (1) deposit magnetic materials (e.g., cobalt, iron) onto particle monolayers by vapor deposition^51,52 and (2) incorporate magnetic nanoparticles within polymeric colloids.⁵³ Anisotropic particles have been demonstrated using both methods by stretching polymeric colloids within a sacrificial matrix;^54,55 however, it remains a challenge to quantify and control the direction of particle magnetization relative to particle shape.

4 Conclusions

We have developed a dynamical model to demonstrate how time-periodic magnetic fields can encode directed migration of ferromagnetic spheres near a planar wall under shear flow. By tuning the frequency, magnitude, and waveform of the applied field, we achieve programmable rheotactic responses—including enhanced downstream migration, upstream motion against the flow, and cross-stream drift perpendicular to it. In addition to analyzing simple fields (e.g., static, oscillating, rotating), we design more complex, time-periodic fields to optimize specific rheotactic behaviors. This design process reveals a fundamental trade-off: high-performance responses such as rapid upstream rheotaxis tend to be sensitive to system parameters, whereas more robust, orientation-agnostic designs yield modest performance gains.

Realizing programmable rheotaxis in experimental systems will require strategies to (i) enhance performance within the linear response regime, and/or (ii) improve the robustness of optimized designs operating in the nonlinear regime. Toward this goal, we are exploring the use of non-spherical particles, such as ellipsoids, to strengthen the coupling between magnetic and shear-induced torques. If successful, these advances could lay the groundwork for self-guided microrobots capable of autonomous navigation through complex microstructured environments, such as the human vasculature, where local flows and boundaries naturally guide particle motion for biomedical applications. Such real world environments introduce new challenges to investigate and overcome such as (i) the effects of time-periodic shear rates, (ii) the impact of non-Newtonian rheology, and (iii) the need for rotational invariance in three dimensions rather than two demonstrated herein.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information: Derivation of the dynamical model; assessment of the model assumptions; multiple-scale perturbation analysis and its numerical solution; analysis of particle dynamics in an oscillating field; downstream particle velocity in a rotating field; perturbation results for migration in a precessing field; and implementation details of both the direct optimization approach and the robust optimization based on perturbation theory. See DOI: https://doi.org/10.1039/d5sm00682a

Data for this article, including Mathematica and Python notebooks for computing field-driven particle dynamics, are available on github at https://doi.org/10.5281/zenodo.15792129.

Acknowledgements

This material is based upon work supported by the National Science Foundation under Grant No. CBET-2153202. A. R. received additional support from the Barnard Summer Research Institute (SRI).

Notes and references

1. K. Dhatt-Gauthier, D. Livitz, Y. Wu and K. J. M. Bishop, JACS Au, 2023, 3, 611–627.
2. R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone and J. Bibette, Nature, 2005, 437, 862–865.
3. A. Ghosh and P. Fischer, Nano Lett., 2009, 9, 2243–2245.
4. P. Tierno, R. Golestanian, I. Pagonabarraga and F. Sagués, Phys. Rev. Lett., 2008, 101, 218304.
5. T. O. Tasci, P. S. Herson, K. B. Neeves and D. W. M. Marr, Nat. Commun., 2016, 7, 10225.
6. F. Martinez-Pedrero and P. Tierno, Phys. Rev. Appl., 2015, 3, 051003.
7. T. Yang, T. O. Tasci, K. B. Neeves, N. Wu and D. W. M. Marr, Langmuir, 2017, 33, 5932–5937.
8. K. Han, C. W. Shields IV, N. M. Diwakar, B. Bharti, G. P. López and O. D. Velev, Sci. Adv., 2017, 3, e1701108.
9. T. O. Tasci, D. Disharoon, R. M. Schoeman, K. Rana, P. S. Herson, D. W. M. Marr and K. B. Neeves, Small, 2017, 13, 1700954.
10. J. Leclerc, H. Zhao, D. Bao and A. T. Becker, IEEE Trans. Rob., 2020, 36, 975–982.
11. Q. Wang, X. Du, D. Jin and L. Zhang, ACS Nano, 2022, 16, 604–616.
12. B. J. Nelson, I. K. Kaliakatsos and J. J. Abbott, Annu. Rev. Biomed. Eng., 2010, 12, 55–85.
13. M. Sun, K. F. Chan, Z. Zhang, L. Wang, Q. Wang, S. Yang, S. M. Chan, P. W. Y. Chiu, J. J. Y. Sung and L. Zhang, Adv. Mater., 2022, 34, 2201888.
14. J. Law, X. Wang, M. Luo, L. Xin, X. Du, W. Dou, T. Wang, G. Shan, Y. Wang and P. Song, et al., Sci. Adv., 2022, 8, eabm5752.
15. J. G. Lee, R. R. Raj, N. B. Day and C. W. Shields IV, ACS Nano, 2023, 17, 14196–14204.
16. Y. Dou and K. J. M. Bishop, Phys. Rev. Res., 2019, 1, 032030.
17. A. T. Liu, M. Hempel, J. F. Yang, A. M. Brooks, A. Pervan, V. B. Koman, G. Zhang, D. Kozawa, S. Yang, D. I. Goldman, M. Z. Miskin, A. W. Richa, D. Randall, T. D. Murphey, T. Palacios and M. S. Strano, Nat. Mater., 2023, 22, 1453–1462.
18. A. M. Brooks, S. Yang, B. H. Kang and M. S. Strano, Adv. Intell. Syst., 2024, 6, 2400002.
19. M. P. Kummer, J. J. Abbott, B. E. Kratochvil, R. Borer, A. Sengul and B. J. Nelson, IEEE Trans. Rob., 2010, 26, 1006–1017.
20. J. Jiang, Z. Yang, A. Ferreira and L. Zhang, Adv. Intell. Syst., 2022, 4, 2100279.
21. Y. Hong, N. M. Blackman, N. D. Kopp, A. Sen and D. Velegol, Phys. Rev. Lett., 2007, 99, 178103.
22. C. Lozano, B. Ten Hagen, H. Löwen and C. Bechinger, Nat. Commun., 2016, 7, 12828.
23. J. Palacci, S. Sacanna, A. Abramian, J. Barral, K. Hanson, A. Y. Grosberg, D. J. Pine and P. M. Chaikin, Sci. Adv., 2015, 1, e1400214.
24. K. J. M. Bishop, S. L. Biswal and B. Bharti, Annu. Rev. Chem. Biomol. Eng., 2023, 14, 1–30.
25. Y. Wu, A. Ramos and K. J. M. Bishop, ACS Appl. Eng. Mater., 2024, 2, 2479–2487.
26. Y. Dou, P. M. Tzelios, D. Livitz and K. J. M. Bishop, Soft Matter, 2021, 17, 1538–1547.
27. T. Pedley and J. O. Kessler, Annu. Rev. Fluid Mech., 1992, 24, 313–358.
28. J. S. Guasto, R. Rusconi and R. Stocker, Annu. Rev. Fluid Mech., 2012, 44, 373–400.
29. W. E. Uspal, M. N. Popescu, S. Dietrich and M. Tasinkevych, Soft Matter, 2015, 11, 6613–6632.
30. L. Ren, D. Zhou, Z. Mao, P. Xu, T. J. Huang and T. E. Mallouk, ACS Nano, 2017, 11, 10591–10598.
31. J. Katuri, W. E. Uspal, J. Simmchen, A. Miguel-López and S. Sánchez, Sci. Adv., 2018, 4, eaao1755.
32. R. Baker, J. E. Kauffman, A. Laskar, O. E. Shklyaev, M. Potomkin, L. Dominguez-Rubio, H. Shum, Y. Cruz-Rivera, I. S. Aranson, A. C. Balazs and A. Sen, Nanoscale, 2019, 11, 10944–10951.
33. P. Sharan, Z. Xiao, V. Mancuso, W. E. Uspal and J. Simmchen, ACS Nano, 2022, 16, 4599–4608.
34. G. Jing, A. Zöttl, É. Clément and A. Lindner, Sci. Adv., 2020, 6, eabb2012.
35. F. P. Bretherton and N. M. V. Rothschild, Proc. R. Soc. London, Ser. B, 1961, 153, 490–502.
36. S. Kim and S. J. Karrila, Microhydrodynamics: Principles and selected applications, Dover, 2005.
37. A. J. Goldman, R. G. Cox and H. Brenner, Chem. Eng. Sci., 1967, 22, 653–660.
38. G. B. Jeffery, Proc. London Math. Soc., 1915, 2, 327–338.
39. W. Dean and M. O'Neill, Mathematika, 1963, 10, 13–24.
40. M. E. O'Neill, Mathematika, 1964, 11, 67–74.
41. H. Brenner, Chem. Eng. Sci., 1961, 16, 242–251.
42. A. J. Goldman, R. G. Cox and H. Brenner, Chem. Eng. Sci., 1967, 22, 637–651.
43. J. Diebel, Matrix, 2006, 58, 1–35.
44. W. Fei, M. M. Driscoll, P. M. Chaikin and K. J. M. Bishop, Soft Matter, 2018, 14, 4661–4665.
45. W. Fei, P. M. Tzelios and K. J. M. Bishop, Langmuir, 2019, 36, 6977–6983.
46. D. Livitz, K. Dhatt-Gauthier and K. J. M. Bishop, Soft Matter, 2023, 19, 9017–9026.
47. M. H. Kroll, J. D. Heliums, L. V. Mclntire, A. I. Schafer and J. L. Moake, Blood, 1996, 88, 1525–1541.
48. A. Kaiser, A. Snezhko and I. S. Aranson, Sci. Adv., 2017, 3, e1601469.
49. R. Storn and K. Price, J. Global Optim., 1997, 11, 341–359.
50. S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering, CRC Press, 2018.
51. S. K. Smoukov, S. Gangwal, M. Marquez and O. D. Velev, Soft Matter, 2009, 5, 1285–1292.
52. B. Ren, A. Ruditskiy, J. H. Song and I. Kretzschmar, Langmuir, 2012, 28, 1149–1156.
53. A. Spatafora-Salazar, D. M. Lobmeyer, L. H. Cunha, K. Joshi and S. L. Biswal, Soft Matter, 2021, 17, 1120–1155.
54. J. G. Lee, A. Al Harraq, K. J. M. Bishop and B. Bharti, J. Phys. Chem. B, 2021, 125, 4232–4240.
55. H. M. Gauri, R. Patel, N. S. Lombardo, M. A. Bevan and B. Bharti, Small, 2024, 20, 2403007.

---