# Generic non-uniqueness of minimizing harmonic maps from a ball to a sphere

Antoine Detaille  Universite Claude Bernard Lyon 1, ICJ UMR5208, CNRS, Ecole
Centrale de Lyon, INSA Lyon, Université Jean Monnet, 69622 Villeurbanne,
France<EMAIL_ADDRESS>and  Katarzyna Mazowiecka  Institute
of Mathematics,University of Warsaw, Banacha 2, 02-097 Warszawa, Poland
<EMAIL_ADDRESS>
###### Abstract.

In this note, we study non-uniqueness for minimizing harmonic maps from
$B^{3}$ to ${\mathbb{S}}^{2}$. We show that every boundary map can be modified
to a boundary map that admits multiple minimizers of the Dirichlet energy by a
small $W^{1,p}$-change for $p<2$. This strengthens a remark by the second-
named author and Strzelecki. The main novel ingredient is a homotopy
construction, which is the answer to an easier variant of a challenging
question regarding the existence of a norm control for homotopies between
$W^{1,p}$ maps.

###### Key words and phrases:

Harmonic maps, homotopy theory

###### 2010 Mathematics Subject Classification:

58E20, 46E35

###### Contents

  1. 1 Introduction
  2. 2 Homotopy construction
  3. 3 Proof of the generic non-uniqueness

##  1\. Introduction

Minimizing harmonic maps from $B^{3}$ to ${\mathbb{S}}^{2}$ are defined as
mappings with the least Dirichlet energy

(1.1) $E(u)\coloneqq\int_{B^{3}}|\nabla u|^{2}\,\mathrm{d}x$

among maps $u\in W^{1,2}(B^{3},{\mathbb{S}}^{2})$ with fixed boundary datum
$u\big{\rvert}_{\partial B^{3}}=\varphi\in W^{\frac{1}{2},2}(\partial
B^{3},{\mathbb{S}}^{2})$. Here, we minimize in the class of Sobolev maps with
values in a manifold (in our case, a sphere); for $s>0$ and $p\geq 1$, this
space is defined as

$W^{s,p}(\mathcal{M},{\mathcal{N}})\coloneqq\\{v\in
W^{s,p}(\mathcal{M},\mathbb{R}^{L})\colon v(x)\in{\mathcal{N}}\text{ for a.e.
}x\in\mathcal{M}\\}\text{,}$

where ${\mathcal{N}}\subset\mathbb{R}^{L}$ is a Riemannian manifold embedded
into $\mathbb{R}^{L}$ (in our case, ${\mathcal{N}}={\mathbb{S}}^{2}$) and
$\mathcal{M}$ is a compact Riemannian manifold (in our case,
$\mathcal{M}=B^{3}$ or $\mathcal{M}={\mathbb{S}}^{2}$).

The space $W^{1,2}(B^{3},{\mathbb{S}}^{2})$ is not a linear space, but it is
nevertheless a complete metric space endowed with the metric defined by

${\rm dist\,}{(u,v)}=\left\|{u-v}\right\|_{W^{1,2}(B^{3})}\text{.}$

We emphasize that, although being a subset of it, the class
$W^{1,2}(B^{3},{\mathbb{S}}^{2})$ exhibits some striking qualitative
differences with the linear space $W^{1,2}(B^{3},\mathbb{R}^{3})$. For
example, not every mapping $u\in W^{1,2}(B^{3},{\mathbb{S}}^{2})$ can be
approximated by smooth maps $u_{i}\in C^{\infty}(B^{3},{\mathbb{S}}^{2})$ in
the strong topology of $W^{1,2}$; see [14, Section 4]. However, maps
$\varphi\in W^{1,2}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ can be approximated in
$W^{1,2}$ by smooth maps $\varphi_{i}\in
C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$; see [13, Section 3].

For $\varphi\in W^{\frac{1}{2},2}(\partial B^{3},{\mathbb{S}}^{2})$, we also
define the space

$W^{1,2}_{\varphi}(B^{3},{\mathbb{S}}^{2})\coloneqq\\{v\in
W^{1,2}(B^{3},{\mathbb{S}}^{2})\colon v=\varphi\text{ on }\partial B^{3}\text{
in the trace sense}\\}$

and note that this space is always nonempty. For instance, for a given smooth
boundary datum $\varphi\in C^{\infty}(\partial B^{3},{\mathbb{S}}^{2})$, one
can easily construct an extension $u\in W^{1,2}(B^{3},{\mathbb{S}}^{2})$ of
$\varphi$, simply by considering $u(x)=\varphi\bigl{(}\frac{x}{|x|}\bigr{)}$.
More generally, any boundary map $\varphi\in W^{\frac{1}{2},2}(\partial
B^{3},{\mathbb{S}}^{2})$ admits an extension $u\in
W^{1,2}(B^{3},{\mathbb{S}}^{2})$; see [6, Theorem 6.2]. Once again, we
emphasize that this is not an immediate consequence of the analogue property
of linear Sobolev spaces. For example, there exists a boundary datum
$\varphi\in W^{\frac{1}{2},2}(\partial B^{3},{\mathbb{S}}^{1})$ which has _no_
extension $u\in W^{1,2}(B^{3},{\mathbb{S}}^{1})$; see [6, 6.3].

Minimizing harmonic maps satisfy the following system of Euler–Lagrange
equations

(1.2) $\left\\{\begin{array}[]{rcll}-\Delta u&=&|\nabla u|^{2}u&\mbox{in
$B^{3}$,}\\\ u&=&\varphi&\mbox{on $\partial B^{3}$}\text{.}\end{array}\right.$

It is known that for every non-constant boundary datum, the system (1.2)
admits infinitely many solutions; see [12]. Minimizers of (1.1) are not the
only solutions to (1.2) (see, e.g., [5, Section 3]). However, even in the
class of minimizing harmonic maps, we do not have uniqueness for a given
boundary datum $\varphi\colon B^{3}\to{\mathbb{S}}^{2}$; there are many known
examples. To list a few:

  * •

in [3, Section 3], there is an example of a planar boundary datum which admits
two different minimizers, one with values on the southern hemisphere and the
other one with values on the northern hemisphere;

  * •

in [4, 2.2. Corollary], there is an example of a boundary datum for which
there exists a 1-parameter family of distinct energy minimizing maps;

  * •

in [7, Section 5], there is an example of a boundary map which serves as a
boundary datum for at least two minimizers, one singular and the other one
regular;

  * •

in [1, 5.5 Theorem], there is an example of a boundary datum with mirror
symmetry for which there are at least two different minimizers without the
mirror symmetry.

Nevertheless, in the class of minimizing harmonic maps, we have the following
_generic uniqueness_ result ([1] attributes this theorem to Almgren).

######  Theorem 1.1 ([1, Theorem 4.1]).

Let $\varphi\in W^{1,2}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$. For every
$\varepsilon>0$, there exists $\psi\in
W^{1,2}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ such that
$\left\|{\varphi-\psi}\right\|_{W^{1,2}({\mathbb{S}}^{2})}<\varepsilon$ and
for which there exists exactly one energy minimizer $u\colon
B^{3}\to{\mathbb{S}}^{2}$ having boundary datum $\psi$. Moreover, $\psi$
coincides with $\varphi$ outside of $B_{\varepsilon}(x)\cap{\mathbb{S}}^{2}$,
for some $x\in{\mathbb{S}}^{2}$.

In [11], the second-named author and Strzelecki suspected that _generic non-
uniqueness_ occurs, when taking into account small perturbation of the
boundary datum in the topology of the space $W^{1,p}$ for $p<2$. The main
result of this note is the strengthening of [11, Remark 4.1].

######  Theorem 1.2.

Let $\varphi\in C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$. For every
$\varepsilon>0$, there exists $\psi\in
C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ such that
$\left\|{\varphi-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{2},{\mathbb{S}}^{2})}<\varepsilon$
which serves as a boundary datum for at least two energy minimizing maps from
$B^{3}$ to ${\mathbb{S}}^{2}$ having a different number of singularities.

Otherwise stated, Theorem 1.2 asserts that boundary data for which non-
uniqueness occurs are dense in $W^{1,p}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$.
This strengthens [7, Section 5] and [11, Remark 4.1], which provide existence
of _one_ boundary map for which non-uniqueness occurs. To be precise, as it is
stated, Theorem 1.2 only asserts that boundary data subjected to non-
uniqueness are dense in $C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ with
respect to the $W^{1,p}$ topology. In turn,
$C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ is dense in
$W^{1,p}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ (see e.g. [2, Theorem 1]), which
ensures the density of boundary data for which non-uniqueness occurs in the
whole $W^{1,p}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$.

Both Theorem 1.1 and Theorem 1.2 are in line with the _stability_ results: On
one hand, it is known that small perturbations of boundary data (for which
there is a unique minimizer) in the $W^{1,2}$ norm do not change the number of
singularities for corresponding minimizers (see [7] for perturbations in the
$W^{1,\infty}$ norm, [10] and [8] for perturbations in the $W^{1,2}$ norm). On
the other hand, small perturbations of the boundary datum in the $W^{1,p}$
norm for $p<2$ can change the number of singularities for corresponding
minimizers [11].

We prove Theorem 1.2 in Section 3. To do so, roughly speaking, we follow an
example by Hardt–Lin [7, Section 5]. We start with any smooth boundary datum
and use the construction of a boundary map (homotopic to the original one) of
[11] (see [9] for necessary modifications) for which a _Lavrentiev gap
phenomenon_ occurs. In Section 2, we show that a homotopy between these two
maps can be chosen small in $W^{1,p}$-norm for $p<2$, which is the novelty of
this note, and prove that within this homotopy, there is a boundary datum with
the required properties.

As we explained, our key contribution in this note, which allows the
transition from the existence to the density of boundary data where non-
uniqueness occurs, is the homotopy construction presented in Section 2. We
conclude this introduction with some extra comments concerning this
construction.

Assume that one is given $1\leq p<2$ and two maps $\varphi$ and $\psi\in
C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ that have the same topological
degree. Therefore, there exists a continuous, and even smooth homotopy
connecting $\varphi$ to $\psi$. A natural question is whether or not, knowing
that $\varphi$ and $\psi$ are close with respect to the $W^{1,p}$ distance,
one can choose the homotopy between $\varphi$ and $\psi$ to remain close to
$\varphi$ and $\psi$ all along the deformation. More precisely, one could for
instance expect that there exists a constant $C>0$ depending on $p$ such that
a homotopy $H\in C^{\infty}({\mathbb{S}}^{2}\times[0,1],{\mathbb{S}}^{2})$
between $\varphi$ and $\psi$ can be chosen so that

(1.3) $\left\|{\varphi-H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{2})}\leq
C\left\|{\varphi-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{2})}\quad\text{for every
$0\leq t\leq 1$.}$

Here, $H_{t}$ stands for the map $H(\cdot,t)$. The question is already
interesting if we assume in addition that $\varphi$ and $\psi$ coincide
outside of a small disk. For instance, one could ask whether or not a homotopy
such that (1.3) holds can be found under the additional assumption that
$\varphi=\psi$ outside of a ball of radius $r$, for some $r>0$ sufficiently
small, possibly depending on the map $\varphi$ that would be fixed in advance.

We are not able to solve this question, and a precise statement of the problem
in a more general context is given as 2.3. However, we are able to solve a
weaker version of this problem, which is nevertheless sufficient for our
purposes. Namely, we prove that, if the maps $\varphi$ and $\psi$ coincide
outside of a small ball, then a smooth homotopy between them can be found such
that $\left\|{\varphi-H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{2})}$ is
controlled, not by the distance between $\varphi$ and $\psi$, but by the sum
of their norms on a neighborhood of the region where they differ. This is the
content of the main result of Section 2, Proposition 2.1. This allows us to
deduce that, for a fixed $\varphi$ and a given $\varepsilon>0$, one can choose
the radius $r>0$ sufficiently small such that, for any map $\psi$ sufficiently
close to $\varphi$ such that $\varphi=\psi$ outside of $B_{r}(x)$, a homotopy
$H$ connecting $\varphi$ to $\psi$ can be found such that

