Traveling band formation in feedback-driven colloids

TL;DR

This paper demonstrates how time-delayed feedback in colloidal suspensions induces self-organized traveling bands, with properties tunable by feedback parameters, supported by simulations and theoretical analysis.

Contribution

It introduces a novel feedback mechanism causing traveling band formation in colloids, combining simulation and theory to analyze the phenomenon.

Findings

Traveling bands form due to feedback in colloids.

Band width and speed are tunable by delay and potential range.

Traveling bands persist in finite and confined systems.

Abstract

Using simulation and theory we study the dynamics of a colloidal suspension in two dimensions subject to a time-delayed repulsive feedback that depends on the positions of the colloidal particles. The colloidal particles experience an additional potential that is a superposition of repulsive potential energies centered around the positions of all the particles a delay time ago. Here we show that such a feedback leads to self-organization of the particles into traveling bands. The width of the bands and their propagation speed can be tuned by the delay time and the range of the imposed repulsive potential. The emerging traveling band behavior is observed in Brownian dynamics computer simulations as well as microscopic dynamic density functional theory (DDFT). Traveling band formation also persists in systems of finite size leading, in the case of circularly confined systems, to rotating traveling waves.

Figures

Equations

$$\gamma \frac{\mathrm{d}{\mathbf{r}}_i}{\mathrm{d}t}={\mathbf{f}}_i(t)$$

$$+\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{\mathbf{F}}_{\mathrm{Yuk}}\left({\mathbf{r}}_i(t)-{\mathbf{r}}_j(t)\right).$$

$$\mathbf{F}(\mathbf{r})=-{\nabla }_{\mathbf{r}}{V}_{\mathrm{fb}}\left(r\right)$$

$$V_{\mathrm{fb}}(r)=A\exp (-\frac{r^2}{2b^2}),$$

$${\mathbf{F}}_{\mathrm{Yuk}}(\mathbf{r})=-{\nabla }_{\mathbf{r}}\Phi (r)$$

$$\Phi (r)=\frac{V_0}{r}\exp (-\kappa r)$$

$$v_{\mathrm{s}}=\frac{{\lambda }_{\mathrm{b}}}{2\tau }.$$

$$\displaystyle\mathbf{F}_{i}(t)=\Bigg{[}$$

$$\displaystyle+\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}\mathbf{F}_{\mathrm{Yuk}}\left(\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)\right)\Bigg{]}\,.$$

$${\mathbf{v}}_i(t)=\frac{{\mathbf{F}}_i(t)}{\gamma }$$

$$\mathbf{v}$$

$$⟨B(t^{\prime })⟩=\frac{1}{T}{\int }_{t_0}^{t_0+T}\mathrm{d}{t}^{\prime }B(t^{\prime })$$

$$\rho_{\mathrm{b}}\left(s\right)=\Big{\langle}\sum_{i=1}^{N}\delta\left(s-\mathbf{r}_{i}(t^{\prime})\cdot\hat{\mathbf{e}}_{\mathrm{v}}\right)\Big{\rangle}\,.$$

$$S(\mathbf{k})=\Bigg{\langle}\frac{1}{N}\sum_{i,j=1}^{N}\mathrm{e}^{-i\mathbf{k}\left(\mathbf{r}_{i}(t^{\prime})-\mathbf{r}_{j}(t^{\prime})\right)}\Bigg{\rangle}\,.$$

$$\Delta(t)=\Bigg{\langle}\frac{1}{N}\sum_{i=1}^{N}\left(\mathbf{r}_{i}\left(t^{\prime}+t\right)-\mathbf{r}_{i}\left(t^{\prime}\right)\right)^{2}\Bigg{\rangle}\,.$$

$$\frac{\partial }{\partial t}X(t)=h\left(X(t),X(t-\tau )\right)+\sqrt{2D}\Gamma (t),$$

$$\frac{\partial}{\partial t}w(x,t|x_{\tau},t_{\tau})\Big{|}_{t_{\tau}=t-\tau}=\hat{L}(x,\nabla_{x},x_{\tau})w(x,t|x_{\tau},t_{\tau})$$

$$\hat{L}(x,{\nabla }_x,x_{\tau })=$$

$$\hat{L}(\mathbf{r},\nabla ,{\mathbf{r}}_{\tau })=$$

$$\displaystyle\frac{\partial w(\mathbf{r},t|\mathbf{r}_{\tau},t_{\tau})}{\partial t}\Bigg{|}_{t_{\tau}=t-\tau}$$

$$=\frac{1}{\gamma }\nabla \left[k_{B}T\nabla -\mathbf{F}(\mathbf{r},{\mathbf{r}}_{\tau })\right]w(\mathbf{r},t|{\mathbf{r}}_{\tau },t_{\tau }).$$

$$\displaystyle\frac{\partial w(\mathbf{r},t;\mathbf{r}_{\tau},t_{\tau})}{\partial t}\Bigg{|}_{t_{\tau}=t-\tau}$$

$$=\frac{1}{\gamma }\nabla \left[k_{B}T\nabla -\mathbf{F}(\mathbf{r},{\mathbf{r}}_{\tau })\right]w(\mathbf{r},t;{\mathbf{r}}_{\tau },t_{\tau }).$$

$$\frac{\partial w(\mathbf{r},t)}{\partial t}=$$

$$-\frac{1}{\gamma }\int \mathrm{d}{\mathbf{r}}_{\tau }\nabla \mathbf{F}(\mathbf{r},{\mathbf{r}}_{\tau })w(\mathbf{r},t;{\mathbf{r}}_{\tau },t_{\tau }).$$

$$\displaystyle\frac{\partial w(\mathbf{r}^{N},t|\mathbf{r}_{\tau}^{N},t_{\tau})}{\partial t}\Bigg{|}_{t_{\tau}=t-\tau}$$

$$=\frac{1}{\gamma }\sum _{i=1}^N{\nabla }_i\left[k_{B}T{\nabla }_i-\left(-{\nabla }_i{U}_{\mathrm{tot}}({\mathbf{r}}^N,{\mathbf{r}}_{\tau }^N)\right)\right]w({\mathbf{r}}^N,t|{\mathbf{r}}_{\tau }^N,t_{\tau }).$$

$$U_{\mathrm{tot}}\left({\mathbf{r}}^N,{\mathbf{r}}_{\tau }^N\right)=$$

$$+\sum _{i,j=1}^N{V}_{\mathrm{fb}}(|{\mathbf{r}}_i-{\mathbf{r}}_{{\tau }_j}|).$$

