Dynamics of a bottom-heavy Janus particle near a wall under shear flow

Abstract

In this study, Brownian dynamics simulations are implemented to investigate the motion of a bottom-heavy Janus particle near a wall under varying shear flow conditions and at small Péclet (Pe) numbers. The stochastic motion of the Janus particle impacted by surface forces is described using a set of coupled Langevin equations that takes into account the Janus particle orientation. Interactions arising from surface potentials are found to depend on the separation distance between the Janus particle and the wall, the properties of the surfaces involved, and the thickness of the Janus particle cap. When shear flow is introduced in the system, the dynamical behavior of the Janus particle is also governed by the strain rate. Furthermore, the effect of friction on the dynamical behavior of the Janus particle under shear flow is investigated and reveals that the rotational motion of the Janus particle slows down slightly when the particle is close to the surface. In summary, we demonstrate the ability to utilize Brownian dynamics simulations to capture the rich dynamical behavior of a bottom-heavy Janus particle near a wall and under a range of shear flow conditions, cap thicknesses, and surface charges.

---

1 Introduction

Besides their ubiquity as amphiphilic colloids in nature, Janus colloids have a wide variety of applications in areas such as biomedicine,^1–3 energy storage,⁴ catalysis,⁵ and environmental science.⁶ In many of these anticipated applications Janus particles will encounter various types of boundaries such as fluid/fluid interfaces⁷ and solid walls of blood vessels, channels, and porous structures. As such, understanding and predicting individual particle–wall interactions is important. For example, particle–wall interactions can determine drug delivery rates,⁸ guide particle swimming,^9–11 result in scale formation in pipes,¹² or change the flow in porous media.¹³

There is a long scientific history regarding the description of particle–wall interactions.^14–19 Most of the literature addresses spheres or ellipsoids near boundaries with a few studies focusing on spherical passive Janus particles discussed below. Experimentally, three techniques, atomic force microscopy (AFM), surface force apparatus (SFA), and total internal reflection microscopy (TIRM), have been used to measure particle–wall interactions. AFM and SFA require the immobilization of the probe particle,^20–22 whereas TIRM does not. TIRM²³ is an experimental light scattering technique that indirectly measures particle position above a surface as a function of time resulting in a probability density function (PDF), p(h), in which h is the height of the particle surface from the wall. It then infers the particle–wall potential, U(h), from analysis of the PDF. The particle–wall interaction potential and scattering morphology become more complex functions when Janus particles are involved due to their non-uniform surface charge and potentially heavier cap that influence particle height and orientation in a complex manner. In a recent perspective article, Yan and Wirth²⁴ introduced Scattering Morphology Resolved TIRM (SMR-TIRM) and demonstrated that assessment of anisotropic particle interactions near boundaries with SMR-TIRM is more accurate when combined with Brownian dynamics (BD) simulations. Using BD simulations, Rashidi et al.^25,26 explored the dynamics of a polystyrene Janus particle with a non-uniform surface potential with and without a gold cap in the vicinity of a negatively charged wall and successfully derived the 3D potential energy surface of the particle that can be used in the interpretation of TIRM measurements. In a follow up study, Rajupet et al.²⁷ used the surface element integration (SEI) method to elucidate the impact of van der Waals forces on Janus particle–wall interactions and found that at short distances (relevant in systems of high ionic strength) van der Waals interactions and cap nonuniformity significantly impact the potential landscape experienced by the Janus particle. More recently, the focus has been on understanding the behavior of Janus particles in channels under shear flow, which further complicates the interpretation of TIRM data. Therefore, the development of robust theoretical models and their numerical analysis are a central tool to fully capture the physics of Janus particles, as models are needed to interpret experimental results and, at the same time, might predict novel behaviors to be in turn verified experimentally.

A number of simulation studies have investigated the dynamic behavior of Janus particle suspensions in channels and near solid boundaries with focus on particle–particle and particle–wall interactions under shear flow. Using BD simulations, Mohammadi et al.²⁸ investigated the binding kinetics of Janus particles under shear flow and determined the coagulation rate of a particle pair. Kobayashi et al.²⁹ studied the rheology of Janus nanoparticle suspensions in nanotubes of varying surface chemistry with dissipative particle dynamics and found that flow properties strongly depended on the structure of the colloid phase and the wall chemistry. Nikoubashman et al.^30,31 employed hybrid molecular dynamics simulations to investigate the shear-induced break up of Janus particle micelles within a slit-like channel under shear flow and found that the shear flow can lead to enhancement of micelles and the formation of bigger aggregates some of which show enhanced stability due to their cluster symmetry. In a follow up study,³² they found that increasing particle–wall interactions causes adhesion and cluster formation on the channel walls. DeLaCruz-Araujo et al.³³ showed using BD simulations that exposure of micellar, vesicular, worm-like and lamellar Janus particle aggregates to shear flow induces rearrangement, deformation, and break-up of the aggregates that is more easily controlled by the Peclet number (Pe) than the interaction potential. In a second study,³⁴ they investigated the switch between parallel and perpendicular alignment of lamellar Janus particle aggregates by shear flow and developed a simple scaling argument based on the torque balance on a single Janus particle pair.

Other studies have focused on a single passive Janus particle near a wall. For example, Ramachandran et al.³⁵ used a finite element method to study the dynamics and rheology of a dilute suspension of slip-stick Janus particles in creeping flow. They showed that depending on the ratio of the slip length to the particle radius the rotational motion of the particle varies. In addition, the translational motion of the particle is predicted to be impacted by the initial orientation of the particle. Furthermore, the suspension exhibits non-Newtonian behavior for a specific volume fraction of slip-stick Janus spheres. Troufa et al.³⁶ numerically simulated the dynamics of a slip-stick Janus sphere within a Newtonian fluid confined in a cylindrical micro-channel under Poiseuille flow and found two behavior regimes where the particle either shows periodic oscillations near the wall or migration towards the channel center if near the centerline.

Most recent studies have primarily focused on the steering of active particles near walls,^37–41 though open questions regarding the interactions of passive Janus particles with walls under shear flow persist. Additionally, Janus particles with thick caps (50–100 nm) have been explored, where the cap volume is required to attain the field interaction strength needed to control the Janus particle behavior.^42–47 Understanding how the interplay of the particle's and substrate's properties influences the particles height, rotation speed, and the frequency with which its surface interacts with the substrate is essential for advancing applications in self-assembly, catalysis, drug delivery, and biosensing.