$\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{2})}\leq\varepsilon\quad\text{for every
$0\leq t\leq 1$;}$

see Corollary 2.2. This is sufficient to prove our main result, Theorem 1.2,
but does not solve 2.3, as in our proof the radius $r>0$ of the ball outside
of which the maps $\varphi$ and $\psi$ are required to coincide has to depend
on $\varepsilon$, ruling out the possibility of controlling $\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{2})}$ uniformly in $t$ solely by
$\left\|{\varphi-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{2})}$ with our argument.

Notation. We denote by $B^{3}$ a Euclidean unit ball in $\mathbb{R}^{3}$. We
will write ${\mathbb{S}}^{n}$ for the unit $n$-dimensional sphere. For a point
$x\in{\mathbb{S}}^{n}$ and $r>0$, we will write $B_{r}(x)$ for a geodesic ball
of radius $r$ around $x$. We will write $A\precsim B$ whenever there is a
constant $C$ (independent of all crucial quantities) such that $A\leq CB$.
Throughout this paper, the term _minimizer_ will always refer to an
${\mathbb{S}}^{2}$-valued mapping minimizing the Dirichlet energy with given
boundary datum.

Acknowledgments.

  * •

The project is co-financed by the Polish National Agency for Academic Exchange
within Polish Returns Programme - BPN/PPO/2021/1/00019/U/00001 (KM).

  * •

The project is co-financed by National Science Centre, Poland grant
2022/01/1/ST1/00021 (KM).

##  2\. Homotopy construction

We will assume in this section that ${\mathcal{N}}$ is a (non necessarily
compact) Riemannian manifold. We work on the sphere ${\mathbb{S}}^{n}$, but
the result may be readily extended to an arbitrary domain, either an open
subset of $\mathbb{R}^{n}$ or a Riemannian manifold $\mathcal{M}$ of dimension
$n$. We also always assume that $p<n$.

######  Proposition 2.1.

Let $\varphi\in C^{\infty}({\mathbb{S}}^{n},{\mathcal{N}})$ and $p<n$. For
every $r>0$, for every $x\in{\mathbb{S}}^{n}$, and every $\psi\in
C^{\infty}({\mathbb{S}}^{n},{\mathcal{N}})$ homotopic to $\varphi$ and
satisfying $\varphi=\psi$ on ${\mathbb{S}}^{n}\setminus B_{r}(x)$, there
exists a homotopy $H\in C^{\infty}({\mathbb{S}}^{n}\times[0,1],{\mathcal{N}})$
from $\varphi$ to $\psi$ such that

$\sup_{0\leq t\leq 1}\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{n})}\leq
C\left(\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}+\left\|{\psi}\right\|_{W^{1,p}(B_{2r}(x))}\right)\text{,}$

for some constant $C>0$ depending only on $n$ and $p$.

This proposition can be used in combination with Lebesgue’s lemma to obtain a
homotopy which remains close to $\varphi$ in $W^{1,p}$. Indeed, choosing $r$
sufficiently small, depending on $\varphi$, we may ensure that
$\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}$ is as small as we want,
uniformly with respect to $r$. Since
$\left\|{\psi}\right\|_{W^{1,p}(B_{2r}(x))}\leq\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}+\left\|{\varphi-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{n})}$,
assuming in addition that
$\left\|{\varphi-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{n})}$ is small, we can
make $\sup_{0\leq t\leq 1}\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{n})}$ as small as we want. This yields
the following corollary.

######  Corollary 2.2.

Let $\varphi\in C^{\infty}({\mathbb{S}}^{n},{\mathcal{N}})$ and $p<n$. For
every $\varepsilon>0$, there exists $r>0$ sufficiently small, depending on
$\varphi$, and there exists $\delta>0$ such that, for every
$x\in{\mathbb{S}}^{n}$ and every $\psi\in
C^{\infty}({\mathbb{S}}^{n},{\mathcal{N}})$ homotopic to $\varphi$ and
satisfying $\varphi=\psi$ on ${\mathbb{S}}^{n}\setminus B_{r}(x)$ and
$\left\|{\varphi-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{n})}\leq\delta$, there
exists a homotopy $H\in C^{\infty}({\mathbb{S}}^{n}\times[0,1],{\mathcal{N}})$
from $\varphi$ to $\psi$ such that

$\sup_{0\leq t\leq 1}\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{n})}\leq\varepsilon\text{.}$

###### Proof of Proposition 2.1.

Let $G\in C^{\infty}({\mathbb{S}}^{n}\times[0,1],{\mathcal{N}})$ be any
homotopy connecting $\varphi$ to $\psi$ with $G_{0}=\varphi$ and $G_{1}=\psi$.
Since $\varphi=\psi$ outside of $B_{r}(x)$, we may assume that $G$ is
stationary outside of $B_{r}(x)$, i.e., for each $t\in[0,1]$, we have
$G_{t}=\varphi=\psi$ on ${\mathbb{S}}^{n}\setminus B_{r}(x)$. Consider
$\tau>0$, which will be chosen sufficiently small at a later stage. We are
going to rescale $G$, $\varphi$, and $\psi$ from $B_{r}(x)$ to a smaller ball
$B_{\tau}(x)$, while keeping them unchanged outside of $B_{2r}(x)$. More
specifically, let $\left(\Phi_{t}\right)_{0\leq t\leq 1}$ be a family of
smooth diffeomorphisms of ${\mathbb{S}}^{n}$ such that $\Phi_{t}={\rm id\,}$
outside of $B_{2r}(x)$ and such that, on $B_{2r}(x)$, in the local chart given
by the exponential map around $x$, $\Phi_{t}$ is expressed as

$\begin{cases}\frac{rx}{\left(1-t\right)r+t\tau}&\text{if
$\left|x\right|\leq\left(1-t\right)r+t\tau$,}\\\
\frac{x}{\left|x\right|}\left(\frac{r}{2r-\left(1-t\right)r-t\tau}\left(\left|x\right|-\left(1-t\right)r-t\tau\right)+r\right)&\text{if
$\left(1-t\right)r+t\tau\leq\left|x\right|\leq 2r$.}\end{cases}$

We define $H\in C^{\infty}({\mathbb{S}}^{n}\times[0,1],{\mathcal{N}})$ by

$H_{t}\coloneqq\begin{cases}\varphi\circ\Phi_{3t}&\text{if $0\leq
t\leq\frac{1}{3}$,}\\\ G_{3\left(t-1/3\right)}\circ\Phi_{1}&\text{if
$\frac{1}{3}\leq t\leq\frac{2}{3}$,}\\\
\psi\circ\Phi_{1-3\left(t-2/3\right)}&\text{if $\frac{2}{3}\leq t\leq
1$.}\end{cases}$

Of course, $H$ is a homotopy from $\varphi$ to $\psi$. It remains to show
that, if $\tau>0$ is suitably small, then $H$ satisfies the required estimate.

For $0\leq t\leq\frac{1}{3}$, we note that $\varphi-\varphi\circ H_{t}=0$
outside $B_{2r}(x)$. We readily obtain bounds on the Jacobian and the
derivatives of $\Phi_{t}$, so that the change of variable theorem combined
with $n-p>0$ implies that

$\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{n})}\leq\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}+\left\|{\varphi\circ\Phi_{3t}}\right\|_{W^{1,p}(B_{2r}(x))}\precsim\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}\text{.}$

Similarly, for $\frac{2}{3}\leq t\leq 1$, we have

$\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{n})}\leq\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}+\left\|{\psi\circ\Phi_{3t}}\right\|_{W^{1,p}(B_{2r}(x))}\lesssim\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}+\left\|{\psi}\right\|_{W^{1,p}(B_{2r}(x))}\text{.}$

Concerning $\frac{1}{3}\leq t\leq\frac{2}{3}$, we estimate

$\begin{split}\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{n})}&\leq\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}+\left\|{G_{3\left(t-1/3\right)}\circ\Phi_{1}}\right\|_{W^{1,p}(B_{2r}(x))}\\\
&\lesssim\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}+\left\|{G_{3\left(t-1/3\right)}}\right\|_{W^{1,p}(B_{2r}(x)\setminus
B_{r}(x))}+\tau^{\frac{n-p}{p}}\left\|{G_{3\left(t-1/3\right)}}\right\|_{W^{1,p}(B_{2r}(x))}\text{.}\end{split}$

Since the homotopy $G$ has been assumed to be stationary outside of
$B_{r}(x)$, we know that
$\left\|{G_{3\left(t-1/3\right)}}\right\|_{W^{1,p}(B_{2r}(x)\setminus
B_{r}(x))}=\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x)\setminus B_{r}(x))}$.
On the other hand, by compactness, we have

$\sup_{0\leq t\leq 1}\left\|{G_{t}}\right\|_{W^{1,p}(B_{2r}(x))}\leq C_{1}$

for some possibly large constant $C_{1}>0$. We may assume that either
$\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}\neq 0$ or
$\left\|{\psi}\right\|_{W^{1,p}(B_{2r}(x))}\neq 0$. Indeed, if
$\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}=0=\left\|{\psi}\right\|_{W^{1,p}(B_{2r}(x))}$,
this implies that both $\varphi$ and $\psi$ are identically zero — note that
this may only happen if $0\in{\mathcal{N}}$ — and we may directly conclude by
choosing $H$ to be constantly zero. As $p<n$, we may therefore choose $\tau>0$
sufficiently small, depending on $C_{1}$, so that

$\tau^{\frac{n-p}{p}}\left\|{G_{3\left(t-1/3\right)}}\right\|_{W^{1,p}(B_{2r}(x))}\leq\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}+\left\|{\psi}\right\|_{W^{1,p}(B_{2r}(x))}\quad\text{for
every $\frac{1}{3}\leq t\leq\frac{2}{3}$.}$

Hence, we deduce that

$\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{n})}\lesssim\left\|{\varphi}\right\|_{W^{1,p}(B_{2r}(x))}+\left\|{\psi}\right\|_{W^{1,p}(B_{2r}(x))}\quad\text{for
every $\frac{1}{3}\leq t\leq\frac{2}{3}$.}$

This concludes the proof. ∎

In Corollary 2.2, both the $\delta>0$ controlling
$\left\|{\varphi-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{n})}$ _and_ the $r>0$
depend on $\varepsilon$. A very natural question is whether or not one may
find a homotopy $H$ so that $\sup_{0\leq t\leq 1}\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{n})}$ is controlled only by
$\left\|{\varphi-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{n})}$. More precisely,
we formulate the following open question (cf. [11, Problem, p.11]).

######  Open Problem 2.3.

Let $\varphi\in C^{\infty}({\mathbb{S}}^{n},{\mathcal{N}})$. Does there exist
some $r>0$, possibly depending on $\varphi$, such that for every
$x\in{\mathbb{S}}^{n}$ and every $\psi\in
C^{\infty}({\mathbb{S}}^{n},{\mathcal{N}})$ homotopic to $\varphi$ and
satisfying $\varphi=\psi$ on ${\mathbb{S}}^{n}\setminus B_{r}(x)$, there
exists a homotopy $H\in C^{\infty}({\mathbb{S}}^{n}\times[0,1],{\mathcal{N}})$
from $\varphi$ to $\psi$ such that

$\sup_{0\leq t\leq 1}\left\|{\varphi-
H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{n})}\leq\omega\left(\left\|{\varphi-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{n})}\right)\text{,}$

where $\omega$ is a modulus of continuity satisfying $\omega\left(t\right)\to
0$ as $t\to 0$.