$$\displaystyle\frac{\partial w(\mathbf{r}^{N},t;\mathbf{r}_{\tau}^{N},t_{\tau})}{\partial t}\Bigg{|}_{t_{\tau}=t-\tau}$$

$$=\frac{1}{\gamma }\sum _{i=1}^N{\nabla }_i\left[k_{B}T{\nabla }_i+\left({\nabla }_i{U}_{\mathrm{tot}}({\mathbf{r}}^N,{\mathbf{r}}_{\tau }^N)\right)\right]w({\mathbf{r}}^N,t;{\mathbf{r}}_{\tau }^N,t_{\tau }).$$

$$w({\mathbf{r}}^N,t)=\int \mathrm{d}{\mathbf{r}}_{{\tau }_1}\cdots \int \mathrm{d}{\mathbf{r}}_{{\tau }_N}w({\mathbf{r}}^N,t;{\mathbf{r}}_{\tau }^N,t_{\tau }).$$

$$\int \mathrm{d}{\mathbf{r}}_{{\tau }_1}\cdots \int \mathrm{d}{\mathbf{r}}_{{\tau }_N}{U}_{\mathrm{tot}}({\mathbf{r}}^N,{\mathbf{r}}_{\tau }^N)w({\mathbf{r}}^N,t;{\mathbf{r}}_{\tau }^N,t_{\tau })$$

$$=\int \mathrm{d}{\mathbf{r}}_{{\tau }_1}\cdots \int \mathrm{d}{\mathbf{r}}_{{\tau }_N}{U}_{\mathrm{tot}}({\mathbf{r}}^N,{\mathbf{r}}_{\tau }^N)w({\mathbf{r}}^N,t|{\mathbf{r}}_{\tau }^N,t_{\tau })w({\mathbf{r}}_{\tau }^N,t_{\tau }),$$

$$\frac{\partial w({\mathbf{r}}^N,t)}{\partial t}=$$

Full text

††preprint: APS/123-QED

I Introduction

Nonequilibrium systems subject to a feedback potential have been studied extensively in recent times Granger et al. (2016); Loos et al. (2014); Loos and Klapp (2017, 2019); von Lospichl and Klapp (2018); Lopez et al. (2008); Gernert and Klapp (2015); Popli et al. (2018); Yang and Bevan (2018); Lichtner and Klapp (2010). Due to the feedback, used e.g. to stabilize dynamics Lichtner and Klapp (2010); Pyragas (1992); Hövel and Schöll (2005); Cohen (2005); Jun and Bechhoefer (2012) or structure Popli et al. (2018); Khadka et al. (2018), the system dynamics becomes history-dependent. The feedback can be realized through external programming of a laser trap Popli et al. (2018); Blickle and Bechinger (2011); Hanes et al. (2009); Evers et al. (2013); Bewerunge and Egelhaaf (2016); Bäuerle et al. (2018); Jones et al. (2015); Nishizawa et al. (2017) or, more naturally, may arise in autochemotactic particles, i.e., if the particles themselves are part of the production mechanism of the chemical substance they react to. In particular, examples of the latter include biological systems such as bacteria Adler (1966); Keller and Segel (1971) and army ants Couzin and Franks (2003), as well as synthetic microswimmers such as active colloidal particles Bechinger et al. (2016); Liebchen et al. (2015, 2016, 2017); Saha et al. (2014); Sengupta et al. (2009, 2011) or self-propelling droplets Jin et al. (2017, 2018).

In the context of many-particle systems, the topic of pattern formation Cross and Greenside (2009); Cross and Hohenberg (1993); Gonnella et al. (1997); Zakine et al. (2018); Chacko et al. (2015); Krekhov et al. (2013); Lushi et al. (2018) is of central interest. In particular, the Ginzburg-Landau Aranson and Kramer (2002); Puzyrev et al. (2014); Ciszak et al. (2015); Emmerich et al. (2012) and Swift-Hohenberg equations Swift and Hohenberg (1977) are widely used to study pattern formation. Most of these studies present a coarse-grained treatment using effective continuum theories but do not resolve the individual particles. One pattern which is commonly observed in many different systems is traveling waves or moving bands of particles. Examples include actin-waves formed in the biological actin-myosin systems Gerisch et al. (2004); Le Goff et al. (2016); Schaller et al. (2010), metachronal waves in cilia arrays Elgeti and Gompper (2013), the patterning in systems of active agents under various settings Vicsek et al. (1995); Chaté et al. (2008); Menzel and Löwen (2013); Krekhov et al. (2013); Ophaus et al. (2018); Reinken et al. (2019); Schnyder et al. (2017), the formation of bands in passive colloidal suspensions driven by ac Vissers et al. (2011); Wysocki and Löwen (2009) or dc Dzubiella et al. (2002); Reichhardt and Olson Reichhardt (2007); Ikeda et al. (2012); Reichhardt et al. (2018); Reichhardt and Reichhardt (2018); Dzubiella and Löwen (2002) fields, and phase separating mixtures Okuzono and Ohta (2001, 2003); Krekhov et al. (2013). Recent work on pattern forming systems also considers the effect of time-delayed feedback using continuum theories Gurevich (2013); Tabbert et al. (2017); Gurevich and Friedrich (2013); Ciszak et al. (2015).

In this paper, we present a study of feedback-driven colloidal particles as an example of a feedback system of discrete components considered on the fundamental particle level. In our model, the particles are subjected to a feedback potential driving them away from their previous positions. Using Brownian dynamics computer simulations and dynamical density functional theory Marconi and Tarazona (1999, 2000); Archer and Evans (2004); Español and Löwen (2009); Evans et al. (2016); Wittkowski and Löwen (2011), we show that this repulsive feedback leads to self-organization of the particles into a moving band structure reminiscent of a traveling wave. Remarkably, this ordering takes place despite the absence of any attractive interactions in the system, for which static band formation is known to occur Olson Reichhardt et al. (2010); Zhao et al. (2012). The width of the bands and their propagation speed can be tuned by the delay time and the range of the imposed repulsive delay potential. Finally, we demonstrate that traveling band formation also persists under strong confinement leading, in circularly confined systems, to globally rotating and spiraling bands Wioland et al. (2013).

Our model can be realized in experiments for colloidal suspensions. The suspensions can be exposed to a potential energy landscape using optical fields which are programmed via a feedback loop Bäuerle et al. (2018); Leyman et al. (2018); Lavergne et al. (2019). Typically the colloids are attracted towards the intensity maximum of the optical field. However, by inverting the intensity landscape a repulsion is achieved, in which case the particles are driven away from the dark regions.

