Emergent Non-Eulerian Hydrodynamics of Quantum Vortices in Two Dimensions

TL;DR

This paper introduces a coarse-grained model of quantum vortices in two dimensions, revealing a non-Eulerian fluid behavior with anomalous stresses and compressibility, supported by analytical and numerical results.

Contribution

It develops a novel large-scale description of quantum vortices as a non-Eulerian fluid with unique stress properties, extending classical vortex models.

Findings

Emergent vortex fluid exhibits non-Eulerian, compressible behavior.

Analytic vortex shear flow solution matches numerical simulations.

Identifies anomalous stresses caused by quantum vortex singularities.

Abstract

We develop a coarse-grained description of the point-vortex model, finding that a large number of planar vortices and antivortices behave as an inviscid non-Eulerian fluid at large scales. The emergent binary vortex fluid is subject to anomalous stresses absent from Euler's equation, caused by the singular nature of quantum vortices. The binary vortex fluid is compressible, and has an asymmetric Cauchy stress tensor allowing orbital angular momentum exchange with the vorticity and vortex density. An analytic solution for vortex shear flow driven by anomalous stresses is in excellent agreement with numerical simulations of the point-vortex model.

Tables

Table 1. Table 1: Numerical fitting parameters for N = 9522 𝑁 9522 N=9522 and different L 𝐿 L . For each L 𝐿 L , the results are for a single run (time evolution for a single sample from the initial state distribution). For L = { 1 , 2 , 3 , 4 } × 10 3 𝐿 1 2 3 4 superscript 10 3 L=\{1,2,3,4\}\times 10^{3} the fit is made at t = 2 𝑡 2 t=2 , and for L = { 5 , 6 , 7 , 8 , 9 , 10 } × 10 3 𝐿 5 6 7 8 9 10 superscript 10 3 L=\{5,6,7,8,9,10\}\times 10^{3} the fit is made at t = 50 𝑡 50 t=50 .

$L\left[{10}^3\right]$	$v_0\left(\mathrm{fitting}\right)$	$v_0\left(\mathrm{theory}\right)$	${\theta }_0$	$C$
$1$	$1.4265$	$1.5155$	$0.0152$	$0.0820$
$2$	$0.7811$	$0.7577$	$0.0828$	$0.0043$
$3$	$0.5181$	$0.5052$	$-0.0940$	$-0.0563$
$4$	$0.3778$	$0.3789$	$-0.0296$	$0.0175$
$5$	$0.3163$	$0.3031$	$-0.0094$	$0.0047$
$6$	$0.2332$	$0.2526$	$-0.0572$	$-0.0001$
$7$	$0.2143$	$0.2165$	$-0.0807$	$-0.0050$
$8$	$0.1903$	$0.1894$	$-0.1280$	$0.0180$
$9$	$0.1660$	$0.1684$	$0.2033$	$-0.0106$
$10$	$0.1567$	$0.1515$	$0.1662$	$-0.0029$

Equations

$${\mathcal{D}}_t^u\mathbf{u}=-\frac{1}{nm}\nabla p,\nabla \cdot \mathbf{u}=0,$$

$${\mathcal{D}}_t^u\omega =0.$$

$$H=\frac{nm}{2}\int d^2\mathbf{r}|\mathbf{u}{|}^2=\frac{nm}{2}\int d^2\mathbf{r}\psi (\mathbf{r})\omega (\mathbf{r}),$$

$$\frac{d{z}_i}{dt}=\frac{\partial \mathcal{H}}{\partial p_i},\frac{d{p}_i}{dt}=-\frac{\partial \mathcal{H}}{\partial z_i},$$

$$\mathcal{H}=-\pi {\gamma }^2\sum _{i\ne j}{\sigma }_i{\sigma }_j\log \left|\frac{z_i-z_j}{\ell }\right|.$$

$$\frac{d{\bar{z}}_i}{dt}$$

$${\partial }_t\rho$$

$${\partial }_t\sigma$$

$$J_{\mathrm{c}}$$

$$2\sum _{i\ne j}\frac{{\sigma }_i}{z-z_i}\frac{{\sigma }_j}{z_i-z_j}={\left(\sum _i\frac{{\sigma }_i}{z-z_i}\right)}^2+{\partial }_z\sum _i\frac{1}{z-z_i}$$

$$\rho w=\sigma u-2\eta i{\partial }_z\rho ,$$

$$\rho v=\rho u-2i\eta {\partial }_z\sigma .$$

$$\rho w=\sigma v-i\eta {\rho }^{-1}\left({\partial }_z{\rho }^2-{\partial }_z{\sigma }^2\right),$$

$${\omega }_v-\omega =\nabla \times (\mathbf{v}-\mathbf{u})=\eta \nabla \cdot \left({\rho }^{-1}\nabla \sigma \right).$$

$$\nabla \cdot \mathbf{v}=2\eta i\left[{\partial }_z({\rho }^{-1}{\partial }_{\bar{z}}\sigma )-h.c.\right]=-\eta \nabla \times ({\rho }^{-1}\nabla \sigma )\ne 0,$$

$${\mathcal{D}}_t^v\sigma =0,$$

$${\partial }_t(\rho v)+{\partial }_z{\mathcal{T}}_{z\bar{z}}+{\partial }_{\bar{z}}\mathcal{T}+\rho {\partial }_z(2p)=0,$$

$${\partial }_t(\rho v_a)+{\partial }_b{\mathcal{T}}_{ab}+\rho {\partial }_{a}p=0,$$

$${\Pi }_{ab}=-\sigma {\tau }_{ab}-8{\eta }^2\pi {\sigma }^2{\delta }_{ab}-{\eta }^2\sigma {\partial }_b({\rho }^{-1}{\partial }_a\sigma ),$$

$${\tau }_{xy}$$

$${\tau }_{xx}$$

$${\Pi }_{xy}-{\Pi }_{yx}$$

$$p_m$$

$$\frac{\partial {\mathcal{L}}_{ab}^{\mathrm{v}}}{\partial t}+{\partial }_c{\mathcal{M}}_{abc}={\Pi }_{ba}-{\Pi }_{ab}$$