One may expect $\omega$ to be linear in $t$, but any modulus of continuity
would already be of interest. The question is already interesting for maps
${\mathbb{S}}^{2}\to{\mathbb{S}}^{2}$.

##  3\. Proof of the generic non-uniqueness

###### Proof of Theorem 1.2.

Fix $\varepsilon>0$ and $\varphi\in
C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$. We note first that, by Theorem
1.1 combined with Hölder’s inequality, we may find another mapping
$\varphi_{0}\in C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ which admits
exactly one energy minimizer $u_{0}\colon B^{3}\to{\mathbb{S}}^{2}$ among all
maps having boundary datum $\varphi_{0}$, and such that $\varphi_{0}$ differs
from $\varphi$ only on a set $B_{\frac{\varepsilon}{2}}(x_{0})$ for some
$x_{0}\in{\mathbb{S}}^{2}$ and is such that

(3.1)
$\left\|{\varphi-\varphi_{0}}\right\|_{W^{1,p}({\mathbb{S}}^{2})}<\frac{\varepsilon}{2}\text{.}$

We recall that, combining the regularity result [13, Theorem II] with the
boundary regularity [14, Theorem 2.7] of Schoen–Uhlenbeck, $u_{0}$ can have
only a finite number of singularities; let us denote this number by
$M=\\#\operatorname{sing}u$ (possibly $M=0$).

Next, we apply Corollary 2.2 to $\varphi_{0}\in
C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$. We obtain the existence of a
$\delta=\delta(\varepsilon)>0$ and an $r=r(\varphi_{0},\varepsilon)>0$ such
that for any $\psi\in C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ that
differs from $\varphi_{0}$ only on the set $B_{r}(x_{0})$ and such that
$\left\|{\varphi_{0}-\psi}\right\|_{W^{1,p}({\mathbb{S}}^{2})}<\delta$, there
exists a homotopy $H\in
C^{\infty}({\mathbb{S}}^{2}\times[0,1],{\mathbb{S}}^{2})$ with

(3.2) $\sup_{0\leq t\leq
1}\left\|{\varphi_{0}-H_{t}}\right\|_{W^{1,p}({\mathbb{S}}^{2})}<\frac{\varepsilon}{2}.$

Let $\varepsilon_{1}\coloneqq\min\\{\delta,r,\frac{\varepsilon}{2}\\}$. By [9,
Theorem 2.3.1], we construct $\varphi_{1}\in
C^{\infty}({\mathbb{S}}^{2},{\mathbb{S}}^{2})$ with the properties:

  1. (1)

$\deg\varphi_{0}=\deg\varphi_{1}$;

  2. (2)

$\|\varphi_{0}-\varphi_{1}\|_{W^{1,p}}<\varepsilon_{1}$ and
$\varphi_{0}=\varphi_{1}$ except on $B_{\varepsilon_{1}}(x)$ for some point
$x\in{\mathbb{S}}^{2}$;

  3. (3)

$\varphi_{1}$ admits only one energy minimizer $u_{1}\colon
B^{3}\to{\mathbb{S}}^{2}$ having at least $M+1$ singularities.

To be precise, the statement [9, Theorem 2.3.1] gives only that
$\mathcal{H}^{2}(\\{x\in{\mathbb{S}}^{2}\colon\varphi_{0}(x)\neq\varphi_{1}(x)\\})<\varepsilon_{1}$,
but following the lines of the proof, we may deduce that
$\varphi_{0}=\varphi_{1}$ except on $B_{\varepsilon_{1}}(x)$ for some point
$x\in{\mathbb{S}}^{2}$.

Now, let us take the homotopy $H_{t}$ between $\varphi_{0}$ and $\varphi_{1}$
constructed in Corollary 2.2. Let

$\begin{split}\tau\coloneqq\sup\\{t\in[0,1]\colon&\text{ each energy minimizer
with boundary datum $H_{t}$}\\\ &\text{ has at most $M$ singular points in
}B^{3}\\}\text{.}\end{split}$

We argue like in [11, Remark 4.1] (which is a modified argument from [7,
Section 5]). For the convenience of the reader, we state here the main lines
of the reasoning. First, we note that from the Stability Theorem [7], see also
[10, Theorem 8.9], we have $\tau\in(0,1)$.

Now take $s_{i}\nearrow\tau$ and a sequence of minimizing harmonic maps
$u_{i}\in W^{1,2}(B^{3},{\mathbb{S}}^{2})$ with $u_{i}\big{\rvert}_{\partial
B^{3}}=H_{s_{i}}$ and $\\#\operatorname{sing}u_{i}\leq M$. Let us also take
$t_{i}\searrow\tau$ and a sequence of minimizing harmonic maps $v_{i}\in
W^{1,2}(B^{3},{\mathbb{S}}^{2})$ with $v_{i}\big{\rvert}_{\partial
B^{3}}=H_{t_{i}}$ and $\\#\operatorname{sing}v_{i}>M$. Since
$\sup_{i}\left([H_{s_{i}}]_{W^{1,2}({\mathbb{S}}^{2})}+[H_{t_{i}}]_{W^{1,2}({\mathbb{S}}^{2})}\right)<\infty$,
we may deduce from the strong convergence of minimizers, see [1, Theorem 1.2
(4)] (see also [10, Theorem 6.1 (3)]), that up to a subsequence we have

$\begin{split}u_{i}&\to u\quad\text{ strongly in
}W^{1,2}(B^{3},{\mathbb{S}}^{2})\text{,}\\\ v_{i}&\to v\quad\text{ strongly in
}W^{1,2}(B^{3},{\mathbb{S}}^{2})\text{,}\end{split}$

and both $u$ and $v$ are energy minimizers with $u\big{\rvert}_{\partial
B^{3}}=v\big{\rvert}_{\partial B^{3}}=H_{\tau}$. We claim that
$\\#\operatorname{sing}u\leq M$. Indeed, assume on the contrary that
$\\#\operatorname{sing}u>M$. Then, by [1, Theorem 1.8 (2)] (see also [10,
Theorem 2.10]), we would obtain that for each $y\in\operatorname{sing}u$ and
for sufficiently large $i$, there would exist
$y_{i}\in\operatorname{sing}u_{i}$ with $y_{i}\to y$ as $i\to\infty$, a
contradiction.

Moreover, $\\#\operatorname{sing}v>M$. To see this, let us again assume by
contradiction that $\\#\operatorname{sing}v\leq M$. Let now
$z_{i,j}\in\operatorname{sing}v_{i}$ for $j\in\\{1,\ldots,M+1\\}$ be distinct
singular points of $v_{i}$. Now let us observe that for sufficiently large
$i$, we know that that $H_{t_{i}}$ and $H_{\tau}$ are close in $C^{\infty}$.
Hence, by uniform boundary regularity [1, Theorem 1.10 (2)] (see also [10,
Theorem 7.4]), there is a uniform neighborhood of the boundary $\partial
B^{3}$ which contains no singularities of $v$ and $v_{i}$, say ${\rm
dist\,}(z,\partial B^{3})\geq\lambda>0$ for any
$z\in\bigcup_{i}\operatorname{sing}v_{i}\cup\operatorname{sing}v$. Since
singular points converge to singular points, we deduce from [1, Theorem 1.8
(1)] (see also [10, Theorem 2.5]) that for each $j$, we have $z_{i,j}\to
z_{j}$ as $i\to\infty$ and $z_{j}\in\\#\operatorname{sing}v$. The only
possibility for $\\#\\{z_{1},\ldots,z_{M+1}\\}<M+1$ is that two singularities
of $v_{i}$ converge to the same singularity of $v$. This, however, is
impossible, because by the uniform distance between singularities [1, Theorem
2.1] (see also [10, Theorem 2.12]), there exists a universal constant $C$
(independent of the minimizer) such that no singularity can occur next to
$z_{i,j}$ at a distance $C\,{\rm dist\,}(z_{i,j},\partial B^{3})\geq
C\lambda$.

Hence, $H_{\tau}\colon{\mathbb{S}}^{2}\to{\mathbb{S}}^{2}$ serves as a
boundary condition for at least two minimizers $u$ and $v$ having a different
number of singularities. Combining (3.2) with (3.1), we obtain

$\left\|{\varphi-
H_{\tau}}\right\|_{W^{1,p}({\mathbb{S}}^{2})}\leq\left\|{\varphi-\varphi_{0}}\right\|_{W^{1,p}({\mathbb{S}}^{2})}+\left\|{\varphi_{0}-H_{\tau}}\right\|_{W^{1,p}({\mathbb{S}}^{2})}<\frac{\varepsilon}{2}+\varepsilon_{1}\leq\varepsilon.$

This finishes the proof. ∎

## References

  * [1] F. J. Almgren, Jr. and E. H. Lieb.  Singularities of energy minimizing maps from the ball to the sphere: Examples, counterexamples, and bounds.  Ann. of Math. (2), 128(3):483–530, 1988.
  * [2] F. Bethuel and X. Zheng.  Density of smooth functions between two manifolds in Sobolev spaces.  J. Funct. Anal., 80(1):60–75, Sept. 1988.
  * [3] R. Hardt and D. Kinderlehrer.  Mathematical questions of liquid crystal theory.  In Nonlinear partial differential equations and their applications. Collège de France Seminar, Vol. IX (Paris, 1985–1986), volume 181 of Pitman Res. Notes Math. Ser., pages 276–289. Longman Sci. Tech., Harlow, 1988.
  * [4] R. Hardt, D. Kinderlehrer, and F.-H. Lin.  The variety of configurations of static liquid crystals.  In Variational methods (Paris, 1988), volume 4 of Progr. Nonlinear Differential Equations Appl., pages 115–131. Birkhäuser Boston, Boston, MA, 1990.
  * [5] R. Hardt, D. Kinderlehrer, and M. Luskin.  Remarks about the mathematical theory of liquid crystals.  In Calculus of variations and partial differential equations (Trento, 1986), volume 1340 of Lecture Notes in Math., pages 123–138. Springer, Berlin, 1988.
  * [6] R. Hardt and F.-H. Lin.  Mappings minimizing the $L^{p}$ norm of the gradient.  Comm. Pure Appl. Math., 40(5):555–588, 1987.
  * [7] R. Hardt and F.-H. Lin.  Stability of singularities of minimizing harmonic maps.  J. Differential Geom., 29(1):113–123, 1989.
  * [8] S. Li.  Stability of minimising harmonic maps under $W^{1,p}$ perturbations of boundary data: $p\geq 2$.  J. Differential Equations, 296:279–298, 2021.
  * [9] K. Mazowiecka.  Singularities of harmonic and biharmonic maps into compact manifolds.  PhD thesis, University of Warsaw, 2017.
  * [10] K. Mazowiecka, M. Miśkiewicz, and A. Schikorra.  On the size of the singular set of minimizing harmonic maps.  Mem. Amer. Math. Soc., accepted.
  * [11] K. Mazowiecka and P. Strzelecki.  The Lavrentiev gap phenomenon for harmonic maps into spheres holds on a dense set of zero degree boundary data.  Adv. Calc. Var., 10(3):303–314, 2017.
  * [12] T. Rivière.  Infinités des applications harmoniques à valeur dans une sphère pour une condition au bord donnée.  Published in Applications harmoniques entre variétés.  PhD thesis, Université Paris VI, 1993.
  * [13] R. Schoen and K. Uhlenbeck.  A regularity theory for harmonic maps.  J. Differential Geom., 17(2):307–335, 1982.
  * [14] R. Schoen and K. Uhlenbeck.  Boundary regularity and the Dirichlet problem for harmonic maps.  J. Differential Geom., 18(2):253–268, 1983.


# Hydrodynamically Induced Aggregation in Two-Dimensional Active Systems