The paper is organized as follows: In the following section we introduce the underlying Langevin equation including a delay term Ohira and Yamane (2000); Ohira and Hosaka (2005); Ohira and Milton (2009); Trimper and Zabrocki (2004); Le Berre et al. (1990); Guillouzic et al. (1999); Atay (2010); Loos and Klapp (2019) describing the dynamics of the system. We continue with presenting our simulation results in Sec. III. A prediction of the observed traveling wave formation is derived from dynamic density functional theory in Sec. IV. Subsequently, we consider confinement effects in Sec. V. Finally, we conclude with a summary of our main findings and an outlook to possible extensions of the system in Sec. VI.

II Model and Brownian dynamics computer simulations

The Brownian dynamics of $N$ colloidal particles in two spatial dimensions is described by their time-dependent positions $\mathbf{r}_{i}(t)$ ($i=1,\ldots,N$) and governed by the following Langevin equation

[TABLE]

which can be viewed as a force balance equation. The left-hand side of Eq. (1) contains the Stokes drag force with $\gamma$ denoting the friction coefficient. The Gaussian random force $\mathbf{f}_{i}(t)$ mimics the collision of the particle with solvent molecules. This stochastic force is characterized by its first two moments $\langle\mathbf{f}_{i}(t)\rangle=0$ and $\langle\mathbf{f}_{i}(t)\otimes\mathbf{f}_{j}(t^{\prime})\rangle=2D\gamma^{2}\mathds{1}\delta\left(t-t^{\prime}\right)\delta_{ij}$, where $D$ is the short-time diffusion coefficient of the particles, $\delta\left(t\right)$ is the Dirac delta function and $\delta_{ij}$ denotes the Kronecker delta. The important new ingredients in Eq. (1) are the feedback forces $\mathbf{F}\left(\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t-\tau)\right)$. These forces are evaluated at distances between the actual position $\mathbf{r}_{i}(t)$ of particle $i$ and the former positions $\mathbf{r}_{j}(t-\tau)$ of the other particles $j$ (where the special case $i=j$ is included). Here, $\tau$ is the time difference, which we refer to as the delay time of the feedback.

We derive $\mathbf{F}(\mathbf{r})$ from a potential $V_{\mathrm{fb}}(r)$ as

[TABLE]

and assume for simplicity a Gaussian form

[TABLE]

characterized by an energy amplitude $A$ and a range $b$. The Gaussian potential form is a good approximation for optical systems such as optical tweezers and occurs naturally for autochemotactic particles Sengupta et al. (2009). Here, we confine ourselves to the case of repulsive feedback potentials such that the energy amplitude $A>0$ is positive and the special case $A=0$ serves as an equilibrium reference case. For $A>0$, all particles are driven away from the past positions of all particles including their own.

Finally, the equations of motion include direct particle-particle interaction forces

[TABLE]

via a repulsive Yukawa pair potential

[TABLE]

involving an inverse range $\kappa$ and an amplitude $V_{0}$.

We perform Brownian dynamics simulations with a square simulation box of length $L$ and periodic boundary conditions with $N=6400$ particles. Some of the simulations were repeated in a rectangular box, in order to obtain stable traveling bands. This was necessary because at the onset of the formation of traveling bands, the band stability is highly dependent on the commensurability of the box size and the preferred wavelength. Possible wavelengths in the finite system are restricted to those being commensurate with the periodic boundaries, which thus requires the system length to be adjusted. The equation for the particle positions, Eq. (1), is integrated using an explicit Euler scheme with a finite time step of $\Delta t=10^{-4}\tau_{0}$, where $\tau_{0}=b^{2}/D$ denotes the Brownian time scale.

In the following, lengths are normalized to the feedback potential range $b$ and times to the Brownian time $\tau_{0}=b^{2}/D$. Energies are given in terms of the thermal energy $k_{B}T\equiv D\gamma$. In order to keep the set of parameters limited, we maintain $V_{0}=60\,bk_{B}T$, $\kappa=4.5/b$ and $\tau=0.25\,\tau_{0}$ constant in our units. Furthermore, we use a number density $\rho_{0}=N/L^{2}=1/b^{2}$. The units are dropped hereafter for ease of notation. We use the feedback amplitude $A$ as a control parameter and investigate the change in the system structure and dynamics as a function of it.

Our simulation protocol is as follows: First, the system is equilibrated without any feedback potential (corresponding to a two-dimensional pure Yukawa system Löwen (1992)), after which the positions are recorded for updating the feedback potential. We define $t=0$ as the time at which the feedback potential is first introduced into the system. Subsequently, the relaxation of the system is monitored for a long time, several hundreds of time units.

III Simulation Results

III.1 Band formation

Figure 1 shows a typical instance of system relaxation in the case of strong feedback potentials. The initial equilibrated homogeneous fluid state in Fig. 1(a) spontaneously separates into two regions which are either empty or crowded (i.e., exhibiting a high density of particles). This demixed state coarsens further as a function of time into a configuration of system-spanning straight bands at long times, resembling those observed in mixtures of particles subjected to bidisperse Magnus forces Reichhardt and Reichhardt (2019). The width of the formed bands depends on the specific parameters of the particle repulsion and the feedback potential as well as the delay time $\tau$. In particular, no band formation is observed for very small delay times, in which case the feedback potential can be understood as modified direct interactions between the particles, while larger delay times lead to an increase in the time scale on which band formation takes place. The orientation of the bands is tilted relative to the quadratic box with the angle depending on the initial configuration as well as the commensurability of the wavelength and the box size. The commensurability condition leads to a finite set of possible orientations for the bands while the initialization determines which of these is realized.

The emerging bands are observed to move collectively along the normal of their interfaces. The empty regions are found at the former particle positions, i.e., the positions where the potential is inserted, suggesting that particles try to effectively avoid regions where they have been a time $\tau$ before. In more detail, the occurrence of the separation into empty and crowded regions can be explained in the following way:

If in a disordered system as in Fig. 1(a) there are by chance more particles at a particular position at time $t$, the strong repulsive potential imposed at time $t+\tau$ leads to fewer particles at this position, which in turn leads to a small potential and more particles at time $t+2\tau$. The feedback potential thus leads to a self-ordering and the particle distribution effectively changes between one state and its negative image with period $2\tau$. The easiest way to achieve this, namely, the one with the fewest collisions between the particles, is a collective movement into one direction as given for the moving lamellar phase seen in Fig. 2.