$$\mathcal{H}[\rho ,\sigma ]$$

$$\sigma (x,y)={\sigma }_0\sin (\alpha y),\mathbf{v}=(v_0\cos (\alpha y),0),$$

$$\rho w$$

$$=\frac{1}{\pi }{\partial }_{\bar{z}}\sum _{i=1}^N\frac{{\sigma }_i}{z-z_i}\sum _{j,j\ne i}\frac{-i\gamma {\sigma }_j}{z_i-z_j}$$

$$=\frac{-i\gamma }{2\pi }{\partial }_{\bar{z}}\left[{\left(\sum _{i=1}^N\frac{{\sigma }_i}{z-z_i}\right)}^2+{\partial }_z\sum _{i=1}^N\frac{1}{z-z_i}\right]$$

$$=\sigma u-2\eta i{\partial }_z\rho ,$$

$${\rho }_{+}\equiv \sum _{i=1}^{N_{+}}\delta (\mathbf{r}-{\mathbf{r}}_i^{+}),\text{and}{\rho }_{-}\equiv \sum _{i=1}^{N_{-}}\delta (\mathbf{r}-{\mathbf{r}}_i^{-}),$$

$$u=-\sum _{i=1}^N\frac{i\gamma {\sigma }_i}{z-z_i(t)}=u^{+}+u^{-},$$

$$u^{+}=-\sum _{i=1}^{N_{+}}\frac{i\gamma {\sigma }_i^{+}}{z-z_i^{+}(t)}\text{and}{u}^{-}=-\sum _{i=1}^{N_{-}}\frac{i\gamma {\sigma }_i^{-}}{z-z_i^{-}(t)}$$

Full text

Introduction.— Topological defects play a key role in the two-dimensional (2D) superfluid transition Kosterlitz and Thouless (1972, 1973), and in many out of equilibrium phenomena including phase ordering dynamics after a quench Bray (1994); Mondello and Goldenfeld (1990); Bray et al. (2000) and turbulence in Bose-Einstein condensates (BECs) Barenghi et al. (2001); Vinen (2006); Barenghi et al. (2014); Tsatsos et al. (2016); Navon et al. (2016); Anderson (2016); Tsubota et al. (2017); Henn et al. (2009). Under planar confinement foo (a) the collective motion of many interacting quantum vortices in a superfluid leads to the emergent complexity of 2D quantum turbulence (2DQT) Numasato et al. (2010); Neely et al. (2013); Reeves et al. (2012); Bradley and Anderson (2012); Chesler et al. (2013); Neely et al. (2010); Kwon et al. (2014); Stagg et al. (2015); Reeves et al. (2013); Billam et al. (2015a); Skaugen and Angheluta (2016). Experimental advances now allow plane-confined BECs Chomaz et al. (2015); Gauthier et al. (2016), and simultaneous detection of quantum vortex positions and circulations Seo et al. (2017), opening the way to studies of 2DQT and Onsager vortices Onsager (1949); Eyink and Sreenivasan (2006); Joyce and Montgomery (1973); Montgomery and Joyce (1974); Billam et al. (2014); Simula et al. (2014); Yu et al. (2016); Groszek et al. (2016).

In the hydrodynamic regime at low temperature, superfluid vortex motion is much slower than the rate of sound propagation, and a system of many quantum vortices evolves as an almost isolated subsystem: a vortex fluid. The point-vortex model (PVM), a central model for studying 2D classical incompressible turbulent flows Onsager (1949); Eyink and Sreenivasan (2006); Novikov (1975); Aref et al. (1999), describes the dynamics of quantum vortices Fetter (1967); Haldane and Wu (1985), provided they are well-separated. In this regime the vortex core structure is unimportant, and the coupling to acoustic modes is relatively weak Lucas and Surówka (2014). Many studies of collective dynamics of vortices rely on large-scale numerical simulations of the discrete PVM Numasato and Tsubota (2009); Siggia and Aref (1981); Skaugen and Angheluta (2017); Reeves et al. (2017). An alternate approach was initiated by Wiegmann and Abanov Wiegmann and Abanov (2014), who formulated a hydrodynamic description of well-separated vortices of the same circulation, providing a rigorous starting point for studying rotating fluids, vortex clusters with definite sign of vorticity Smith (1989); Smith and O’Neil (1990), and the connection between vortex fluids and quantum Hall liquids Wiegmann (2013). A general 2D turbulent flow or a phase ordering process involves many vortices and anti-vortices, motivating a hydrodynamic theory of the binary vortex system.

In this Letter we develop a hydrodynamic formulation of systems involving a large number of vortices and anti-vortices, providing a framework for describing their emergent collective dynamics at large scales. A system containing many vortices with separation much bigger than the vortex core size $\xi$ can be treated as a fluid on a scale much larger than $\xi$, and smaller than the system size. We generalize the coarse-graining procedure proposed for chiral vortex systems Wiegmann and Abanov (2014) to the binary vortex system by introducing two hydrodynamical velocity fields via vortex number and charge currents. The binary fluid obeys a compressible hydrodynamic equation containing an asymmetric Cauchy stress tensor and an anomalous stress that is analogous to viscous stress while conserving energy. A vortex fluid shear flow driven by the anomalous stress is found, and numerical simulations of the PVM show excellent agreement with the analytical solution. Variation of the coarse-graining scale also demonstrates convergence of many-body dynamics of the PVM to non-Eulerian hydrodynamics. Dissipation due to the interaction between vortices and thermal cloud generates uphill diffusion of the vorticity at macroscopic scales. Chiral flow dynamics Wiegmann and Abanov (2014) is recovered as a special case of the binary fluid.

Two-Dimensional Hydrodynamics.— The dynamics of non-viscous incompressible classical fluids in 2D is described by the Euler equation

[TABLE]

