Active chiral particles under confinement: surface currents and bulk accumulation phenomena

TL;DR

This study explores how chiral active particles behave under confinement, revealing surface currents and bulk accumulation phenomena influenced by chirality and confinement geometry, with analytical and numerical insights.

Contribution

It introduces an extended active Ornstein-Uhlenbeck model for chiral particles under confinement, analyzing their distribution and currents both analytically and numerically.

Findings

Chirality reduces wall accumulation compared to non-chiral particles.

Chiral particles generate currents parallel to confining walls.

Bulk accumulation increases with certain confinement conditions.

Abstract

In this work, we study the stationary behavior of an assembly of independent chiral active particles under confinement by employing an extension of the active Ornstein-Uhlenbeck model. The chirality modeled by means of an effective torque term leads to a drastic reduction of the accumulation near the walls with respect to the case without handedness and to the appearance of currents parallel to the container walls accompanied by a large accumulation of particles in the inner region. In the case of two-dimensional chiral particles confined by harmonic walls, we determine the analytic form of the distribution of positions and velocities in two different situations: a rotationally invariant confining potential and an infinite channel with parabolic walls. Both these models display currents and chirality induced inner accumulation. These phenomena are further investigated by means of a more realistic description of a channel, where the wall and bulk regions are clearly separated. The corresponding current and density profiles are obtained by numerical simulations. At variance with the harmonic models, the third model shows a progressive emptying of the wall regions and the simultaneous enhancement of the bulk population. We explain such a phenomenology in terms of the combined effect of wall repulsive forces and chiral motion and provide a semiquantitative description of the current profile in terms of an effective viscosity of the chiral gas.

Equations

$$\gamma {\dot{\mathbf{x}}}_i$$

$${\dot{\mathbf{e}}}_i$$

$$⟨{\mathbf{e}}_i(t)\cdot {\mathbf{e}}_j(0)⟩={\delta }_{ij}{e}^{-D_{r}t}\cos (\Omega t),$$

$${\dot{\mathbf{u}}}_i=-\frac{{\mathbf{u}}_i}{\tau }+\Omega \hat{\mathbf{z}}\times {\mathbf{u}}_i+\frac{\sqrt{2D_a}}{\tau }{\mathbf{\eta }}_i,$$

$$⟨u_{\alpha }(t)u_{\alpha }(0)⟩=\frac{D_a}{\tau }{e}^{-\frac{t}{\tau }}\cos (\Omega t),$$

$$v_0^2=\frac{D_a}{\tau }d.$$

$$\dot{\mathbf{x}}=\mathbf{v}+\sqrt{2D_t}\mathbf{\xi },$$

$$\dot{\mathbf{v}}=-\frac{1}{\tau }\left(\frac{1}{\gamma }\nabla \varphi +\mathbf{\Gamma }(\mathbf{x})\cdot \mathbf{v}\right)+\Omega \hat{\mathbf{z}}\times \left(\mathbf{v}+\frac{1}{\gamma }\nabla \varphi \right)$$

$$+\frac{\sqrt{2D_a}}{\tau }\mathbf{\eta }+\frac{\mathbf{I}-\mathbf{\Gamma }(\mathbf{x})}{\tau }\sqrt{2D_t}\mathbf{\xi },$$

$${\Gamma }_{\alpha \beta }(\mathbf{x})={\delta }_{\alpha \beta }+\frac{\tau }{\gamma }{\nabla }_{\alpha }{\nabla }_{\beta }\varphi (\mathbf{x}).$$

$$\frac{\partial P(\mathbf{x},\mathbf{v},t)}{\partial t}=\frac{D_a}{{\tau }^2}\sum _{\alpha }\frac{{\partial }^{2}P(\mathbf{x},\mathbf{v},t)}{\partial v_{\alpha }^2}-\sum _{\alpha }{v}_{\alpha }\frac{\partial P(\mathbf{x},\mathbf{v},t)}{\partial x_{\alpha }}+$$

$$\sum _{\alpha }\frac{\partial }{\partial v_{\alpha }}\frac{1}{\tau }\left(\sum _{\beta }({\delta }_{\alpha \beta }+\frac{\tau }{\gamma }{\varphi }_{\alpha \beta }-\Omega \tau {ϵ}_{\alpha \beta })v_{\beta }P(\mathbf{x},\mathbf{v},t)\right)$$

$$+\sum _{\alpha }\frac{\partial }{\partial v_{\alpha }}\frac{1}{\tau }\left(\sum _{\beta }({\delta }_{\alpha \beta }-\Omega \tau {ϵ}_{\alpha \beta })\frac{{\varphi }_{\beta }}{\gamma }P(\mathbf{x},\mathbf{v},t)\right),$$

$${\mathbf{div}}_{\alpha }\cdot {\mathbf{I}}_{\alpha }=\sum _{\alpha }\frac{\partial }{\partial x_{\alpha }}(I_{\alpha }{)}_{\mathbf{x}}+\frac{\partial }{\partial v_{\alpha }}(I_{\alpha }{)}_{\mathbf{v}}.$$

$$(I_{\alpha }^{rev}{)}_{\mathbf{x}}=v_{\alpha }{P}_{st}(\mathbf{x},\mathbf{v}),$$

$$(I_{\alpha }^{rev}{)}_{\mathbf{v}}=-\frac{1}{\tau \gamma }\left(\sum _{\beta }({\delta }_{\alpha \beta }-\Omega \tau {ϵ}_{\alpha \beta }){\varphi }_{\beta }\right)P_{st}(\mathbf{x},\mathbf{v})=0,$$

$$(I_{\alpha }^{irr}{)}_{\mathbf{x}}=0,$$

$$(I_{\alpha }^{irr}{)}_{\mathbf{v}}=-\frac{D_a}{{\tau }^2}\frac{\partial P_{st}(\mathbf{x},\mathbf{v})}{\partial v_{\alpha }}-\frac{1}{\tau }\sum _{\beta }({\delta }_{\alpha \beta }+\frac{\tau }{\gamma }{\varphi }_{\alpha \beta }-\Omega \tau {ϵ}_{\alpha \beta })v_{\beta }{P}_{st}(\mathbf{x},\mathbf{v}).$$

$$P_{st}(\mathbf{x},\mathbf{v})=\mathcal{N}\exp \left(-\frac{\tau }{2D_a}\left(v_x^2+v_y^2\right)\right).$$

