Chiral Edge Modes in Helmholtz-Onsager Vortex Systems

The point vortex model consists of a system of $N$ point vortices in a 2D incompressible, inviscid fluid with constant density. Each vortex has strength $\lambda_{a}$ and position $\mathbf{x}_{a}(t)=(x_{a}(t),y_{a}(t))$ in a domain $\Omega$, which corresponds to the following vorticity field

$\displaystyle\omega=\sum_{a}\lambda_{a}\delta\left(\mathbf{x}-\mathbf{x}_{a}(t)\right)$

The flow field $\mathbf{u}(t,\mathbf{x})$ for this vortex configuration is obtained from the streamfunction $\psi(t,\mathbf{x})$ as [14]

$\displaystyle\mathbf{u}=\nabla\psi\times\mathbf{e}_{z},\qquad\qquad\nabla^{2}\psi=-\omega$

The streamfunction can be written in terms of Green’s functions

$\displaystyle\psi=-\sum_{a=1}^{N}\lambda_{a}G(\mathbf{x},\mathbf{x}_{a}(t))$

where the Green’s function $G$ satisfies

$\displaystyle\nabla^{2}G(\mathbf{x},\mathbf{y})$

$\displaystyle=\delta(\mathbf{x}-\mathbf{y}),\qquad\qquad G|_{\partial\Omega}=0$

The boundary condition for the Green’s function arises from the fact that the boundary of the domain must be a streamline, $\psi|_{\partial\Omega}=\text{const.}$ For simply connected domains, we can set this constant to 0. Using this streamfunction as an ansatz for the inviscid vorticity equation [1, 31, 9]

$$\frac{\partial\omega}{\partial t}+\mathbf{u}\cdot\nabla\omega=0$$

yields a Hamiltonian system where $x_{a}$ and $y_{a}$ are canonically conjugate

$\displaystyle\lambda_{a}\dot{x}_{a}=\frac{\partial\mathcal{H}}{\partial y_{a}},\qquad\qquad\lambda_{a}\dot{y}_{a}=-\frac{\partial\mathcal{H}}{\partial x_{a}}$

(1)

The Hamiltonian is given by

$\displaystyle\mathcal{H}$

$\displaystyle=-\sum_{a<b}\lambda_{a}\lambda_{b}G(\mathbf{x}_{a},\mathbf{x}_{b})-\frac{1}{2}\sum_{a}\lambda_{a}^{2}g(\mathbf{x}_{a})$

(2)

where

$$g(\mathbf{x}_{a})=\lim_{\mathbf{x}\rightarrow\mathbf{x}_{a}}\left(G(\mathbf{x},\mathbf{x}_{a})-\frac{1}{2\pi}\log|\mathbf{x}-\mathbf{x}_{a}|\right)$$

is the regularized Green’s function. The Hamiltonian coincides with the regularized kinetic energy of the fluid [32]

$\displaystyle\mathcal{K}=\frac{1}{2}\int_{\Omega}d^{2}\mathbf{x}\,|\mathbf{u}|^{2}=\frac{1}{2}\int_{\Omega}d^{2}\mathbf{x}\,\nabla\psi\cdot\nabla\psi=\frac{1}{2}\int_{\Omega}d^{2}\mathbf{x}\,\psi\omega=\mathcal{H}$

The quantities defined only depend on length $L$ and time $T$: $[\lambda]=L^{2}/T,\;[\mathcal{H}]=L^{4}/T^{2}$. Choosing a reference length scale $L_{0}$, we can define a timescale $T_{0}=L_{0}^{2}/(N|\lambda|)$. Since the fluid density plays no role in the system, this also gives energy scale $E_{0}=L_{0}^{4}/T_{0}^{2}$ and angular momentum scale $J_{0}=L_{0}^{4}/T_{0}$. In particular

$$\frac{E_{0}}{N^{2}}=|\lambda|^{2},\qquad\qquad\frac{J_{0}}{NL_{0}^{2}}=|\lambda|$$

In the remainder, we consider systems of identical vortices with positive circulation,

$$\lambda_{a}\equiv\lambda>0$$

We choose units so that $L_{0}=1$ and $\lambda=1$. This fixes the timescale, $T_{0}=1/N$, reflecting the fact that the fluid rotates faster as the number of point vortices increases. Below, we will focus on two different domains, the unit disk $D$ with radius $L_{0}=1$ and a bean-shaped domain $B$ obtained from $D$ by a conformal mapping specified below.

In complex coordinates, $z=x+iy$, the Green’s functions $G$ and $g$ for the disk $D$ are [33, 34, 35]

$\displaystyle G^{D}(z,w)$