Roee Bashan and Naomi Oppenheimer  School of Physics and Astronomy and the
Center for Physics and Chemistry of Living Systems, Tel Aviv University

###### Abstract

We investigate a system of co-oriented active particles interacting only via
hydrodynamic and steric interactions. We offer a new method of calculating the
flow created by any active particle in a 2D fluid, focusing on the dynamics of
flow fields with a high-order spatial decay, which we analyze using a
geometric Hamiltonian. We show that when orientational degrees of freedom are
quenched, and the flow has a single, odd power decay, such many-particle
systems lead to stable, fractal-like aggregation, with the only exceptions
being the force dipole. We discuss how our results can easily be generalized
to more complicated force distributions and to other effective two-dimensional
systems.

Out-of-equilibrium ensembles, both biological and artificial, are often
suspended in a fluid. Many times, these active systems have lower dimensions
[1, 2, 3, 4]. A prominent example is a cell membrane — an effective two-
dimensional (2D) fluid in which both passive and active proteins are immersed
[5, 6]. Active membrane proteins act as either rotary proteins, such as ATP
synthase [7, 8, 9], or as shakers — particles that are not self-propelled but
apply active forces on the membrane due to conformational changes,
polymerization, or reorganization [10]. Artificial 2D systems are also
widespread, such as active colloids driven by a chemical reaction [11, 12,
13], by light [14, 15] or by an external magnetic field [16]. Particle
dynamics can be dominated by complex interactions (e.g. by electrostatics,
capillary forces or an external field) [17, 18], or by purely hydrodynamic
flows [19, 20, 21, 22]. In all the above systems, understanding the connection
between flow and structure is crucial.

The flow created by an active particle can be described by a multipole
expansion (similar to electrostatics). Far away from the particle, the leading
order of the flow, the stokeslet, is given by the total force the particle
applies. Closer to the particle, higher multipole contributions start to
dominate. First, the force dipole, which could be divided to the rotlet [7,
23, 24, 25, 9], and the stresslet [26, 27, 28, 29]. Closer still, the
quadrupolar term is significant, then the octopole, and so forth [30, 31, 32].
Thus, by first finding the effect of a single point-force, the so-called
Green’s function of the problem, we can derive the general flow response. This
procedure applies even when there is no net force, in which case the higher
orders dominate at larger distances.

In this work, we explore the dynamics of suspended particles in a 2D viscous
fluid with quenched orientational degrees of freedom, interacting only via the
flow they create at their boundary and steric interactions. We review how the
multipole terms mentioned above are often calculated from the stokeslet, but
then, we focus on improving the common procedure, making higher-order terms
easier to use and calculate. We show that, in many cases, namely for all even
terms of the multiple expansion, the dynamics could be described by a
Hamiltonian. This Hamiltonian is geometric in nature, with the conjugate
variables being $x_{i}$ and $y_{i}$, the position of the $i^{\rm th}$
particle. Phase space in such cases corresponds to configurational space.
Thus, by limiting possible paths in phase-space we can determine the steady-
state distributions such active particles can take. We prove that, in many
cases, particles aggregate. In particular, in all far-field limits, which are
dominated by a single multipole term with an odd power decay. The only
exceptions are the rotlet and the stresslet. We corroborate our findings by
simulating particles that interact by an octopolar force distribution. Figure
(1) illustrates our key finding: a pure multipole force distribution leads to
particle aggregation. Moreover, we claim that our findings can be applied to
more general 2D cases — both for calculating high-order interactions in more
complex fluids, and for understanding their resulting dynamics in many-
particle ensembles.

Figure 1: Snapshots of a simulation of 500 active particles with an octopole
force distribution (i.e., the velocity decays as $1/r^{3}$). Particles are
initialized randomly in a $60\times 60$ box; Their orientation is quenched.
The particles collide quickly only due to hydrodynamic interactions and form
small clusters. At long times, clusters also collide, creating larger
clusters. Finally, The particles form a single, fractal-like, cluster
organized in a $45^{\circ}$ angled square lattice.

The Reynolds number of a flow, $\text{Re}=\frac{\rho UL}{\mu}$, is the ratio
between the inertia of the fluid and the viscous forces in its flow, where
$\rho$ is the fluid’s density, $\mu$ is its viscosity, $U$ is the
characteristic velocity and $L$ the characteristic length. For microswimmers
in biological systems, the Reynolds number is small, as the characteristic
lengths are microscopic. In the low Reynolds number limit, the inertia is thus
negligible, and the flow velocity v(r) (where ${\bf r}=(x,y)$ is the 2D
position vector) is governed by the Stokes equations,

$\begin{cases}\nabla p=\mu\nabla^{2}\textbf{v}+\textbf{f}\\\
\nabla\cdot\textbf{v}=0,\end{cases}$ (1)

where f is an external force. From now on, we will work in unitless notation
where $\mu=1$. The flow due to a point force, the stokeslet, also termed the
Oseen tensor in 2D [33], is given by

$G^{\rm
2D}_{ij}=\frac{1}{4\pi}\left(-\log\frac{r}{r_{0}}\delta_{ij}+\frac{r_{i}r_{j}}{r^{2}}\right),$
(2)

where $r=|\textbf{r}|$, and $r_{0}$ is a normalizing constant, e.g. the size
of the system (it is insignificant for our purposes since it does not appear
in higher order terms). This solution is the Green’s function of the Stokes
equations and can be used to calculate the flow due to any force distribution.
However, an exact calculation is generally difficult, and, when particles are
far apart, unnecessary. Instead, one often approximates the flow in powers of
$r$ using a Taylor expansion of $G^{\rm 2D}$, giving a solution of the form,

$\textbf{v}=\bm{F}\cdot\bm{G}+\bm{D}_{2}\bm{:}\nabla\bm{G}+\cdots+\bm{D}_{N}\raisebox{0.5pt}{\raisebox{-.9pt}
{N}⃝}\nabla^{\otimes N-1}\bm{G}+\cdots.$ (3)

where $\bm{:}$ is the second order inner product, \raisebox{-.9pt} {N}⃝ is the
$N^{th}$ order inner product, $\otimes$ is the tensor product (which, when
found in the exponent, refers to the times to apply the operator), F is the
total force, and $\textbf{D}_{n}$ is the $n^{\rm th}$ force coefficient given
by

$\bm{D}_{n+1}=\frac{(-1)^{n}}{n!}\int_{S}\bm{f}\otimes\bm{r}^{\otimes n}dA.$
(4)

The integration is over the area where the forces are applied. Using the
monopole or dipole terms for solutions far away from the particle is often
enough. Calculations of higher-order terms are messy at best. To evaluate the
$N^{\rm th}$ term, first one needs to calculate the mixed $N-1$ order
derivatives, then calculate $N$ coefficients of $\textbf{D}_{N}$, and finally
contract these two organized in $N-1$ rank and $N$ rank tensors in a slaving
procedure. The contraction result is often complicated in Cartesian
coordinates but is much simpler in polar coordinates. A similar formulation is
also used in 3D [31], with the equivalent $G^{\rm 3D}$. In Sec. I, we will
show a different route to calculating the flow using the biharmonic equation.

Let us now consider a fluid in which an active particle is immersed. The
particle is propagated by the flow only (we will later add steric
interactions), and so its position is given by a simple equation of motion
$\dot{\textbf{r}}=\textbf{v}(\textbf{r})$. Since the flow is incompressible
and has no sources of mass (Eq. 1), it can be described by a potential $\psi$
called the stream function. In 2D,
$\nabla\times(\psi\hat{z})=\textbf{v}\equiv\nabla^{\perp}\psi$, where
$\nabla^{\perp}=-\hat{z}\times\nabla=\hat{x}\partial_{y}-\hat{y}\partial_{x}$.
In the case of a flock of similar particles, all with the same orientation,
each with strength $S_{i}$, the stream function created by all particles is a
simple superposition,
$\psi(\textbf{r})=\sum_{i}S_{i}\psi(\textbf{r}-\textbf{r}_{i})$. The dynamics
of the particles could be described by a geometric Hamiltonian

$H=\frac{1}{2}\sum\limits_{i\neq
j}S_{i}S_{j}\psi(\textbf{r}_{i}-\textbf{r}_{j}).$ (5)

The equations of motions are given by Hamilton equations of $H$, with
$\sqrt{S_{i}}x_{i}$ and $\sqrt{S_{i}}y_{i}$ being conjugate variables. Such a
Hamiltonian description has been useful for point vortices in an ideal fluid
[34, 35] and also, more recently, in viscous systems [36, 37, 27, 38, 39, 40].
However, this formalism only holds when $\psi$ is an even function of r. An
odd component of $\psi$ will cancel out in the (symmetric) summation, and $H$
would be identically zero. In the rest of this work, we will mostly use the
case where the particles are identical, i.e. $S_{i}=1$.

The rest of the paper is organized as follows. In Sec. I we will derive the 2D
flow created by an active particle with any force distribution acting on its
surface. In Sec. II, we will derive the dynamics of two similar particles with
quenched orientational degrees of freedom. As an example, we will show that
when the force distribution is octopolar, dominated by ${\bf D_{4}}$,
particles will collide at a finite time. We will then prove that, except for a
force-dipole, all interactions of a single, even, multipole term lead to
aggregation in the two-particle system. We will also discuss more general
interactions. Our arguments will also be applicable to other similar systems
with flow fields not described by the equations of flow above. In Sec. III, we
will claim that the two-particle case also applies to an ensemble of many
particles and show a simulation resulting in aggregation.

##  I flow due to a single particle

Here, we will derive the stream function generated by a general active
particle in a 2D viscous fluid. For simplicity, we will assume a circular
particle with a radius $R$, applying a force on the fluid at its boundary.
However, the choice that forces are only applied on a circle does not limit
the result of the stream function, and all the solutions in Eq. (3) will
appear here as well (see Supplementary Information for more details). Taking
the curl of the first equation in Eqs. (1), and writing the velocity in terms
of $\psi$, we arrive at an equation for the stream function,

$\nabla^{2}\nabla^{2}\psi=-\nabla^{\perp}\cdot\textbf{f}=\hat{z}\cdot(\nabla\times\textbf{f}).$
(6)

The force is only present at the surface of the particle. We will assume all
the forces are parallel, giving a force of the form
$\textbf{f}=\frac{1}{R}\delta(r-R)F(\theta)\textbf{f}_{0}$ where
$\textbf{f}_{0}$ is a constant and $F(\theta)$ is an angular density
distribution. The choice that forces are parallel creates some restrictions,
but they are easily resolved for a force distribution with multiple
orientations (see SI). The homogeneous biharmonic equation has been used to
study the flow created by active swimmers with varied boundary conditions [41,
42]. Here, we focus on the coupling between the force applied by a swimmer and
the flow, which is given by the non-homogeneous biharmonic equation (6).

In order to get a solution equivalent to Eq. (3), we match the boundary
conditions on the flow at infinity and at $r=0$ (despite the particle having a
small radius $R$), similar to the derivation of the Oseen Tensor [33]. It is
possible to solve Eq. (6) by decomposing the force and the stream function in
a Fourier series as they are periodic in $\theta$,

$\psi(r,\theta)=\sum\limits_{n\in\mathbb{Z}}B_{n}(r)e^{in\theta}\ \ \ ,\ \ \
F(\theta)=\sum\limits_{n\in\mathbb{Z}}C_{n}e^{in\theta}.$ (7)

We will assume no force monopole, thus $C_{0}=0$. Using Eq. (7) in Eq. (6), we
get (see SI for more details),

$\hat{O}_{n}^{2}B_{n}(r)=i\frac{\eta}{R}C_{n-1}\left[\frac{d}{dr}\delta(r-R)-\frac{n-1}{r}\delta(r-R)\right],$
(8)

