This article describes the simulated Brownian motion of a sphere comprising hemispheres of unequal zeta potential (i.e., “Janus” particle) very near a wall. The simulation tool was developed and used to assist in the methodology development for applying Total Internal Reflection Microscopy (TIRM) to anisotropic particles. Simulations of the trajectory of a Janus sphere with cap density matching that of the base particle very near a boundary were used to construct 3D potential energy landscapes that were subsequently used to infer particle and solution properties, as would be done in a TIRM measurement. Results showed that the potential energy landscape of a Janus sphere has a transition region at the location of the boundary between the two Janus halves, which depended on the relative zeta potential magnitude. The potential energy landscape was fit to accurately obtain the zeta potential of each hemisphere, particle size, minimum potential energy position and electrolyte concentration, or Debye length. We also determined the appropriate orientation bin size and regimes over which the potential energy landscape should be fit to obtain system properties. Our simulations showed that an experiment may require more than 106 observations to obtain a suitable potential energy landscape as a consequence of the multivariable nature of observations for an anisotropic particle. These results illustrate important considerations for conducting TIRM for anisotropic particles.

Nanometer to micrometer scale colloidal particles are ubiquitous to consumer products, coatings, cosmetics, and large and small scale industrial processes. Colloidal particles dispersed in a liquid interact via surface forces to form microstructures that will have a profound impact on the mechanics and rheology of that liquid. Surface forces arise from physiochemical properties of the particle and liquid medium, such as surface charge, solution pH, and salinity, and the presence of other entities, such as a surfactant or polymer dispersed in the liquid.1 Multiple experimental techniques, including Atomic Force Microscopy (AFM), the Surface Force Apparatus (SFA), and Total Internal Reflection Microscopy (TIRM), have been developed over the past few decades to probe these interactions, which could be as weak as a few ∼kT. SFA and AFM, modified with a colloidal probe, are robust techniques that have been widely adopted to measure the surface forces between the probe and a substrate.2,3 Similarly, TIRM has been developed to measure the surface force experienced by a single colloidal particle immersed in a liquid nearby, but not adhered to, a substrate. The essence of TIRM is the measurement of stochastic “Brownian” fluctuations of the colloidal particle normal to the nearby substrate.4–6 TIRM tends to be more sensitive than AFM and SFA because TIRM uses a “thermal” energy scale (i.e., ∼kT), whereas AFM and SFA operate on a mechanical energy scale (i.e., ∼kSX2/2, kS is the cantilever spring constant).

TIRM, employing spherical particles with uniform surface chemistry and composition, has been used to accurately measure both equilibrium and non-equilibrium ∼kT scale colloidal interactions, including electrostatic double layer repulsion,7 van der Waals attraction,8 steric interactions,9 depletion interactions,10 electrophoretic forces,11 hydrodynamic interactions,12 and even the Casimir interaction.13 Recently, the TIRM technique has even been developed to measure interactions of colloidal particles and oil droplets near a water-oil interface, thereby extending the technique to liquid-liquid interfaces from conventional liquid-solid interfaces.14 These studies illustrate the utility of TIRM as a tool, which complements SFA and AFM, for surface force measurement. One current limitation of the TIRM technique is that only optically isotropic particles (i.e., sphere of uniform composition) have been used for measurements, with the exception of one example where the particle was small enough to assume that anisotropy did not impact scattering.15 Only recently, there have been efforts to develop the experimental technology and methodology for the application of TIRM to geometrically or chemically anisotropic colloidal particles.16 

Surface forces experienced by anisotropic particles are relevant to many natural and synthetic systems, for instance, carbon nanotubes, graphene sheets, red blood cells, and clay particles, which are all anisotropic. Anisotropic particles can be locally hard or soft, may have complicated shapes, or have a non-uniform surface charge distribution.17 Many groups are now able to synthesize model anisotropic particles, such as spherical and ellipsoidal “Janus” particles.18–20 A significant amount of work has focused on the directed self-assembly of Janus spheres or ellipsoids with different shapes and dissimilar properties.21–24 Surface forces, and associated particle dynamics, are of central importance to the directed self-assembly and ultimate utility of these materials. Non-uniform or non-symmetric physical and chemical properties typically induce spatially and orientationally dependent colloidal interactions, thereby complicating the measurement, prediction, and control of such interactions among anisotropic particles. Researchers have proposed techniques for tracking the rotational motion of a cluster of colloidal particles,25 copper oxide nanorods,26 and micrometric wires by using holographic video microscopy and 2D video microscopy. Recent work has even demonstrated how to account for particle curvature when calculating surface forces for a non-spherical particle.27 However, there is still a need to develop a technique and methodology for sensitively measuring the surface forces between anisotropic particles and a nearby substrate.

We are working to fill this need by developing the experimental technology, methodology, and design for conducting TIRM on anisotropic particles, including spheroids of arbitrary shape and spheres of non-uniform surface chemistry. In this article, we present the results from Brownian dynamics simulations of a Janus sphere very near a boundary at conditions similar to those of a TIRM experiment. We chose to initially work with Janus spheres because they are the closest departure from an isotropic sphere. As was done previously for isotropic particles,28 the goal of these simulations was to generate in silico data similar to those obtained during a TIRM experiment, allowing us to probe and optimize TIRM parameters that would be implemented for a Janus particle with hemispheres of unequal zeta potential, develop a methodology for calculating potential energy landscapes from multivariable trajectory data, and determine how the interpretation of conservative forces may be impacted by non-equilibrium dynamics. We combined a well-established Brownian dynamics stepping algorithm for the uncoupled rotational-translational Langevin equation with a meshing technique for application to a Janus sphere. We then applied this meshing technique to simulate the position and rotational trajectories of a spherical Janus particle with unequal zeta potential on each hemisphere. We utilized the simulated position and orientation observations to assemble a 3D potential energy landscape from these multivariable data, a process required for a TIRM experiment. Finally, we used these energy landscapes to calculate system properties, including hemisphere zeta potential, particle size, and the solution’s Debye length. The described simulation tool, associated results, and data analysis is the first step in methodology development for applying TIRM to an anisotropic particle.