Here, a BD model is introduced that captures the interaction of a bottom-heavy Janus particle exhibiting two distinctly charged surfaces with a charged wall with and without shear flow. The BD model is based on the Hogg, Healy and Furstenau (HHF) sphere–sphere interaction model⁴⁸ modified to apply to a sphere-infinite plate interaction by setting one of the sphere radii to infinity. Note, the model is limited to thin double layers (i.e., κR ≫ 1). The BD model is used to capture the non-trivial interplay between surface properties of the particle, i.e., surface potential and cap thickness, and the nearby wall in the presence of shear flow, where surface forces and hydrodynamic interactions (using Goldman et al.'s dimensionless frictional coefficients)⁴⁹ govern the dynamics and determine the particle height.

2 Simulation model

A hard sphere of radius R with two different hemispheres (1,2) serves as a model particle for a system in which each face is uniformly (assuming constant surface potential), but independently charged, creating the Janus particle (Fig. 1). The Janus particle is suspended in an aqueous solution of ionic strength I, which throughout this study is considered to be deionized (DI) water (Millipore, resistivity 18.2 MΩ cm and viscosity μ = 1 × 10⁻³ Pa s at 25 °C, I = 1 × 10⁻⁶ M). The Janus particle is bounded by a wall, i.e., a substrate made of SiO₂, which is uniformly charged. The normal to the particle (n_p), that passes through the center of the particle and through the cap's center of mass (COM_cap) is employed to characterize the rotational behavior of the particle, shown in Fig. 1. Finally, the Janus particle is subjected to a shear flow with varying strain rates [small gamma, Greek, dot above].

Fig. 1 Schematic representation of a bottom-heavy Janus particle capped on face 2 with a metal of thickness δ and a cap center of mass COM_cap, bounded with a wall w and subject to shear flow of strain rate [small gamma, Greek, dot above]. θ, indicates the particle orientation with θ = 0 indicating a cap/face 2 up orientation, while θ = π indicates a cap/face 2 down orientation. Other relevant parameters are also depicted (see the text).

Anisotropic interactions between the Janus particle and the wall are modeled by utilizing the total interaction energy, U_tot, as the sum of electrostatic double layer (EDL) interactions and the gravitational interaction. Hence, U_tot is given by eqn (1):^50,51

U_tot(z, θ, ψ₁, ψ₂, ψ_w) = U_EDL(z, θ, ψ₁, ψ₂, ψ_w) + U_gravity(z), (1)

where z denotes the height of the Janus particle, i.e., the distance of the Janus particle from the wall, θ is defined as the angle between the normal to the particle (i.e., unit vector n_p) and the normal to the wall (i.e., unit vector n_w), and ψ₁, ψ₂, and ψ_w are the surface potentials of the two faces and the wall, respectively. The contribution from the EDL interaction between a uniformly charged particle and a uniformly charged wall (based on the modified HHF model) is given by eqn (2):where A is introduced as a dimensionless parameter, eqn (3a), and w_EDL(z̄, ϕ_p, ϕ_w) as the dimensionless EDL interaction, eqn (3b):

A: = πRε₀ψ₀²β, (3a)

with z̄ = κz, where z is the height of the particle non-dimensionalized by measuring the length scale in units of Debye screening length, κ⁻¹. In eqn (3a), ψ₀ = 1/qβ is the thermal electrostatic potential used to non-dimensionalize the wall and surface potentials as ϕ_i = ψ_i/ψ₀ with i = w and p, respectively, where q denotes the electron charge (in absolute value) and β = 1/(k_bT), where k_bT is the thermal energy. Note the same notation is used in the ESI,† Section S1, which provides additional information on the derivation. Use of the modified HHF model for the sphere-infinite plate interaction potential assumes a small constant surface potential (<25 mV), i.e., the Debye Hückel approximation, and small double layer thickness compared to the particle size, i.e., the Derjaguin approximation. HHF showed in their work that the model can be extended to surface potentials of 50–60 mV.⁴⁸

The Janus particle assumes various configurations in proximity to the wall and each hemisphere of the Janus particle contributes its own particle–wall EDL interaction to U_EDL. The anisotropic interaction is a function of the Janus particle orientation with respect to the wall as specified by θ (see Fig. 1). The orientation of the cap is incorporated into the model using the orientation-dependent factors yielding U_tot for a Janus particle interacting with a wall, eqn (4):

The governing equations that describe both the translational and rotational motion of the Janus particle near the wall, known as Langevin equations, are:

• Translation of particle center, O, along the x direction:

• Translation along the z direction:

• Rotation around the y axis (passing through particle center, O):

In the above equations, C_d = 6πμR and ξ = 8πμR³ are the translational and rotational drag coefficients, respectively, and f_x and f_θ are Goldman et al.'s⁴⁹ dimensionless frictional coefficients for translational and rotational motion. Note that expressions for C_d and ξ represent a perfect sphere and that the contribution from the cap to the particle shape is neglected here. V_x(z_c) = [small gamma, Greek, dot above]z_c denotes the velocity of the external shear flow that is being exerted on the position of the center of the Janus particle, O (see Fig. 1), by the externally applied strain rate, [small gamma, Greek, dot above]. [small script l]_xη_x(t) and [small script l]_zη_z(t) are the Brownian forces in the x and z directions, respectively, while Eη_r(t) is the rotational Brownian force. F_g = m*g and N_g = −Ξ*sin(θ) are the gravitational force and torque, respectively, where m* and Ξ* correspond to the effective mass of the particle and effective weight of the cap, respectively. Note that in the absence of shear flow, [small gamma, Greek, dot above] = 0, one clearly has V_x = 0. Finally, F_EDL and N_EDL are the contributions from the electric double layer interaction and defined as:

The conservative wall–particle interaction forces are computed through the wall–particle interaction energy as F_tot = −∇U_tot, which includes the electrostatic double layer interaction and gravitational potential energy. Eqn (5)–(9) describe the most general system under analysis in this work and is referred to as model IV. In particular, they account for different surface potentials of the two sides of the Janus particle, an asymmetry in the mass distribution as implied by the cap of thickness δ, the torque induced by the shear flow, and the hydrodynamic friction. Note that an imbalance in the surface potentials and the mass asymmetry induce a force, F_EDL and F_g, respectively, as well as a torque, N_EDL and N_g, with the forces acting also on symmetric systems while the torques are zero.

Before considering model IV, we first consider two simpler models of Janus particles that assess the role of various parameters appearing in the governing equations (eqn (5)–(9)). In Section 4, results from model 0 and model I are presented. Model 0 is not truly a Janus particle: the cap is absent (δ = 0, homogeneous mass distribution) and the surface potential of the two faces is identical (i.e., ψ₁ = ψ₂), making it impossible to distinguish between them. Model I differentiates the role played by inhomogeneity in the surface potentials (ψ₁ ≠ ψ₂) of the particle and is the simplest representation of a Janus particle considered in this work. Note that both model 0 and model I represent particle behavior in the absence of shear flow, i.e., [small gamma, Greek, dot above] = 0 and therefore V_x = 0. In Section 5 the cap thickness is introduced, i.e., we consider δ > 0, while still keeping [small gamma, Greek, dot above] = 0 (hence V_x = 0) and we refer to this formulation as model II. In model II as well as in the subsequent models III and IV, where δ > 0, the geometry of the cap is assumed to be thickest at the pole and becomes gradually thinner as it approaches the equator of the particle, which agrees well with recent experimental observations.⁵² The cap volume is estimated as the volume of a cylinder with a diameter of 2R of thickness δ, V_cap = πR²δ. Typically, Janus particles are coated with a thin layer of metal, e.g., platinum, which is denser than the SiO₂ core of the particle (e.g., ρ_Pt = 21.45 vs. ρ_SiO2 = 2.64 g cm⁻³). Consequently, both the density mismatch and the geometry of the cap render the particle “bottom-heavy”. Model III differs from model II as the effect of shear flow (i.e., [small gamma, Greek, dot above] ≠ 0) is included in the Langevin equation through a non-zero velocity term, V_x. The system Pe number at the highest shear flows studied here ([small gamma, Greek, dot above] = 8 s⁻¹) is small (Pe < 0.006) using the Debye length (κ⁻¹) as the geometric length scale. Note even when the particle diameter is used as the geometric length scale, the system Pe < 1. For the particle, the ratio of the advection to Brownian forces for the particle is expressed by the particle Péclet number, Pe_p = [small gamma, Greek, dot above]R²/D₀, where D₀ is the translational diffusivity of the isolated Brownian particle in the bulk (ESI,† Section S1). Depending on the particle size and strain rate, Pe_p numbers vary over the range of 0 ≤ Pe_p ≤ 0.5 × 10⁴. Subsequently, the effect of the hydrodynamic friction is evaluated by inclusion of frictional parameters, f_x and f_θ, in the translational and rotational terms in model IV. Table 1 summarizes the model parameters for radius (R), surface/wall potentials (ψ₁, ψ₂, and ψ_w), cap thickness (δ), and shear rate ([small gamma, Greek, dot above]).

Table 1 Parameters used in this study

Parameter	Symbol	Value
Particle radius	R	0.75, 1, 2, 4 μm
Potential on face 1	ψ 1	−20, −40 mV
Potential on face 2	ψ 2	−20, −40 mV
Wall potential	ψ w	−20, −50 mV
Cap thickness	δ	0–100 nm
Shear rate	[small gamma, Greek, dot above]	0–8 s^−1

---

Note that, in all equations the particle height (z) and position (x), time (t), force (F), and the torque (N) are non-dimensionalized by the Debye length (κ⁻¹), the characteristic time , thermal force (k_bT/κ⁻¹), and thermal energy (k_bT), respectively. Non-dimensionalized variables are indicated by a bar, e.g., z̄.

3 Methods

A BD simulation is conducted to evaluate the rotational and translational trajectories of a Janus particle by integrating the Langevin equations forward in time employing the explicit Euler time integration scheme. Thermal fluctuations are incorporated in the integration by sampling from the normal distribution N(0,1). Note the standard white noise assumption for the Brownian force applies for the system studied here as the ratio of [small gamma, Greek, dot above] over the bath frequency (≈10¹² s⁻¹) is less than 10⁻¹¹. For larger ratios, i.e., very large [small gamma, Greek, dot above] or more viscous fluids, a non-white noise more accurately describes the random fluctuations.⁵³ Each state point is integrated from time t₀ up to time t₁ and N_p independent particles are simulated. Specific settings of numerical simulations are detailed in Section S2 of the ESI.†

Analysis of the trajectories results in PDFs that show the frequency with which the particle is found in a particular orientation and height at various conditions. In order to assess and minimize bias in the behavior of the Janus particle with respect to the initial position and orientation, random initial positions and orientations are utilized. Equilibrium behavior of the Janus particle is verified to be independent of the initial condition and the early-time stages of the dynamics are disregarded from the data analyzed, removing data from the beginning of the simulation up to the equilibration time t_eq when necessary. See the ESI,† Section S2 for details on the initial condition and on the equilibration time. This approach ensures that the results obtained studying different trajectories of the same Janus particle are statistically equivalent.

When the shear flow is included in the model, i.e., model III, the focus is particularly on the rotational degree of freedom of the Janus particle. To characterize this, the average angular velocity of each trajectory is measured by temporally averaging the velocity as

where the instantaneous angular velocity v_θ(t̄) is simply estimated as (θ(t̄ + dt̄) − θ(t̄))/dt̄. The ensemble average of the angular velocity is then computed by averaging over different simulations, obtaining a constant angular velocity value, ω. The ω value is used as the order parameter of a dynamical transition between a rotating state, where ω > ω₀, and a non-rotating state, where ω ≤ ω₀. The two states are distinguished by setting the threshold ω₀ = 0.05 and verifying the robustness of the results by considering similar values of the threshold.

The trajectories, θ(t̄), exhibit a rich behavior, which is classified using statistical testing. Specifically, time is binned and an average value of the instantaneous velocity v_θ(t̄) is computed for each bin to obtain a time-dependent average angular velocity ω(t̄) (see the ESI,† Section S2), whose behavior is fitted to two functions: a constant and a sinusoidal. A constant behavior of ω(t̄) indicates that the particle rotates with a velocity that is independent of its orientation and possibly does not rotate at all, if ω = 0; a sinusoidal behavior indicates that the rotation is accelerated over time, an effect that is the consequence of the complex interplay between shear forces, gravitational attraction between the particle and the wall and electrostatic repulsion between the wall and the two sides of the Janus particle (which in general have different surface potentials and hence experience different repulsion strengths). The best model among the constant and the sinusoidal is selected using the Akaike Information Criterion,⁵⁴ which appropriately accounts for the larger number of fitting parameters (four) that the sinusoidal model has over the constant (which has one). In this way, each trajectory is classified as being sinusoidal or constant, providing a number of classifications equal to the number of trajectories simulated per state point. The more populated class is selected as the most representative of the state point. Noting that the classification results in two possible classes for each state point, and its reliability is further tested as follows. The class is regarded as a binary variable and the null hypothesis that the outcome of the classification results from the sampling of a binomial distribution with parameter 1/2 and with N_p trials is formulated. The null hypothesis is tested setting the p value to 0.01 and computing the probability of the classification outcome under the null hypothesis. The classification is regarded as reliable if this probability is found to be smaller or equal to p. In this way, each state point is classified as sinusoidal, constant or is not classified at all (i.e., when the classification is regarded as not reliable).

4 Dynamics of a Janus particle near a wall

The behavior of a homogeneously charged (i.e., constant surface potentials: ψ₁ = ψ₂) “Janus” particle of R = 1 μm without a cap (δ = 0) near a wall (SiO₂) with negative surface potential of ψ_w = −50 mV is studied first to test the formulation and validity of the BD simulation, i.e., model 0, and compared to a Janus particle without a cap (δ = 0), but with differing surface potentials (ψ₁ ≠ ψ₂), i.e., model I, in Fig. 2.

Fig. 2 Impact of surface potential, ψ, on the probability density function (PDF) for a Janus particle of R = 1 μm without a cap (δ = 0), with surface potentials, ψ₁ and ψ₂, interacting with a wall of surface potential, ψ_w = −50 mV, using model 0: (A) PDF and (B) equilibrium height for ψ₁ = ψ₂ = −20 mV and model I: (C) PDF for ψ₁ = −20 mV and ψ₂ = −40 mV. (D) Boltzmann distribution prediction of the orientation distribution for the particle simulated in (C).

Fig. 2(A) shows the expected PDF with equal probability for all orientations for a “Janus” particle with uniform surface potential (ψ₁ = ψ₂ = −20 mV, model 0) near a wall with ψ_w = −50 mV. The PDF is uniform because the driving force to rotate each face of the particle toward the wall is equal for both particle halves due to the identical surface potential. The particle rotates around its center, O, at an equilibrium height of z̄ = 6.19 ± 0.45 as shown in the height trajectory of the particle in Fig. 2(B).

The PDF shown in Fig. 2(C) is obtained for a Janus particle where the surface potential of the two halves differ (model I) and are specified as ψ₁ = −20 mV and ψ₂ = −40 mV, while the wall surface potential is kept at ψ_w = −50 mV. Comparing the PDFs in Fig. 2(A) and (C), it is apparent that the particle stays at a similar equilibrium height but shows a bias in its configuration, i.e., the particle tends to orient the face of lower surface potential, face 1, toward the wall (θ = 0) due to the greater surface potential gradient. Calculation of the equilibrium height reveals z̄ = 6.32 ± 0.50 slightly increased from the particle with even surface potential. Careful inspection of the height distribution reveals that the particle is closer to the surface when face 2 is rotated away from the wall than when rotated toward the wall (z̄_max,0 = 6.07 vs. z̄_max,π = 6.70) explaining the increased average equilibrium height value.

The dynamics of a Janus particle with the total potential interaction energy of U_tot(z,θ) should follow the Boltzmann distribution prediction for which the probability distribution of various configurations is given as with A_p chosen such that . Fig. 2(D) shows the Boltzmann distribution prediction for the same parameters used in the simulation of the PDF shown in Fig. 2(C). Comparison of Fig. 2(C) and (D) shows agreement thereby inferring the accuracy of the BD simulation methodology used here. Note data shown in the ESI† indicate that this holds true also for the plain particle model (see the ESI,† Section S3).

5 Dynamics of a bottom-heavy Janus particle near a wall

After evaluating the accuracy of models 0 and I for a uniformly charged “Janus” particle with equal surface potentials (ψ₁ = ψ₂) and the simplest Janus particle with differing surface potentials (ψ₁ ≠ ψ₂), respectively, model II introduces the effect of the cap, here made of platinum (ρ_Pt = 21.45 g cm⁻³), creating a bottom-heavy Janus particle. The weight of the cap and its geometry are incorporated in the simulations as described in Section 2. In this context, the cap of the particle (i.e., its bottom-heaviness) contributes to the rotational dynamics of the particle and tends to rotate the cap toward the wall to antiparallel align the particle director (n_p) with the wall's director (n_w) via the corresponding torque of N_g.

Fig. 3 shows the effect of incorporating a Pt cap (δ > 0) in the dynamical behavior of a particle with ψ₁ = ψ₂ and a Janus particle (ψ₁ ≠ ψ₂) with varying δ. Generally speaking, the bottom-heavy Janus particle stays at an equilibrium height and fluctuates with different preferred configurations near the wall due to an interplay of the effects of the cap and surface potentials. The PDF in Fig. 3(A) shows that for the case in which ψ₁ = ψ₂ (the wall potential is kept at ψ_w = −50 mV), the particle prefers rotating the face with the cap, face 2, toward the wall (θ = π). Furthermore, according to the height trajectory (Fig. 3(B)), the bottom-heavy Janus particle fluctuates at an equilibrium height of z̄ = 6.09 ± 0.46 slightly smaller (and therefore closer to the wall) compared to the “Janus” particle with uniform surface potential in Fig. 2(B) due to the additional gravitational force of the cap.

Fig. 3 Impact of cap thickness and surface potential on the probability density function (PDF) of a bottom-heavy Pt-capped Janus particle of R = 1 μm with cap thickness, δ, and surface potentials, ψ₁ and ψ₂, interacting with a wall of surface potential, ψ_w = −50 mV using model II: (A) PDF and (B) equilibrium height z̄ for ψ₁ = ψ₂ = −20 mV and δ = 10 nm. (C) Same as (A) with ψ₂ = −40 mV. (D) Same as (C) with δ = 30 nm.

Fig. 3(C) displays the PDF for the case in which the cap weight and the more negative surface potential on the same particle face, face 2, exert competing demands on the particle orientation; the cap forces face 2 towards the wall (θ = π) due to the gravitational torque, while the more negative surface potential forces face 2 away from the wall (θ = 0). For δ = 10 nm, the surface potential gradient dominates over the cap torque, and the particle exhibits a preference for rotation of face 2 away from the wall. The average height is found to be z̄ = 6.37 ± 0.51, with the preferred height for the cap-up configuration (θ = 0) peaking at z̄ = 5.99, while the particle in cap-down configuration (θ = π) is found most often at z̄ = 6.70, in good agreement with the stronger repulsion between the wall and face 2. An increase in cap thickness to δ = 30 nm leads to the dominance of the gravitational cap torque that quenches the rotational motion of the particle to a mostly cap-down configuration (θ = π) as shown in Fig. 3(D). The average equilibrium height of the particle during the simulation is z̄ = 6.50 ± 0.45 with cap-up and cap-down configuration heights peaking at z̄ = 5.80 and 6.51, respectively. It is clear that the dynamic rotational behavior and the height of a bottom-heavy Janus particle depend strongly on the interplay of cap and surface properties. The complexity of this interaction is increased further by the introduction of shear flow and particle radius variations discussed in the next section revealing an intricate Janus particle flow behavior near a wall.

6 Dynamics of a bottom-heavy Janus particle near a wall under shear flow

In model III, the bottom-heavy Janus particle is subjected to shear flow with strain rate [small gamma, Greek, dot above] acting on the particle center, O. Initial simulations are performed with a strain rate of [small gamma, Greek, dot above] = 6 s⁻¹ for sets of three Janus particles with R = 1 μm, δ = 86 nm, ψ₁ = −20 mV, ψ₂ = −40 mV, and ψ_w = −50 mV. Data for translation in x̄, θ, and z̄ as a function of t̄ are provided for the three particles in the ESI,† Section S4 (Fig. S7, top left). The three traces shown for x̄(t̄) indicate that the Janus particles move forward in the direction of the flow. The θ(t̄) traces show that each Janus particle starts in a unique, random orientation. Over time, θ values fluctuate between 0 and 2π in a cyclic fashion implying that the Janus particle rotates from a cap-down orientation (θ = 0) to a cap-up orientation (θ = π) back to a cap-down orientation (θ = 2π). Note that the particle orientation θ = 0 is equivalent to θ = 2π resulting in the abrupt change in θ at the end of each cycle. The z̄(t̄) traces closely follow the θ(t̄) traces with the particles at a high z̄ when the particles are cap-down and a lower z̄ when the particles are cap-up. A change in particle radius to R = 4 μm at the same [small gamma, Greek, dot above] = 6 s⁻¹ and cap thickness, and for the same simulation length shows the expected longer distance x̄(t̄) and overall smaller z̄ values with the oscillations in z̄ more clearly matching the particle rotation (ESI,† Fig. S7, bottom left). Decrease of [small gamma, Greek, dot above] to 1 s⁻¹ (ESI,† Fig. S7, top and bottom right) leads to a different θ behavior for both Janus particle sizes, where the Janus particle rotates until it assumes a constant θ value, i.e., the strain rate is not strong enough to induce rotation (ω < ω₀) and the particle slides along the wall with constant θ. The distinct behavior observed in the cap rotation, i.e., θ(t), between the two particle sizes and as a function of [small gamma, Greek, dot above] prompts a careful search of the [small gamma, Greek, dot above],δ-parameter space for R = 0.75, 1, 2, and 4 μm with ψ₁/ψ₂/ψ_w combinations of −20/−20/−20 mV, −20/−20/−50 mV, −40/−20/−20 mV, −40/−20/−20 mV, −20/−40/−20 mV, and −20/−40/−50 mV (see the ESI,† Section S5 and Fig. S8–S12). The findings from this parameter scan are summarized here.

Fig. 4 displays the angular velocity ω (see Section 3) of a Janus particle with the same ψ₁, ψ₂ and ψ_w values considered in Section 5 (Fig. 3) for the same particle radius R = 1 μm (Fig. 4(A) and (B)) and for a larger particle with radius R = 4 μm and surface potentials ψ₁ = −20 mV, ψ₂ = −40 mV and ψ_w = −50 mV (Fig. 4(C)) for cap thicknesses δ = 0–100 nm and strain rates [small gamma, Greek, dot above] = 0–8 s⁻¹. All systems display the transition between a rotating (ω > ω₀) and a non-rotating behavior (ω<ω₀) described above as marked by the white line in Fig. 4. Above the transition line, ω increases with increasing [small gamma, Greek, dot above] at a constant cap thickness δ, whereas at constant [small gamma, Greek, dot above] and increasing δ, the change in ω is non-trivial. As the distance between the particle center, O, and the top of the cap, i.e., the lever arm, grows linearly with R at a given δ, the torque acts more effectively on large particles, leading to the higher values of ω observed at the larger R (see colors in Fig. 4(A), (B)vs.Fig. 4(C)).

Fig. 4 Angular velocity maps for a bottom-heavy Janus particle under varying shear flow ([small gamma, Greek, dot above] = 0–8 s⁻¹) and cap thickness (δ = 0–100 nm) using model III for: (A) R = 1 μm, ψ₁ = −20 mV, ψ₂ = −20 mV and ψ_w = −50 mV. (B) As in (A), but ψ₂ = −40 mV. (C) As in (B), but R = 4 μm. The color represents the average angular velocity as displayed by the color bar on the right. The transition line (white) is drawn at ω = 0.05 s⁻¹.

The onset of a non-negligible angular velocity at δ = 0 follows the expected trends by increasing from Fig. 4(A)–(B) due to the higher surface potential asymmetry and decreasing from Fig. 4(B) and (C) due to the increase in particle radius. However, once the cap is introduced (δ > 0) and increased in thickness, the transition line shows a distinct behavior for each system, ranging from a monotonic increase (Fig. 4(A)) to a slight decrease followed by a slow increase (Fig. 4(B)) and sharp decrease followed by a linear increase. The transition line, a locus of points in the diagram of the form (δ,[small gamma, Greek, dot above]), fully characterizes the transition: above (below) the locus the system is in the rotating (non-rotating) state. For example, if two points of the locus (δ₁,[small gamma, Greek, dot above]₁) and (δ₂,[small gamma, Greek, dot above]₂) are chosen such that the system is at the transition at those points, it is true for (almost) all systems that if δ₁ < δ₂, then [small gamma, Greek, dot above]₁ < [small gamma, Greek, dot above]₂. In other words, it is (almost) always true that increasing the cap thickness requires a stronger shear flow for the rotating phase to be entered (see Fig. 4(A) and the ESI,† Section S5). The only exception to this behavior is observed for Janus particles whose heavy side (face 2) is more repelled (ψ₂ = −40 mV) by the wall then the light side (face 1, ψ₁ = −20 mV) as shown in Fig. 4(B) and (C) (see also the ESI,† Section S5, Fig. S8 and S9). In this case, it is observed that δ₁ < δ₂ results in [small gamma, Greek, dot above]₁ ≥ [small gamma, Greek, dot above]₂ for sufficiently small δ (lower left corner). This behavior is persistent across different particle sizes and is more discernible for particles with R ≥ 2. A different way of describing this phenomenon is to see the locus of points as a set of values of [small gamma, Greek, dot above] that is a function of δ. In this perspective, it can be stated that [small gamma, Greek, dot above] is (almost) always a growing function of δ, with an exception at small δ when the surface potential asymmetry (ψ₁ less negative than ψ₂) makes the heavy face, face 2, more repelled by the wall. A minimum in the transition line then represents the cap thickness at which the opposing cap torque and surface potential gradient balance each other.

Surprisingly, the location of the non-rotating region depends only weakly on the particle radius, surface potential balance, and cap thickness. This behavior is explained by taking into account the complex interplay between particle radius, surface potential, and cap thickness, which together determine – for a given strain rate – the average distance z̄ between the particle and the wall (see the ESI,† Section S5, and Fig. S10). For example, the strain rate is less effective when a more negative ψ₂ than ψ₁ results in a low z̄, but is more effective when a large R and/or a strong asymmetry in the particle's mass distribution due to a large δ also result in a low z̄.

At fixed R (Fig. 4(A)vs.Fig. 4(B)), an increase in electrostatic repulsion only slightly shifts the transition line to lower [small gamma, Greek, dot above] values, as it pushes the bottom-heavy particles farther from the wall, thus favoring their rotation. On increasing R (Fig. 4(B)vs.Fig. 4(C)), the lever arm grows, favoring the rotating behavior at a given [small gamma, Greek, dot above], but this kind of increase is compensated by an increase in particle weight, which brings the particle closer to the wall, thus disfavoring the particle rotation due to the electrostatic repulsion. The previous argument helps to clarify the aforementioned complex dependence of ω on δ at constant [small gamma, Greek, dot above]. In particular, a non-monotonic behavior is observed when the repulsion strength between the wall and the heavy side of the particle is strong. This is because, as a consequence of the electrostatic repulsion, the particle is further away from the wall when the heavy side is closer to it. This effect interplays with the previously noted ones and triggers the rotation of particles with thin caps (including the δ = 0 case) at values of [small gamma, Greek, dot above] that would not be sufficient to trigger it in the absence of the surface potential asymmetry: the surface potential asymmetry pushes the particle upward enough to allow [small gamma, Greek, dot above] to act more effectively on it, hence triggering the rotational behavior.

The phenomenology is, however, richer than what can be deduced from Fig. 4 only. A visual inspection of individual trajectories (Fig. 5) for the same system under different conditions shows evident differences. All curves are for R = 4 μm, ψ₁ = −20 mV, ψ₂ = −40 mV and ψ_w = −50 mV, i.e., they correspond to state points in the angular velocity map shown in Fig. 4(C). In Fig. 5(A), the curve for δ = 86 nm and [small gamma, Greek, dot above] = 1.0 s⁻¹ (blue) is an example of a trajectory sampled in the non-rotating phase and indeed the particle orientation simply experiences Brownian fluctuations. On increasing the shear rate to 6.0 s⁻¹ while reducing the cap thickness to 16 nm (orange), the particle enters the rotating phase. The trajectory shows a growth of θ with time that is well approximated by a linear function, indicating a constant angular velocity. On increasing the cap thickness back to 86 nm (green), the behavior of θ(t̄) changes significantly: the first half of the rotation, for 0 ≤ θ < π is much faster than the second half. This change in angular velocity happens because the cap is heavy so gravity tends to keep it close to the wall and more work is required for the torque to fully rotate a particle in the opposite orientation, where the heavy side is aligned with the wall normal n_w (θ = 0). This behavior is ascribed to an oscillating trend of the angular velocity, which is higher in the first part of the rotation and lower in the second. The observation motivates the classification introduced in Section 3. Fig. 5(B) illustrates a difficulty encountered by the classification. The curves in Fig. 5(B) show a mixed behavior despite corresponding to different simulations of the same state point in the angular velocity map. One of the curves (red) shows a single rotation over the whole temporal window of the numerical integration (t̄ ≈ 47), the purple curve also shows a single rotation, but almost at the end of the integration (t̄ ≈ 52) and the cyan curve does not show a rotation at all, despite the fact that a rotation would be observed if a longer numerical integration time were used. However, there is not a length of the numerical integration time that guarantees observation of at least one rotation for systems that exhibit rotations and not to observe rotation at all for systems that do not undergo rotations. When the model parameters are such that the rotational behavior exists, but the angular velocity is small, two different scenarios can be observed: very slow rotations with constant angular velocities, or the limit behavior of the sinusoidal angular velocity previously described.

Fig. 5 Representative angular trajectories for a system with R = 4 μm, ψ₁ = −20 mV, ψ₂ = −40 mV, and ψ_w = −50 mV. (A) Angular trajectories in a short-time window for state points specified in the legend shown in (B). (B) Angular trajectories in a long-time window for systems whose parameters are given in the legend. Note that the three trajectories shown in (B) correspond to different simulations of the same system. The dashed vertical line is the equilibration time after which the data are analyzed.

Fig. 5(B) is an example of this kind of limit behavior: as shown by the red curve, the time it takes to complete the first half of the rotation is extremely short, due to the heavy cap driving face 2 toward the wall, while the time it takes to complete the second half of the rotation is exceedingly large as shear flow and surface potential repulsion must overcome the gravitational torque from the cap. This phenomenon can be described by the first passage time distribution as a parametric function of the orientation θ:⁵⁵ this kind of function describes the probability that θ reaches a certain (parametric) value in a time t*, which is the random variable of the distribution. Stochastic processes near a phase transition are typically characterized by a fat tailed distribution of the first passage time⁵⁶ and the results presented here suggest that the system under analysis makes no exception: for a large value of the parameter θ, the distribution would likely appear heavy tailed, indicating that no matter how long the numerical integration is, the probability to not observe a full rotation would still be non-zero.

Note the regime of constant angular velocity can also be seen as a purely rotational regime, where the translational motion of the particle is strongly correlated with its orientation, as one expects for a ball rotating on a plane, while the variable angular velocity, induced by surface potential imbalance and asymmetry of the mass distribution, can be seen as a type of motion in which the correlation between orientation and position along the plane is increasingly weaker as the linearity in θ(t) is lost, due to a “slipping” behavior in which the particles translate while rotating independently of the translating motion.

The results of the classification are shown as state diagrams in Fig. 6 and in the ESI,† Section S6. The transition line between the rotating and the non-rotating behavior is identified by the white transition line. Almost all the state points below it are ascribed to constant angular velocity (green state point) with a small region of unclassified states (red state point) inside the non-rotating phase only for large particle radius, small cap thickness, and small shear flow. Just above the transition line, a region of unclassified states points is systematically observed. This is due to the phenomenon exemplified in Fig. 5(B) and explained in the previous paragraph: near the transition, some realizations appear as constant (cyan curve) and some are clearly not well fitted by a constant function (red and purple curves), resulting in state points that are not classified at all. Further above the transition line, at large cap thickness and strong shear flow, the system is systematically characterized by a sinusoidal angular velocity (blue state point). On progressively reducing the cap thickness, the system enters first a region of unclassified state points in which the distinction between the two behaviors is not simple as the difference in slopes between the two regimes of rotation is extremely small, but not small enough for the state points to be systematically ascribed to a constant behavior. On further reducing the cap thickness, a region of constant angular velocity is entered. This region displays an interesting behavior as the model parameters change: on reducing R, it systematically shifts to the left (Fig. 6(C)vs.Fig. 6(B) and ESI,† Fig. S12), it is observed to shift to the right (Fig. 6(A)vs.Fig. 6(B) and ESI,† Fig. S11, S12) if the surface potential of the heavy side of the particle is increased (in absolute value, at fixed R) and it shifts upward if the other side is overcharged with respect to the wall (ESI,† Fig. S11). See the ESI,† Section S6 for a more detailed analysis of the results.

Fig. 6 State diagram of a bottom-heavy Janus particle under varying shear flow ([small gamma, Greek, dot above] = 0–8 s⁻¹) and cap thickness (δ = 0–100 nm) for: (A) R = 1 μm, ψ₁ = −20 mV, ψ₂ = −20 mV and ψ_w = −50 mV. (B) Same as in (A), but ψ₂ = −40 mV. (C) Same as in (B), but R = 4 μm The white line is the same transition line drawn in Fig. 4, and colors represent the classification of each state point as described by the legend on the right. Symbols in the legend stand for sinusoidal (S), i.e., non-constant angular velocity, constant (C), i.e., constant angular velocity, and unclassified (U).

Despite being apparently intricate, this behavior is relatively simple to understand. Constant but sinusoidal angular velocity is observed when the shear rate is strong enough to induce rotational motion (above the transition line), but the surface potential and density imbalances are not so strong that one side is more prone to be closer to the wall than the other. Hence, at comparable repulsion between the two sides and the wall (ψ₁ = ψ₂), only the cap induces asymmetry and this can be compensated when the cap is small and the shear flow is large. Reducing ψ₂ to more negative potentials means that the light side of the particle is more likely to stay close to the wall, which is balanced by the cap pushing the particle in the opposite orientation, hence increasing the width of the constant regime. When the light side is strongly repelled by the wall, the region shifts upward as a given asymmetry in the mass distribution requires a strong shear rate to compensate the repulsion at small cap sizes. Furthermore, when the constant behavior region is observed at large cap thickness, it is also seen that a new sinusoidal region exists on its left, when both cap thickness and shear rate are small (ESI,† Fig. S11 and S12, bottom row).

Summarizing the behavior of bottom-heavy Janus particles under shear flow, it can be stated that outside of the rotating region the system is characterized by a sinusoidal angular velocity if the cap is thick enough to induce a strong asymmetry. Reducing the cap thickness at a given shear rate implies reducing the variability in ω(t̄) up to completely balancing the speed of rotation in the two half periods and, possibly, on further reducing the cap thickness a sinusoidal angular velocity is observed again (especially for large particles), as the surface potential imbalance is not compensated by the weak flow. Data provided in the ESI† (Sections S5 and S6) shows that this behavior is persistent across different model settings.

Last but not least, in model IV, the impact of hydrodynamic friction is considered taking into account the Janus particle's closeness to the wall for the R = 1 μm and R = 4 μm systems discussed at the beginning of this section (δ = 86 nm, ψ₁ = −20 mV, ψ₂ = −40 mV, and ψ_w = −50 mV). Based on Goldman et al.'s report,⁴⁹ frictional factors f_x and f_θ, which are functions of z̄, are determined and introduced into eqn (5) and (7), respectively. For a bottom-heavy Janus particle of R = 1 μm and cap thickness of 86 nm at an average height of z̄ = 6 ± 1, f_x = 0.9887 and f_θ = 0.9882 indicating that these particles experience very little hydrodynamic friction. The situation changes when the particle size is increased to R = 4 μm and the particles display an average z̄ = 3 ± 0.5 resulting in an f_x = 0.8342 and f_θ = 0.8428. Inclusion of these friction factors in the simulation reduces the angular velocity of the R = 4 μm Janus particle resulting in approximately 20% fewer rotations as shown in the ESI,† Section S7, Fig. S13, while the R = 1 μm Janus particle is too far away from the surface to be substantially impacted by friction. Note that Goldman et al.'s⁴⁹ friction coefficients were calculated for a uniform, neutrally buoyant sphere near a plane wall under simple shear flow. Geonzon et al.⁵⁷ experimentally verified Goldman et al.'s⁴⁹ friction coefficients for unmodified SiO₂ spheres (R = 1.25 μm) under shear flow using optical tweezers. More recently, the translational and rotational drags of a Janus particle close to a wall without shear flow were determined using fluorescently labeled Pt-coated MF beads of R = 1.24 ± 0.05 μm and were also found to be in good agreement with Goldman et al.'s values.⁵⁸ Experimental verification of the effect of shear flow on the hydrodynamic drag force for Janus particles near a wall would be desirable. Unfortunately, Janus particles with metal caps experience thermophoresis in optical traps due to the heating of the metal cap making such measurements non-trivial.⁵⁹

7 Conclusion

The behavior of a bottom-heavy, charged Janus particle near a charged wall under shear flow is fully rationalized, providing insights into how both confinement and shear flow can be utilized to manipulate the particle's translational and rotational motion.

The dynamical behavior of a Janus particle is described using a model based on gravitational and electric double layer interactions. The simplest Janus particle model, model I, shows that the Janus particle without a cap of differing density tends to orient toward the wall with the face of lower surface potential because of the greater potential gradient. For the case of a uniformly charged bottom-heavy Janus particle, i.e., a particle with a cap of mismatched density but uniform surface potential, model II, the cap quenches the rotational motion of the particle leading to a preference for a cap-down orientation adjacent to the wall in good agreement with existing literature.^25,26 Adjustment of the caps surface potential leads to a competition between gravitational torque and electrostatic repulsion in which the position of the particle above the wall and its rotational behavior become a function of radius, surface potential, and cap thickness.

The relationship is further complicated by the introduction of shear flow at small Pe_p in model III. Thorough analysis of the particle behavior for a range of shear flows (0–8 s⁻¹), cap thicknesses (0–100 nm), radii (0.75–4 μm), and surface potential combinations (ψ_w = −20, −50 mV and ψ_1,2 = −20, −40 mV) reveals complex angular velocity profiles for the bottom-heavy Janus particle with regions of no rotation, constant rotation, and time-varying rotation. Use of statistical analysis leads to state diagrams that reveal δ,[small gamma, Greek, dot above]-combinations for which the passive, bottom-heavy Janus particle shows constant angular velocity with and without rotation and a sinusoidal angular velocity. The latter state can be described as the particle reaching a cap-down orientation followed by a sliding event before completing its rotation where the sliding time can be engineered and therefore is of potential interest when patch–wall interactions need to be timed in applications such as drug delivery or environmental remediation. Finally, we note that hydrodynamic friction, model IV, effectively reduces the particle's ability to rotate by ≈20% when the particle is close to the wall, which occurs for large radii and/or thick caps at high shear rates, but can be neglected for small particles and/or thin caps.

For future work, it is worth noting that the choice of orientation-dependent factors balances simplicity and numerical efficiency in the BD model employed here. A more rigorous formulation could be developed based on the Janus sphere–sphere potentials presented by Popov and Hernandez.⁶⁰ A detailed discussion on the selection of weighting factors and their influence on the torque is provided in Section S8 of the ESI.†

Author contributions

Z. J., U. C., and I. K. discussed and conceptualized the project. Z. J. formulated the theoretical model, developed the initial code and provided a first draft of the manuscript as part of her PhD thesis.⁶¹ I. K. revised the theoretical model and finalized the code. Numerical simulations were performed by Z. J. I. K. and D. N. D. N. conceptualized, implemented, and performed the statistical analysis needed for the angular velocity maps and state diagrams. D. N., E. B., and I. K. discussed the interpretation of angular velocity maps and state diagrams and revised the manuscript during I. K.'s sabbatical stay in E. B.'s group. All authors have read and approved the final version.

Data availability

Data supporting the statements and conclusions made in the manuscript are included as part of the ESI.† Additional data is available upon reasonable request from the corresponding authors.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Science Foundation (NSF) under award CBET-1705565. I. K. and U. C. also acknowledge support from NSF grant number EES-2112550. I. K. thanks Eric McPherson and Justin Pinca for the helpful discussions. D. N. and E. B. acknowledge that this research was funded in part by the Austrian Science Fund (FWF) [Project no. Y-1163-N27]. For open access purposes, the authors have applied a CC BY public copyright license to any author accepted manuscript version arising from this submission.

References

1. H. Su, C.-A. Hurd Price, L. Jing, Q. Tian, J. Liu and K. Qian, Mater. Today Bio, 2019, 4, 100033.
2. J. Jeevanandam, K. X. Tan, J. Rodrigues and M. K. Danquah, Pharmaceutics, 2023, 15, 1614.
3. K. X. Tan, M. K. Danquah, J. Jeevanandam and A. Barhoum, Pharmaceutics, 2023, 15, 423.
4. A. Enayati Gerdroodbar, R. Damircheli, S. N. Eliseeva and M. Salami-Kalajahi, J. Electroanal. Chem., 2023, 948, 117831.
5. C. Chen, L. Zhang, N. Wang, D. Sun and Z. Yang, Macromol. Rapid Commun., 2023, 44, 2300280.
6. C. Marschelke, A. Fery and A. Synytska, Colloid Polym. Sci., 2020, 298, 841–865.
7. E. L. Correia, N. Brown and S. Razavi, Nanomaterials, 2021, 11, 374.
8. M. Cooley, A. Sarode, M. Hoore, D. A. Fedosov, S. Mitragotri and A. Sen Gupta, Nanoscale, 2018, 10, 15350–15364.
9. S. Das, A. Garg, A. I. Campbell, J. Howse, A. Sen, D. Velegol, R. Golestanian and S. J. Ebbens, Nat. Commun., 2015, 6, 1–10.
10. J. Simmchen, J. Katuri, W. E. Uspal, M. N. Popescu, M. Tasinkevych and S. Sanchez, Nat. Commun., 2016, 7, 1–9.
11. Z. Xiao, M. Wei and W. Wang, ACS Appl. Mater. Interfaces, 2018, 11, 6667–6684.
12. J. P. Jean-Pierre Minier, Particles in Wall-Bounded Turbulent Flows: Deposition, Re-Suspension and Agglomeration, Springer, Cham, 1st edn, 2017.
13. Z.-G. Feng, J. Gatewood and E. E. Michaelides, J. Fluids Eng., 2021, 143, 081106.
14. B. Derjaguin and L. Landau, Prog. Surf. Sci., 1993, 43, 30–59.
15. E. J. W. Verwey, J. T. G. Overbeek and K. van Nes, Theory of the stability of lyophobic colloids, Elsevier Pub. Co., New York, 1948.
16. S. Asakura and F. Oosawa, J. Chem. Phys., 1954, 22, 1255–1256.
17. J. N. Israelachvili, Q. Rev. Biophys., 1973, 6, 341–387.
18. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
19. S. Rajupet, Langmuir, 2021, 37, 13208–13217.
20. J. M. Thwala, J. Mathunjwa, S. K. Mkhonta, G. C. Bwembya and B. Mamba, Colloids Surf., A, 2023, 670, 131315.
21. A. Specht, D. Krämer, N. Helfricht and G. Papastavrou, Langmuir, 2025, 41(10), 6515–6527.
22. J. Israelachvili, Y. Min, M. Akbulut, A. Alig, G. Carver, W. Greene, K. Kristiansen, E. Meyer, N. Pesika, K. Rosenberg and H. Zeng, Rep. Prog. Phys., 2010, 73, 036601.
23. D. C. Prieve, Adv. Colloid Interface Sci., 1999, 82, 93–125.
24. J. Yan and C. L. Wirth, J. Appl. Phys., 2022, 131, 150903.
25. A. Rashidi and C. L. Wirth, J. Chem. Phys., 2017, 147, 224906.
26. A. Rashidi, S. Razavi and C. L. Wirth, Phys. Rev. E, 2020, 101, 042606.
27. S. Rajupet, A. Rashidi and C. L. Wirth, Phys. Rev. E, 2021, 103, 032610.
28. M. Mohammadi, E. D. Larson, J. Liu and R. G. Larson, J. Chem. Phys., 2015, 142, 024108.
29. Y. Kobayashi and N. Arai, Langmuir, 2017, 33, 736–743.
30. E. Bianchi, A. Z. Panagiotopoulos and A. Nikoubashman, Soft Matter, 2015, 11, 3767–3771.
31. A. Nikoubashman, E. Bianchi and A. Z. Panagiotopoulos, Soft Matter, 2015, 11, 3946.
32. A. Nikoubashman, Soft Matter, 2017, 13, 222–229.
33. R. A. DeLaCruz-Araujo, D. J. Beltran-Villegas, R. G. Larson and U. M. Córdova-Figueroa, Soft Matter, 2016, 12, 4071–4081.
34. R. A. DeLaCruz-Araujo, D. J. Beltran-Villegas, R. G. Larson and U. M. Córdova-Figueroa, Langmuir, 2018, 34, 1051–1060.
35. A. Ramachandran and A. S. Khair, J. Fluid Mech., 2009, 633, 233.
36. M. Trofa, G. DAvino and P. L. Maffettone, Phys. Fluids, 2019, 31, 083603.
37. B. R. Si, P. Patel and R. Mangal, Langmuir, 2020, 36, 11888–11898.
38. A. Choudhary, K. V. S. Chaithanya, S. Michelin and S. Pushpavanam, Eur. Phys. J. E: Soft Matter Biol. Phys., 2021, 44, 97.
39. P. Vinze and S. Michelin, Phys. Rev. Fluids, 2024, 9, 014202.
40. P. Bag, S. Nayak and P. K. Ghosh, Soft Matter, 2024, 20, 8267–8277.
41. V. Mancuso, M. N. Popescu and W. E. Uspal, Soft Matter, 2024, 20, 8742–8764.
42. B. Ren, A. Ruditskiy, J. H. K. Song and I. Kretzschmar, Langmuir, 2012, 28, 1149–1156.
43. L. Baraban, R. Streubel, D. Makarov, L. Han, D. Karnaushenko, O. G. Schmidt and G. Cuniberti, ACS Nano, 2013, 7, 1360–1367.
44. R. Huhnstock, M. Reginka, A. Tomita, M. Merkel, K. Dingel, D. Holzinger, B. Sick, M. Vogel and A. Ehresmann, Sci. Rep., 2021, 11, 21794.
45. T. Huang, V. Misko, A. Caspari, A. Synytska, B. Ibarlucea, F. Nori, J. Fassbender, G. Cuniberti, D. Makarov and L. Baraban, Commun. Mater., 2022, 3, 60.
46. Y. Wen, J. Liu, W. Wang, P.-Y. Lai and P. Tong, Soft Matter, 2024, 20, 9208–9218.
47. A. Khoshzaban, A. Magazzú, M. G. Donato, O. M. Maragò, M. B. Unlu, M. N. Cizmeciyan and P. Elahi, ACS Photonics, 2025, 12, 1936–1943.
48. R. Hogg, T. W. Healy and D. W. Fuerstenau, Trans. Faraday Soc., 1966, 62, 1638–1651.
49. A. Goldman, R. Cox and H. Brenner, Chem. Eng. Sci., 1967, 22, 653–660.
50. G. Trefalt, S. H. Behrens and M. Borkovec, Langmuir, 2016, 32, 380–400.
51. S. L. Carnie and D. Y. Chan, J. Colloid Interface Sci., 1993, 161, 260–264.
52. A. Rashidi, M. W. Issa, I. T. Martin, A. Avishai, S. Razavi and C. L. Wirth, ACS Appl. Mater. Interfaces, 2018, 10, 30925–30929.
53. S. Pelargonio and A. Zaccone, Phys. Rev. E, 2023, 107, 064102.
54. K. P. Burnham and D. R. Anderson, Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach, Springer, New York, 2nd edn, 2002.
55. C. W. Gardiner, Handbook of Stochastic Methods: For Physics, Chemistry and the Natural Sciences, Springer, 3rd edn, 2009.
56. S. di Santo, P. Villegas, R. Burioni and M. A. Muñoz, Phys. Rev. E, 2017, 95, 032115.
57. L. C. Geonzon, M. Kobayashi and Y. Adachi, Soft Matter, 2021, 17, 7914–7920.
58. V. Sharma, F. Fessler, F. Thalmann, C. M. Marques and A. Stocco, J. Colloid Interface Sci., 2023, 652, 2159–2166.
59. S. Auschra, A. Bregulla, K. Kroy and F. Cichos, Eur. Phys. J. E: Soft Matter Biol. Phys., 2021, 44, 90.
60. A. Popov and R. Hernandez, J. Phys. Chem. B, 2023, 127, 1664–1673.
61. Z. Jalilvand, Motion of Janus and Patchy Particles Near Various Boundaries, PhD thesis, CUNY Academic Works, 2020.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5sm00229j

---