where
$\hat{O}_{n}=\frac{1}{r}\frac{d}{dr}\left(r\frac{d}{dr}\right)-\frac{n^{2}}{r^{2}}$
and $\eta$ is the force in complex notation
$\eta=|\mathbf{f_{0}}|e^{-i\arg(\textbf{f}_{0})}$. The solution is given only
by the real part of Eq. (17). This equation implies that a Fourier component
$n$ in the stream function $\psi$ is caused by a component $n-1$ of the force.
For example, a force caused by two opposing normal point forces on the surface
has only odd components, which leads to only even components in $\psi$. Every
Fourier component, $B_{n}(r)$, in the stream function $\psi$ is a linear
combination of four power laws in $r$:
$r^{n},r^{n+2},\frac{1}{r^{n}},\frac{1}{r^{n-2}}$, where only negative
exponents will contribute at $r>R$ (see SI for a detailed solution). Outside
the particle, at $r>R$ the solution is a linear combination of terms of the
form $\ln r,\sin(2\theta),\sin(n\theta)/r^{n},\sin(\theta)/r^{n-2}$ (ignoring
phases), where the logarithm is due to a degeneracy for $n=0$. The complete
solution outside the particles is given by

$\displaystyle\frac{4\psi}{iR\eta}=\sum\limits_{n>1}\left\\{C_{n-1}e^{in\theta}\left[\frac{1}{n}\left(\frac{R}{r}\right)^{n}+\frac{1}{(1-n)}\left(\frac{R}{r}\right)^{n-2}\right]\right.$
(9)
$\displaystyle+\left.\frac{C_{-n-1}}{n(n+1)}e^{-in\theta}\left(\frac{R}{r}\right)^{n}\right\\}-2C_{-1}\ln\frac{r}{R}+\frac{1}{2}C_{-2}e^{-i\theta}\frac{R}{r}.$

This solution allows us to calculate the velocity field created by any active
particle applying a force on the fluid in one direction. As mentioned earlier,
this result is easily generalized to any force orientation as outlined in the
SI. Furthermore, the terms in Eq. 9 can apply to active particles that have
different boundary conditions at the surface of the particle. For instance, a
squirmer is an active particle defined by the flow on its surface, which has
no radial component [43, 44]. Its 3D description is well-known. Similarly to
Eq. (9), a general 2D squirmer is given by:

$\psi=B_{0}\ln\frac{r}{R}+\sum\limits_{n>1}B_{n}\left(\frac{R^{n}}{r^{n}}-\frac{R^{n-2}}{r^{n-2}}\right)\sin(n\theta+\phi_{n})$
(10)

For any $B_{n}$ and $\phi_{n}$. These terms are of the same form as the
different multiples of a 2D flow given in Eq. 9, i.e.
$\sin(n\theta)/r^{n},\sin(n\theta)/r^{n-2}$.

##  II two particle dynamics

This section will focus on a particular example of the flow given by a single
force multipole, as given by Eq. (9), and the dynamics of two particles
interacting by such a flow. We will then discuss the dynamics of two particles
interacting by a general solution to Eq. (6), with any desired force
distribution. We will show that when the orientation is quenched, similar
particles driven by a pure, even, multipolar term will always collide, the
only exceptions being the force monopole and force dipole. For some general
force distribution, two particles will always attract in the far-field limit,
given that the total force is zero and there is no contribution from a force
dipole. If there are no local minima and maxima in the stream function, the
particles will also collide.

###  II.1 Case Study: Octopole Term

Figure 2: A realization of a particle whose stream function’s leading term is
$1/r^{2}$, which is the octopole term in Eq. (3), using point forces. We
investigate the interaction between two identical particles with this force
distribution, with a fixed orientation.

We will consider a force distribution such that only the third harmonic of the
force distribution is present. In equation (3), this corresponds to the forth
term, $\mathbf{D}_{4}$ dominating, and the stream function is given by:

$\psi_{3}=f_{0}C_{3}R\left[\left(\frac{R^{2}}{12r^{2}}-\frac{R^{4}}{16r^{4}}\right)\sin
4\theta+\frac{R^{2}}{24r^{2}}\sin 2\theta\right],$ (11)

assuming $C_{3}=C_{-3}$ is real, and the force is in the $x$ direction. In the
far-field limit, taking $R\rightarrow 0$ and $f_{0}\rightarrow\infty$ such
that $R^{3}f_{0}$ remains constant, only the terms $~{}1/r^{2}$ remain, which
are the octopole moment

$\psi_{\text{oct}}=\frac{R^{3}}{12r^{2}}C_{3}f_{0}\left(\frac{1}{2}\sin
2\theta+\sin 4\theta\right).$ (12)

There are many force distributions that have this property. The obvious one is
the pure harmonic written above. A more realistic realization given by point
forces is presented in Fig. 2. A distribution of forces on a disk can always
be constructed so that any term is dominant in the multipole expansion (see
SI).

From this point on, for simplicity, we will focus on one of the terms in Eq.
12 and take $\psi=S\sin(4\theta)/r^{2}$ (This can be realized with a force
distribution at two different directions, see SI), which represents well the
general case for the dynamics caused by a single (even) term. Its
corresponding flow is displayed in Fig. 3. A system of two identical particles
creating such a flow is described by the Hamiltonian in Eq. (5).

In Hamiltonian mechanics, symmetries of the Hamiltonian correspond to
conservation laws. In this case, the Hamiltonian is symmetric to translation
in time and space. Therefore, the Hamiltonian is conserved as well as, what we
term, the “center of activity”, $(\textbf{r}_{1}+\textbf{r}_{2})/2$, which is
analogous to the center of mass. It is therefore enough to consider the
dynamics of the relative difference between particles
$\textbf{d}=\textbf{r}_{1}-\textbf{r}_{2}$, which has a 2D phase space. The
paths d can take are displayed in Fig. 3, which are given by
$r^{2}\propto\sin(4\theta)$ as a result of the conservation of the Hamiltonian
$H=\psi(\textbf{d})$. It is also straightforward to solve the equations of
motion directly, giving

$\begin{cases}\cos(4\theta)=\cos(4\theta_{0})-16\frac{H^{2}}{S}t\\\
r^{4}=r_{0}^{4}-32St\cos(4\theta_{0})-256H^{2}t^{2},\end{cases}$ (13)

where $r$ and $\theta$ are the polar coordinates in the phase space of d.
Theta is constrained to only one of the intervals
$(0,\frac{\pi}{4}),(\frac{\pi}{4},\frac{\pi}{2})\cdots(\frac{7\pi}{4},2\pi)$
which are given by half the period of the harmonic. Also, in general, the
radius does not decrease monotonically with time, however it always decreases
to zero unless $H=0$. The special paths with $H=0$ are the only paths that
will not necessarily converge to the origin. In principle, they can diverge to
infinity, though a divergence is unstable as even a small perturbation to the
Hamiltonian will cause the path to be bound. Therefore, every stable solution
of the pair leads to a collision. We will now prove that, in the far-field
limit, this holds true in general, except for dipolar force terms.

Figure 3: The vector field generated by a stream function
$\psi=S\frac{\sin(4\theta)}{r^{2}}$, where blue corresponds to low strength
and red is high. In purple we mark example paths that a particle can take. The
straight paths are characterized by a zero Hamiltonian, and alternate being
stable or unstable. The leaf-shaped paths are for the non-zero case,
alternating a positive and negative Hamiltonian.

###  II.2 General dynamics of two Particles

Let us consider a case of two identical particles, each creating a stream
function $\psi$ that is described by a single (even) multipole term in the
far-field limit. We will later consider a dipolar flow which is given by
$\psi=\alpha\ln r+\beta\sin(2\theta+\phi)$. Any other flow is given by

$\psi=\left[\alpha\sin(m\theta+\varphi)+\beta\sin((m+2)\theta+\phi)\right]\frac{R^{m}}{r^{m}}$
(14)

for some even $m$, any coefficients $\phi,\varphi$ and any non-trivial
$\alpha$ and $\beta$. We will now show that all such flows (excluding the
dipole terms) will result in a collision of the two particles, similar to the
octopole term described above.

To prove this claim, we will examine the paths that the relative distance
between particles, d, takes in its phase-space, whose motion is governed by
$\psi$, $\dot{\textbf{d}}=2\nabla^{\perp}\psi(\textbf{d})$. The “grad perp”
operator, $\nabla^{\perp}$, has an important geometrical interpretation:
unlike the gradient, which points to the direction where a scalar function is
most increasing, the “grad perp” operator points to the direction
perpendicular to that (anti-clockwise) which is also the direction where the
function is constant. This interpretation is connected to the fact that $H$ is
conserved along paths.

Figure 4: The categories of possible paths for the separation between two-
active particles, d. If bound, the path either encircles the origin, does not
encircle it, or crosses it. If the path diverges, then it either crosses the
origin or not. Each path is described by an equation $\psi(x,y)=H$. In the
text we prove which paths are permitted and under what circumstances.

First, we note that the stream function has no local maxima/minima, as a cross
section of it at a constant $\theta$ gives $\psi\propto r^{-m}$, which never
has local extrema. Thus, a point path is unstable. This is expected since a
stable point would have a negative divergence of the velocity field. Thus,
every path curve must be either closed or diverging to infinity. And, since
$\psi(r\rightarrow\infty)=0$, every path diverging to infinity has a zero
Hamiltonian, and every path with non-zero Hamiltonian is necessarily bound and
closed. The different possible path categories are shown in Fig. 4. We will
now assess whether or not each is possible.

Path a. A closed path that does not contain the origin, such as path a,
encircles an area which is a compact subset of the phase-space, where $\psi$
is continuous. Because the stream function is constant along the edge of the
region, that is, along a path, there must be a local minimum or maximum inside
that region, contradicting our earlier conclusion (this is a higher
dimensional case of Rolle’s theorem [45]). This argument, which disallows
paths like a, does not hold for the closed paths like d and e because their
closed area contains the origin, and so the stream function is not continuous
on it.

Paths b and c. Now consider a path that diverges to infinity, like paths b and
c. As mentioned, such a path has zero Hamiltonian, so every path that crosses
it must also have zero Hamiltonian. Therefore, the curve will split the plane
in two — paths on one side of the curve that have a non-zero Hamiltonian
cannot cross to the other side. For paths like b, there is a side which does
not contain the origin, and so closed paths contained in it must be like path
a. Since we showed such paths are not allowed, it then follows that all paths
in that side have zero hamiltonian, and thus $\psi=0$. We assume that there
cannot be such an area in which the stream function is constant since that
corresponds to an area with no flow. That means every path diverging to
infinity must cross the origin, That is, it must be similar to path c. In
fact, a path similar to path c must exist, because there is some
$\theta_{0}\in[0,2\pi)$ such that $\psi(r,\theta_{0})=0$ for all $r$, and so
the ray $\theta=\theta_{0}$ is a path similar to path c.

Path d. Since a divergent path such as c necessarily exists, every closed path
encircling the origin (path d) crosses it. Thus, it also has a zero
Hamiltonian. Here, again, the path splits the plane in two and forces paths
similar to path a. Therefore, such a path as d cannot exist.

In conclusion, the only possible paths are either closed curves crossing the
origin such as path e or divergent curves with a zero Hamiltonian crossing the
origin such as path c. However, the side of the divergent path that goes to
infinity as time increases is unstable. Consider a small perturbation to the
Hamiltonian, $H=\varepsilon\neq 0$, the particle’s trajectory must now be
bound. Thus, it is unstable and will turn to path e. And so, every stable path
crosses the origin, which corresponds to the particles colliding.

This conclusion only applies to the far-field limit. As the particles get
closer, eventually, higher-order terms of the flow will contribute to the
dynamics. If the stream function $\psi$ has no minima or maxima, even in
close-range, then by the same reasoning as above, the particles will collide
because the only allowed paths are still c and e. However, if the stream
function has an extremum in close range, then all path types are allowed.
Particles will still aggregate up to these short distances.