Consider the particle of radius a and separation distance h shown in Fig. 1(a). The particle is nearby, but not adhered to, the substrate and will remain mobile at conditions where there is a strong electrostatic repulsion between it and the neighboring substrate. The electrostatic repulsive force depends on the solution’s Debye length, the size of the particle, and the particle’s and substrate’s Stern potential, which is typically equated with the zeta potential of each surface. In addition to electrostatic repulsion, we also consider the weight of the particle because experiments typically are conducted on polystyrene or silica particles that have a density not matched with that of water. The colloidal interaction force Fc, which is equal to the negative gradient of the interaction energy, between a spherical particle and a flat plate is given by

Fc=dϕc(h)dh=κBexpκhElectrostatic43πa3ρpρfgGravitational,
(1)
B=64πε0εfakTe2tanheζs4kTtanheζp4kT,
(2)
κ=2e2Cε0εfkT,
(3)

where κ is the solution’s Debye parameter, B is the electrostatic charge parameter, ρp and ρf are density of the particle and substrate, respectively, g is the gravitational acceleration, ε0 is the electric permittivity of vacuum, εf is the relative permittivity of water, e is the charge on an electron, ζs and ζp are the zeta potentials of the surface and particle, respectively, and C is the electrolyte concentration in the bulk. Note that the expressions above can be directly applied to a sphere of uniform zeta potential, but not for the Janus sphere shown in Fig. 1(a) because of the additional dependence of Fc on θ′. We later describe a meshing technique to account for such a dependence.

FIG. 1.

(a) Illustration of a spherical Janus particle with radius a and hemispheres of zeta potential A, ζA, and B, ζB. The particle is separated a distance h from the nearby boundary. θ′ rotation is defined as in the axis parallel to the boundary, while φ rotation is defined as in the axis perpendicular to the boundary. (b) Schematic representation of orientation of a spherical Janus particle. The angle θ is defined between the center of side A and the axis perpendicular to the substrate, such that when the center of side “B” is downward, θ = 0°, and when the center of side “A” is downward, θ = 180°. θ = 90° occurs when the center of the transition is pointed downward. (c) Particle with surface meshing to account for non-uniform zeta potential. The projected surface from each mesh point was used to calculate the interaction force between the sphere and substrate, which is similar to the Derjaguin approximation.

FIG. 1.

(a) Illustration of a spherical Janus particle with radius a and hemispheres of zeta potential A, ζA, and B, ζB. The particle is separated a distance h from the nearby boundary. θ′ rotation is defined as in the axis parallel to the boundary, while φ rotation is defined as in the axis perpendicular to the boundary. (b) Schematic representation of orientation of a spherical Janus particle. The angle θ is defined between the center of side A and the axis perpendicular to the substrate, such that when the center of side “B” is downward, θ = 0°, and when the center of side “A” is downward, θ = 180°. θ = 90° occurs when the center of the transition is pointed downward. (c) Particle with surface meshing to account for non-uniform zeta potential. The projected surface from each mesh point was used to calculate the interaction force between the sphere and substrate, which is similar to the Derjaguin approximation.

Close modal

A TIRM experiment tracks the stochastic fluctuations of a particle in the direction normal to the nearby substrate [z-axis in Fig. 1(a)]. The translational and rotational diffusion coefficients, Dz and Drθ, respectively, of a particle decrease when it approaches a boundary due to an increase in the hydrodynamic drag force. The increase in translational hydrodynamic drag for the slow motion of a sphere approaching a wall is well-known and has been confirmed experimentally.12,29,30 Note that both normal and rotational mobility will decrease due to the wall. The translational diffusion coefficient for a spherical particle in the direction normal to the boundary is

Dz=kTfq(h)=kT6πηaq(h),
(4)
qh=6h2+2ha6h2+9ha+2a2,
(5)

where f is the linear quasi-steady Stokes drag coefficient f=6πηa, η is the fluid viscosity, and q(h) is the wall correction factor.

The exact analytical solution of the Stokes equation for the rotation of a sphere around an axis parallel to a substrate in a creeping flow when the sphere-wall gap is very small has been derived by Dean and O’Neill.31 Liu and Prosperetti32 reported an approximate solution by Maude33 for this analytical problem close to the boundary. Goldman34,35 plotted graphically three different approximations and an exact (corrected Dean and O’Neill) result of the Dean and O’Neill solution. As in our system the dimensionless gap range is between 0.015 and 0.15, the Maude approximation solution is not applicable. To use exact numerical solution (corrected Dean and O’Neill) in our range, we fitted these data in the appropriate range. The fitted expressions are

Dr,θ=kTfr,/qθ(h)=kT8πηa3/qθ(h),
(6)
qθh=0.9641ha0.1815.
(7)

The motion of a Brownian particle very near a boundary is described by a force balance on the particle with the inclusion of a stochastically fluctuating force. At small Reynolds number, the inertial term in the force balance is zero and there is a well-known numerical solution of the Langevin equation by Ermak and McCammon.36 The trajectory of a single particle was predicted via

ht+Δt=ht+dDzdhΔt+DzkTFΔt+HΔt.
(8)

In this algorithm, h is calculated in each forward time step and HΔt is a Gaussian random height displacement variable satisfying H=0. The variance of this noise is H2=2DzΔt. These expressions satisfy the dissipation theorem that defines particle-fluid viscous dissipation as the origin of particle fluctuations in the suspending fluid.37 Note that Eq. (8) allows for the friction coefficient to depend on height. Previous work has shown that if height dependency of the friction coefficient is assumed to be independent of the separation height, there will be a systematic error between the real and simulated potential energy profiles.28 Two conditions must be satisfied to use the stepping algorithm. One is that the time step is sufficiently short such that force, height, and all other calculated properties are constant. The second condition is that the time step is much longer than the momentum relaxation time of the particle.