Based on the previous consideration, a band moves over its full periodicity $\lambda_{\mathrm{b}}$ during twice the delay time $\tau$ such that a scaling expression for the magnitude of the expected band velocity is obtained as

[TABLE]

Remarkably, knowledge of the static property, namely, the wavelength and orientation of the band structure, thus provides an estimate of the dynamics of the system, i.e., the band velocity. Likewise, determining the velocity of the bands for a given feedback time yields an approximative value for the band periodicity $\lambda_{\mathrm{b}}$. Moreover, through determination of the velocity and periodicity, the feedback time can be estimated, which, in particular in biological systems, might not be easily accessible otherwise.

The prediction for the band velocity obtained from the scaling expression can be compared to the simulation results. For the latter case, we use the systematic force

[TABLE]

acting on particle $i$ at time $t$ to define an instantaneous drift velocity

[TABLE]

of this particle. From this expression, the mean global drift velocity $\mathbf{v}$ is then obtained by averaging over all particles as

[TABLE]

where

[TABLE]

denotes an average taken over time $t^{\prime}$ for an observable $B(t^{\prime})$. The time $t_{0}$ is bigger than a typical relaxation time of the system and $T$ is the width of the time window over which the average is performed. Here, we use $t_{0}=500$ and $T=500$. In the ordered band state the mean drift velocity equals the velocity of the bands. The comparison between the two velocities is shown in the next section as a function of the feedback potential amplitude $A$ [see Fig. 7].

III.2 Dependence on the feedback strength

The system shows different structuring depending on the amplitude of the applied feedback potential $A$. The patterns are shown in Fig. LABEL:fig_snapshot_t500 with the color code indicating the drift velocity directions of the individual particles. As a general result, we find that the average of these velocities, the mean drift velocity $\mathbf{v}$, is typically normal to the band direction.

With respect to potential strength $A$, we observe that, while for small $A\lesssim 1$ [Fig. LABEL:fig_snapshot_t500(a)], diffusion prevents any structure formation, higher potential amplitudes lead to patterning [Fig. LABEL:fig_snapshot_t500(b)] and, for even larger potential strength, to the formation of a band structure. The necessary potential amplitude for pattern formation can be estimated by considering at which point the feedback force $\mathbf{F}$ becomes comparable to the interparticle repulsion force $\mathbf{F}_{\mathrm{Yuk}}$. From the condition $\mathbf{F}^{2}=\mathbf{F}_{\mathrm{Yuk}}^{2}$, we find that $A\approx 6$, which is in reasonable agreement with the simulation results which indicate the start of band formation at $5\leq A^{*}\leq 6$. Close to this threshold, the band state can be constituted of two distinct band orientations [Fig. LABEL:fig_snapshot_t500(c)]. Further increasing the potential leads to a single band orientation with a considerable number of particles traveling in the opposite direction [Fig. LABEL:fig_snapshot_t500(d)]. At even higher potentials [Figs. LABEL:fig_snapshot_t500(e) and LABEL:fig_snapshot_t500(f)], the system forms stable bands with all particles moving as part of the bands into the same direction.

We have checked the system for hysteretic behavior via additional simulations. For a finite system, a small hysteresis effect is observed in forming the bands, which appears to be due to the finite system size and the prescribed periodic boundary conditions.

The density distribution within the bands $\rho_{\mathrm{b}}$ can be determined by changing to a co-moving frame. Taking an average along the band tangential, the density only depends on the position $s$ in the drift direction $\hat{\mathbf{e}}_{\mathrm{v}}$,

[TABLE]

One period of these profiles is shown in Fig. 4, revealing an increasing layering of particles for strong feedback amplitudes and hence an increasing internal order of the bands.

In the following, the formed patterns are explored in more detail via the structure factor, defined by

[TABLE]

Figure LABEL:fig_Skxky shows $S(\mathbf{k})$ for the cases presented in Fig. LABEL:fig_snapshot_t500 as a function of the components of the wave vector $\mathbf{k}=\left(k_{x},k_{y}\right)$. The outer black ring corresponds to the mean particle distance. For higher potential, its radius is shifted away from the equilibrium value $2\pi\rho_{0}^{1/2}$ towards higher wave numbers, indicating that the mean distance between particles in the band structure is considerably smaller than the one in the homogeneous system without the feedback potential.

Further, the inner black ring first appearing in Fig. LABEL:fig_Skxky(b) represents the ordering due to the feedback potential. With respect to the feedback potential strength $A$, different stages of ordering are observed at this wave number which we denote by $k^{*}$. First, for small potentials, directionally independent patterning at $k^{*}\approx 0.82$ is found [Fig. LABEL:fig_Skxky(b)], which for stronger amplitudes develops a directional dependence with two preferred band orientations [Fig. LABEL:fig_Skxky(c)]. While medium potentials [Figs. LABEL:fig_Skxky(c) and LABEL:fig_Skxky(d)] still show some remainder of the initial orientationally-independent ordering, indicated by the light gray ring at $k^{*}$, this feature disappears at large feedback potential strength $A$ [Figs. LABEL:fig_Skxky(e) and LABEL:fig_Skxky(f)]. For strong feedback $A$, only a single band orientation is found. For large potential amplitudes [Figs. LABEL:fig_Skxky(d)–(f)], the nonsinusoidal form of the density profile, visualized in Fig. 4, is reflected in higher harmonics which lead to additional peaks at multiples of $k^{*}$ in the structure factor. For easier comparison, the azimuthal average $S(k)$ is shown in Fig. 6 for different feedback strengths $A$, illustrating the shift in the mean particle distance as well as the growth of the new wave number $k^{*}$ and its higher harmonics.

Extracting the wavelength of the band pattern from the structure factor $S(k)$, allows for a comparison between the estimated value of the band velocity obtained from the scaling expression introduced in the previous section and the simulation results. The prediction for the band velocity according to the scaling expression $v_{\mathrm{s}}$, given by Eq. (6), and the simulation results for the magnitude $v$ of the global drift velocity $\mathbf{v}$, defined by Eq. (9), are shown in Fig. 7. Above a threshold value $5\leq A^{*}\leq 6$, $v$ increases sharply. For high feedback amplitudes $A$ it saturates to a value that depends on the specific choice of the remaining potential parameters as well as the delay time $\tau$. The agreement between the theoretical prediction and the simulation results is acceptable but not exact. The reason for the discrepancy is that the bands filled with particles and the empty spaces in between are not exactly equal in width, a consequence of the feedback force pushing the particles from behind. The difference in width is also visible in Fig. 2. The length traveled by the bands is thus overestimated when using half the wavelength of the band structure instead of the width of the bands, leading to a slightly higher predicted velocity for the present parameters.