$$⟨(x_{\alpha }(t)-x_{\alpha }(0){)}^2⟩=2{\int }_0^{t}d{t}_1{\int }_0^{t_1}d{t}_2⟨u_{\alpha }(t_1)u_{\alpha }(t_2)⟩$$

$$\displaystyle=2D_{\Omega}\Bigl{(}t+\frac{(1-\Omega^{2}\tau^{2})\tau}{1+\Omega^{2}\tau^{2}}\left(e^{-t/\tau}\cos(\Omega t)-1\right)$$

$$\displaystyle-\frac{2\Omega\tau^{2}}{1+\Omega^{2}\tau^{2}}e^{-t/\tau}\sin(\Omega t)\Bigr{)},$$

$$(I_x^{irr}{)}_{\mathbf{v}}=-\Omega v_y{P}_{st}(\mathbf{v})\ne 0,$$

$$(I_y^{irr}{)}_{\mathbf{v}}=\Omega v_x{P}_{st}(\mathbf{v})\ne 0.$$

$$\varphi (\mathbf{x})=k\frac{(x^2+y^2)}{2}.$$

$$P(\mathbf{x},\mathbf{v})=\mathcal{N}\exp \left[-\frac{\tau }{2D_a}\Gamma \left((v_x+\frac{\Gamma -1}{\Gamma }\Omega y{)}^2+(v_y-\frac{\Gamma -1}{\Gamma }\Omega x{)}^2\right)\right]$$

$$\times \exp \left[-\frac{1}{2D_a\tau }\frac{\Gamma -1}{\Gamma }\left({\Gamma }^2+{\Omega }^2{\tau }^2\right)(x^2+y^2)\right],$$

$$⟨\mathbf{v}{⟩}_{\mathbf{x}}=\frac{\int d\mathbf{v}P(\mathbf{x},\mathbf{v})\mathbf{v}}{\int d\mathbf{v}P(\mathbf{x},\mathbf{v})},$$

$$⟨v_x{⟩}_{\mathbf{x}}$$

$$⟨v_y{⟩}_{\mathbf{x}}$$

$$⟨v_x^2⟩-⟨v_x{⟩}^2=⟨v_y^2⟩-⟨v_y{⟩}^2=\frac{1}{\Gamma }\frac{D_a}{\tau }$$

$$\rho (x,y)={\mathcal{N}}^{\prime }\exp \left(-\frac{1}{2D_a\tau }\frac{\Gamma -1}{\Gamma }\left({\Gamma }^2+{\Omega }^2{\tau }^2\right)(x^2+y^2)\right).$$

$$2p_{v}A=⟨\nabla \varphi \cdot \mathbf{x}⟩,$$

$$p_{v}A=\frac{D\gamma }{\Gamma +\frac{{\Omega }^2{\tau }^2}{\Gamma }}.$$

$$P(\mathbf{x},\mathbf{v})\approx \mathcal{N}\exp \left[-\frac{\tau }{2D_a}\left(\Gamma v_x^2+v_y^2+2\frac{\Gamma -1}{\Gamma }\Omega \tau v_x{v}_y\right)\right]$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: a* Gran Sasso Science Institute (GSSI), Via. F. Crispi 7, 67100 L’Aquila, Italy. E-mail: lorenzo.caprini@gssi.it*††footnotetext: *b**School of Sciences and Technologies, Università di Camerino, via Madonna delle Carceri, 62032, Camerino, and INFN Perrugia, Italy. E-mail: umberto.marinibettolo@unicam.it *

1 Introduction

The dynamics of self-propelled particles is an exciting new area of research at the frontier between physics, biology and bioengineering 1, 2, 3, 4.These particles are ubiquitous in nature or can be artificially created and employed in a near future in tasks such as targeted drug delivery and nanosurgery. Examples of the first category are bacteria, ciliates, Synechococcus a type of blue-green alga, sperm cells, while artificial self-propellers are represented by Janus particles catalytically driven or magnetic microparticles subject to a magnetic field. Each particle, whose size may range between 10 $\mu m$ up to several microns, consumes energy from the environment in order to move and is fueled either either by metabolic processes or by chemical reactions: E. coli move forward by a rotational motion of their spiral-shaped flagella while synthetic Janus colloids are driven by a catalytic chemical reaction on their surface.

It is sufficient a small departure from the right-left symmetry relative to the propulsion axis of self-propelled particles or an external magnetic field to determine circle swimming in two dimensions or helical swimming in three dimensions, a phenomenon termed as chirality or handedness. Indeed, in nature one observes spiral-like swimming trajectories in E. coli bacteria 5 and spermatozoa in bulk suspensions 6 or circle-like motion near a planar substrate 7 and clockwise treadmilling in FtsZ proteins on membranes 8.

Analogous chiral trajectories are produced by artificial microswimmers, for instance, L-shaped particles 9, 10. In nature, magnetotactic bacteria are widespread motile prokaryotes, endowed with an organelle called magnetosome which causes the cell to behave like a miniature compass and swim parallel to the magnetic field lines  11, 12. The direction of the magnetic field is so crucial that this kind of bacteria die when taken to the opposite hemisphere of Earth.

As discussed by ten Hagen et al. 13 a system of self-propelled particles can be described at two different levels: i) a fine-grained level fully taking into account the self-propulsion mechanism, the degrees of freedom of the microswimmers, such as moving internal elements, and the solvent hydrodynamics 14, ii) a coarse-grained level where the motion of microswimmers is modeled by using the Brownian overdamped equations and the self-propulsion is represented in terms of effective forces and torques 15. In the following, we have found convenient, following the literature in the field, to adopt the second level of description which allows a straightforward application of the tools of statistical mechanics and computer simulation to large collections of active particles. Among these tools there are two popular theoretical approaches to the study of active particles: Active Brownian particle (ABP) model 16, 17, and the Active Ornstein-Uhlenbeck particle (AOUP) model 18. Both models represent the microswimmer motion through a viscous solvent by means of Langevin overdamped dynamics and include the distinguishing feature of active particles, namely the self-propulsion mechanism as an effective force term. Such a self-propulsion force has the following properties: a) it is typically noisy and its orientation is random, b) active forces acting on different self-propelled particles are mutually independent and c) the direction of driving persists over time scales of the order of the microseconds. As briefly discussed in the next section, the ABP and AOUP differ by the way the self-propulsion force is modeled, but in spite of that many important aspects of their nonequilibrium behavior are similar. So far, the mechanism responsible for the active chiral behavior has been included only in the ABP by adding an effective constant torque term in the corresponding Langevin equation 15. Such a torque induces a systematic drift in the orientation of active particles, i.e. a rotation with a constant angular speed of the active velocity, leading to helical trajectories in three dimensions and to circular ones in two dimensions. While the effect of chirality on the bulk behavior of self-propelled particles has been investigated 10, the properties of self-propelled chiral particles under confinement are not so well known and understood. Confining surfaces may induce surface accumulation, sliding motion and selection of active particles according to their chirality 19, 20, 21 and therefore it is interesting to investigate under which conditions these effects can be enhanced or suppressed.