The rotational trajectory was predicted with an expression analogous to Eq. (8). The rotational stepping algorithm is38,39

θt+Δt=θt+Dr,θTexinΔtkT+GΔt.
(9)

The GΔt term is stochastic Gaussian standard Brownian rotational displacement. This random noise is defined by the following variance G2=2Dr,θΔt. The Texin term accounts for the time-dependent internal and external torques. External torques may include contributions from both electrostatics and gravity. Electrostatic torque will arise because of the mismatch in electrostatic repulsion very near the boundary between hemispheres A and B. Gravitational torque will arise because of the mismatch in densities of the base particle and cap material. Common examples of cap materials on polystyrene (ρ = 1.055 g/cm3) or silica (ρ = 2.65 g/cm3) particles are gold (ρ = 19.32 g/cm3),40 silica,41 and carbon.42 These phenomena will induce torques that will make the middle term on the right-hand side of Eq. (9) non-zero. We conducted simulations to isolate the effects of either electrostatic or gravitational torque. We found that electrostatic torque had a negligible effect, with an average root mean squared error of 0.0642 kT when compared with results neglecting the effect. However, we did find that gravitational torque had a significant impact on the rotational sampling of a Janus particle and, consequently, the potential energy landscape. Thus, the cap weight and thickness (even for very thin caps) will affect the potential energy landscape. Herein, we assumed that the time dependent internal and external torques are zero and, thus, our results are applicable only to systems where the cap density matches that of the base particle. Given the very large parameter space, including cap density and thickness, we plan on exploring the influence of cap weight on rotational dynamics of a Janus particle in a future contribution. All simulations presented in this article were for a polystyrene particle with a cap of matching density.

Note from Fig. 1(a) that the force felt by the Janus particle will be axisymmetric about φ, the axis perpendicular to the substrate. We utilized this symmetry to make the simulation more efficient by only stepping orientation fluctuations in θ′ and ignoring rotation of the particle in the φ direction. Equations (8) and (9) were executed for either 2.4 × 106 or 4.8 × 106 time steps for a single simulation, with each time step Δt = 5 ms, unless otherwise noted. A single simulation provides trajectories in both h and θ′, which were then used to prepare histograms as described later in this article. Subsequently, potential energy profiles were calculated from these histograms. The potential energy profiles and landscapes presented later in this article were obtained from averaging ten simulations. At each bin height, the mean and standard deviation of potential energy were calculated and compared with the analytical potential energy profile. The convergence criteria for each height is satisfied if the standard deviation and difference between the meshing method simulation (mean point of 10 sets of simulations) and true potential energy profile are less than 0.3 kT. Note that this convergence criterion is similar to that recommended by Sholl et al.28 

The stepping algorithms described above were implemented for a spherical Janus particle of different zeta potentials on each hemisphere. The sphere was meshed into small parts in both the polar angle θ′ and azimuthal angle φ [see Fig. 1(c)] to account for the change in zeta potential upon Janus particle rotation in the θ′ direction. After projecting the small curved surfaces parallel to the substrate, the electrostatic double layer force between each projected flat surface and the substrate was calculated. In the case of monovalent salt, the electrostatic double layer repulsion pressure of two flat parallel surfaces is

PdL(i)=64CkTtanheζs4kTtanheζmpi4kTexp(κhdL(i)),
(10)

where ζs and ζmp are the zeta potential of the substrate and flat projected surface of the mesh point i. The electric double layer repulsion force between each mesh point and the substrate is proportional to the product of the pressure and the projected surface area for each mesh point,

FdL(i)=PdL(i)A(i).
(11)

The effective height (hdL) in the pressure equation was equal to the mid-point of the mesh surface height,

hdL=hp+L/2.
(12)

Finally, the total electrostatic double layer repulsion force between the particle and the substrate was calculated by summing the forces at each mesh point,

FdL_total=FdL(i).
(13)

Once implemented, the meshing method together with the Brownian Dynamics simulation allows for a direct prediction of particle height and orientation evolving in time very near a boundary for a Janus particle of unequal zeta potential.

The optimum mesh size was determined by comparing the electrostatic portion of the potential energy profile as calculated from Eq. (1) (i.e., non-meshed particle) to that calculated from Eq. (13) (i.e., meshed particle) for different mesh sizes, but with uniform zeta potential. The gravitational portion of the force calculation was unnecessary for this comparison because the gravitational force depends on neither the zeta potential of the particle nor the height of the mesh point. We found that increasing the number of mesh points decreased the error in the comparison between a non-meshed and meshed particle. A particle with 30 × 30 meshing (900 total points with 30 points in θ′ and 30 points in φ) had nearly a 22% error when compared with a non-meshed particle with the same properties. A particle with 100 × 100 meshing had less than 1% error when compared with a non-meshed particle with the same properties. Table I summarizes these results. We chose a particle with 100 × 100 meshing to have acceptable numerical error for our purposes as using more mesh points severely decreased computational speed, yet did not appreciably increase accuracy.

TABLE I.

Error associated with varying discretization in meshing. Comparison of double-layer repulsion force for optimum mesh size calculations.