$\displaystyle=\frac{1}{2\pi}\log\left|\frac{z-w}{z^{*}w-1}\right|\qquad\qquad g^{D}(z)=-\frac{1}{2\pi}\log|1-z^{*}z|$

This gives a compact expression for the Hamiltonian

$\displaystyle\mathcal{H}$

$\displaystyle=-\frac{1}{2\pi}\sum_{a<b}\lambda^{2}\log\left|\frac{z_{a}-z_{b}}{z_{a}^{*}z_{b}-1}\right|+\frac{1}{4\pi}\sum_{a}\lambda^{2}\log|1-z_{a}^{*}z_{a}|$

and the velocity

$\displaystyle\lambda\dot{z}_{a}=-2i\frac{\partial\mathcal{H}}{\partial\bar{z}_{a}}$

$\displaystyle=\frac{i\lambda^{2}}{2\pi}\sum_{b\neq a}\frac{1}{z_{a}^{*}-z_{b}^{*}}-\frac{i\lambda^{2}}{2\pi}\sum_{b}\frac{1}{z_{a}^{*}-(1/z_{b})}$

The first sum captures the flow generated by point vortices at the locations $z_{b}$, in the the absence of boundaries. Analogously, the second sum, which includes $b=a$, is due to image vortices located at the points $1/z_{b}^{*}$ [Fig. 1(a)]. The velocity of the whole fluid, $\mathbf{u}(t,z)=(u(t,z),v(t,z))$, is

$\displaystyle u(t,z)+iv(t,z)=\frac{i\lambda}{2\pi}\sum_{a}\left[\frac{1}{z^{*}-z_{a}^{*}}-\frac{1}{z^{*}-(1/z_{a})}\right]$

The Hamiltonian for the unit disk can be split up into a bulk energy and a boundary energy, arising from the interaction between real vortices and image vortices [Figs. 1(a) and 1(b)]

$\displaystyle\mathcal{H}$

$\displaystyle=-\frac{1}{4\pi}\sum_{a\neq b}\lambda^{2}\log\left|z_{a}-z_{b}\right|+\frac{1}{4\pi}\sum_{a,b}\lambda^{2}\log|1-z_{a}z_{b}^{*}|$

The singularities of $\mathcal{H}$ correspond to the locations of these real and image vortices. The energy is singular when the $a$’th vortex approaches the real vortex at $z_{b}$ or the image vortex at $1/z_{b}^{*}$ [Figs. 1(a) and 1(b)].

A bean-shaped domain $B$ can be obtained from the unit disk by a conformal mapping. Consider the image of $D$ under a map of the form $f^{-1}(z)=(b^{\prime}z-1)^{2}/a^{\prime}$. When $b^{\prime}=1$, this mapping produces a nonconvex cardioid domain, which posseses a cusp singularity at the boundary. By taking $b^{\prime}=0.9$, we instead obtain a smooth nonconvex bean-shaped domain, $B$ [Fig. 1(c)]. We set $a^{\prime 2}=2b^{\prime 2}(b^{\prime 2}+2)$ so that $B$ has area $\pi$, equal to the unit disk. The Green’s function in $B$ transforms with the conformal map $f$

$\displaystyle G^{B}(z,w)$

$\displaystyle=G^{D}(f(z),f(w))=\frac{1}{2\pi}\log\left|\frac{f(z)-f(w)}{f(z)^{*}f(w)-1}\right|$

and the regularized Green’s function is

	$\displaystyle g^{B}(z)$	$\displaystyle=\lim_{w\rightarrow z}\left(\frac{1}{2\pi}\log\left|\frac{f(z)-f(w)}{f(z)^{*}f(w)-1}\right|-\frac{1}{2\pi}\log|z-w|\right)$
		$\displaystyle=-\frac{1}{2\pi}\left(\log|1-f(z)^{*}f(z)|-\log|f^{\prime}(z)|\right)$

As before, the form of the Hamiltonian indicates the presence of image vortices

$\displaystyle\mathcal{H}$

$\displaystyle=-\frac{1}{4\pi}\sum_{a\neq b}\lambda^{2}\log\left|f(z_{a})-f(z_{b})\right|+\frac{1}{4\pi}\sum_{a,b}\lambda^{2}\log|1-f(z_{a})f(z_{b})^{*}|-\frac{1}{4\pi}\sum_{a}\lambda^{2}\log|f^{\prime}(z_{a})|$

(3)