For the individual particles, the directed drift motion becomes visible in the mean-squared displacement (MSD), which we define by

[TABLE]

The MSD changes qualitatively with the onset of the traveling wave instability: Increasing $A$ over the threshold value for band formation changes the long time behavior from diffusive ($\propto t$) to a directed drift motion ($\propto t^{2}$), similar to what is found for active (self-propelled) Brownian particles Howse et al. (2007); Bechinger et al. (2016) (see Fig. 8). Thus, the feedback potential effectively provides a source of self-propulsion.

IV Dynamical density functional theory

The previously observed values for the potential strength $A^{*}$ and wave number $k^{*}$ characterizing the onset of the pattern formation instability are now derived from microscopic dynamic density functional theory (DDFT). The theory requires the particle interactions as the only input.

DDFT is based on the Smoluchowski equation Dhont (1996) corresponding to the Langevin equation (1). The equivalent time-delayed Smoluchowski equation is known for the one-particle case Frank (2005a, b); Loos and Klapp (2017, 2019); Guillouzic et al. (1999); Atay (2010) and is briefly discussed in Sec. IV A. Subsequently, this single-particle case is extended to the case of many particles (in Sec. IV B) and the approximations used to obtain a self-contained DDFT equation are introduced in Sec. IV C. Finally, the system’s stability against traveling-waves is investigated within DDFT in Sec. IV D and the obtained predictions are compared to the simulation results in Sec. IV E.

IV.1 One-particle time-delayed Smoluchowski equation

To determine the time-delayed Smoluchowski equation equivalent to Eq. (1), we first revisit the case of just a single feedback-driven particle. In Frank (2005a), the Smoluchowski equation corresponding to a time-delayed Langevin equation was derived. Equations (1), (12), and (14) of this work are relevant here and are reproduced below for the case of a constant noise amplitude. Specifically, in Frank (2005a) it was shown that a time-delayed Langevin equation of the form

[TABLE]

for a general state variable $X(t)$ subject to a Gaussian noise $\Gamma(t)$, defined by $\langle\Gamma(t)\rangle=0$ and $\langle\Gamma(t)\Gamma(t^{\prime})\rangle=\delta(t-t^{\prime})$, is equivalent to the Smoluchowski equation

[TABLE]

with

[TABLE]

Here, the conditional probability $w(x,t|x_{\tau},t_{\tau})$ gives the probability for the system to be in state $x$ at time $t$, under the condition that it was in state $x_{\tau}$ at time $t_{\tau}=t-\tau$. The time derivative on the left-hand side of Eq. (15) only acts on $t$ but not $t_{\tau}$. The one-particle equivalent of our Eq. (1) is contained in this solution by identifying the actual and the time-shifted system states $x$ and $x_{\tau}$ with the particle positions $\mathbf{r}$ and $\mathbf{r}_{\tau}$. We set $h(\mathbf{r},\mathbf{r}_{\tau})=F(\mathbf{r},\mathbf{r}_{\tau})/\gamma$ for which

[TABLE]

where $\Delta$ denotes the Laplace operator. The Smoluchowski equation then reads

[TABLE]

From this equation, the joint probability $w(\mathbf{r},t;\mathbf{r}_{\tau},t_{\tau})=w(\mathbf{r},t|\mathbf{r}_{\tau},t_{\tau})w(\mathbf{r}_{\tau},t_{\tau})$ is obtained by multiplication with $w(\mathbf{r}_{\tau},t_{\tau})$ as

[TABLE]

Finally, integration over the past position $\mathbf{r}_{\tau}$ yields an expression for the probability $w(\mathbf{r},t)=\int\mathrm{d}\mathbf{r}_{\tau}w(\mathbf{r},t;\mathbf{r}_{\tau},t_{\tau})$ to find a particle at $\mathbf{r}$ at time $t$ without specifying its past position Frank (2005a, b); Loos and Klapp (2017, 2019); Guillouzic et al. (1999); Atay (2010); i.e., it leads from the joint probability to have a particle at position $\mathbf{r}$ at time $t$ and at position $\mathbf{r}_{\tau}$ at time $t_{\tau}$ to one for the joint probability to have a particle at position $\mathbf{r}$ at time $t$ and at any position at time $t_{\tau}$. The resulting probability $w(\mathbf{r},t)$ is then no longer dependent on the past time $t_{\tau}$ and its time evolution is given by

[TABLE]

IV.2 Extension to the many-particle case

We now return to the original problem of finding a Smoluchowski equation equivalent to Eq. (1). In this case, the system state $X(t)$ is given by a set of particle positions which should still obey Eq. (15). The many-body conditional probability density $w(\mathbf{r}^{N},t|\mathbf{r}^{N}_{\tau},t_{\tau})$ to have $N$ particles at positions $\mathbf{r}^{N}=\mathbf{r}_{1},\cdots,\mathbf{r}_{N}$ at time $t$ under the condition that the previous corresponding positions were $\mathbf{r}^{N}_{\tau}=\mathbf{r}_{\tau_{1}},\cdots,\mathbf{r}_{\tau_{N}}$ at time $t_{\tau}=t-\tau$ is thus given by

[TABLE]

Here, the gradient or divergence $\nabla_{i}$ indicates differentiation with respect to $\mathbf{r}_{i}$. The forces on the particles are expressed as potential gradients $-\nabla_{i}U_{\mathrm{tot}}(\mathbf{r}^{N},\mathbf{r}^{N}_{\tau})$ with the derivative in this term only intended to act on $U_{\mathrm{tot}}$ and not on $w(\mathbf{r}^{N},t|\mathbf{r}^{N}_{\tau},t_{\tau})$. The potential includes the direct pair interaction potential between the particles, $\Phi$, and a contribution due to the feedback potential, $V_{\mathrm{fb}}$:

[TABLE]

Analogously to the one-particle case, multiplication of Eq. (21) with the probability $w(\mathbf{r}_{\tau}^{N},t_{\tau})$ to have had $N$ particles at the positions $\mathbf{r}_{\tau}^{N}$ at time $t_{\tau}$ leads to an equation for the joint probability $w(\mathbf{r}^{N},t;\mathbf{r}_{\tau}^{N},t_{\tau})=w(\mathbf{r}^{N},t|\mathbf{r}_{\tau}^{N},t_{\tau})w(\mathbf{r}_{\tau}^{N},t_{\tau})$, given by

[TABLE]