Non-meshed particleMeshed particle
Mesh sizeelectrostatic force (N)electrostatic force (N)Error%
20 × 20 9.6898 × 10−16 1.4158 × 10−15 46.114 
30 × 30 9.6898 × 10−16 9.6898 × 10−16 22.013 
60 × 60 9.6898 × 10−16 1.0029 × 10−15 3.504 
80 × 80 9.6898 × 10−16 9.8595 × 10−16 1.751 
100 × 100 9.6898 × 10−16 9.7851 × 10−16 0.983 
120 × 120 9.6898 × 10−16 9.7458 × 10−16 0.577 
140 × 140 9.6898 × 10−16 9.7224 × 10−16 0.335 
Non-meshed particleMeshed particle
Mesh sizeelectrostatic force (N)electrostatic force (N)Error%
20 × 20 9.6898 × 10−16 1.4158 × 10−15 46.114 
30 × 30 9.6898 × 10−16 9.6898 × 10−16 22.013 
60 × 60 9.6898 × 10−16 1.0029 × 10−15 3.504 
80 × 80 9.6898 × 10−16 9.8595 × 10−16 1.751 
100 × 100 9.6898 × 10−16 9.7851 × 10−16 0.983 
120 × 120 9.6898 × 10−16 9.7458 × 10−16 0.577 
140 × 140 9.6898 × 10−16 9.7224 × 10−16 0.335 

Finally, the orientation and height trajectories of the particle were used to calculate the potential energy landscape. The probability of finding a particle at each height depends on the potential energy at that location. Boltzmann’s equation is given by

ph=Aexp(ϕhkT),
(14)

where phdh is the probability of finding the particle between h and h + dh and ϕh is potential energy at that height. The shape of the histogram of particle observations n(h) is the same as the shape of the probability density function ph.5,6 In other words, n(h) is directly proportional to p(h) and the potential energy profile can be deduced by inverting Eq. (14) to give the following equation:

ϕchϕchmkT=lnnhmnh,
(15)

where ϕchm is the potential energy of the particle in the most probable height, nhm is the number of particle observation at the most probable height and nh is the number of particle observation at that height. Comprehensive details of analyzing data associated with TIRM can be found elsewhere.5 We included the effect of orientation in assembling the potential energy landscape by calculating a potential energy profile based on height observations (Eq. 15) at each theta interval. In other words, the potential energy profile at each theta was computed assuming a Boltzmann distribution of height at each of those theta. The potential energy profiles at each theta were then assembled to obtain the potential energy landscape for a Janus particle.

We validated our method by conducting Brownian dynamic simulations of a meshed spherical Janus particle with zeta potential across each hemisphere assumed to be equal (ζA = ζB = ζS = −50 mV) and subsequently compared the potential energy profile to the known analytical profile of an isotropic sphere. Figure 2 shows the comparison between the potential energy profile obtained by dynamics and the analytical profile. The simulated data shown in Fig. 2 are those which converged to >6 kT based on the criteria described in Sec. II. These converged data demonstrate the validity of our meshing and stepping algorithms.

FIG. 2.

Potential energy profiles of a meshed spherical Janus particle with equal zeta potential across each hemisphere (red circles) and the known analytical potential profile from Eq. (1) of an isotropic sphere (blue line). All data reported from our simulations have converged to within 0.3 kT of the analytical values. The converged data reach a maximum potential energy >6 kT, which is the threshold typically used for TIRM experiments. Note that only every fourth point of the simulated profile was reported for clarity.

FIG. 2.

Potential energy profiles of a meshed spherical Janus particle with equal zeta potential across each hemisphere (red circles) and the known analytical potential profile from Eq. (1) of an isotropic sphere (blue line). All data reported from our simulations have converged to within 0.3 kT of the analytical values. The converged data reach a maximum potential energy >6 kT, which is the threshold typically used for TIRM experiments. Note that only every fourth point of the simulated profile was reported for clarity.

Close modal

Once validated, we used our method to simulate the dynamics of a Janus particle of unequal zeta potential and subsequently calculate the potential energy landscape in both h and θ. A Janus particle of radius 3 μm, ζA = −5 mV, and ζB = −50 mV was simulated for 4.8 × 106 time steps. The raw simulated data consisted of the height and orientation of the particle at each time step. Figure 1(b) shows the convention we use to define hemisphere orientation with respect to the angle θ. In all of the simulations, the initial orientation of the particle was set to π/2. Potential energy landscapes were calculated from 10 simulations.

1. Janus particle rotational and height fluctuations

Figure 3 shows the rotational fluctuations experienced by the Janus particle for a representative simulation. The particle undergoes stochastic fluctuations because of rotational diffusion, rotating between ∼0 and ∼20 radians. The time scale for rotational diffusion of a micrometer scale particle scales with ∼a−3 (rather than ∼a−1 as with Dz) and is quite slow, thereby producing rotational fluctuations over only ∼3 revolutions in 6 h of simulated time.

FIG. 3.

Rotational fluctuations from a representative simulation of a Janus particle with two different zeta potentials (ζA = −5 mV, ζB = −50 mV) for 4.8 × 106 time steps.

FIG. 3.

Rotational fluctuations from a representative simulation of a Janus particle with two different zeta potentials (ζA = −5 mV, ζB = −50 mV) for 4.8 × 106 time steps.

Close modal

Figure 4 shows the 3D histogram, as well as cut away panels of the 2D histograms, for height and orientation. The bin size in height was equal to Δh = 1 nm, and the bin size in orientation was equal to Δθ = 2°. The histogram of orientation comprising bins between 0° and 180° utilizes the rotational symmetry of the particle in two ways. First, all rotational fluctuations were converted to a confined range, between 0° and 360°, because of the 2π orientation repetition of the particle; for instance, θ = 90° is identical to θ = 450°. Second, the particle has rotational symmetry such that it is only necessary to consider orientation positions sampled between θ = 0° and θ = 180°. There is a peak in height at each orientation, which is the maximum number of observations or the location of the most probable height at that orientation. The location of the most probable height can also be seen in Fig. 4(a). Note that the observable peak in Fig. 4(a) is not sharp (as compared to an isotropic particle) because the most probable height at each orientation, which differs, is being displayed in 2D.

