Ferrofluidic Manipulator: Theoretical Model for Single-Particle Velocity

Theoretical models have crucial importance in the design of algorithms for feedback control of robotic micromanipulation platforms. More importantly, theoretical models provide an understanding of the limits of each manipulation approach along with the identification of the key parameters that influence motion performance. Here, we provide a mathematical framework and numerical implementation of the velocity field for a single diamagnetic particle pinned at the air–ferrofluid interface when it is actuated by one or more inclined electromagnets in a system previously introduced as the ferrofluidic manipulator (Cenev et al., 2021). The theoretical model uses a magnetic dipole approximation for a dipole located between the tip of a solenoid and the end of the coil. The model reveals that the forces due to gravitation, the applied magnetic field, and the capillary action have decreasing contributions to the overall velocity field, respectively. The model assumes overdamped dynamics, and therefore, it is time-independent. The model is valid for an infinite number of electromagnetic solenoids. The theoretical predictions are in good agreement with the estimations from experimental data realized for one and two actuated solenoids.


Ferrofluidic Manipulator: Theoretical
Model for Single-Particle Velocity Zoran M. Cenev , Ville Havu, Jaakko V. I. Timonen , and Quan Zhou Abstract-Theoretical models have crucial importance in the design of algorithms for feedback control of robotic micromanipulation platforms.More importantly, theoretical models provide an understanding of the limits of each manipulation approach along with the identification of the key parameters that influence motion performance.Here, we provide a mathematical framework and numerical implementation of the velocity field for a single diamagnetic particle pinned at the air-ferrofluid interface when it is actuated by one or more inclined electromagnets in a system previously introduced as the ferrofluidic manipulator (Cenev et al., 2021).The theoretical model uses a magnetic dipole approximation for a dipole located between the tip of a solenoid and the end of the coil.The model reveals that the forces due to gravitation, the applied magnetic field, and the capillary action have decreasing contributions to the overall velocity field, respectively.The model assumes overdamped dynamics, and therefore, it is timeindependent.The model is valid for an infinite number of electromagnetic solenoids.The theoretical predictions are in good agreement with the estimations from experimental data realized for one and two actuated solenoids.
Index Terms-Microelectromechanical systems, smart material-based devices, soft robotics systems.