The positions of real and image vortices again follow from the singularities of the Hamiltonian. The first term is singular when $z_{a}$ approaches the real vortex at $z_{b}$, and the second term is singular when $z_{a}$ approaches the image vortex at $f^{-1}\left(1/f(z_{b})^{*}\right)$ [Fig. 1(c)]. Furthermore, the first and third terms of the Hamiltonian in (3) are bounded below, unlike the second term. To see this, observe that the boundedness of the domain means that the expressions $-\log|f(z_{a})-f(z_{b})|$ in the first term are bounded below. Similarly, $-\log|f^{\prime}(z_{a})|$ is bounded below whenever $f^{\prime}(z_{a})$ is bounded, as is the case for the mappings we consider. On the other hand, the second term of (3) approaches negative infinity when real vortices approach image vortices, although this divergence is not independent of divergences in the first term. Nevertheless, this heuristic argument indicates that an effective image vortex attraction at low energies will give rise to vortex clustering at the boundary of the domain. This image vortex dipole mechanism (Fig. 1) is responsible for the formation of edge modes.

In addition to the energy $\mathcal{H}$, an incompressible, inviscid fluid in 2D also has conservation laws for linear momentum $\mathcal{P}$, angular momentum $\mathcal{A}$, vorticity $\mathcal{V}$, and enstrophy $\mathcal{E}$

$\displaystyle\mathcal{P}=\int_{\Omega}d^{2}\mathbf{x}\,\mathbf{u},\qquad\mathcal{A}=\int_{\Omega}d^{2}\mathbf{x}\,\left(\mathbf{x}\wedge\mathbf{u}\right),\qquad\mathcal{V}=\int_{\Omega}d^{2}\mathbf{x}\,\omega,\qquad\mathcal{E}=\int_{\Omega}d^{2}\mathbf{x}\,\omega^{2}$

where the pseudoscalar $\mathbf{x}\wedge\mathbf{u}=xv-yu$ corresponds to the $z$-component of the 3D cross-product. The linear momentum vanishes in a bounded domain, which can be shown componentwise using the streamfunction $\mathbf{u}=(\psi_{y},-\psi_{x})$

$\displaystyle\mathcal{P}\cdot\mathbf{e}_{x}=\int_{\Omega}d^{2}\mathbf{x}\,\nabla\cdot\left(\psi\mathbf{e}_{y}\right)=\int_{\partial\Omega}d\mathbf{n}\cdot\mathbf{e}_{y}\psi=0$

where $d\mathbf{n}$ is the normal curve element. The second equality follows from the divergence theorem and the final equality uses the fact that $\psi|_{\partial\Omega}$ is constant and $\partial\Omega$ is a closed curve. The vorticity and enstrophy have straightforward expressions in terms of the vortex strengths, however the enstrophy is singular

$\displaystyle\mathcal{V}=\sum_{a}\lambda_{a}=N\lambda,\qquad\qquad\mathcal{E}=\sum_{a,b}\lambda_{a}\lambda_{b}\delta(\mathbf{x}_{a}-\mathbf{x}_{b})=\lambda^{2}\sum_{a,b}\delta(\mathbf{x}_{a}-\mathbf{x}_{b})$

In the systems of identical point vortices we consider here, the conservation of $\mathcal{V}$ is equivalent to the conservation of particle number, $\mathcal{V}=N$. The conversation law for angular momentum depends on the pressure field $p(t,\mathbf{x})$, which can be found in terms of the velocity field, $\mathbf{u}$, through the Euler equation for an inviscid fluid

$\displaystyle-\nabla p=\frac{\partial\mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}$

Taking the cross product of the Euler equation with $\mathbf{x}$ gives an equation for angular momentum transport; integrating this transport equation and using the divergence theorem yields the conservation law for angular momentum in terms of the boundary advection of angular momentum and the boundary torque due to pressure

$\displaystyle\frac{d}{dt}\mathcal{A}$

$\displaystyle=-\int_{\partial\Omega}d\mathbf{n}\cdot\mathbf{u}\,\left(\mathbf{x}\wedge\mathbf{u}\right)+\int_{\partial\Omega}d\mathbf{n}\wedge\mathbf{x}\,p=\int_{\partial\Omega}d\mathbf{n}\wedge\mathbf{x}\,p$

The advection term vanishes since $\mathbf{u}$ has no normal component on the boundary. Angular momentum is globally conserved if the boundary torque vanishes, which occurs when the domain $\Omega$ has the appropriate symmetry [14, 24]. Using Stokes’ theorem, the angular momentum of the fluid can be expressed as

$\displaystyle\mathcal{A}=\frac{1}{2}\int_{\Omega}d^{2}\mathbf{x}\,\left[\nabla\wedge\left(\mathbf{u}|\mathbf{x}|^{2}\right)-\omega|\mathbf{x}|^{2}\right]=\frac{1}{2}\int_{\partial\Omega}d\mathbf{r}\cdot\mathbf{u}\,|\mathbf{x}|^{2}-\frac{1}{2}\mathcal{J}$