FIG. 4.

Histogram of position and orientation observations. (a) 2D view of the histogram that shows just the height and number of observations. (b) 3D histogram in both height and orientation. (c) 2D view of the histogram that shows just the orientation and number of observations.

FIG. 4.

Histogram of position and orientation observations. (a) 2D view of the histogram that shows just the height and number of observations. (b) 3D histogram in both height and orientation. (c) 2D view of the histogram that shows just the orientation and number of observations.

Close modal

A distinct difference between Figs. 4(a) and 4(c) is that there is a clear height preference, while there is no clear orientation preference. One may expect there to be an orientation preference because of the change in mobility associated with lower heights sampled at small zeta potentials. For instance, rotational diffusion is further hindered when the particle rotates towards the lower zeta surface and samples lower heights. We tested this hypothesis by tabulating the orientation preferences for varying zeta conditions. Over a set of ten simulations, there was no clear preference between the hemispheres with a smaller zeta potential, as compared to the hemisphere with a larger zeta potential. However, after averaging those ten simulations, there was a very small preference for the hemisphere with the smaller zeta potential to be orientated towards the substrate. Table II summarizes the preference in zeta potential values. Given that there is no direct dependence of orientation fluctuations on the potential energy [via Eq. (9)], one would not expect there to be a most probable orientation. However, the very weak dependence of rotational diffusion on height could be the origin of the minor preference shown in our results summarized in Table II. Note that we expect the presence of a cap with a density not matching that of the base particle to profoundly impact these rotational sampling results.

TABLE II.

Number of particle observation in various positions.

No. of observationNo. of observation
ζA (mV)ζB (mV)A downB down
20 2 444 485 2 355 515 
50 2 424 742 2 375 258 
10 60 2 524 228 2 275 772 
100 2 480 467 2 319 533 
No. of observationNo. of observation
ζA (mV)ζB (mV)A downB down
20 2 444 485 2 355 515 
50 2 424 742 2 375 258 
10 60 2 524 228 2 275 772 
100 2 480 467 2 319 533 

2. Calculating and interpreting the potential energy landscape

Simulated observations of height were used to compute the potential energy profile of the particle at each individual orientation [Eq. (15)]. At each orientation, the maximum number of observations occurs at the most probable height, which appears as nhm in Eq. (14). The bin size of the height was set Δh = 1 nm, and the bin size in orientation was equal to either Δθ = 1° or 2°. Figures 5(a) and 5(b) show the potential energy landscape calculated from simulated data. Note that these reported data satisfy two conditions for convergence: (1) the standard deviation from 10 sets of potential energy landscapes should be less than 0.3 kT at each point and (2) the difference between the average of 10 sets of potential energy landscapes and the analytical potential energy landscape should be less than 0.3 kT at each point. The analytical potential energy landscape was calculated by fixing the orientation in a specific position, meshing the particle, and then using Eqs. (1) and (10)–(13). This process was repeated for all desired heights to obtain the landscape. An analytical potential energy landscape for standard simulation conditions is shown in Figs. 5(c) and 5(d). This procedure provided minimum and maximum heights for which the simulation converged in each orientation.

FIG. 5.

Potential energy landscape for a Janus particle with ζA = −5 mV and ζB = −50 mV. (a) 3D and (b) 2D views of potential energy landscape. The height and orientation bin sizes were 1 nm and 2°, respectively. The potential energy landscape displays the expected expo-linear profile, but with a transition at θ = 90° at the location of the boundary between ζA and ζB. (c) 3D and (d) 2D views of analytical potential energy landscape for a Janus particle with ζA = −5 mV and ζB = −50 mV.

FIG. 5.

Potential energy landscape for a Janus particle with ζA = −5 mV and ζB = −50 mV. (a) 3D and (b) 2D views of potential energy landscape. The height and orientation bin sizes were 1 nm and 2°, respectively. The potential energy landscape displays the expected expo-linear profile, but with a transition at θ = 90° at the location of the boundary between ζA and ζB. (c) 3D and (d) 2D views of analytical potential energy landscape for a Janus particle with ζA = −5 mV and ζB = −50 mV.

Close modal

Figure 6(a) shows that the potential energy profile at 0°θ°<2° is shifted to larger heights in comparison to orientation at 178°θ°<180° in Fig. 6(c). The reason for the shift to larger heights at 0°θ°<2° is that at this orientation hemisphere B, with ζB = −50 mV, is wholly exposed to the surface with no contribution from hemisphere A. As a consequence of the relatively stronger electrostatic repulsive force, the particle samples larger heights with a larger most probable height. Conversely, the electrostatic repulsion force at 178°θ°<180° is orientated such that hemisphere A, with ζA = −5 mV, is wholly exposed with no contribution from hemisphere B. This orientation experiences a relatively weaker electrostatic repulsive force, and the particle samples smaller heights at this orientation.

FIG. 6.

The potential energy deduced from simulation data of a 3 μm radius polystyrene Janus sphere with ζA = −5 mV and ζB = −50 mV in 1 mM NaCl electrolyte (open blue circles). The simulation conditions are 4.8 × 106 iterations and Δt = 5 ms. The selected orientation bin size in this figure is Δθ = 2°; thus, there are 90 orientation positions. The solid red line is a fit to the data as described in the main text. Panel (a) is for θ bin 0°θ<2°, panel (b) is for θ bin 90°θ<92°, and panel (c) is for θ bin 178°θ<180°.

FIG. 6.

The potential energy deduced from simulation data of a 3 μm radius polystyrene Janus sphere with ζA = −5 mV and ζB = −50 mV in 1 mM NaCl electrolyte (open blue circles). The simulation conditions are 4.8 × 106 iterations and Δt = 5 ms. The selected orientation bin size in this figure is Δθ = 2°; thus, there are 90 orientation positions. The solid red line is a fit to the data as described in the main text. Panel (a) is for θ bin 0°θ<2°, panel (b) is for θ bin 90°θ<92°, and panel (c) is for θ bin 178°θ<180°.

