Reconfigurable Vortex-like Paramagnetic Nanoparticle Swarm with Upstream Motility and High Body-length Ratio Velocity

Drug delivery systems with high-targeted doses can minimize excipients, reduce side effects, and improve efficacy. Human blood circulation is a complex circulatory system, and the motion control of microrobots in the static flow field in vitro is completely different from in vivo. How to achieve precise counterflow motion for targeted drug delivery without vascular blockage and immune rejection is the biggest challenge for micro-nano robots. Here, we propose a control method that enables vortex-like paramagnetic nanoparticle swarm (VPNS) to move upstream against the flow. By mimicking the clustering motion of wild herring schools and the rolling of leukocytes, VPNS are incredibly stable even when subjected to high-intensity jet impacts in the blood environment, can travel upstream, anchor at the target location, and dissipate when the magnetic field is withdrawn, which greatly reduces the risk of thrombosis. VPNS can also upstream along the vessel wall without an additional energy source and has a marked targeted therapeutic effect on subcutaneous tumors.

Where ρ is the fluid density, μ is hydrodynamic viscosity, υ is kinematic viscosity. For example, when a particle with a radius of 40μm is rotating at an angular velocity of 100 rad/s in water with a kinematic viscosity of 1.007cSt (20℃), the Reynolds number is 2.56×10 -10 , which conforms to ≪ 1, which means the influence of inertia on particle motion can be ignored. Thus, equation (1) can be simplified as: By substituting variables, equation (3) can be reduced to: 4 Where -, , / , ⁄ . Equation (4) formally conforms to nonuniform oscillators. Fig. S3 shows the vibration state of the nonuniform oscillators system as the ratio / changes.
Assuming that transient decay to zero as well as the system reaches steady-state. When / 1, equation (4) approaches to the fixed point / , therefore in a stable state with =0. When / 1, average velocity within a period is calculated: The vibration period is as follows: T 1 2 6 Thus, the average angular velocity of the particle is: Where is the critical angular velocity of a particle, corresponding to the critical frequency is (Fig. S4).
The above theory can be summarized as follows: regardless of inertia and Brownian forces, the rotation of a spherical magnetic particle is controlled by the interaction of magnetic force moment and viscous moment and exhibits two different modes of motion: synchronous and asynchronous rotation.

When
, the phase lag between the particle and the magnetic field is constant because the particle rotates synchronously with the magnetic field because of the balance between magnetic force moment and vicious moment.

When
, the binding of the particle by the viscous force exerted by the surrounding liquid environment exceeds magnetic driving force, and the phase lag between particle and magnetic field increases with the increase of frequency, indicating that particle no longer rotates synchronously with the magnetic field and its angular velocity decreases with the increase of magnetic field frequency.
The Mason number is defined as the dimensionless ratio of the magnetization to the viscous force of the particles, that is, . . When there are only two particles, the time t is defined as the shear time required for one particle to migrate under the influence of the other, and denotes the shear time of a particle moving at an angular frequency . When ( 1), the magnetic force acting on the particles exceeds the viscous force; hence, the particles aggregate into short chains. As Mn increases, the magnetic force operating on the particles tends to decrease, and chains of a shorter length are formed, especially when ( 1). The shearing time of the particles is less than the time required for the particles to move under the magnetic field, which results in no stable chains being established. Neglecting the magnetization relaxation effect of MNPs, the particles' angular frequency of the particles is approximately equal to the magnetic field rotation frequency, and Mn is described as the most relevant function of and n (the number of magnetic fields particles). The specific number of particles in the short-chain was measured from 20 to 40, and the median n = 30 was selected for the analytical treatment. The Mason number and magnetic field frequency curves are shown in Fig. S8A. Theoretically, if the magnetic field frequency exceeds approximately 20 Hz, the particle chain is destabilized, consistent with the experimental results (Fig. S8B).
Wall effect When a spherical particle on the bottom surface rolls, it moves horizontally because of the role of horizontal flow in a particle's physical imbalance; the force analysis is shown in Fig. SA. Furthermore, velocity distribution function P(V) as shown in Fig. S5B around 4 mm/s at peak. The size of the peak depends on the driving magnetic field frequency. Moves on the nonslip boundary, particle's speed limit is . Fig. S5C shows that the rolling velocity distribution of particles varies with frequency, and the velocity distributes around 0.2V0 ~ 0.3 V0.
Two main factors can cause velocity to decline. First, fluid can act as a lubricant, causing partial slip between particle and boundary. Second, a rotating particle moving horizontally near the solid surface experience a lift force , which is opposite to the direction of gravity. According to Stokes' law, it can be deduced that the lift force of rotating particle is:

8
Where is particle translation velocity. When is equal to gravityG 4/3π , the particle has no contact with the bottom surface.

Using
, the critical condition ( / 3 /4 1) for the particle to escape from the bottom surface can be obtained. Aranson et al.[53] observed that for a nickel particle with a radius of 60~100 μm ( / 9), the lift force will exceed gravity while 60~100 Hz, which means that the particle cannot contact the bottom for a long time above this frequency. Accordingly, it should move parallel to the surface at speed less than .
Fang et al. [54] research shows that when a particle rotates near the boundary, the horizontal viscous force of the particle is:

9
Where: 3 4 ℎ 1 3 8ℎ 10 When a particle moves horizontally, the drag force in the horizontal is:

11
Where: In addition, the bottom surface of a particle is subjected to the friction of the fluid.

13
Where is the fluid friction coefficient.
The curves of the magnetic particle viscous force and fluid drag force vary with distance from the bottom surface are shown in Fig. S6.

Note S2. Interpretation of the transformation from MNP suspension to vortex
Under the action of the external magnetic field, the spinning particle generates a small vortex field with itself as the vortex core. Spin-induced hydrodynamic interaction and magnetic dipole force attract nearby particles to get close to each other, while the short-range repulsion prevents the particles from getting too close together. As more particles continue to gather, a vortex is formed.

Emergence of vortex
The Formation Process Vortex exerts long-term attraction on adjacent systems as a dynamic gathering system. When two vortices are far apart, that is, the center distance d of two vortices is greater than radius of a vortex, vortices rotate slowly. When / exceeds a particular critical value, the attraction between vortices leads to vortex fusion. The phenomenon of vortex merging largely depends on the initial value of ( / ). When the ratio is lower than 0.3, the primary vortex forms a large stable vortex by swallowing small vortices one by one along a spiral trajectory.
This process can be summarized into three stages exchange, fusion, and stabilization, as shown in Fig. S9.
First, once two vortices are close enough, the vortices move towards each other driven by vortex advection (long-range attraction interaction), and the two closest parts of vortices begin to exchange members.
Then, when two vortices contact each other, all members enter the fusion zone, and an irregularly shaped vortex forms quickly. Meanwhile, the vorticity field merges into one through the diffusion process. Next, the vortex merges into an ellipse and spits out two strong vortex filaments.
Eventually, a stable circular symmetrical vortex is formed, fusing all members of these two vortices, while filaments gradually curl up and dissipate around the vortex core.
Surface fraction Since density-dependent diffusion in non-equilibrium systems leads to phase separation and a large number of fluctuations, and the surface fraction has a significant effect on the dynamic assembly of vortices. Research shows that in a given period, the surface fraction is lower levels (φ 2.08%) observed only a few small vortexes (less than 5 particles), and big vortex present at higher levels (φ 11.99%) [25].
Furthermore, the critical value R/d can estimate the minimum surface fraction required for vortex formation[27].

Parameter of collective dynamics
Diameter Magnetic force between two magnetic particles can be expressed as Where vector and are the magnetic moment of particle 1 and particle 2, vector is the center distance between particles 1 and 2, μ0 is vacuum permeability [55].
Model-1 When the particle cluster is regarded as a rigid disk (Fig. S10), the particle at edge conforms to the law of circular motion, and the component of magnetic force along the central line ( ) provides the centripetal force.

15
After a certain simplification of Equation (15) ( ), the cluster radius varied with its rotational speed can be obtained.