where $n$ is the constant density, $\mathbf{u}$ is the fluid velocity, ${\cal D}^{u}_{t}\equiv\partial_{t}+\mathbf{u}\cdot\nabla$ is the material derivative with respect to $\mathbf{u}$, $m$ is the atomic mass, and $p$ is the fluid pressure determined by $\nabla\cdot(\mathbf{u}\cdot\nabla\mathbf{u})=-(nm)^{-1}\nabla^{2}p$. Taking the curl of Eq. (1), the Helmholtz equation for the vorticity $\omega\equiv[\nabla\times\mathbf{u}]_{z}$ is

[TABLE]

The kinetic energy of the fluid reads foo (b)

[TABLE]

where $\psi$ is the stream function, and ${\mathbf{u}}=\nabla\times(\psi\mathbf{\hat{z}})$ and $-\nabla^{2}\psi=\omega$. Here $\mathbf{\hat{z}}$ is a unit normal vector to the fluid plane.

For a 2D flow, it is convenient to use complex coordinates $z=x+iy$, $\partial_{z}=(\partial_{x}-i\partial_{y})/2$ and the complex velocity $u=u_{x}-iu_{y}$. In terms of complex notation $\nabla\cdot\mathbf{u}=\partial_{z}\bar{u}+\partial_{\bar{z}}u$, $\mathbf{u}\cdot\nabla=\bar{u}\partial_{z}+u\partial_{\bar{z}}$, $\omega=[\nabla\times\mathbf{u}]_{z}=i(\partial_{\bar{z}}u-\partial_{z}\bar{u})=2i\partial_{\bar{z}}u$, and $u=2i\partial_{z}\psi$. We also use subscripts $a,b$ to denoted the Cartesian components of vectors, and use vector, complex, or component notation where convenient.

Point-Vortex System.— A superfluid containing vortices with separation larger than the core size $\xi$ is nearly incompressible, away from vortex cores foo (c). For a BEC described by a macroscopic wavefunction $\Psi$, the associated Gross-Pitaevskii equation (GPE) governing time evolution can be mapped to the form Eq. (1) in the incompressible regime (constant density $n=|\Psi|^{2}$) Pitaevskii and Stringari (2016). The single valuedness of the wave function $\Psi$ requires that the circulation of a vortex excitation must be quantized in units of circulation quantum $\kappa\equiv 2\pi\hbar/m$, and the vorticity has a singularity at the position of the vortex core $\mathbf{r}_{i}$: $\omega(\mathbf{r})=\kappa\sigma_{i}\delta(\mathbf{r}-\mathbf{r}_{i})$ with the sign $\sigma_{i}=\pm 1$. Hereafter we set $nm=1$ for convenience.

We consider a system containing $N_{+}$ singly-charged quantum vortices and $N_{-}$ anti-vortices. The total number of vortices is $N=N_{+}+N_{-}$ and the vortex with sign $\sigma_{i}$ is located at $\mathbf{r}_{i}$. The fluid velocity generated by these quantum vortices far from the fluid boundary is completely determined by the vorticity: $u=2i\partial_{z}\psi=(2\pi i)^{-1}\int d^{2}\mathbf{r}^{\prime}\;\omega(\mathbf{r}^{\prime})/(z-z^{\prime})=-\sum^{N}_{j=1}i\gamma\sigma_{j}/(z-z_{j})$, where the stream function of the fluid is $\psi(\mathbf{r})=-\gamma\sum_{i}\sigma_{i}\log\left|(\mathbf{r}-\mathbf{r}_{i})/\ell\right|$, and the vorticity is $\omega(\mathbf{r})=2\pi\gamma\sum_{i}\sigma_{i}\delta(\mathbf{r}-\mathbf{r}_{i})$. Here $\ell$ is a length-scale introduced to ensure the correct dimension, and $\gamma=\kappa/2\pi$ is a convenient unit of circulation. As shown by Helmholtz, the above fluid velocity $u$ is a singular solution of Eq. (2). A quantum vortex generates a flow in the bulk fluid, while the vortex core is a point-like particle and has its own dynamics driven by the flow generated by the other vortices. The dynamics of the vortex cores has a Hamiltonian structure:

[TABLE]

with canonical momentum $p_{i}=i\pi\gamma\sigma_{i}\bar{z}_{i}$, and Hamiltonian

[TABLE]

The formal solutions of Eq. (4) are the Kirchhoff equations Kirchoff (1877)

[TABLE]

Using the vortex degrees of freedom, the kinetic energy of the fluid $H={\cal H}+E_{\rm self}$, where ${\cal H}$ is energy of interaction between vortices, and $E_{\rm self}=N\pi\gamma^{2}\log(\ell/\xi)$ is the total self-energy. The PVM can be seen as a limiting case of the vortex method Chorin (1994); Marchioro and Pulvirenti (2012) and well-approximates incompressible classical fluids with $\kappa$ determined by the injection scale.

Vortex Fluid Hydrodynamics.— For large $N$, the emergent collective dynamics of the discrete vortex system Eq. (6) can be described by a few hydrodynamic variables. By coarse-graining microscopic vortex distributions over patches containing many vortices, we derive a hydrodynamic formulation of the PVM, describing the vortex dynamics on scales greater than the patch scale. The core size $\xi$ is much smaller than the patch scale, and serves as a natural ultraviolet cut-off.

Under Hamiltonian evolution, conservation laws ensure the following continuity equations

[TABLE]

where vortex number density $\rho(\mathbf{r})\equiv\sum_{i}\delta(\mathbf{r}-\mathbf{r}_{i})$, vortex charge density $\sigma(\mathbf{r})\equiv\sum_{i}\sigma_{i}\delta(\mathbf{r}-\mathbf{r}_{i})=(2\pi\gamma)^{-1}\omega$, and

[TABLE]

are the currents for charge and number respectively. The coarse-grained vortex charge velocity field $w$ and vortex velocity field $v$ are defined according to the hydrodynamic relations $J_{\rm c}\equiv\rho w$ and $J_{\rm n}\equiv\rho v$ foo (d). Using the identity

[TABLE]

with Eq. (6), we can rewrite $J_{\rm c}$, to obtain an important relation between the vortex charge velocity $w$, and the fluid velocity $u$, given by

