Theory of the vortex-clustering transition in a confined two-dimensional quantum fluid
Xiaoquan Yu, Thomas P. Billam, Jun Nian, Matthew T. Reeves, Ashton S. Bradley
TTheory of the vortex-clustering transition in a confined two-dimensional quantum fluid
Xiaoquan Yu, ∗ Thomas P. Billam,
2, 3
Jun Nian,
4, 5
Matthew T. Reeves, and Ashton S. Bradley † Department of Physics, Centre for Quantum Science,and Dodd-Walls Centre for Photonic and Quantum Technologies, University of Otago, Dunedin, New Zealand. Joint Quantum Centre (JQC) Durham-Newcastle, School of Mathematics and Statistics,Newcastle University, Newcastle upon Tyne, NE1 7RU, United Kingdom. Joint Quantum Centre (JQC) Durham-Newcastle, Physics Department,Durham University, Durham, DH1 3LE, United Kingdom. Institut des Hautes ´Etudes Scientifiques, Bures-sur-Yvette, France. C. N. Yang Institute for Theoretical Physics, Stony Brook University, Stony Brook, U.S.A. (Dated: August 4, 2016)Clustering of like-sign vortices in a planar bounded domain is known to occur at negative temperature, aphenomenon that Onsager demonstrated to be a consequence of bounded phase space. In a confined superfluid,quantized vortices can support such an ordered phase, provided they evolve as an almost isolated subsystemcontaining su ffi cient energy. A detailed theoretical understanding of the statistical mechanics of such states thusrequires a microcanonical approach. Here we develop an analytical theory of the vortex clustering transition ina neutral system of quantum vortices confined to a two-dimensional disc geometry, within the microcanonicalensemble. The choice of ensemble is essential for identifying the correct thermodynamic limit of the system,enabling a rigorous description of clustering in the language of critical phenomena. As the system energy in-creases above a critical value, the system develops global order via the emergence of a macroscopic dipolestructure from the homogeneous phase of vortices, spontaneously breaking the Z symmetry associated withinvariance under vortex circulation exchange, and the rotational SO(2) symmetry due to the disc geometry.The dipole structure emerges characterized by the continuous growth of the macroscopic dipole moment whichserves as a global order parameter, resembling a continuous phase transition. The critical temperature of thetransition, and the critical exponent associated with the dipole moment, are obtained exactly within mean fieldtheory. The clustering transition is shown to be distinct from the final state reached at high energy, known as su-percondensation. The dipole moment develops via two macroscopic vortex clusters and the cluster locations arefound analytically, both near the clustering transition, and in the supercondensation limit. The microcanonicaltheory shows excellent agreement with Monte Carlo simulations, and signatures of the transition are apparenteven for a modest system of 100 vortices, accessible in current Bose-Einstein condensate experiments. I. INTRODUCTION
In two-dimensional (2D) classical fluids, giant coherentvortex structures can emerge from microscopic vortex exci-tations [1–3], as end states of turbulent fluid dynamics involv-ing an inverse energy cascade [4, 5]. The Great Red Spotin Jupiter’s atmosphere [6–8] is a well-known example. On-sager explained this spontaneous formation of large-scale vor-tices from an equilibrium statistical mechanics point of viewby studying a point-vortex model in a bounded domain [1, 9].The phenomenon stems from the bounded phase space, whichsupports negative temperature states that favour the sponta-neous clustering of like-sign vortices. As an equilibrium sta-tistical mechanics problem, cluster formation in 2D vortexsystems has attracted much attention [5, 10–20]. Buildingon Joyce and Montgomery’s formulation of vortex statisticalmechanics [10], Smith and O’Neil [17, 18] investigated thesingle charge 2D plasma using the point-vortex model in arotating container and showed that the formation of a non-axisymmetric cluster resembles a second-order phase transi-tion. The analogy to a phase transition is a powerful tool andprovides deep insights into the nature of vortex cluster forma- ∗ [email protected] † [email protected] tion. To the best of our knowledge, the nature of the transi-tion to large-scale clustered states in bounded neutral vortexsystems remains an open problem and fundamental questionsremain unanswered: What is the exact transition temperature,and the universality class of the transition? What is the precisemechanism of cluster formation? A rigorous theory capable ofaddressing these questions may expand our knowledge of crit-ical phenomena [21] and pattern formation [22], and lead toa deeper understanding of collective vortex states in a broadclass of 2D systems [23].In this work we develop an analytical treatment withinmean field theory [10, 18] of the clustering transition in a neu-tral point-vortex model confined to a uniform disc geometry.The theory provides a systematic description of the negativetemperature phenomena in terms of inverse temperature β ( ε ),a function of the system energy ε parametrizing microcanon-ical states. We analytically find β c < a r X i v : . [ c ond - m a t . qu a n t - g a s ] A ug at β s < β c [4], involving point-like concentration of vorticeswith the same circulation, and divergence of the system en-ergy. The supercondensation does not break any symmetry (itoccurs within the symmetry broken phase) and is an end pointof a continuous process as energy increases. We test the mi-crocanonical theory against Monte Carlo (MC) sampling ofthe microcanonical ensemble for finite vortex numbers, andfind good agreement even for a modest system containing atotal of N =
100 vortices. To probe the role of finite vor-tex number further, we compute the energy dependence of themacroscopic dipole moment, for a range of N from 50 to 1000.The critical exponent determining the order parameter scalingemerges clearly even for N =
100 vortices.
II. BACKGROUNDA. Negative Temperature Vortex States
As a measure of motion, absolute temperature is alwayspositive. However under certain conditions, negative temper-ature states can occur for an almost isolated subsystem. Theconcept of temperature is applicable to such a subsystem ifit is in a long-lived metastable state. Since the coordinates( x , y ) of a vortex are conjugate variables, the crucial featureof 2D point-vortex systems in a bounded domain is that thephase space volume is finite. Let the phase space volumebe Ω ( ε ) for a given energy ε , then the thermodynamic en-tropy is k B ln Ω ( ε ), where k B is the Boltzmann constant. Forthe extreme situation in which all the vortex dipoles (vortex-antivortex pairs) collapse, we then have Ω ( ε = −∞ ) =
0. Theopposite limit occurs when all the like-sign vortices concen-trate into a point, in which case Ω ( ε = ∞ ) =
0. It is clearthat Ω ( ε ) must reach its maximum at an intermediate energy −∞ < ε m < ∞ . Then the inverse temperature1 T ≡ k B ∂ ln Ω ( ε ) ∂ε (1)is thus negative for ε > ε m . High energy equilibrium vor-tex states at negative temperature correspond to large-scaleclusters [9]. Negative temperature states can arise in systemswith finite and discrete degrees of freedom, for instance lo-calized spin systems [24–26]. Recently such negative temper-ature states have been realized for motional degrees of free-dom [27], where the finite band width sets an upper bound forthe energy spectrum [28]. B. Motivation from Bose-Einstein condensates
An interesting scenario for potential realization of Onsagervortex states may occur in planar Bose-Einstein condensates(BECs). In atomic BECs the point-vortex model provides anexcellent description of vortex dynamics in a hydrodynamicregime [29–32]. Furthermore, experimental observations ofsmall-scale clustering [33, 34], and Gross-Pitaevskii simula-tions demonstrating clustering in forced [35, 36] and decaying[37] homogeneous systems, in harmonic traps [38], and disc geometries [39] have raised the prospect that large-scale clus-tered states may be realized, amid increasing theoretical [40–49] and experimental [32, 50–53] interest.The scenario of negative temperature vortex states withina low positive temperature superfluid may appear implausi-ble as the superfluid in a BEC is compressible and the vortex-sound coupling must be taken into account. In general, vortex-sound coupling can play an important role in a compressiblesuperfluid [23, 54, 55], however, there are also regimes wherethe acoustic loss rate can be much slower than the dynami-cal timescales of vortex evolution. Several theoretical workshave considered the viability of negative temperature vortexstates using the Gross-Pitaevskii equation (GPE) [37, 39, 40].They show that the compressibility of the quantum fluid doesnot prevent the clustering of like-sign vortices, and moreoverit can even provide a mechanism to drive clustering, via anevaporative heating process [39], whereby annihilations oflow-energy vortex-antivortex pairs “heat up” the remainingvortices and the end state shows a clustered configuration oflike-sign vortices. Experimental observations that like-signvortices can aggregate into long-lived multiply charged struc-tures provides further evidence that such a regime may be ac-cessible [33, 34]. These experiments can be seen as first stepstowards realizing large-scale, negative temperature, clusteredstates of quantum vortices.Several works have considered the possibility of an inverseenergy cascade in two-dimensional BEC, suggesting that thecompressibility of the fluid may be fatal for large-scale clus-tering, instead causing transport of energy to small scalesdriven by vortex dipole annihilation [56–58]. However, theseworks considered specific scenarios that immediately prohib-ited clustering, due to either starting from initial conditionsdominated by acoustic energy [56, 57], or evolving the sys-tem according to an over-damped equation of motion [58]. Asshown in a systematic study of energy transport [44], a clearregime exists where energy is transported to large scales via avortex clustering process. The regime is one of weak dissipa-tion, and large vortex number. Recent analytical calculationsfurther support the possibility that the sound-vortex interac-tion is not strong enough to prevent energy transport to largescales [42], consistent with GPE simulations [35–39, 46].In summary, although the temperature of the bulk liquid isalways positive, rendering the negative temperature states ul-timately unstable, if the temperature of a BEC is low enoughand the mean distance between vortices is much larger thanthe healing length, the clustered phase is expected to be long-lived. Indeed, GPE simulations indicate that a regime existswhere a separation of timescales allows clustered states to belong-lived. The cause of this timescale separation appearsto have two sources: (i) the vortex-antivortex recombinationrate is rather slow in 2D (stemming from weak coupling tothe sound field; this weak coupling also holds for single vor-tices), (ii) clustering of same sign vortices suppresses closeapproaches of vortices and antivortices. Thus, in a long timewindow, a system of 2D quantum vortices may evolve as analmost isolated subsystem of a confined BEC at su ffi cientlylow temperature. Such quasi-equilibrium states may also haveconnections with non-thermal fixed point behavior [59, 60].While much evidence points to a weak coupling regime, a de-tailed microscopic picture of the mechanism remains an openproblem and such states have yet to be observed in the labo-ratory, posing a exciting future challenge for both theoreticaland experimental work. C. Point-Vortex Regime
We consider a total of N quantum vortices in a planar su-perfluid, with core size determined by the healing length ξ .In a BEC at su ffi ciently low temperature, coupling to the ther-mal cloud is negligible, and provided that the average distancebetween vortices greatly exceeds the healing length, vorticesare only weakly coupled to phonon modes. The point-vortexmodel then applies, and the other degrees of freedom can beneglected. In this regime the point-vortex model can also berigorously obtained from the GPE [29, 37] and the incom-pressible kinetic energy spectrum of quantum vortices in theGPE is well approximated by the point-vortex energy spec-trum [35, 37].Let us now consider N = N + + N − quantum vortices con-sisting of N + vortices and N − = N + anti-vortices confinedto a uniform disc with radius R , by a hard-wall boundary.Hereafter we set R =
1, giving the point-vortex Hamilto-nian [61, 62] H = − (cid:88) i (cid:44) j κ i κ j ln | r i − r j | + (cid:88) i , j κ i κ j ln (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( r i − ¯ r j ) | r j | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (2)where κ i = ± r i . Theimage terms are given by the second summation in Eq. (2) and¯ r j = r j / | r j | is the position of the j th image vortex. The imageterms ensure that the velocity normal to the boundary van-ishes; without losing generality we choose images such thatthe stream function is zero on the boundary. The Hamiltonian(2) has Z symmetry (invariant under κ i → − κ i ) and rotationalSO(2) symmetry due to the disc geometry. In a dilute gasBEC the Hamiltonian Eq. (2) is measured in the energy unit E = ρ m a κ / π , where ρ is the superfluid density, κ = h / m a is the circulation quantum and m a is the atomic mass. Eq. (2)describes well-separated ( ξ √ N (cid:28)
1) quantum vortices in 2Dincompressible superfluids. Note that the second summationof Eq. (2) depends on the shape of the domain and is absent foran unbounded domain. The disc geometry is chosen not onlybecause it can be realized in current BEC experiments [63–65] but also the simple form of the domain leads to the pos-sibility of an exact analytical treatment. In general the point-vortex model also describes vortices in classical inviscid, in-compressible fluids [5, 15, 16, 66], and guiding-center plasmadynamics [10, 12].
D. Scaling and Continuous Hamiltonian
The model Eq. (2) exhibits negative temperature states athigh energies. In the standard thermodynamic limit, namely N → ∞ and V → ∞ with H / N and N / V held fixed, where V is the area, negative temperature states do not exist. Inorder to support negative temperature states, a bounded do-main is necessary. A well-defined thermodynamic limit forthe vortex system in a bounded domain is found by a carefulchoice of scaling [1, 67, 68]. In the clustered phase, the energy H ∼ N is due to the sum over all vortex pairs. Hence H / N is finite-valued as N → ∞ . Moreover the entropy per vortexis bounded by the area of the domain. A formulation of thetheory with a well-defined thermodynamic limit that can de-scribe negative temperature states is thus given by the rescaledHamiltonian H ≡ H / N , (3)equivalent to introducing the rescaling κ i → κ i / N ± = κ i / N .The validity of this rescaling for the bounded vortex systemhas been verified rigorously [67, 68].A course-grained formulation of Eq. (3) valid for N (cid:29) n ± ( r ) ≡ N ± (cid:88) i δ ( r − r ± i ) , (4)satisfying the normalization condition (cid:82) d r n ± =
1, where r ± i is the position of the vortex i of circulation ± κ i respectively,then the vorticity field is σ ( r ) ≡ n + ( r ) − n − ( r ) . (5)At large N ± , the continuous coarse-grained Hamiltonian reads H e ff = (cid:90) d r d r (cid:48) σ ( r ) φ ( r , r (cid:48) ) σ ( r (cid:48) ) , (6)which gives a well-defined continuum limit of Eq. (2) as N → ∞ and in the clustered phase H e ff ∼ O (1) [1, 67, 68].The potential φ ( r , r (cid:48) ) is the Green’s function of the Laplaceoperator ∇ φ ( r , r (cid:48) ) = − πδ ( r − r (cid:48) ) , (7)where φ ( r , r (cid:48) ) = φ ( r (cid:48) , r ) , and φ ( r , r (cid:48) ) ∼ − | r − r (cid:48) | as | r − r (cid:48) | → ψ ( r ) ≡ (cid:90) d r (cid:48) φ ( r , r (cid:48) ) σ ( r (cid:48) ) (8)satisfies the Poisson equation ∇ ψ = − πσ ( r ) , (9)with the boundary condition ψ ( r = , θ ) = C . Here θ is theazimuthal angle and C is a constant chosen to be zero, whichis equivalent to including image terms in Eq. (2). Recall thatthe radial velocity v r = r − ∂ψ/∂θ , and this boundary condi-tion ensures that there is no flow across the boundary of thecontainer.In the disc container, the energy H e ff = E = (cid:90) d r σψ = π (cid:90) d r |∇ ψ | ≥ L = (cid:90) d r r σ (11)are conserved quantities. E. Statistical Mechanics
Let us briefly recall the statistical mechanics of the systemdeveloped by Joyce and Montgomery [10], who derived self-consistent equations for the vorticity field via a maximum en-tropy principle.The most probable density distribution of vortices is ob-tained by maximizing the entropy per vortex S = − (cid:90) d r ( n + ln n + + n − ln n − ) , (12)with N + , N − , E and L constrained to fixed values in the micro-canonical ensemble. The variational equation δ S − βδ E − αδ L − µ + δ N + / N + − µ − δ N − / N − = , (13)where β , α and µ ± are Lagrange multipliers, gives the maxi-mum entropy state n ± ( r ) = exp (cid:104) ∓ βψ ( r ) ∓ α r + γ ± (cid:105) , (14)with γ ± = − µ ± −
1. Here β = β ( E ) is the dimensionless inversetemperature describing the vortex subsystem at fixed energy, α is proportional to the rotation frequency of the container, and γ ± set the normalization. Eq. (14) together with Eq. (9) pro-vides a mean-field description of the distribution of vorticesin a disc [1]. III. CLUSTERED PHASES
In the following we investigate possible stable large-scalecoherent structures described by Eqs. (9), (14). Non-trivialsolutions for the vortex density exist only for β <
0. We notethat some exact solutions of Eqs. (9), (14) in rectangular do-mains are known [70–73], however, we are unaware of anyexact solutions for the disc geometry. Our approach is to ex-tend the bifurcation theory of Smith and O’Neil [17, 18], orig-inally developed by for the single charge vortex system, to theneutral vortex system confined to a disc.
A. Constrained Variation of the Maximum Entropy State
Let us start at a solution n ± of Eqs. (9), (14) at E and L , andconsider a nearby solution n ± + δ n ± at E + δ E and L + δ L . The corresponding changes of the constraints are0 = (cid:90) d r δ n ± , (15) δ E = (cid:90) d r ψδσ + (cid:90) d r δψδσ, (16) δ L = (cid:90) d r r δσ. (17)The variation of Eq. (9) is(4 π ) − ∇ δψ = − [ δ n + ( r ) − δ n − ( r )] = n ( ψδβ + βδψ + r δα ) − n + δγ + + n − δγ − , (18)where n = n + + n − , and δγ ± , δβ and δα can be expressed aslinear functions of δ E , δ L and δψ . In contrast to the singlecharge vortex system [17, 18], the macroscopic dipole forma-tion in the neutral vortex system does not require finite rota-tion frequency of the container, and energy (equivalently tem-perature) is the only control parameter of the transition forgiven angular momentum L .We now focus on the non-rotating container, setting α = n ± = n / = /π is the solutionof Eqs. (9), (14), at which σ = ψ = E = L =
0. Evaluating the variations of the constraints Eqs.(15), (16),(17) from the homogeneous state to the leading order in δψ gives δγ + + δγ − = , (19) δ E = , (20) δγ + = β n (cid:90) d r δψ. (21)Plugging Eqs. (19), (20), (21) into Eq. (18), we obtain the zeromode ( δ E =
0) equation (cid:104) ∇ − πβ n (1 − ˆ I ) (cid:105) δψ = , (22)where the integral operator is defined asˆ I δψ ≡ π (cid:90) d r δψ, (23)and the solutions are required to satisfy the boundary con-dition δψ ( r = , θ ) =
0. The onset of clustering occurs ifEq. (22) has non-zero solutions for the fluctuation δψ . Notethat in our mean field approach the value of β is undefined inthe homogeneous phase, and its precise value at the transitionis determined by the specific unstable mode that develops at E = B. Vortex Clustering Transition
We now give a formulation of the clustering transition asa symmetry-breaking perturbation of the homogeneous neu-tral vortex configuration, and investigate the system depen-dence on energy near the transition. It is natural to decomposeEq. (22) into azimuthal Fourier harmonics ψ l ( r , θ ) = ˆ ψ l ( r ) e il θ (24)with the integer mode number l . The normalization conditionfor ψ l is chosen to be − (cid:90) d r ψ l ∇ ψ l = π (25)for convenience. Let us now consider a small variation alongthe mode l , δψ = δψ l ≡ (cid:15)ψ l , where (cid:15) is a small amplitude.In this work we focus on l > Z symmetry. Break-ing the Z symmetry while preserving the SO(2) symmetryinvolves the l = L . This polar angle θ -independent clustered phase is not pursued further here. For l >
0, ˆ I δψ =
0, and consequently δγ + = δγ − = δ L = (cid:34) r ddr r ddr − l r (cid:35) ψ l = − λψ l , (26)with the Dirichlet boundary condition ψ l ( r = , θ ) = λ = j l , m , where j l , m is the m th positive zero of the Bessel function of the firstkind J l ( r ). Due to the dipolar form of the dominant symmetrybreaking perturbation, the root j , (cid:39) .
832 plays an importantrole in what follows. The onset of the θ -dependent clusteredstable phase occurs when the l = m = β c = − j , / π n (cid:39) − . , (27)and the corresponding normalized non-zero solution is ψ ( r , θ ) = √ (cid:2) | J ( j , ) | j , (cid:3) − J ( j , r ) cos( θ ) . (28)Here we choose the positive y -axis as the dipole directionwithout loss of generality. Higher values of l and m corre-spond to additional ordered phases with lower entropy at finiteenergy, and hence have less statistical weight. For zero angu-lar momentum L = E (cid:44)
0, the maximumentropy state arises from l = m = T c = − π n ( j , ) E k B N ± (cid:39) − . π (cid:126) ρ Nm a k B . (29)The clear role of the boundary condition distinguishes theclustering transition from supercondensation, and uniquelydetermines the numerical factor.As a consistency check, the critical temperature (29) canalso be found by varying the entropy (12) along the mode l = δ S = S − S = − (cid:32) n β c (cid:90) d r ψ (cid:33) δ E + o ( δ E ) , (30) where S = − (cid:82) d r ( n /
2) ln( n /
2) is the entropy of theuniform phase, and in the last step we used Eq. (A9).The energy change along the mode l = δ E = (cid:90) d r δσδψ (cid:39) − (cid:15) n β c (cid:90) d r ψ = (cid:15) , (31)where the normalization condition Eq. (25) is used.We then immediately have ∂ S ∂ E (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) E = = − n β c (cid:90) d r ψ = β c , (32)consistent with the foregoing analysis.In the dipole phase, the macroscopic dipole moment D ≡ | D | ≡ (cid:12)(cid:12)(cid:12)(cid:12)(cid:88) κ i r i (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) d r r σ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (33)serves as the order parameter, which becomes non-zero spon-taneously as the system crosses the transition point. Impor-tantly here D does not refer to local order or o ff -diagonal long-range order but rather describes the global coherent structureof the clustered phase. Near the transition, 2 δ E = (cid:15) , andhence D = D ( E − E c ) ν + O ( δ E ) = ˜ D | β − β c | ν + O ( δβ ) (34)with the critical exponent ν = / , (35)and E c =
0. The critical exponent ν takes the same value as thecritical exponent for the condensation field in the mean fieldXY model. The coe ffi cients are D = (cid:112) − J ( j , ) / J ( j , ) = D = C c /β c , where C c ≡ − β c ∂ E ∂β (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) β = β c (cid:39) .
91 (36)is the microcanonical specific heat at the vortex clusteringtransition point (see Appendix A). The positiveness of the mi-crocanonical specific heat in the neutral dipole phase suggeststhat ensemble equivalence may hold for the confined neutralvortex system. This behavior is in stark contrast with thesingle charge dipole phase where the microcanonical specificheat is negative [18]. Note that D and C c do not depend onthe chosen normalization of ψ .Clear evidence of the emergence of global order in the vor-tex clustering transition can be seen in real space. In the dipolephase, | σ | is symmetric under x → − x , and ψ (0 , y ) =
0. Di-viding the disc into two regions A ± = { x ≷ } , we define theaverage position of net vortices in region A ± as r ± ≡ (cid:82) A ± d r r | σ | (cid:82) A ± d r | σ | . (37)The separation between vorticity centres can be measuredby d = | d | ≡ | r + − r − | . (38)Evaluating Eq. (37) close to the clustering transition gives d = d c + O ( δ E ) (39)with d c = H − ( j , ) (cid:39) . , (40)and H ( x ) is the first-order Struve function [74]. C. Supercondensation
Far above the vortex clustering transition each clustershrinks into a rather small region centred at r ± = | r ± | ≡| r ∓ d / | . In the supercondensation regime the angular varia-tion of the vortex distributions in the neighbourhood of r ± = n ± = A (cid:48) (2 − πβ A (cid:48) r ± ) , (41)with the smoothness condition n (cid:48)± ( r ± = =
0, where A (cid:48) isfixed by the normalization condition. For high energy, A (cid:48) → A ( β ) with A ( β ) ≡ π (1 − β/β s ) , (42)and β → β s = −
2. The inverse temperature β s defines theuniversal supercondensation point, at which vortices concen-trate spatially into point-like clusters, irrespective of eitherthe boundary condition or the net vortex charge [18]. InBEC units, supercondensation occurs at T s = − π (cid:126) ρ N / m a k B ,consistent with an estimate based upon the free energy argu-ment [5] and with canonical MC sampling [39]. Viewing theclusters as points, the equilibrium condition for the two pointvortices with opposite circulation yields d → d s , where d s = √ − / (cid:39) .
972 (43)is independent of the cluster charges. This equilibrium con-figuration balances the force of interaction between the twoclusters with the force arising interaction with the container;the latter is generated by the images in Eq. (2). The limitingforms Eq. (40) and Eq. (43) give the cluster separation for en-ergies near and well above the transition, and clearly the loca-tion of the clusters is only rather weakly dependent on systemenergy.The main contribution to the energy comes from vorticeswith the same circulation, hence E (cid:39) − β − (cid:2) ln (1 − β/β s ) − β/ (cid:3) . (44) As β → β s , |(cid:104) D (cid:105) − (cid:104) D (cid:105) max | ∼ β − β s with (cid:104) D (cid:105) max = d s , and themicrocanonical specific heat has asymptotic form C s = − β ( ∂ E /∂β ) ∝ ( β − β s ) − > . (45)The power law divergence of the microcanonical specific heatis determined by the asymptotically exact vortex density dis-tribution Eq. (41). In any physical realization of the point-vortex model short-range repulsion becomes important in thesupercondensation limit. Hence, instead of collapsing to apoint associated with infinite energy, a tight crystal latticestructure would occur [39]. IV. MONTE CARLO SIMULATIONS
Starting from the uniform phase at minimum energy E = L = A. Microcanonical Monte Carlo
Our microcanonical Monte Carlo simulations follow thegeneral scheme of Refs. [17, 18], based on the Demon algo-rithm of Creutz [75], modified to e ffi ciently sample the phasespace of the neutral point-vortex system. Specifically, we cou-ple the system of point vortices to a demon with energy andangular momentum degrees of freedom E D and L D , and writethe total energy and angular momentum as E Total = E + E D , (46) L Total = L + L D , (47)where E and L are the point-vortex energy and angular mo-mentum. To sample the microcanonical ensemble for desiredenergy and angular momentum E (cid:48) and L (cid:48) , we fix E Total = E (cid:48) and L Total = L (cid:48) and perform a random walk subject to theconstraints | E D | < E max and | L D | < L max . Typical con-straints on the Demon used here are E max = . × − and L max = × − .Our random walk consists of two alternating types of MonteCarlo step. The first Monte Carlo step is comprised of N pair-wise random moves of two randomly-chosen vortices, suchthat r i → r i + δ r i , where i = , δ r i are chosen uniformly randomly froma square of side length ∆ centered on the origin. We update ∆ to maintain an e ffi cient success rate for the random walk(typically such that 50% of pairwise moves succeed). ThisMonte Carlo step is very similar to that used in Ref. [18] for single-circulation point-vortex systems.To more e ffi ciently explore the possible values of the dipolemoment in the microcanonical ensemble of neutral point-vortex systems, we follow each Monte Carlo step of the pari-wise type described above with a single Monte Carlo step of asecond type. This second Monte Carlo step consists of a simi-lar pairwise move, but in this case applied to a pair consistingof; (a) the center of circulation of all positive vortices, and (b)the center of circulation of all negative vortices. Such a move D N ( ∆ D ) | d | E d s d c -1.00.01.0 ˜ σ E − E c = . E − E c = . E − E c = . (a)(b)(c)(d) FIG. 1. (Color online) A comparison between microcanonical the-ory and MC sampling for N =
100 point-vortices on the clusteringtransition. (a) shows the dipole moment D (solid red) for increasingenergy E ; also shown is the mean field prediction D = ( E − E c ) / (dotted black). E c (cid:39) − × − — here representing the critical en-ergy for N =
100 — is estimated by fitting to the numerical data[see Sec. IV C]. The inset in (a) shows the variance of dipole mo-ment for N =
50 (dash dotted green), N =
100 (solid red), and N =
200 (dashed blue). (b) shows the separation between vorticitycentres d ( E ), Eq. (38); the numerical data (solid line) is shown along-side asymptotic mean field predictions for high energies d s [double-dashed line, Eq. (43)] and energies close to the transition d c [dot-dashed line, Eq. (40)]. (c) and (d) show, respectively, results fromMC-sampling and mean field theory for the scaled vortex density ˜ σ (see text). To present the large energy range on one color map, wedefine ˜ σ ≡ σ for E − E c = . , .
1, and ˜ σ ≡ σ/
10 for E − E c = . is achieved by moving all the positive (negative) vortices by avector δ r + ( δ r − ), such that1 N + (cid:88) κ i > r i → N + (cid:88) κ i > r i + δ r + , (48)1 N − (cid:88) κ i < r i → N − (cid:88) κ i < r i + δ r − . (49)The random moves δ r ± are chosen from the same distributionas the pairwise moves above, using the same value of ∆ . Theaddition of this second type of step leads to greater ergodicityin the values of dipole moment obtained during the randomwalk. N E c [10 − ] ∆ E c [10 − ]50 − − − − . . . E c ( N ) and its uncertainty ∆ E c ( N ) ob-tained by our numerical fitting procedure [see Sec. IV C]. B. Vortex Clustering Transition for N = A comparison between the analytical predictions of meanfield theory and the microcanonical MC sampling for N = ν and D are in good agreement with the mean fieldpredictions. Before the transition (cid:104) D (cid:105) has a small value whichis almost independent of energy, due to the finite value of N .Here the brackets refer to the ensemble average. The vari-ance of the order parameter ( ∆ D ) ≡ (cid:104) D (cid:105) − (cid:104) D (cid:105) is evaluatednumerically and it shows a peak around the transition point,where the dipole structure becomes unstable (see Fig. 1(a) in-set). Moreover, N ( ∆ D ) grows with increasing N near thetransition point while it has almost no N -dependence else-where. These non-Gaussian fluctuations of the dipole momentindicate non-trivial correlations between vortices due to crit-icality. As a check of our MC-sampling, we independentlyverified the results in Fig. 1 (a) by sampling using the Wang-Landau algorithm [76].The values of d c and d s are also confirmed by MC simu-lations as shown in Fig. 1(b). Note that d (cid:44) N . In crossing the transition r ± jumps from zero to a nearlyenergy-independent constant. Positive (negative) vortices ac-cumulate at the point r ± gradually as energy increases, whichyields continuous growth of the dipole moment (cid:104) D (cid:105) .In Fig. 1, (c) and (d) we compare MC-sampling and meanfield theory for the density di ff erence σ at di ff erent ener-gies. For E − E c = . , .
1, the mean field result for σ isevaluated perturbatively via Eq. (14) by plugging in the zeromode Eq. (28). Agreement between the spatial distributionsis very good near the clustering transition where a single un-stable mode is dominant. At the high energy E = .
02, wecompare MC-sampling with the supercondensation asympoticform [Eq. (41)] that ignores the interaction between clusters.The MC-sampling reveals additional elliptical distortion ofeach cluster.
C. Order Parameter Dependence on E for Increasing N We also extract the average dipole moment D as a func-tion of energy for a range of N . We first estimate the criticalenergy E c for each N using a simple fitting procedure in theregion close to the transition. Given discrete samples of thedipole moment D i = D ( E i ), and its variance, obtained from -3.0-2.5-2.0-1.5-1.0-0.50.00.5-6 -5 -4 -3 -2 -1 0 1 N = N = N = N = N = l n ( D ) ln( E − E c ) D E FIG. 2. (Color online) Logarithmic plot of average dipole moment D against the energy above the critical energy, E − E c . Solid lineshows the predicted mean-field scaling | E − E c | / . Inset shows anon-logarithmic plot of D against E . MC simulations conducted at energies E i , we minimize χ = (cid:88) i ∈R (cid:16) D i − √ E i − E c (cid:17) ( ∆ D i ) , (50)where R is the set of i such that 0 < E i − E c < .
1. Thischoice of energy range leads to at least 7 points in the set R for every value of N . This procedure yields the values of E c reported in Table I. We can see that E c ( N ) →
0, the meanfield prediction, as N increases. We calculate an estimate ofthe uncertainty in E c ( N ) from the χ + ∆ E c ( N ).The error in E c decreases with increasing N , corresponding toincreasing goodness of fit [77].Having estimated E c independently for each N , in Fig. 2 weplot dipole moment D against E − E c on a log scale. The scal-ing in the transition region becomes very close to | E − E c | / as N increases, with good agreement for the scaling alreadyevident for N = D observedbelow the transition energy systematically decrease in ampli-tude with increasing N [Fig. 2 (inset)]. Also note that the low-est accessible energy of the finite-size system increases with N ; in the thermodynamic limit there are no states with E < V. DISCUSSION AND CONCLUSIONSA. Phases of the (Extended) Point-Vortex Model
In Fig. 3 we summarize the equilibrium phases of the sys-tem of point-vortices confined to a disc. As an end point ofa continuous process as energy increases, the system enters the supercondensation regime at β s = − β c , with β s < β c < ffi cult todistinguish [39]. The microcanonical approach reveals thatthe apparent coincidence of the transitions is due to numericalproximity of β s and β c . In the microcanonical ensemble, statesare parametrized by energy, giving a clear distinction betweenthe vortex clustering transition and supercondensation.Our mean field analysis gives analytical expressions for thestructure of the clustered phase and critical behavior of thedipole moment near β (cid:39) β c . For β c < β <
0, Eqs. (9), (14)only have the trivial solution, namely n ± ( r ) = n and ψ = β c < β < β > ff ective field for-mulation of the point-vortex model is the Sine-Gordon fieldtheory. At β pc > − Kosterlitz − Thouless(BKT) transition at β BKT [81–83].
B. Symmetries and Order Parameter
One should distinguish the vortex clustering transition dis-cussed in this paper from the BEC transition. The fundamen-tal degrees of freedom in our system are the vortices but notthe bosonic atoms of the host fluid. The clustering transitionis not described by a local order parameter. The macroscopicdipole moment instead describes the global spatial structureof the system. The clustered phase breaks the underlyingSO(2) symmetry as well as the Z symmetry, and we empha-size that the SO(2) symmetry breaking occurs at zero angularmomentum. There are also clustered phases which break the Z symmetry while preserving the SO(2) symmetry, carryingfinite angular momentum. The macroscopic dipole moment (cid:104) D (cid:105) = (cid:104)| (cid:80) i κ i r i |(cid:105) functions as a global order parameter as (cid:104) D (cid:105) provides a good description of the transition from the uniformphase (cid:104) D (cid:105) = (cid:104) D (cid:105) (cid:44)
0. The crit-ical exponent for the order parameter, Eq. (35), has the samevalue as the mean field XY model. However an anomalousscaling correction might be present due to a remnant of somescreening e ff ects close to the transition, which is beyond thescope of the mean field theory approach taken here. Near thetransition we expect that N ( ∆ D ) ∼ ( E − E c ) − γ as N → ∞ .While our current level of MC sampling is su ffi cient to deter-mine the scaling of (cid:104) D (cid:105) , identifying γ would require signifi- BKT c s pc Super- Condensation Clustering Transition Pair- Collapse BKT Transition
FIG. 3. (Color online) Schematic of phases of a neutral system ofpoint vortices confined to the disc geometry. Energy per vortex de-creases from left to right. In dimensionless units the supercondensa-tion temperature occurs at the universal point β s = −
2. The vortexclustering transition occurs at β c (cid:39) − .
835 [Eq. (27)]. In the posi-tive temperature regime, the pair-collapse limit is reached at β pc . Theaddition of short-range repulsion between vortices allows the point-vortex model to extend further into the positive temperature regimeto encompass the BKT transition at β BKT > β c (see text). cant further investigation.The clustering transition is not described by a local orderparameter and hence one can not apply the standard scenarioof phase transitions straightforwardly. However, the cluster-ing transition does have features of a second-order phase tran-sition: 1) it spontaneously breaks the underlying symmetries;2) there is an order parameter associated with the broken sym-metry that grows continuously as a function of the tuning pa-rameter crossing the transition; 3) the fluctuation of the orderparameter has a peak near the transition and the peak growsas the number of vortices N increases.The important role of symmetry in this system warrantssome discussion of the role of the SO(2) symmetry of the trap.We note that a very high degree of cylindrical trap symme-try was essential for preserving su ffi cient angular momentumto achieve rotating Bose-Einstein condensation [84]. Further-more, recent advances in digital-micromirror device technol-ogy [65] could allow the construction of highly circular traps.Another promising avenue may be o ff ered by using specificGauss-Laguerre modes to induce radial confinement [85], aslaser fields can be e ffi ciently shaped into specific modes thatare SO(2) symmetric. C. Conclusions
We study clustering phenomena in the neutral point-vortexmodel confined to a disc, in the microcanonical ensemble.Utilizing concepts from the standard theory of phase transi-tions, we give a clear picture of the clustering transition ofquantum vortices, including the transition temperature, thepower law scaling of the order parameter near the transition,and the distinction between clustering and supercondensa-tion. The clustering transition breaks the underlying Z andSO(2) symmetries, exhibiting features of a mean field XY-type second-order phase transition. The mean field theory isfound to agree well with microcanonical Monte Carlo simu- lations, consistent with the dominance of long-range interac-tions in the clustered regime. By uncovering the nature of theclustering transition in a bounded domain, our theory gives atreatment of the transition relevant for experimental realiza-tion in Bose-Einstein condensates, and suggests future direc-tions for testing ensemble equivalence in systems with long-range interactions. Our approach could be extended to arbi-trary confining geometries, and provides a starting point forstudying vortex clustering phenomena in quantum fluid sys-tems that are strongly forced and damped [86, 87]. ACKNOWLEDGMENTS
We acknowledge M. M¨uller, P. B. Blakie, T. P. Simula, S.A. Gardiner, S. Flach, S. Denisov, T. Li and A. M. Mateo foruseful discussions. TPB acknowledges support from the JohnTempleton Foundation via the Durham Emergence Project.ASB is supported by a Rutherford Discovery Fellowship ad-ministered by the Royal Society of New Zealand.
Appendix A: Specific Heat Near the Transition
In order to investigate the microcanonical specific heat atthe transition point, we need to study higher-order perturba-tion theory in (cid:15) . Similar analyses for single charge vorticescan be found in Ref. [18]. In the following we always usethe fact that at the transition point σ = α = δα = δψ for l = δψ (cid:39) (cid:15)ψ + ψ (2)1 + ψ (3)1 , where ψ ( i ) ( i = ,
3) is of order O ( (cid:15) i ).To the O ( (cid:15) ), ∇ δψ = − π [ δ n + ( r ) − δ n − ( r )] = π (cid:104) n ( ψδβ + βδψ + r δα ) − n + δγ + + n − δγ − (cid:105) (A1)becomes L ψ (2)1 = π (cid:18) − n + δγ + + n − δγ − (cid:19) , (A2)where L ≡ ∇ − π n β . Multiplying ψ ( r , θ ) = ˆ ψ ( r ) cos( θ ) toEq. (A2) and integrating over space, the right-hand side of theequation is trivially zero and the left-hand side of the equationis also zero because (cid:82) d r L ψ (2) ψ = (cid:82) d r ψ (2) L ψ = η which satisfies (cid:32) r ddr r ddr − π n β c (cid:33) η = π n (A3)with boundary conditions η ( r = = η (cid:48) ( r = =
0. Thesolution of Eq. (A3) reads η ( r ) = [2 β c J ( j , )] − J ( j , r ) − / β c . Then ψ (2)1 can be expressed as ψ (2)1 = η ( − δγ + + δγ − ) . (A4)0Conservation of the number of vortices leads to (cid:90) d r ( δ n + − δ n − ) = , (A5) (cid:90) d r ( δ n + + δ n − ) = . (A6)At O ( (cid:15) ), the constraint Eq. (A5) becomes (cid:90) d r ∇ δψ = (cid:90) d r ∇ ψ (2)1 = , (A7)and then we have (cid:90) π d θ ∂ψ (2)1 ∂ r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r = = , (A8)which is satisfied automatically. Then Eq. (A6) takes the fol-lowing form (cid:90) d r (cid:18) n + δγ + + n − δγ − + n β (cid:15) ψ (cid:19) = . (A9)At O ( (cid:15) ), the constraint 0 = δ M = − π − (cid:82) d r ψ (2) requiresthat δγ − = δγ + , and therefore ψ (2)1 =
0. In order to obtain the β − E relation, we need to expandEq. (A1) to O ( (cid:15) ): L ψ (3)1 = π n (cid:15)ψ (cid:20) δβ + β c δγ + + β c δγ − + (cid:15) β c ψ (cid:21) . (A10)Multiplying ψ ( r , θ ) to Eq. (A10) and integrating overspace, the left-hand side of the equation is zero because (cid:82) d r L ψ (3) ψ = (cid:82) d r ψ (3) L ψ =
0. Hence we have (cid:15) (cid:90) d r ψ (cid:2) δβ + β c ( δγ + + δγ − ) (cid:3) + (cid:15) β c (cid:90) d r ψ = . (A11)Combining Eq. (A9) and Eq. (A11), we obtain δβ = ˜ C δ E ,where˜ C = β c (cid:20) n (cid:18) (cid:90) d r ψ (cid:19) − (cid:90) d r ψ (cid:21)(cid:18) (cid:90) d r ψ (cid:19) − (cid:39) − . . (A12)Hence we arrive at Eq. (36). [1] G. L. Eyink and K. R. Sreenivasan, Rev. Mod. Phys. , 87(2006).[2] E. D. Siggia and H. Aref, Phys. Fluids , 171 (1981).[3] G. Bo ff etta and R. E. Ecke, Annu. Rev. Fluid Mech. , 427(2012).[4] R. H. Kraichnan, J. Fluid Mech. , 155 (1975).[5] R. H. Kraichnan and D. Montgomery, Rep. Prog. Phys. , 547(1980).[6] B. Turkington, A. Majda, K. Haven, and M. DiBattista, P. Nat.Acad. Sci. , 12346 (2001).[7] F. Bouchet and J. Sommeria, J. Fluid Mech. , 165 (2002).[8] D. S. Choi and A. P. Showman, Icarus , 597 (2011).[9] L. Onsager, Il Nuovo Cimento , 279 (1949).[10] G. Joyce and D. Montgomery, J. Plasma Phys. , 107 (1973).[11] D. Montgomery and G. Joyce, Phys. Fluids , 1139 (1974).[12] S. F. Edwards and J. B. Taylor, Proc. Roy. Soc. A: Math., Phys.and Eng. Sci. , 257 (1974).[13] J. H. Williamson, J. Plasma Phys. , 85 (1977).[14] D. L. Book, S. Fisher, and B. E. McDonald, Phys. Rev. Lett. , 4 (1975).[15] Y. B. Pointin and T. S. Lundgren, Phys. Fluids , 1459 (1976).[16] L. J. Campbell and K. O’Neil, J. Stat. Phys. , 495 (1991).[17] R. A. Smith, Phys. Rev. Lett. , 1479 (1989).[18] R. A. Smith and T. M. O’Neil, Phys. Fluids B , 2961 (1990).[19] Y. Yatsuyanagi, Y. Kiwamoto, H. Tomita, M. Sano, T. Yoshida,and T. Ebisuzaki, Phys. Rev. Lett. , 054502 (2005).[20] J. G. Esler and T. L. Ashbee, J. Fluid Mech. , 275 (2015).[21] J. Cardy, Scaling and Renormalization in Statistical Physics (Cambridge University Press, 1996); N. Goldenfeld,
Lec-tures on Phase Transitions and the Renormalization Group (Addison-Wesley, 1992); M. Kardar,
Statistical Physics ofFields (Cambridge University Press, Cambridge, 2009).[22] M. Cross and P. Hohenberg, Rev. Mod. Phys. , 851 (1993). [23] A. J. Chorin, Vorticity and Turbulence , Applied MathematicalSciences, Vol. 103 (Springer, New York, NY, 1994); L. M. Pis-men,
Vortices in Nonlinear Fields , 1st ed., From Liquid Crys-tals to Superfluids, from Non-equilibrium Patterns to CosmicStrings (Oxford University Press, New York, 1999).[24] E. M. Purcell and R. V. Pound, Phys. Rev. , 279 (1951).[25] A. S. Oja and O. V. Lounasmaa, Rev. Mod. Phys. , 1 (1997).[26] P. Medley, D. M. Weld, H. Miyake, D. E. Pritchard, and W. Ket-terle, Phys. Rev. Lett. , 195301 (2011).[27] S. Braun, J. P. Ronzheimer, M. Schreiber, S. S. Hodgman,T. Rom, I. Bloch, and U. Schneider, Science , 52 (2013).[28] L. D. Carr, Science , 42 (2013).[29] A. L. Fetter, Phys. Rev. , 100 (1966).[30] S. J. Rooney, P. B. Blakie, B. P. Anderson, and A. S. Bradley,Phys. Rev. A , 023637 (2011).[31] T. Aioi, T. Kadokura, T. Kishimoto, and H. Saito, Phys. Rev. X , 021003 (2011).[32] G. Moon, W. J. Kwon, H. Lee, and Y.-i. Shin, Phys. Rev. A ,051601 (2015).[33] T. W. Neely, E. C. Samson, A. S. Bradley, M. J. Davis, andB. P. Anderson, Phys. Rev. Lett. , 160401 (2010).[34] T. W. Neely, A. S. Bradley, E. C. Samson, S. J. Rooney, E. M.Wright, K. J. H. Law, R. Carretero-Gonz´alez, P. G. Kevrekidis,M. J. Davis, and B. P. Anderson, Phys. Rev. Lett. , 235301(2013).[35] A. S. Bradley and B. P. Anderson, Phys. Rev. X , 041001(2012).[36] M. T. Reeves, T. P. Billam, B. P. Anderson, and A. S. Bradley,Phys. Rev. Lett. , 104501 (2013).[37] T. P. Billam, M. T. Reeves, B. P. Anderson, and A. S. Bradley,Phys. Rev. Lett. , 145301 (2014).[38] A. White, C. F. Barenghi, and N. P. Proukakis, Phys. Rev. A , 013635 (2012). [39] T. P. Simula, M. J. Davis, and K. Helmerson, Phys. Rev. Lett. , 165302 (2014).[40] M. T. Reeves, T. P. Billam, B. P. Anderson, and A. S. Bradley,Phys. Rev. A , 053631 (2014).[41] G. W. Stagg, N. G. Parker, and C. F. Barenghi, J. Phys. B: At.Mol. Opt. Phys. , 095304 (2014).[42] A. Lucas and P. Sur´owka, Phys. Rev. A , 053617 (2014).[43] A. T. Powis, S. J. Sammut, and T. P. Simula, Phys. Rev. Lett. , 165303 (2014).[44] T. P. Billam, M. T. Reeves, and A. S. Bradley, Phys. Rev. A ,023615 (2015).[45] G. W. Stagg, A. J. Allen, N. G. Parker, and C. F. Barenghi,Phys. Rev. A , 013612 (2015).[46] M. T. Reeves, T. P. Billam, B. P. Anderson, and A. S. Bradley,Phys. Rev. Lett. , 155302 (2015).[47] A. Skaugen and L. Angheluta, Phys. Rev. E , 032106 (2016).[48] A. Skaugen and L. Angheluta, Phys. Rev. E , 042137 (2016).[49] H. Salman and D. Maestrini, arXiv:1607.00092 (2016).[50] W. J. Kwon, G. Moon, J.-y. Choi, S. W. Seo, and Y.-i. Shin,Phys. Rev. A , 063627 (2014).[51] K. E. Wilson, Z. L. Newman, J. D. Lowney, and B. P. Anderson,Phys. Rev. A , 023621 (2015).[52] W. J. Kwon, S. W. Seo, and Y.-i. Shin, Phys. Rev. A , 033613(2015).[53] J. H. Kim, W. J. Kwon, and Y. Shin, arXiv:1607.00092 (2016).[54] N. Parker, N. P. Proukakis, and C. F. Barenghi, Phys. Rev. Lett. , 160403 (2004).[55] N. G. Parker, A. J. Allen, C. F. Barenghi, and N. P. Proukakis,Phys. Rev. A , 013631 (2012).[56] R. Numasato and M. Tsubota, J. Low Temp. Phys. , 415(2010).[57] R. Numasato, M. Tsubota, and V. S. L’vov, Phys. Rev. A ,063630 (2010).[58] P. M. Chesler, H. Liu, and A. Adams, Science , 368 (2013).[59] B. Nowak, S. Erne, M. Karl, J. Schole, D. Sexty, and T. Gasen-zer, arXiv:1302.1448 (2013).[60] T. Langen, T. Gasenzer, and J. Schmiedmayer, J. Stat. Mech. , 064009 (2016).[61] G. B. Hess, Phys. Rev. , 189 (1967).[62] P. K. Newton, The N-Vortex Problem , Applied MathematicalSciences, Vol. 145 (Springer New York, New York, NY, 2001).[63] K. Henderson, C. Ryu, C. MacCormick, and M. G. Boshier,New J Phys , 043030 (2009).[64] A. L. Gaunt, T. F. Schmidutz, I. Gotlibovych, R. P. Smith, and Z. Hadzibabic, Phys. Rev. Lett. , 200406 (2013).[65] G. Gauthier, I. Lenton, N. M. Parry, M. Baker, M. J. Davis,H. Rubinsztein-Dunlop, and T. W. Neely, arXiv:1607.00092(2016).[66] H. Aref, P. L. Boyland, M. A. Stremler, and D. L. Vainchtein, Fundamental Problematic Issues in Turbulence , edited byA. Gyr, W. Kinzelbach, and A. Tsinober (Birkh¨auser, Basel,1999).[67] G. L. Eyink and H. Spohn, J. Stat. Phys. , 833 (1993).[68] E. Caglioti, P. L. Lions, C. Marchioro, and M. Pulvirenti,Comm. Math. Phys. , 229 (1995).[69] C. C. Lin, P. Nat. Acad. Sci. , 570 (1941).[70] A. C. Ting, H. H. Chen, and Y. C. Lee, Phys. Rev. Lett. ,1348 (1984).[71] A. C. Ting, H. H. Chen, and Y. C. Lee, Physica D , 37 (1987).[72] B. N. Kuvshinov and T. J. Schep, Phys. Fluids , 3282 (2000).[73] D. Gurarie and K. W. Chow, Phys. Fluids , 3296 (2004).[74] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series,and Products , 4th ed. (Academic Press, New York, 1980).[75] M. Creutz, Phys. Rev. Lett. , 1411 (1983).[76] F. Wang and D. P. Landau, Phys. Rev. Lett. , 2050 (2001).[77] I. Hughes and T. Hase, Measurements and their Uncertainties:A practical guide to modern error analysis (Oxford UniversityPress, Oxford, 2010).[78] E. H. Hauge and P. C. Hemmer, Phys. Norv. , 209 (1971).[79] F. Cornu and B. Jancovici, J. Stat. Phys. , 33 (1987).[80] F. Cornu and B. Jancovici, J. Chem. Phys. , 2444 (1989).[81] J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys. , 1181 (1973).[82] L. P. Kadano ff , J. Phys. A-Math. Gen. , 1399 (1978).[83] B. Nienhuis, J. Stat. Phys. , 731 (1984).[84] P. Haljan, I. Coddington, P. Engels, and E. Cornell, Phys. Rev.Lett. , 210403 (2001).[85] A. Ramanathan, K. C. Wright, S. R. Muniz, M. Zelan, W. T.Hill, C. J. Lobb, K. Helmerson, W. D. Phillips, and G. K.Campbell, Phys. Rev. Lett. , 130401 (2011).[86] K. G. Lagoudakis, M. Wouters, M. Richard, A. Baas, I. Caru-sotto, R. Andre, L. E. S. I. Dang, and B. Deveaud-Pledran, Nat.Phys. , 706 (2008).[87] G. I. Harris, D. L. McAuslan, E. Sheridan, Y. Sachkou,C. Baker, and W. P. Bowen, Nat. Phys. (2016),10.1038 //