where $d\mathbf{r}$ is the tangential curve element and

$$\mathcal{J}=\sum_{a}\lambda|\mathbf{x}_{a}|^{2}$$

If the domain has circular symmetry, $\mathcal{A}$ is conserved. In particular, when $\Omega$ is the unit disk $D$, the circulation term can be simplified

$\displaystyle\mathcal{A}=\frac{1}{2}\sum_{a}\lambda_{a}-\frac{1}{2}\sum_{a}\lambda_{a}|\mathbf{x}_{a}|^{2}=\frac{1}{2}\mathcal{V}-\frac{1}{2}\mathcal{J}$

where $\mathcal{J}$ is the conserved vortex angular momentum. The statistical mechanics of the system can be described in terms of the nonzero and nonsingular conserved quantities $\mathcal{H},\mathcal{A}$ and $\mathcal{V}$.

The different statistical regimes of the point vortex system can be identified through the density of states. For the disk $D$, conservation of both $\mathcal{H}$ and $\mathcal{J}$ implies there is a joint density of states $W(E,J)$ [Fig. 2(a)]. Defining $\xi=(x_{1},y_{1},x_{2},y_{2},...,x_{N},y_{N})$, we have

$\displaystyle W(E,J)$

$\displaystyle=\frac{1}{\pi^{N}}\int_{D^{N}}d\xi\,\delta\left(\mathcal{H}(\xi)-E\right)\delta\left(\mathcal{J}(\xi)-J\right)$

The marginal densities $W_{\mathcal{H}}(E)$ and $W_{\mathcal{J}}(J)$ over energy and angular momentum follow from the joint density by [Fig. 2(b) and 2(c)]

	$\displaystyle W_{\mathcal{H}}(E)$	$\displaystyle=\frac{1}{\pi^{N}}\int_{D^{N}}d\xi\,\delta\left(\mathcal{H}(\xi)-E\right)$
	$\displaystyle W_{\mathcal{J}}(J)$	$\displaystyle=\frac{1}{\pi^{N}}\int_{D^{N}}d\xi\,\delta\left(\mathcal{J}(\xi)-J\right)$

Flows where $\mathcal{J}$ is large in magnitude must involve a concentration of vortices at the boundary, and thus $\mathcal{J}$-conservation may be thought of as an angular momentum barrier. However, due to the correlation between $\mathcal{H}$ and $\mathcal{J}$ [Fig. 1(a)], it is hard to disentangle this effect from the low-energy boundary attraction described above.

As is typical of a generic bounded domain, the bean-shaped domain has no additional conserved quantities, so there is only a density of states over energy [Fig. 2(d)]. This density, $W_{\mathcal{H}}(E)$, has a maximum at some energy, $E=E^{*}$, which arises from the fact that the phase space is bounded [Figs. 2(b) and 2(d)]. Here, we focus on the low-energy regime, $E<E^{*}$.

To identify the low-energy regime for large $N$ and understand the corresponding statistical states, we return to the density of states for a general bounded domain

$\displaystyle W_{\mathcal{H}}(E)$

$\displaystyle=\frac{1}{\pi^{N}}\int_{\Omega^{N}}d\xi\,\delta\left(\mathcal{H}(\xi)-E\right)$

The function $W_{\mathcal{H}}(E)$ can also be understood as the probability density function (pdf) of the energy when $N$ points are chosen uniformly at random in the domain. Let $Z_{a}=X_{a}+iY_{a}$ be a collection of such uniform random variables for $a=1$ to $N$. The energy is now a random variable $\hat{H}$

$\displaystyle\hat{H}=-\frac{1}{4\pi}\sum_{a\neq b}G(Z_{a},Z_{b})+\frac{1}{4\pi}\sum_{a}g(Z_{a})$

where the first sum contains $N(N-1)$ terms, the second sum contains $N$ terms and the strength of each vortex $\lambda$ has been set to 1 as before. Previous work [36] has identified the limiting distribution of $\hat{H}$ in bounded domains. For the identical vortex systems considered here, this limiting distribution is Gaussian. An informal derivation of this fact follows from the theory of $U$-statistics [37, 38]. Set $\bar{H}=\hat{H}/N^{2}$, and consider

$$\sqrt{N}(\bar{H}-\mu)=\sqrt{N}\left(-\frac{1}{4\pi N^{2}}\sum_{a\neq b}G(Z_{a},Z_{b})-\mu\right)+\frac{1}{4\pi N\sqrt{N}}\sum_{a}g(Z_{a})$$