[TABLE]

with anomalous kinetic coefficient $\eta=\gamma/4$; here we have used $\partial_{\bar{z}}(1/z)=\pi\delta(\mathbf{r})$ and $[\partial_{z},\partial_{\bar{z}}](1/z)=0$. A fundamental relation linking the vortex velocity field $v$ to the fluid velocity $u$ can also be derived by decomposing $J_{\rm n}$ SM (1), to give

[TABLE]

The two velocity fields are related by

[TABLE]

and the vorticity of the vortex velocity field $v$, $\omega_{v}\equiv i(\partial_{\bar{z}}v-\partial_{z}\bar{v})$, has the anomalous correction

[TABLE]

In the hydrodynamic formulation the quantities $\rho$, $\sigma$, $J_{\rm c}$ and $J_{\rm n}$ represent averages of the corresponding microscopic quantities over patches, giving smooth coarse-grained quantities on the patch scale SM (2).

The relation Eq. (12) links the superfluid velocity field that is irregular at a vortex core, to a vortex fluid velocity field that is regular. In other words, the velocity of a vortex at position $\mathbf{r}$ is the fluid velocity excluding the flow generated by the vortex itself at $\mathbf{r}$. The regularization involves subtracting the singular term, namely the pole at the vortex core. For a single vortex at the origin the superfluid velocity $u=-i\gamma\sigma_{i}/z$ and $2i\eta\partial_{z}\sigma/\rho=4i\eta\sigma_{i}/z$, Eq. (12) yields $v=u-2i\eta\partial_{z}\sigma/\rho=0$. The correction cancels out the superfluid velocity field due to the local vortex, giving the physical result that a vortex does not move under the action of its own velocity field. Eq. (11) has a similar interpretation.

The binary vortex fluid is compressible:

[TABLE]

as also seen from $\nabla\cdot\mathbf{w}=\bar{u}\partial_{z}\left(\sigma/\rho\right)+u\partial_{\bar{z}}\left(\sigma/\rho\right)\neq 0$. In the chiral limit, $\sigma=\rho$, $w=v$ and $\nabla\cdot\mathbf{v}=0$ Wiegmann and Abanov (2014). The chiral vortex fluid is rigid, as the energy cost to compress a vortex fluid containing $N$ vortices scales as $N^{2}$ due to repulsive interactions between like-sign vortices. For a binary vortex fluid, the presence of two opposite sign vortices softens the vortex fluid such that gapless excitations can occur, making the vortex fluid compressible.

Anomalous Euler Equation.— Using Eqs. (8), (13), (15), $\sigma$ satisfies the vortex fluid Helmholtz equation

[TABLE]

and vortex charge is conserved, moving with the vortex velocity; use of Eq. (12) shows consistency with Helmholtz’ equation for the fluid, Eq. (2). A straightforward calculation SM (3) yields the anomalous Euler equation of the vortex velocity field $v$

[TABLE]

where ${\cal T}_{z\bar{z}}=\rho v\bar{v}+16\eta^{2}\pi\sigma^{2}+4\eta^{2}\sigma\partial_{\bar{z}}(\rho^{-1}\partial_{z}\sigma)$ and ${\cal T}=\rho vv+4\eta^{2}\sigma\partial_{z}(\rho^{-1}\partial_{z}\sigma)-4\eta i\sigma\partial_{z}v$ are complex components of the momentum flux tensor. Eq. (17) does not explicitly contain $w$ and Eqs. (7) (16) and  (17) fully describe the binary vortex fluid. In Cartesian coordinates Eq. (17) becomes

[TABLE]

where the momentum flux tensor ${\cal T}_{ab}=\rho v_{a}v_{b}-\Pi_{ab}$ Landau and Lifshitz (2013), with the emergent Cauchy stress tensor

[TABLE]

reflecting the macroscopic effects of the topological nature of quantum vortices. The anomalous stress

[TABLE]

does not cause energy dissipation and is formally identical to that of the chiral vortex fluid Wiegmann and Abanov (2014). The anomalous stress $\tau_{ab}$ vanishes in uniform rotation with angular velocity $\Omega$, where $v=-\Omega i\bar{z}$. Although there is no frictional viscosity in the binary vortex fluid, the shear stress ($\Pi_{ab}\neq 0$ for $a\neq b$) is non-vanishing, induced by gradients of $\rho$ and $\sigma$ and the anomalous stress $\tau_{ab}$. A similar situation can be found in the GPE when quantum pressure is important Winiecki et al. (2000).

A conspicuous feature of the Cauchy stress tensor $\Pi_{ab}$ for the binary vortex fluid is that it is asymmetric, with non-trivial commutator linking to the compressibility

[TABLE]

Note that under the transformation $x\leftrightarrow y$, the velocity $v_{x}\rightarrow-v_{y}$ and $v_{y}\rightarrow-v_{x}$, and hence $\nabla\cdot\mathbf{v}\rightarrow-\nabla\cdot\mathbf{v}$. For the chiral fluid $\sigma=\rho$, $\nabla\cdot\mathbf{v}=0$ and $\Pi_{ab}$ is symmetric Wiegmann and Abanov (2014). The local mechanical pressure of the binary vortex fluid is given by the normal stress:

[TABLE]

Angular Momentum.— The canonical angular momentum of the point-vortex system which is associated with rotational symmetry reads $L^{\rm c}\equiv\sum^{N}_{i}\mathbf{r}_{i}\times\mathbf{p}_{i}=-\pi\gamma\sum_{i}\sigma_{i}r^{2}_{i}$, and is equivalent to the fluid angular momentum $L_{\rm f}=\int d^{2}\mathbf{r}\;\mathbf{r}\times\mathbf{u}=-1/2\int d^{2}\mathbf{r}\;r^{2}\omega$. $L^{\rm c}$ is conserved as long as the system has rotational symmetry. We can also consider the orbital angular momentum (OAM) of vortices $L^{\rm v}\equiv\sum^{N}_{i}\mathbf{r}_{i}\times\mathbf{v}_{i}$. For a chiral vortex system ($\sigma_{i}=1$), $L^{\rm v}=2\eta N(N-1)$. In terms of the hydrodynamic field $v$, $L^{\rm v}=\int d^{2}\mathbf{r}{\cal L}^{\rm v}$ with ${\cal L}^{\rm v}=\mathbf{r}\times(\rho\mathbf{v})$. In Cartesian coordinates ${\cal L}^{\rm v}_{ab}=\rho(x_{a}v_{b}-x_{b}v_{a})$. Using Eq. (18), we obtain the continuity equation

