Maximilian R.
Bailey
ab,
C. Miguel
Barriuso Gutiérrez
b,
José
Martín-Roca
bc,
Vincent
Niggel
a,
Virginia
Carrasco-Fadanelli
d,
Ivo
Buttinoni
d,
Ignacio
Pagonabarraga
ef,
Lucio
Isa
a and
Chantal
Valeriani
*bg
aLaboratory for Soft Materials and Interfaces, Department of Materials, ETH Zürich, Zürich, Switzerland
bDepartamento de Estructura de la Materia, Física Térmica y Electrónica, Universidad Complutense de Madrid, Madrid, Spain. E-mail: cvaleriani@ucm.es
cDepartamento de Química Física, Facultad de Química, Universidad Complutense de Madrid, Madrid, Spain
dDepartment of Physics, Institute of Experimental Colloidal Physics, Heinrich-Heine University, Düsseldorf, Germany
eDepartament de Física de la Matèria Condensada, Universitat de Barcelona, Barcelona, Spain
fUniversitat de Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona, Barcelona, Spain
gGISC - Grupo Interdiplinar de Sistemas Complejos, Madrid, Spain
First published on 29th December 2023
The underlying mechanisms and physics of catalytic Janus microswimmers is highly complex, requiring details of the associated phoretic fields and the physiochemical properties of catalyst, particle, boundaries, and the fuel used. Therefore, developing a minimal (and more general) model capable of capturing the overall dynamics of these autonomous particles is highly desirable. In the presented work, we demonstrate that a coarse-grained dissipative particle-hydrodynamics model is capable of describing the behaviour of various chemical microswimmer systems. Specifically, we show how a competing balance between hydrodynamic interactions experienced by a squirmer in the presence of a substrate, gravity, and mass and shape asymmetries can reproduce a range of dynamics seen in different experimental systems. We hope that our general model will inspire further synthetic work where various modes of swimmer motion can be encoded via shape and mass during fabrication, helping to realise the still outstanding goal of microswimmers capable of complex 3-D behaviour.
The rational design of chemical microswimmers displaying tailored motion in 3-D would greatly profit from models which can capture empirical observations.16 However, their experimental simplicity masks a complex system of chemical and mass transfer relationships, the underlying mechanisms and physics of which are still the subject of ongoing debate.23,24 It is generally accepted that a more complete description of such “chemically active colloids” requires the full solution of their phoretic fields, including details of the colloids, the substrate, and the solution composition.25–27 Unfortunately, such detailed descriptions call for a high level of technical expertise, are system-specific due to phoretic mobility parameters, have only been solved for spherical and ellipsoidal structures, and do not easily allow the inclusion of thermal fluctuations, which are critical when describing the dynamics of micron-scale objects.
To avoid the consideration of chemical phoretic fields and only account for hydrodynamic flows, the “squirmer” model is frequently invoked.28 This model was first proposed for microorganisms such as Paramecia and Volvox,29,30 and is now used as a generic description for various active systems. Squirmers swim due to a self-generated, usually stationary31 and axi-symmetric velocity field which is typically evaluated as tangential across its surface. Considered as a spherical rigid body, the squirmer can be defined using two modes describing its swimming velocity and force-dipole (B1 and B2 respectively32). Promisingly, it has been shown that a bottom-up model of Janus self-diffusiophoretic microswimmers – accounting for the unique interaction potentials of the different chemical species with the separate hemispheres of a Janus particle – will produce flow fields characteristic of spherical squirmers.33 Recent experimental studies investigating tracer flows around chemical microswimmers34 and their directed motion under flow (“rheotaxis”)35,36 have also indicated that the flow fields generated by such active agents are characteristic of squirmer-type systems. Therefore, “coarse-graining” Janus chemical microswimmers as squirmers provides a viable approach to model their behaviour despite neglecting the (albeit important) contributions arising from chemical fields26,27 and an asymmetric surface.33
A number of numerical strategies have been implemented to simulate the coarse-grained dynamics of (chemical) microswimmers, amongst which a multi-particle-collision dynamics (MPCD) description of the solvent is perhaps the most frequently invoked due to its ability to include thermal fluctuations at a lower computational cost than e.g. the Lattice Boltzmann (LB) method.32,37–46 Nevertheless, we note particularly relevant studies where LB approaches were utilised to study the behaviour of raspberry-type and other shape anisotropic microswimmers.47,48 Like MPCD, dissipative particle dynamics (DPD)49 coarse-grains the solvent as point-like fluid particles (packets of fluid molecules), and thus inherently includes the effects of thermal noise while solving the Navier–Stokes equations.50 By utilising a particle-based approach to model the solvent, the computational cost is significantly reduced, allowing the simulation of hydrodynamic interactions in systems where advection dominates over diffusion (high Péclet number, Pe).51 In DPD, the solvent particles themselves act as local thermostats due to their competing stochastic and dissipative pair-wise terms, conserving momentum and thus accounting for hydrodynamics and thermal fluctuations.52 Additionally, DPD makes use of softer inter-particle potentials, allowing greater timesteps compared to MPCD, thereby enabling simulations of larger systems at longer timescales.53,54 Therefore, DPD emerges as a suitable candidate to model microswimmers as squirmers, as it can handle hydrodynamics in high Pe systems (where directed motion dominates over diffusion) with (relatively) large time-steps, while properly dealing with thermal fluctuations. For these reasons, a DPD framework capturing the hydrodynamic interactions of microswimmers in the presence of confining boundaries – by applying tangential solvent forces around the swimmer – was recently developed by some of the authors.55 As DPD is already implemented in the open-access molecular dynamics program Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS),56 this provides the opportunity to exploit a range of in-built functions to extend previous studies, with the goal of mimicking the behaviour of chemical microswimmers above a substrate.
Here, we further develop this model to consider the influence of mass and shape asymmetries on the dynamics of spherical microswimmers in the presence of a bounding substrate, mirroring common experimental conditions for chemical active colloids.55,57–59 The strength of our approach lies in its modularity, which allows the simple inclusion of these asymmetries that may have otherwise presented significant complications when considering other numerical schemes. We find that the interplay between hydrodynamic interactions, gravity, bottom-heaviness, and shape is sufficient to qualitatively capture the 2- and 3-D physics of a range of catalytic (and photo-catalytic) Janus microswimmer systems,17,18,20i.e. while neglecting contributions from chemical59 and light60,61 gradients. Specifically, there is a competition between the hydrodynamic attraction to the substrate that grows with the swimming speed and the active force required to overcome gravity, which determines the ability of the microswimmer to enter the bulk. This balance can be furthermore adjusted by introducing mass or shape asymmetry to the particles. Our coarse-grained approach thus opens the door to the informed design of chemical microswimmers whose dynamics can be encoded via shape or mass during fabrication, helping to realise the still outstanding goal of active colloids capable of truly 3-D dynamics.
Fig. 1 (a) Schematic representation of a 2D section of the colloid using a raspberry model with cap-front (CF) mass imbalance and shape asymmetry. The colloid is built from 18 DPD filler particles (red and blue), distributed on the surface of a sphere radius Rin = 1, each with a DPD cut-off for interactions with other DPD particles (e.g. solvent particles), resulting in an overall effective microswimmer radius of R = 2. Shape asymmetry is implemented by introducing an additional particle Pshape (green) displaced a distance Rshift from an off-axis filler particle. Here, the equatorial 2D slice containing the shape asymmetry is shown, thus only depicting the 8 filler particles within this section. The thruster particle (not visible) generates the force field (light purple arrows), here consistent with a pusher-type squirmer, emanating from the center of mass of the colloid and fixed with its internal frame of reference. (b–d) Snapshot of the full simulation box (b), as well as the definition of the angles θx,y,z (c) and the azimuthal angle Φ (d) Graphics presented here and elsewhere were generated in part using the visualisation software Ovito.63 |
The microswimmers motility originates from one central “thruster” particle, which generates the squirmer-type force field experienced by the solvent particles in the spherical shell between the outer surface of the microswimmer (dashed black circle) and RH=4 (dashed black circle and dotted black circle respectively, see Fig. 1a). These solvent particles in turn apply an equal and opposite reaction force to their nearest colloid “filler” particle (18–20 DPD particles, depending on whether or not shape asymmetry is introduced to the microswimmer, see Fig. 1a), resulting in the net self-propulsion force of the microswimmer . The mass of specific filler particles can be adjusted, allowing us to introduce mass asymmetry, and thus model microswimmers equipped with heavy “caps” (e.g. Pt-coated catalytic microparticles). We define mass asymmetry as (where mhemisphere and mcap are the mass of a particle hemisphere and the heavy cap respectively), noting that the motion can thus be defined with the heavier cap at the back (CB) or at the front (CF) of the swimmer with respect to the swimming direction. Mass imbalance is thus introduced by increasing the mass of the particles constituting the cap (blue particles, see Fig. 1). In all cases, we use a squirmer active stress parameter β = B2/B1 = −2.5, i.e. a weak pusher, based on qualitative matching of simulated trajectories to the behaviour described in.20 We note that this squirmer parameter corresponds to that experimentally determined by Campbell et al.34 for a similar Pt catalytic system as that studied in ref. 20, which provides further evidence for the suitability of our selection. Solvent conditions are selected to ensure that all microswimmers with radius R studied remain in the low Reynolds (Re) number regime (,62 also see ESI, Fig. S1†), where vp is the particle speed, and specifically by controlling for the kinematic viscosity ν via the solvent parameters.49 To map different simulations to experiments, we modulate the ratio of swimming velocity Vswim to sedimentation velocity Vgravity. A minimum of 2000 time steps are used for equilibrating solvent conditions in all simulations, as were periodic boundary conditions. The number of solvent DPD particles depends on the simulation dimensions, ranging from 47002 for a (20 × 20 × 20, xyz) box, to 57197 solvent DPD particles for a (16 × 16 × 40) box.
Fig. 2 Rotational dynamics of CB microswimmers as a function of scaled time, simulated to reproduce the properties of the chemical active colloids studied in ref. 64. (a) Mean-squared-angular-displacement (MSAD) of the passive colloid (no activity). (b) Introducing activity to the microswimmer dampens its rotational diffusivity, particularly about the x- and y-axes, which saturate due to the hydrodynamic coupling to the substrate. We note the scaling of the MSAD with time here is with , and is therefore quadratic. Error bars depict the standard error of the mean from 2401 frames (sub-sampled from 50000 simulation steps). |
We then map our squirmer-inspired DPD model to the experimental findings presented by Carrasco-Fadanelli and Buttinoni.20 Specifically, we adjust Vswim:Vgravity to reproduce the different conditions studied in their work, while using the same value for mass asymmetry of their Janus polystyrene spheres coated with a Pt thin film (2.8 μm with a coating of 4 nm, masymm = 1.22) with CB motion. When Vswim < Vgravity, the gravitational torque from the heavier cap aligns the particle's internal orientation axis away from the substrate (see Fig. 3b), however the particle is unable to leave the substrate due to the force of gravity (see Fig. 3a and c). As Vswim is increased and becomes larger than Vgravity, the particle is able to leave the substrate when its internal orientation axis is directed upwards (see Fig. 3d–f, and ESI, Fig. S2†), assisted by its bottom-heaviness.59 Notably, we find that we can reproduce the structure of the microswimmer dynamics previously seen experimentally17–19 (see ESI, Fig. S3†), without requiring the previously hypothesised self-shadowing effects.60,61 The two scenarios qualitatively capture the behaviour described in ref. 20, where the gravitational torque due to the mass asymmetry introduced by the denser Pt cap favours an anti-parallel orientation of the microswimmer with gravity, competing with hydrodynamic interactions and thermal fluctuations. Interestingly, by further increasing Vswim, the simulated microswimmer is once again confined to 2-D motion in the xy plane (see Fig. 3g–i). However, here the confinement arises due to the internal orientation of the microswimmer, rather than its inability to overcome the applied gravitational force. For fast swimmers, the flow fields created by the microswimmer lead to a strong hydrodynamic attraction to the substrate, effectively “trapping” the particle in 2-D (compare insets in Fig. 3a and g). This is similar to the “sliding-state" behaviour proposed by Uspal and coworkers as a result of hydrodynamic flows from self-generated phoretic gradients,10 and mirrors the experimental observations of many chemical microswimmer systems, including that discussed in Fig. 2.
Fig. 3 Dynamics of CB microswimmers, simulated to reproduce the properties of the chemical active colloids studied in.20z/R refers to the height that a microswimmer reaches above the substrate with respect to its radius R (a, d, and g). Φ is the azimuthal angle of the microswimmer (b, e, and h) as defined in Fig. 1, and v0 is the time-series of the microswimmer's instantaneous velocity (c, f, and i, with Vswim = 〈v0〉/R). The insets in (a) and (g) show the solvent velocity flow fields around the (fixed) microswimmer at the different self-propulsion forces. Here, we progressively increase the swimming to sedimentation velocity (Vswim:Vgravity) from top to bottom to qualitatively map onto the different systems studied by Carrasco-Fadanelli and Buttinoni. The first two cases (a–c and d–f) are in good agreement with the findings presented in ref. 20, while the final case (g–i) recaptures the typical dynamics observed for chemical active colloids, where phoretic particles are essentially confined to 2-D (also see Fig. 2b). |
To summarise, we find that hydrodynamic interactions and gravitational forces on a mass-asymmetric spherical microswimmer are sufficient to qualitatively capture the physics of chemical microswimmers above a substrate under different conditions. We also confirm the importance of gravity on both the translational and rotational dynamics of active colloids suspended just above a planar surface.
Fig. 4 Swimming dynamics of CF microswimmers with respect to their cap orientation, in the absence gravity and bounding substrates, as a function of introduced shape asymmetry Rshift (colour coded from blue to red with increasing Rshift). (a) Increasing divergence 〈θ〉 in the internal body axis and the observed swimming direction is measured with increasing Rshift. We note that the values of 〈θ〉 are a function of the time-step evaluated (here, , where DT is the constant obtained for a passive sphere with Rshift = 0). Inset: linear growth in 〈θ〉 with ω values fitted from the microswimmers’ MSAD (see ESI, Fig. S6†). Graphical inset: angle θ between the swimming direction and internal body axis, as defined in the text. The particle depicted has Rshift = 0.5. (b) Rotational Pèclet number as a function of Rshift for Δt = 1. Error bars either indicate the standard error of the mean, or for the fits of ω and DR are obtained from the co-variance matrix, accounting for error propagation where relevant, from 1001 frames (sub-sampled from 50000 simulation steps). |
To quantify the effect that increasing Rshift has on the dynamics of our bulk microswimmers, we calculate the angle, θ, between the vectors describing the particle's internal orientation axis, A, (governing the direction of the self-propulsion force, see Fig. 1a), and the particles displacement vector, B, (), where A is defined by 3 particles along the microswimmers body axis at time t and B = [xt+Δt − xt, yt+Δt − yt, zt+Δt − zt]. We note that the value of 〈θ〉 is dependent on the Δt over which the displacement vector is evaluated. We find a positive correlation between Rshift and 〈θ〉 (coefficient of determination = 0.953, see Fig. 4a), which we attribute to the increasing torque experienced by the microswimmer due to the solvent forces applied to Rshape at the distance Rshift. We hypothesise that the non-zero value for 〈θ〉 at Rshift = 0 arises due to the presence of thermal fluctuations from the solvent and their effect on microswimmer motion (see ESI, Fig. S4† for further discussion). The source of the growing divergence between the internal orientation of the particles and their swimming velocity can also be observed in the MSAD of the microswimmers as the trends become increasingly ballistic with larger Rshift (see ESI, Fig. S6†). Notably, the fitted angular velocity ω (ballistic component of the MSAD, MSAD = 2DRΔt + ω2Δt2, where DR is the rotational diffusion coefficient) grows with shape asymmetry, and we determine a strong linear relationship between ω and 〈θ〉 (see Fig. 4a, bottom inset – coefficient of determination = 0.981). We furthermore note a linear growth of the “rotational” Péclet () with Rshift (see Fig. 4b, coefficient of determination = 0.982), demonstrating the increasingly deterministic orientational motion of the microswimmer as shape asymmetry becomes more pronounced. Finally, we note that introducing Pshape along axes resulting in a reduced shape asymmetry with respect to A have a negligible effect on θ, indicating the importance of appropriate shape selection (see ESI; Fig. S5, S7, and S8†).
After establishing the effect of shape asymmetry on microswimmer motion in the absence of external fields and confining boundaries, we then introduce gravity and a substrate to simulate the experimental conditions in ref. 17 and 18 using masymm = 1.25. We define Vswim:Vsediment ∼ 1.5 to prevent the microswimmers from escaping too far into the bulk. Under these conditions, microswimmers with no or low shape asymmetry are unable to leave the bottom substrate due to the joint forces of gravity and hydrodynamic attraction (see Fig. 5a). However, beyond a threshold value of Rshift > 0.25, the microswimmers display 3-D motion and begin to loop into the bulk, demonstrated by the emerging tail in P(z/R). This ability to leave the substrate is reflected in the angles which the internal orientation of the microswimmer makes with the substrate, Φ. Specifically, with increasing shape asymmetry, we observe an increase in positive values in the distribution P(Φ) (see Fig. 5b), which indicates that the particle points upwards and away into the bulk, thus overcoming gravity.
Fig. 5 Probability distributions for CF microswimmers describing their out-of-plane motion P(Δz/R) and orientation P(Φ), as a function of their introduced shape asymmetry Rshift, indicated by the colour code in the legend. (a) Above a threshold (Rshift > 0.25), the microswimmers are able to leave the substrate and enter the bulk. For visualisation, Δz where Δz/R < 0.05 is set to 0, and the presence of the wall is indicated. The particle depicted at the wall has an orientation Φ ∼ −π/8. (b) The ability to leave the substrate is reflected in the emergence of a shoulder in P(Φ) for Φ > 0, as Vswim > Vgravity and therefore the microswimmer can overcome its sedimentation velocity and enter the bulk. (c) Example quasi-3D trajectory of a microswimmer with Rshift = 0.5. Note that the motion is predominantly in the xy plane, with looping out-of-plane (z) segments as described in ref. 17–19. |
We thus propose that the angular velocity ω, introduced by the shape asymmetry of the microswimmers, drives its internal swimming orientation away from the substrate, competing with the effects of gravity and hydrodynamic wall interactions (as well as thermal fluctuations). Above a certain threshold shape asymmetry – governed by the dynamics of the system – the microswimmer is able to escape the substrate and move out-of-plane.
In summary, we demonstrated that shape asymmetry is an important design parameter to control the dynamics of chemical microswimmers, in particular to promote 3-D motion. Specifically, the presence of these introduced asperities explains the unconventional swimming patterns observed in some experimental systems.17,18
Footnote |
† Electronic supplementary information (ESI) available: Further details of numerical results and extracted parameters. See DOI: https://doi.org/10.1039/d3nr03695b |
This journal is © The Royal Society of Chemistry 2024 |