Energy spectra of vortex distributions in two-dimensional quantum turbulence

We theoretically explore key concepts of two-dimensional turbulence in a homogeneous compressible superfluid described by a dissipative two-dimensional Gross-Pitaeveskii equation. Such a fluid supports quantized vortices that have a size characterized by the healing length $\xi$. We show that for the divergence-free portion of the superfluid velocity field, the kinetic energy spectrum over wavenumber $k$ may be decomposed into an ultraviolet regime ($k\gg \xi^{-1}$) having a universal $k^{-3}$ scaling arising from the vortex core structure, and an infrared regime ($k\ll\xi^{-1}$) with a spectrum that arises purely from the configuration of the vortices. The Novikov power-law distribution of intervortex distances with exponent -1/3 for vortices of the same sign of circulation leads to an infrared kinetic energy spectrum with a Kolmogorov $k^{-5/3}$ power law, consistent with the existence of an inertial range. The presence of these $k^{-3}$ and $k^{-5/3}$ power laws, together with the constraint of continuity at the smallest configurational scale $k\approx\xi^{-1}$, allows us to derive a new analytical expression for the Kolmogorov constant that we test against a numerical simulation of a forced homogeneous compressible two-dimensional superfluid. The numerical simulation corroborates our analysis of the spectral features of the kinetic energy distribution, once we introduce the concept of a {\em clustered fraction} consisting of the fraction of vortices that have the same sign of circulation as their nearest neighboring vortices. Our analysis presents a new approach to understanding two-dimensional quantum turbulence and interpreting similarities and differences with classical two-dimensional turbulence, and suggests new methods to characterize vortex turbulence in two-dimensional quantum fluids via vortex position and circulation measurements.