Close modal

Figure 7 is a 2D height and orientation view of Fig. 5. The smallest and largest height is a potential energy slice at 3 kT. This figure illustrates the limiting and transitional behavior of the Janus particle as it rotates from large zeta (−50 mV) to small zeta potential (−5 mV). The particle experiences limiting behavior as θ approaches 0° and 180°, with each of these limits having potential energy profiles closely matching particles with a uniform zeta potential of either –50 mV (for 0°) or −5 mV (for 180°). The transition region between 80° and 100° is identifiable by the boundary found at the smallest heights, which is where the potential energy is equal to 3 kT. The boundary at 3 kT transitioned from approximately 95 nm to 80 nm, while the most probable height transitioned from approximately 110 nm to 90 nm (summarized in Table III). This signature results from the change in magnitude of the surface zeta potential as a consequence of rotation, as illustrated in the bottom panel of Fig. 7. A change in zeta potential magnitude will induce a chance in the magnitude of the electric double layer force. The resolution of the transition will be an important aspect of data analysis for TIRM experiments. These simulations show that even for a perfectly discrete boundary, there is a diffuse transition between hemisphere A and hemisphere B. Further, simulations with 3 μm polystyrene particles (results not shown) demonstrated that the shape of the transition was not a strong function of particle size. The potential energy landscape was also used to determine the appropriate time step for such a simulation. The time step was systemically changed for a Janus particle with −5 mV and −50 mV zeta potentials (rows 5–10 in Table III). Our results showed that landscapes obtained from simulations with time steps of 5 ms, 6 ms, 7 ms, and 8 ms have good agreement with analytical potential energy landscapes, while step times <5 ms and >9 ms have inaccurate results.

FIG. 7.

(a) 2D Orientation–height view of potential energy landscape. As can be seen, the position and shape of the Janus spherical particle in this simulation in the ranges of 0°θ<2°, 90°θ<92°, and 178°θ<180° are shown. The potential is sliced at 3 kT and the upper and bottom boundaries represent this slice. There is a distinct transition when rotation through past 90°. (b) The 2D view of the analytical potential energy profile for standard simulation conditions and bin size equal to 2°. The potential is sliced at 3 kT and the upper and bottom boundaries represent this slice.

FIG. 7.

(a) 2D Orientation–height view of potential energy landscape. As can be seen, the position and shape of the Janus spherical particle in this simulation in the ranges of 0°θ<2°, 90°θ<92°, and 178°θ<180° are shown. The potential is sliced at 3 kT and the upper and bottom boundaries represent this slice. There is a distinct transition when rotation through past 90°. (b) The 2D view of the analytical potential energy profile for standard simulation conditions and bin size equal to 2°. The potential is sliced at 3 kT and the upper and bottom boundaries represent this slice.

Close modal
TABLE III.

Parameter fits for various conditions.

Real dataFitting parameters
ζA (mV)ζB (mV)Δt (ms)Orientation bin size (deg)Hm(A) (nm) error%Hm(B) (nm) error%k (nm) error%r (μm) error%ζA error%ζB error%
20 79.89 102.3 0.1071 3.038 2.81 23.38 
     0.12 0.33 2.88 1.2 40 17 
10 60 95.73 112.07 0.1049 3.05 11.6 65.89 
     0.3 0.35 0.7 1.6 16 
100 89.16 115.38 0.1062 3.035 6.34 141 
     0.4 0.26 1.9 1.1 20 41 
50 88.61 110.4 0.1074 3.071 7.73 80.23 
     0.21 0.18 3.31 2.3 54.6 60.4 
50 89.07 110.81 0.1045 3.046 5.86 52.48 
     0.3 0.54 0.4 1.5 17 4.8 
50 88.18 95.05 0.1071 2.96 10.18 85.87 
     0.6 13.74 2.9 1.2 116 71.7 
50 89.22 110.78 0.1042 3.038 5.64 51.72 
     0.48 0.53 0.12 1.29 12.8 3.4 
50 89.3 110.82 0.1044 3.035 5.41 52.28 
     0.57 0.56 0.29 1.17 8.3 4.5 
50 89.52 110.94 0.1027 3.045 4.87 45.67 
     0.8 0.67 1.3 1.5 2.5 8.6 
10 50 10 89.59 111.9 0.102 3.039 4.7 41.56 
     0.9 0.8 2.06 1.33 5.8 16.8 
Real dataFitting parameters
ζA (mV)ζB (mV)Δt (ms)Orientation bin size (deg)Hm(A) (nm) error%Hm(B) (nm) error%k (nm) error%r (μm) error%ζA error%ζB error%
20 79.89 102.3 0.1071 3.038 2.81 23.38 
     0.12 0.33 2.88 1.2 40 17 
10 60 95.73 112.07 0.1049 3.05 11.6 65.89 
     0.3 0.35 0.7 1.6 16 
100 89.16 115.38 0.1062 3.035 6.34 141 
     0.4 0.26 1.9 1.1 20 41 
50 88.61 110.4 0.1074 3.071 7.73 80.23 
     0.21 0.18 3.31 2.3 54.6 60.4 
50 89.07 110.81 0.1045 3.046 5.86 52.48 
     0.3 0.54 0.4 1.5 17 4.8 
50 88.18 95.05 0.1071 2.96 10.18 85.87 
     0.6 13.74 2.9 1.2 116 71.7 
50 89.22 110.78 0.1042 3.038 5.64 51.72 
     0.48 0.53 0.12 1.29 12.8 3.4 
50 89.3 110.82 0.1044 3.035 5.41 52.28 
     0.57 0.56 0.29 1.17 8.3 4.5 
