Cybloids – creation and control of cybernetic colloids

Abstract

Colloids play an important role in fundamental science as well as in nature and technology. They have had a strong impact on the fundamental understanding of statistical physics. For example, colloids have helped to obtain a better understanding of collective phenomena, ranging from phase transitions and glass formation to the swarming of active Brownian particles. Yet the success of colloidal systems hinges crucially on the specific physical and chemical properties of the colloidal particles, i.e. particles with the appropriate characteristics must be available. Here we present an idea to create particles with freely selectable properties. The properties might depend, for example, on the presence of other particles (hence mimicking specific pair or many-body interactions), previous configurations (hence introducing some memory or feedback), or a directional bias (hence changing the dynamics). Without directly interfering with the sample, each particle is fully controlled and can receive external commands through a predefined algorithm that can take into account any input parameters. This is realized with computer-controlled colloids, which we term cybloids – short for cybernetic colloids. The potential of cybloids is illustrated by programming a time-delayed external potential acting on a single colloid and interaction potentials for many colloids. Both an attractive harmonic potential and an annular potential are implemented. For a single particle, this programming can cause subdiffusive behavior or lend activity. For many colloids, the programmed interaction potential allows to select a crystal structure at wish. Beyond these examples, we discuss further opportunities which cybloids offer.

---

1 Introduction

The behaviour of a many-body system is determined by the interactions between its constituents, which can range from atoms, (bio)molecules and colloids to biological cells, animals and humans or even to planets, stars and galaxies. To establish a link between the interactions and the observed behaviour is one of the central tasks of statistical physics.^1,2 For a systematic and quantitative experimental test of theoretical predictions, it is crucial that the interactions can be tuned. In atomic systems, the available set of interactions is restricted by the number of atomic species. In contrast, in colloidal systems the interactions can in principle be tuned.^3–5 To change the colloid–colloid interactions, however, one usually needs to interfere with the sample, e.g., change the chemical composition or structure of the particles or add ingredients such as salts or polymers.⁶ There are only few systems that allow for external control of the interactions. For example, the interactions of microgel particles and paramagnetic particles are susceptible to changes in temperature^7–9 and an external magnetic field,^10–15 respectively. However, the range of possible modifications remains quite limited. Irrespective of these limitations, colloids continue to play an important role in fundamental science and have had a particularly large impact on statistical physics. They have helped us to understand collective phenomena, ranging from glass formation^16–18 to phase transitions^19–23 as well as the clustering of active Brownian particles.^24,25 Furthermore, this understanding has a large impact also on other areas of sciences, ranging from life science, e.g. the crowding in biological cells²⁶ and the swarming of bacteria,^27–30 to medicine, e.g. the mechanism leading to cataract,³¹ and to material sciences, e.g. the rational design of materials such as the spontaneous assembly of hierarchical structures.^32,33 Since the success of colloidal systems originates from the specific properties (both chemical and physical) of the colloidal particles, the availability of particles with the appropriate characteristics is of prime importance.

Here we present a principle to create colloids with freely selectable properties. The key idea is to expose colloids to an externally created laser-optical potential which is programmed at wish. This allows to program the forces acting on each colloidal particle arbitrarily and thus gives rise to huge flexibility for generating novel colloidal properties. These properties might depend on several parameters, for example, previous configurations (hence introducing some memory or feedback), on the presence of other colloids (hence mimicking specific pair or many-body interactions) or a directional bias (hence changing the dynamics). Each individual particle is fully controlled and can receive commands through a predefined algorithm that can take into account any input parameter without directly interfering with the physical or chemical properties of the particles. We term these computer-controlled colloids cybloids as a portmanteau of cybernetic colloids.

The desired properties are implemented by external programming which does neither require to change nor to directly interfere with the sample. This significantly extends the range of properties that can be realised, controlled, investigated and explored to, e.g., design and fabricate previously inaccessible states. Therefore, this method provides new possibilities ranging from novel model systems to new materials.

How can colloidal particles be programmed in detail? This is realised using a control circuit (Fig. 1) using time-delayed feedback.^34–41 (The details of the experimental setup are given in Section 2.) The particles are imaged using an optical microscope and their positions, r_i(t), are followed.⁴² The actual and previous particle positions and possibly other parameters are used as input parameters to calculate the total potential, U_prog(r,t), that a particle at position r should experience at time t. The dependence of U_prog(r,t) on the input parameters can be freely chosen and determines the characteristics of the cybloids. The potential U_prog(r,t) is imposed on the particles by an appropriate light pattern I(r,t).⁴³ This procedure exploits the susceptibility of colloids to external potentials. Here, electromagnetic radiation is applied, as in optical tweezers.⁴⁴ Whereas tightly focused light beams are used in optical tweezers, extended light patterns are applied to control cybloids.⁴⁵ Modern optical methods allow us to create almost any time-varying light pattern I(r,t)⁴⁶ and hence to program a large variety of particle properties, including almost any particle–particle or particle-independent potential, cooperation rule or mode of activity.

Fig. 1 Implementation of cybloids. First, particle positions r_i are determined from the microscope image. These positions (sensor data) then serve as input parameters for the programmed potential U_prog({r_i},t) (controller) which is applied back onto the particles (system to be controlled). A (small but finite) time delay τ naturally arises, due to the time necessary for particle localization and potential calculation. Note that, while here, we constrain ourselves to the simple case of a potential dependent on the particle positions, in principle, any quantity (local density, particle trajectories, etc.) can be used in the programmed potential.

In Section 3, we illustrate this approach by first implementing a feedback potential for a single colloidal particle and a simple programmed potential that only takes the (past) particle position as an input. If the cybloid is coupled to an attractive harmonic potential centered around its past position, this leads to a transient subdiffusive regime as the particle is retracted to its past position. However, if the harmonic potential is inverted, the particle is kicked away from its previous position resulting in an intermediate ballistic dynamical regime. In our experiment, we realize such a repulsive harmonic potential by an annular ring potential⁴⁷ centered around the previous particle position. In Section 4, we consider groups of cybloids and program their interaction potential as single or multiple attractive rings around their (past) positions. By tuning the distance between different rings one can tune the crystalline lattice and generate exotic crystallites such as ones dominated by a square structure in two dimensions which is relatively uncommon for single-component systems. We compare our experimental results with theory and Brownian dynamics computer simulations and we find good agreement.