Let us now consider the different paths that the dipolar terms allow. If the
stream function contains a dipolar term of the form $\sin 2\theta$, then paths
that diverge to infinity can have a non-zero Hamiltonian, and so paths like b
are allowed. If the stream function contains the term $\ln r$, no paths can
diverge to infinity, and so paths like d are allowed, while paths like b and c
are not allowed. Examples for both are depicted in Fig. 5. Note that paths
like path a are never allowed in the absence of extrema.

Figure 5: Phase space (configurational space) examples where dipole terms are
allowed showing paths b and d exist in such cases. In blue, $\psi=2\ln
r+\frac{\sin 4\theta}{r^{2}}$ with $H=1.2$ and in red $\psi=\sin
2\theta+\frac{\cos 4\theta}{r^{2}}$ with $H=0.8$. These paths are also stable,
as a perturbation to the Hamiltonian gives a similarly shaped path.

Lastly, we will address the dynamics when $\psi$ is not even. If the stream
function is odd (the far-field limit always has a defined parity), then no
Hamiltonian can be constructed. However, it is easy to see that now the
relative distance, d, is conserved, while the “center of activity” changes
with velocity given by $\textbf{v}(\textbf{d})$. This results in the pair
moving together at a constant velocity, and is true even in close range
interaction. In general, if $\textbf{v}=\textbf{v}_{O}+\textbf{v}_{E}$ which
are even and odd respectively, then

$\begin{cases}\dot{\textbf{r}}_{\text{cm}}=\textbf{v}_{E}(\textbf{d})\\\
\dot{\textbf{d}}=2\textbf{v}_{O}(\textbf{d})\end{cases}$ (15)

and so the relative distance between particles, d, has the same dynamics as
described earlier for an even stream function, as long as $\textbf{v}_{O}\neq
0$ (note that v is odd when $\psi$ is even and vice versa).

##  III many particle dynamics

A system of many particles can no longer be solved analytically. However, for
an even stream function it is possible to construct a Hamiltonian as in Eq.
(5), from which follows the conservation of the center of activity
$\sum_{i}\textbf{r}_{i}/N$ as well as $H$ itself. We integrate over the
system’s equations of motion using the python library scipy.integrate.DOP853,
which is an 8th order Runge-Kutta method. Since the particles are expected to
collide, we also add steric interactions to the system, given by
$\Delta\textbf{v}_{\text{steric}}=k_{\text{steric}}(l_{\text{steric}}-|\textbf{r}_{ij}|)\hat{\textbf{r}}_{ij}$
if $|\textbf{r}_{ij}|<l_{\text{steric}}$ and zero otherwise, where
$\textbf{r}_{ij}$ is the distance from the i-th and j-th particles. The steric
interaction gives a strong repelling force when particles overlap. In our
simulation we used $l_{\text{steric}}=0.5$ and $k_{\text{steric}}=10^{4}$. The
new interactions break the Hamiltonian description, and so $H$ is no longer
conserved. However, the center of activity should remain constant as the
interaction is still symmetric. In addition, between collisions the
Hamiltonian is still conserved (see Fig. 6). The steric interactions are a
highly simplified description of the close-range forces between active
particles that are not described by stokes flow.

Figure 6: The Hamiltonian of a 4-particle system during a collision. A swift
change in the Hamiltonian is observed at time $t=8.5$ when a collision occurs.

We initialize hundreds of randomly positioned particles. Each particle creates
the velocity field discussed in Sec. III. Namely, we take an octopolar flow
with the stream function $\psi=\sin 4\theta/r^{2}$. Since we have established
the similarity of all far-field flows, excluding dipolar flow, as well as
flows created by stream functions with no extrema, we can use this simulation
to infer the behavior of ensembles interaction via such stream functions. The
result of the simulation for different times is displayed in Fig. 1.

Immediately pairs of particles start attracting and eventually collide. This
behavior is a direct result of the two particle dynamics, because the stream
function falls fast with $r$, $\psi\sim r^{-2}$, and so the majority of a
particle’s motion is determined by its closest neighbor, while the weak effect
of other particles forces it to take stable paths only. The colliding
particles quickly align at perfect $45^{\circ}$ lines, corresponding to the
angles where the velocity field is radial only. The stability results from the
steric interactions — once the particles collide, the steric force cancels any
velocity in the radial direction, and the particles slide on each other until
the tangential velocity is also zero. Since the velocity field increases near
the particle, the sliding motion is very fast compared to the global dynamics.
Thus, a collision is characterized by a swift change to the Hamiltonian since
the pair’s contribution to $H$ changed from relatively constant (in time) to
zero. We clearly see this rapid change in Fig. 6, which shows a simulation for
four particles. A collision of two particles is clearly observed by the
decrease in the Hamiltonian at $t=8.5$. In simulations of hundreds of
particles, collisions are too frequent to see the effect clearly.

The two-particle behavior described above repeats on a global scale: all the
particles collide into clusters, and all the clusters collide into bigger
clusters. This is because, at large distances, a cluster functions as a
particle whose strength is the combined strength of its constituents.
Therefore, the dynamics of clusters can be described by a Hamiltonian in Eq.
(5), now with varying strengths. The center of activity is now given by
$\sum_{i}S_{i}\textbf{r}_{i}/\sum_{i}S_{i}$ and is still conserved. In the
case of two-cluster dynamics, the vector
$\textbf{d}=\textbf{r}_{1}-\textbf{r}_{2}$ is sufficient to describe the
dynamics, and its equation of motion is
$\dot{\textbf{d}}=(S_{1}+S_{2})\nabla^{\perp}\psi(\textbf{d})$. Therefore the
two-cluster dynamics also lead to a collision, resulting in the clusters
mimicking the behavior of the particles. The collisions continue until all the
particles form a single stable cluster in the form of an incomplete square
lattice, see Fig. 1. This final result can be deduced from the Hamiltonian,
since, in each collision, the formed cluster must have zero Hamiltonian (i.e.
its self interaction).

The crystallization is thus a general behavior of any ensemble of identical
particles with quenched orientation, as long as the two-particle dynamics lead
to a collision, in the manner discussed in Sec.II.2. Particles will tend to
stabilize in incomplete crystal clusters, with an angle approximately given by
the lowest power term of the stream function, which is dominant at close-range
interactions. More precisely, the stability angles are such that the velocity
field at the particles’ radius ($l_{\text{steric}}/2$) at those angles is
radial only. Of course, this generalized conclusion is true if the system does
not find dynamic stability.

In cases where the stream function is not a pure multipole and possesses
extremum points, the dynamics do not necessarily result in the collision of
two particles. The particles still draw closer up to distances where the next-
order multipole term starts to dominant. Two examples of such cases are given
in the SI. One still results in crystallization, with particles tending to
stay at minimal points of the stream function. The second example is of
squirmers. The particles draw closer, but within the aggregate, there are
still chaotic dynamics.

##  IV Discussion

In this work, we introduced a new method to characterize the flow created by
2D active particles, working directly with a stream function instead of the
velocity field. We used this method to investigate the dynamics of a two-
particle system, which we then used to infer the behavior of a multi-particle
ensemble. We found that in all cases where the stream function does not have a
local extremum and does not have a contribution from the monopole or the
dipole terms, two particles will always collide. This includes all cases where
the flow field is given by a pure even force moment (such as in the far-field
limit). Our simulations suggest that many-body systems of such particles
aggregate into an ordered cluster. The formed clusters appear to be fractal-
like and are somewhat reminiscent of diffusion-limited aggregation [46] but
with two differences — (1) particle motion is driven by active flows, not
diffusion. In fact, the only random component is in the initial positions of
the particles. (2) There is short-range repulsion between particles, not
attraction. Particles stick together due to the balance between steric
interactions and active flows.

The steps we took in this paper are not limited to a simple description of 2D
fluid dynamics, but can be generalized to other effective 2D fluids with more
complex flow equations. Such systems are commonly encountered in biological
membranes [47, 48, 8, 10, 4], or in experimental systems when particles
sediment close to the bottom or float to the interface [14, 49, 50, 51]. Such
cases are associated with Green’s functions that differ from Eq. (2) and flow
equations that differ from Eq. (6). Solutions to their flow equations will
produce all the force multipoles of those systems.

In addition, the arguments and results laid out in Sec.II.2, regarding
possible paths that two particles can take, are also relevant to these more
complex systems, in describing particles whose orientation is frozen. The
requirements for particle dynamics to result in a collision are: (a) The
stream function has no extrema, e.g. by the Hessian test (and twice-
differentiability everywhere except the origin). (b)
$\psi(r\rightarrow\infty)=0$, and (c) there is at least one path diverging to
infinity.

For instance, in a membrane, at small distances, the fluid is conserved in the
membrane, and the flow created by an active protein is given by the dipole
terms in Eq. (9). At large distances, momentum is exchanged with the outer
three-dimensional fluid [47, 6]. An active particle in a membrane at large
distances creates a flow with a stream function that dies out at infinity, as
$\psi=S\frac{\sin 2\theta}{r}$ [10]. The ray $\theta=0$ is a path diverging to
infinity. The Hessian is negative,
$\det\textbf{H}(\psi)=-2S(5+3\cos^{2}2\theta)/r^{6}<0$. Thus, two oriented
active proteins interacting via this stream function are bound to aggregate.

We have pointed out a connection between two-particle dynamics and multi-
particle dynamics. It is tempting to say that as long as $\psi\sim r^{-1}$ or
lower, such that two particles must collide (per our discussion in Sec.II.2),
then larger systems will exhibit the same dynamics because the interactions
are dominated by the nearest neighbors. However, it is unclear if this
assessment is necessarily correct. Based on our findings, we hypothesize that
such a conclusion is correct, but it remains for future work to assess its
validity.

We have focused on interactions that lead to collisions between two particles
— those of pure multipolar flows and cases where there are no extremum points
in the stream function. In other types of dynamics, particles attract up to a
certain distance, but there is no static final state. Ensembles with such
interactions can have chaotic dynamics, but also global structures and order,
as mentioned in the SI. These systems require closer examination.

We have limited our discussion to far-field hydrodynamic interactions in the
low Reynolds limit. However, it is clear that in close range, the interaction
of the particles can be influenced by other effects. For Instance, Yoshinaga
and Liverpool [52] showed the importance of lubrication forces in dense
interaction between active swimmers, which may affect the crystallization and
stability of a many-particle ensembles. Similar effects are expected to change
the final state of our system, but in principle, the characteristics of the
far-range interactions should remain the same, meaning the particles are still
expected to aggregate into large clusters.

Acknowledgments
We thank David B. Stein for his helpful comments.