50 89.52 110.94 0.1027 3.045 4.87 45.67 
     0.8 0.67 1.3 1.5 2.5 8.6 
10 50 10 89.59 111.9 0.102 3.039 4.7 41.56 
     0.9 0.8 2.06 1.33 5.8 16.8 

3. Determining particle and solution properties from potential energy landscape

One outcome of a TIRM measurement is the determination of particle and solution properties from the position and orientation observations and ultimately the potential energy landscape. The simulation described herein is a tool for developing a methodology for determining particle properties from a TIRM measurement on an anisotropic particle. Herein, we fit the potential energy landscapes for a variety of conditions with a non-linear least square method to determine particle and solution properties, which are summarized in Table III. The potential energy landscape was used to determine particle and solution properties, including particle radius a, hemisphere A zeta potential ζA, hemisphere B zeta potential ζB, the solution’s Debye length κ1, and the most probable height hm, as would be done in a TIRM experiment.

The potential energy profile at every orientation (including in the transition region) was fit to determine particle radius and Debye length, while the potential energy profiles only at the limiting conditions were used to calculate the most probable height and zeta potentials. Determining particle and solution properties from a fit of the potential energy landscape requires treating different regions of the landscape because of variations in sensitivity of parameters to the fit. Particle radius and the solution’s Debye length are robust to the fitting process because small changes in each of those parameters will produce significant changes in the potential energy profile. However, the particle’s zeta potential appears only in the electrostatic parameter B, which itself is a pre-factor to the exponential electrostatic repulsion. Although large changes in zeta potential at our conditions will induce large changes in B, these large changes in B do not produce significant changes in the potential energy profile. This issue was exacerbated at the transition region, where there are contributions from both hemispheres to the electrostatic repulsion. Consequently, the electrostatic part of the interaction is known to be difficult to determine, even for isotropic particles.5 Thus, once the transition region was identified, we fit the potential energy profile only at the limiting conditions as θ → 0° or θ → 180°. For example, among 90 orientation conditions (with Δθ = 2°), the average of first 40 orientations was used to calculate the most probable height and zeta potential of the B hemisphere and the average of last 40 orientations was used to calculate hemisphere A parameters.

Based on our results summarized in Table III, we determined that it would be possible to measure properties of the Janus particle as shown by the good agreement between fitting parameters and real system parameter values (rows 1-3 and 5 in Table III). Figure 6 shows a sample of fitting results for all three orientations. Figure 6(a) shows the orientation 0°θ<2° corresponding to hemisphere B (ζB = −50 mV) of the Janus particle on the bottom side, and Fig. 6(c) shows the orientation 178°θ<180° corresponding to hemisphere A (ζA = −5 mV) of the Janus particle is on the bottom side.

The accuracy of the parameter fits was impacted by the choice of bin size in orientation. Rows 4 and 5 of Table III show that the zeta potential fits are more accurate when choosing Δθ = 2° as compared to Δθ = 1°. The origin of this difference is in the number of observations in each bin that is necessary for an accurate fit. Reducing the orientation bin size from Δθ = 2° to Δθ = 1° reduces the number of observations binned into each orientation. The reduction in observations at each orientation reduces the accuracy of the potential energy profile calculated for that orientation. Although operating with larger orientation bin sizes will assist in determining parameters, larger orientation bins reduces the resolution of the transition region, which may be necessary for determining the patterning on the surface and will become especially important for particles with more elaborate patterning or diffuse boundaries. Ultimately, these results illustrate that there is a balance between fit parameters and potential energy landscape resolution when attempting to resolve non-uniformities along the particle surface, which is one significant difference between a TIRM experiment with an isotropic particle versus one with an anisotropic particle. Further, far more observations will be required given the need to resolve potential energy in both h and θ for an anisotropic particle, as compared with just h for an isotropic particle. A typical TIRM experiment on an isotropic particle requires ∼105 observations, but our simulation results suggest that ∼106–107 observations will be required to calculate a potential energy profile comprising data of similar accuracy.

This article describes the method and results from a Brownian dynamics simulation suitable for a Janus sphere of unequal zeta potential. Our effort was motivated by the need for guidance and methodology for analyzing data from a TIRM experiment. For this purpose, we have developed a method for meshing a sphere of arbitrarily patterned zeta potential that integrates with a Brownian dynamics simulation. Following validation of the meshing method, we simulated the position and rotational trajectories of the particle near a boundary and subsequently calculated potential energy landscapes for the Janus particles. We derived for the first time a 3D potential energy profile landscape for a Janus particle as would be obtained in a TIRM system. This simulation and landscape can predict the behavior of the Janus particle and help us in the optimization of experimental parameters in an experimental TIRM system for a Janus particle. Results show that the potential energy landscape of a Janus sphere has a transition region at the location of the boundary between the two Janus halves, which depends on the relative zeta potential magnitude. The potential energy was accurately fit to obtain parameters including zeta potential magnitude in each hemisphere, particle size, minimum potential energy position and electrolyte concentration, or Debye length. Time step and orientation bin size are important parameters in the simulation and data analysis process. Additionally, simulations show that an experiment may require more than 106 observations to obtain a suitable potential energy landscape. These results demonstrated the utility of Brownian dynamics simulations as a tool to probe TIRM for anisotropic particles.

This work was supported by the Cleveland State University Office of Research Startup Fund and NSF No. 1344954.