I. INTRODUCTION
A CTUATION and manipulation of a single particle or a particle collective at the fluid-fluid interface have received particular attention due to the rich phenomenology occurring at the interface.Techniques utilizing light-controlled interfacial flows were used for transporting oil droplets [2], [3], liquid marbles [4], and steel beads [5].Micropost perturbed interfaces [6], [7] and bubble-induced capillary forces [8] were used to transport objects by tuning the shape of the fluid-fluid interfacial contact area.
Studies on the manipulation of microscale matter at the airliquid interface by magnetic means can be also found in the literature.For example, Snezhko et al. [9], [10], [11] studied active and dynamic self-assembled structures, such as snakes, stars, and pulsating clusters induced by a vertical alternating magnetic field in a population of magnetic particles suspended on a fluid-fluid (air-water and oil-water) interfaces.Grosjean et al. [12], [13], [14] reported magnetocapillary self-assemblies that are able to move along the air-water interface when powered by oscillatory, uniform magnetic fields.Dkhil et al. [15] developed a closed-loop controller for positioning a magnetic particle at the air-water interface.Dong and Sitti [16] reported the formation and cooperative behavior of a swarm of magnetic agents at the air-water interface, and Wang et al. [17] demonstrated embedding information into structures of spinning magnetic microdisks at the air-water interface.Droplets and liquid marbles pinned at the air-paramagnetic liquid interface were actuated (or manipulated) by the application of a nonuniform magnetic field arising from a permanent magnet by Vialetto et al. [18].
Similarly, we have experimentally demonstrated and theoretically described the physical mechanism behind the motion and trapping of micro-and millimeter-sized particles pinned at the air-paramagnetic liquid interface in the presence of a nonuniform magnetic field [19].The paramagnetic liquid was an aqueous solution of paramagnetic salts such as manganese dichloride (MnCl 2 ) or holmium (III) nitrate pentahydrate (Ho(NO 3 )5H 2 O).This study examined the interaction among the particle, the liquid, and the magnet in an axisymmetric case.The motion of a particle on an air-paramagnetic liquid interface is a sum of the gravitational potential energy, the capillary energy from the interface deformation created by the nonuniform magnetic field, and the magnetic energy from the particle and the liquid.Consequently, a particle can be pulled toward and eventually trapped or pushed away from the magnet based on the particle magnetization, density, and contact angle with the liquid.
Other types of magnetic liquids similar to salt-based paramagnetic liquids are ferrofluids.Ferrofluids represent colloidal dispersions of magnetic nanoparticles within a carrier liquid with numerous applications recently reviewed in [20].Recent developments also involve ferrofluid droplet formation under varying magnetic energy landscapes [21], [22], bubble manipulation on ferrofluid-infused surfaces [23], and manipulation of ferrofluid droplets in a two-dimensional (2-D) environment [24] to name a few.
We have also built the ferrofluidic (FF) manipulator, which successfully was used in the automatic manipulation of individual polyethylene (PE) and polystyrene spherical particles, a silicon chip, and poppy and sesame seeds, all pinned at the air-ferrofluid interface [1].The operation principle of the FF manipulator is illustrated in Fig. 1(a).When a direct current is supplied to a solenoid, it generates magnetic field, which deforms the ferrofluid causing an interface deformation (a bump).The interface deformation acts as a slope onto which a single particle slides down.Note that the particle must have a density similar to the fluid's density to be pushed.Fig. 1(b) shows a 3-D cut-view of the FF manipulator with solenoids arranged in a circular pattern with 45°offset between each other and 45°offset with respect to the air-ferrofluid interface.Fig. 1(c) shows a photograph of the device.A water-like density particle such as PE is placed on the air-ferrofluid interface and can be manipulated in a quasi-2D workspace of about 8 mm, as shown in Fig. 1(d).There are eight solenoids in total, four short solenoids (denoting core neck length 1-2 mm) enumerated with odd numbers and four long solenoids (denoting core neck length 4-6 mm), enumerated with even numbers.In the scope of this work, we consider the particle to be at an initial distance from solenoid 1 along the Y-direction.
The motion of a single diamagnetic particle in the case of the FF manipulator cannot be simply described by the same energy minimization equations as in [19] or the velocity equation as in [18] since these equations are applicable to an axis-symmetric case only.Describing the motion of a single particle in the case of the FF manipulator necessities conversion to the Cartesian coordinates since the inclination of the solenoids should be considered.Such conversion is not a trivial problem because gravitational and capillary energies will depend on the nonlinear properties of the interface deformation and its derivatives.Moreover, the open-loop analysis in [1] characterized the velocity of a particle by studying the actuation velocity, i.e., the mean particle velocity during the first second.This characterization was sufficient to provide knowledge for building the control algorithm based on heuristically obtained and mathematically optimized parameters.However, the open-loop characterization does not quantify the overall velocity profile of a single particle and therefore limits the understanding of its velocity profile.
Here, we provide a theoretical framework that models the velocity of a single diamagnetic particle pinned at the air-ferrofluid interface when actuated by one or more inclined electromagnets as in the case of the FF manipulator.The validity of the theoretical model was verified against the velocity estimated from the experimental data for one and two actuated electromagnetic solenoids.The theoretical model is built on first principles and does not involve fitting parameters or extracting information from the experimental data of the motion of the particles.