In the present work, we focus on the properties of confined chiral particles and motivated by the simpler mathematical structure of the AOUP, which facilitates the theoretical treatments, we extend it to account for chiral effects. As in the case of the ABP, this goal is achieved by adding a torque term to the self-propulsion.

The paper is organized as follows: In Sec. 2 we introduce the chiral active Ornstein-Uhlenbeck particle (CAOUP) model, in Sec. 3 we study the special case of CAOUP confined in a two-dimensional harmonic trap, in Sec. 4 we study the behaviour of the model in the case of confinement in a parabolic channel and in a slit. Finally, in Sec 5 we present the conclusions and discuss some future directions.

2 Planar chiral active motion

In the following, we shall describe an assembly of $N$ mutually independent active chiral particles moving in the two dimensional $(x,y)$ plane. Each particle at position ${\mathbf{x}}_{i}$ is subject to four kinds of drivings: an active or self-propulsion force, ${\mathbf{f}}^{a}_{i}$, a drag force, $-\gamma\dot{\mathbf{x}}_{i}$ due to the friction with the solvent, a white noise force $\gamma\sqrt{2D_{t}}\boldsymbol{\xi}_{i}$ representing the thermal agitation of the solvent, whose intensity depends on the thermal diffusion coefficient, $D_{t}$, and a force, $-\nabla\phi({\mathbf{x}}_{i})$ due to a potential, $\phi$ which confines their motion in a restricted region of the two-dimensional space. The active force, ${\mathbf{f}}^{a}_{i}=\gamma v_{0}{\mathbf{e}}_{i}$ has a constant magnitude proportional to a velocity $v_{0}$ and a direction ${\mathbf{e}}_{i}\equiv(\cos\theta_{i},\sin\theta_{i})$ dependent on an angle $\theta_{i}(t)$ undergoing unbiased rotational diffusion with rotational diffusion coefficient $D_{r}$. Inertial effects are neglected because of the low Reynolds number regime as well as hydrodynamic effects and the solvent reaction. Van Teffelen et al. 15 have extended the standard ABP model to account for the chiral motion of microswimmers by adding a new ingredient: they imposed a constant angular drift of amplitude $\Omega$ to the dynamics of the angle $\theta_{i}(t)$ to represent an effective constant torque uniformly applied to the particles. The resulting chiral ABP model is described by the following stochastic equations:

[TABLE]

where $\hat{\mathbf{z}}$ is the unit vector normal to such a plane, $\boldsymbol{\xi}_{i}$ and $w_{i}$ are independent Gaussian noises with $\delta$-correlated components, unit variance and zero mean. The torque turns the standard exponential form of the autocorrelation function of the orientation vector into a damped oscillatory behavior

[TABLE]

where $t>0$. Notice that when $|\Omega t|>\pi/2$ the force autocorrelation function becomes negative.

A series of interesting studies of the chiral ABP have been presented by various authors 19, 21, 22, 23, and mainly rely on numerical simulations and/or the analysis of the low noise, quasi-deterministic behavior of equations (1) and (2).

Hereafter, motivated by the success of the active Ornstein-Uhlenbeck model in reproducing and predicting the main behaviors of non-chiral microswimmers 24, 25, 26, we consider its chiral extension. Moreover, the AOUP is considered to be a valid alternative tool to investigate the properties of active particles because of the feature which makes it analytically more treatable than the ABP, namely the property that the fluctuations of the self-propulsion force are Gaussian. We introduce the chiral version of the AOUP model by assuming the same governing equation (1) for ${\mathbf{x}}_{i}$ as in the ABP, but writing the dynamics of the active force ${\mathbf{f}^{a}_{i}}=\gamma{\mathbf{u}_{i}}$ as:

[TABLE]

where $\mathbf{\eta}_{i}$ is a Gaussian noise and each component of $\mathbf{u}_{i}$ evolves according to an Ornstein-Uhlenbeck process of characteristic time, $\tau$ and strength $D_{a}$ and is subject to a tangential drift at a fixed frequency $\Omega$ around an axis orthogonal to the plane of motion. It is very simple to show that, in the case of freely moving particles ($\phi=0$), the two-time correlation functions of the chiral versions of ABP and AOUP are similar. A brief calculation gives the following two-time the average:

[TABLE]

where Greek indices denote Cartesian components and the particle index will be omitted from now. Although, it is known that even in the non-chiral case there are some important differences between the two models (e.g. the stationary distribution of $f^{a}_{\alpha}$ is circular in ABP and Gaussian in AOUP), previous studies 27, 28 have shown that they share many important aspects of their nonequilibrium behavior. It is not too unlikely to assume that this holds true even in the the case of angular drift. Thus, neglecting autocorrelation functions of the self-propulsion ${\mathbf{f}^{a}_{i}}(t)$ beyond second order and noticing that Eqs. (3) and (5) have the same functional form, we establish the following mapping between the parameters of the two models by making the correspondence:

[TABLE]

For theoretical work it is convenient to transform Eqs (1) and (4) into a system of equations involving the effective velocity, $\mathbf{v}$ defined as $\mathbf{v}=\mathbf{u}-\frac{\nabla\phi}{\gamma}$ and generalizing the transformation of 29, 30 we obtain a new set of equations (see appendix A):

[TABLE]

where $\mathbf{I}$ is the identity matrix and $\boldsymbol{\Gamma}$ is an effective friction tensor given by:

[TABLE]