(4)

where

$$\mu=-\frac{1}{4\pi}\mathbb{E}\left[G(Z_{a},Z_{b})\right]$$

The second term in (4) is a sum of only $N$ independent random variables, and thus converges to $0$. A central limit type theorem [37] for the sequence of $N(N-1)$ random variables $G(Z_{a},Z_{b})$ shows that $\sqrt{N}(\bar{H}-\mu)$ is asymptotically normal, where $\mu$ and the variance $\sigma^{2}$ are given by [37, 39, 36]

$\displaystyle\mu=-\frac{1}{4\pi}\mathbb{E}G(Z_{a},Z_{b}),\qquad\qquad\sigma^{2}=\lim_{N\rightarrow\infty}N\,\text{Var}(\bar{H})=\lim_{N\rightarrow\infty}\frac{\text{Var}(\hat{H})}{N^{3}}$

Furthermore, $U$-statistics theory [37] indicates that $\sigma^{2}>0$ holds whenever $\mathbb{E}(G(z,Z_{a}))$ is not constant as a function of $z$. Formalizing the above argument would require showing that $G$ and $g$ meet suitable regularity conditions [37]. The scaling of the variance is $\text{Var}(\bar{H})\sim O(1/N)$, analogous to the central limit theorem. This scaling may be understood intuitively as a consequence of the fact that there are only $N$, and not $N^{2}$, independent points. Although $\mathbb{E}(\hat{H}/N^{2})\sim O(1)$, as $N\rightarrow\infty$, there is an $O(1/N)$ contribution to the mean coming from the boundary terms, $\mathbb{E}(g(Z_{a}))$. For finite $N$, the pdf of $\hat{H}/N^{2}$ can therefore be better approximated by a Gaussian with mean

$$\mathbb{E}\left(\frac{\hat{H}}{N^{2}}\right)=\mu+\frac{\mu_{b}}{N}$$

where $\mu_{b}=\mathbb{E}(g(Z_{a}))/4\pi$. The limit allows us to define the low-energy regime in terms of a rescaled energy $\epsilon$, which measures the distance from the mean energy in units of standard deviation

$\displaystyle\epsilon=\frac{\sqrt{N}}{\sigma}\left[\frac{\mathcal{H}}{N^{2}}-\left(\mu+\frac{\mu_{b}}{N}\right)\right]$

(5)

Similarly, for the disk, the angular momentum becomes a random variable $\hat{J}$

$\displaystyle\hat{J}=\sum_{a}|Z_{a}|^{2}$

where $\mathbb{E}(|Z|^{2})=1/2$ and $\text{Var}(|Z|^{2})=1/12$. Using the central limit theorem, the appropriately rescaled angular momentum is therefore

$\displaystyle\ell=\sqrt{12N}\left(\frac{\hat{J}}{N}-\frac{1}{2}\right)$

Statistical edge modes persist in the disk at low energy for large vortex numbers $N$. When $N$ is sufficiently large, the densities of state approach their limiting values [Figs. 5(a) and 5(b)], and the rescaled energy $\epsilon$ defines the low-energy regime. The robustness of edge modes at low $\epsilon$ [Figs. 5(c) and 5(d)], can be understood analytically through the mean-field approach of maximizing the entropy of the vortex system at fixed energy, particle number, and angular momentum if the domain has circular symmetry [40, 41, 32]. It is convenient to explicitly state factors of the vortex circulation $\lambda$ in the mean-field picture. We set $\rho=\omega L_{0}^{2}/(N\lambda)$ to be the dimensionless vortex density, so $\rho$ integrates to $L_{0}^{2}=1$. Similarly, we rescale the streamfunction, $\tilde{\psi}=\psi L_{0}^{2}/N$ so that $\nabla^{2}\tilde{\psi}=-\lambda\rho$. The constrained entropy functional is [21]

	$\displaystyle\mathcal{I}[\rho]$	$\displaystyle=S-\beta\lambda^{2}T_{0}^{2}\mathcal{H}-\alpha\int_{\Omega}d^{2}\mathbf{x}\,\rho-\frac{\beta\gamma L_{0}^{2}}{N}T_{0}\mathcal{J}$		(6)
		$\displaystyle=-\int_{\Omega}d^{2}\mathbf{x}\,\rho\log\rho-\beta\lambda\int_{\Omega}d^{2}\mathbf{x}\,\tilde{\psi}\rho-\alpha\int_{\Omega}d^{2}\mathbf{x}\,\rho-\beta\gamma\lambda\int_{\Omega}d^{2}\mathbf{x}\,r^{2}\rho$