I. INTRODUCTION
Turbulence in three-dimensional (3D) classical fluids is associated with a cascade of energy from large length scales defined by the details of an energy-forcing mechanism, to small length scales where viscous damping removes kinetic energy from the fluid. This range of length scales, and the range of associated wavenumbers k, define the inertial range of energy flux [1]. As shown by Kolmogorov in 1941 [2], the energy cascade corresponds to a kinetic energy spectrum that is proportional to k −5/3 in the inertial range. Turbulence in a 3D fluid is also often associated with the decay of large patches of vorticity into ever smaller regions of vorticity; this Richardson cascade of vorticity provides an important visual picture of the fluid dynamics involved in 3D turbulence [3].
Remarkably, two-dimensional (2D) incompressible classical fluids exhibit very different turbulent flow characteristics due to the existence of an additional inviscid invariant: in the absence of forcing and dissipation, the total enstrophy [4] of a 2D fluid is conserved in addition to the fluid's kinetic energy [5][6][7][8][9]. The fluid dynamics during forced 2D turbulence are highly distinctive when compared with 3D flows: rather than decaying into smaller patches, vorticity aggregates into larger coherent rotating structures [10] (see [11] for a more detailed picture in terms of turbulent stress imposed on smallscale vortices). Accompanying these 2D fluid dynamics is an inverse energy cascade in which energy moves from a small forcing scale to progressively larger length scales, defining an inertial range for energy transport with energy flux in a direction opposite that of 3D turbulence. Eventually, energy is transported into flows characterized by length scales that are on the order the system size [12], for which dissipation may occur. Additionally, there is an enstrophy cascade, in which enstrophy is transported from the forcing scale to progressively smaller scales [13]. Thus in 2D turbulence, the kinetic energy distribution contains at least these two distinctly different spectral regimes.
Quantum turbulence [14] involves chaotic flow in a superfluid [15][16][17][18][19] and is often associated with a random vortex tangle in 3D [15]. In general, the quantization of circulation strongly constrains the velocity fields allowed in quantum turbulence, which must be irrotational everywhere within the fluid, yet inertial ranges with k −5/3 spectral dependence are still found in 3D quantum turbulence [20]. In an incompressible superfluid (such as HeII), the vortex core diameter can be neglected for all practical purposes, inspiring the study of point-vortex models of superfluid dynamics. Such a model was used by Onsager to first predict the aggregation of vortices within inviscid 2D fluids, and was the context for his prediction of the quantization of vortex circulation in a superfluid [21]. Despite the historical importance of this approach in stimulating advances in 2D classical turbulence [22], characteristics of 2DQT remain little known, due in part to the difficulty of achieving the necessary 2D confinement for incompressible superfluids. The increasing relevance of 3D turbulence concepts to dilute-gas Bose-Einstein condensate (BEC) experiments [23][24][25] and recent theoretical work on 2DQT [19,[26][27][28][29][30][31] have highlighted the need for a treatment of turbulence in 2D superfluid systems that incorporates the concept of compressibility from the outset. Motivated by recent experimental demonstrations of the confinement needed to study 2DQT in a dilute-gas BECs [32,33] our aim in the present paper is to present a new approach to solving some of the open problems of 2DQT in the context of such a system.
In a BEC, the vortex core size is non-negligible, and stems from the healing length ξ, a scale of fundamental importance in BEC dynamics that is typically about two orders of magnitude smaller than the system size [34]. Compressibility also allows for a rich array of physical phenomena in these superfluids; in particular, a vortex dipole [32] (comprised of two vortices of opposite sign of circulation) can recombine, releasing vortex energy as a burst of acoustic waves. The opposite process of vortex dipole generation from sound may also efficiently occur.
Recent theoretical studies of decaying quantum turbulence in 2D BECs have shown that when the vortex dipole annihilation process is dominant it sets up a direct cascade of energy over the scales associated with the dipole decay, suggesting that this annihilation mechanism could prohibit an inverse energy cascade from occurring in a compressible superfluid [29,30]. Moreover, enstrophy in a quantum fluid is associated with the number of vortex cores; if vortices annihilate, then enstrophy may not be conserved, bringing into question the existence of energy and enstrophy cascades in 2DQT, and the universal nature of 2D turbulence and its correspondence with 2DQT.
The general characteristics of 2D quantum turbulence in compressible quantum fluids, including the capacity for these systems to show an inverse energy cascade, enstrophy conservation, and vortex aggregation have thus remained largely unknown. However, a recent study of the formation of vortex dipoles during the breakdown of superfluid flow around an obstacle in a highly oblate BEC experimentally and numerically observed aggregation of like-sign vortices into larger-scale coherent structures [32], and found time scales over which vortex number and hence enstrophy may remain constant. The vortex clustering effect inhibits the dipole-decay mechanism by keeping vortices distant from antivortices (vortices of opposite circulation), and suggests that an inverse cascade might be observed under the right conditions of forcing. System dynamics consistent with the existence of an inverse energy cascade were indeed found in a recent study of forced 2DQT in a BEC [33].
In this article we address 2D quantum turbulence in a compressible quantum fluid from an analytical perspective. We determine the kinetic energy spectra of vortex distributions in a homogeneous compressible superfluid in a quasi-exact manner via an analytical treatment of the physics of the vortex core. We are thus able to study the properties of vortex configurations and their resulting spectra in BEC. We develop a technique to sample spatially localized vortex distributions with power-law behavior over a well defined scale range. We are thus able to identify the conditions for an inertial range in fully polarized and neutral systems. A polarized cluster is sampled using a specific exponent for the vortex locations relative to the cluster center, which is size and scale dependent. The specific radial exponent is shown to determine the velocity distribution in the classical limit and we thus identify an expanding inertial range with a steepening velocity distribution.
Making use of the universality of the spectral region generated by the vortex core, we identify an analytical form of the Kolmogorov constant that we test against dynamical simulations of the damped GPE. The derivation of the Kol-mogorov constant occurs for a highly idealized vortex distribution. Thus the complex flows generated by real forcing require that we introduce a new parameter called the clustered fraction, and modify our ansatz to account for imperfect clustering, based on the universality of the Kolmogorov constant. The modified ansatz agrees well with the numerical simulations of grid turbulence, supporting our analytical identification of the Kolmogorov constant.

II. BACKGROUND
The starting point for much of BEC theory is the Gross-Pitaevskii equation (GPE), which provides a capable description of trapped Bose-Einstein condensates at zero temperature [34]. Our model, outlined below, consists of a damped GPE (dGPE) description of a finite-temperature BEC which can be derived from the stochastic GPE theory [35,36]. In this section we develop a link between the dGPE and the classical Navier-Stokes equation, identifying a quantum viscosity arising from the damping. The corresponding Reynolds number is defined in direct analogy with classical fluids. We then state some key properties of a single quantum vortex, and define the decomposition of kinetic energy into its compressible and incompressible components.

A. Damped Gross-Pitaevskii theory
The damped Gross-Pitaevskii equation of motion for the quantum fluid wave function ψ(r, t) has been obtained phenomenologically [37], within ZNG theory [38], and via a microscopic reservoir theory [35,39,40], and we will consider it within the context of the latter framework, for which the full equation of motion is the Stochastic Projected Gross-Pitaevskii equation (SPGPE). The SPGPE is derived by treating all atoms with energies above an appropriately chosen energy cutoff ǫ cut as thermalized, leading to a grandcanonical description of the atoms below ǫ cut . A dimensionless temperature-dependent rate γ describes Bose-enhanced collisions between thermal reservoir atoms and atoms in the BEC. Neglecting the noise, we obtain the equation of motion for the condensate wave function (in the frame rotating with the chemical potential µ) For atoms of mass m in an external potential V(r, t), the operator L gives the GPE evolution: where the interaction parameter is g = 4π 2 a/m, for s-wave scattering length a. This equation of motion has been used extensively in previous studies of vortex dynamics [38,39,41,42] and provides a capable description of dynamical BEC phenomena. In general the damping parameter is small (γ ≪ 1), and it is typically much smaller than any other rates characterizing the evolution. Defining the Gross-Pitaevskii Hamiltonian and condensate atom number the equation of motion (1) evolves the grand-canonical Hamil- The stationary solution minimizing K C is the ground state satisfying µψ 0 (r) ≡ Lψ 0 (r). This is a consequence of the nonlinear form of the damping in (1). The damping term arises from collisions between high energy atoms that lead to a Boseenhanced growth of the matter wave field, with instantaneous energy determined by L. The equation of motion thus describes a system coupled to a thermal reservoir in the chosen frame of reference.
The SPGPE provides a rigorous framework for the dGPE derivation, originating from a microscopic treatment of the reservoir interaction. In particular, γ can be calculated explicitly [36] for a system with well-defined reservoir parameters µ, T, and ǫ cut , i.e. a system close to thermal equilibrium. In essence it is computed via a reduced Boltzmann collision integral that accounts for all irreversible s-wave interactions that can change the condensate population by interacting with the thermal cloud. If the thermal cloud is 3D (i.e. β −1 ≡ k B T is greater than the potential well mode spacing in each spatial dimension) the damping takes the explicit form where Φ[z, s, α] is the Lerch transcendent, and with λ dB ≡ 2π 2 /mk B T the thermal deBroglie wavelength. The dimensionless rate γ 0 provides a useful estimate of the full damping strength when the cutoff ǫ cut is unknown. Equation (6) is independent of position, and valid over the region V(r, t) ≤ 2ǫ cut /3, provided the potential can be treated semiclassically [36]. The summation gives Bose-enhancement corrections due to the Bose-Einstein distributed reservoir atoms, and is typically of order 1-20 in SPGPE simulations with a consistently determined energy cutoff [43]. Typicaly γ ∼ 5 × 10 −4 in 87 Rb experiments [42,44].

B. Heuristic derivation of a quantum Reynolds number
In this section we consider the role of dissipation within the dGPE description, and show how to recover the celebrated Navier-Stokes equation (NSE). In doing so we find an explicit expression for the viscosity which has a microscopic quantum origin, stemming from s-wave scattering of incoherent reservoir particles with a coherent superfluid. While not offering a practical reformulation (the GPE and its generalizations are capable numerical workhorses), this indicates a connection between the dGPE and the NSE in the hydrodynamic regime, allowing the identification of a parameter analogous to the kinematic viscosity of classical fluids.
The fluid dynamics interpretation of the Gross-Pitaevskii equation is based on the Madelung transformation, which we now apply to the damped GPE (1), writing ψ(r, t) = ρ(r, t) exp [iΘ(r, t)], where ρ(r, t) is the number density of the superfluid (number of atoms per unit volume), and Θ(r, t) is the macroscopic phase of the quantum fluid. The velocity is then given by v(r, t) = ∇Θ(r, t)/m. The resulting equations of motion (with implicit t and r dependence) for density and velocity are then given by where an effective potential U eff is defined as The last term is called the quantum pressure, which is very small except where ρ changes sharply, such as near vortex cores. By neglecting this term in the absence of dissipation we are considering the so-called hydrodynamic regime. We now consider the γ term in (9) Note that v · ∇ρ ≡ 0 for an isolated quantum vortex. In the absence of acoustic energy this will also be a good approximation for a system of vortices provided their cores are well separated, since the density gradient of each vortex is localized to a region where the velocity is dominated by the single-vortex velocity field. It should thus be a reasonable approximation to neglect the second term in (11). In a superfluid the curl term in the expansion ∇(∇ · v) = ∇ × (∇ × v) + ∇ 2 v may also be consistently neglected away from vortex cores; similarly we neglect the curl term in when taking the gradient of (10). We then find that (9) reduces to a quantum Navier-Stokes equation for the velocity field: where the kinematic quantum viscosity is in analogy with classical fluids. In this regime, (12) is coupled to the continuity equation with hydrodynamic potential The source term in (14) drives the system towards particlenumber equilibrium with the reservoir. In the Thomas-Fermi regime µ − U H (r) ≈ 0, restoring approximate particle-number conservation.
Making use of (7), we can give an order-of-magnitude estimate for the viscosity We can also estimate a quantum Reynolds number as for BEC flow with characteristic speed U and length scale L. We can write the quantum Reynolds number as We note that temperature only enters the expression through the deBroglie wavelength of the matter wave field, in the ratio λ dB /a, which is typically very large for a BEC. Note that strong scattering corresponds to strong damping, and hence low Reynolds number. A large deBroglie wavelength corresponds to a relatively cold system, which is hence expected to be weakly damped and have a high Reynolds number. We thus have a dimensionless ratio in the form where /L is interpreted as the quantum momentum associated with the transverse length scale of the flow. A concrete example is provided by a recent experimental study of 2DQT generated by stirring a highly oblate, toroidally confined BEC [33]. The initial system consists of ∼ 2.6 × 10 6 atoms of 87 Rb at a temperature of ∼ 100nK. Using these numbers in our analysis the dimensionless damping parameter γ 0 ∼ 6 × 10 −4 gives a kinematic quantum viscosity ν 0 q ∼ 6 × 10 −2 µm 2 s −1 . The trapping potential confines the flow to an annular channel of width L ∼ 30 µm. The nominal flow speed can be estimated from numerical simulations of the dGPE [33], giving a value U ∼ 5 µm s −1 as the peak value occurring in the bulk flow during the stirring sequence. These values give an estimate Re 0 q ∼ 600. We can alternatively estimate a Reynolds number at the scale of the forcing in the experiment, which is of order of the size of the vortex dipoles nucleated, d ∼ 10 ξ, with healing length ξ ∼ 0.5 µm. Such dipoles have a characteristic speed v d ∼ 146 µms −1 , for which we estimate the Reynolds number of the forcing scale as ∼ v d d/ν 0 q = 6.2×10 3 . These large values suggest that turbulent flow in a finite-temperature BEC may exist across a wide range of length scales if the analogy is made with the classical Reynolds numbers that correspond to turbulent flow [1]. This interpretation is broadly consistent with the experimental and numerical observations of chaotic vortex dynamics [33].
We emphasize that the quantum Reynolds number estimate proposed here is applicable to a finite-temperature weakly interacting superfluid and may provide a general condition in analogy with classical fluids that is independent of dimension. However, taking the zero-temperature limit gives an infinite value for Re 0 q , and in this regime the superfluid dissipation stems from vortex reconnections (or annihilation in two dimensions) and coupling to the sound field [45][46][47]. The detailed description of criteria for superfluidity in the cross-over from the zero-temperature to high-temperature regimes is an open problem [48]. Furthermore, due to the reversed direction of energy transfer in 2DQT, the scales of interest have to be reexamined; we do not pursue this here. Our aim is to establish a conceptual link between the dGPE and the NSE, given by Eq. (12). In doing so we have shown how to identify the equivalent viscosity in a finite temperature BEC.

C. Two-dimensional vortex wavefunction
In the remainder of this work we limit our analysis to homogeneous compressible quantum fluids in two dimensions, and redefine our spatial and velocity coordinates accordingly: r = (x, y) = r (cos θ, sin θ) and v = (v x , v y ). We thus confine our attention to the regime of an effective 2D GPE, with modified interaction parameter. While 2D BEC systems can be created through extremely tight confinement in one dimension, a regime of effective 2D vortex dynamics can also be accessed in less oblate systems, giving a 2D analysis wider applicability. For example, although the BECs of References [32] and [33] were three dimensional, the confinement along one dimension was strong enough to limit vortex motion to a plane and suppress vortex bending and tilting away from the tighttrapping direction. Aspects of BEC dimensionality in regards to vortices and Kelvin waves were analyzed in [44], further indicating that sufficiently oblate 3D BECs may be considered 2D as far as vortex dynamics and turbulence are concerned. At the same time, such systems can remain far enough away from the quasi-2D limits in which a Berezhkinski-Kosterlitz-Thouless (BKT) transition has been observed [49], and BKT physics may thus be neglected.
For our analysis of kinetic energy spectra, we require certain properties of a quantized vortex, namely the asymptotic character of the wavefunction for large and small length scales. The Gross-Pitaevskii equation describing the homogeneous (V = 0) 2D Bose gas is obtained from (1) by taking γ = 0 and using an interaction parameter g 2 = g/l where l is the characteristic thickness of the 3D system [30]: For example, in a system with harmonic trapping in the zdirection characterized by trapping frequency ω z , the length scale is l = √ 2πl z where l z = /mω z is the z-axis harmonic oscillator length, and the confinement is assumed sufficient to put the wavefunction into the z-direction single-particle ground state.
For solutions with chemical potential µ containing a single vortex at the origin (with circulation normal to the plane of the quantum fluid) we can write [50] where ξ = /mc is the healing length for speed of sound c = µ/m, and n 0 = µ/g 2 is the 2D particle density for r ≫ ξ and is taken to be a constant. The vortex radial amplitude function The boundary conditions are χ(0) = 0, and the derivative χ ′ ≡ dχ/dσ evaluated at σ = 0 must be chosen such that it is consistent with χ(∞) = 1 and χ ′ (∞) = 0. The value is determined numerically to be Λ = 0.8249 . . . . The state (21) has the velocity field of a quantum vortex v(r) = mr (∓ sin θ, ± cos θ).

D. Kinetic energy decomposition
We make use of the decomposition of the kinetic energy into compressible and incompressible parts [28,51]. The 2D case of the Gross-Pitaevskii energy functional (3) can be de- Respectively, these define the components of energy that can be attributed to kinetic energy, potential energy, interaction energy, and quantum pressure. We are interested in the kinetic energy, E K . We define a density-weighted velocity field u(r, t) ≡ ρ(r, t)v(r, t), then decompose this into u(r, t) = u i (r, t) + u c (r, t), where the incompressible field u i satisfies ∇ · u i = 0, and the compressible field u c satisfies ∇ × u c = 0. We can further decompose the kinetic energy as where the portion of E K attributed to compressible or incompressible kinetic energy is defined as The compressible component is attributed to the kinetic energy contained in the sound field, while the incompressible part gives the contribution from quantum vortices. Our analysis below only involves E i . Because we focus on vortex configurations at instants in time, we drop the explicit time dependence from the remainder of our expressions.
In k−space, the total incompressible kinetic energy E i is given by The one-dimensional spectral density in k-space is given in polar coordinates by integrating over the azimuthal angle to give which, when integrated over all k, gives the total incompressible kinetic energy

A. Incompressible kinetic energy spectrum of a vortex
We now consider the kinetic energy spectrum of a single quantum vortex in a 2D BEC. For an arbitrary wavefunction the decomposition into compressible and incompressible parts must be performed prior to carrying out the transformation to the spectral representation. However, for a quantum state containing a single vortex and no acoustic energy [i.e. the single vortex wavefunction ψ 1 (21)] we note that the wavefunction is automatically incompressible, i.e. the compressible part is identically zero: The first term vanishes due to the orthogonality of the density gradient and velocity of a vortex, and the second due to the form of (24). Thus the incompressible spectrum is the entire spectrum for a single quantum vortex.
For a single vortex we can thus ignore the incompressible decomposition and cast the kinetic energy spectrum in terms of the properties of the radial amplitude function χ(σ) = ρ(σξ)/n 0 obtained from (22). We have where J 0 is the zeroth-order Bessel function of the first kind. Similar analysis gives F y (k) = −F x (k). We can thus find the one-vortex spectrum [see (32)] where we define the dimensionless core spectral function and we have introduced the unit of enstrophy giving Ωξ 3 as the natural unit for the kinetic energy density. The core spectral function has the small-kξ asymptotic form For kξ ≫ 1, J 0 (kξσ) is highly oscillatory except at σ = 0 where it is unity, and the Taylor expansion of χ ′ (σ) can be truncated at zeroth order to give We thus have the asymptotic spectra for a single vortex The kξ ≪ 1 regime arises purely from the irrotational velocity field of a quantum vortex, while the kξ ≫ 1 regime is a property of the core of a compressible quantum vortex. The kξ ≫ 1 regime explicitly depends on the slope of the wavefunction at the core of a vortex. The cross-over between these regions occurs in the vicinity of kξ ≈ 1, hence we take kξ = 1 as distinguishing the infrared (kξ < 1) and ultraviolet (kξ > 1) regimes in the remainder of our analysis. The scale kξ = 1 thus serves to define an important length scale of the problem, namely l v ≡ 2πξ. In Fig. 1 we see that at this distance from the vortex core the deviation of the amplitude from the background value is very small. This is the scale beyond which the details of the core structure are no longer important in characterizing the wavefunction, or equivalently that the fluid density has approximately reached its bulk value. The irrotational velocity field in Eq. (24) is the only remaining signature of a vortex at this range from its center and beyond. We note that our derivation of the k −3 power-law stemming from the quantum vortex core structure is consistent with recent analysis of the Kelvin-wave cascade in 3D [52]. Now that we have identified the properties of a single vortex, it is natural to ask whether a unit of enstrophy can be attributed to a single quantum vortex, and to compare this with the quantity defined in (37). The point-vortex model suggests that this can be done, but gives a singular result, which is nevertheless known to be proportional to the number of vortices [30]. The problem is also evident if we attempt to evaluate the enstrophy of a single vortex from the spectrum (35). Multiplying by k 2 to produce an enstrophy spectral density and integrating this over the ultraviolet regime kξ > 1, we are faced with the singular integral ∞ 1 dk/k. In this work we therefore define a new asymptotic quantity with units of enstrophy as This quantity plays a fundamental role for a compressible superfluid because it completely specifies the large-k region of the incompressible kinetic energy spectrum. Because the spectrum in this region of k-space is determined by the core structure of quantized vortices, we call this unit the onstrophy to both recall Onsager's contribution to our understanding of quantized vorticity in a superfluid and emphasize the difference between enstrophy in classical and quantum fluids. For a single quantum vortex we find, using (35), the onstrophy which differs from (37) by the factor Λ 2 = 0.6805 . . . , a property of the vortex core in a compressible superfluid.

B. Vortex wavefunction ansatz and kinetic energy spectrum
To study 2D kinetic energy spectra we will make extensive use of an algebraic ansatz for the wavefunction of a single vortex in a homogeneous superfluid. Numerical evaluation of the exact core function (36) is not straightforward due to the highly oscillatory integrand, and the need to determine the vortex amplitude χ(σ) extremely accurately over a large range of length scales. In order to accurately represent the spectrum it will be crucial that our ansatz have the correct asymptotic properties for small and large length scales described immediately above (23). Making use of the slope at the origin com-puted for the exact solution in (23), we use the ansatz wavefunction: The general form of this ansatz has been previously used to describe the shape of a vortex core [50], but here we use a length scale Λ −1 ξ that enforces matching the slope of the ansatz density distribution to the exact value at the center of the core. The state (44) has the irrotational velocity field of a quantum vortex specified in (24) and reproduces the asymptotic slope of the exact solution near the origin, as shown in Figure 1.
We now compute the kinetic energy spectrum for a single vortex by evaluating (31) using the form (44).
where I j and K j are modified Bessel functions of the first and second kind, respectively, of order j. Since |F y | 2 =|F x | 2 , we find the incompressible energy spectrum of a single vortex where and where we define The function f (z) has the following asymptotics: for z ≪ 1 The function F Λ (kξ) thus has the asymtotics which are identical to those of F(kξ). The two functions are very similar, with only small differences evident in the crossover region kξ ∼ 1, as seen in Figure 2. We use F Λ (kξ) instead of F(kξ) for describing the kinetic energy spectrum for a vortex core in the remainder of this work as it is numerically expedient and does not alter any of the physical consequences of our analysis. Towards the end of this paper we will compare !" the asymptotic results of our analysis with spectra determined numerically from forced dGPE dynamics. The spectrum of a single vortex is shown in Figure 3, and compared with the spectrum of a vortex-antivortex pair (a vortex dipole), and that for two vortices of the same sign (a vortex pair). These two-vortex spectra are analyzed in the following section.

C. Two-vortex spectra
Extension of the discussion in Section III A leads us to conclude that a wavefunction that only contains vortices (i.e., no sound field) separated by more than a few healing lengths will thus be approximately incompressible according to the decomposition. The approximation breaks down through the non-orthogonality of v and ∇ ρ(r) near a vortex core due to the velocity field induced by the other vortices. However, in the close vicinity of a vortex core, where ∇ ρ(r) is significant, the velocity is dominated by the velocity field of that vortex core. An arrangement of vortices separated by more than a few healing lengths will thus be approximately incompressible. In the following analytical treatment we will neglect any compressible part that arises from an assembly of vortices described by the ansatz (44).
A two-vortex state in a homogeneous system with no boundaries has density-weighted velocity field ρ(r)[v 1 (r) + v 2 (r)], where v j is the velocity field around vortex j = 1, 2 taken separately. If the vortex cores are separated by d ≫ ξ it will also be a very good approximation to write ρ(r) = ρ 1 (r)ρ 2 (r)/n 0 , with The density weighted two-vortex velocity field can then be written as The final term is which is only significant when considering the velocity field of one vortex in the close vicinity of the other vortex core. K 12 is therefore a negligible correction to the spectrum for widely separated vortices. This approximation may be trivially generalized to arbitrary numbers of vortices provided the cores do not overlap appreciably. This approximation is central to our treatment, as it allows recovery of familiar point-like vortex physics in the infrared regime kξ ≪ 1. Strict validity is limited to the regime where the intervortex spacing d is bounded below by ∼ l v = 2πξ, the scale at which the core structure becomes evident (see Fig. 1). We require the transform where the superscript d denotes the case of a vortex dipole, and subscript j = x, y indicates the x and y components of the density-weighted velocity fields of each vortex. As above, the subscripts 1 and 2 denote vortices 1 and 2, and the superscript i denotes that it is the incompressible, or divergence-free portion of the density-weighted field that is of interest here. To account for the opposite signs of circulation for the two vortices, without loss of generality we choose vortex 1 as positively charged, located at r 0 = (d/2)x, so that v 1 (r) = v(r − r 0 ) where v(r) is the central vortex velocity field (24). Vortex 2 has velocity field v 2 (r) = −v(r + r 0 ). We then have For a vortex dipole we can then write where F is the spectrum of a single vortex. Using (45), and the fact that F y (k) = −F x (k), we find for a vortex dipole The spectrum of a pair of vortices of the same circulation separated by d is calculated similarly to be The spectra (59) and (60) are shown in Figure 3. It is clear that for scales less than the vortex separation distance d there is interference in k-space, leading to oscillations in the spectrum. The difference between the dipole and pair is that the interference fringes are offset, and the infrared asymptotics are different, a feature we discuss further below. The spectrum of the vortex pair is clearly similar to that of the single vortex in the far infrared region, but the additional kinetic energy of the vortex pair state is observed throughout the spectrum.

D. Kinetic energy spectrum of N-vortex configurations
Extending the above analysis, for a general system of N singly quantized vortices with circulation signs κ p = ±1 located at r p the kinetic energy spectrum is given by We note the resemblance to point-vortex results which also have the Bessel function dependence [29,53]. The function F Λ (kξ) gives the incompressible limit for small k [J 0 (kr) approaches unity for small k and finite r], and introduces the physics of compressible superfluids for 1 kξ. We can write (61) as where is a purely configurational function involving summation of M = N(N − 1)/2 distinct intervortex distances. This function has the limits where the total dimensionless circulation is defined by for any contour C enclosing all N + positive and N − = N − N + negative vortices. We then find that the onstrophy for the N-vortex system is and consequently This is one of our central results: the ultraviolet regime kξ ≫ 1 has a universal asymptotic form that is independent of the vortex configuration, and that resembles the ultraviolet spectrum of classical 2D turbulence that is identified with a direct cascade of enstrophy. If we try to evaluate the classical definition of enstrophy, the result is singular, yet the onstrophy definition (42) gives a well-defined additive quantity that is singularity free and depends only on the total number of vortices in the system.

E. Infrared behavior
When Γ 0 the far infrared limit (65) gives This configuration-independent k −1 power law arises from the far-field velocity distribution of a collection of point vortices, which becomes equivalent to that of a single vortex of charge Γ at sufficiently large scales.
When Γ ≡ 0 we use the small-argument expansion J 0 (z) ≃ 1 − z 2 /4, and the asymptotic form (49), to find the kξ ≪ 1 behavior determined by the configurational information contained in the intervortex distances |r p − r q |: The simplest case involves a single vortex dipole and has only one length scale, namely the vortex separation, and the low-k form E d i (k) ≃ Ωξ 2 d 2 k/2, as shown in Fig. 3. In general, when Γ = 0 the infrared region of the spectrum is sensitive to the vortex configuration, but approaches a power-law for low-k that has a configuration-independent exponent. The linear decay of kinetic energy as k → 0 stems from the cancellation of the far-field velocity profiles for length scales greatly exceeding the largest intervortex separation in any neutral configuration of vortices.

IV. KOLMOGOROV SPECTRUM
In the previous section we obtained an explicit expression (62) for the incompressible kinetic energy spectrum that incorporates the compressible nature of individual vortex cores through the function F Λ (kξ) (derived via an ansatz for the vortex core profile), which captures the essential physics of the corresponding exact solution F(kξ) defined in (36). For both functions, point-vortex physics is recovered at large lengthscales (kξ ≪ 1). If the dynamical evolution is such that an inertial range associated with an inverse energy cascade develops, we should expect a Kolmogorov power law E i (k) ∝ k −5/3 over the inertial range. It is clear from the form of (62) that this law can only depend on the spatial configuration of the vortices. We now seek to understand the simplest situations that may show evidence for the existence of such an inertial range. We consider forcing occurring via vortex and energy injection at a forcing scale k F ∼ ξ −1 , and describe vortex configurations that do and do not lead to a Kolmogorov law for k < ξ −1 .
We now assume an idealized case in which the infrared spectrum is continuous with the universal k −3 law of the ultraviolet spectrum at the scale kξ ≈ 1. This constraint imposes a strong restriction on the infrared spectrum, completely determining its form in the case that it satisfies a power law. In this respect the universal ultraviolet region has significant physical consequences. The power-law approximation to the universal ultraviolet region based on (67) has the form The number of vortices determines the N-vortex onstrophy ζ N (67), from which the power-law approximation to the ultraviolet energy spectrum is completely determined. This power law is a very good approximation, as will be seen by sampling different vortex configurations below. The infra-red or configurational regime is then given by the kξ ≪ 1 regime of F Λ (kξ): At this point we consider the consequences of assuming that a turbulent system will have a k −5/3 law in the configurational regime, and that this power law is continuous with (71) at kξ = 1. We suppose that E N i,C (k) ∝ k −5/3 . Continuity at k = 1/ξ then requires E N i,C (1/ξ) = E N i,U (1/ξ), and gives the infrared spectrum Thus the constraint that the universal regime is continuous at the cross-over scale k = ξ −1 completely constrains the form of the configurational spectrum. Physically this may correspond to an inertial range that extends upwards from the smallest scale of the infrared region given forcing at a wavenumber k F ∼ ξ −1 . Note that this expression (73) has no reference to the signs of the vortex circulations, the degree of circulation polarization, or vortex clustering. By assuming continuity at kξ = 1 with a UV spectrum that has a universal N-vortex form, we have implicitly assumed that all N vortices are involved in determining the spectrum of the inertial range. We might expect that this will give a very good description for a completely polarized system exhibiting fully developed turbulence. When there is clustering in mixtures of different sign vortices, the spectrum may well still approach a Kolmogorov law, but there is no reason to expect that it will cross over so smoothly. We return to this problem when we compare our analysis with numerical simulations of forced turbulence in Sec. V C.
It is useful at this point to give a simplified reiteration of Novikov's argument for the power law for the vortex distribution being −1/3 in Kolmogorov turbulence [53]. To obtain power-law behavior we must consider the spectrum for a vortex distribution involving many length scales. For simplicity we assume all vortices have the same sign of circulation, κ p ≡ κ. The configuration function (63) has M = N(N − 1)/2 terms in the summation, and can be written as in terms of an average over distinct vortex separations s i . We introduce the intervortex distance distribution P(s) such that P(s)ds is the fraction of intervortex distances in the range [s, s + ds). In the continuum limit We seek a distribution P(s) that will generate a Kolmogorov law from the N-vortex spectrum (62) for scales larger than the vortex core, k ≪ ξ −1 , in the large-N regime. We then find from (72) that The scale invariance of turbulence naturally leads to the assumption that the intervortex separation distribution is a power-law P(s) ∼ s −α over the scale range of interest. The requirement of a power law in the kinetic energy then gives the scaling relation The integral is convergent for −1/2 < α < 1, allowing In particular, the universal Kolmogorov law β = 5/3 occurs for as obtained by Novikov [53] for point-vortices. We will test this scaling argument for the exponent in numerical sampling of localized vortex configurations in the following sections. Testing if this vortex separation power-law holds in simulations and experiments may give a quantitative measure of fully developed 2D turbulence in a compressible superfluid, and a way to identify the inertial range as the scale range over which this power-law can be identified. In 2D classical turbulence a k −3 region of the kinetic energy spectrum is often associated with a direct enstrophy cascade. We note that this exponent β = 3 is ruled out by (79). Hence, within this continuum analysis the k −3 power-law spectrum cannot occur in the configurational region for 2D quantum turbulence, as long as the vortex distribution follows a simple power law. This result suggests that if a direct enstrophy cascade were to occur in the configurational region of the spectrum, a different type of vortex distribution would be necessary. This makes intuitive sense, since direct enstrophy cascades may be associated with the stretching of patches of vorticity in the 2D plane. Furthermore, direct enstrophy cascades and energy spectra proportional to k −3 have been noted in simulations of superfluid helium thin films [54,55]. Nevertheless, since we have also shown that the kξ > 1 range is determined entirely by the core structure, this gives a strong indication that a direct enstrophy cascade cannot occur in this region in compressible 2DQT.
It is clear that configurations containing one or a few characteristic length scales, such as a vortex dipole or a vortex lattice, cannot lead to a power law spectrum for E N i,C (k). In the case of a vortex lattice the intervortex distance distribution has many discrete peaks [36]. The vortex dipole and a vortex pair each have a single length scale and this leads to characteristic interference fringes in the energy spectrum seen in Fig. 3.

A. Sampling spatial vortex distributions
We now test our analysis of the spectrum by numerically sampling several vortex distributions {r p , κ p } N p=1 and evaluating (62). A straightforward test of the statistical argument for the Kolmogorov power law to occur involves sampling the power law (80) and evaluating (74). Indeed, it is easily verified that this generates a k −5/3 spectrum in the configurational region. However, the connection of such a sampling to particular spatial vortex distributions is not clear, and in fact the mapping is not unique. To make this connection concrete we require a way of sampling finite, localized, spatial vortex distributions for which the vortex separations are power-law distributed.
In an ideal, infinitely extended vortex configuration exhibiting the power law (80), the system is translationally invariant and the coordinate origin can be placed at any particular vortex, yielding the same power law for the radial distribution of vortices from the origin. In a finite system the scale invariance can only persist up to scales of order the largest vortex separation. Furthermore, the vortices must be separated by a minimum distance to satisfy the assumptions used in deriving the energy spectrum from the point-vortex model. In practice it is necessary to use a self-consistent sampling scheme in order to generate the correct power-law distributions for localized finite configurations.
Our sampling scheme for a configuration of N vortices is: 1. Sample the radial distance r p of each vortex from the coordinate origin according to a power law probability distribution ∝ r −ᾱ p . The exponentᾱ is distinct from α due to the finite system size and localization of the distribution.
3. For a given N, compute the kinetic energy spectrum, averaging over n s = 100 samples of vortex position data found via the foregoing routine. Iterate until the spectral power law of interest is found.
We will sample a power-law distribution of vortex distances from the origin. In practice there is a lower (r min ∼ ξ) and upper (r min ∼ R = system size) cutoff for power-law scaling. We thus wish to sample the distribution that is normalized on the interval r min ≤ r < r max . We sample r values from this distribution using uniform random variates x ∈ [0, 1) via the transformation [56]: Note that whenᾱ < 1, as is always the case in this work, some kind of ultraviolet cutoff is required for the distribution to be normalizable. Here we have made a choice that gives a power-law distribution over a well-defined scale range (See e.g. Ref. [56] for other common choices).

B. Classical velocity distribution of a large cluster
A fully polarized configuration of vortices with a given radial power law forms a quantum analogue of the coherent vortices of forced 2D turbulence in classical fluids. As an arrangement of many vortices, the velocity distribution must approach a classical limit, according to Bohr's correspondence principle, in much the same way that a rotating Abrikosov lattice generates a velocity field that approaches that of a rotating rigid body [50]. In what follows we find the classical velocity field, and identify the physical significance of the radial exponentᾱ.
To determine the velocity field we first compute the fraction of vortices enclosed by a circular contour around the origin, with radius r, for the distribution (82). Taking r min = ξ, r max = R, and considering scales ξ ≪ r ≪ R we have where we used the fact that where n = f r N is the number of vortices enclosed by the contour of radius r. Using (84), we obtain the velocity profile This inertial cluster has a power-law velocity profile determined by the specific radial exponentᾱ. In contrast, the velocity profile of an Abrikosov vortex lattice rotating at frequency ω approaches that of a rigid body v φ (r) = ωr, and a single quantum vortex has profile v φ (r) = /mr. The inertial cluster velocity profile (86) is compared with sampled distributions (see below) in Figure 4.

C. Scale expansion of a large cluster
We now illustrate the role of the radial power law exponentᾱ and the classical velocity distribution by sampling a large cluster. We vary our choice of scale range for the vortices (r max in Eq. 82), and investigate how the range of k −5/3 changes.
A characteristic wave number measuring cluster size for a sample involving n c vortices in a given cluster is given bȳ Wherer = 1 n c n c p=1 r p . In the figure we plot the scalek to give an indication of the range of k −5/3 scaling. We also plot k max corresponding to the largest vortex separation scale in the system For k < k max the velocity field approaches that of a single vortex.
In Figure 4 vortex distributions are sampled for N = 100 vortices of the same sign using the sampling scheme (81), (83), for r min = ξ and different values of r max . Individual samples are shown to indicate the spread of vortices, and the velocity profiles and kinetic energy spectra are computed by averaging over n s = 100 samples. The mean azimuthal velocity compares well with Equation (86), showing that the specific radial exponentᾱ in (82) also determines the power law of the azimuthal velocity field, as seen in (86). For scales larger !!""" " !""" !!""" " !""" Eq. (61) !!""" !#""" " #""" !""" " than r max , v φ (r) returns to the r −1 scaling for a charge N vortex [clearly seen in Fig. 4 (a)]. N = 100 vortices distributed up to r max = 60πξ gives an inertial range in the kinetic energy spectrum of approximately one decade, forᾱ = 0.4 (as the number of vortices in a given scale range increases,ᾱ → 1/3). Increasing the upper scale cutoff to r max = 400πξ, 4000πξ gives α = 0.65, 0.83, with 1.5 and 2 decades of inertial range respectively. Thus expanding the scale range of power-law behavior expands the inertial range, but requires the azimuthal velocity profile to steepen. In the UV region of the spectrum the k −3 law always holds [Eq. (71)], while the inertial range holds fork k ξ −1 , and single vortex behavior is apparent for k < k max .

D. Clustering in a neutral distribution
We now consider the role of increasing clustering in expanding the inertial range of the kinetic energy spectrum.
First, in n s = 100 samples, we distribute N + = N − = 100 vortices over a uniform periodic domain; one such sample is shown in the inset of Figure 5 (a). The corresponding kinetic energy spectrum shows the correct UV-region spectrum given by Eq. (71), and also approaches the E i (k) ∝ k form for k < k max , given by Eq. (70). For k max k ξ −1 the spectrum is less steep than k −5/3 and the system lacks an inertial range. In Figure 5 (b) the vortices are sampled as 40 clusters of n c = 5 vortices of the same sign according to (81) and (83), withᾱ = 0.8, and r min = ξ, r max = 10πξ. The 10 +ve and 10 -ve cluster centers are uniformly distributed as in Fig. 5 (a), as seen in the sample (inset). This distribution yields a k −5/3 power-law kinetic energy spectrum over ∼ 1 decade of wave numbers. By further expanding the scale of clustering while reducing the number of clusters to preserve N + and N − , we find the inertial range can be extended. In Figure 5 (c) samples consists of n c = 20 vortices in each cluster, distributed between r min = ξ −1 and r max = 100πξ, withᾱ = 0.75, and giving ∼ 2 decades of inertial range. We note that, compared with Fig. 4,k, shown in Fig. 5 (b) and (c) does not correspond so well with the lower bound on the inertial range, presumably because of the significant space between clusters in the neutral system.

V. FORCED 2D QUANTUM TURBULENCE
A. Scenario of forced turbulence 6. Illustration of an inertial range (non-shaded region) for the incompressible portion of kinetic energy in forced compressible 2DQT. The E(k) ∝ k −3 region arises from the structure of the vortex core and thus is not a signature of vortex configurations and vortex turbulence. This ultraviolet region can thus not support energy cascades, nor does this region correspond to enstrophy cascades. Net energy injected at k F ∼ ξ −1 in the form of vortices can only move towards the infrared. The Kolmogorov law E(k) ∝ k −5/3 occurs in the inertial range of fully developed turbulence. The far-infrared region is given by E(k) ∝ k for a system with no net vorticity, and is evident for k ≪ L −1 where L is the largest intervortex distance. For forcing at smaller wavenumbers, the spectrum may be more complex, possibly involving other forms of energy and enstrophy flux.

in a compressible 2D superfluid
The canonical model of 2D classical turbulence consists of a velocity field described by the 2D Navier-Stokes equation The density ρ of the incompressible fluid is held constant by the pressure field p, ν is the kinematic viscosity, f v is a forcing term and λ represents linear frictional damping arising from irreducible 3D aspects of the system in which the 2D flow resides. If the fluid is subjected to suitable forcing it will develop an inverse energy cascade and a direct enstrophy cascade, with associated k −5/3 and k −3 power laws respectively [13]. An inverse energy cascade induced by small-scale forcing can be steady because the −λv term damps energy at large length scales [9]. For a homogeneous compressible superfluid subject to forcing from an external potential, (12) can be written as where the forcing f v ≡ −∇V(r, t)/m is assumed to be spatially localized. The lack of a −λv frictional damping term means, in the classical case, that if an inverse energy cascade develops as a result of steady forcing, it is not expected to be stationary. In the superfluid case, the compressibility of the fluid allows vortex-antivortex annihilation, which couples energy into the sound field. This interaction between the sound and vorticity fields renders the calculation of energy fluxes in compressible superfluid systems particularly difficult and somewhat ambiguous [30]. As we have shown above, in contrast with the classical Kraichnan scenario utilizing a 2D Navier-Stokes analysis, a k −3 spectrum for k > ξ −1 for a 2D quantum fluid is not caused by a direct enstrophy cascade but is rather a consequence of the vortex core structure, and thus should not be interpreted in terms of vortex configuration dynamics (note that vortex core shape excitations can be neglected, since they constitute a component of the sound field). Given forcing at a wavenumber k F ∼ ξ −1 , and minimal vortex-antivortex annihilation, the incompressible kinetic energy can only move toward the infrared. This scenario is shown schematically in Figure 6.
It has been shown that dipole recombination provides a route for a direct energy cascade to develop in 2D GPE dynamics [30]. This mechanism provides a means for oppositesign vortices to approach zero distance, coupling vortex energy to the sound field during vortex annihilation. However, if the forcing leads to significant clustering of like-sign vortices faster than recombination occurs, or prior to recombination occurring, dipole-decay will be strongly inhibited. This suggests that under the right conditions of forcing an inverse energy cascade can become the dominant mechanism of energy transport between distinct length scales.

B. Kolmogorov constant and clustered fraction
By making use of the universal onstrophy and the condition of continuity at kξ ≈ 1 we have found that the k −5/3 power-law given by (73) describes the spectrum of numerically sampled vortex configurations that exhibit a s −1/3 power law for the vortex separation data. While individual spectra and configurations do not give information about dynamics, in particular, the direction in k-space of any energy cascades, the power law suggests the existence of an inertial range comprised of vortices. In a cascade, such a configuration will transfer incompressible energy between scales while conserving energy. Assuming the infrared portion of our double-power-law analysis, namely (73), will also describe such a cascade, we can cast it as a statement about the Kolmogorov constant in terms of the one-vortex onstrophy and the slope of the radial wavefunction at the vortex core.
To write (73) in standard form, we introduce the unique Nvortex quantity with dimensions energy/mass/time that can be constructed from Ωξ 2 , m, and /Ωξ 2 : We then find where the remaining quantities have been absorbed into the dimensionless Kolmogorov constant: µ = 2 /mξ 2 in the homogeneous system, and the bar notation distinguishes the quantum system. In the dilute Bose gas the 2D interaction parameter is µ/n 0 = g 2 = 4π 2 a/ml, where l is the characteristic thickness of the three-dimensional system [30,57]. In terms of this length we find We emphasize that the physical input needed to arrive at this form of the Kolmogorov constant is (i) accounting for the structure of a compressible quantum vortex in determining the ultraviolet spectrum, and (ii) imposing continuity of the ultraviolet spectrum at kξ = 1 to a Kolmogorov powerlaw in the infrared. In classical turbulence, C 2D ≃ 7 [58]. To give an example of howC 2D may be evaluated for a compressible superfluid exhibiting 2DQT, we consider a 87 Rb BEC that is homogeneous in the x − y plane, and that is harmonically trapped in the z-dimension with trap frequency ω z = 2π×5000 Hz. We use m = 1.44 × 10 −25 kg, a = 5.8nm, for which l = 2π /mω z = 0.38 µm, and g 2 = 0.197 2 /m. For these values,C 2D = 0.212. Note also that by defining the configurational rate constant (91), we have confined this discussion to a scale invariant distribution involving N vortices. This expression suggests that such a configuration can support an inverse energy cascade at the rate ǫ N .
The foregoing discussion involves an ideal distribution of N vortices configured with the α = 1/3 power law. It is clear from Fig. 5 that the vortices do not have to all have the same sign of circulation, but they must be configured into clusters of vortices with the same sign. The universality of C 2D in classical turbulence in incompressible fluids leads us to postulate that the condition that all N vortices are power-law clustered can be further relaxed. For fully developed quantum turbulence involving N vortices, we interpret N as the participation number, representing the number of vortices in a scale-free turbulent configuration, which in the case of a fully polarized cluster is maximal. Imperfect clustering involves fewer vortices in power-law cluster configurations, and an effective participation number that is the number of clustered vortices N c < N, namely, the number with nearest-neighbors of the same sign. Making the replacement N → N c in (73), we propose the ansatz spectrum as a more general definition for systems that have incomplete clustering in the inertial scale range. We test this hypothesis in the next section in dynamical simulations of the forced dGPE.
It is important to note that the condition of continuity at kξ = 1 is no longer exactly met, since only the value N c ≡ N will produce an infrared spectrum with k −5/3 that is continuous with the ultraviolet power law approximation at k = 1/ξ. We also note that a more general measure of clustering, namely the polarization index (P) was introduced in Ref. [59] measuring the degree and type of spatial clustering of like-sign vortices in 3DCT. The Kolmogorov k −5/3 spectrum was found to correspond to the partially polarized value P = 1/3, while other scaling laws yield differing polarizations. The clustered fraction used in the present work is a simpler (global) measure of polarization as it does not contain information about the spatial distribution of vortices, but it is only relevant for 2DQT.
We can also write down an expression forC 2D for a general energy spectrum that may be computed numerically from simulation data, by making use of the ansatz (95). This is equivalent to using (92) with N → N c , from which we can define the function In a region where the spectrum is approximately k −5/3 ,C * 2D (k) will be approximately constant, and it may be compared with the prediction (93). We test this numerically in the next section.
We note that in 3DQT an energy bottleneck has been predicted for the direct energy cascade [59,60], and also observed in GPE simulations [61]. It occurs due a mismatch between the rates of energy transport at large length scales (hydrodynamic regime) and small length scales (Kelvin wave cascade). The mismatch causes energy to pile up at the length scale where the two cascades meet. This raises the possibility of a bottleneck in 2D, although Kelvin waves are disallowed in 2DQT, so this particular mechanism would not be relevant. However, for a given forcing mechanism, it is possible that the rate of transporting energy to large length scales may not be high enough to remove all of the vortex energy introduced at the forcing scale. Thus a bottleneck could still occur, and our assumption of continuity of the spectrum at kξ ≃ 1 may not hold in general. We return to this question in the next subsection, where we find some indication of an energy bottleneck at the forcing scale in numerical simulations.

C. Damped Gross-Pitaevskii Dynamics
We now consider a simulation of the forced dGPE that generates significant clustering of vortices of the same sign. The system consists of a homogeneous superfluid with periodic boundary conditions, stirred by dragging four Gaussian obstacle beams through it at a constant speed [32], thus modeling grid turbulence in a BEC. When an obstacle is dragged through a superfluid sufficiently rapidly, superfluidity cannot be maintained. For slowly moving obstacles, the superfluid will adapt to the forcing, and vortices are not formed. Above a critical velocity v c [62] vortex dipoles are periodically formed in the wake of the obstacle, injecting linear momentum into the superfluid. Sufficiently rapid motion (v ≫ v c ) causes many vortices to be nucleated behind the obstacles in a chaotic fashion [63], involving clustering of like sign vortices. Our choice of obstacle speed puts the system dynamics in the latter category.
We work in units of µ, ξ, and ξ/c for energy, length, and time, respectively. In these units the specific parameters we choose (see the previous subsection) are g 2 = 0.197µξ 2 , corresponding to homogeneous density n 0 = µ/g 2 = 5.26ξ −2 , and N tot = 5.5 × 10 6 particles in a homogeneous 2D system of side length L = 1024ξ. The Gaussian potentials each have fixed 1/e 2 width of w 0 = √ 8ξ, and height V 0 = 100µ, and are initially located at x = −L/2 + 8ξ, y = ±L/8, ±3L/8. Numerically, we proceed by first finding a ground state of (1), for a homogeneous system with periodic boundary conditions, subject to the localized obstacle beams. We then transform into a frame translating at v 0 = 0.8c, and maintain the obstacle locations relative to this frame, creating a dragging grid of stirring beams. A small amount of initial noise is added to the wave function to break the reflection symmetries of the system. We thus evolve the system according to (1) for the same potential, but with the Galilean transformed nonlinear operator L → L + i v 0 ∂ x . During evolution the dimensionless damping rate is set to γ = 0.003.
The time evolution of the system is shown in Figure 7 at three times, at approximately (1/3, 2/3, 1)L/v 0 , so as to avoid any periodic flow effects in the x-direction. The four obstacle beams generate many vortices (up to N ∼ 10 3 ), and significant clustering (N c /N 0.6). The incompressible kinetic energy spectrum shows a wide region that is well described by the k −5/3 form of Equation (95). For later times [ Figure 7 (b), (c)] the spectrum shows a significant pile up around the forcing scale k F ∼ ξ −1 , suggesting a mismatch between the rates of injection and transport of incompressible kinetic energy.
In Figure 8 we compare the functionC * 2D (k) [Eq. (96)], as numerically computed from our simulation data, with the analytical prediction of the Kolmogorov constantC 2D [Eq. (94)]. The region of k −5/3 appears as a broad flat region that is in close agreement withC 2D = 0.212 pertaining to our simulation parameters.

VI. CONCLUSIONS
To summarize, we have investigated relationships between the concepts of 2D turbulence in classical fluids and the emerging topic of 2D quantum turbulence of vortices, specifically as it relates to Bose-Einstein condensates. We established a link between the hydrodynamic limit of the damped GPE and the Navier-Stokes equations, providing an estimate of a quantum Reynolds number for superfluid flows in BECs. We have given a theoretical treatment of the incompressible kinetic energy spectrum that explicitly incorporates the vortex core structure in a compressible superfluid. The incompressible kinetic energy spectrum for a compressible superfluid is deconstructed in terms of single-vortex contributions determining a unique ultraviolet power-law where the energy spectrum scales as k −3 , and a contribution that depends on the configuration of vortices within the fluid that determines the infrared region of the spectrum. For the configurational regime we find: 1. The spectrum only depends on the distribution of vortex separations and the sign of the circulation of each quantum vortex. If the distribution of vortex separation s for a system of vortices of the same sign is a power-law ∝ s −α with exponent α = 1/3, the kinetic energy spectrum will take the universal Kolmogorov form ∝ k −5/3 , as shown for point vortices [53]. Localized clusters of N vortices of the same circulation with this power law distribution can be constructed by sampling using a specific radial exponentᾱ that depends on the number of vortices and the scale range over which they are distributed.
2. The azimuthal velocity field of a large cluster is determined byᾱ. By inflating the scale range of a cluster we find thatᾱ increases, the velocity field is steepened and the inertial range expands to larger scales. In a neutral system the inertial range can be extended by increasing the size of clusters while decreasing their number.
3. The universal form of the UV region of the kinetic energy spectrum imposes a strong constraint. If the Kolmogorov power law occurs in the infrared region, then the postulate of continuity between the infrared and ultraviolet regions completely determines the spectrum when the ultraviolet and infrared regions are approximated as power laws. Physically, this corresponds to the inertial range extending down to the smallest configurational scale of the system ∼ ξ. We note that the postulate of continuity may not be relevant for all systems or forcing mechanisms.
4. We infer an analytical value for the Kolmogorov constant [Eq. (93)] under the conditions of spectral continuity at the cross-over scale for a system of vortices of the same sign. To assess the validity of this inference for dynamical situations we compare our analytical results with spectra from a numerical simulation of the forced dGPE for the specific case of a dragging a grid of obstacles through an otherwise homogeneous BEC. We find reasonable agreement provided we introduce the concept of a clustered fraction N c /N ≤ 1, which is the fraction of vortices that have same-sign nearest neighbors. This measure discounts all vortex dipoles from the configurational analysis. We then observe good agreement between our Kolmogorov ansatz [Eq. (95)], and the spectrum calculated from the dGPE data. We also find that the predicted value of the Kolmogorov constant is in close agreement with the numerical simulations [ Figure 8].
We note that while our analysis indicates that vortex positions and circulations are enough to determine an approximate incompressible kinetic energy spectrum, the reverse is not necessarily true: a Kolmogorov spectrum does not carry information about specific vortex distributions. Nevertheless, our analysis does indicate that the number of vortices in a quantum fluid can in principle be directly determined from the ultraviolet energy spectrum. Moreover, the concept of a cascade in turbulence implies system dynamics and energy transport, yet aside from our numerical simulation example, our analytical approach is an instantaneous measure. Importantly, one must determine means of characterizing vortex motion and relate such dynamics to the cascade concept. The field of 2D quantum vortex turbulence is relatively new, compared with the much longer histories of 3D superfluid turbulence, 2D classical turbulence, and even dilute-gas Bose-Einstein condensation. Point-vortex models have been extensively used in descriptions of superfluid dynamics as well as in 2D classical turbulence, although point-vortex distributions can only serve as approximate models of real 2D classical flows. Our approach merges concepts from each of the above subjects in order to develop a new understanding of 2D quantum turbulence. By considering the compressibility of a dilute-gas BEC, we find an analytical expression for the ultraviolet incompressible kinetic energy spectrum and an Nvortex equivalent of enstrophy in a quantum fluid, which we term the onstrophy. For the infrared region, point vortex models are sufficient, and vortex configurations serve to identify spectra as summarized above. Taken together, the primary new outcome of our analysis is a link between vortex distributions, vortex core structure, and power-law spectra for 2D compressible quantum fluids.
Future work on 2D quantum vortex turbulence will involve numerical simulations and comparisons with our analytical results, extension of this analysis to confined and inhomogeneous density distributions, characterization of vortex dynamics and the time-dependence of energy spectra particularly in relation to clustering [64], inverse-energy cascades, and the nonthermal fixed point [31,65], and investigation of connections with weak-wave turbulence in BEC [16,18,19,[66][67][68][69]. We also believe that observing vortex distributions such as the α = 1/3 power-law for localized clusters may provide a new means of quantitatively characterizing 2D quantum vortex turbulence through direct experimental observations of vortex locations in a forced 2D superfluid, and we are working towards realizing such experimental observations.