II. EXPERIMENTAL DATA
A 550 ± 50 μm PE particle (WPMS-1.25,Cospheric LLC, USA) was carefully placed on the air-liquid interface of an inhouse developed ferrofluid [1] and brought into the workspace.The particle was displaced by supplying electrical power to solenoid 1 only or simultaneously to solenoids 1 and 8.
The experimental data were generated for the purpose of [1].In the context of this work, the raw data has been processed differently.The particle position data were loaded and smoothed with a smoothing spline for X and Y axes separately to minimize the fluctuations and consequently ease the computation of the derivative for each axis.The axial derivative essentially represents the velocity along each axis.Finally, the total velocity was calculated as a vector sum of the velocity components v x and v y .Each velocity component was also smoothed with a smoothing spline to minimize the fluctuations.
Here, we assume that the initial velocity was 0 mm/s.During the experimental work, however, a feedback control was used to bring the particle to a specific location, stop the controller for ∼1 s, and then actuate a solenoid or two to quantify the open-loop behavior.Also note that here we show the data where the sum of gravitational, magnetic, and capillary forces is nonzero; and we dismiss the data associated with a random movement that usually occurs at the end of each experimental trial.We also want to emphasize that we disregard the torque induced to the particle when the magnetic field gradient comes at an angle to the particle.We assume that the magnitude is negligibly small since we have not observed any meaningful rotation of a particle under manipulation in the experimental setup, the analyzed data for open-loop actuation, and the data from the automatic manipulation [1].Appendixes A and B comment on the data acquisition and processing.The first comment refers to the calibration of the data and the second to the limits of the estimation of the particle velocity from the experimental data.
We conducted two different types of experiments to examine the velocity profiles of the single nonmagnetic PE particle when actuated by a single solenoid.Each experimental trial was repeated five times.The first experiment examined the dependence of the particle velocity when the initial position of the particle was set to 2.23, 2.94, and 4.36 mm away from solenoid 1, respectively.In the absolute Y coordinate, the initial distances are -1.77,-1.06, and 0.36 mm.The velocity profiles estimated from the experimental data as a function of time are shown in Fig. 2(a).In all cases, the estimated velocity increased slightly in the beginning, and peaked at about 1 s time point, followed by a monotonic decay.The second experiment examined the particle velocity versus the actuation current with actuation currents of 0.95, 1.19, and 1.43 A applied to solenoid 1.The estimated velocities from the experimental data are plotted in Fig. 2(b).Similar to the results with varying initial distances, the estimated velocity initially slightly increased, and peaked at about 1 s, followed by a monotonic decay.
Finally, a third experiment studied the dependence of the particle velocity when it was simultaneously actuated by one short (solenoid 1) and one long solenoid (solenoid 8) at 1.43 A excitation current in each solenoid.As evident from Fig. 2(c), and similar to the first two experiments, the estimated velocity from the experimental data initially slightly increased, and peaked at about 1 s timepoint, following a monotonic descent.

III. MATHEMATICAL MODEL
This mathematical model encompasses the magnetic, the capillary, and the gravitational energy of a particle floating at the surface of a ferrofluid.In the first step, we discuss how the interface deformation arises from the nonuniform magnetic field generated by an inclined dipole at a distance z above the interface.Then we consider a spherical particle trapped at this interface, followed by performing force decomposition.Finally, we determine the velocity.We consider the general case where the magnetic field has no symmetry and arises from several solenoids at arbitrary positions and orientations.

A. Magnetic Stress Tensor
The solenoids, numbered by i = 1, . . ., n, are at positions R i = X i e x + Y i e y + Z i e z , in respect to the origin O(0, 0, 0) where e i are the unit vectors and Z i is the height of the center of the magnetic dipole above the ferrofluid interface.X i and Y i are the in-plane coordinates of the magnetic dipole of the ith solenoid.The magnetic dipoles are oriented according to (1) Then, the magnetic field is the sum of the dipolar fields (2) where μ 0 is the magnetic permeability of vacuum and r = xe x + ye y + ze z is the position vector of interest (e.g., the center of gravity of a spherical particle) with respect to the origin in Cartesian coordinates x, y, z.The difference r = r − R i is the vector from the magnetic dipole to the position of the particle.The scalar product of two vectors m and r is with the magnitude of the vector The interfacial magnetic pressure Π is defined through the square of the magnetic field z , evaluated at the interface z = 0, same as in [19].Hence, we obtain where μ ff and μ 0 are the magnetic permeabilities of the ferrofluid and free space, respectively.