where $\beta$, $\alpha$, $\gamma$ are Lagrange multipliers enforcing the constraints of constant energy, particle number, and angular momentum, respectively, with units $[\beta]=T^{2}/L^{6}$, $[\gamma]=L^{2}/T,[\alpha]=0$. If $\Omega$ does not have circular symmetry, then $\gamma=0$. By analogy with the first law of thermodynamics, $\beta$ and $\alpha$ are sometimes thought of as an inverse temperature and a chemical potential. Extremizing $\mathcal{I}$ yields

$\displaystyle\frac{\delta\mathcal{I}}{\delta\rho}$

$\displaystyle=-\log\rho-1-\beta\lambda\tilde{\psi}-\beta\lambda\gamma r^{2}-\alpha=0$

and thus

$\displaystyle\rho=e^{-1-\alpha}e^{-\beta\lambda\left(\tilde{\psi}+\gamma r^{2}\right)}=\frac{1}{Z}e^{-\beta\lambda\left(\tilde{\psi}+\gamma r^{2}\right)}$

where $Z=e^{1+\alpha}$. Using the streamfunction-vorticity relation, $\nabla^{2}\tilde{\psi}=-\lambda\rho$, gives the mean-field equation

$\displaystyle\nabla^{2}\tilde{\psi}=-\frac{\lambda}{Z}e^{-\beta\lambda\left(\tilde{\psi}+\gamma r^{2}\right)},\qquad\qquad\tilde{\psi}|_{\partial\Omega}=0$

(7)

where

$\displaystyle Z=\int_{\Omega}d^{2}\mathbf{x}\,e^{-\beta\lambda(\tilde{\psi}+\gamma r^{2})}$

(8)

By setting $\phi=-\beta\lambda\tilde{\psi}$ and $\beta_{1}=\beta/Z$, the mean-field equation can be expressed as a nonlinear eigenvalue problem

$\displaystyle\nabla^{2}\phi$

$\displaystyle=\beta_{1}\lambda^{2}e^{-\beta\gamma\lambda r^{2}}e^{\phi},\qquad\qquad\phi|_{\partial\Omega}=0$

where

$\displaystyle\beta$

$\displaystyle=\beta_{1}\lambda^{2}\int_{\Omega}d^{2}\mathbf{x}\,e^{-\beta\lambda\gamma r^{2}}e^{\phi}$

In the disk, we can simplify this equation by assuming axisymmetry, $\phi=\phi(r)$, and neglecting angular momentum, $\gamma=0$, based on the strong correlation between angular momentum and energy [Fig. 5(a)]

$\displaystyle\phi^{\prime\prime}+\frac{1}{r}\phi^{\prime}-\beta_{1}\lambda^{2}e^{\phi}=0,\qquad\phi^{\prime}(0)=0,\quad\phi(1)=0$

The boundary condition at 0 comes from the requirement that $\phi$ is smooth at the origin. This equation can be solved exactly to give a vortex density $\rho$ with one free parameter [32]

$\displaystyle\rho(r)=\frac{8}{8\pi+\beta}\frac{1}{\left(1-\frac{\beta r^{2}}{8\pi+\beta}\right)^{2}},\qquad\beta>-8\pi$

(9)

This solution agrees quantitatively with simulations of Hamiltonian dynamics [Figs. 5(e) and 5(f)], confirming the validity of the mean-field theory at low energies. Furthermore, the mean-field solution predicts a wide range of low-energy edge modes. Whenever $\beta>0$, the solution (9) is maximized on the boundary. For large $\beta$, the solution describes an increasingly concentrated edge mode. The role of image vortices in this edge mode solution can be visualized by mapping the system onto a sphere [Figs. 1 and 6].

Operators of the form $\nabla^{2}-m^{2}$ have a topological interpretation in certain contexts, including in the theory of chiral hydrodynamics [42]. In our system, this operator appears in the linearization of the mean-field equations (7) after neglecting angular momentum ($\gamma=0$)

$\displaystyle\left(\nabla^{2}-\frac{\beta}{Z}\lambda^{2}\right)\tilde{\psi}=-\frac{\lambda}{Z}$

(10)

This linearization is valid when $|\beta\lambda\tilde{\psi}|$ is small, which holds close to the boundary, see Eq. (7). Although this equation has exponentially decaying solutions, the assumption that $|\beta\lambda\tilde{\psi}|$ is small means that these solutions are not valid everywhere. The true decay behavior is given by exact solutions to the full nonlinear equation, such as (9). Setting $\psi_{0}=\tilde{\psi}-1/\beta\lambda$ removes the constant term in equation (10)