## References

  * Lauga  _et al._ [2006] E. Lauga, W. R. DiLuzio, G. M. Whitesides, and H. A. Stone, Swimming in circles: motion of bacteria near solid boundaries, Biophysical journal 90, 400 (2006).
  * Lushi  _et al._ [2014] E. Lushi, H. Wioland, and R. E. Goldstein, Fluid flows created by swimming bacteria drive self-organization in confined suspensions, Proceedings of the National Academy of Sciences 111, 9733 (2014).
  * Miles  _et al._ [2019] C. J. Miles, A. A. Evans, M. J. Shelley, and S. E. Spagnolie, Active matter invasion of a viscous fluid: Unstable sheets and a no-flow theorem, Physical Review Letters 122, 098002 (2019).
  * Bagaria and Samanta [2022] S. Bagaria and R. Samanta, Dynamics of force dipoles in curved fluid membranes, Physical Review Fluids 7, 093101 (2022).
  * Robertson [1981] J. D. Robertson, Membrane structure, The Journal of cell biology 91, 189 (1981).
  * Oppenheimer and Diamant [2009] N. Oppenheimer and H. Diamant, Correlated diffusion of membrane proteins and their effect on membrane viscosity, Biophysical Journal 96, 3041 (2009).
  * Lenz  _et al._ [2003] P. Lenz, J.-F. Joanny, F. Jülicher, and J. Prost, Membranes with rotating motors, Physical review letters 91, 108104 (2003).
  * Oppenheimer  _et al._ [2019] N. Oppenheimer, D. B. Stein, and M. J. Shelley, Rotating membrane inclusions crystallize through hydrodynamic and steric interactions, Physical Review Letters 123, 148101 (2019).
  * Oppenheimer  _et al._ [2022] N. Oppenheimer, D. B. Stein, M. Y. B. Zion, and M. J. Shelley, Hyperuniformity and phase enrichment in vortex and rotor assemblies, Nature communications 13, 804 (2022).
  * Manikantan [2020] H. Manikantan, Tunable collective dynamics of active inclusions in viscous membranes, Physical review letters 125, 268101 (2020).
  * Howse  _et al._ [2007] J. R. Howse, R. A. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh, and R. Golestanian, Self-motile colloidal particles: from directed propulsion to random walk, Physical review letters 99, 048102 (2007).
  * Vutukuri  _et al._ [2020] H. R. Vutukuri, M. Lisicki, E. Lauga, and J. Vermant, Light-switchable propulsion of active particles with reversible interactions, Nature communications 11, 2628 (2020).
  * Palacci  _et al._ [2014] J. Palacci, S. Sacanna, S.-H. Kim, G.-R. Yi, D. J. Pine, and P. M. Chaikin, Light-activated self-propelled colloids, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20130372 (2014).
  * Ben Zion  _et al._ [2022] M. Y. Ben Zion, Y. Caba, A. Modin, and P. M. Chaikin, Cooperation in a fluid swarm of fuel-free micro-swimmers, Nature Communications 13, 184 (2022).
  * Modin  _et al._ [2023] A. Modin, M. Y. Ben Zion, and P. M. Chaikin, Hydrodynamic spin-orbit coupling in asynchronous optically driven micro-rotors, Nature Communications 14, 4114 (2023).
  * Dreyfus  _et al._ [2005] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, and J. Bibette, Microscopic artificial swimmers, Nature 437, 862 (2005).
  * Aranson [2013] I. S. Aranson, Active colloids, Physics-Uspekhi 56, 79 (2013).
  * Genkin  _et al._ [2017] M. M. Genkin, A. Sokolov, O. D. Lavrentovich, and I. S. Aranson, Topological defects in a living nematic ensnare swimming bacteria, Physical Review X 7, 011029 (2017).
  * Yoshinaga and Liverpool [2017] N. Yoshinaga and T. B. Liverpool, Hydrodynamic interactions in dense active suspensions: From polar order to dynamical clusters, Physical Review E 96, 020603 (2017).
  * Marchetti  _et al._ [2013] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Reviews of modern physics 85, 1143 (2013).
  * Saintillan and Shelley [2013] D. Saintillan and M. J. Shelley, Active suspensions and their nonlinear models, Comptes Rendus Physique 14, 497 (2013).
  * Rallabandi  _et al._ [2019] B. Rallabandi, F. Yang, and H. A. Stone, Motion of hydrodynamically interacting active particles, arXiv preprint arXiv:1901.04311 (2019).
  * Lushi and Vlahovska [2015] E. Lushi and P. M. Vlahovska, Periodic and chaotic orbits of plane-confined micro-rotors in creeping flows, Journal of Nonlinear Science 25, 1111 (2015).
  * Yeo  _et al._ [2015] K. Yeo, E. Lushi, and P. M. Vlahovska, Collective dynamics in a binary mixture of hydrodynamically coupled microrotors, Physical review letters 114, 188301 (2015).
  * Samanta and Oppenheimer [2021] R. Samanta and N. Oppenheimer, Vortex flows and streamline topology in curved biological membranes, Physics of Fluids 33 (2021).
  * De Graaf and Stenhammar [2017] J. De Graaf and J. Stenhammar, Lattice-boltzmann simulations of microswimmer-tracer interactions, Physical Review E 95, 023302 (2017).
  * Shoham and Oppenheimer [2023] Y. Shoham and N. Oppenheimer, Hamiltonian dynamics and structural states of two-dimensional active particles, Physical Review Letters 131, 178301 (2023).
  * Kos and Ravnik [2018] Ž. Kos and M. Ravnik, Elementary flow field profiles of micro-swimmers in weakly anisotropic nematic fluids: Stokeslet, stresslet, rotlet and source flows, Fluids 3, 15 (2018).
  * Lauga  _et al._ [2021] E. Lauga, T. N. Dang, and T. Ishikawa, Zigzag instability of biased pusher swimmers, Europhysics Letters 133, 44002 (2021).
  * Kim and Karrila [2013] S. Kim and S. J. Karrila,  _Microhydrodynamics: principles and selected applications_ (Courier Corporation, 2013).
  * Spagnolie and Lauga [2012a] S. E. Spagnolie and E. Lauga, Hydrodynamics of self-propulsion near a boundary: predictions and accuracy of far-field approximations, Journal of Fluid Mechanics 700, 105 (2012a).
  * Spagnolie and Lauga [2012b] S. E. Spagnolie and E. Lauga, Hydrodynamics of self-propulsion near a boundary: predictions and accuracy of far-field approximations, Journal of Fluid Mechanics 700, 105 (2012b).
  * Lisicki [2013] M. Lisicki, Four approaches to hydrodynamic green’s functions–the oseen tensors, arXiv preprint arXiv:1312.6231 (2013).
  * Onsager [1949] L. Onsager, Statistical hydrodynamics, Il Nuovo Cimento (1943-1954) 6, 279 (1949).
  * Newton and Platzer [2002] P. K. Newton and M. Platzer, N-vortex problem: Analytical techniques, Appl. Mech. Rev. 55, B15 (2002).
  * Zöttl and Stark [2012] A. Zöttl and H. Stark, Nonlinear dynamics of a microswimmer in poiseuille flow, Physical review letters 108, 218104 (2012).
  * Chajwa  _et al._ [2019a] R. Chajwa, N. Menon, and S. Ramaswamy, Kepler orbits in pairs of disks settling in a viscous fluid, Physical Review Letters 122, 224501 (2019a).
  * Lenz  _et al._ [2004] P. Lenz, J.-F. Joanny, F. Jülicher, and J. Prost, Membranes with rotating motors: Microvortex assemblies, The European Physical Journal E 13, 379 (2004).
  * Chajwa  _et al._ [2019b] R. Chajwa, N. Menon, and S. Ramaswamy, Kepler orbits in pairs of disks settling in a viscous fluid, Physical Review Letters 122, 224501 (2019b).
  * Bolitho and Adhikari [2022] A. Bolitho and R. Adhikari, Overdamped lattice dynamics of sedimenting active cosserat crystals, arXiv preprint arXiv:2202.06102 (2022).
  * Crowdy and Or [2010] D. G. Crowdy and Y. Or, Two-dimensional point singularity model of a low-reynolds-number swimmer near a wall, Physical Review E 81, 036313 (2010).
  * Crowdy and Davis [2013] D. G. Crowdy and A. M. Davis, Stokes flow singularities in a two-dimensional channel: a novel transform approach with application to microswimming, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 469, 20130198 (2013).
  * Lighthill [1952] M. J. Lighthill, On the squirming motion of nearly spherical deformable bodies through liquids at very small reynolds numbers, Communications on pure and applied mathematics 5, 109 (1952).
  * Blake [1971] J. R. Blake, A spherical envelope approach to ciliary propulsion, Journal of Fluid Mechanics 46, 199 (1971).
  * Furi and Martelli [1995] M. Furi and M. Martelli, A multidimensional version of rolle’s theorem, The American mathematical monthly 102, 243 (1995).
  * Witten Jr and Sander [1981] T. A. Witten Jr and L. M. Sander, Diffusion-limited aggregation, a kinetic critical phenomenon, Physical review letters 47, 1400 (1981).
  * Saffman and Delbrück [1975] P. Saffman and M. Delbrück, Brownian motion in biological membranes., Proceedings of the National Academy of Sciences 72, 3111 (1975).
  * Hosaka  _et al._ [2017] Y. Hosaka, K. Yasuda, R. Okamoto, and S. Komura, Lateral diffusion induced by active proteins in a biomembrane, Physical Review E 95, 052407 (2017).
  * Lopez and Lauga [2014] D. Lopez and E. Lauga, Dynamics of swimming bacteria at complex interfaces, Physics of Fluids 26 (2014).
  * Tsang and Kanso [2014] A. C. H. Tsang and E. Kanso, Flagella-induced transitions in the collective behavior of confined microswimmers, Physical Review E 90, 021001 (2014).
  * Wioland  _et al._ [2013] H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler, and R. E. Goldstein, Confinement stabilizes a bacterial suspension into a spiral vortex, Physical review letters 110, 268102 (2013).
  * Yoshinaga and Liverpool [2018] N. Yoshinaga and T. B. Liverpool, From hydrodynamic lubrication to many-body interactions in dense suspensions of active swimmers, The European Physical Journal E 41, 1 (2018).

##  V Supplementary Material

###  V.1 Solving for the Streamfunction

Here we develop the equation for the Fourier components of $\psi$ denoted by
$B_{n}(r)$ from the flow equation given in the main text, and solve it. The
right-hand-side is given by

$\displaystyle\nabla^{\perp}\cdot\textbf{f}$
$\displaystyle=\frac{1}{R}\textbf{f}_{0}\cdot\nabla^{\perp}(\delta(r-R)F(\theta))$
(16) $\displaystyle=\frac{1}{R}\Re\left\\{\eta
e^{i\theta}\left[\frac{1}{r}\delta(r-R)\frac{\partial}{\partial\theta}-i\frac{\partial\delta}{\partial
r}(r-R)\right]F(\theta)\right\\},$

where we used
$\textbf{f}_{0}=f_{0}(\hat{x}\cos\varphi+\hat{y}\sin\varphi)=\Re\left\\{(\hat{r}+i\hat{\theta})\eta
e^{i\theta}\right\\}$. Plugging in the Fourier sums, the left-hand-side is
identified with $-\hat{O}^{2}_{n}B_{n}$, as a result of the Laplacian in polar
coordinates. Thus, the resulting equation is given by

$\hat{O}_{n}^{2}B_{n}(r)=\frac{1}{R}i\eta
C_{n-1}\left[\frac{d}{dr}\delta(r-R)-\frac{n-1}{r}\delta(r-R)\right],$ (17)

We ignore the $\Re$ operator in Eq. (1.1) because the solution is now given by
the real part of $\psi$. To solve this equation, one must first solve the
homogeneous equation $\hat{O}_{n}^{2}B_{n}(r)=0$. This is an Euler equation
whose solutions are always given by monoms. Here, after plugging in
$B_{n}=r^{m}$ we get $m=\pm n,2\pm n$. Thus, the solution to eq. (17) is given
by:

$B_{n}(r)=\frac{1}{R}\begin{cases}x_{n}\left(\frac{r}{R}\right)^{|n|}+y_{n}\left(\frac{r}{R}\right)^{|n|+2}&r<R\\\
z_{n}\left(\frac{r}{R}\right)^{-|n|}+w_{n}\left(\frac{r}{R}\right)^{-|n|+2}&r>R\end{cases},$
(18)