[TABLE]

with angular momentum flux tensor ${\cal M}_{abc}=x_{a}{\cal T}_{bc}-x_{b}{\cal T}_{ac}$, indicating that the OAM is conserved if and only if $\Pi_{ab}=\Pi_{ba}$. Since $\Pi_{ab}$ is asymmetric for the binary vortex fluid, the OAM is not conserved regardless of the symmetry, consistent with the property of the corresponding discrete point-vortex system  SM (4). For the incompressible chiral vortex fluid ($\nabla\cdot\mathbf{v}=0$), the OAM is conserved. The divergence $\nabla\cdot\mathbf{v}$ is a source term in Eq. (23), and $\eta$ can be seen as a rotational viscous coefficient mediating the exchange between vortex fluid OAM and the internal vorticity and density degrees of freedom. The conserved OAM of the binary point-vortex system can be constructed by considering the sign-weighted OAM $L^{\rm w}\equiv\sum^{N}_{i}\mathbf{r}_{i}\times(\sigma_{i}\mathbf{v}_{i})=2\eta[(\sum_{i}\sigma_{i})^{2}-N]$. In the chiral case $L^{\rm w}=L^{\rm v}$ and $L^{\rm w}=-2\eta N$ for the neutral system; in terms of $w$, $L^{\rm w}=\int d^{2}\mathbf{r}\ \mathbf{r}\times(\rho\mathbf{w})$.

Vortex Fluid Hamiltonian.— Since $\sigma_{i}v_{i}\propto dp_{i}/dt$ and $\sum_{i}\delta(\mathbf{r}-\mathbf{r}_{i})\partial{\cal H}/\partial z_{i}=\sum_{i}\delta(\mathbf{r}-\mathbf{r}_{i})\sigma_{i}v_{i}$, the vortex charge velocity $w$ satisfies the canonical equation $-i\pi\gamma\rho w=\sigma\partial_{z}(\delta{\cal H}[\rho,\sigma]/\delta\sigma)+\rho\partial_{z}(\delta{\cal H}[\rho,\sigma]/\delta\rho)$, which gives the Hamiltonian of the binary vortex fluid foo (e)

[TABLE]

Here $H[\sigma]$ is the fluid Hamiltonian Eq. (3), and the second term is the self-energy of vortices foo (f). The Hamiltonian Eq. (24) is the hydrodynamic formulation to the discrete point-vortex Hamiltonian Eq. (5). It can be decomposed as ${\cal H}[\rho,\sigma]=H_{\rm kin}+H_{\rm int}+H_{\rm so}$, with vortex fluid kinetic energy $H_{\rm kin}=1/2\int d^{2}\mathbf{r}|\mathbf{v}|^{2}$, “internal energy” $H_{\rm int}=\eta^{2}/2\int d^{2}\mathbf{r}(\rho^{-2}|\nabla\sigma|^{2}-16\pi\rho\log\rho)$ and a “spin-orbit coupling” term $H_{\rm so}=-\eta\int d^{2}\mathbf{r}\rho^{-1}\mathbf{v}\cdot(\hat{z}\times\nabla\sigma)$ that couples vortex distribution energy to vortex fluid kinetic energy.

Vortex Shear Flows.— We consider a static vortex flow with $v_{y}=0$, $\partial_{x}v_{x}=0$ and $\partial_{x}\sigma=0$ in a neutral vortex system of constant density $\rho=\rho_{0}$. For such a static flow Eq. (18) reduces to $\partial_{b}{\cal T}_{ab}=\partial_{b}(\rho_{0}v_{a}v_{b}-\Pi_{ab})=\partial_{b}\Pi_{ab}=0$ foo (g), where the only non-trivial component reads $\left(\partial_{y}\sigma\partial_{y}+\sigma\partial^{2}_{y}\right)v_{x}+\eta\left[16\pi\sigma\partial_{y}+\rho^{-1}_{0}(\partial_{y}\sigma\partial^{2}_{y}+\sigma\partial^{3}_{y})\right]\sigma=0$, indicating that $\mathbf{v}$ is completely determined by the emergent Cauchy stress tensor including the anomalous stress $\tau_{ab}$. An exact shear flow solution to Eq. (18) can then be found:

[TABLE]

where $v_{0}=\alpha\eta\sigma_{0}(8\pi\alpha^{-2}-\rho^{-1}_{0})$, and $\alpha$ is a real parameter. Numerical simulations of the PVM with $N=9522$ foo (h) in a doubly-periodic square box with side length $L=6\times 10^{3}\xi$ SM (5) show excellent agreement with the prediction of Eq. (25) (see Fig. 1). Good agreement is also seen for $N=450$ and $L=300\xi$, being nearly within reach of current BEC experiments Kwon et al. (2014); Chomaz et al. (2015); Gauthier et al. (2016). There may also be more exotic vortex flows with enhanced anomalous term in Eq. (12), as may be caused by a gradient discontinuity in $\sigma$.