$\displaystyle\left(\nabla^{2}-\beta_{1}\lambda^{2}\right)\psi_{0}=0$

(11)

where $\beta_{1}=\beta/Z$ as before.

The topological interpretation of this simplified equation comes from taking square root of the differential operator

$\displaystyle M=\begin{pmatrix}-i\partial_{y}&&-\partial_{x}+\lambda\sqrt{\beta_{1}}\\ \partial_{x}+\lambda\sqrt{\beta_{1}}&&i\partial_{y}\end{pmatrix},\qquad M^{2}=-\left(\nabla^{2}-\beta_{1}\lambda^{2}\right)I_{2}$

where $I_{2}$ is the $2\times 2$ identity matrix. The matrix $M$ is the 2D Dirac Hamiltonian, and has topological properties [42]. In the Fourier representation, $-i\partial_{x}\mapsto q_{x}$ and $-i\partial_{y}\mapsto q_{y}$, we can write the matrix operator $M$ as

$\displaystyle M(q_{x},q_{y},\lambda)=\begin{pmatrix}q_{y}&&-iq_{x}+\lambda\sqrt{\beta_{1}}\\ iq_{x}+\lambda\sqrt{\beta_{1}}&&-q_{y}\end{pmatrix}$

(12)

When viewed as a function on the parameter space $(q_{x},q_{y},\lambda)$, a topological quantity known as a Chern number [42, 43] can be associated with $M$ (Appendix). This invariant measures the winding of the eigenvectors of $M$ across the $(q_{x},q_{y})$ plane, and depends only on the sign of $\lambda$. As a consequence, regions in the physical $(x,y)$ space where $\lambda$ changes sign are special, and host localized steady-state flows [42] (Appendix). This corresponds exactly to our physical picture of vortex dynamics. Our systems of positive vortices have $\lambda=1$, and are surrounded by negative image vortices with $\lambda=-1$ (Fig. 6). At the boundary, where $\lambda$ changes sign, we observe statistical edge modes. This argument requires $\beta<0$, and so is not incompatible with the high energy regime, where edge modes are not observed. Although this description [42, 44] is not the standard picture of topological edge modes [43], it presents a possible avenue for exploring connections with other topological phenomena.

Low-energy edge modes also survive in nonconvex domains at large $N$. As in the disk, the rescaled energy $\epsilon$ defines low-energy states [Fig. 7(a)] which display statistical edge modes [Figs. 7(b) and 7(c)]. Expressed as a nonlinear eigenvalue problem, the mean-field equation without angular momentum is

$\displaystyle\nabla^{2}\phi=\beta_{1}\lambda^{2}e^{\phi},\qquad\qquad\phi|_{\partial\Omega}=0$

(13)

where $\rho$ is obtained as

$\displaystyle\rho=\frac{1}{\beta}\nabla^{2}\phi,\qquad\qquad\beta=\beta_{1}\lambda^{2}\int_{\Omega}d^{2}\mathbf{x}\,e^{\phi}$

The numerical solution of this equation agrees with time-averaged simulations of the Hamiltonian dynamics [Fig. 7(c)-7(f)]. Mean-field theory thus remains valid for irregular, nonconvex domains at low energy, and produces statistical states with edge modes. In particular, any smooth solution of (13) with $\beta_{1}>0$ gives a vortex density which is maximized on the boundary. To see this, observe that $\phi$ is subharmonic, $\nabla^{2}\phi>0$. By the maximum principle for subharmonic functions [32], $\phi$ must be maximized on the boundary of its domain. Since $\phi$ is constant on the boundary, we have

$$\phi(\mathbf{x}_{b})\geq\phi(\mathbf{x})\qquad\forall\mathbf{x}_{b}\in\partial\Omega,\;\mathbf{x}\in\Omega$$

The vortex density $\rho$ is an increasing function of $\phi$,

$\displaystyle\rho=\frac{1}{\beta}\nabla^{2}\phi=\frac{\beta_{1}\lambda^{2}}{\beta}e^{\phi}$

Thus $\rho$ is also maximized on the boundary.