Such a transformation maps the original problem of an active particle onto the dynamics of a fictitious passive particle, immersed in a heat bath of amplitude $\sqrt{2D_{a}}/\tau$. The fictitious particle experiences: i) a deterministic force proportional to the potential gradient (the first term in the r.h.s of Eq. (8)), ii) a Stokes force dependent on the second derivative of the potential (the second term in the r.h.s) , iii) an effective Lorentz force proportional to the torque $\Omega$ and orthogonal to the velocity $\mathbf{v}$ and iv) a term proportional to $\Omega$ and orthogonal to the potential gradient. Hereafter, the thermal diffusion coefficient, $D_{t}$ is set equal to zero, not only for simplifying the analytical work but also because it is often much smaller than $D_{a}$ 31. For the following analysis, it is convenient to write the Fokker-Planck equation (FPE) 32 for the $P(\mathbf{x},\mathbf{v},t)$ distribution corresponding to the dynamics (7) and (8). In the limit of vanishing thermal noise we write:

[TABLE]

where for notational convenience we have adopted an explicit Cartesian representation, using Greek indices to denote two dimensional vector components, introduced the symbols $\phi_{\alpha}$ and $\phi_{\alpha\beta}$ for the first and second partial derivatives of the potential, respectively, and the symbol $\epsilon_{\alpha\beta}$ for the antisymmetric $2\times 2$ tensor such that $\epsilon_{yx}=-\epsilon_{xy}=1$.

2.1 Detailed balance

At variance with the standard AOUP model ($\Omega=0$) in order to obtain the steady state solution of the FPE (10) we ought to solve a second order partial differential equation corresponding to the vanishing of the phase-space divergence of the probability current 33 ${\bf I}_{\alpha}=\left((I_{\alpha})_{\mathbf{x}},(I_{\alpha})_{\mathbf{v}}\right)$, where $\alpha=x,y$:

[TABLE]

It is convenient to split the probability current into a reversible and an irreversible contribution , i.e. ${\bf I}_{\alpha}={\bf I}^{rev}_{\alpha}+{\bf I}^{irr}_{\alpha}$, according to the their parity under time reversal, i.e. the components $(\mathbf{x},\mathbf{v})$ of the reversible (irreversible) part of the current are transformed under time-reversal in the same (opposite) way as the time derivative of $(\mathbf{x},\mathbf{v})$, respectively. Explicitly, the components of the probability current are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Let us remark that the $x$-component ($y$-component) of the current depends on the $y$-component ($x$-component) of the gradient of the potential.

The system is microscopically reversible when the detailed balance condition holds 34, 35. Such a condition is satisfied if ${\bf I}^{irr}_{\alpha}=0$ and ${\bf div}_{\alpha}\cdot{\bf I}^{rev}_{\alpha}=0$ in the steady state.

2.2 Unconfined chiral active motion

In the free case $\phi=0$, the time-independent solution of Eq. (10) is the steady state distribution, $P_{st}$, and is uniform in space and has the form:

[TABLE]

We notice that while the two-time properties of the velocity distribution depend on $\Omega$, as already seen in Eq. (5), the presence of chirality does not affect the shape of the Maxwell-like velocity distribution (16) with respect to the case $\Omega=0$, a situation which is modified when the particles are confined as we shall show below.

We, now, consider how the torque $\Omega$ modifies the diffusive properties of the particles and find that for both ABP and AOUP the mean square displacement is given by the formula 15:

[TABLE]

where $D_{\Omega}=\frac{D_{a}}{1+\Omega^{2}\tau^{2}}$ and the thermal contribution $2D_{t}t$ to the mean square displacement has been neglected. One sees that the long-time diffusion coefficient is decreased by a factor $(1+\Omega^{2}\tau^{2})^{-1}$ with respect to the non-chiral case 36, an effect which is explained by noticing that the particles perform cycloid trajectories.

An easy check of the lack of detailed balance is provided by the free case whose distribution function is given by (16). At variance with the AOUP with $\phi=0$, the condition that the components of the irreversible current must vanish is violated:

[TABLE]

3 Spontaneous circulation of totally confined active particles

We begin by considering a CAOUP moving in two dimensions and confined in a harmonic trap 37:

[TABLE]

The linearity of the Ornstein-Uhlenbeck process allows to obtain the exact form of the nonequilibrium steady state, whose phase-space distribution reads:

[TABLE]

where $\Gamma=(1+\frac{\tau}{\gamma}k)$ with $k>0$ and ${\cal N}$ is a normalization constant.

At variance with passive systems, the positional and velocity coordinates are correlated and the steady state is characterized by currents linearly increasing with $\Omega\neq 0$. In fact, the average velocity field at fixed position $\mathbf{x}$, defined as

[TABLE]

is given by the formulas:

[TABLE]

Such a field has the structure of a vortex centered at the minimum of the potential. The confining force pins to the origin the trajectories of different particles with the result of producing a single coherent macroscopic vortex, which on the contrary is absent in the free-case as evident from the distribution (16). Finally, the velocity variance

[TABLE]

is reduced with respect to the free case being $\Gamma>0$, but is torque independent.

We also consider the corresponding marginal (configurational) distribution, $\rho(x,y)=\int d\mathbf{v}P(\mathbf{x},\mathbf{v})$:

[TABLE]

We remark that the chirality determines a more concentrated distribution of particles near the bottom of the potential well, thus effectively increasing its stiffness. In other words, the chirality acts as a centripetal force whose strength is proportional to $\Omega^{2}\tau^{2}$. The non-chiral limit $\Omega\to 0$ is smoothly recovered and one finds the well known AOUP distribution in the quadratic trap.

By inserting the exact solution (20) in Eqs. (12)-(15) one finds that also in the case of harmonically confined chiral particles at variance with the non chiral AOUP the detailed balance condition is violated due to the presence of circulating currents. Due to this lack of detailed balance condition the unified colored noise approximation (UCNA) 38 fails even in this simple case of a chiral AOUP confined in a harmonic trap 24.

3.1 Virial Pressure

The virial pressure formula is obtained by equating the pressure exerted by the particles on the walls to the force per unit length that the walls exert on the particles 39, 40. Thus, we get

[TABLE]

where $A$ is the area and the average is performed using the steady state distribution. For the rotationally invariant harmonic trap we find using Eq. (24) the following expression of the virial pressure

[TABLE]

Such a result shows that the pressure decreases as $\Omega$ increases since the particles tend to be more localised near the bottom of the well with respect to the $\Omega=0$ case.

4 Confinement in a channel