Integration over the past positions $\mathbf{r}_{\tau}^{N}$ yields an equation for the probability density solely dependent on the set of current positions $\mathbf{r}^{N}$ given by

[TABLE]

Performing the integration, the potential term in Eq. (23),

[TABLE]

cannot be traced back to a simple expression due to the fact that both $U_{\mathrm{tot}}(\mathbf{r}^{N},\mathbf{r}_{\tau}^{N})$ and $w(\mathbf{r}^{N},t|\mathbf{r}_{\tau}^{N},t_{\tau})$ depend on the previous particle positions $\mathbf{r}_{\tau}^{N}$. Intuitively, this can be understood as being due to the coupling of $U_{\mathrm{tot}}(\mathbf{r}^{N},\mathbf{r}_{\tau}^{N})$ to the particle trajectories. This renders the reduction of the potential to a single value at the current position impossible. The many-particle equivalent of Eq. (20) obtained from integration of Eq. (23) thus reads

[TABLE]

Next, we define the one- and two-particle densities, $\rho$ and $\rho^{(2)}$, as well as the time-shifted two-particle density, $\rho^{(2)}_{\mathrm{s}}$, via

[TABLE]

Here, the instantaneous densities $\rho\left(\mathbf{r}_{1},t\right)$ and $\rho^{(2)}\left(\mathbf{r}_{1},\mathbf{r}_{2},t\right)$ do not depend on $t_{\tau}$ due to integration over the past positions $\mathbf{r}_{\tau}^{N}$. Furthermore, note that $\rho^{(2)}_{\mathrm{s}}(\mathbf{r}_{1},t;\mathbf{r}_{\tau},t_{\tau})$ includes a self-term (first part) as well as a term due to other particles (second part) such that the resultant two-particle density is defined as the probability to find one particle at $\mathbf{r}_{1}$ at time $t$ and any particle (i.e., the same or another) at position $\mathbf{r}_{\tau}$ at time $t_{\tau}$.

To obtain an equation for the evolution of the one-particle density, Eq. (26) is integrated with $N\int\mathrm{d}\mathbf{r}_{2}\cdots\int\mathrm{d}\mathbf{r}_{N}$. Using the above definitions, this leads to the exact equation

[TABLE]

It is worth mentioning that the last term, i.e., the interaction via the feedback potential, is similar to the preceding one, meaning that it takes the same form as a pair potential. The two terms differ in that the feedback does not act between positions at the same time but positions at time $t$ and previous ones at time $t_{\tau}$. Also, the last term includes a self-term which is not present in the direct interactions.

IV.3 Derivation of the DDFT equation

Starting from the many-body time-delayed Smoluchowski equation [Eq. (30)], which is stochastically equivalent to the Langevin equation [Eq. (1)], we construct a dynamical density functional theory (DDFT) Archer and Evans (2004). The DDFT is then used in the next section to investigate the stability of the homogeneous density state against wave formation.

Equation (30) is exact but not self-contained as it includes the instantaneous and the time-shifted two-particle density. To obtain a closed equation, we first rewrite the direct interaction potential using the adiabatic approximation of DDFT,

[TABLE]

which expresses the particle interaction forces via the direct correlation function $c^{(1)}(\mathbf{r})$. The latter is related to the excess free energy, i.e., the nonideal part of the free energy, $\mathcal{F}_{\text{exc}}$, by

[TABLE]

in equilibrium. The adiabatic approximation assumes the same relation in nonequilibrium.

Second, it is necessary to focus on the time-shifted two-particle density. In fact, even for the noninteracting particle case ($\Phi=0$) this term leads to a hierarchy of equations whose solution is not at all trivial and subject to current research Frank (2005a, b); Loos and Klapp (2017, 2019); Guillouzic et al. (1999); Atay (2010). Here, we use the mean-field-like approximation

[TABLE]

At this point, we are thus neglecting the correlations between the present and past positions of the particles. The description is then equivalent to having particles move in an effective external potential

[TABLE]

which is obtained through a convolution of the time-shifted one-particle density with the feedback potential $V_{\mathrm{fb}}$. This is a crude approximation for systems where the different sets of actual particle paths are important and the one-particle density is not sufficient to describe the system. It should, however, be decent if we are close to a homogeneous state due to a small amplitude of the feedback potential or if we consider only a small perturbation to this state, as done in the next section. Such a homogeneous reference state is expected to lose its memory of the specific realization of past positions on the time scale $1/D\rho_{0}$, i.e., the time needed for particles to diffuse the mean interparticle distance. However, this approximation is in particular not well justified for dilute feedback systems or very short delay times $\tau\ll 1/D\rho_{0}$.

Finally, inserting the above approximations into Eq. (30), a self-contained approximative equation for the one-particle density is obtained:

[TABLE]

IV.4 Linear stability analysis

Next, the effect of a small perturbation $\rho_{\mathrm{p}}(\mathbf{r},t)$ on a homogeneous state of constant density $\rho_{0}$ is investigated with

[TABLE]

being the space- and time-dependent density. We expand the direct correlation function $c^{(1)}(\mathbf{r})$ about the bulk fluid value $\rho_{0}$ up to linear order,

[TABLE]

with the direct pair-correlation function

[TABLE]

and insert this expression into Eq. (35). Linearization of the result in the perturbation density $\rho_{\mathrm{p}}(\mathbf{r},t)$ and nondimensionalization yields the equation

[TABLE]

where $\tilde{\rho}\left(\mathbf{k},t\right)$ indicates the Fourier transform of $\rho_{\mathrm{p}}(\mathbf{r},t)$, $\tilde{c}(k;\rho_{0})$ the one of $c^{(2)}\left(|\mathbf{r}-\mathbf{r}_{2}|;\rho_{0}\right)$, and $\tilde{V}_{\mathrm{fb}}$ the one of the feedback potential.

Taking the ansatz for a wave solution of wave vector $\mathbf{k}$ and amplitude $\varepsilon$,

[TABLE]

and inserting it into Eq. (39), we obtain the dispersion relation

[TABLE]

The solution for $\lambda$ can be separated into a real part $\alpha$ describing the growth of the perturbation and an imaginary part $\omega$ giving the angular frequency of the traveling wave as

[TABLE]

This gives the two equations

[TABLE]

which implicitly define the solution. Since $\tilde{c}(k;\rho_{0})$ for the equilibrium bulk fluid is required here as an input, we use a reference simulation without the feedback potential to obtain $\tilde{c}(k;\rho_{0})$ which is related to the structure factor of the system by Hansen and McDonald (1990)

[TABLE]