Dissipation and Annihilation.— To model dissipation due to interaction between superfluid vortices and a non-condensed thermal cloud Jackson et al. (2009), we consider a dissipative point-vortex model Rica and Tirapegui (1990); Törnkvist and Schröder (1997); Billam et al. (2015b); Moon et al. (2015) of the form $d\bar{z}_{i}/dt=v_{i}+i\eta^{\ast}\sigma_{i}v_{i}$ foo (i), where the dimensionless dissipation rate $\eta^{\ast}$ measures energy damping foo (j). In the presence of dissipation, the continuity equations (7) and (8) still hold subject to the substitutions $\mathbf{J}_{n}\rightarrow\mathbf{J}^{\prime}_{n}=\mathbf{J}_{n}-\eta^{\ast}\mathbf{\hat{z}}\times\mathbf{J}_{c}$ and $\mathbf{J}_{c}\rightarrow\mathbf{J}^{\prime}_{c}=\mathbf{J}_{c}-\eta^{\ast}\mathbf{\hat{z}}\times\mathbf{J}_{n}$. The dissipative Helmholtz equation is ${\cal D}^{\hat{v}}_{t}\sigma=-\eta^{\ast}\left(8\pi\eta\rho\sigma+\eta\nabla^{2}\sigma-\mathbf{v}\times\nabla\rho\right)$, where $\hat{\mathbf{v}}\equiv\mathbf{v}-\eta^{\ast}\eta\nabla\log\rho$, the terms $-8\pi\eta\eta^{\ast}\rho\sigma$ and $\eta^{\ast}\mathbf{v}\times\nabla\rho$ describe damping and transverse damping respectively. The negative sign in front of the diffusion term induces uphill diffusion of $\sigma$, concentrating vorticity, and may contribute to inverse energy cascades and vortex clustering processes  Billam et al. (2015a); Skaugen and Angheluta (2017).

Various forms of vortex number loss can modify the Hamiltonian theory, including vortex annihilations due to collisions, and dissipative vortex dipole decay or boundary loss. For a low temperature BEC containing many vortices the dissipation rate is small, $\eta^{*}\simeq 10^{-4}$ Bradley and Anderson (2012), and annihilations due to collisions are the dominant limitation for the Hamiltonian PVM approach. However, their influence is regime-dependent foo (k) and can be negligible for highly vorticity-polarized states, examples of which are the vortex shear flow (see  Fig. 1) foo (l), the enstrophy cascade in decaying 2DQT Reeves et al. (2017) and Onsager clustered states Reeves et al. (2014); Billam et al. (2014); Skaugen and Angheluta (2016). In the vortex dipole gas regime, dipole-dipole collisions are frequent and the theory is valid for times shorter than the dipole lifetime. The characteristic quantities are the collision cross-section, $\sigma_{\rm cs}\sim\xi$, the root-mean-square dipole velocity $v_{\rm{rms}}\sim\kappa d^{-1}$ (for $d\lesssim 10\xi$ the Jones-Roberts soliton Jones and Roberts (1982); Jones et al. (1986) forms), the mean free path ${\ell}_{m}\sim(\sigma_{\rm cs}\rho_{\rm d})^{-1}$, and the mean free time $\tau={\ell}_{m}/v_{\rm rms}\sim 10/(\kappa\rho_{\rm d})$, where $\rho_{\rm d}=\rho/2$ is the dipole density. The average total collision frequency per unit area is then $\tau^{-1}\rho_{\rm d}$, and assuming that every collision annihilates two vortices, we find the rate equation $dN/dt=-\Gamma_{2}N^{2}$, where the two-body decay rate $\Gamma_{2}=2\tau^{-1}\rho_{\rm d}A$ and $A$ is the system area. Such a decay behavior has recently been reported Kwon et al. (2014); using their experimental parameters we estimate $\Gamma_{2}\simeq 0.054\ s^{-1}$, close to the rate observed for the highest temperature. For the same parameters, the Hamiltonian PVM is a valid description up to the mean free time (lifetime) of a dipole, $\tau\sim 300{\rm ms}$ SM (6). For regimes at higher vortex energy than the dipole gas limit, the decay rate per vortex can be strongly suppressed Billam et al. (2014). Including vortex-antivortex annihilations into the hydrodynamic theory systematically is an important direction to explore in the future.

Conclusion.— Our hydrodynamical formulation of the point-vortex model reveals that the collective motion of many vortices emerges as a binary vortex fluid with rich phenomenology including an asymmetric Cauchy stress tensor and compressible flow. Examples of the former are relatively rare, and are associated with interactions between internal degrees of freedom and the bulk fluid flow, as occurs in liquid crystals de Gennes and Prost (1995). Excellent agreement between the analytic solution for vortex shear flow and the numerical simulations demonstrates the validity of the coarse-grained approach to collective vortex motion. The compressible vortex fluid may support density waves, however global plane-waves appear as excitations on top of non-trivial static background flows and are thus likely to be transient in general. The hydrodynamic theory suggests many areas for exploration including local density and vorticity excitations, regular vortex flows, connections with classical flows, and states of fully-developed two-dimensional quantum turbulence.

Acknowledgments

We acknowledge M. T. Reeves, B. P. Anderson, P. B. Blakie, Giulio Giusteri, L. A. Toikka, L. Williamson, and J. M. Floryan for useful discussions. We are grateful to P. Wiegmann and A. G. Abanov for helpful feedback on an earlier version of the manuscript. ASB is supported by a Rutherford Discovery Fellowship administered by the Royal Society of New Zealand.

I Supplemental Material

I.1 Hydrodynamic velocities

We first demonstrate Eq. (11) in the main manuscript. Combining Eqs. (6), (9), (10) in the main manuscript, we obtain

[TABLE]

where $\eta=\gamma/4$. In the last step we have used $\partial_{\bar{z}}(1/z)=\pi\delta(\mathbf{r})$ and $[\partial_{z},\partial_{\bar{z}}](1/z)=0$.

Let us now consider the relation between the vortex velocity field $v$ and the fluid velocity $u$. For this aim we introduce

[TABLE]

satisfying $\sigma=\rho_{+}-\rho_{-}$ and $\rho=\rho_{+}+\rho_{-}$. Here $\mathbf{r}^{+}_{i}$ is the position of a positive vortex and $\mathbf{r}^{-}_{i}$ is the position of an anti-vortex.

The fluid velocity has a decomposition:

[TABLE]

where

[TABLE]

are the fluid velocities generated by vortices with positive circulation and with negative circulation respectively. Here $\sigma^{\pm}_{i}=\pm 1$.