Before delving into the study of the confined chiral gas, let us briefly summarize some recent results concerning the case $\Omega=0$: the self-propelled particles accumulate in the proximity of repulsive walls 41, 42. Such a phenomenon is more evident when the persistence time, $\tau$, increases and $D_{a}/\tau$ is large with respect to the thermal noise amplitude, $\gamma D_{t}$. In fact, a particle pushed by the self-propulsion force along the normal to the wall remains trapped there and escapes only when the propulsion direction changes, i.e. typically after a time $\tau$. This mechanism leads to the formation of a peak in the proximity of each wall and is reflected in Eq. (8) by the presence of a Stokes force term, $-\frac{1}{\tau}\Gamma(\mathbf{x})\cdot\mathbf{v}$, which opposes the motion and traps the particle near the walls. In the region where the potential gradients are negligible, instead, the particles undergo an underdamped motion and their density is almost uniform or smoothly varying in the presence of thermal noise ($D_{t}\neq 0$) 30.

Intuitively, we expect a reduction of the accumulation of the particles at the walls with chirality. In fact, the time, $t_{w}$, a particle spends in front a wall is a decreasing function of $|\Omega|$, as we argue by considering the form of the two-time velocity autocorrelation function Eq. (5), characterized by the two characteristic time-scales $\tau$ and $1/\Omega$. In the non-chiral limit, $\Omega=0$, the residence time at the wall is $t_{w}\approx\tau$, because after this time a particle typically inverts its self-propulsion and goes back to the bulk. On the other hand, in the case $|\Omega|>0$, the first value where the autocorrelation changes its sign is $t_{\Omega}=\pi/(2|\Omega|)$, the smallest zero of the cosine function in Eq. (5). Based upon this remark, we identify two regimes: I) $t_{\Omega}\gg\tau$, where the role of the chirality is negligible and does not affect the distribution of the particles in the channel, hence $t_{w}\approx\tau$; II) $t_{\Omega}\leq\tau$, where the correlation (5) changes sign for times shorter than $\tau$. In this case, we can estimate the typical time needed to change orientation and leave the wall as $t_{w}\sim t_{\Omega}$. Following Lee 43, we propose a coarse-grained description where the particles belong to two different populations: a bulk populations of $n_{b}$ members and a wall population of $n_{w}$ elements, characterized by residence times $t_{b}$ and $t_{w}$, respectively. Due to the permanent injection of energy, there is a continuous exchange of particles between the two populations, so that to achieve the steady state we must have $n_{w}/t_{w}=n_{b}/t_{b}$. We may conclude that the wall population decreases as $t_{w}$ decreases, that is when $\Omega$ grows.

In the following, we employ two different types of set-ups in order to assess quantitatively the effects of chirality in confined systems. We begin with the study of an infinite parabolic channel where the particles are confined only in the $x$-direction and free to move along the $y$-direction. We shall determine the exact full steady state distribution function and show the existence of steady momentum currents induced by the chirality.

We then consider the more realistic slit case, i.e. a model where the wall and bulk regions are clearly separated and the walls are modeled by means of truncated repulsive harmonic potentials. In this case of non-constant potential curvature, it is possible to observe the interplay between wall accumulation and wall depletion due to the competition between wall attraction and a chiral effective force pushing the particles away from the boundaries.

4.1 Parabolic channel

A one dimensional quadratic well $\phi(\mathbf{x})=k\frac{x^{2}}{2}$ mimics a channel limited by harmonic walls along the $x$-direction and unbounded along the $y$-direction. In this case, the distribution depends on the two components of the velocity $\mathbf{v}$, but only from the $x$-coordinate transversal to the channel. Its representation is a multivariate Gaussian distribution containing diagonal terms proportional to $(v_{x}^{2},v_{y}^{2},x^{2})$, but also three cross terms of type $(v_{x}v_{y},xv_{y},xv_{x})$. The exact expression of each of the six proportionality constants, featuring in the Gaussian, is a function of the control parameters, $k,\tau,\gamma,\Omega,D_{a}$ and is reported in appendix B. In the main text, to ease the presentation we write the more transparent form of the distribution obtained by neglecting contributions to the coefficients beyond the linear order in $\Omega$:

[TABLE]

By integrating over $v_{x}$ and $v_{y}$ we obtain the $x$-dependent average velocity along the $y$-direction as a function of $\Omega$

[TABLE]

The result (27) is valid to all orders in $\Omega$ and for $\Omega>0$ predicts the existence of a current parallel to the walls and directed along the downward $y$-direction for $x<0$ and upwards for $x>0$, while reversing the sign of $\Omega$ the current in the two halves of the channel changes sign. Instead, there is no net average current ($\langle v_{x}\rangle_{x}$ along the $x$-direction. Finally, using the exact form of the distribution given by Eq. (46) we obtain associated marginalized positional distribution function

[TABLE]

which shows that the density near the center is enhanced with respect to the $\Omega=0$ case as if an extra effective potential $\frac{\Gamma-1}{\Gamma}\Omega^{2}\tau^{2}\ x^{2}/2$ was pushing the particles towards the midpoint. In Fig.  1 we display two typical trajectories for two different values of $\Omega$, which become more and more localised as the chirality increases.

Unlike the rotationally invariant case, the variances of the velocity of the particles weakly depend on the value of $\Omega$ as shown by the exact solution of appendix B:

[TABLE]

[TABLE]

and

[TABLE]

where the subscript means that the average refers to the midpoint $x=0$.

4.2 The slit

We turn, now, to study a different type of channel where the stiffness of the walls and the width of the channel can be varied independently. A two-dimensional collection of independent chiral active particles is confined between two parallel repulsive soft walls, exerting a force piece-wise linear, characterized by an elastic constant $k$ according to the formula:

[TABLE]

where $\Theta$ is the Heaviside step-distribution. We choose a large value of $k$ so that the penetration inside the wall is negligible. The space between the walls extending from $x=-L$ to $x=L$, instead, forms a force-free region. The setup, recently studied in the case of non-chiral particles, $\Omega=0$, provides a more realistic description of a straight capillary because we may clearly distinguish a boundary, potential region from a bulk-like region 30. In the numerical simulations, in order to model an infinite vertical channel we assume periodic boundary conditions along the $y$-direction. At variance with the non chiral case, where the problem becomes effectively one-dimensional, in the $\Omega\neq 0$ case, as already found in the previous subsection 4.1, both $x$ and $y$ directions matter and the system develops a steady vertical momentum current and the bulk region becomes overpopulated with respect to the $\Omega=0$ case.

In Fig. 2, we plot the density at the walls, $n_{w}$, as a function of $\Omega$, for two different parameter configurations: we fix the ratio $D_{a}/\tau=1$ and perform numerical simulations for $\tau=1$ and $\tau=10$. We evaluate $n_{w}$ by counting the number of particles in the regions $x\leq-L$ and $x\geq L$, for the left and right wall, respectively. In the regime I) identified in section 4 the decrease of $n_{w}$ is slower than $1/|\Omega|$, since the effect of the chirality does not affect the wall population being $t_{\Omega}>\tau$. On the other hand, in regime II) the scaling $1/\Omega$ fairly agrees with the data for both choices of $\tau$.