IV.5 DDFT results and comparison to simulations

Results for $\alpha(k)$ and $\omega(k)$ are shown in Fig. 9. Figure 9(a) displays the range of unstable wave numbers defined via a positive $\alpha(k)$ and puts these into relation to the structure factor $S(k)$. Indeed, the maximally unstable wave number, i.e., the wave number at which $\alpha(k)$ has its maximum, coincides with the structural bandwidth characterized by $k^{*}$.

Figure 9(b) shows the full dispersion $\alpha(k)$ and $\omega(k)$ for various $A$. The range of unstable wave numbers is zero below the transition ($A<4.8$) and increases with the feedback amplitude $A$ above it. From the frequency $\omega(k)$, close to the transition, the phase velocity of the imposed wave is determined by $\omega(k)/k$. At $k=k^{*}$ this phase velocity should be close to the band velocity $v$ such that the band velocity predicted by the theory is

[TABLE]

According to Eqs. (6) and (45) we obtain an approximative phase shift of $\omega(k^{*})\tau=\pi$ within the delay time. In fact, Fig. 9(b) reveals that the phase shift $\omega(k^{*})\tau$ is close to $\pi$ but deviations appear due to the approximative nature of Eq. (6). The phase shift $\omega(k^{*})\tau$ is close to the ideal value $\pi$ for any $k$ close to the maximally unstable wave vector with discrepancies indicating that the wave does not change exactly between a pattern and its negative image after a delay time $\tau$.

The theoretical prediction of the band velocity is compared to the simulation data in Fig. 9(c) and reasonable agreement is found with respect to onset and magnitude of the velocity as a function of the feedback potential amplitude $A$.

Figure 9(d) illustrates the stability in the plane defined by the wave number $k$ and the feedback amplitude $A$. There is a separatrix distinguishing a stable [$\alpha(k)<0$] from an unstable [$\alpha(k)>0$] regime. The instability occurs beyond a critical amplitude $A^{*}\approx 4.8$ with an instability wave number $k^{*}=0.87$. The DDFT prediction for the critical amplitude $A^{*}\approx 4.8$ compares favorably with the simulation results of $5\leq A^{*}\leq 6$ which documents the predictive power of the microscopic theory. Deviations may be due to the approximative pair correlation function entering as an input for the nonequilibrium steady state or due to the adiabatic approximation.

In the unstable regime, $\alpha(k)>0$, the maximally unstable wave number, i.e., the one with the fastest growth rate, is shown by the blue dashed line in Fig. 9(d). This prediction can be compared to the actual band wave number $k^{*}$ found in the simulations shown as green dots. While at the onset of the transition these two wave numbers are very close [see also Fig. 9(a)], they exhibit a different trend for increasing $A$ indicating that nonlinear effects play a larger role away from criticality.

V Circular confinement

In the case of periodic boundary conditions as discussed in Sec. III the commensurability of the box size and the (projected) wavelength effectively constrains the possible band orientations. Inversely, this constraint may affect the actually realized wavelength as certain combinations of orientation and wavelength that fit well with the box size may be preferred. In the following, we aim to investigate this effect in more detail by confining the system to a prescribed geometry. We choose a circular confinement which is imposed onto the system by an external potential $\Phi_{\mathrm{c}}(r)$ which is radially symmetric. The potential is soft but diverges at a distance $r=R$ from the center such that a circular confinement is realized. We choose the external potential form

[TABLE]

and fix the radius to $R\approx 45\,b$ such that the density $\rho_{0}=N/\pi R^{2}=1/b^{2}$ with $N=6400$ is the same as in the previous case.

Representative snapshots of the system for different feedback potential amplitudes $A$ are shown in Fig. LABEL:fig_circularconfined. Band formation is found as for the nonconfined case, with the pattern reaching further inwards if the feedback potential strength is increased. This is because the confining boundaries set a preferred band orientation and thus facilitate the formation of a traveling wave along the boundary. Due to the circular form of the confining potential, the bands are now forced into a rotational motion. Moreover, the circular geometry does not allow parallel wavefronts with a constant separation running from the boundary to the center. Thus, defects are unavoidable. At low and intermediate feedback potentials the formed pattern is approximately rotationally symmetric with respect to the center but not very well ordered. In contrast, at high values of $A$, the waves are well-established and have an about constant wavelength. However, to allow for this arrangement, a separation into domains of different band orientation with a complicated topology at the center is required and observed.

To further examine the system response in this case, we consider the time-averaged radial density profile

[TABLE]

which by symmetry depends only on the radial distance $r$ to the center. Next, we also define a time-averaged angular drift velocity profile $\omega_{\mathrm{c}}(r)$:

[TABLE]

where $\hat{\mathbf{e}}_{z}$ is the unit vector perpendicular to the plane of motion. The profile $\omega_{\mathrm{c}}(r)$ likewise depends on the radial distance to the center only. Figure 11 shows results for the density profile $\rho_{\mathrm{c}}(r)$ and the angular velocity $\omega_{\mathrm{c}}(r)$.

For growing feedback potential strength $A$, an accumulation of particles near the system boundary is observed as signaled by a peak in the density profile close to the wall accompanied by a (relative) depletion of particles near the center. This effect has also been found for self-propelled particles in circular confinement Smallenburg and Löwen (2015); Jamali and Naji (2018) and can qualitatively be understood here on similar grounds as suggested by the MSD (Fig. 8) which resembles active Brownian particles. The difference between the two cases is that here, the individual particle velocity is typically directed tangentially to the wall. Therefore, in the feedback case the accumulation is driven by wall curvature and larger $A$ enhances the wall accumulation effect. In contrast, in the active particle case it is driven by persistence in the otherwise unconstrained propulsion direction.

Considering the angular velocity, we find that for small feedback potentials it is almost zero at all distances $r$ from the center. Increasing $A$ leads to the formation of a two-layer state with the inner part revolving in the opposite direction to the outer one, signaled by a sign change in $\omega_{\mathrm{c}}$. The continuous change in $\omega_{\mathrm{c}}$ with distance $r$ is consistent with a motion spiraling outwards. At even higher values of $A$, this feature disappears and all particles rotate like a rigid body with a constant joint angular velocity apart from some remaining distance-dependence in the central part of the profile. Thus, the slope of $\omega_{\mathrm{c}}$ is negative, implying a motion which is spiraling inwards. Furthermore, close to the wall, particles move in the opposite direction to the bands or at least significantly more slowly than the band. This is especially observed for intermediate feedback strengths, signaled by a dip in the angular velocity to values below zero. This effect is interpreted as an escape of particles from the ideal bands induced by wall curvature, leading to counterpropagating particles which are visible in the system snapshots of Figs. LABEL:fig_circularconfined(b) and (c). It vanishes again for strong feedback strengths [see Fig. LABEL:fig_circularconfined(d)].

