Collective Surfing of Chemically Active Particles
Hassan Masoud
1,2
,*
and Michael J. Shelley
1
,†
1
Applied Mathematics Laboratory, Courant Institute of Mathematical Sciences, New York University, New York, New York 10012, USA
2
Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, New Jersey 08544, USA
(Received 17 December 2013; published 26 March 2014)
We study theoretically the collective dynamics of immotile particles bound to a 2D surface atop a 3D
fluid layer. These particles are chemically active and produce a chemical concentration field that creates
surface-tension gradients along the surface. The resultant Marangoni stresses create flows that carry the
particles, possibly concentrating them. For a 3D diffusion-dominated concentration field and Stokesian
fluid we show that the surface dynamics of active particle density can be determined using nonlocal 2D
surface operators. Remarkably, we also show that for both deep or shallow fluid layers this surface
dynamics reduces to the 2D Keller-Segel model for the collective chemotactic aggregation of slime mold
colonies. Mathematical anal ysis has established that the Keller-Segel model can yield finite-time, finite-
mass concentration singularities. We show that such singular behavior occurs in our finite-depth system,
and study the associated 3D flow structures.
DOI:
10.1103/PhysRevLett.112.128304 PACS numbers: 47.57.J−, 47.20.Ma, 68.03.Cd
Active fluids or suspensions have attracted much atten-
tion for their interesting, often unexpected dynamics [1].
Examples of active fluids include suspensions of micro-
swimmers such as bacteria
[2–4], or of chemically or
optically driven particles [5–10], and complex networks of
biopolymers and molecular motors
[11,12]. These active
soft materials differ from passive ones as they continuously
consume energy from the surrounding environment to do
work and are far from thermodynamic equilibrium, even in
steady state conditions. This energy may come from
chemical energy and the work done on the fluid can be
used for self-propulsion. Most studies of active fluids have
focused on suspensions of motile particles
[1] whose active
stresses, produced by swimming, and spontaneous flows
can enhance mixing
[13,14], and affect chemotactic aggre-
gation (see Ref.
[15] and references therein).
Here, we investigate a new type of activ e fluid with a
different source of active stress. Consider immotile (non-
swimming), but chemically active, particles that are bound to a
flat surface sitting above a fluidlayer
[16–20]. Theseparticles’
activity creates (or depletes) a spatially diffusing chemical
concentration field. At the surface, this chemical changes the
local surface tension, and any consequent gradients in surface
tension will produce Marangoni shear stresses. These “acti v e”
Marangoni stresses will produce fluid flows that move the
particles. We call this surfing
[21], and study it in the case of
negligible fluid inertia and a chemical species whose transport
is dominated by diffusion. Our theoretical analyses and
simulations suggest that, if particle concentration raises sur-
face tension (throughthe induced concentrationfield), surface
flows of chemical surfers can yield large aggre gations
associated with vortical flows in the bulk. These might be
harnessed for microfluidic manipulations
[22,23] or flow
assisted self-assembly [24,25].
Very surprisingly, we also show that for both sufficiently
deep or shallow fluid layers the surface dynamics of
particle density recovers the 2D parabolic-elliptic Keller-
Segel (KS) model. Originally conceived to describe the
chemotactically driven aggregation of motile slime molds,
the KS model
[26] is a canonical model of mathematical
biology
[27,28]. It describes the collective chemotactic
dynamics of motile organisms that secrete and respond to a
diffusing chemoattractant. KS dynamics can lead to aggre-
gation and, under easily met analytical conditions, chemo-
tactic collapse, which is an infinite pointwise density in
finite time
[27,28]. In the fluid dynamical context this
collapse is associated with divergent velocity gradients and
intense vortical structures.
That the dynamics of our system is, in certain limits,
identical to the KS system is not obvious. However, the result
derives from a very direct Fourier transform argument. This
reveals an unexpected connection between self-attractiv e
chemotaxis (termed “autotaxis”)oforganismsandapurely
physical fluid system driven by Marangoni stresses. In the
following, we first describe how the dynamics of active
particles is modeled followed by the derivation of the fluidic
analogy to the KS model, and give results.
Consider a flat free surface, at z ¼ 0, that sits above a
layer of Newtonian fluid, of viscosity μ, bounded by an
impermeable solid wall at z ¼ −H [see Fig.
1(a)]. Let
ψðx
;tÞbe the number density field of active particles on the
free surface. The particles are carried passively by the
surface fluid velocity Uðx;tÞ and diffuse along the surface
with diffusion constant D
p
. Then, ψ satisfies
ψ
t
þ ∇
2
· ðUψ Þ¼D
p
Δ
2
ψ (1)
in the z ¼ 0 plane where the subscript 2 denotes spatial
derivatives in the x ¼ðx; yÞ plane. These particles are
PRL 112, 128304 (2014)
PHYSICAL REVIEW LETTERS
week ending
28 MARCH 2014
0031-9007=14=112(12)=128304(5) 128304-1 © 2014 American Physical Society
chemically active and either deplete or produce a chemical
species that diffuses into the bulk. Thus, the concentration
field Cðx;z;tÞ of the chemical species satisfies
C
t
þ u · ∇C ¼ D
c
ΔC; (2)
where D
c
is a diffusion constant, Δ is the 3D Laplacian in x
and z, and uðx;z;tÞ is the 3D fluid velocity. The flux
boundary conditions of Eq.
(2) are
D
c
C
z
ðx; 0;tÞ¼
_
mψðx;tÞ and C
z
ðx; −H; tÞ¼0; (3)
where
_
m is the rate of production (
_
m>0) or consumption
(
_
m<0) of the chemical species per active particle.
The inertialess incompressible fluid flow with velocity
uðx;z;tÞ¼ðu; v; wÞ, pressure p, and stress tensor σ ¼
−pI þ μ½∇u þð∇uÞ
T
is driven by Marangoni shear
stresses produced by surface gradients in C. We make
the standard assumption that the surface tension γ depends
linearly on the surface concentration of the chemical
species
[29,30], so that γðx;tÞ¼γ
0
þ αCðx; 0;tÞ, where
γ
0
and α are constants. To determine the surface velocity
U ¼ðu; vÞj
z¼0
, we must solve the 3D Stokes equations
∇ · σ ¼ −∇p þ μΔu ¼ 0 and ∇ · u ¼ 0 (4)
with the boundary conditions
μðu
z
;v
z
Þj
z¼0
¼ ∇
2
γ ¼ α∇
2
C; wðx; 0;tÞ¼0; (5)
and no slip at z ¼ −H.
We consider the coupled Eqs.
(1)–(5) in a L × L × H
domain V that is L periodic in the x and y directions. We
make these equations dimensionless by scaling particle
density, length, chemical concentration, time, and velocity
with, respectively, the average conserved particle density
¯
ψ, L,
_
m
¯
ψ L=D
c
, advection time scale τ ¼ μD
c
=
_
mα
¯
ψ,
and L=τ.
Assuming Cðx;z;0Þ¼0, the integration of Eq.
(2) over
the volume and applying the boundary conditions, Eqs. (3) ,
yields
¯
CðtÞ¼ð
R
CdVÞ=
¯
V ¼ t=ðδPe
c
Þ, where
¯
V is the
domain volume, δ ¼ H=L, and Pe
c
¼ τ
c
=τ with τ
c
¼
L
2
=D
c
being the time scale of diffusion of the chemical
species. Letting Cðx;z;tÞ¼ϕðx;z;tÞþ
¯
CðtÞ, Eq.
(2)
becomes
Pe
c
ðϕ
t
þ u · ∇ϕÞ¼Δϕ − δ
−1
for − δ ≤ z ≤ 0: (6)
Assuming Pe
c
≪ 1 and so neglecting the left-hand side, ϕ
satisfies a quasisteady Poisson equation (after an initial
transient), which can be solved via 2D Fourier transform in
x (see the Supplemental Material
[31]). This yields
ˆ
ϕð0;zÞ¼z
2
=ð2δÞþz þ δ =3 and
ˆ
ϕðk;z;tÞ¼ðcothkδ cosh kz þ sinh kzÞ
ˆ
ψðk;tÞ=k (7)
for k ≠ 0, where k ¼ð2πn
1
; 2πn
2
Þ is the 2D wave vector
with k ¼jkj and n
1
, n
2
are integers. The surface gradient
d
∇
2
ϕj
z¼0
¼ðik=kÞcothðkδÞ
ˆ
ψðk;tÞ¼
c
R
δ
½ψ is a general-
ized Riesz transform
[32].
Though slightly more complicated, Eq.
(4) with boun-
dary conditions (5) can again be solved via Fourier trans-
form in x (see the Supplemental Material
[31]), which
allows us to relate the surface velocity to the density of
active particles,
ˆ
Uðk;tÞ¼ðik=k
2
ÞΩðkδÞ
ˆ
ψðk;tÞ (8)
for k ≠ 0, where
ΩðλÞ¼cothλðsinh
2
λ − λ
2
Þ=ðsinh 2λ − 2λÞ: (9)
Note that
ˆ
Uðk ¼ 0;tÞ¼0. Ω has a simple structure: for λ
small, Ω ≃ 1=4 þ Oðλ
2
Þ whereas for λ large, Ω approaches
1=2 exponentially fast [Fig.
1(b)]. That Ω remains finite as
δ → 0 is interesting. While viscous dissipation increases
with narrowing of the gap between the free surface and
solid wall, the surface gradient of C, and therefore the
driving Marangoni stress, increases with a decrease in δ
[see Eq.
(7)]. These two effects cancel, which gives rise to a
finite Ω even when δ is very small.
Hence, in both limits, Eq.
(8) reduces to
β
ˆ
Uðk;tÞ¼ðik=k
2
Þ
ˆ
ψðk;tÞ; (10)
where β ¼ 2 (deep layer) or 4 (shallow layer). In real space
we can write Eq.
(10) as
βU ¼ −Δ
−1
2
∇
2
ψ; (11)
λ
λ
z
(a)
(b)
L
y
FIG. 1 (color online). (a) Schematic illustrating chemically
active particles (dark circles) bound to a flat fluid surface sitting
above a fluid layer of depth H. The gray-scale map represents the
concentration field of the chemical species produced by the active
particles. (b) Variation of Ω as a function of λ.
PRL 112, 128304 (2014)
PHYSICAL REVIEW LETTERS
week ending
28 MARCH 2014
128304-2
and Eq. (1) becomes
ψ
t
− ∇
2
· ½ðψ =βÞðΔ
−1
2
∇
2
ψÞ ¼ Pe
−1
p
Δ
2
ψ; (12)
where Pe
p
¼ τ
p
=τ with τ
p
¼ L
2
=D
p
the time scale of
diffusion of active particles. Surprisingly, rescaling of
Eq.
(12) when
_
mα > 0 recovers the 2D parabolic-elliptic
KS model for autotactic aggregation. This model is given
by equations for the organismal density φðx;tÞ and the
collectively produced chemoattractant concentration
ρðx;tÞ,
φ
t
þ ∇
2
· ðχφ∇
2
ρÞ¼Δ
2
φ and Δ
2
ρ ¼ −φ: (13)
Read as a kinetic equation for conservation of species
number, the first equation of Eqs.
(13) states that the
organismal speed is χ∇
2
ρ with “chemotactic” strength χ.
The Poisson equation states that the chemoattractant is
produced locally at a rate proportional to the organismal
density, and is rapidly diffused. Its solution can be written
formally as ρ ¼ −Δ
−1
2
φ. Thus, the evolution equation for φ
has the same form as Eq.
(12). However, though ψ and φ
are governed by the same equation, surface distributions of
C and ρ are different since Δ
2
C ¼ ∇
2
· R
δ
½ψ ≠−ψ. Much
is known about Eqs.
(13). For instance, given a sufficient
mass of organisms in the plane, the 2D KS model leads to
chemotactic collapse in finite time
[27,28,33,34]. The
collapse singularity is approximately self-similar, i.e.,
φðx;tÞ ≈ ζ
−2
Φðjxj=ζÞ for some function Φ and a scale ζ
whose dominant algebraic part is
ffiffiffiffiffiffiffiffiffiffiffi
t
c
− t
p
, with t
c
the
collapse time
[35]. In addition, ζ → 0 and φ becomes a
Dirac δ function in x as t → t
c
. See Refs.
[27,28] for
comprehensive reviews.
Consistent with this, linear stability analysis of Eqs. (1)
and (8) shows that the system is unstable to 2D surface
flows if Pe
p
> 2=Ωð
ffiffiffi
2
p
δÞ as Pe
p
scales linearly with
¯
ψ
(see the Supplemental Material
[31]). In this scenario,
particle activity locally increases the surface tension
(α
_
m>0) and Marangoni stresses produce flows that
concentrate the surface density of particles, leading to
yet higher surface tension. Unlike the 2D KS problem, the
surface flow in our system is associated with fully 3D fluid
flow and structures.
We simulate Eqs.
(1) and (8) using a semi-implicit,
second-order in time, Fourier pseudospectral method for
_
mα > 0. We set ψðx; 0Þ¼1 þ 0.1 cos ð2πxÞcos ð2πyÞ and
fix Pe
p
¼ 400 which gives rise to instability independent
of δ. We observe a rapid accumulation of active surface-
bound particles at the center and corners of the domain,
where the initial concentration is peaked. Figure
2(a) shows
a snapshot of ψðx;tÞ near the collapse time when
ψ
max
¼ ψð1=2; 1 =2Þ ≈ 10. Due to the similarity of patterns
across δ, we only show the results for the case of δ ¼ 1=2.
We see an analogous, though broader, distribution for C on
the surface [Fig.
2(b)]. Note that C is one derivative
smoother than ψ [see Eq. (7)], and on the surface blows
up more rapidly than ρ in the KS model as ψ collapses.
The increased active particle density is accompanied by
fluid flow towards the blowup points [Fig. 2(c)]. This
surface flow generates a 3D flow in the bulk. Figure
2(d)
shows the in-plane vorticity field at z ¼ 0 overlaid with the
contour map of its magnitude, which highlights the 3D
nature of the velocity field. The out-of-plane vorticity is
zero since U is a 2D gradient [Eq.
(8)]. To better understand
the bulk flow, consider the δ → ∞ case (β ¼ 2). Applying
the incompressibility condition at z ¼ 0 gives [see Eq.
(11)]
∇
2
· U ¼ −w
z
j
z¼0
¼ð1 − ψÞ=2: (14)
This suggests that the out-of-plane rate of strain directly
inherits the structure of ψ, and so diverges with any
collapse. Figure
3(a) demonstrates such a behavior where
the surface flow creates a vortex ring just below the free
surface. Though with a different strength and dimension, a
0 0.25 0.5 0.75 1
0
0.25
0.5
0.75
1
(a)
(b)
(c)
(d)
0
10
170.4
171
0
3.5
FIG. 2 (color online). For δ ¼ 1=2, contours of (a) particle density ψ, (b) chemical surface concentration C, and (c) the surface velocity
field from Eq. (8). (d) The in-plane surface vorticity field ω overlaid with a color map of vorticity magnitude ω. The data correspond to
t=t
c
≈ 0.9 where t
c
is the estimated collapse time. Simulation parameters are Δt ¼ 1=400 , Δx ¼ Δy ¼ 1=256, Pe
c
¼ 0.1, Pe
p
¼ 400,
and ψ ðx; 0Þ¼1 þ 0.1 cos ð2πxÞcos ð2πyÞ. Vorticity is scaled by 1=τ.
PRL 112, 128304 (2014)
PHYSICAL REVIEW LETTERS
week ending
28 MARCH 2014
128304-3
vortex ring emerges in shallower fluid layers as well [see
Figs.
3(b) and 3(c)].
Finally, we examine the effect of fluid layer depth on the
putative collapse time. Figure
4 shows the divergence of
ψ
max
for different values of δ. We see that t
c
increases with
decreasing δ, but remains finite for both deep or shallow
layers. The collapse time is very well described by the
shallow depth KS model for δ ≲ 0.01 and by the infinite
depth KS model for δ ≳ 1. While the form of singularity in
the asymptotic limits of δ appears to be consistent with KS
collapse, the precise nature of the apparent singularity for
finite δ is beyond the scope of this Letter.
Can our results be realized experimentally? The values
for Pe
p
and Pe
c
used in our simulations can be realized in
an aqueous system of size L ∼ 100 μm and characteristic
fluid velocity L=τ ∼ 1 μm =s, using chemically active par-
ticles of radius R ∼ 100 nm that produce a chemical species
with a diffusion constant D
c
∼ 10
−9
m
2
=s. Surface tension
in this system is strong enough to keep the free surface flat.
A challenge in developing such a system is designing a
chemical reaction whose product increases the surface
tension or whose reactants decrease the surface tension.
Alternatively, an endothermic chemical reaction could
increase surface tension by lowering the local temperature.
Also, our results near the collapse time are not expected to
be realized quantitatively since collapse leads to violation
of many of our modeling assumptions such as finite
velocity gradients, linearity of surface tension dependence
on C, and large dynamic range of the surface tension
(which, in reality, is small).
In this work, we considered active particles of isotropic
shape that are much smaller than the characteristic size of
the system. The use of anisotropic, or larger, active particles
[16] could lead to other types of instabilities with different
emerging flow patterns. Note that Marangoni stresses and
associated 3D fluid flow in our system could be used for
microfluidic manipulations
[22] and directed self-assembly
[24,25]. In fact, Marangoni stresses, induced by a chemical
reaction, were utilized in a related system to create self-
propelling liquid microdroplets
[36]. It has been also
reported that bacteria in biofilms exploit Marangoni
stresses for dispersal
[21], so it would be interesting to
see whether any surface bound organisms take advantage of
the Marangoni stresses for aggregation (say by consuming
surfactant).
Lastly, we showed analytically that for sufficiently deep
or shallow fluid layers the collective surfing of active
particles is described by the iconic Keller-Segel model and
we used the existing knowledge about its behavior to
enhance our understanding of singularities in flow.
(b)
7362.2
7362.35
2362
.3
2363
(c)
(a)
170.3
171
0
0.25
0.5
0.75
1
0
0.25
0.5
0.75
1
0
0.25
0.5
0.75
1
-
-
-
-
FIG. 3 (color online). The 3D velocity field overlaid with a color map of the chemical concentration C, at the time when
ψ
max
¼ ψ ðL=2;L=2Þ ≈ 10. δ ¼ 1=2, 1=10, 1=20 in (a), (b), and (c), respectively. Simulation parameters are the same as for Fig.
2. The z
axis is scaled by δ. The black arcs in (a) highlight vortical flows in the bulk.
1/
ψ
max
t
5
20
35
50
0.0001 0.001 0.01 0.1 1 10
t
c
δ
Ω = 1/2 (KS)
δ
= 1/2
δ
= 1/10
δ
= 1/20
δ
= 1/40
Ω = 1/4 (KS)
0
0.25
0.5
0.75
1
1.25
1.5
0 1020304050
FIG. 4 (color online). Time variation of ψ
max
¼ ψð1=2; 1=2Þ
for different values of δ and asymptotic cases of Ω ¼ 1=2, 1=4 for
which the evolution of ψ obeys the 2D KS model. The inset
shows the variation of the estimated collapse time t
c
as a function
of δ. Crosses show the asymptotic values of t
c
for Ω ¼ 1=2, 1=4.
Simulation parameters are the same as for Fig.
2.
PRL 112, 128304 (2014)
PHYSICAL REVIEW LETTERS
week ending
28 MARCH 2014
128304-4
Alternatively, it might prove fruitful to explore aspects of
biological systems described by the KS model and its
variants using chemically active particles that create active
surface stresses for motion.
We thank J. Bedrossian, J. W. M. Bush, S. Childress, E.
Nazockdast, L. Ristroph, and J. Zhang for useful con-
versations. H. M. and M. J. S. were supported by DOE
Grant No. DE-FG02-88ER25053. H. M. acknowledges
support from NSF Grant No. DMR-0844115 and the
Institute for Complex Adaptive Matter. M. J. S. acknowl-
edges support from NSF (NYU MRSEC DMR-0820341).
*
Hassan.Masoud@courant.nyu.edu
†
shelley@cims.nyu.edu
[1] D. Saintillan and M. J. Shelley, C.R. Phys. 14, 497 (2013).
[2] A. Sokolov, I. S. Aranson, J. O. Kessler, and R. E. Goldstein,
Phys. Rev. Lett. 98, 158102 (2007).
[3] H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler, and
R. E. Goldstein,
Phys. Rev. Lett. 110, 268102 (2013) .
[4] L. H. Cisneros, J. O. Kessler, S. Ganguly, and R. E. Goldstein,
Phys. Rev. E 83, 061907 (2011).
[5] W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. St
Angelo, Y. Y. Cao, T. E. Mallouk, P. E. Lammert, and
V. H. Crespi, J. Am. Chem. Soc. 126, 13424 (2004).
[6] Y. Hong, N. M. K. Blackman, N. D. Kopp, A. Sen, and
D. Velegol,
Phys. Rev. Lett. 99, 178103 (2007).
[7] D. Takagi, A. B. Braunschweig, J. Zhang, and M. J. Shelley,
Phys. Rev. Lett. 110, 038301 (2013).
[8] J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine, and
P. M. Chaikin,
Science 339, 936 (2013).
[9] R. Golestanian, Phys. Rev. Lett. 108, 038303 (2012).
[10] T. Bickel, G. Zecua, and A. Würger, arXiv:1401.7833.
[11] T. Shinar, M. Mana, F. Piano, and M. J. Shelley, Proc. Natl.
Acad. Sci. U.S.A. 108, 10508 (2011)
.
[12] M. Sheinman, C. P. Broedersz, and F. C. MacKintosh,
Phys. Rev. Lett. 109, 238101 (2012).
[13] D. Saintillan and M. J. Shelley, Phys. Rev. Lett. 100, 178103
(2008)
; Phys. Fluids 20, 123304 (2008).
[14] A. Sokolov, R. E. Goldstein, F. I. Feldchtein, and I. S.
Aranson, Phys. Rev. E 80, 031903 (2009).
[15] E. Lushi, R. E. Goldstein, and M. J. Shelley, Phys. Rev. E
86, 040902(R) (2012)
.
[16] H. Zhang, W. Duan, L. Liu, and A. Sen,
J. Am. Chem. Soc.
135, 15734 (2013).
[17] E. Lauga and A. M. J. Davis, J. Fluid Mech. 705, 120
(2012)
.
[18] H. Masoud and H. A. Stone, J. Fluid Mech. 741, R4 (2014).
[19] H. A. Stone and A. Ajdari, J. Fluid Mech.
369, 151
(1998).
[20] H. A. Stone,
Phys. Fluid s A 2, 111 (1990).
[21] T. E. Angelini, M. Roper, R. Kolter, D. A. Weitz, and
M. P. Brenner,
Proc. Natl. Acad. Sci. U.S.A. 106 , 18109
(2009).
[22] H. A. Stone, A. D. Stroock, and A. Ajdari, Annu. Rev. Fluid
Mech. 36, 381 (2004).
[23] H. Masoud and A. Alexeev,
Chem. Commun. (Cambridge)
1 (2011) 472.
[24] M. T. F. Rodrigues, P. M. Ajayan, and G. G. Silva, J. Phys.
Chem. B 117, 6524 (2013)
.
[25] B. A. Grzybowski, H. A. Stone, and G. M. Whitesides,
Nature (London) 405, 1033 (2000).
[26] E. F. Keller and L. A. Segel,
J. Theor. Biol. 26, 399 (1970);
30, 225 (1971); 30, 235 (1971).
[27] D. Horstmann, Jahresber. Dtsch. Math.-Ver. 105, 103
(2003).
[28] T. Hillen and K. J. Painter,
J. Math. Biol. 58, 183 (2009).
[29] W. E. Acree, J. Colloid Interface Sci. 101, 575 (1984).
[30] W. F. Paxton, A. Sen, and T. E. Mallouk, Chem. Eur. J. 11,
6462 (2005).
[31] See Supplemental Material at
http://link.aps.org/
supplemental/10.1103/PhysRevLett.112.128304
for a de-
tailed derivation of Eqs. (6–8) and the origin of the stability
condition.
[32] G. F. Carrier, M. Krook, and C. E. Pearson, Functions of a
Complex Variable: Theory and Technique (McGraw-Hill,
New York, 1966).
[33] S. Childress and J. K. Percus,
Math. Biosci. 56, 217 (1981).
[34] S. Childress, Chemotactic Collapse in Two Dimensions
(Springer-Verlag, Berlin, 1984), pp. 61–66.
[35] M. A. Herrero and J. J. L. Velazquez,
J. Math. Biol. 35, 177
(1996).
[36] S. Thutupalli, R. Seemann, and S. Herminghaus, New J.
Phys. 13, 073021 (2011).
PRL 112, 128304 (2014)
PHYSICAL REVIEW LETTERS
week ending
28 MARCH 2014
128304-5