In Fig.3 we show the marginalized distribution $\rho(x)$ for a system of active chiral particles subject to the dynamics of Eqs. (7)-(8) in the presence of the wall force (30) . We consider the dependence of the profile on the internal torque by varying $\Omega$ and keeping $D_{a}$ and $\tau$ fixed to values such that there is an appreciable accumulation at the walls. Upon introducing a small chirality the accumulation at the walls is reduced, as shown in Fig.3 (a). In an intermediate regime of $\Omega$, the role of the chirality is more consistent: on one hand, particles accumulate in front of the wall (as in the case $\Omega=0$), on the other hand, the profile in the bulk is no longer flat-like. In particular, we observe an emptying phenomenon in a layer near the wall favoring the accumulation in the inner region, as shown in Fig.3 (b). A further increase of $\Omega$ depletes the density near the walls and enhances it in the bulk region until the wall region remains almost completely empty. For these values of $\Omega$, the situation is completely inverted with respect to the case $\Omega=0$ and the wall behaves as if the effective wall-potential were repulsive.

We stress that, at variance with the harmonically confined CAOUP of Sec. 4.1, it is possible to observe a transition from a situation where particles accumulate near the walls to a situation where they accumulate in the bulk.

We, now, explain the process which leads to the emptying of the wall region and to the enrichment of the central region in the case of the channel when $\Omega\neq 0$ and to this purpose we rewrite Eq. (8) under the explicit form:

[TABLE]

One sees that on one side the chirality acts as an effective magnetic field of strength $\Omega$ periodically rotating the velocity direction without changing its magnitude and the other side, for $|x|>L$, in combination with the wall repulsion it produces a tangential force always along the $y$ direction. As a result, the particle accelerates and gains "kinetic energy". Such a force, proportional to $k\Omega/\tau$ only acts in the boundary region and generates a vertical momentum current having opposite directions at the two walls. When a particle penetrates the potential region, its velocity components $v_{x}$ is strongly damped by the large Stokes force $-\Gamma v_{x}/\tau$ and thus it remains trapped there. At the same time, the $v_{y}$ component of the velocity rapidly increases under the action of the tangential force. This process continues for a time $t_{\Omega}$ after which the torque rotates the velocity vector and transfers the accumulated "kinetic energy" to the $x$ component of the velocity pointing towards the bulk. In the regime of large $|\Omega|\tau\gg 1$, the tangential field acts on a shorter time scale than the scale of the dissipative Stokes force and the velocity is nearly unaffected by the friction. The particle leaves the wall with a velocity proportional to $\Omega$ and enters the potential-free region where it performs a spiral motion with a large initial radius, due to the energy accumulated at the wall. This radius continuously shrinks due to the dissipation caused by the Stokes bulk force, as shown in Figure 4. After a time $\sim\tau$ particles the dissipation becomes relevant and $\langle|v|\rangle\sim\sqrt{D_{a}/\tau}$, i.e. the steady state typical velocity. The process here described goes on forever since a fluctuation of the self-propulsion can drive a bulk particle to reach again the wall region to start a new cycle. It is clear, that the removal of particles from the walls due to $\Omega$ has the effect of depleting the wall accumulation and increasing the bulk population.

To test the existence of large surface currents leading to the formation of an intermediate emptying region we have studied the profile of the vertical velocity. In the two panels Fig. 5 we display the average velocity $\langle v_{y}\rangle_{x}$ as a function of the distance from the center of the slit $x$ for two different values of $\Omega$. We notice the monotonic behavior of the velocity, which can be described by a hyperbolic sine function since the average velocity in the left-hand side has the opposite direction with respect to velocity in the right-hand side. As discussed above, the mechanism which causes such a behavior is the interplay between the wall force and the chiral force, however, in this case, the non-uniformity of the Stokes force in the $x$-direction leads to a richer structure both in the density and in the velocity profiles than the one observed in the parabolic channel. In fact, the Stokes force opposing the motion along the $x$ direction is large only in the wall region $|x|>L$, while the centripetal force $\Omega\hat{\mathbf{z}}\times\mathbf{v}$ is spatially uniform.

4.3 Momentum profile in the potential-free region $|x|<L$ and effective viscosity

As shown in the two panels of Fig. 5, the vertical component of the average velocity field varies in the region bounded by the walls. Hereafter, we derive the equations for the current by projecting the Fokker-Planck equation onto an appropriate hydrodynamic space. By integrating the FPE (10) over $v_{x},v_{y}$ we obtain the continuity equation:

[TABLE]

with

[TABLE]

In the steady state $\frac{\partial j_{x}(x,t)}{\partial x}=0$ and by symmetry the particle current ${\bf j}$ is a function of $x$ only. From the continuity equation and the wall boundary conditions it follows that $j_{x}(x)=0$, while the vertical component $j_{y}(x)$ varies with $x$. Next, we define the pressure tensor ${\mathcal{P}}$ as:

[TABLE]

multiply the stationary FPE by $v_{y}$ and $v_{x}$, respectively, and integrate over the velocity. In the potential-free region, $|x|<L$, we obtain the following equations relating the components of the pressure tensor to the current $j_{y}$ :

[TABLE]

[TABLE]

In order to determine $j_{y}(x)$ from Eq. (32), we must relate it to ${\mathcal{P}}_{xy}(x)$, a task which will be pursued hereafter using a simple kinetic argument.