1.
J.
Berg
,
An Introduction to Interfaces and Colloids
, 1st ed. (
World Scientific
,
2010
).
2.
H.-J.
Butt
,
B.
Cappella
, and
M.
Kappl
,
Surf. Sci. Rep.
59
,
1
(
2005
).
3.
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.
73
,
36601
(
2010
).
4.
D. C.
Prieve
and
N. A.
Frej
,
Langmuir
6
,
396
(
1990
).
5.
D. C.
Prieve
,
Adv. Colloid Interface Sci.
82
,
93
(
1999
).
6.
M. A.
Bevan
and
S. L.
Eichmann
,
Curr. Opin. Colloid Interface Sci.
16
,
149
(
2011
).
7.
S. G.
Flicker
and
S. G.
Bike
,
Langmuir
9
,
257
(
1993
).
8.
M. A.
Bevan
and
D. C.
Prieve
,
Langmuir
15
,
7925
(
1999
).
9.
M. A.
Bevan
and
P. J.
Scales
,
Langmuir
18
,
1474
(
2002
).
10.
T. D.
Edwards
and
M. A.
Bevan
,
Langmuir
28
,
13816
(
2012
).
11.
C. L.
Wirth
,
P. J.
Sides
, and
D. C.
Prieve
,
J. Colloid Interface Sci.
357
,
1
(
2011
).
12.
M. A.
Bevan
and
D. C.
Prieve
,
J. Chem. Phys.
113
,
1228
(
2000
).
13.
C.
Hertlein
,
L.
Helden
,
A.
Gambassi
,
S.
Dietrich
, and
C.
Bechinger
,
Nature
451
,
172
(
2008
).
14.
L.
Helden
,
K.
Dietrich
, and
C.
Bechinger
,
Langmuir
32
,
13752
(
2016
).
15.
S. L.
Eichmann
,
B.
Smith
,
G.
Meric
,
D. H.
Fairbrother
, and
M. A.
Bevan
,
ACS Nano
5
,
5909
(
2011
).
16.
C.
Bolton
and
R. R.
Dagastine
, in
American Chemical Society Colloid and Surface Science Symposium
,
2017
.
17.
S. C.
Glotzer
and
M. J.
Solomon
,
Nat. Mater.
6
,
557
(
2007
).
18.
Z.
He
and
I.
Kretzschmar
,
Langmuir
29
,
15755
(
2013
).
19.
L.
Hong
,
S.
Jiang
, and
S.
Granick
,
Langmuir
22
,
9495
(
2006
).
20.
A. A.
Shah
,
B.
Schultz
,
K. L.
Kohlstedt
,
S. C.
Glotzer
, and
M. J.
Solomon
,
Langmuir
29
,
4688
(
2013
).
21.
A. A.
Shah
,
B.
Schultz
,
W.
Zhang
,
S. C.
Glotzer
, and
M. J.
Solomon
,
Nat. Mater.
14
,
117
(
2014
).
22.
Y.
Liu
,
W.
Li
,
T.
Perez
,
J. D.
Gunton
, and
G.
Brett
,
Langmuir
28
,
3
(
2012
).
23.
W. L.
Miller
and
A.
Cacciuto
,
Phys. Rev. E
80
,
021404
(
2009
).
24.
Q.
Chen
,
J. K.
Whitmer
,
S.
Jiang
,
S. C.
Bae
,
E.
Luijten
, and
S.
Granick
,
Science
331
,
199
(
2011
).
25.
G. L.
Hunter
,
K. V.
Edmond
,
M. T.
Elsesser
, and
E. R.
Weeks
,
Opt. Commun.
19
,
17189
(
2011
).
26.
F. C.
Cheong
and
D. G.
Grier
,
Opt. Express
18
,
6555
(
2010
).
27.
I.
Torres-Díaz
and
M. A.
Bevan
,
Langmuir
33
,
4356
(
2017
).
28.
D. S.
Sholl
,
M. K.
Fenwick
,
E.
Atman
, and
D. C.
Prieve
,
J. Chem. Phys.
113
,
9268
(
2000
).
29.
P.
Holmqvist
,
J. K. G.
Dhont
, and
P. R.
Lang
,
J. Chem. Phys.
126
,
044707
(
2007
).
30.
M. D.
Carbajal-Tinoco
,
R.
Lopez-Fernandez
, and
J. L.
Arauz-Lara
,
Phys. Rev. Lett.
99
,
138303
(
2007
).
31.
W. R.
Dean
and
M. E.
O’Neill
,
Mathematika
10
,
13
(
1963
).
32.
Q.
Liu
and
A.
Prosperetti
,
J. Fluid Mech.
657
,
1
(
2010
).
33.
A. D.
Maude
,
Br. J. Appl. Phys.
14
,
894
(
1963
).
34.
A. J.
Goldman
,
R. G.
Cox
, and
H.
Brenner
,
Chem. Eng. Sci.
22
,
653
(
1967
).
35.
A. J.
Goldman
, “
Investigations in Low Reynolds Number Fluid-Particle Dynamics
” Ph.D. thesis, New York University, 1966; http://adsabs.harvard.edu/abs/1966PhDT........93G.
36.
D.
Ermak
and
J.
McCammon
,
J. Chem. Phys.
69
,
1352
(
1978
).
37.
M.
Hütter
and
H.
Christian Öttinger
,
J. Chem. Soc., Faraday Trans.
94
,
1403
(
1998
).
38.
B.
ten Hagen
,
S.
van Teeffelen
, and
H.
Löwen
,
Condens. Matter Phys.
12
,
725
(
2009
).
39.
E.
Dickinson
,
J. Chem. Soc., Faraday Trans. 2
81
,
591
(
1985
).
40.
S.
Gangwal
,
O. J.
Cayre
,
M. Z.
Bazant
, and
O. D.
Velev
,
Phys. Rev. Lett.
100
,
058302
(
2008
).
41.
Y.
Hotta
,
P. C. A.
Alberius
, and
L.
Bergström
,
J. Mater. Chem.
13
,
496
(
2003
).
42.
I.
Buttinoni
,
J.
Bialké
,
F.
Kümmel
,
H.
Löwen
,
C.
Bechinger
, and
T.
Speck
,
Phys. Rev. Lett.
110
,
238301
(
2013
).