Introducing the velocity of a vortex with positive (negative) circulation

[TABLE]

we can write

[TABLE]

where

[TABLE]

We can further decompose ${\cal I}_{+}$ into

[TABLE]

where

[TABLE]

and

[TABLE]

In order to see the last step of Eq. (I.1) is true, let us integrate Eq. (I.1) over an arbitrary region $\Omega$ which contains $M$ positive vortices. On one hand, we obtain

[TABLE]

On the other hand side, we obtain

[TABLE]

showing that Eq. (I.1) is correct.

Similarly for ${\cal I}_{-}$ we have

[TABLE]

Collecting all the terms above we finally obtain Eq. (12) in the main manuscript:

[TABLE]

I.2 Orbital Angular Momentum of Point Vortex Systems

We show that the orbital angular momentum is not conserved for a binary vortex system with $N=3$. Firstly, let us consider $N=2$ with $\sigma_{1}=-\sigma_{2}=1$. The velocities of the two vortices read

[TABLE]

The orbital angular momentum of the two opposite vortices is

[TABLE]

and then

[TABLE]

We now consider $N=3$ with $\sigma_{1}=\sigma_{2}=-\sigma_{3}=1$, and then

[TABLE]

It is straightforward to find

[TABLE]

and in general

[TABLE]

For example for $z_{1}=1$, $z_{2}=2+i$ and $z_{3}=1+2i$, $\operatorname{Im}\left(\sum^{3}_{j=1}\frac{d\bar{v}_{j}}{dt}\bar{z}_{j}\right)=1$. Therefore, in general

[TABLE]

It is expected that Eq. (46) holds for a binary point vortex system with $N>3$. This is consistent with the corresponding coarse-grained vortex flow, where the Cauchy stress tensor is asymmetric.

I.3 Derivation of the Anomalous Euler equation

Here we outline the derivation of Eq. (17) in the main manuscript. Combining the Euler equation

[TABLE]

the continuity equation

[TABLE]

the Helmholtz equation

[TABLE]

and the relation

[TABLE]

we obtain

[TABLE]

This can be simplified as follows. Using Eqs. (48) and  (49), we have

[TABLE]

and

[TABLE]

Substituting Eqs. (52) and  (I.3) in to Eq. (51) and expressing $u$ in terms of $v$ and $\sigma$ using Eq. (50), many terms on the right hand of Eq. (51) are canceled, and we obtain

[TABLE]

Introducing the complex momentum flux tensors

[TABLE]

and

[TABLE]

Eq. (54) becomes

[TABLE]

which is Eq. (17) in the main manuscript.

I.4 Macroscopic Quantities

In this section we show how to construct macroscopic quantities from the corresponding microscopic quantities following a similar procedure in electromagnetism Jackson (2007). The microscopic definition of charge density reads

[TABLE]

Considering a patch with characteristic length scale $\ell$ containing $M$ point vortices, the average charge density in the patch reads

[TABLE]

Fluctuations on the scale smaller than $\ell$ are averaged out in Eq. (I.4), giving a relatively smooth macroscopic quantity on the scale of the patch size $\ell\gg\xi$. In this coarse-grained picture, a patch is viewed as a point $\mathbf{r}$ and the macroscopic fields describe collections of vortices at the patch. Similarly, for the vortex density

[TABLE]

The coarse-grained, macroscopic vortex number current reads

[TABLE]

The relation $\langle J_{\rm n}\rangle=\langle\rho\rangle v$ gives the vortex velocity field

[TABLE]

which is used to calculate the average vortex velocities for different bin sizes, showing that the average vortex velocities $M^{-1}\sum_{i}v_{i}$ approaches the continuous field $v$ with increasing the bin size (see Fig. 1 in the main manuscript).

I.5 Numerical Simulations of the Point vortex model and the Vortex Shear Flow

In the simulations we use the dissipative PVM

[TABLE]

We consider a neutral point vortex system in a square box with side length $L$. Since the vortex shear flow described by Eq. (25) in the main manuscript is periodic in the $y$ direction and homogeneous in the $x$ direction, we impose periodic boundary conditions. In a doubly-periodic square box, the velocity of the $i$th vortex in the presence of the $j$th, $\mathbf{v}_{i}^{(j)}$, reads Weiss and McWilliams (1991)

[TABLE]

where $\mathbf{r}_{ij}^{\prime}=(\mathbf{r}_{i}-\mathbf{r}_{j})(2\pi/L)$. To simulate the dynamics of large number of vortices in a doubly-periodic square box, we adopted the computational scheme developed in Ref. Reeves et al. (2017). We measure length $L$, velocity $v$ and time $t$ in units of $\xi$, $\kappa/\xi$ and $\xi^{2}/\kappa$ respectively.

The annihilation of vortex–antivortex pairs is modelled by removing vortex dipoles with separation less than the healing length $\xi$. To take into account the energy dissipation due to the sound wave radiation of accelerating vortices, we increase the dissipation $\eta^{*}$ for the same-sign vortex pairs when their separation approaches the healing length $\xi$ Pismen (1999). In our simulations we use the same dissipation as in Ref. Reeves et al. (2017):

[TABLE]

and we choose $r_{2}=\xi$ and $r_{1}=0.1\xi$ without losing generality. When computing the evolution of vortex $i$, the background dissipation rate $\eta^{*}=10^{-4}$ in Eq. (63) is replaced with the local dissipation rate $\eta^{*}_{i}$. During our simulations we find that vortex number loss is weak and the highest loss ratio is around $3\%$.

The initial vortex configuration is obtained by sampling vortex coordinates according to the distribution $\sigma=\sigma_{0}\sin(2\pi y/L)$, choosing $\sigma_{0}=\rho_{0}$. We sample positions of point vortices (antivortices) according to the distributions

[TABLE]

respectively, where $y\in[-L/2,L/2]$. We then make replicas of these vortices and distribute them uniformly along the $x$-axis. In our sampling the constant vortex density is anisotropic. Being anisotropic is irrelevant as long as the densities along the $x$-axis and the $y$-axis are defined at similar scales.