In the situation depicted in Fig. 5, let us consider the average vertical momentum flux $J_{y}^{+}$ per unit time and per unit length crossing a segment at $x=x_{0}$ (with $-L<x_{0}<L$) and originating from the region $x<x_{0}$. We assume that: a) the particles passing from the one side to the other of the unit vertical segment at $x_{0}$ move at the constant average velocity, $|\overline{v}_{x}|$ b) the physical space can be divided into cells of linear size $\lambda_{\Omega}$, the smallest distance over which the mean values of the physical observables vary. Under these hypotheses, the magnitude of $J_{y}^{+}$ can be written as:

[TABLE]

Similarly, the momentum flux of the particles coming from the right is

[TABLE]

where the geometrical factor $1/4$ takes into account the different possible directions of the velocities of the particles under the assumption of isotropy of the velocity distribution. The result is a net transfer of the $y$-component of momentum across the segment $x=x_{0}$ given by:

[TABLE]

The rate of change of momentum per unit length, $(J_{y}^{+}-J_{y}^{-})$, due to the diffusion of the particles in the horizontal direction is equivalent to the shear force per unit length exerted by the particles with $x<x_{0}$ on the particles with $x>x_{0}$, that is to the off-diagonal component of the pressure tensor:

[TABLE]

with $\nu_{\Omega}=\frac{1}{2}\lambda_{\Omega}|\overline{v}_{x}|$. The same kind of elementary kinetic argument, this time replacing the transported quantity $j_{y}(x)$ by $\rho(x)$, predicts that a gas of tracer chiral active particles of non uniform density $\rho(x)$ generates a net flux $J_{N}$, defined as the number of particles per unit time crossing a unit vertical segment according to the formula:

[TABLE]

The ratio between $-J_{N}$ and $-\frac{\partial\rho(x)}{\partial x}$ defines the the self-diffusion coefficient of chiral particles $D_{\Omega}=\frac{1}{2}\lambda_{\Omega}|\overline{v}_{x}|$. One concludes that the kinematic shear viscosity and the diffusion coefficient are equal: $\nu_{\Omega}=D_{\Omega}$, a result well known in kinetic theory of dilute gases. By comparing this result with the exact calculation of the long-time diffusion coefficient $D_{\Omega}=\frac{D_{a}}{1+\Omega^{2}\tau^{2}}$ given by Eq. (17) we can now directly determine the value of $\nu_{\Omega}$ and thus obtain

[TABLE]

Finally, substituting in Eq. (32) we obtain

[TABLE]

whose solution is

[TABLE]

where the constant $A$ is fixed by the odd boundary conditions at $x=\pm L$. Using Eq. (33) we can obtain ${\mathcal{P}}_{xx}$:

[TABLE]

where the constant ${\mathcal{P}}_{0}$ is determined by the boundary conditions. Let us remark that $j_{y}(x)$ is proportional to $\langle v_{y}\rangle_{x}$ in the regime where the density $\rho(x)$ varies slower than $\langle v_{y}\rangle_{x}$ itself. Under this condition it is legitimate to compare directly the solution $j_{y}(x)$ of Eq. (41) with the profile $\langle v_{y}\rangle_{x}$ of Fig. 5. The present estimate of the shear kinematic viscosity is of the right order and we find a quantitative agreement between the prediction of of Eq. (41) and the numerical simulation by a simple rescaling of the characteristic length: $\sqrt{\nu_{\Omega}\tau}\to\sqrt{2\nu_{\Omega}\tau}$.

5 Conclusions

In this paper, we have introduced an extension of the active Ornstein-Uhlenbeck model to study the behavior of a gas of active chiral particles under confinement. The handedness of particles’ motion has been accounted for by means of an effective torque term in the self-propulsion forcing. The presence of the torque has a profound influence on the properties of the AOUP as it breaks the detailed balance condition even in the unconfined case. Under confinement, we find that the torque is responsible for the appearance of steady momentum currents and of reduction of the accumulation of the particles at the container’s boundaries. Explicit illustrations of these phenomena are shown by means of the exact solution of the stationary Fokker-Planck equation in two cases of parabolic confinement. In a more realistic case of confinement by stiff walls, beside the emptying of the boundary region and the enrichment of the bulk region, we find a structure the velocity field akin to the one observed in sheared viscous fluids. This behavior is explained in terms of a simple kinetic argument and the velocity profile is reproduced.

Since the linear propulsive and rotational behavior can independently be tuned, by changing $\tau$ and $\Omega$ respectively, one could drive more efficiently the motion of active particles for instance by partly suppressing their diffusivity, or obtain more efficient harmonic traps. Active rotation does also have an impact on the properties of interacting systems: Liao and Klapp have found that it generally opposes motility-induced clustering and phase separation, as demonstrated by a narrowing of the coexistence region of two-dimensional chiral ABP upon increase of the propulsion angular velocity 44. A similar behaviour should also be observablein systems of interacting CAOUP.

Finally, concerning the practical interest of the systems studied, one could envisage the possibility of selecting chiral active particles on the basis of their handedness in pharmaceutical or biotechnological applications where chiral levogyre or dextrogyre properties correspond to different functionalities 19, 45.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

We thank Andrea Puglisi and Claudio Maggi for illuminating discussions.

Appendix A Transformation of variables

In this appendix, we derive the governing equations for the transformed variables $(x_{i},v_{i})$. Let us consider the AOUP equation:

[TABLE]

and define the "velocity", $v_{i}$: $v_{i}=\dot{x}_{i}-\sqrt{2D_{t}}\xi_{i}(t)$. We differentiate (43) and eliminate $u_{i}$ with the help of Eq. (4) and use again (43). Adopting Stratonovich convention for the derivatives of $\phi$ we finally get:

[TABLE]

where we used the Einstein repeated indexes convention.

Appendix B Exact distribution in the parabolic channel

Hereafter, we report the exact phase-space distribution in the case of a parabolic channel. Its form as discussed in the main text depends on six constants $A,K_{1},K_{2},K_{3},M_{1},N_{2}$ according to the formula:

[TABLE]

For the sake of notational simplicity we define:

[TABLE]

and write the following six equations for the coefficients in terms of the independent control parameters:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The average velocity along the y is

[TABLE]

and is proportional to $\frac{\Omega\tau}{\gamma}k$.

Notes and references

• Ramaswamy 2010

S. Ramaswamy, The Mechanics and Statistics of Active Matter, 2010, 1, 323–345.

• Bechinger et al. 2016