In the systems of positive point vortices considered here, time-averaged dynamical simulations are found to agree with vortex distributions obtained by Monte Carlo sampling at fixed energy, even for large $(N=500)$ vortex numbers (Fig. 8). Although there are mathematical questions which remain unresolved, it is believed that this observed ergodicity is an appropriate assumption for the point vortex model under certain conditions [55, 21]. In other cases, the system is known to be non-ergodic [56, 57]. Numerical studies indicate evidence for ergodicity at intermediate energies for mixed vortex systems (with $\lambda_{a}=\pm 1$) in bounded [26, 36] and periodic [36] domains. On the sphere, the ergodic hypothesis appears to hold for a range of energies [58]. At very low energies, however, dipole formation in mixed vortex systems can produce quasi-stationary states [59, 36]. Furthermore, it is possible that the relaxation time scales unfavorably with the vortex number $N$ in certain scenarios [21]. Since we only consider positive vortices, we do not expect dipole formation to obstruct equilibration at very low energies. However, relaxation timescale issues could be responsible for the discrepancy we observe in the vortex distributions obtained from dynamical simulations, Monte Carlo sampling and mean-field theory, near the center ($r\to 0$) of the disk [Fig. 8(c)].

The edge mode obtained for $N=500$ vortices at rescaled energy $\epsilon=-2.5$ (Fig. 8) is shallower than that obtained for $(N,\epsilon)=(80,-2.5)$ (Figs. 5 and 7). In particular, the mean-field solution for the $(N,\epsilon)=(80,-2.5)$ edge mode has $\beta=15$, wheras the $(N,\epsilon)=(500,-2.5)$ edge mode has $\beta=5.2$. However, in both cases, the mean-field solution (9) accurately describes the vortex density

$\displaystyle\rho(r)=\frac{8}{8\pi+\beta}\frac{1}{\left(1-\frac{\beta r^{2}}{8\pi+\beta}\right)^{2}},\qquad\beta>-8\pi$

This suggests that, by choosing an $\epsilon$ which corresponds to a larger $\beta$, a sharper edge mode can be found for $N=500$.

The mean-field equations (7) can be derived using the canonical ensemble. Following the presentation in Ref. [60], the canonical partition function for $N$ vortices in a domain $\Omega$ is

$\displaystyle Z(\beta)=\int_{\Omega^{N}}d\xi\,e^{-\beta\mathcal{H}(\xi)}$

where $\xi=(x_{1},y_{1},...,x_{N},y_{N})$ is a point in $\Omega^{N}$. This can be rewritten in terms of the dimensionless vortex density $\rho=\omega T_{0}$. Let $\Gamma[\rho]$ be the functional corresponding to the total vortex density

$\displaystyle\Gamma[\rho]=\int_{\Omega}d^{2}\mathbf{x}\,\rho(\mathbf{x})=1$

Setting $W[\rho]$ to be the number of microstates $\{\xi\}$ corresponding to the macrostate $\{\rho(\mathbf{x})\}$, allows us to write $Z(\beta)$ as an integral over fields $\rho(\mathbf{x})$

$\displaystyle Z(\beta)=\int_{\Omega^{N}}d\xi\,e^{-\beta\mathcal{H}(\xi)}=\int\mathcal{D}\rho\,W[\rho]\,\delta\left(\Gamma[\rho]-1\right)e^{-\beta\mathcal{H}[\rho]}$

The energy is

$\displaystyle\mathcal{H}[\rho]=-T_{0}^{2}\int_{\Omega}d^{2}\mathbf{x}\,d^{2}\mathbf{x}^{\prime}\,\rho(\mathbf{x})G(\mathbf{x}-\mathbf{x}^{\prime})\rho(\mathbf{x}^{\prime})=T_{0}^{2}\int_{\Omega}d^{2}\mathbf{x}\,\tilde{\psi}\rho$

where the second equality uses $\tilde{\psi}=\psi T_{0}$ from (6). By a counting argument [32, 60], the entropy $S[\rho]$ is related to $W[\rho]$ by

$\displaystyle S[\rho]=-\int_{\Omega}d^{2}\mathbf{x}\,\rho\log\rho=\log W[\rho]$

The partition function $Z(\beta)$, and the probability $P[\rho]$ of the density $\rho$ are therefore

	$\displaystyle Z(\beta)$	$\displaystyle=\int\mathcal{D}\rho\,\delta\left(\Gamma[\rho]-1\right)e^{S[\rho]-\beta\mathcal{H}[\rho]}$
	$\displaystyle P[\rho]$	$\displaystyle\propto\delta\left(\Gamma[\rho]-1\right)e^{S[\rho]-\beta\mathcal{H}[\rho]}$

The most probable distribution follows from maximizing $S-\beta\mathcal{H}$ subject to $\Gamma[\rho]=1$. In other words, we must maximize the functional

$\displaystyle\mathcal{I}[\rho]=S-\beta\mathcal{H}-\alpha\int_{\Omega}d^{2}\mathbf{x}\,\rho$

Upon rescaling $\beta$, this is exactly the maximization in (6) without angular momentum ($\gamma=0$). As before, this maximization yields the mean-field equations (7) with $\gamma=0$.