For high feedback potential amplitude, the emergence of domains with different band orientations is observed. These domains are separated by grain boundaries which, interestingly, are neither static nor co-rotating with the same mean angular speed as the particles. Instead, they rotate with a considerably slower angular speed albeit in the same direction. This is documented in Fig. LABEL:fig_snapshots_phases, where typical system snapshots are shown during one full but slow rotation of the system grain boundaries. Four boundaries separating different domains are indicated in Fig. LABEL:fig_snapshots_phases.

The angular position of the grain boundaries as a function of time is shown in Fig. 13. Their angular positions change approximately linearly in time. A linear fit yields an average rotation speed of the grains of $\omega_{\mathrm{g}}\approx-0.04$ for the system of Fig. LABEL:fig_circularconfined(d), i.e., for a feedback amplitude of $A=20$. This is one order of magnitude smaller than the mean particle angular velocity which is $\omega_{\mathrm{c}}\approx-0.5$ in the outermost part of the cavity [see Fig. 11(b)]. We speculate that the slower grain boundary speed is related to the group velocity of propagating bands in the bulk which is also much smaller than the phase velocity. In fact, from Fig. 9(b), the group velocity $d\omega(k)/dk|_{k=k^{*}}$ is only 13% of the phase velocity $\omega(k^{*})/k^{*}$ at the wave number $k^{*}$ for the given potential strength ($A=20$).

To summarize, the feedback-driven suspension forms moving bands also under circular confinement. Due to the impenetrability of the walls, the band velocity has to be tangential to the wall; i.e., the band lamellae must be oriented perpendicular to the wall if the bands move along the normal of their interface as in the unconstrained case. While confinement to a straight (i.e., uncurved) slit geometry supports band formation perpendicular to the walls by imposing a preferred band orientation (data not shown), wall curvature destabilizes such bands due to the incompatibility with straight bands. Since the constraint of perpendicularity cannot be sustained for a straight band along a curved wall, the bands must either bend or change their width as a function of radial distance to the center such that $\lambda_{\mathrm{b}}\equiv\lambda_{\mathrm{b}}(r)$. Furthermore, as in the periodic boundary case, the fixed length of the circumference also poses a commensurability constraint for the band wavelength. However, this second constraint is weakened due to the grain boundaries. For the large feedback strengths considered here, the system exhibits tilted band formation. Dynamically this implies that spiraling waves are formed for which both inward and outward orientation are observed. In particular, particles near the curved wall cannot accommodate their individual motion with the global band motion and depending on $A$ are either left behind or propagate faster [see Fig. 11(b)]. As a further important consequence, the convex curvature of the boundary leads to a significant particle accumulation near the wall and a depletion from the system center. The latter constitutes a natural source of disorder as topologically the cavity can never fulfill the constraints of curved bands meeting in the center. Instead the center becomes a nucleation point for grains with different band orientations. For strong feedback strengths, these grains co-rotate with much smaller angular velocities than that of the individual bands.

VI Conclusions

We have shown that subjecting colloidal particles to a repulsive feedback potential dependent on their previous positions can lead to pattern formation. For a sufficiently strong repulsive potential, we have found a transition to a traveling band phase which can be predicted by a linear stability analysis of dynamical density functional theory. Under circular confinement, the transition persists but becomes more complex, exhibiting wall accumulation of particles, spiraling patterns, and creation of banded domains. Our findings can be verified in experiments by using colloids in an external light field Hanes et al. (2009); Evers et al. (2013); Jones et al. (2015); Bewerunge and Egelhaaf (2016); Bäuerle et al. (2018) or employing autochemotactic particles Jin et al. (2017, 2018).

Future work can be performed along the following directions: First, it is worth considering an attractive feedback potential with a negative amplitude $A<0$. For attractive feedback, particles tend to stay at their previous positions. Starting from an initial homogeneous system, this leads to the formation of particle clusters accompanied by a slowing down of the dynamics and to subsequent coarsening. The emerging structures are expected to be similar to those found in phase separation Tanaka (2000). It would thus be interesting to identify and compare the scaling laws for the mean cluster size as a function of time in different regimes for various system parameters.

Second, while our scenario will not qualitatively change for other types of soft repulsion between the particles, a long-ranged attraction Olson Reichhardt et al. (2010); Zhao et al. (2012) in the interparticle forces can give rise to equilibrium gas-liquid phase separation Vliegenthart and Lekkerkerker (2000) with a critical point. If a repulsive feedback potential is applied, there is competition between bulk phase separation and band formation which may lead to interesting new structures.

Third, if one considers higher particle densities or lower temperature than in this paper, the equilibrium two-dimensional Yukawa system exhibits freezing into a hexatic or hexagonal crystalline phase Kapfer and Krauth (2015). Applying a repulsive feedback is then expected to lead to the formation of traveling crystalline bands (“solids”), the detailed structure of which still has to be worked out. The density increase in the bands induced by the feedback will additionally support and enhance traveling crystal formation. Thus, it is expected that the phase boundary will depend on the amplitude $A$. Dynamical density functional theory can in principle be applied to describe crystallization in nonequilibrium van Teeffelen et al. (2008); Zimmermann et al. (2016).

Fourth, active particles in feedback potentials present a rich playground to investigate further complex collective effects due to competition between activity and feedback forces. Even in the absence of feedback the collective behavior of active particles is rich Klamser et al. (2018) and it remains to be explored how swarming Elgeti et al. (2015) and motility-induced phase separation Cates and Tailleur (2015) compete with feedback potentials. Again dynamical density functional theory can be employed for an appropriate instability analysis Menzel et al. (2016).

Last, possible future work should also include extending the system presented here to three dimensions and considering the effect of more complex confining geometries Beppu et al. (2017), for example, nonconvex walls and moving boundaries. Likewise, multiple (competing) Ahlborn and Parlitz (2005) or more complex Ciszak et al. (2015); Pyragas and Pyragas (2018) feedback terms can be considered. The latter may also take a density-dependent form, thus describing quorum sensing Bäuerle et al. (2018).

Acknowledgements.

We thank W. Zimmermann, B. Liebchen, C. Hoell and M. Tarama for helpful discussions. We gratefully acknowledge financial support by the DFG Grants No. LO 418/19-1 and No. EG 269/6-1.