as long as $n\neq 0,\pm 1$ which will be resolved later. The condition on the
coefficients $x_{n},y_{n},z_{n},w_{n}$ are given by the behavior at $r=R$. $B$
and $B^{\prime}$ must be continuous at $r=R$, because their discontinuity
would result in higher derivatives of $\delta(r-R)$ in the equation. To obtain
discontinuity constraints on $B^{\prime\prime}$ and $B^{\prime\prime\prime}$
we first integrate over the equation from $R-\epsilon$ to $R+\epsilon$ and get
one condition, and second, we multiply both sides by $(r-R)$ and integrate
over the same interval, getting the second condition

$\left.rB^{\prime\prime}_{n}(r)\right|_{R-\epsilon}^{R+\epsilon}=i\eta
C_{n-1},\left.r^{2}B^{\prime\prime\prime}_{n}(r)\right|_{R-\epsilon}^{R+\epsilon}=-i\eta
C_{n-1}(n+1).$ (19)

After plugging in Eq. (18), we get equations for the coefficients

$\left(\begin{array}[]{cccc}1&1&-1&-1\\\ |n|&|n|+2&|n|&|n|-2\\\
|n|(|n|-1)&(|n|+1)(|n|+2)&-|n|(|n|+1)&-(|n|-1)(|n|-2)\\\
(|n|-1)(|n|-2)&(|n|+1)(|n|+2)&(|n|+1)(|n|+2)&(|n|-1)(|n|-2)\end{array}\right)\left(\begin{array}[]{c}x_{n}\\\
y_{n}\\\ z_{n}\\\ w_{n}\end{array}\right)=\left(\begin{array}[]{c}0\\\ 0\\\
-1\\\ \frac{n+1}{|n|}\end{array}\right)i\eta C_{n-1}R.$ (20)

The solution is

$\displaystyle z_{n}=\frac{i\eta
C_{n-1}R}{4(1+|n|)}\left(\frac{1}{|n|}+\Theta(n)\right)$ $\displaystyle
w_{n}=\frac{i\eta C_{n-1}R}{4(1-|n|)}\Theta(n),$ (21)

where $\Theta$ is the step function. When $n=-1$ the solution can be gained
from a limit in n, since $w_{n}=0$ for negative $n$. As mentioned in the main
text, we don’t concern ourselves with the case $n=1$, and so $n=0$ is the only
special case. Its homogeneous solutions have a degeneracy, and so the general
homogeneous solution is given as a linear combination of $1,r^{2},\ln
r,r^{2}\ln r$. Applying the same continuity conditions gives the solution in
the main text, where a global constant is ignored.

Figure 7: (a) The velocity field created by the stream function
$\psi=(\frac{R^{2}}{3r^{2}}-\frac{R^{4}}{4r^{4}})\sin(4\theta)+\frac{R^{2}}{2r^{2}}\sin
2\theta$. The red circle marks the particle. The green circle has twice its
radius, noting the closest distance another particle can approach before
steric interactions take effect. The red crosses mark the extremum points. (b)
The final state of a simulation of particles with this stream function. The
particles clearly form lines tilted at the angles of the extrema. In addition,
a few particles collide and form diamonds, similar to the clusters presented
in the main text. It is possible that at even longer times, all the particles
will form a single line.

###  V.2 Constraints on the Force Distribution

As mentioned in the main text, limiting the force distribution to a circle and
to a constant direction does not constrain the solution of the flow given by
the series solution. The coefficients in the series solution are given by:

$D_{ji_{1}i_{2}\dots i_{n}}=\frac{(-1)^{n}}{n!}\int_{S}f_{j}r_{i_{1}}\dots
r_{i_{n}}dA$ (22)

where $S$ is the region the force is applied. The series solution is given by
a sum of the contraction between $\textbf{D}_{n+1}$ and
$\nabla^{n}\textbf{G}$, in the from $D_{ji_{1}\cdots
i_{n}}\partial_{i_{1}}\cdots\partial_{i_{n}}G_{jk}$. Now, we will show that
the possible leading terms in a solution are independent of $S$ in most cases,
and indeed independent in our case, a circle. Because the permutation of
$i_{p},i_{q}$ doesn’t change the component of either tensor, D and
$\nabla^{n}\textbf{G}$, the contraction is of the form:

$\left(\begin{array}[]{c}{n\choose 0}D_{i11\cdots 11}\\\ {n\choose
1}D_{i11\cdots 12}\\\ \vdots\\\ {n\choose n}D{i22\cdots
22}\end{array}\right)\cdot\left(\begin{array}[]{c}\partial_{1}^{n}\partial_{2}^{0}G_{ij}\\\
\partial_{1}^{n-1}\partial_{2}^{1}G_{ij}\\\ \vdots\\\
\partial_{1}^{0}\partial_{2}^{n}G_{ij}\end{array}\right)$ (23)

where in the first line all $i_{p}=1$, in the second line a single index is
equal to $2$ and so on. We now denote $D_{i\alpha}$ the components of
$\textbf{D}_{n}$ with $\alpha$ being the amount of $i_{p}=1$. Since the
derivative of the green function does not depend on $S$, to show the result
space of Eq. 23 does not depend on $S$, it is enough to show that the possible
space of $D_{i\alpha}$ does not depend on $S$ (not a necessary but a
sufficient condition). The components $D_{i\alpha}$ are defined as

$D_{i\alpha}=\frac{(-1)^{n}}{n!}\int_{S}f_{i}x^{\alpha}y^{n-\alpha}dA=\langle
f_{i},b_{\alpha}\rangle_{S},$ (24)

where $\langle\cdot,\cdot\rangle_{S}$ is the standard inner product for
function vector spaces, defined over $S$, and
$b_{\alpha}=\frac{(-1)^{n}}{n!}x^{\alpha}y^{n-\alpha}$. First assume the
direction of the force is fixed, say $f_{2}=0$. As long as
$\mathcal{B}=\\{b_{\alpha}|\alpha=1\cdots n\\}$ are linearly independent under
S, for every $\textbf{D}_{n}$ there exists $f_{1}\in\text{span}\mathcal{B}$
such that Eq. (24) is satisfied. In particular, it is given by
$f_{1}=\sum_{\alpha}D_{1\alpha}b^{*}_{\alpha}$, where $b^{*}_{\alpha}$ are the
dual basis of $\mathcal{B}$. In particular, when $S$ is the unit circle,
$b_{\alpha}=\cos^{\alpha}\theta\sin^{n-\alpha}\theta$ are linearly
independent. Therefore, the possibilities of $\textbf{D}_{n}$ are not limited
by the choice that $S$ is the unit circle, and therefore, the result of the
series solution is not constrained by this choice. An example where the choice
of $S$ does matter is when $\mathcal{B}$ is linearly dependent, for instance,
if $S$ is a section of the line $y=x$.

However, the orientation of the force does matter. The solutions for the force
that give a specific $D_{i\alpha}$ is $f_{i}=\sum D_{i\alpha}b^{*}_{\alpha}$.
Because this force does not have a constant direction (in general), it is not
possible to use the argument above here.

However, a solution for a general force distribution is given by the sum of
two familiar solutions — one considering only the $x$-component of the force,
and second only with the $y$-component. Therefore, the most general force
distribution is given as a sum of the same terms discussed in the main text
($\psi_{nm}$ in Sec. IIIB), just with different coefficients. As a result, the
discussion on the dynamics in the main text applies to any force distribution.
An example is given in the following section.

###  V.3 Constructing Multipoles with Point Forces

Here we highlight how to construct force distributions using point forces
only, such that a specific multipole can be dominant. As mentioned in the main
text, the term $r^{-n}$ ($n>0$) couples to the harmonics $\pm n,n+2$ in the
streamfunction. These harmonics couple to the harmonics $\pm n-1,n+1$ of the
force distribution. Therefore if the lowest (positive) harmonic in a force
distribution is $n+1$ (which is accompanied by $-n-1$), then $r^{-n}$ will be
the leading term in the stream function.

$\psi=\frac{R^{n+1}}{4r^{n}}\Re\left\\{-\frac{i\eta
C_{n+1}}{n+1}e^{i(n+2)\theta}+\frac{i\eta
C_{-n-1}}{n(n+1)}e^{-in\theta}\right\\},$ (25)

where $C_{-n-1}=C_{n+1}^{*}$ because the force distribution is real. In order
to accomplish Eq. (25) with point forces, a force distribution is given by

$F(\theta)=\sum\limits_{k=0}^{2n+1}(-1)^{k}\delta\left(\theta-\theta_{0}+\frac{\pi
k}{n+1}\right).$ (26)

Intuitively, this function has a period of $2\pi/(n+1)$ and so its smallest
harmonic must be $n+1$, and one can verify that using the Fourier expansion of
the delta function, as well as the fact that the roots of
$f(z)=1-z+z^{2}-\cdots-z^{2n+1}$ are all the $2n+2$ roots of unity except
$-1$. This distribution is used in the main text in Sec. IIIA with $n=2$.

As mentioned in the main text, for $n=2$ this distribution creates both second
and forth harmonics in $\psi$. As mentioned earlier, using a force
distribution with a varying direction, it is possible to achieve more varied
linear combination of the terms. For example, we can sum two force
distributions, each in a constant direction, in order to isolate only one of
the harmonics. The distribution in Eq. (26) with $n=2$ gives (as noted in the
main text with $\theta_{0}=0,\textbf{f}_{0}=f_{0}\hat{x},C_{3}=\frac{3}{\pi}$)

$\psi_{\text{oct}}=\frac{R^{3}}{4\pi r^{2}}f_{0}\left(\sin
4\theta+\frac{1}{2}\sin 2\theta\right).$ (27)

Now, to isolate only the forth harmonic, we rotate the force distribution (and
$\textbf{f}_{0}$) by $\pi/2$, creating a distribution with
$\textbf{f}_{0}||\hat{y}$. These two distributions can be summed to create a
distribution with a varying force direction. Summing Eq. (27) with itself
after the transformation $\theta\rightarrow\theta+\frac{\pi}{2}$ gives

$\psi=\frac{R^{3}}{2\pi r^{2}}f_{0}\sin 4\theta.$ (28)

Using this decomposition, any term we have mentioned $\psi_{nm}$ can be
dominant.

###  V.4 Stream Functions with Extrema

As mentioned in the main text, a simple example of a stream function with
minima and maxima is given by a squirmer with
$\psi=(\frac{R^{2}}{r^{2}}-\frac{R^{4}}{r^{4}})\sin(4\theta)$, which has
extrema at $r=\sqrt{2}R$ and
$\theta=\frac{\pi}{8},\frac{3\pi}{8}\cdots\frac{15\pi}{8}$. In close range,
the particles seem to exhibit mostly chaotic dynamics, with weak stability
around the stream function’s extrema.

Another example is given by the stream function
$\psi=(\frac{R^{2}}{3r^{2}}-\frac{R^{4}}{4r^{4}})\sin(4\theta)+\frac{R^{2}}{2r^{2}}\sin
2\theta$, generated by a particle with only a third harmonic in its force
distribution, but with a varying force direction. Its velocity field is
displayed in Fig.7. In a two-particle interaction, it is clear that most paths
that $\mathbf{d}$ takes lead to the $y$-axis, and very few paths (that begin
very close the the particle) result in collision and static stability. Once
$\mathbf{d}$ reaches the $y$-axis, small perturbations will change its paths
into orbits around the extrema points. Thus, particles are likely to be stable
around the extrema in multi-particle interactions.

Indeed, simulating a multi-particle interaction, in a similar manner to the
main text, shows dynamical stability in the formation of rods, angled
according to the angle of the extrema, at $\approx\pm 9^{\circ}$ measured from
the $y$-axis. These structures, while relatively stable, are still much easier
to change than the structures created by attraction of two particles at their
boundary. In fact, while the rods can scatter and destroy each other, the rare
crystallizations in this simulation remain uninterrupted.

These examples show that when extrema are present, and crystallization is not
forced, the dynamics can be very ordered or very chaotic, and a closer look is
needed to understand close-range dynamics in general.