16
Where is mass of particle 2, and , represent radial component along the central line, is the sum of the vector product related to the magnetic moment.
If magnetic moments of particle 1 and particle 2 are assumed to be parallel, which means the direction of magnetic force is along the centerline, then the equation (17) can be simplified as 2 , and its radial component is itself. Total number of particles in the cluster is / . Thus, equation (16) can be further simplified to: 18 Equation (18) shows that the cluster radius is related to the magnetic moment m2， density , particle radius of a single particle and magnetic field angular velocity Ω. Comparing theoretical result and experimental data in Fig. 2E, the overall trend is followed with frequency increases while the cluster size decreases correspondingly, which confirm that the theoretical model is valid. However, theoretical and experimental values have a significant deviation in lowfrequency.
Model-2 When cluster rotates vertically, considering the influence of gravity (a component of gravity along the central line ), equation (15) can be modified as:

19
After simplification, we can get the modified cluster size. Where cos is a component of gravity along the direction of the centerline, and θ is the angle between gravity and the direction of the centerline. Therefore, is a function with . Takes the average value of during 0~2 , and a specific cluster size can be obtained. Moreover, comparing these two models shown in Fig. S11, model-2 shows that the change rate of cluster size at low frequency is larger than model-1 after considering gravity; in other words, gravity has a more significant impact on cluster size than the magnetic force at low frequency.
When the cluster moves horizontally near the wall, fluid resistance at the edge inhibits the cluster's size. However, since fluid resistance is related to the overall translational velocity of the cluster, the model involved is too complex, so the influence of fluid resistance on the size of the cluster is not considered in this work.

Note S3. Flow field distribution model for microfluidic flow channels with a square crosssection
Considering the low-Re fluid environment in the experimental conditions, we combine the Navier-Stokes equations with the continuity equation for viscous incompressible peristaltic flow as follows: where ⃗ is the body force per unit mass. We assume that a constant pressure gradient exists within the fluid, that is:

∆ . 23
The hydraulic conductivity[30] is: where (ratio of area to the perimeter) represents the morphological factor, and denotes the ratio of the side lengths for a square section , 1.
We define ∆ . Corresponding to the volumetric flow input used in this study: . 25 Thus, deriving the fluid velocity : The graph of , was plotted using MATLAB (Mathworks), shown in Fig. S21.
Note S4. Dynamics model of simplified VPNS upstream rolling in the Poiseuille flow Rigid disc model The VPNS converts rotational motion into translational motion near the wall, owing to the slip-free boundary of the hydrodynamics. This means that the fluid may slide on the particle surface but not on the plain wall. This hinders the rolling of the VPNS near a surface with mismatched fluid pressure at its top and bottom, thereby inducing different velocities at each side. The snapshot from the experimental side view demonstrates that there is a highly narrow lubrication distance h between the VPNS and the wall when it is rotating upstream, which originates from the pressure difference generated by the high-speed rotation of the VPNS in the viscous fluid that makes the integral being held up by the lubricating film. Thus, it is suitable to apply the lubrication theory to analyze the forces on the VPNS. We simplify the VPNS as a microdisk of radius R and thickness t. Because motion in a super-viscous environment is overdamped, we assume a total force balance on the rotating VPNS in the Stokes state; the propulsive and drag forces in the translational direction are balanced. Moreover, the gravitational, lift, and buoyancy forces perpendicular to the propulsive direction are also balanced. Specifically, the propulsive force generated by wet friction can be expressed as:

, 29
where , ⁄ is the boundary friction coefficient, which is inversely proportional to the radius of the rotating disc, and , is a function of the material and surface of the target object. where ⁄ is the shear Reynolds number that characterizes the fluid inertia, k is the fluid shear rate, and is the rotational angular velocity. Then, the dynamic equilibrium equation is established according to the force analysis: , 35 or: In the vertical direction, clusters experience Saffman lift force caused by cluster rotation, support force provided by liquid, buoyancy force comes from liquid, Magnus lift force arises from shear flow and cluster gravity. Therefore, we deduce the upstream velocity of the VPNS : , 84.63 1⁄ 2 2 ℎ ⁄ 2 t . 37 The curve of the upstream migration velocity is plotted by taking the diameter of VPNS as the independent variable, as shown in Fig. S25. The gradient of the translational velocity decreases and tends towards a constant value as the swarm diameter increases. This may be owing to the faster fluid velocity at the edge of larger diameter swarms, which hinders the increase in migration velocity.
The gap between cluster and boundary As we assume above, the cluster is regarded as a homogeneous magnetic disk, whose thickness is , the radius is , and total magnetic moment is . According to the distribution of particles in the cluster, the magnetic moment of each layer is / , is the magnetic moment of each particle, and the number of cluster layers is /2 , so the total magnetic moment is Where , represents fluid velocity at the upper and lower surfaces of the cluster, respectively. For example, the velocity distribution of a one-dimensional circular tube is 1 2 40 Where r is the distance between the cluster center and the center axis of the pipe, by solving equations (36) and (39), we can obtain the higher-order equation about h. Mathematical software can only obtain an approximate numerical solution since there is no analytic solution for a higherorder equation.
For example, 370 μm/s , 300 μm, and frequency 20~200 Hz are set for onedimensional circular pipe flow. Moreover, the height between the lower surface of the cluster and the wall is shown in Fig. S12B. With the increase of frequency, the height increases gradually from 0.1352 μm to 0.2946 μm, and the corresponding cluster velocity increases from -12.65 μm/s to 119.03 μm/s. The inset shows the cluster velocity when different height values are set. When the height increases, the same diameter's velocity also increases, indicating that height cannot be ignored while calculating the cluster velocity.

Note S5. Manipulation of clusters
According to equation (16), the cluster velocity can be derived. As we know from the above, since ℎ is a numerical solution, so is . The analytical expression of cluster velocity is listed here to show cluster velocity factors more intuitively.
We assume that the vessel center velocity is 2 m/s and the vessel diameter is 2 cm, so its velocity distribution is shown in Fig. S13A. As the frequency increases, the cluster diameter decreases rapidly to a flat level, and the cluster velocity also increases rapidly to a flat level. When the frequency reaches 2500 Hz, the cluster velocity is still negative, indicating that the cluster cannot move against the flow and can only flow with the liquid at a slow speed (Fig.  S13B).
Considering the critical frequency / 2448 Hz of a single particle, the critical frequency of cluster synchronous movement with the magnetic field is smaller than that of a single particle. The premise of equation (19) is that the cluster rotation is synchronized with the magnetic field; in other words, the magnetic field frequency is equivalent to the cluster frequency. When the increasing magnetic field frequency exceeds the critical cluster frequency, this equation is no longer accurate, which means the result obtained by increasing the magnetic field frequency to improve the cluster velocity is invalid.   S2. Force on a single magnetic particle in a magnetic field. Phase lag (φ) between magnetic moment (m) for a spheroidal particle and magnetic field (B), is the angular velocity of the magnetic field, θ is the rotational angle of the particle. 1 the system stops completely, and a semi-stable fixed point / appears at π/2. (Ⅲ) / 1 the system is unstable and showing periodic changes where fastest at -π/2 and lowest at π/2.

Fig. S4. Particle rotation frequency ( ) and magnetic field frequency ( ). When
, . is the critical frequency of a particle. When , decays rapidly.  Fig. S6. Viscous force and fluid drag force vary with distance from the bottom surface. As h increases, normalized Decreases, which means that the horizontal force due to rolling decreases. Normalized becomes larger and tends to 1, which is close to the fluid resistance of a spherical particle in the unbound fluid.   is the magnetic force between two magnetized particles, G is the gravity of particle 2, , are viscous force and shear force, respectively. Components of Moreover, G along the centerline provides the centripetal force.

Fig. S25. Analytical relationship between the VPNS upstream velocity and diameter. (A)
Plotted surfaces of VPNS upstream velocity versus its diameter and thickness. The maximum thickness of the VPNS reaches 4 μm, and the image shows that the thickness has almost no effect on the upstream velocity. (B) The powerful injection affects the correspondence between upstream velocity and the diameter of VPNS. The increase in vortex diameter induced the moderated increase in upstream velocity.