B. Young-Laplace Equation
The Young-Laplace equation is a mathematical description of the interfacial capillary pressure difference between two static fluids.Therefore, the Young-Laplace equation for the FF manipulator can be described as follows: where γ is the surface tension of the ferrofluid, u(x, y) is the interface deformation with H(x, y) being its mean curvature, g is the gravitational acceleration constant, ρ ff is the density of the ferrofluid, and Π(x, y) is the magnetic pressure subjected to the air-ferrofluid interface defined in (5).Since the pressure Π is given as an explicit function of the coordinates x and y , it is convenient to express the gradient of u as ∇u (x, y) = e x u x (x, y) + e y u y (x, y) (7) such that u x = ∂u/∂x and u y = ∂u/∂y.Here, we denote the gradient operator with ∇.Similarly, let us define u xx = ∂ 2 u/∂x 2 and u yy = ∂ 2 u/∂y 2 , then the mean curvature has the following form: with its gradient ∇H (x, y) = e x H x (x, y) + e y H y (x, y) ( where H x = ∂H/∂x and H y = ∂H/∂y.The analytical expressions for the solutions of the interface deformation u(x, y) and its derivatives are provided in Appendix C.

C. Forces Acting on a Particle
The effective mass of a spherical particle floating on the airferrofluid can be expressed as in [19], [25] where ρ p and ρ ff denote the density of the particle and the ferrofluid, respectively, and V p is the volume of the particle.The immersion volume V imm represents the volume of the replaced ferrofluid by the particle and can be expressed as [19], [25] with a p being the particle radius, and θ denoting the wetting angle between the particle and the ferrofluid.When a nonuniform magnetic field is subjected to a particle floating on the air-ferrofluid interface, the particle experiences the following forces: the magnetic force exerted on the particle F M , the gravitational force F G from the deformed airferrofluid interface, the capillary force F C at the interface, and the drag force opposing the motion of the particle F D .Hence, we obtain The magnetic force has the following form: where χ ff = μ ff μ 0 − 1 and χ p denote the magnetic susceptibility of the ferrofluid and the particle.The magnetic force has a positive sign because here the ferrofluid has much higher magnetic susceptibility than the particle (i.e., χ ff χ p ).The gravitational force along the air-ferrofluid interface has the following form: Similarly, the capillary force reads where l = γ/gρ ff is the capillary length of the ferrofluid.Both, the gravitational and the capillary force have a positive sign.In the former case, the reason is that the effective mass m ef f is positive due to the similarity of the densities of the particle and the ferrofluid, and in the latter case due to the positive curvature of the deformed interface.
When the particle is moving on the interface, it will experience a drag opposing the movement.This drag force could be approximated as Here, f d = 2Cπηr 0 represents the friction factor and v p is the velocity of the particle; C is a drag constant ranging from 1 to 3, η is the viscosity of the ferrofluid, and r 0 = a p sin θ is the radius of the contact line.
The expanded form of the force could be rewritten as or

D. Velocity Exerted on a Single Particle
According to Newton's second law, the total force acting on the particle exerts acceleration dv p /dt; hence, we obtain However, assuming that the motion of the particle is overdamped (usually 1 or maximum 2 body lengths per second), the acceleration dv p /dt is equated to 0. Therefore, we obtain the velocity of the particle (20) Note that when a particle resides at the air-liquid interface it exhibits a random movement due to locally induced airflow (e.g., from movements of the operator nearby) as well as thermocapillary convective flows resulting from the heat exposure (e.g., from illuminating source and variations in the ambient temperature) of the interface.These contributions are not accounted in the presented model.

IV. NUMERICAL IMPLEMENTATION
COMSOL Multiphysics 6.0 (COMSOL Group, Sweden) was used for solving the partial differential equation (PDE) (6) for the interface deformation u and its first and second derivatives for obtaining the mean curvature H (8). The data for u and its derivatives were transferred to MATLAB (v2021) where the gradient of the curvature was computed.Finally, the velocity was computed with (20) and the resulting velocity profiles were compared to the estimated velocity profiles from the experimental data.The density (1071 kg/m 3 ), viscosity (0.9 mPs), and surface tension (74.75 mN/m) of the ferrofluid are taken from [1].The magnetic susceptibility (0.016) was measured with a vibrating sample magnetometer mode in Quantum Design PPMS equipment (Quantum Design, USA).The contact angle of the particle with the ferrofluid was estimated to ∼75°.

A. Calibration to Experimental Data
The magnitude of the magnetic flux density (or magnetic field) was numerically computed by a finite-element model (FEM) implemented in COMSOL and experimentally measured in [1].The magnetic dipole position and strength were adjusted in a 2-D axisymmetric study against the FEM data for two different cut lines [Fig.3(a)], one along the easy axis of the solenoid [Fig.3(b)] and one off-centered line resembling the position of the air-ferrofluid interface [Fig.3(c)].For an actuation current of 1.43 A, the identified magnetic dipole strength was 0.0021 Am 2 for a short solenoid at a location of 15 mm from the tip of the solenoid core.For a long solenoid, the magnetic dipole strength was 0.0012 Am 2 at the location of 5 mm from the tip.The position of magnetic dipoles in global Cartesian coordinates (i.e., with respect to the coordinate system of the FF manipulator) were (0, -5.41, 1.50) mm and (3.82, -3.82, 1.50) mm for the short solenoid and the long solenoid, respectively.

B. Solving for the Young-Laplace Equation and Obtaining the Interface Deformation u and Its Derivatives
A 2-D stationary COMSOL model was defined for solving a coefficient defined PDE with diffusion coefficient c = γ, absorption coefficients a = ρ ff g, and source term f = Π.The geometry consisted of two circles, one with a diameter of 8 mm, approximating the dimensions of the workspace; and the other one with a diameter of 57 mm, approximating the size of the Petri dish containing the ferrofluid.The solution for u was numerically obtained.The first and second derivatives of u, the mean curvature H, and the gradient of the mean curvature, ∇H, were computed by numerical differentiation.Finally, the velocity was computed using (20).

A. Numerical Results for the Interface Deformation and the Mean Curvature
Fig. 4(a) illustrates the experimental procedure with varying initial distances whose results are shown in Fig. 2(a).The particle is initially located at a particular initial distance (denoted with initDist on the figure) from the actuated solenoid (solenoid 1).The particle is displaced, and its position change has been recorded.Fig. 4(b) shows a surface plot of the computed interface deformation u when a short solenoid (solenoid 1) was actuated with 1.43 A. The interface deformation along the Y-axis is shown in Fig. 4(c).The peak of the interface deformation is ∼35 μm.This magnitude of the interface deformation has a reasonable magnitude even though a precise measurement is lacking.The interface deformation is a bump that monotonically rises, reaches a peak, and decays.The numerically computed peak of the interface deformation is located at -498 mm from the origin; hence, it is outside of the workspace of the FF manipulator and at about 0.4 mm offset distance from the location of the magnetic dipole.
The first derivative of u along X-and Y-axes and the vector field of the negative gradient of u are shown in Fig. 8(a) and (b), respectively.Note that the gravitational contribution in the velocity (14) depends on the ∇u.The second derivative of u along X-and Y-axis is shown in Fig. 8(c)(i and ii).The mean curvature along the Y-axis is shown in Fig. 4(d).The mean curvature contains two local peaks and one valley.The left peak is 0.818 m -2 and the right is 0.760 m -2 .The difference of ∼7% signifies asymmetricity caused by the 45°inclination of the solenoid with respect to the air-ferrofluid interface.The 3-D surface plot of the mean curvature is shown in Fig. 8(d).The gradient of the mean curvature along the X and Y axes is shown in Fig. 8(e)(i and ii).The simulated results show that similarly to the gravitational force the capillary force points away from the peak of the deformation adding to the velocity of the particle.

B. Velocity Profile. Comparison of the Particle Velocity Estimations From the Theoretical Model and From the Experimental Data
The surface plot of the particle velocity is shown in Fig. 5(a) and has a cylinder cone shape.The cylinder cone shape features a steep conical hill and a hollow cylinder, or crater, at the peak of the conical hill, resembling the shape of a volcano.The peak surface (red data points) has a tilt toward the positive Y-axis.The surface plots of v x and v y are shown in Fig. 8(f).The vector field of the particle velocity is shown in Fig. 5(b)(i and ii), and it indicates that the particle velocity dominantly is directed away from the peak of the interface deformation.The velocity field also indicated symmetricity along the Y-axis.Fig. 5(c) shows the velocity profile along the Y-axis.The velocity has two peaks and one valley close to the location of the origin of the magnetic dipole and monotonically decreases at the tails.The velocity profile in Fig. 5(c) predicts a local minimum where a PE particle could be trapped.Indeed, this scenario is possible, but such a scenario would be energetically very unfavorable.One analogy of this scenario could be to balance a ball on top of another ball, which is theoretically possible, but practically extremely rare to occur.
Fig. 5(d) shows the individual contributions to the particle velocity.The most dominant contribution is from the gravitational force, followed by the magnetic force contribution.The least dominant is the capillary force (or curvature) contribution.
Finally, Fig. 5(e) shows the comparison between the theoretically estimated velocity and the estimated velocity from the experimental data.The data is for a short solenoid (solenoid 1) actuated at 1.43 A. The theoretical and experimental estimated particle velocities follow the same trend and order of magnitude.The discrepancy in the magnitude ranges from ∼0% (overlapping curves) to a maximum of 40%.The root-mean-square errors are 170.7,156.6, and 191.2 μm/s, for initial distances of 2.23, 2.94, and 4.36 mm, respectively.Fig. 6 shows the particle velocity for different actuation currents.Fig. 6(a) illustrates the theoretically estimated particle velocities for the three actuation currents 1.43, 1.19, and 0.95 A at the initial distance of 2.48 mm.The estimated velocity from the experimental data in Fig. 2(b) depicts the velocity with respect to time, whereas Fig. 6 depicts the velocity with respect to the position of the particle.Fig. 6(b) provides a comparison between theoretical velocity and estimated particle velocity from the experimental data for the three actuation currents.The theoretical velocities match the decaying trend very well.However, the velocity magnitude features an error rate of factor 1.3 up to factor 2 at the initial phase of the movement.The largest error is obtained for the highest actuation current (1.43 A), and the smallest error for the lowest actuation current (0.95 A).The theoretical predictions fit the estimation from the experimental data better for distances further away from the magnetic dipole.The root-mean-square errors are 59.1, 165.8, and 332.5 μm/s for actuation currents of 1.43, 1.19, and 0.95 A, respectively.Fig. 7 shows the particle velocity for the simultaneous actuation of two solenoids (solenoids 1 and 8).Fig. 7(a) illustrates the actuation scenario and Fig. 7(b) shows the surface plot of the interface deformation.The maximum magnitude is ∼45 μm and it is located below the magnetic dipole of solenoid 1.A second local peak does also exist under the magnetic dipole of the solenoid 8. Fig. 7(c) shows the surface plot of the velocity.The surface plots of v x and v y are shown in Fig. 8(g).The surface plot of the particle velocity shows two craters and one dominant peak.Similar to the case of a single solenoid, the particle velocity is largely directed away from the magnetic dipoles such that the stronger decay is present for the stronger magnetic dipole (solenoid 1).The vector field of the particle velocity is shown in Fig. 7(d)(i and ii).Note that the arrows of the vector on the left-hand side of the workspace point toward the north-west direction, but the ones on the right-hand side point toward the north.Finally, Fig. 7(e) shows the comparison between the theoretically estimated particle velocity and estimated velocity from the experimental data when a short solenoid (solenoid 1) and a long solenoid (solenoid 8) have been actuated with  The reasons for the mismatch of the theoretical prediction [in Figs.5(e), 6(b), and 7(e)] and the estimations from the experimental data could be attributed to experimental errors due to the constant presence of thermal noise (e.g., thermocapillary convective flows) and its increase in dominance the further the particle goes from the actuated solenoid(s); experimental errors due to the ongoing evaporation and change of the magnetic susceptibility of the ferrofluid; the estimation errors arising from the image acquisition frequency, image processing, and data analysis including two subsequent averaging steps and fitting procedures, see Appendixes A and B; inaccuracy arising from the assumed constant wetting contact angle and triple contact line of the PE particle while residing on the deformed interface that lead to the estimation of the deformed interface u; and the assumption of overdamped dynamics (disregarding the initial acceleration) in the theoretical model.
One additional reason for the discrepancy in Fig. 7(e) is the inaccuracies in the location of the magnetic dipole of the long solenoid (solenoid 8).

VI. SUMMARY AND CONCLUSION
This article reported a theoretical model and implementation of the velocity field for a single diamagnetic particle pinned at the air-ferrofluid interface.The theoretical model is valid for any number of electromagnetic solenoids, even though we have demonstrated its validity only with one particle and up to two solenoids.Our theoretical model correctly estimates the magnitude and the trend of the interface deformation as well as the magnitude and trend of the velocity.The theoretical velocity features a root-mean-square error of ∼15%, 18%, and 30% of the maximum velocity estimated from the experimental data in each case for the variable initial distances.For the case of actuation with variable current, we obtained a root-mean-square error of ∼11%, 24%, and 43% of the maximum velocity estimated from the experimental data.For the case of two-solenoid actuation, the root-mean-square error was ∼25%.
We have also theoretically analyzed the individual components in the velocity profile and not just matched the experimental results per se.We found that all forces (gravitational, magnetic, and capillary) are contributing to the overall velocity such that the biggest dominance comes from the gravitational force, then from the magnetic force, and the least dominant contribution from the capillary force.We also found that in the magnetic dipole approximation, the location of the dipole lies somewhere between the tip of the solenoid and the end of the coil.The theoretical model is time-independent, i.e., assumes overdamped dynamics so the initial acceleration is not captured.The changing magnetic field in the workspace of the FF manipulator induces currents in the solenoids, which further generates magnetic field from the induced currents.However, we anticipate that the generated magnetic field from the induced currents will have negligible magnitude than the magnetic field generated by the supplied direct currents, thus such contributions are not accounted for.
The presented model provides an understanding of the underlying physics of the manipulation principle in the ferrofluid manipulator.It also can serve as a means for the development of control algorithms.The particle velocity is based on the constant parameters for the particle and the ferrofluid and only the magnetic field remains to be the only varying quantity, even though the magnetic field gradient ∇(B 2 ) is the one present in (20).Note that the gradient of the deformation interface ∇u and the curvature of the deformation interface H are still dependent on the magnetic field since (A4) and (A6) show that the ∇u and H are proportional to the magnetic pressure Π, which in turn is proportional to the square of the magnitude of the magnetic field B. One infers that the velocity of the particle is strictly proportional of the applied magnetic field, which can be controlled by modulating the (direct) current supplied to each solenoid.
The theoretical model and the manipulation principle of the ferrofluid manipulator are fundamentally different from the magnetic actuation principles in the established wireless/untethered magnetic methods in robotics [26].The proposed model may facilitate the development of other theoretical models and provide design considerations for future micromanipulation platforms that will use ferrofluids for indirect manipulation at the air-liquid interface.
The utilization of ferrofluids to manipulate nonmagnetic media (diamagnetic or weakly paramagnetic), specifically on the air-ferrofluid interface, could be extended in several new avenues.Distinct patterns may be generated if a multiparticle system is considered.Such a system could be interesting from the perspective of active matter physics, but also from a robotic perspective because feedback control can be used for phase shifting or motion control of entire clusters [27], [28].An alternative approach could be the utilization of living matter (cells/microorganisms) or oil droplets instead of synthetic particles.

A. Comment on Experimental Data Recalibration
The experimental data has been obtained for characterizing the actuation velocity (velocity during the first second) in an open loop regime for demonstrating dependencies to particle velocity from the solenoid to particle distance, actuation current, and the number of actuated solenoids in [1].The pixel-to-mm calibration of the camera was performed by capturing an image of a ruler with a resolution of 100 μm.Here, we recalibrate the data by taking the particle diameter (datasheet value: 550 μm) as ground truth, to which we adjust the pixel-to-mm ratio.

B. Comment on Estimation of Particle Velocity From the Experimental Data
The camera used for recording the particle motion was used for data acquisition and feedback control of the particle simultaneously.The acquisition frequency was set to 10 frames per second.The feedback control was used to position the particle in the initial position and then actuate it in order to examine its open-loop behavior.
We estimate the velocity from the position data of the particle.The particle position is obtained by using an object tracking algorithm followed by a smoothing step to the position data (first averaging step).A derivation is applied to the smoothed position data to obtain the particle velocity.The obtained velocity data is also smoothed (second averaging step) to obtain a reasonable estimate since the fluctuating velocity is full of noise.Note that we are applying two smoothing steps in the processing of the data.The velocity obtained from the experimental data has an associated error, which is difficult to estimate.Hence, it should not be taken as ground truth.

Fig. 1 .
Fig. 1.Principle of operation of the FF manipulator (adapted and modified from [1]).(a) Diamagnetic particle with water-like density is pinned at the air-ferrofluid interface and is pushed away from a solenoid when it is actuated (electric current supplied to it).The actuation by single solenoid observed by (i) front view; (ii) side view; and (iii) top view.(b) Cut-view of a computer-aided design (CAD) rendered image of the FF manipulator's mechanical assembly consisting of a holding block enclosing the eight solenoids located on top of a Petri dish.The Petri dish contains a ferrofluid and a particle pinned at the air-ferrofluid interface.(c) Photograph of the FF manipulator showing the holding block with the solenoids and the dish with the ferrofluid.(d) CAD image of the top-view of the FF manipulator with the eight solenoids enumerated from 1 to 8 and a single particle (e.g., PE particle) in the workspace at an initial distance from Solenoid 1 along the Y-axis.

Fig. 2 .
Fig. 2. Estimated particle velocity (µm/s) from experimental data over time (s) at different actuation conditions.(a) Actuation of a single particle at three different initial distances 2.23, 2.94, and 4.36 mm of the particle from solenoid 1 at actuation current of 1.43 A. (b) Actuation of a single particle at three different actuation currents 0.95, 1.19, and 1.43 A of solenoid 1 for the initial distance of 2.48 mm.(c) Simultaneous actuation of a single particle with short + long (SL) solenoids (solenoids 1 and 8) for actuation current of 1.43 A in each solenoid and initial distance of 4.16 mm from solenoid 1.In all cases, each actuation experiment has been repeated for five times.The error bars denote 95% confidence interval.

Fig. 3 .
Fig. 3. Numerical implementation and calibration of the magnetic dipole in a short solenoid.(a) Short solenoid geometry for axisymmetric numerical study and definition of cut lines.(b) Calibration of the magnetic dipole against the FEM model of the solenoid line along the axis of symmetry (cut line * ).(c) Calibration of the magnetic dipole against the FEM model of the solenoid line along the off-centered line (cut line * * ).

Fig. 4 .
Fig. 4. Numerical results for the interface deformation and mean curvature when the solenoid 1 (short solenoid) is actuated at 1.43 A. (a) Illustration of the top-view of the FF manipulator when the solenoid 1 is actuated.(b) 2-D surface deformation u (in meters) for the given actuation input.(c) Interface deformation u and (d) mean curvature H along the ordinate.

Fig. 5 .
Fig. 5. Numerical results for the particle velocity when the solenoid 1 (short solenoid) is actuated at 1.43 A at three different initial distances.(a) Magnitude of the theoretical particle velocity.(b) Vector field of the theoretical particle velocity (i) around the workspace, and (ii) close-in to the origin of the interface deformation.The red circle denotes the workspace.(c) Theoretical particle velocity along Y-axis across the whole dish.(d) Individual contributions to the theoretical particle velocity.(e) Comparison of the theoretical versus estimated particle velocity from the experimental data.

Fig. 6 .
Fig. 6. Results for particle velocity at different actuation currents.(a) Theoretically estimated particle velocity for actuation currents 1.43, 1.19, and 0.95 A for the initial distance of 2.48 mm.(b) Comparison between theoretical velocity and velocity estimated from experimental data for actuation currents of (i) 1.43 A, (ii) 1.19 A, and (iii) 0.95 A. Theoretical velocities are represented with circle/triangle/asterisk-marked line and estimated velocities from experimental data with lines with error bars.1.43 A per solenoid.The largest error (∼60%) is obtained at the beginning of the actuation and with time (or position change) the error is decreased up to a point to meet the estimates from the experimental data.The trend in both cases is matching.The root-mean-square error is 187.3 μm/s.The reasons for the mismatch of the theoretical prediction [in Figs.5(e), 6(b), and 7(e)] and the estimations from the experimental data could be attributed to experimental errors due to the constant presence of thermal noise (e.g., thermocapillary convective flows) and its increase in dominance the further the particle goes from the actuated solenoid(s); experimental errors due to the ongoing evaporation and change of the magnetic susceptibility of the ferrofluid; the estimation errors arising from the image acquisition frequency, image processing, and data analysis including two subsequent averaging steps and fitting procedures, see Appendixes A and B; inaccuracy arising from the assumed constant wetting contact angle and triple contact line of the PE particle while residing on the deformed interface that lead to the estimation of the deformed interface u; and the assumption of overdamped dynamics (disregarding the initial acceleration) in the theoretical model.One additional reason for the discrepancy in Fig.7(e) is the inaccuracies in the location of the magnetic dipole of the long solenoid (solenoid 8).

Fig. 7 .
Fig. 7. Two-solenoid actuation.(a) Illustration (to scale) of the top-view of the FF manipulator when the solenoids 1 and 8 are actuated with 1.43 A. (b) Surface deformation in meters.(c) Surface plot of the magnitude of the particle velocity.(d) Vector field of the particle velocity (i) around the workspace, and (ii) close-in to the origin of the interface deformation.The red circle denotes the workspace.(e) Comparison of the theoretical particle velocity versus the particle velocity estimated from the experimental data.

Fig. 8 .
Fig. 8. Intermediate numerical results.Actuation of solenoid 1 (short solenoid) at 1.43 A. (a) Gradient of the interface deformation along (i) X-axis and (ii) Y-axis.(b) Vector field of the negative gradient of the interface deformation.The red semicircle represents part of the workspace.(c) Second derivative of the interface deformation along (i) X-axis and (ii) Y-axis.(d) Surface plot of the mean curvature H. (e) Gradient of the mean curvature along (i) X-axis and (ii) Y-axis.(f) Surface plot of the particle velocity along (i) X-axis and (b) Y-axis for actuation of solenoid 1 (short solenoid) at 1.43 A. (g) Surface plot of the particle velocity along (a) X-axis and (b) Y-axis for simultaneous actuation of solenoids 1 and 8 at actuation current of 1.43 A per solenoid.