C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Reviews of Modern Physics, 2016, 88, 045006.

• Marchetti et al. 2013

M. Marchetti, J. Joanny, S. Ramaswamy, T. Liverpool, J. Prost, M. Rao and R. A. Simha, Reviews of Modern Physics, 2013, 85, 1143.

• Romanczuk et al. 2012

P. Romanczuk, M. Bär, W. Ebeling, B. Lindner and L. Schimansky-Geier, The European Physical Journal Special Topics, 2012, 202, 1–162.

• Lauga et al. 2006

E. Lauga, W. R. DiLuzio, G. M. Whitesides and H. A. Stone, Biophysical journal, 2006, 90, 400–412.

• Friedrich and Jülicher 2008

B. Friedrich and F. Jülicher, New Journal of Physics, 2008, 10, 123025.

• DiLuzio et al. 2005

W. R. DiLuzio, L. Turner, M. Mayer, P. Garstecki, D. B. Weibel, H. C. Berg and G. M. Whitesides, Nature, 2005, 435, 1271.

• Loose and Mitchison 2014

M. Loose and T. J. Mitchison, Nature cell biology, 2014, 16, 38.

• Kümmel et al. 2013

F. Kümmel, B. ten Hagen, R. Wittkowski, I. Buttinoni, R. Eichhorn, G. Volpe, H. Löwen and C. Bechinger, Physical review letters, 2013, 110, 198302.

• Löwen 2016

H. Löwen, The European Physical Journal Special Topics, 2016, 225, 2319–2331.

• Blakemore 1975

R. Blakemore, Science, 1975, 190, 377–379.

• Lefèvre and Bazylinski 2013

C. T. Lefèvre and D. A. Bazylinski, Microbiology and Molecular Biology Reviews, 2013, 77, 497–526.

• Ten Hagen et al. 2015

B. Ten Hagen, R. Wittkowski, D. Takagi, F. Kümmel, C. Bechinger and H. Löwen, Journal of Physics: Condensed Matter, 2015, 27, 194110.

• Najafi and Golestanian 2004

A. Najafi and R. Golestanian, Physical Review E, 2004, 69, 062901.

• van Teeffelen and Löwen 2008

S. van Teeffelen and H. Löwen, Physical Review E, 2008, 78, 020101.

• Ebeling et al. 1999

W. Ebeling, F. Schweitzer and B. Tilch, BioSystems, 1999, 49, 17–29.

• Cates and Tailleur 2013

M. Cates and J. Tailleur, EPL (Europhysics Letters), 2013, 101, 20010.

• Szamel 2014

G. Szamel, Physical Review E, 2014, 90, 012111.

• Mijalkov and Volpe 2013

M. Mijalkov and G. Volpe, Soft Matter, 2013, 9, 6376–6381.

• Volpe et al. 2014

G. Volpe, S. Gigan and G. Volpe, American Journal of Physics, 2014, 82, 659–664.

• Li et al. 2014

Y. Li, P. K. Ghosh, F. Marchesoni and B. Li, Physical Review E, 2014, 90, 062301.

• Ai et al. 2015

B.-q. Ai, Y.-f. He and W.-r. Zhong, Soft Matter, 2015, 11, 3852–3859.

• Ao et al. 2015

X. Ao, P. K. Ghosh, Y. Li, G. Schmid, P. Hänggi and F. Marchesoni, EPL (Europhysics Letters), 2015, 109, 10003.

• Marconi and Maggi 2015

U. M. B. Marconi and C. Maggi, Soft matter, 2015, 11, 8768–8781.

• Das et al. 2018

S. Das, G. Gompper and R. G. Winkler, New Journal of Physics, 2018, 20, 015001.

• Caprini et al. 2018

L. Caprini, U. M. B. Marconi and A. Puglisi, arXiv preprint arXiv:1810.12652, 2018.

• Fily and Marchetti 2012

Y. Fily and M. C. Marchetti, Physical Review Letters, 2012, 108, 235702.

• Farage et al. 2015

T. F. Farage, P. Krinninger and J. M. Brader, Physical Review E, 2015, 91, 042310.

• Marconi et al. 2016

U. M. B. Marconi, N. Gnan, M. Paoluzzi, C. Maggi and R. Di Leonardo, Scientific reports, 2016, 6, 23297.

• Caprini and Marconi 2018

L. Caprini and U. M. B. Marconi, Soft matter, 2018.

• Howse et al. 2007

J. R. Howse, R. A. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh and R. Golestanian, Physical review letters, 2007, 99, 048102.

• Risken 1984

H. Risken, Fokker-Planck Equation, Springer, 1984.

• Puglisi and Marini Bettolo Marconi 2017

A. Puglisi and U. Marini Bettolo Marconi, Entropy, 2017, 19, 356.

• Cates 2012

M. Cates, Reports on Progress in Physics, 2012, 75, 042601.

• Marconi et al. 2017

U. M. B. Marconi, A. Puglisi and C. Maggi, Scientific reports, 2017, 7, 46496.

• ten Hagen et al. 2011

B. ten Hagen, S. van Teeffelen and H. Löwen, Journal of Physics: Condensed Matter, 2011, 23, 194119.

• Dauchot and Démery 2018

O. Dauchot and V. Démery, arXiv preprint arXiv:1810.13303, 2018.

• Hanggi and Jung 1995

P. Hanggi and P. Jung, Advances in Chemical Physics, 1995, 89, 239–326.

• Marconi et al. 2016

U. M. B. Marconi, C. Maggi and S. Melchionna, Soft matter, 2016, 12, 5727–5738.

• Marini Bettolo Marconi et al. 2017

U. Marini Bettolo Marconi, C. Maggi and M. Paoluzzi, The Journal of chemical physics, 2017, 147, 024903.

• Wensink and Löwen 2008

H. Wensink and H. Löwen, Physical Review E, 2008, 78, 031409.

• Elgeti and Gompper 2013

J. Elgeti and G. Gompper, EPL (Europhysics Letters), 2013, 101, 48003.

• Lee 2013

C. F. Lee, New Journal of Physics, 2013, 15, 055007.

• Liao and Klapp 2018

G.-J. Liao and S. H. Klapp, Soft matter, 2018, 14, 7873–7882.

• Meinhardt et al. 2012

S. Meinhardt, J. Smiatek, R. Eichhorn and F. Schmid, Physical review letters, 2012, 108, 214504.