We keep track of the coordinates of vortices initially located at $x=0$. From the trajectories of these vortices we can extract the vortex velocity field distribution. There is no motion along the $y$ direction. We fit our results using the fitting function

[TABLE]

The fitting parameters $v_{0}$, $\theta_{0}$ and $C$ are shown in Table 1. The parameter $v_{0}$ barely changes if $\theta_{0}$ and $C$ are set to be zero for each $L$.

For each $L$ the fitting parameter $v_{0}$ remains nearly constant over a long time interval (see Fig. 2). Fig. 3 shows the coarse-grained vortex velocities for various bin sizes for a much smaller point vortex system and the results still show good agreement with the theory.

I.6 Vortex number Loss

The proposed theory describes emergent properties of well-separated quantum votices at large scales and is valid whenever the PVM is applicable and the vortex number loss is low. In various regimes of interest such conditions can be fulfilled and vortex losses are only a weak modification of the Hamiltonian vortex motion. In general, the collision rate depends on the vortex distribution and varies from case to case. For the vortex shear flow described in the main manuscript, our simulations of the dissipative PVM show that the vortex number loss is below $3\%$. In the recent study on the enstrophy cascade based on the dissipative PVM, the annihilation rate is also shown to be negligible (less than $1\%$) Reeves et al. (2017).

In a BEC vortex annihilation events involve physics outside the scope of the Hamiltonian PVM. The short range approach of vortices is associated with the formation of Jones-Roberts soliton Jones and Roberts (1982); Jones et al. (1986), a metastable density dip that occurs at the onset of the gradual process of complete annihilation into sound waves. The Jones-Roberts soliton forms at $d\approx 10\xi$ Jones et al. (1986); SM (7) marking the limit of validity of the PVM for BEC experiments. We estimate the annihilation rate for situations where vortex collision events are most frequent. In the absence of energy dissipation the vortex number loss is dominated by vortex dipole annihilations Stagg et al. (2015); Karl and Gasenzer (2016), being significant when the system is in the homogeneous dipole gas regime, where collisions happen between vortex dipole pairs. In this regime the vortex collision cross-section $\sigma_{\rm cs}\sim\xi$. The root-mean-square vortex velocity $v_{\rm{rms}}\sim\kappa d^{-1}\sim\kappa\ (10\xi)^{-1}$, the mean free path ${\ell}_{m}\sim(\sigma_{\rm cs}\rho_{\rm d})^{-1}$, and the mean free time $\tau={\ell}_{m}/v_{\rm rms}\sim 10/(\kappa\rho_{\rm d})$, where $\rho_{\rm d}=\rho/2$ is the dipole density. Hence the average collision frequency for one vortex dipole is $Z=\tau^{-1}=\kappa\rho_{\rm d}/10$ and the average total collision frequency per unit area is $Z_{t}=Z\rho_{\rm d}$. We assume that for every dipole-dipole collision there are two vortices being annihilated, then the vortex number loss during $\delta t$ is $\delta N(t)=N(t+\delta t)-N(t)=-2Z_{t}A\delta t$, where $A$ is the total area of the system.

Hence the decay rate equation reads

[TABLE]

leading to $\rho(t)\sim t^{-1}$. Such vortex decay behavior has been observed in the recent experiment Kwon et al. (2014) and is consistent with Refs. Schole et al. (2012); Stagg et al. (2015). Note that the pre-factor in Eq. (68) is over-estimated as some of the dipole-dipole scatterings may not evolve annihilation Schole et al. (2012). Taking the parameters in the experiment on ${}^{23}\rm Na$ Kwon et al. (2014), where the atomic mass $m_{\rm a}=23\times 1.7\times 10^{-27}\rm{kg}$ and the radius of Thomas-Fermi density profile $R_{\rm TF}=7\times 10^{-5}m$, we estimate the two-body decay rate $\Gamma_{2}=(1/20)\ \kappa A^{-1}\simeq 0.054\ s^{-1}$. The result is comparable to the hottest experiment and is around $10$ times of the value for the coldest experiment in Ref. Kwon et al. (2014). Here we used $A=\pi R^{2}_{\rm TF}$ and $\kappa=h/m_{\rm a}$ with $h\simeq 6.6\times 10^{-34}\rm{m^{2}kgs^{-1}}$. For the same parameters the mean free time (lifetime) of a dipole is $\tau\sim 300{\rm ms}$, setting a time limit up to when the Hamiltonian PVM is a valid description.

In above estimations, we assume that dipole-dipole scattering processes are dominant. Vortex decay rates will be modified if three-body scattering processes (a vortex dipole scattering with a free vortex) are important Karl and Gasenzer (2016), which could be the case when close to the Berezinskii-Kosterlitz-Thouless transition. It is worth mentioning that for a finite system, drift-out vortex number loss could be also relevant, while less important if energy dissipation is low  Stagg et al. (2015).

I.7 Experimental Accessibility

Experiments have realized large atom-number condensates of up to $N_{a}\sim 10^{8}$ atoms Streed et al. (2006); van der Stam et al. (2007), with atomic number densities $n_{0}\sim 10^{14}$ cm*-3*, using 23Na Streed et al. (2006). The healing length $\xi\approx 0.6\mu\mathrm{m}$. If such a condensate can be realized in a uniform quasi-2D trap with side length $L$ and thickness $\sim 6.6\xi$ Kwon et al. (2014), then this would allow that $L\sim 1000\xi$ and the vortex number $N\sim 450$. Since here $u_{\mathrm{rms}}/c\lesssim 0.3$, such a system is within the incompressible regime. Experiments in quasi-2D have achieved $N\sim 60$ vortices for $L\sim 120\xi$ in harmonically trapped systems Kwon et al. (2014), and hard-wall confinement with $L\gtrsim 200\xi$ Gauthier et al. (2016). With a lot of scope for larger condensates, and increasing experimental control, it is reasonable to expect conditions for which the hydrodynamic theory applies could be achieved in near-future BEC experiments.