2 Materials and methods

2.1 Sample preparation and experimental methods

We use colloidal suspensions of sulfate latex beads (IDC) with radius R = 1.4 μm in heavy water (D₂O). Due to the larger density of heavy water (1.11 g ml⁻¹), the particles (1.05 g ml⁻¹) are pushed to the top by gravity as well as by radiation pressure. The suspension is stirred together with ion exchange resin to reduce the salt content. Although no further effort is taken to keep the salt concentration low, this is sufficient to avoid particles sticking to the glass surface. Subsequently, the suspension is loaded into a sample cell. The cell is constructed from a microscope slide and a cover glass (number #1) separated by two cover glasses (number #0) used as spacers that are fixed using UV glue.⁴⁸ This results in a narrow capillary which is sealed with glue after filling.

The colloidal particles are observed using an inverted microscope (Eclipse TE2000-U, Nikon) equipped with a 60× objective. Images are recorded using a CMOS camera with USB interface (MAKO U-130, AVT) and the particle locations {r⃑_i(t)} are determined in real-time with a LabVIEW subroutine (IMAQ particle analysis report VI) that is based on a standard particle tracking algorithm. Light patterns I(r,t) are created with a holographic optical tweezers setup⁴³ that is combined with the inverted microscope. It is based on a SLM i.e. spatial light modulator (LCR-2500, Holoeye), 2-D galvanometer-mounted mirrors (Quantum Scan 30, Nutfield Technology Inc.) and a diode-pumped solid-state laser (wavelength λ = 532 nm, Ventus 532-1500, Laser Quantum). A control software based on LabVIEW and MATLAB handles the image acquisition, particle tracking and hologram generation. The holograms are calculated using an iterative algorithm, the Gerchberg–Saxton algorithm⁴⁹ and uploaded to the SLM. Finally, the light field is obtained by the speckle pattern created by the SLM. The particle tracking and light field generation naturally leads to a time delay τ between the measured particle positions {r⃑_i(t)} and the introduction of the corresponding light field I(r,t). The light field is kept the same while the next one is calculated, leading to discrete potential update times nτ (where n = 0, 1, 2,…). Here, the delay time τ = 0.2 s is much smaller than the Brownian time τ_B = R²/D₀ ≈ 40 s where D₀ ≈ 0.05 μm² s⁻¹ = is the diffusion coefficient without an external potential but in the presence of the confining glass plates. This implies that the particle is still very close to its past position when the potential is introduced.

2.2 Langevin equation of motion and computer simulations

We model the time-dependent trajectory r(t) = (x(t),y(t)) of a single colloid in two spatial dimensions using Brownian dynamics. We describe the particle motion by the time-delayed over-damped Langevin equation^50,51where γ denotes the friction coefficient and F_prog(r − r_p) is the systematic force the particle experiences at its position r(t) due to the feedback potential centered around its past position r_p. f(t) is a Gaussian random force that describes the thermal motion of the particle and is characterized by its first two moments 〈f(t)〉 = 0 and 〈f(t)f(t′)〉 = 2dD₀γδ(t − t′), where 〈⋯〉 indicates an average over different realizations of the noise, D₀ is the diffusion coefficient of the particles and d = 2 gives the spatial dimension. The feedback force

F_prog(r) = −∇U_prog(r) (2)

is derived from the programmed potential U_prog(r). We use two types of potentials: a harmonic attractive potentialand a Gaussian ring potentialAs in the experiment, we define times nτ with n = 0, 1, 2,… at which the position of the particle is evaluated to update the past position r_p entering into the potential. The programmed potential U_prog(r) is then given by eqn (3) and (4) respectively and centered around the past particle position

r_p = r((k − 1)τ) (5)

for the time interval kτ ≤ t < (k + 1)τ, such that the last multiple of τ smaller than time t is kτ and the feedback potential introduced at this point corresponds to the particle positions a delay time τ prior to this which is (k − 1)τ. The feedback potential is kept the same until the next potential update at time (k + 1)τ. Consequently, the past position used in the feedback term has a time shift of between τ and 2τ to the present time t.

For the case of multiple particles, the programmed potential U_prog is composed of contributions of the form of eqn (4) centered around all past particle positions r_p,j. Here, r_p,j = r_j((k − 1)τ) indicates the past position of the jth particle. Further, we introduce an additional force to describe the direct repulsive steric interactions between the particles. This part is modelled via a smooth Weeks–Chandler–Anderson (WCA) pair-potential which takes the form

where σ is the particle diameter, r is the distance between the two particles and ε sets the interaction strength. Thus, the equation of motion in the multiple-particle case readssuch that every particle i is experiencing a fluctuating Brownian force f_i, direct interaction via the WCA potential with all particles j but itself (j ≠ i) and feedback forces including the self-term j = i.

Brownian dynamics simulations are performed based on eqn (1) and (7) with an explicit Euler forward integration scheme with a finite time step Δt = 10⁻⁴ s (all other parameters are given in the subsequent figure captions explicitly). For the simulations, the system is first equilibrated without the programmed potential (pure diffusive motion). Then, the feedback potential is introduced and the dynamics again equilibrated for ∼100 s.

3 Feedback-driven dynamics of a single colloidal particle

First we investigate the effect of a feedback potential on the dynamics of a single colloidal particle. Here, we study simple forms of the programmed potential, given by a harmonic trap (eqn (3)) and an annular ring potential (eqn (4)). In general, any interaction potential can be implemented. We shall first present an analytic solution for a harmonic trap and then turn to experiment and simulation.

3.1 Analytic solution for a harmonic feedback trap

In the following we derive an analytic expression for the mean squared displacement (MSD) of a single particle in a two-dimensional harmonic feedback potential. For the harmonic potential, the programming force is given by

F(r − r_p) = −κ(r − r_p) (8)

with the actual and past particle positions r = r(t) and r_p = r((k − 1)τ) respectively. Here, k ∈ 0, 1, 2,… is determined by the relation kτ ≤ t < (k + 1)τ, i.e. kτ is the last update time of the feedback potential before time t. For the harmonic potential, the x and y directions of motion decouple. In the following, we thus continue on to solve the one-dimensional equation of motionwith x_p = x((k − 1)τ). From this the two-dimensional MSD, follows as

〈(r(t₀ + t̃) − r(t₀))²〉 = 〈(x(t₀ + t̃) − x(t₀))²〉 + 〈(y(t₀ + t̃) − y(t₀))²〉

= 2〈(x(t₀ + t̃) − x(t₀))²〉. (10)

Here, we define the MSD as the mean squared distance the particle travels within a time t̃ starting from a fixed reference time t₀. It is thus given as the mean squared distance traveled between the fixed times t₀ and t₀ + t̃, while averaging over different realizations of the thermal noise with 〈⋯〉 denoting the thermal noise average.

We solve eqn (9) on the intervals between feedback potential updates kτ ≤ t < (k + 1)τ. Within these time intervals x_p does not change, so we use the ansatz to obtain the solution for the particle position

where the integration constant is set by x(t = kτ) = x(kτ). This solution still depends on the previous positions of the particle. For t → (k + 1)τ, we obtain the recursive formula for the particle position at the feedback potential update times

Now, from eqn (12), we can see that the thermal noise within the last delay time (from kτ to (k + 1)τ) enters with the term

while previous positions (that depend on the noise integrals at earlier times) enter with additional constant prefactors. We thus use the ansatzfor the particle position at the update times kτ, where c_n are (parameter-dependent) constants. Here, we cut the sum at k − 1, i.e. the last interval entering into the sum is for t = 0 to t = τ. This demands that the recursion according to eqn (12) ends atWe thus prescribe x(0) = 0 and x(−τ) = 0 at these two times, i.e., the particle starts from the position x(0) = 0 and the first feedback potential is centered at this position.

Next, we plug eqn (13) into both sides of eqn (12). As the noise is δ-correlated, the integrals for different time intervals are completely independent of each other and should as such have matching prefactors on both sides of the equation. We thus compare the respective prefactors of the integrals to find

Finally, we solve this recurrence equation for c_nvia the characteristic root technique to obtainand arrive at the full solution for the particle position

The mean squared displacement (MSD) for k̃τ ≤ t₀ < (k̃ + 1)τ and l̃τ ≤ t₀ + t̃ < (l̃ + 1)τ (and k̃ ≤ l̃), cf.Fig. 2, i.e. the mean squared distance the particle travels between two fixed times t₀ and t₀ + t̃ while averaging over different realizations of the thermal noise, is then given by

Here, c_n = 0 for n < 0 and τ_s = t₀ − k̃τ, τ_e = t₀ + t̃ − l̃τ are the time offsets to the respective last potential update (τ_s, τ_e ∈ [0,τ)). Finally, the abbreviation A(i,j) is defined aswhich can be solved as a geometric series asEqn (18) gives the mean squared distance that a particle travels between the times t₀ and t₀ + t̃. For sufficiently large t₀, i.e., large k̃, the effect of the initial condition disappears. In this limit, the two terms and in eqn (20) vanish (as long as ). Note, however, that eqn (18) not only depends on the time difference t̃ as well as the number of potential renewals l̃ − k̃ during that time but also on the offsets τ_s = t₀ − k̃τ and τ_e = t₀ + t̃ − l̃τ that the times t₀ and t₀ + t̃ have to the respective last renewal point of the potential. This is due to the discrete potential updates which cause the system to lose its invariance under translation in time (it is invariant only for multiples of the delay time τ). For measurements of the MSD, the resulting curve may thus vary depending on the protocol used. In particular, to obtain a full time average over t₀, the sampling of the particle position must be more frequent than the feedback update. Otherwise, i.e. for low sampling rates (= number of measurements of the particle position per delay time τ), the resulting MSD curve will not only be “rough” due to the small number of data points but actually differ from the full average over t₀, see Fig. 3.

Fig. 2 Timeline of the system. The feedback is turned on at t = 0 with the particle at x = 0 and the feedback potential centered at this position. The MSD is calculated between two later times t₀ and t₀ + t̃ that lie in the intervals k̃τ ≤ t₀ < (k̃ + 1)τ and l̃τ ≤ t₀ + t̃ < (l̃ + 1)τ with k̃, l̃ ∈ 0, 1, 2,… determined by these relations (where k̃ ≤ l̃). The feedback potential is updated at multiples of the delay time τ. The equation of motion is solved piece-wise for the intervals kτ ≤ t < (k + 1)τ, i.e. between two feedback potential update points.

Fig. 3 Analytic result for the MSD of a single particle in an attractive harmonic feedback potential and different sampling rates f_sampl (= number of measurements of the particle position per delay time τ) of the particle position. We find a decrease in the MSD compared to free diffusion (as indicated by the black line 2D₀t̃ at short times) and oscillations whose period is set by 2τ which arise as the particle is being pulled back to its past position. At long times the particle shows diffusion with a reduced diffusion constant D_{[small script l]} < D₀. Discrepancies appear in the MSD for small sampling rates as these do not cover the full average over t₀. A high sampling rate is thus necessary to obtain the “correct” MSD.

Further, we have not made any assumptions with respect to the sign of κ. The above expressions hold for both attractive (κ > 0) and repulsive potentials (κ < 0). However, for large negative values of κ, the particle continues to be accelerated away from its past position with increasing forces and the MSD diverges exponentially. An easy way to see this is from eqn (12). Subtracting x(kτ) from this equation, we find that the distance the particle travels between two potential updates

grows with a scaling factor of (on average, the additional noise term averages to zero). If this scaling factor is larger than 1 the distance increases exponentially for repeated iterations. This implies that should be smaller than 2 to ensure convergence. We therefore constrain ourselves to potential strengths that fulfill . In this case, at long times the particle dynamics becomes diffusive with the long-time diffusion coefficientwhich only depends on the strength of the external potential κ and the delay time τ.

For an attractive potential (κ > 0), the feedback potential leads to a decrease in the MSD compared to free diffusion as the particle is pulled back towards its past position, see Fig. 3 and 7. Additionally, oscillations of period 2τ appear that are linked to the potential update: after a potential update at time t = 0, the particle diffuses in a non-moving harmonic potential until the next potential update a delay time τ later. The new potential (at t = τ) is then centered around the position the particle had at the previous potential update (at t = 0), i.e. the particle experiences a force drawing it back towards its past position, leading to a reduction in the MSD. At the next potential update (at t = 2τ), the particle is then drawn towards its position at time t = τ, driving the particle away from the initial starting position and thus leading to an increase in the MSD (see also Fig. 7B). The repetition of this process leads to the observed oscillations while the additional spatial diffusion leads to overall diffusive behaviour at long times with diffusion coefficient D_{[small script l]} < D₀. This long time diffusive motion can be understood as the diffusive motion of the feedback potential center.

For a weakly repulsive potential (κ < 0, ), we find a ballistic (∝t̃²) regime at medium times between two diffusive regimes with diffusion constants D₀ and D_{[small script l]} > D₀, see Fig. 4. At short times, the particle dynamics is dominated by diffusion. On intermediate time scales, the introduced feedback potential leads to an effective particle propulsion, as the particle is pushed away from its past position. Due to spatial diffusion, the angle between current and past position changes, causing the motion to become diffusive again on long time scales but now with a larger diffusion coefficient D_{[small script l]} > D₀. Getting closer to the divergence point (), the transition from propulsion to long-time diffusion shifts to longer times. This is accompanied by the divergence of the long-time diffusion coefficient as

where .

Fig. 4 Analytic result for the MSD in repulsive harmonic feedback potentials of varying strength (average taken over t₀). The repulsive feedback potential, pushing the particle away from its past position, leads to propulsion of the particle as indicated by the ballistic ∝t̃² regime of the MSD at intermediate times. Spatial diffusion additionally leads to changes in the angle between current and past position, causing the motion to become diffusive again on long time scales.

It should also be noted that, as the memory of previous displacements decays with a factor over one decay time (see eqn (21)), this decay becomes increasingly slow towards the divergence point and the dynamics becomes highly non-local in time. An increasingly long equilibration time is needed for the effect of the initial condition to disappear. So for a system with long memory , what happens if, instead of an equilibrated feedback system, we start from t₀ = 0, i.e. the particle sitting in the centre of the feedback potential? In this case, the resulting MSD shows super-ballistic ∝t̃³ behaviour instead of ballistic motion, see Fig. 5. From eqn (18), we find that for the MSD in the limit the highest order in time is indeed given by

where we used l̃ = (t̃ + t₀ − τ_e)/τ and the limit valuesandMathematically, the t̃³ appears from the A(i,j) terms, specifically A(l̃ − 1,l̃ − 1), A(l̃ − 2,l̃ − 1), A(l̃ − 1,l̃ − 2), A(l̃ − 2,l̃ − 2) all scale with l̃³ as their highest order. Physically, the ∝t̃³ behaviour can be understood in the following way: starting from the coarse-grained particle velocity v̄_k := (x(kτ) − x((k − 1)τ))/τ, that follows from eqn (13) aswe find that the velocity average 〈v̄_k〉 = 0 is zero, while its mean squared value changes asFor long memory, i.e., close to the divergence point , the coarse-grained velocity becomeswith mean 〈v̄_k〉 = 0, while its mean square is given byand thus grows linearly in k (where k ≈ t/τ), i.e. ∝t, leading to a 〈v̄²〉t̃² ∝ t̃³ behaviour in the MSD. This is similar to the t̃³-scaling of an active particle in shear flow.⁵² The velocity of an individual particle varies over time. The change in the coarse-grained velocity is given by the noise terms that are added for increasing k in eqn (29). These terms may be positive or negative depending on the realization of the thermal noise. The super-ballistic behaviour observed in the MSD is therefore not caused by the speed-up of an individual particle but rather by the accumulation of noise in the particle velocities leading to a broadening of the velocity distribution after the feedback potential is introduced. For appropriate choices of and t₀, the particle shows both ballistic and super-ballistic motion with the proportionality in the MSD changing as t̃–t̃²–t̃³–t̃, see Fig. 6. The t̃² regime corresponds to motion at constant velocity, given by the particle velocity at t₀ according to eqn (30), while the t̃³ regime is caused by the growing variance of the particle velocity due to the long memory in the system.

Fig. 5 Analytic result for the MSD in moderately strong repulsive harmonic feedback potentials of varying strength close to the divergence point and starting with the particle in the potential centre (i.e. t₀ = 0 fixed). The repulsive feedback potential, pushing the particle away from its past position, leads to an effective propulsion of the particle. Due to the growing (variance of the) particle velocity over time, the particle shows super-ballistic ∝t̃³ at intermediate times. The long-time diffusion coefficient remains the same as for the equilibrated system. The magenta-colored line corresponds to the highest order term of the MSD given in eqn (24), i.e. in the limit . [For reference: ln(2) ≈ 0.6931472.]

Fig. 6 Analytic result for the MSD in repulsive harmonic feedback potentials of varying strength starting from a not fully equilibrated feedback system (average taken over t₀ ∈ [100τ, 101τ) after first introduction of the feedback potential). The initial particle velocity (given by eqn (30) after 100τ) leads to a ∝t̃² regime. The cyan-colored line corresponds to a particle that travels at this constant velocity. The additional speed-up with time leads to a ∝t̃³ regime after multiple potential renewals. The magenta-colored line is the same as in Fig. 5. [For reference: ln(2) ≈ 0.6931472.]

In our calculations, we have derived a full analytic solution for the trajectory a single particle in a harmonic feedback potential. We have shown that the feedback potential can lead to both reduced and increased long-time diffusion as well as particle propulsion at intermediate times. Furthermore, care must be taken when defining the mean squared displacement (MSD) as the system is not completely symmetric under shifts in time and strong memory effects can lead to the retaining of the initial condition for long times. In the following, we drop the tilde from t̃ in the MSD, implying that the system is equilibrated long enough for the initial condition to be lost, and that the average over t₀ is executed properly such that only the relative time offset is relevant.

3.2 Experimental realization and simulation results

In the experiments, two types of feedback potentials are realized: first, an attractive harmonic potential is implemented. Second, a repulsive harmonic potential should be considered. However, as the colloidal particles are attracted to regions of high light intensity, construction of a repulsive potential is difficult. While in principle, it is possible to invert the intensity profile and illuminate the whole sample while decreasing the light intensity around the particle, this distributes the light intensity over a large area, requiring strong laser powers and potentially causing a significant heat-up in the sample. Therefore, in our work for the second potential, we chose to program an attractive ring instead. This potential pulls the particle towards a bright ring centered around the past particle position, thus leading to an effective repulsion.

3.2.1 Attractive harmonic potential. An attractive harmonic potential with stiffness κ > 0 is imposed onto the particle centered around its previous position. While the actual light beam has a Gaussian shape, the particle stays close to the potential centre such that the potential can be approximated as harmonic (Fig. 7A). Due to the time delay τ of the feedback loop, the harmonic potential centre is located at the previous (rather than the actual) particle position and thus the particle tends to move back towards it (Fig. 7B). The experimental setup was realized in three different ways: first, the spatial light modulator (SLM) was used as a mirror only and the Gaussian shape of the beam was exploited to create the desired light pattern. The beam location was controlled using galvanometer-mounted mirrors (GMM) (Fig. 7C). Second, the SLM was used to create a harmonic potential of different stiffnesses κ and the beam location was again controlled using the GMM (Fig. 7D). Third, the SLM was used to create a harmonic potential with different κ and to steer the beam (not shown in this figure, but discussed for other examples, see Fig. 9 and 10). All three realizations yielded equivalent results and perfectly agree with our analytic calculations.

Fig. 7 (A) One dimensional cut through the applied attractive potential. (B) The harmonic potential is centered around the position that the particle occupied a time τ earlier. Trajectories within the respective potential are indicated by the same color. The potential is kept the same for a duration τ while the next potential is calculated. As a result, the starting point of the coloured trajectory becomes the centre of the next programmed potential and the particle is pulled back to this previous position. (C) and (D) MSD as a function of time t obtained while applying a static harmonic potential (colored dotted lines (experimental data), dashed gray (analytic results)) and a harmonic feedback potential (coloured solid lines (experimental data), black lines (analytic results)) which was created by (C) a Gaussian beam and (D) the SLM. [Parameters of the analytic curves: D₀ = 0.026R² s⁻¹, τ = 0.2 s, values of κ/γ given in legend.]

In particular, we measured the two-dimensional mean squared displacement (MSD) 〈(r(t + t₀) − r(t₀))²〉, as a function of time t (with the average taken over t₀) first for a static harmonic potential (no feedback) and then for the feedback case. The static harmonic potential case is used to calibrate the potential strength, i.e. determine κ. By fitting the MSD to the known result for a harmonic trap,⁵³ we obtain both the diffusion coefficient D₀ within the potential as well as the potential stiffness κ, see dotted and dashed lines in Fig. 7C and D. Next, the feedback loop is introduced, i.e. the same harmonic potential is applied with its centre at the position the particle was at a time τ = 0.2 s earlier. As predicted by our analytic calculations, at short times, t < 0.1 s, diffusive motion is observed; at intermediate times, 0.1 s < t < 3 s, the delayed feedback results in oscillations; at long times, t > 3 s, diffusion is re-established with a significantly reduced diffusion coefficient. The long time diffusive behaviour then corresponds to the diffusive motion of the position of the potential centre. The experiments perfectly agree with our theoretical predictions (Fig. 7C and D).

3.2.2 Attractive annular (= effectively repulsive) potential. Next, we applied an annular ring potential, see Fig. 8A and B. Due to the form of the potential, the particle is pulled towards the ring surrounding its past position and is thus effectively pushed away from its past position. We use a ring potential of radius 1.4R (where R is the particle radius). In this case, the circular minimum typically is not reached by the particle within the delay time τ before the potential is relocated (see Fig. 8A and B). The light pattern is created using the SLM.

Fig. 8 (A) One dimensional cut through the applied annular potential. (B) The potential is centered around the position that the particle occupied a time τ earlier. Trajectories within the potential are indicated by the same color. The potential is kept the same for a duration τ while the next potential is calculated. The starting point of the colored trajectory then becomes the centre of the next programmed potential. The repeated relocation of the potential continuously pushes the particle away from the centre. Additionally, due to the time delay τ, the particle is typically off-centre in the refreshed potential. During two successive steps, the particle (as well as the potential) hence tends to advance in the same direction. (C) and (D) Mean squared displacement (MSD) and exponent μ as a function of time t in the experiment (C) and in the simulations (D). The programmed feedback potential leads to a ∝t² behaviour at intermediate times (μ ≈ 2) while the particle dynamics is diffusive at both short and long times (μ ≈ 1). The experimental results show the transition from ballistic to diffusive motion. (The initial diffusive regime is not covered.) The simulations reveal additional oscillatory behaviour with period of about τ resulting from the discretized potential update. [Experimental setup in (C): potential created with the SLM, potential strength is slightly larger for the red curve; simulation parameters for (D): D₀ = 0.026R² s⁻¹, τ = 0.2 s, b² = 0.1R², r_ring = 1.4R, potential amplitude given in figure legend, MSD averaged over 20 000 s.]

Three regimes of the MSD are observed (Fig. 8C): at short times, the particle freely diffuses with diffusion coefficient D₀. These times are typically shorter than the delay time τ = 0.2 s and hence are not recorded in the experiments (Fig. 8C) but covered in the simulations (Fig. 8D). At intermediate times, 0.03 s < t < 1 s, the motion is superdiffusive, almost ballistic. This is also indicated by the dynamical exponent μ = d ln〈(r(t₀ + t) − r(t₀))²〉/d ln t, which reaches a value of almost 2 in the experiments. At these times, the repeated relocation of the potential continuously pushes the particle away from the potential centre causing an effective propulsion. As the particle is typically off-centre in the refreshed potential, it tends to proceed in the same direction during successive steps. This propulsion direction is retained for long times (≈10 s). Since the particle also undergoes spatial diffusion, the propulsion direction slowly changes, such that at very long times diffusion is re-established, with a diffusion coefficient higher than the short-time one. The long-time diffusive motion can be understood as the diffusive motion of the potential centre. Again, the experimental findings qualitatively agree with simulations, in particular the transition from ballistic ∝t², to diffusive, ∝t motion (Fig. 8C and D). The simulations show additional oscillations of period τ in μ that are caused by the discretized potential updates: The effective propulsion is caused by the feedback force whose strength varies with the distance of the particle from its past position. The force thus varies periodically with period τ.

Note that the super-ballistic ∝t³ behavior observed in the simulations does not appear for this potential form as the former necessitates repulsion over a large spatial region.

4 Feedback-driven dynamics of groups of colloidal particles

In addition to individual particles, we also studied an ensemble of particles to explore their cooperative behaviour. Here, we program a potential landscape made up of annular rings centered around all previous particle positions.

For multiple particles and no feedback potential in a dilute system, the particles are found in a disordered fluid state. By imposing an attractive interaction potential, however, ordered structures can be induced and the formation of “living” chains and small clusters (Fig. 9) was observed (see ESI,† Movies S1 and S2). Through collisions, these chains grow and only very rarely break. Similar to the individual particles, the chains undergo superdiffusion. In very dilute samples, dimers exist for long times before forming trimers. They can be observed to not only perform translational motion but also rotate and slightly change their separation (Fig. 10 and Movie S3, ESI†). Similar chains and “rotators” were observed in computer simulations of the system (see ESI,† Movies S4 and S5). It is remarkable that the directed motion results from an instability-like situation: using an initial fluctuation the system is set into directed motion, while the subsequent iterations amplify the stochastic feedback systematically.

Fig. 9 Cooperative motion induced through programming. (a) Particles with an attractive annular potential with radius 2R and delay time τ = 0.2 s assemble in ‘living’ chain-like structures in the experiment. (b) Similar chain-like structures are obtained in the simulation. [Simulation parameters: N = 13, D₀ = 0.026R² s⁻¹, τ = 0.2 s, V₀/(γR²) = −0.35 s⁻¹, b² = 0.02R², r_ring = 2R, ρ = 0.02/R².] (Movies can be found in the ESI† as Movies S1 and S2 (experiment, same parameters (but different realization) for Movie S1, stronger potential for Movie S2) and S4 (simulation, same data set as shown above), respectively.)

Fig. 10 (top) Translation and rotation of dimers assembled through the feedback potential. Attractive annular potential with radius 4R and delay time τ = 0.2 s. The micrographs are taken in intervals of 0.2 s. (bottom) Simulation snapshots for an attractive annular potential. [Simulation parameters: N = 13, D₀ = 0.026R² s⁻¹, τ = 0.2 s, V₀/(γR²) = −20.0 s⁻¹, b² = 0.02R², r_ring = 2.0R, ρ = 0.02/R².] (Movies can be found in the ESI† as Movies S3 (experiment, same parameters but different realization) and S5 (simulation, same data set as shown above), respectively.)

Further, we programmed an attractive annular interaction potential with radius 2R and width R/7 and found crystallization of the particles into a triangular lattice (Fig. 11A). The triangular arrangement optimises the particle separation (to the radius of the interaction potential) and maximizes the number of neighbours with this separation. During the assembly process, the particles are observed to rotate around each other (similar to the “rotators” in Fig. 10) thus keeping their ideal separation while more neighbours approach (similar to the chain formation in Fig. 9). This separation, i.e. the size of the crystal unit cell, can be tuned through the radius of the annular interaction potential (Fig. 11A–C).

Fig. 11 Particle clusters formed for different parameters of the feedback potential. (A)–(E) Experimental results. (top) Schematic representation of the respective particle–particle interaction potential, which determines the total potential (middle) imposed on the particles (schematically represented in red), and the resulting particle arrangement (bottom). (F)–(H) Corresponding simulation snapshots. “Single ring” feedback potentials (A)–(C) lead to triangular lattices of corresponding spacing. “Double ring” feedback potentials (with ratio of radii ) (D) and (E) can be used to stabilize a square lattice. [Simulation parameters: D₀ = 0.026R² s⁻¹, τ = 0.2 s, (F)–(H) V₀/(γR²) = −10.0 s⁻¹, b² = 0.02R², (I) and (J) V₀/(γR²) = −5.0 s⁻¹, b² = 0.25R², N = 100, ρ = 0.04/R², (F) N = 100, ρ = 0.04/R², r_ring = 2R, (G) N = 50, ρ = 0.02/R², r_ring = 4R, (H) N = 40, ρ = 0.016/R², r_ring = 6R, (I) r_ring = 2R, , (J) r_ring = 2.2R, . All simulations were started from randomly placed particles except for (H) which was started from a triangular lattice.]

An interaction potential with two concentric attractive annuli of ratio (e.g. radii 2R (as before) and ) and identical width R/7 (as before) leads to a square arrangement of the particles (Fig. 11D). The larger distance along the diagonal of the unit cell is favoured by the second annulus which overcompensates the effect of the reduced number of neighbours and hence stabilises the square lattice. Again, an increase in the radii of the annuli leads to larger unit cells (Fig. 11E). The crystal formation was also found for computer simulations of the respective systems (Fig. 11F–J). Note that these crystal structures are not the result of a fixed template that favors certain spatial positions, but is a result of the self-organization of the particles according to the programmed interaction potential. The formed crystals are therefore invariant under translation and rotation and hence retain all possible zero energy modes.

5 Conclusions

We have introduced cybloids, i.e. cybernetic colloids, whose interaction potential is computer-controlled through a programmed feedback potential. Here, we have considered a time-delayed feedback potential that takes the past particle positions as an input parameter. For an experimental setup there will always be a slight delay between the input parameter and the potential output, due to the time needed to determine the parameter, calculate the potential and construct the corresponding light field. We have then systematically and quantitatively investigated the characteristics of these cybloids, including their arrangement, dynamics and collective behaviour (Fig. 7–11).

We have shown that the dynamics of the cybloids can be tuned by the programmed potential and delayed feedback. Through the programming, a cybloid particle is made to diffuse on short and long times with the diffusion coefficients given by the particle's and the potential's diffusion coefficients, respectively, while for intermediate times transient sub- or superdiffusion appears. This demonstrates that a broad range of dynamics is accessible and can be comprehensively investigated using this approach. Further, for groups of particle, we have shown the formation of “living” chains and small clusters as well as crystal lattices. These examples show that programming can also be exploited to induce cooperative behaviour, similar to the recently reported implementation of specific cooperation rules.⁵⁴

The dynamics of a single cybloid particle shown in Fig. 7 and 8 is a feature of the time delay. A particle is attracted to or repelled from its past position. As the constructed feedback potentials are symmetric around their center, without time delay the particle would not feel any force. Contrarily, for the formation of static structures such as the crystals shown in Fig. 11, we expect that the results do not differ significantly from particles with the corresponding instantaneous particle–particle interaction potential.

With the programmed potential, crystal formation can be induced and the crystal structure and the unit cell size controlled. Here, this is achieved through an externally imposed interaction potential between the particles, not a stationary particle-independent potential, i.e. fixed template, for example, an array of tweezers. Programmed particles hence provide the opportunity to systematically and quantitatively investigate the relation between particle–particle interactions and particle arrangements, e.g. different crystal symmetries.

In this work, we have chosen harmonic and Gaussian ring potentials to illustrate the potential of cybloids. We used the harmonic potential as a well-studied standard case for which the single particle dynamics is solvable analytically. As shown in Fig. 7C, experimentally, this case can be realized even without the use of an SLM. Likewise, the used ring-shaped potential can in principle also be realized with an axicon lens.

Beyond the examples we discussed here, programming can be exploited to impose pair and multi-body interactions, to spur drifts, bias and swarming, to induce memory effects or to couple two otherwise unrelated samples and to “clone” samples. Many further dependencies, situations and couplings can be imagined. In any case, the properties of each particle can be prescribed and hence the sample is controlled on all relevant length scales, without the need to change or directly interfere with the sample nor the need that such a sample can actually be prepared. This ideally complements the already very popular possibility to observe and follow each particle^21,55–59 and hence characterize the system on all relevant length scales. Therefore, cybloids open a huge range of new opportunities which can be exploited to create previously inaccessible model systems, explore the rational design of materials or mimic natural and industrially relevant situations.

Author contributions

DS performed the experiments, ST did the theoretical calculations and performed the simulations. SUE and HL supervised the project. DS, ST, SUE and HL discussed the results and contributed to writing the paper.

Data availability

Simulations were performed with a forward Euler scheme based on eqn (1) and (7) in the manuscript. All necessary simulation parameters are provided in Section 2.2 of the manuscript and in the respective figure captions (Fig. 7–11). Fig. 3–6 are plots of the analytic result which is given in full in eqn (18) of the manuscript. Experimental data used to generate Fig. 7 and 8 are available as ESI.†

Conflicts of interest

The authors declare no competing interests.

Acknowledgements

We mourn the loss of our colleague Stefan Egelhaaf, who left us far too early. Many of the ideas in this paper originated from him. We are grateful for his creativity, his thoroughness, his persistence in pursuing questions and questioning data until they were consistent and understood, and for the many joint discussions that ultimately led to this publication. This work was financially supported by a joint project of SUE and HL of the German Research Foundation (DFG) with the project numbers EG 269/6 and LO 418/19. DS thanks Dennis Rohrschneider for his help with the experiments.

References

1. J. Hafner, From Hamiltonians to Phase Diagrams: The Electronic and Statistical-Mechanical Theory of sp-Bonded Metals and Alloys, Springer textbook, 1987, vol. 70.
2. H. Löwen, Phys. Rep., 1994, 237, 249.
3. A. Yethiraj and A. van Blaaderen, Nature, 2003, 421, 513.
4. J. Dzubiella, J. Chakrabarti and H. Löwen, J. Chem. Phys., 2009, 131, 044513.
5. F. Schmidt, A. Callegari, A. Daddi-Moussa-Ider, B. Munkhbat, R. Verre, T. Shegai, M. Käll, H. Löwen, A. Gambassi and G. Volpe, Nat. Phys., 2022, 19, 271.
6. M. Oh and C. A. Mirkin, Nature, 2005, 438, 651.
7. M. Rey, M. A. Fernandez-Rodriguez, M. Karg, L. Isa and N. Vogel, Acc. Chem. Res., 2020, 53, 414.
8. M. J. Bergman, N. Gnan, M. Obiols-Rabasa, J.-M. Meijer, L. Rovigatti, E. Zaccarelli and P. Schurtenberger, Nat. Commun., 2019, 9, 5039.
9. J. Harrer, M. Rey, S. Ciarella, H. Löwen, L. M. C. Janssen and N. Vogel, Langmuir, 2019, 35, 10512.
10. K. Zahn, R. Lenke and G. Maret, Phys. Rev. Lett., 1999, 82, 2721.
11. U. Gasser, C. Eisenmann, G. Maret and P. Keim, ChemPhysChem, 2010, 11, 963.
12. P. Dillmann, G. Maret and P. Keim, Eur. Phys. J.-Spec. Top., 2013, 222, 2941.
13. L. Assoud, F. Ebert, P. Keim, R. Messina, G. Maret and H. Löwen, Phys. Rev. Lett., 2009, 102, 238301.
14. S. Deutschländer, T. Horn, H. Löwen, G. Maret and P. Keim, Phys. Rev. Lett., 2013, 111, 098301.
15. S. Huang, G. Pessot, P. Cremer, R. Weeber, C. Holm, J. Nowak, S. Odenbach, A. M. Menzel and G. K. Auernhammer, Soft Matter, 2016, 12, 228.
16. P. Pusey, E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon and M. E. Cates, Philos. Trans. R. Soc., A, 2009, 367, 4993.
17. A. Ivlev, H. Löwen, G. Morfill and C. P. Royall, Complex plasmas and colloidal dispersions: particle-resolved studies of classical liquids and solids, World Scientific, Series in Soft Condensed Matter, 2012, vol. 5, pp. 1–336.
18. G. L. Hunter and E. R. Weeks, Rep. Prog. Phys., 2012, 75, 066501.
19. V. J. Anderson and H. N. W. Lekkerkerker, Nature, 2002, 416, 811.
20. M. Dijkstra, R. van Roij, R. Roth and A. Fortini, Phys. Rev. E, 2006, 73, 041404.
21. E. Allahyarov, K. Sandomirski, S. U. Egelhaaf and H. Löwen, Nat. Commun., 2015, 6, 7110.
22. B. Li, D. Zhou and Y. Han, Nat. Rev. Mater., 2016, 1, 15011.
23. M. Chaudhuri, E. Allahyarov, H. Löwen, S. U. Egelhaaf and D. A. Weitz, Phys. Rev. Lett., 2017, 119, 128001.
24. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
25. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
26. F. Höfling and T. Franosch, Rep. Prog. Phys., 2013, 76, 046602.
27. V. Lobaskin and M. Romenskyy, Phys. Rev. E, 2013, 87, 052135.
28. P. Grančič and F. Štěpánek, J. Phys. Chem. B, 2013, 117, 8031.
29. G. Ariel, A. Rabani, S. Benisty, J. D. Partridge, R. M. Harshey and A. Beer, Nat. Commun., 2015, 6, 8396.
30. E. F. W. Heffern, H. Huelskamp, S. Bahar and R. F. Inglis, Proc. R. Soc. B, 2021, 288, 20211111.
31. N. Dorsaz, G. M. Thurston, A. Stradner, P. Schurtenberger and G. Foffi, Soft Matter, 2011, 7, 1763.
32. S. C. Glotzer and M. J. Solomon, Nat. Mater., 2007, 6, 557.
33. V. N. Manoharan, Science, 2015, 349, 1253751.
34. S. Tarama, S. U. Egelhaaf and H. Löwen, Phys. Rev. E, 2019, 100, 022609.
35. M. Ángel Fernández-Rodríguez, F. Grillo, L. Alvarez, M. Rathlef, I. Buttinoni, G. Volpe and L. Isa, Nat. Commun., 2020, 11, 4223.
36. J.-J. Liao and F.-J. Lin, Acta Phys. Sin., 2020, 69, 220501.
37. L. Alvarez, M. A. Fernandez-Rodriguez, A. Alegría, S. Arrese-Igor, K. Zhao, M. Kröger and L. Isa, Nat. Commun., 2021, 12, 4762.
38. H. Massana-Cid, C. Maggi, G. Frangipane and R. Di Leonardo, Nat. Commun., 2022, 13, 2740.
39. R. A. Kopp and S. H. L. Klapp, Phys. Rev. E, 2023, 107, 024611.
40. R. A. Kopp and S. H. L. Klapp, EPL, 2023, 143, 17002.
41. F. Pakpour and T. Vicsek, Phys. A, 2024, 634, 129453.
42. J. C. Crocker and D. G. Grier, J. Colloid Interface Sci., 1996, 179, 298.
43. R. D. L. Hanes, M. C. Jenkins and S. U. Egelhaaf, Rev. Sci. Instrum., 2009, 80, 083703.
44. A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm and S. Chu, Opt. Lett., 1986, 11, 288.
45. R. F. Capellmann, A. Khisameeva, F. Platten and S. U. Egelhaaf, J. Chem. Phys., 2018, 148, 114903.
46. Y. Yang, A. Forbes and L. Cao, Opto-Electron. Sci., 2023, 2, 230026.
47. M. C. R. Bell-Davies, A. Curran, Y. Liu and R. P. A. Dullens, Phys. Rev. E, 2023, 107, 064601.
48. M. C. Jenkins and S. U. Egelhaaf, J. Phys.: Condens. Matter, 2008, 20, 404220.
49. R. Gerchberg and R. W. Gerchberg, Optik, 1972, 35, 237.
50. U. Küchler and B. Mensch, Stochastics and Stochastics Reports, 1992, vol. 40, p. 23.
51. L. Giuggioli, T. J. McKetterick, V. M. Kenkre and M. Chase, J. Phys. A, 2016, 49, 384002.
52. B. ten Hagen, R. Wittkowski and H. Löwen, Phys. Rev. E, 2011, 84, 031105.
53. B. Lukić, S. Jeney, I. C. V. Sviben, A. J. Kulik, E.-L. Florin and L. Forró, Phys. Rev. E, 2007, 76, 011112.
54. F. A. Lavergne, H. Wendehenne, T. Bäuerle and C. Bechinger, Science, 2019, 364, 70.
55. V. Prasad, D. Semwogerere and E. R. Weeks, J. Phys.: Condens. Matter, 2007, 19, 113102.
56. C. P. Royall, A. A. Louis and H. Tanaka, J. Chem. Phys., 2007, 127, 044507.
57. S.-H. Lee, Y. Roichman, G.-R. Yi, S.-H. Kim, S.-M. Yang, A. van Blaaderen, P. van Oostrum and D. G. Grier, Opt. Express, 2007, 15, 18275.
58. M. Jenkins and S. Egelhaaf, Adv. Colloid Interface Sci., 2008, 136, 65.
59. M. Leocmach and H. Tanaka, Soft Matter, 2013, 9, 1447.

---

Footnotes

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

‡ These authors contributed equally to this work.

§ Present address: Laboratory of Biological Computation, College of Life Sciences, Department of Bioinformatics, Ritsumeikan University, 1 Chome-1-1 Nojihigashi, Kusatsu, Shiga 525-0058, Japan.

---