Angular Momentum and the Absence of Vortices in the Cores of Fuzzy Dark Matter Haloes

Scalar Field Dark Matter (SFDM), comprised of ultralight (& 10 eV) bosons, is distinguished from massive (& GeV), collisionless Cold Dark Matter (CDM) by its novel structure-formation dynamics as Bose-Einstein condensate (BEC) and quantum superfluid with wave-like properties, described by the Gross-Pitaevski and Poisson (GPP) equations. In the free-field (“fuzzy”) limit of SFDM (FDM), structure is inhibited below the de Broglie wavelength λdeB, but resembles CDM on larger scales. Virialized haloes have “solitonic” cores of radius∼ λdeB that follow the ground-state attractor solution of GPP, surrounded by CDM-like envelopes. As a superfluid, SFDM is irrotational but can be unstable to vortex formation; outside of vortices it remains vorticity-free. We previously showed that halo cores can form vortices, from angular momentum expected during structure formation, if a strong enough repulsive self-interaction (SI) is present, which inhibits structure below a second length scale λSI, with λSI > λdeB, suggesting FDM cores could not. FDM simulations later found vortices, but only outside halo cores, consistent with our suggestion. Extending our analysis now to FDM, we show explicitly that vortices should not arise in solitonic cores from angular momentum, modelling them as either Gaussian spheres or compressible, (n = 2)-polytropic, irrotational Riemann-S ellipsoids. For typical halo spin parameters, angular momentum per particle is below ~, the minimum required for one singly-quantized vortex in the center. Even for larger angular momentum, vortex formation is not energetically favoured.


INTRODUCTION
In recent years, ultralight bosonic dark matter has attracted increasing attention in the community as an alternative to the standard cold dark mater (CDM) paradigm. While CDM is modelled astrophysically as a collisionless "gas" without internal particle self-interactions, except for gravity, ultralight bosons can exhibit a richer particle phenomenology, affecting structure formation and galactic dynamics. This way, the search for signature effects on galactic and cosmological scales helps to distinguish ultralight bosons from CDM. Moreover, the theoretical predictions of excessive clustering in standard CDM face issues, when confronted with observations. These problems have been known for a while as "cusp-core problem", "missingsatellite-problem", "too-big-too-fail-problem" and possibly further issues, see e.g. Weinberg et al. (2015). On the other hand, ultralight bosons with m ∼ (10 −25 − 10 −18 ) eV, bring about a substantial cutoff in the small-scale structure formation, potentially alleviating the above problems, since structure formation is inhibited below a length scale which depends upon the boson parameters, either the de Broglie wave length λdeB discussed below, or a second length scale λSI > λdeB, if a strongly repulsive particle self-interaction (SI) is present. Ultralight bosons encompass a large family of models. The common thing they share, however, is that they are described by a ⋆ E-mail: sonja.schobesberger@protonmail.ch † E-mail: tanja.rindler-daller@univie.ac.at ‡ E-mail: shapiro@astro.as.utexas.edu scalar field, hence also called "scalar field dark matter" (SFDM) and the further requirement that they shall be "cold" (i.e. non-relativistic) as of some point in their cosmic evolution, in order to be able to make up for the present dark matter (DM) energy density. Notwithstanding the small-scale problems, bosons of higher mass are also considered and they remain viable DM candidates. The predecessor of them all is the QCD axion, a fundamental (pseudo-) scalar particle with a mass of about m ∼ 10 −5 eV, introduced to resolve the CP problem of QCD (Peccei & Quinn (1977)). It has been suggested as a DM particle in Weinberg (1978) and Wilczek (1978). Later, ultralight particles have been considered in string theories and other extra-dimensional models, see e.g. Frieman et al. (1995), Günther & Zhuk (1997), Svrcek & Witten (2006), Arvanitaki et al. (2010), or Fan (2016. In this work, we will consider the simplest model of SFDM with a Lagrangian that contains kinetic energy and a (rest-)mass term 1 . Furthermore, we will only consider the non-relativistic description of such DM particles, which is suitable for the physics of DM haloes. In this context, it has been also known as Bose-Einstein-Condensate (BEC) dark matter (or BEC-DM, or BEC-CDM) in order to emphasize the idea that the ultralight bosons undergo Bose-Einstein condensation, which motivates the use of a scalar field in the first place. SFDM without SI has been also called fuzzy DM, ψDM, wave DM or free SFDM. In this work, we will call it fuzzy DM (FDM), in observance of the early paper by Hu et al. (2000), wherein that term was coined.
However, the "fuzziness" of FDM in the sense introduced by these authors actually requires more than just the lack of SI. It basically demands that the de Broglie wavelength of the bosons, for given mass m and collective velocity v, becomes of the order of galactic scales, λdeB 2π = mv = 1.92 kpc 10 −22 eV m 10 km s −1 v , which implies that the corresponding velocity should be of order the virial velocity of the gravitationally bound object. In this model, the de Broglie length is of same order of magnitude as the characteristic Jeans length, or the size of the smallest gravitationally bound object expected in BEC-DM. BEC-DM as FDM can cure the small-scale-crisis, once the scale that fits the smallest (sub)halo will host the smallest observed types of galaxies. In other words, if we believe that the smallest galactic DM structures should not undershoot a size of about 1-2 kpc, then we require a boson mass of order m 10 −22 eV. As a result, this regime of very small masses has received considerable amount of interest in the literature.
In contrast, the QCD axion and generally other bosons of much higher mass, m ≫ 10 −22 eV, have de Broglie lengths λdeB ≪ 1 kpc, i.e. invisible on any galactic scales of interest, hence they would not be examples for FDM in the strong sense, for their wave nature is not potent on galactic scales anymore, implying they would fail to solve the small-scale problems of CDM. Yet, they remain valid models, as long as the nature of DM is unsettled and as long as no other reasons have been found to exclude them definitely.
The quantum nature of SFDM as a BEC and quantum superfluid distinguishes its dynamics from that of a collisionless gas like CDM, in more ways than just the suppression of structure below the de Broglie wavelength of the former. For example, BECs are generally irrotational, i.e. vorticity-free. However, as known from laboratory experiments, quantum vortices can arise, either from the rotation of a single BEC, or from the merging of multiple BECs which initially had no net angular momentum. Vortices result that represent topological defects, outside of which the rest of the system remains vorticityfree. Like "tornadoes", in which the density drops to zero towards the centre, while velocity diverges, such quantum vortices can be a significant departure from the smooth background, with dynamical and structural consequences. As a quantum superfluid, SFDM, or BEC-DM, must also be irrotational -free of vorticity. However, even with irrotational initial conditions, BEC-DM can also become unstable to vortex formation, as found for BECs in the lab.
In general, vorticity in quantum fluids is discrete and quantized. Also, vortices are fundamental building blocks of quantum turbulence, because quantum turbulence exemplifies itself in the form of large disordered vortex tangles, where individual vortices interact via vortex reconnections. Moreover, these turbulent flows may also exhibit partial spatial polarization by organizing into vortex bundles enabling large-scale energy flows. The field of quantum turbulence is far from being understood (see e.g. Tsatsos et al. (2016) for lab examples). This paper will focus on a well-posed problem, namely the question of whether vortices can form within rotating FDM ground-state structures in gravitational equilibrium. We will study the appearance of a single vortex, we do not study quantum turbulence.
Several authors have previously pursued related questions. Brook & Coles (2009) consider uniform objects, and estimate the respective forces to determine when a vortex within this background should arise and be stable. Silverman & Mallett (2002) postulated vortices in galactic haloes by comparing the critical angular velocity for vortex formation with the rotation rate of M31, although the formulae they used apply only to systems with strong SI, but without self-gravity. Subsequently, Yu & Morgan (2002) and Zinner (2011) have studied the influence of a vortex lattice on the velocity profile of a spherical galactic halo. In the end, all these papers are very heuristic in nature.
Our analysis here aims to study the question of vortex formation with analytical rigour, and as such, it is inspired by the work of Rindler-Daller & Shapiro (2012) (henceforth abbreviated RS12); earlier results were published in Rindler-Daller & Shapiro (2010). There, the problem was analyzed whether the typical amount of angular momentum, parameterized by the cosmological spin-parameter from large-scale structure formation, applied to BEC-DM haloes (and halo cores) was sufficient to make vortex formation energetically favoured, and for which parameters of the BEC-DM boson mass and SI coupling strength that would happen, if so. The analysis of RS12 was carried out in the regime of strongly repulsive SI, or the so-called Thomas-Fermi (TF) regime, where gravity is balanced by the SI pressure. The analysis applied to rotating systems in equilibrium (i.e. rotation is added to balance gravity), which are ground state solutions of the underlying equations of motion.
It was shown that vortex formation is favoured for a large part of the boson parameter space. In fact, the parameter space of favoured vortex formation overlaps the parameter space of the TF regime to a great extent. In addition, approximate halo models for this regime were introduced in RS12, which were shown to be viable models for rotating BEC-DM haloes, having different amount of angular momentum. In particular, rotating haloes (or halo cores) without vortices were shown to be well represented by so-called irrotational Riemann-S ellipsoidal figures, whose equation of state is an (n = 1)polytrope, as appropriate for the TF regime. Now, when RS12 found that vortices were likely to form in BEC-DM haloes with SI, their analysis showed that, even within the TF regime, a minimum SI coupling strength was required to make vortex formation energetically favoured. It was reasoned, therefore, that, if SI were absent altogether, as in FDM, vortices would not form at all, i.e. in halo cores supported against gravity only by FDM quantum pressure and the amount of rotation expected from largescale structure formation. The purpose of the present paper is to revisit this suggestion by RS12 that no vortices are expected for FDM, by performing a new, detailed analysis along the lines of that which RS12 performed in the TF regime 2 . As we will show by this new analysis, the expectation expressed in RS12 was correct and in accordance with 3D simulations of FDM halo formation which came later. In particular, while vorticity has been reported in FDM structure formation simulations in some locations, it is entirely absent from the halo cores. Our results help to explain this exclusivity.
Haloes of significant size will be not only composed of the core, but will have an envelope region, whose properties on average should resemble CDM, if structure formation of BEC-DM shall be similar to CDM on large scales. In fact, we suggested in Rindler-Daller & Shapiro (2014) that the polytropic, SI pressure support that set the size of halo cores in the TF regime would be supplemented on larger scales by the support of wave motion, generated by the wave nature of BEC-DM and its quantum pressure during the virialization of haloes that assemble from infall and mergers, making it possible for haloes to be much larger than their polytropic cores in which only SI dominates.
Meanwhile, FDM has been studied in more detail, including simulations of halo formation that report all virialized haloes have solitonic cores of the size of the de Broglie wavelength (as evaluated inside haloes), supported against gravity by quantum pressure, see Schive, Chiueh & Broadhurst (2014a), Schwabe, Niemeyer & Engels (2016) and Mocz et al. (2017). And all haloes also show a wave-supported envelope outside this solitonic core, with a profile that resembles that in CDM haloes, but in which wave motions provide the random internal motions responsible for virial equilibrium, instead of random particle orbits, as first described by Rindler-Daller & Shapiro (2014). These same simulations find that vorticity and signs of quantum turbulence are generated during structure formation (from vorticity-free initital conditions), but only outside of the solitonic cores. The origin of this vorticity has not been well-studied, but its absence from solitonic cores is consistent with our suggestion in RS12 that vortex formation by instability in the presence of angular momentum requires a sufficiently strong repulsive SI in order to be energetically favoured. However, the analysis in RS12 was limited to the TF regime and, so, our purpose here is to extend that analysis to the case of FDM.
The fuzzy regime presents an important challenge for analytical methods, in that all length scales of interest, in particular the characteristic size of the system under question and perturbations in that system like vortices, are of similar order of magnitude. This is in stark contrast to the TF regime studied in RS12, where length scales are separated by many orders of magnitude, allowing to focus on leadingorder considerations. We shall show that, unlike the TF regime, for solitonic cores in FDM, in which gravity is balanced by quantum pressure and rotation, vortex formation cannot be triggered by angular momentum as it is in the TF regime. For vortex formation to be triggered by angular momentum, the specific angular momentum must first satisfy a necessary condition that it exceeds the minimum value that gives each particle an angular momentum of . If this necessary condition is satisfied, then it is further required that vortex formation be energetically favoured, in order to establish that vortex formation will take place. In the TF regime, both conditions can be met for the typical amounts of specific angular momentum for cosmological haloes, for a large range of particle mass and SI strength. However, for FDM, we will show here that the necessary (minimum) condition is generally not met for typical amounts of halo angular momentum. We further show that, even for angular momentum large enough to meet the necessary condition, vortex formation is nevertheless not energetically favoured. This is consistent with and can explain the fact that simulations of structure formation in the FDM model do not find vortices in the solitonic cores of FDM haloes.
Our current analysis will share certain assumptions of RS12: we will also confine to a DM-only analysis; we will limit the consideration to rotating, equilibrium systems, i.e. regions whose size is comparable to the characteristic length scale of hydrostatic objects, and whose density profile corresponds to the (approximate) ground state of the underlying equations of motion, the Gross-Pitaevskii-Poisson system of equations. In FDM, considered in this paper, the characteristic scale is the de Broglie length of bosons, which are distributed according to a numerically calculated solitonic profile, while in the TF regime the characteristic scale is proportional to the radius of an (n = 1)-polytrope of the underlying density profile. The (approximate) ground-state solutions considered in RS12 and here are viable models for the central core region of a big halo, or of the entire region of a small halo which consists only of the core. The reader shall keep this in mind, when we talk about "haloes" and "halo cores" in our work. However, the (approximate) ground-state solutions that we study here also apply to systems composed of DM bosons of higher mass, not just ultra-light. This paper is organized as follows. In section 2, we present the fundamental equations. In section 3, we review the basic concepts of SFDM in the fuzzy limit (FDM) and present models of non-rotating solitonic cores. The latter are the building blocks of our models with rotation, which are introduced in section 4. section 5 includes the derivation of our main results, namely the analysis of the conditions of vortex formation, where we show that vortices due to angular momentum are not expected in FDM cores. Finally, our conclusions and discussions are presented in section 6. Also, some more technical background is provided in two appendices.

FUNDAMENTALS OF GRAVITATIONALLY BOUND BEC-DM STRUCTURES
In the non-relativistic limit, BEC-DM objects under self-gravity are generally described by the Gross-Pitaevskii (GP) equation of motion (Pitaevskii (1961), Gross (1961), Ruffini & Bonazzola (1969)) for the complex scalar wavefunction ψ( r, t) of the bosons in the Bose-Einstein condensate (BEC), which is coupled to the Poisson equation (GPP), where |ψ| 2 ( r, t) = n( r, t) corresponds to the number probability density, m denotes the boson mass and Φ the gravitational potential. We assume that all N particles comprising a given object of volume V are in the condensed state described by ψ, hence The last term in Eq. (3) describes an effective 2-particle SI with a coupling 3 constant g = 4π 2 as/m. Its sign is determined by the sign of the scattering length as. Within this framework, m and g are the fundamental particle parameters. In this work, we will set g = 0, but in order to compare our analysis with previous results we include it in the description of the basic equations. Equ.
(3) admits a quantum fluid description upon introducing hydrodynamical variables via the Madelung transformation (Madelung (1927)): where ρ( r, t) = m|ψ| 2 is the mass density. Through the current and the identification of we can define the bulk (or flow) velocity field as the phase gradient v = m ∇S.
Assuming particle number is conserved, we identify an Euler-like equation of motion, and a continuity equation, The so-called Bohm quantum potential, defined as gives rise to what is often referred to as 'quantum pressure'. In addition, SI gives rise to a pressure of polytropic form, with n = 1 and the polytropic constant Kp depends upon the DM particle parameters m and g.
In the spirit of RS12, we will also restrict our analysis to stationary systems and their corresponding energy, i.e. to the time-independent GP equation including the chemical potential µ, ∆ψs( r) + (mΦ + g|ψs| 2 )ψs( r) .
The time-independent GP equation can be obtained from the timedependent equation (3) by inserting the state Stationary states have this form of wavefunction, which evolves harmonically in time and yields the time-independent density ρ = m|ψs| 2 and hence a time-independent gravitational potential. Of course, ψs itself can be decomposed analogously to ansatz (6) as where from now on the subscript s will be omitted. The GP energy functional is given by By means of decomposition (17) we can write the total energy, as a sum of the total kinetic energy the gravitational potential energy, and the internal energy, KQ is the part of the kinetic energy accounting for the quantum-like phenomena; it is absent in classical galactic dynamics. T describes the kinetic energy of the system due to bulk motions, like rotation or internal motions. USI = PSIdV is determined by the SI pressure given in (14). Since we are interested in g = 0, we will have that USI = 0. However, in the course of our analysis we will see that, under certain assumptions, KQ can be approximated by an energy which formally looks like the internal energy of a polytrope: the internal energy which arises from any polytropic pressure, P = Kpρ 1+1/n , is given as Here, we still leave the general dependency on n. Using all these energy contributions, we can write the equilibrium scalar virial theorem as We will assume that gravitationally bound FDM halo cores in equilibrium fulfil (approximate) virial equilibrium.

Vorticity
The GP framework was originally designed to account specifically for quantized vorticity in BECs. Indeed, the Madelung transformation unveils the superfluid character of this dissipation-free mean-field formulation. A zero-viscosity (i.e. dissipation-free) fluid has a conservative velocity field, i.e. a gradient flow v ∝ ∇S, which implies irrotationality, i.e. ω = ∇ × ∇S = 0 . However, it is only at first glance, that the formulation leaves no room for fluid vorticity. Wherever the mass density ρ = m|ψ| 2 vanishes, the phase S and hence the velocity flow v become discontinuous and ill-defined. The implication of irrotationality would only hold true if the phase function S had continuous first and second derivatives everywhere. However, this is not any longer the case, once a vortex line appears, where v diverges along its centre. (Since the density goes to zero towards the centre, there are no particles moving with infinite speed inside the vortex core.) It turns out that the phase function along vortex lines has a non-trivial winding and, as a result, the vorticity there does not vanish: Since the wavefunction is required to be singly-valued, a circulation along a contour C enclosing a vortex may not change ψ. Hence, S can vary at most by 2πd, where d is the winding number, also called the vortex charge. As a consequence, this circulation is an integer multiple of the quantum of circulation κ = h/m, The wavefunction of an axially-symmetric vortex aligned with the z-direction at the origin in cylindrical coordinates is a stationary solution and has the form The velocity field around this vortex has a form given by the directed along the azimuthal direction. It turns out that the above vortex wavefunction is an angular momentum eigenstate, where the z-component of the angular momentum is given by lz = d and thus yields the total angular momentum LQM denotes the angular momentum required to sustain a singlyquantized (or singly-charged) vortex within an object composed of N bosons. There are many perspectives on such a vortex. It can be seen as a perturbation, topological defect, or excitation with a potentially higher energy than the quantum ground state. For simplicity, we will focus in this paper on singly-quantized vortices with d = 1, which tend to be energetically more favoured than vortices with d > 1.

THE FUZZY REGIME
Considering the three terms on the right-hand side of the Euler-like equation of motion (11), two regimes can be distinguished. On the one hand, there is the regime where SI, or in other words the scattering of the DM particles is entirely neglected and hence solely quantum pressure works against gravity and prevents gravitational collapse. This regime is referred to as fuzzy regime, in conjunction with the requirement of Eq.
(2). In RS12, this regime was called BEC-CDM of type I. On the other hand, there is the regime of strongly interacting particles, where SI pressure dominates and balances gravity. This regime is the TF regime (or type II BEC-CDM in RS12), and its characteristic length scale is usually also required to be of order 1 kpc, in order to resolve the CDM small-scale crisis. The regime of interest in this work is the fuzzy regime, so it is important to review a few of its properties that are necessary in order to justify our choice of approximate models in the forthcoming.
The fuzzy regime effectively amounts to the limit of no SI (g = 0). Thus, we are interested in solutions of the time-independent GPP system in (15) and (4) without SI, namely This system of equations has been also known in the literature as the Schrödinger -Poisson (SP) or Schrödinger -Newton equations, see e.g. Moroz, Penrose & Tod (1998), Tod & Moroz (1999), provided we identify the eigenenergy E with the chemical potential µ. Without loss of generality, Tod & Moroz (1999) assume that ψ is real and they introduce the non-dimensionalized functions 4 Σ and V by setting 5 This change of variables, together with the assumption of spherical symmetry, simplifies the equations to a new set, The prime denotes differentiation with respect to the spherical radial coordinate r. SP admits a specific scaling invariance: given a solution (Σ(r), V (r)), for any arbitrary real γ, there is another solution of the form 4 Tod & Moroz (1999) use the notation S for Σ, which we avoid, because S in our paper refers to the phase function. 5 Notice that their function U has the dimension of gravitational energy and corresponds to mΦ.
Their analysis yields that there exists a discrete family of finite, smooth, normalizable solutions where the jth solution (j ∈ N) has j − 1 zeros and that the energy eigenvalues increase monotonically towards 0 with increasing j. In the case of spherical symmetry, two relevant asymptotic forms are derived: The solutions admit power series expansions near r = 0, Σ0V0r 2 +O(r 4 ) and V = V0 − 1 6 Σ 2 0 r 2 +O(r 4 ). (37) On the other hand, for large r and real constants A,B and k, there are solutions looking like In later sections of their paper, they show that the bound state solutions -those are the ones we are interested in -, have exactly these asymptotic forms.
In principle, the above differential equations (31)-(32) or (34)-(35), respectively, can be solved by looking for regular and finite solutions for the variables Σ (resp. ψ) and V (resp. Φ) of the SP system, reducing to an eigenvalue problem which can be solved numerically. Many authors, such as Kaup (1968), Ruffini & Bonazzola (1969), Membrado, Pacheco & Sañudo (1989), Seidel & Suen (1994), Guzmán & Ureña-López (2004), or Hui et al. (2017, have numerically calculated eigenstates of the SP system in spherical symmetry, e.g. assuming isolated systems, i.e. that ψ and Φ approach 0 as r goes to infinity and that they are regular at the origin. The scaling invariance means that a solution level j forms a one-parameter family which can be specified by the total mass M . Then, the central density is given by where the dimensionless constant ρj (depending on the eigenstate label j) is calculated numerically. The densest state is the ground state, j = 1, from where the central density strongly decreases with increasing level number j. The j = 1 "soliton" state is a longterm attractor for bound, isolated FDM objects, since excited states decay to the soliton state, through dispersion of probability density to infinity, a process known as "gravitational cooling" (Seidel & Suen (1994)). The total mass M = N m of the soliton is conserved and finite, however it has no compact support and thus has to be cut off artificially, if we want to have a finite size. It is customary to pick a radius which includes 99% of the mass. The resulting massradius relation of the bound ground-state has been calculated in Membrado, Pacheco & Sañudo (1989), where "S" stands for the numerically calculated "soliton". This result leads up to a general discussion of characteristic length scales, before we present our models. The analysis of vortex formation in RS12 focused on the TF regime. While conditions were established which determine when a halo (core) is described in either one of the regimes, the focus was targeted on haloes in the TF regime, i.e., the quantum-pressure term, eq.(13), was neglected at the scales of the halo (core), but not at the scale of vortices within such haloes! There is a characteristic length scale introduced in the classic literature on BECs, for example Pethick & Smith (2008), namely the so-called healing length 6 ℓ. This length describes the distance over which the wavefunction tends to its background value when subjected to a localized perturbation and is the result of balancing the quantum kinetic term with SI, withρ = 3M/(4πR 3 ), M and R are the total mass and radius of the object and ψ has dimension [Length −3/2 ]. To put it another way, it takes the BEC this distance to "heal" a local disturbance. By estimating the forces associated with the quantum-kinetic term and the SI term in the same way, RS12 conclude that in the TF regime the characteristic size of the system -the (halo) core size -is much larger than both, the healing length and the de Broglie wavelength, i.e.
Moreover, it is shown that λdeB ∼ 4.3 ℓ. However, in the fuzzy regime we have no SI. By replacing the SI energy by the gravitational energy, we introduce the gravitational healing length, ℓgrav, and derive it by setting the quantum kinetic energy equal to the gravitational energy, see the GP energy functional (18), The kinetic energy is of order 2 /(2mℓ 2 grav ) and the gravitational energy is of order mGM/(2ℓgrav), which yields where, again, M and R are the total mass and radius of the object. Approximating the spatial dimension in Φ with ℓgrav, i.e. with the same length scale used to approximate the Laplace operator on the left-hand side, corresponds to the notion of a Jeans analysis, where ℓgrav can be understood as the smallest length scale for bound structures, and also as the length scale of local perturbations. (Another possibility would have been to choose the gravitational energy to be of order mGM/(2R), which represents global properties as opposed to the local sensitivity of the Laplace operator.). The forthcoming energy analysis in the fuzzy regime gives rise to a quantity with the dimension of mass, which we define to be the characteristic particle mass, mc, whose dimensional combination is As a result, we have Previous works used a slightly different definition of the characteristic particle mass, e.g. mH in RS12. The comparison yields In BEC-DM without SI (i.e. FDM), any gravitationally bound, 6 The literature, incl. RS12, use the notation ξ for the healing length, which we avoid, because later we will use the notation ξ in a different context.
ground-state solution of equ.(31)-(32) has a spatial extent of order λdeB, and this is also true for the solitonic cores found in simulations of FDM. If we require where we use vcore ≈ vcirc = (GM/R) 1/2 , we have On the other hand, using R99,S of equ.(40) in the definition (46) yields m mc = √ 9.946 ≈ 3.154.
However, we can explicitely see again, what we have pointed out before, namely that (40) and (45) are within a factor 10 of similar order of magnitude, in stark contrast to (43) of the TF regime. Therefore, it follows from (47) that m has to be of the same order as mc, or in other words, there is only a relatively small range allowed for m/mc in the fuzzy regime. In order to factor in the imprecision in all these estimates, and in light of our models discussed below, we will consider a fiducial choice of range for that ratio in the forthcoming analysis, namely There are two interpretations for our characteristic particle mass mc.
On the one hand, it is the mass which results from requiring (It is precisely the fact that all length scales in this regime are of similar order, which makes the fuzzy limit of the GPP framework a challenging ground for analytic theory.) On the other hand, just as in RS12, we note another meaning for mc by observing that m = (2/ √ 3)mc = mH , if the characteristic gravitational angular frequency equals the angular frequency of a uniformly rotating object with mass M and angular momentum where LQM is given by (30). We stress that equ.(53) applies only to the ratio m/mc, and is not a bound per se on the particle mass m. As long as the ratio stays within that range, it doesn't matter whether we consider ultralight bosons, or bosons with the mass of the QCD axion, because R99,S in (40) decreases for increasing m.

Approximations to the FDM ground state density distribution
We have already mentioned that the SP system can be solved exactly only by numerical means, yielding the respective wavefunction ψ and the gravitational potential Φ coupled to it, and that the ground state solution, the soliton, is an attractor of a system in isolation. However, during FDM halo formation, multiple mergers of solitons can form bigger haloes, whose envelopes are the result of complicated wave dynamics (Schive, Chiueh & Broadhurst (2014a); Schwabe, Niemeyer & Engels (2016); Mocz et al. (2017)). Still, the centres of these bigger haloes remain soliton-like and that central profile has been "empirically" fit to the outcome of merger simulations. Schive et al. (2014b) introduced a function of the form with the central density ρc, the core radius c, and r the distance from the centre, in order to fit that central core region made up by the "central soliton". In fact, the core size parameterized by c is of order λdeB, if the virial (≈ circular) velocity of the core is used to evaluate λdeB.
Our analysis described in subsequent sections dedicates its focus to a sound analytical treatment of the question of vortices in selfgravitating FDM cores, thus we want to apply analytical models of their structure. To this end, we will use two different approximate models. Each one incorporates a test function for the mass density ρ = m|ψ| 2 , i.e. an approximate density profile of the FDM ground state, and reduce the search for a solution to solving the Poisson equation for Φ with given ρ in the first model, or relying on global quantities for the second model, respectively. While the profile in (58) has been used extensively in the literature, we decided not to use it for our analysis here. After all, it is an empirical fit to numerical data and, in our opinion, it still lacks a clear physical foundation. Instead, we will use a Gaussian profile to model the halo core, which we introduce below.

The Gaussian sphere
Our first approximate density model is inspired by the well-known Gaussian "wave-packet" of quantum mechanics. It is described as where σ 2 is the variance of the Gaussian, see Figure 1 for example plots. This model will serve as one of our approximations for the density profile of a vortex-free, bound FDM halo core in spherical symmetry, where ρ0 denotes the (vortex-free) matter density, ρc is the central density of the system and rs is the radial distance in spherical coordinates. Similar to the exact, numerical profile, this Gaussian approximation has no compact support. However, it is normalized as which yields The Gaussian profile (or ansatz) finds its justification in various previous analysis, where the Gaussian has been used in variational calculations, see e.g. Baym & Pethick (1996), Chavanis (2011), or Schiappacasse & Hertzberg (2018), as well as a fit to the exact soliton, e.g. in Marsh & Pop (2015) or Guzmán & Avilez (2018). Also, the asymptotic behaviour of Σ in (38), derived by Tod & Moroz (1999), shows that the exponential drop-off is very much faster than formula (58) would provide, i.e. the sharper fall-off is much better captured by the Gaussian. More recently, the usefulness and appropriateness of the Gaussian ansatz as an approximate model for the numerical soliton solution has been further shown in Padilla et al. (2020).
In the forthcoming, we will cut off the Gaussian profile at the radius R = R99,G, which includes 99% of its mass. This procedure is motivated by (40), and we will neglect the error introduced by effectively replacing M by M99, i.e. 99% of the total mass of the Gaussian. Given the finite size and in order to highlight the geometry that we will adopt, we will refer to this model as the "Gaussian sphere" (we will often use "R", "R99,G", or "RG" synonymously for the radius of the Gaussian sphere). Now, the mass-radius relationship of the Gaussian sphere has been derived in Chavanis (2011) and Padilla et al. (2020), and it follows easily, once we apply the virial condition, equ. (26), implying which differs only by a factor of ≈ 1.835 from (40). Of course, R99,G is not the same as R99,S (the former includes 99% of the total mass of the Gaussian, the latter includes 99% of the total mass of the numerical soliton). Now, if we insert R99,G into (46), we have m mc = 2.328 , thus smaller than (52), but see (53).
In later parts of this paper, it will be useful to express spatial variables of the Gaussian model in units of σ, and the dimensionless variables which result carry a tilde, i.e.
x ≡ x σ , where x stands for any spatial variable. In particular, for the cutoff radius R99,G we havẽ

The (n = 2)-polytrope
The second approximate density model to substitute for the exact, numerical soliton is inspired by considerations of Chavanis (2019), where he presents an argument as to why the density profile of an (n = 2)-polytrope is a particular solution to the hydrostatic version of Eq. (3) in the fuzzy limit, and the Poisson equation, Before introducing the conventional dimensionless variables θ and ξ and considering spherically symmetric configurations, the Lane-Emden equation (see appendix A) takes the form Setting n = 2 in Eq. (68), dividing it by √ ρ, applying the Laplacian operator ∆ and then substituting ∆ √ ρ on the right-hand side again through Eq. (68) itself, yields  (67), provided that However, this result should be handled with care since Eq. (68) implies (69), but (69) does not imply (68) due to the nature of operations between them. Owing to this non-equivalency, the polytrope of index n = 2 with fixed polytropic constant in (70), in particular its density profile is a valid approximation for our purpose, but it is not an exact solution of the GPP system of equations in the fuzzy limit. (Appendix A gives an overview on polytropic spheres and on the numerical approach to solve the Lane-Emden equation.) Figure 1 shows a plot of the (n = 2)-polytropic density profile (71). Note that we use again the notation ρc for the central density of the FDM core profile, like with the Gaussian sphere or Schive profile before. In addition to the above argument, we remind the reader that the mass-radius relationship in Eq. (40) is -up to prefactors-the same as the one for an (n = 2)-polytropic sphere. Denoting its radius as R0 and using equ.(70), we arrive at (see also eq.(120) or appendix A, eq.(A10)). The prefactor differs only by ≈ 1.033 from (64), thus in this respect the Gaussian sphere and the (n = 2)-polytropic sphere are very similar to each other. A closer look at the energies of a self-gravitating FDM object with the polytropic density profile (71) and the energy of an (n = 2)-polytropic sphere shows the following difference, see Chavanis (2019). We have seen that the quantum kinetic energy term is given by (21). Integration by parts yields The surface term vanishes since (71) indicates ρ ′ |r=R ∝ θθ ′ | ξ=ξ 1 and therefore ρ ′ (R) = 0 at the surface of the complete polytrope. By inserting (68) for n = 2, the quantum kinetic term takes the form The last equality follows from (25), the internal energy that arises from a polytropic pressure of index n = 2. In other words, we know that in the fuzzy limit there is no SI which would lead to an internal energy arising from a polytropic SI pressure. Moreover, we know that the total energy of a standard polytropic sphere, i.e. a spherical system in hydrostatic equilibrium with a density profile according to the solutions of the Lane-Emden equation (see appendix A), does not incorporate a "quantum" kinetic energy term. However, we have seen that certain operations and identifications allow to associate a polytropic density profile in the form of (71) with an FDM halo core density profile. Hence, one might guess that the quantum kinetic term corresponds to an internal energy of the form of eq.(25) and, according to Chavanis (2019), this is true up to a factor of 3/2. Unlike the exact, numerical "soliton" or the Gaussian profile, the density profile (71) has a compact support from the outset ("complete polytrope"), i.e. the density becomes zero at a finite radius. Furthermore, unlike e.g. the polytrope of index n = 1, used in the TF regime, the (n = 2)-polytrope cannot be represented in a closed-form expression, which seems counterproductive at first sight. However, we can benefit from previous calculations on approximate equilibrium solutions for uniformly and non-uniformly rotating polytropes carried out by Lai, Rasio & Shapiro (1993) (henceforth abbreviated LRS93) when constructing our second halo core model. Thus, we can exploit the sole fact that our second halo core model is based upon an (n = 2)-polytrope, in order to derive global energy expressions of uniformly (in the sense of constant angular velocity Ω) and nonuniformly (in the sense of superposed velocity fields) rotating FDM halo cores. This is described in the next section.

Preliminaries
From a cosmological point of view, we expect that angular momentum plays a crucial role in halo dynamics. The forming large-scale structure of the Universe imparts tidal torques onto growing haloes, which helps them to acquire angular momentum in their early phases.
Since BEC-DM is supposed to exhibit structure formation similarly to standard CDM over scales much larger than the de Broglie wavelength (see e.g. Schive, Chiueh & Broadhurst (2014a)), we adopt the same premise as in RS12, that we can make use of the so-called spin parameter which is the dimensionless ratio where E denotes the total energy, L the angular momentum and M the total mass of the object under question. Cosmological N-body simulations of standard CDM have been analysed, using this spinparameter in order to quantify the degree of rotational support of a CDM halo with net angular momentum L. Typical values for λ are in the range [0.01, 0.1] (see e.g. Antonuccio-Delogu et al. (2010)). In contrast, a self-gravitating, rotationally-supported disk has λ ∼ 0.4. Now, if the FDM object in question is a halo, or halo core, it is reasonable to expect that its value of spin parameter is similar to a typical CDM halo, as given in the range above. In particular, this implies that FDM haloes and their cores should not be rotationally supported.
Nevertheless, given our ignorance, we consider a somewhat broader range in λ than is typical for CDM, namely λ ∈ [0.01, 0.2]. In addition, and related to the idea that cores should be in approximate virial equilibrium, we require that the so-called t-parameter, which measures the importance of gravitational potential energy W over the energy of bulk motions T , stays below a certain threshold for our models, namely t ≡ T /|W | < 0.3. Now, the way we will assign a bulk angular momentum to our equilibrium models is to let them rotate with a constant angular velocity Ω, and the associated angular momentum will be used in eq.(75). BEC-DM haloes experience superfluid currents once Ω = 0, which will manifest themselves by a non-trivial phase function S (see also RS12). As long as no defects like vortices appear in the velocity flow, S is a smooth function which we will denote S0 and the respective velocity is v = ∇S0/m (see Eq. (10)).
One key point is that we study configurations in equilibrium. In other words, we consider stationary solutions given by (16) and (17), i.e. quantum states with time-independant observables which are eigenstates of the Hamiltonian, and add the distinction between different reference frames to our considerations. Rotating configurations at constant Ω in the rest frame of the object correspond to stationary solutions in the co-rotating frame. In this co-rotating frame, the bulk velocity is given by where primed quantities and variables denote those in the co-rotating frame.
In the next section, we will consider this co-rotating frame, assuming without loss of generality that the objects rotate about the z-axis, i.e. Ω = (0, 0, Ω). We will analyse whether and under what conditions vortex formation in these rotating, self-gravitating objects is energetically favoured. Our analysis requires a complete model of the equilibrium object, i.e. a geometry, density and velocity profile consistent with the demands of the GPP framework. In the following two subsections, we present the two complete models upon which our investigation is based, adopting a rotating Gaussian sphere in the first case, and a rotating, (n = 2)-polytropic, irrotational Riemann-S ellipsoid in the second case, and in the course of that we calculate the angular momenta L of objects which rotate at constant Ω, and derive useful relations which will be required in our energy analysis of the next section.

FDM cores as rotating Gaussian spheres
This model approximates the object under question as a sphere with radius R = R99,G = 2.576σ, rotating with constant angular velocity Ω = (0, 0, Ω) about the z-axis. The density profile is chosen according to the first density model, introduced in Eq. (59)-(60). This spherical object has a bulk rotation in the rest frame, Ω × r, hence the Gaussian sphere shows no net velocity in the co-rotating frame, Now, the total angular momentum of a uniformly rotating system whose axis of rotation coincides with an axis of symmetry, is generally given by where I denotes the moment of inertia. For a continuous body rotating about a specified axis, it can be written as if we decompose the position vector into a component parallel to the axis of rotation and perpendicular to the axis of rotation. The radial distance between each mass element and the axis of rotation is given by r ⊥ , provided that the centre of the mass distribution is located at the origin of the coordinate system. Thus, in spherical coordinates The total angular momentum of a uniformly rotating sphere filled with matter distributed according to our Gaussian density profile (59)-(60) is then given by with spherical coordinates (rs, θ, φ). This yields where is an expression with dimension [Length 5 ] and Erf(x) denotes the Gauss error function given by the integral

Comparison to the λ-spin parameter
We will use the spin parameter λ, Eq. (75), in order to derive meaningful angular velocities for our FDM halo cores for given λ ∈ [0.01, 0.2]. To this aim, we will express λ as a function of the angular velocity and the DM particle mass, λ = λ(Ω, m). We use global quantities in the rest frame of the rotating object, see eq.(19-23). The total energy is and Given the Gaussian profile in (59) and the definition in Eq. (83), the quantum kinetic energy term KQ amounts to The rotational kinetic energy incorporates the square of the bulk velocity in the rest frame, i.e.
v 2 = ( Ω × r) · ( Ω × r) = Ω 2 r 2 − ( Ω · r) 2 = Ω 2 r 2 ⊥ . This yields where the moment of inertia I, identified by expression (79), is given on the right-hand side (divided by Ω) of (82). Now, we need to calculate the gravitational potential energy W in (87), for which we need to determine the gravitational potential Φ0 first. Given our Gaussian ansatz for the density, we are required to solve the Poisson equation (4), thus Φ0 will be the solution of the following boundary value problem in spherical coordinates: where A is a constant. Imposing regularity, Eq. (90b), is effectively equivalent to requiring the absence of gravitational force at the centre. The calculation yields as the solution to system (90). Requiring a bound, isolated configuration, i.e.
sets A = −2πGρc/a and yields Inserting this solution along with the Gaussian density in (87) results in Collecting the expressions of (82), (88), (89) and (94), we get for the spin-parameter, The square of the spin-parameter can be written as with the gravitational angular velocity 7 defined as and M and R are the mass and radius of the rotating Gaussian sphere andR in (66). Furthermore, denotes a dimensionless quantity and B is given in Eq. (83). Inserting We can see that λ is a cumbersome function of (Ω,R, m/mc).
Since the final goal is to derive values for Ω from given λ-values in the range [0.01, 0.2], we calculated the roots of Eq. (101) for given λ 2 and m/mc. Since Eq. (101) is a non-linear function inΩ, several values forΩ may correspond to one λ-value, depending on m/mc. We deal with this multitude of solutions by selecting only those whose t-parameters fulfill t < 0.3. As an illustration for the numbers we get, let us consider e.g. m/mc = 2π, see (51)

FDM cores as (n = 2)-polytropic, irrotational Riemann-S ellipsoids
DM haloes presumably acquire angular momentum by tidal torques caused by large-scale structure, during their formation. Nonspherical equilibrium configurations are thus expected. In general, ellipsoidal shapes are expected for any object with non-vanishing angular momentum, yet spherical models are often employed due to their simplicity. Indeed, the rotating Gaussian sphere of our first model has the great advantage of being an analytical model for the ground-state, and this will be of substantial help in the next section, when we employ the energy analysis of the perturbed system with vortex. However, the Gaussian sphere suffers not only that is has a spherical shape, but also that it is not strictly irrotational in the rest frame. Therefore, we consider a further model, namely the irrotational Riemann-S ellipsoid, which is strictly irrotational in the rest frame (prior to any possible vortex formation), because on top of a uniform rotation, an internal velocity field is superposed, which combine to yield a vanishing net vorticity in the rest frame. Exact solutions exist for Riemann-S ellipsoids only for the case of uniform density, but fortunately, as we shall see below, LRS93 developed compressible (non-uniform) generalizations of Riemann-S ellipsoids (and other classical figures of rotation) by using their "ellipsoidal approximation", see appendix B. In fact, their approximate solutions for the compressible cases agree well with the true equilibria.
Irrotational Riemann-S ellipsoids have been applied to haloes in RS12 for the first time. RS12 showed that, in general, a rotating ellipsoidal halo cannot be both axisymmetric and irrotational, if it is nonsingular at the origin, and it was this insight which prompted their application of the irrotational Riemann-S ellipsoid. As an equation-ofstate for this Riemann-S ellipsoid, RS12 used an (n = 1)-polytrope, which is appropriate for the TF regime. The results of LRS93 were then used in order to calculate the global energies of such ellipsoids.
We have already established in the previous section, that the (n = 2)-polytropic density profile is an appropriate approximation to the actual "soliton" in the fuzzy regime, hence we paved the way for applying the results of LRS93 for (n = 2)-polytropic, irrotational Riemann-S ellipsoids, in this work. Nevertheless, we require some care, regarding the quantum-kinetic energy KQ of the ellipsoid, because this term is absent in LRS93 and RS12. Chandrasekhar (1969) revisits and elaborates on the incompressible Riemann ellipsoids -homogeneous fluid masses which maintain their ellipsoidal configuration under rotation and under their own gravity at all times with a velocity field that is given by the sum of a uniform rotation with angular velocity Ω and internal motions with uniform vorticity ζ ′ . In the case of so-called Riemann-S ellipsoids, the vectors Ω and ζ ′ are both oriented along the same principal axis and one can define sequences along which the ratio where Ω = | Ω| and ζ ′ = | ζ ′ |, is constant. The geometry of this non-axisymmetric object is given by its three semi-axes (a1, a2, a3) along Cartesian coordinates (x, y, z), or equally by its eccentricities Also, these ellipsoids fulfil a1 ≥ a3 ≥ a2 , i.e. they are all prolate bodies (see also Figure B2 in appendix B). Now, by applying an energy variational method and their ellipsoidal approximation (see appendix B for a summary), LRS93 developed generalized Riemann-S ellipsoids in the sense that they considered density and pressure profiles according to those of a polytrope of index n. As a result, the approximate internal and gravitational potential energy of polytropic Riemann-S ellipsoids can be written as and respectively, where the constants k1, k2 (depending upon n) and the dimensionless ratio f (depending upon the geometry) are given in appendix B, see eq.(B5, B6) and (B1). Here, RR denotes the mean radius of the ellipsoid 8 , given by The geometry and velocity field of this non-axisymmetric Riemann-S ellipsoid in equilibrium are closely related as follows. Let us denote unit vectors e1, e2, e3 along (x, y, z). One starts with an object that rotates rigidly with constant angular velocity Ω = Ω e3 about the z-axis. Then, one superposes an internal velocity field with uniform vorticity parallel to Ω, specified by the requirement that the resulting velocity vector at any point shall be tangent to the isodensity surface at that point. From this follows that the fluid velocity in the rest frame can be written as 8 The subscript R stands for "Riemann" and emphasizes that the system is being described by the mean radius defined in Eq. (106), as opposed to the Gaussian sphere. where In addition, there is a relation between the angular frequency of the internal motions Λ and the vorticity, By means of the above relations, and assuming polytropes of index n, LRS93 find for the angular momentum L and rotational kinetic energy T , and respectively. The moment of inertia I is given by and the definition of the constant κn and its values for n ∈ (1, 2) can be found in Eq. (A13). A bound BEC object described by the GPP framework is irrotational, whenever vortices are not present. In fact, the construction of Riemann-S ellipsoids allows to guarantee irrotationality of the object in its rest frame, given a certain choice of the parameter fR in eq.(102). This can be shown by introducing the circulation along the equator, and the vorticity in the rest frame Along the so-called irrotational Riemann-S sequence the ratio fR = −2 yields ζ = Γequator = 0, as required. Figure B1 in appendix B illustrates how the defining internal velocity field of the irrotational Riemann-S ellipsoid with n = 2 guarantees zero vorticity in the rest frame. RS12 provide similar visualisations of the velocity fields (their figure 2) for polytropic index n = 1. Two crucial equilibrium conditions can be obtained upon extremizing the total energy E = U + T + W of the Riemann ellipsoid with respect to all variations of the central value of the polytropic density profile and the axis ratios. For fR = −2, these conditions read and Here, we defined where the gravitational angular velocity of the ellipsoid is also defined as and RR is the mean radius in (106). The expressions for B12 and qn are given in the appendices, equ.(B4) and (A14), respectively. In other words, in case of equilibrium, one axis ratio determines the other one, thereby immediately fixing the geometry and furthermorẽ Ω. Moreover, equilibrium yields an important relation between the radius of the non-rotating spherical polytrope of same mass M , n and Kp as the Riemann ellipsoid, which is (see also (A10)), and the mean radius in (106). The relationship between the two radii reads A third equilibrium condition is the virial relation Eq. (26) which yields the total equilibrium energy see LRS93 for details. Inserting n = 2, the mean radius of the ellipsoid is where g(e1, e2) ≡ f (e1, e2)(1 − 2T /|W |), since f can be written as a function of the axis ratios, or equivalently as a function of the eccentricities. For increasing angular momentum (parameterized using λ below), the mean radius RR in (123) increases with respect to the spherical radius R0 (see also Table B1).
On the other hand, the total equilibrium energy for n = 2 is

Comparison to the λ-spin parameter
Now, we have everything in place to make the connection between Riemann-S ellipsoid and spin parameter. We can write λ as where follows from inserting (105) and (111). A division of Eeq by the gravitational potential energy amounts to and again, t ≡ T /|W |. Rewriting and multiplying the expressions (126) and (127) yields .
Thus, specifying the spin-parameter λ allows us to determine both axis ratios self-consistently by solving a system of equations consisting of relation (116) and (128). Appendix B demonstrates how setting a value for the spin-parameter λ effectively fixes the geometry of the ellipsoid and furthermore several dimensionless global quantities such asΩ in Eq. (117). Thereby, λ = (0.01, 0.05, 0.1, 0.15, 0.2) correspond to angular velocities Ω = (0.55513, 0.55659, 0.55266, 0.54376, 0.53113), respectively (see Table B1 for the corresponding t-parameters and other quantities). Note that these numbers vary not as much as the corresponding numbers Ω/Ωgrav of the first model at the end of subsection 4.2; moreover, they are not monotonic, due to the nature of the irrotational Riemann-S ellipsoid.

STABILITY OF ROTATING FDM HALO CORES TO VORTEX FORMATION
This section includes the most important parts of this paper, namely the study of the conditions for vortex formation within rotating FDM halo cores. We study the necessary and sufficient conditions for vortex formation, and we will show that vortices will not form for any FDM parameters. The failure of fulfilling the necessary condition excludes already some part of the FDM parameter space, while the rest is excluded by the failure of fulfilling the sufficient condition. We study these conditions for two different kinds of equilibrium models, and our conclusions are the same: solitonic cores are not subject to vortex formation. We will discuss in section 6 why this result has been anticipated and that it is in accordance with (and not in contradiction to) simulation results in the literature.

Necessary condition for vortex formation: L ≥ LQM
It is known from laboratory BECs and superfluids that a vortex will arise only, once an applied rotation surpasses a critical value. In our context, there is a minimum amount of angular momentum necessary for one singly-quantized vortex to appear, namely LQM of equ.(30). As in RS12, we first derive a relationship between the angular momentum L of our models and LQM . If we divide L in equ.(82) by LQM , we get the following relationship for the rotating Gaussian sphere where withΩ and Ωgrav defined in (100) and (97), respectively. The value of l(R) = l(2.576) = 0.69655 is fixed by our choice of R99,G in (66). Therefore, the ratio of L/LQM is determined solely byΩ, or equivalently by the product of mΩ/(mcΩgrav). Requiring implies a lower bound onΩ, namely For instance, if we compare the above values forΩ for m/mc = 2π at the end of subsection 4.2 to the lower bound in (133), we can see that the case λ = 0.01 does not fulfill the necessary condition for vortex formation for this choice of m/mc. Indeed, the corresponding angular momenta are L/LQM = (0.24138, 1.21484, 2.48472, 3.90358, 5.81385), and the first entry is smaller than one. There is a general trend as follows: the higher m/mc, the higher is L/LQM for a fixed λ, and the necessary condition for vortex formation is increasingly fulfilled. For the Riemann-S ellipsoid, we can proceed in an analogous manner: combining the relations (111) (1 − e 2 1 ) 1/3 (1 − e 2 2 ) 1/3 (2 − e 2 1 ) , withΩ in (117)(118). It shows that for fixed eccentricities -which also fixΩ -, and fixed polytropic index (here n = 2), the amount of angular momentum (in units of LQM ) depends only on the DM particle mass. Here, we expressed the DM mass in units of mc,R, where RG and mc,G refer to the radius and particle mass of the Gaussian values before, in order to compare on an equal footing (i.e. mc,G = mc, but we highlight now the different models). Different to the Gaussian sphere, the geometry of the Riemann-S ellipsoid is fixed by λ and therefore, only the DM particle mass is left to determine L/LQM . The latter grows for increasing m/mc,R or m/mc,G, respectively, i.e. a sufficiently high particle mass ratio is required to sustain at least one vortex. At higher masses, even more angular momentum can be provided. Of course, this general trend is the same as in the case of the (n = 1)-polytropic Riemann-S ellipsoid of RS12. In fact, if we were to plot equ. (134) for different values of L/LQM = (1, 10, 100) with m/mc,G as a function of λ, it would look almost indistinguishable from figure 3 (left-hand panel) in RS12, because L/LQM depends only very weakly on the polytrope index n. Thus, a range of high m/mc,G would correspond to a range of high m/mH in their notation. While the analysis in RS12 is valid only for high values of m/mH (or m/mc,G, respectively), this is not the case anymore for the FDM regime, which we study here. We have to content to small range in m/mc,G, see (53). In Figure 2, we show the ratio m/mc,G as a function of λ for which L/LQM = 1, according to equ. (134-135). The curve which results yields a lower bound on the particle mass, for a given λ-value of the Riemann-S ellipsoid, above of which the necessary condition for vortex formation is fulfilled. In short, if the size of the halo core is of the order of the equilibrium radius in the absence of rotation, Eq.(40), as we expect for small λvalues of interest for FDM haloes, then m/mc is limited to the lower end of the range in (53), and, in that case, L/LQM < 1 and the necessary condition for vortex formation is not met. Indeed, the smaller λ is, the larger is the value of m/mc required to make L = LQM , and, hence, further and further above the range of (53). The detailed thresholds depend upon the model (Gaussian sphere versus Riemann-S ellipsoid), but the trend is the same. This behaviour is a fundamental reflection of the nature of FDM cores and distinguishes them from the cores studied in RS12 for the TF regime, for which the size was related, instead, to the equilibrium radius for the (n = 1)-polytropes when SI pressure supports the cores against gravity, instead of the quantum pressure of FDM. In the TF case, therefore, the corresponding ratio m/mc was not restricted to a small range like (53) for FDM.

Energy analysis
Our energy analysis is based upon the comparison between the total energy of the unperturbed system (without a vortex) and the system with a vortex, all in the co-rotating frame. More precisely, the question will be whether the amount of angular momentum quantified by those spin-parameter values, which fulfill the necessary condition of vortex formation, is sufficient to make vortex formation also energetically favoured. Again, we stress that the vortex is required to be energetically favoured in the co-rotating frame of reference, where the wavefunction with vortex is stationary.
In the co-rotating frame, the angular momentum is given by and the GP energy functional in the fuzzy limit (g = 0) is In the course of the following energy analysis, we will use this expression for the energy in the co-rotating frame. For brevity, the primes on variables indicating the co-rotating frame will be omitted, except on phase functions and energies.

Energy splitting and vortex ansatz
For this model, we construct a wavefunction which accounts for a vortex in an otherwise unperturbed Gaussian sphere, calculate the total energy arising from this wavefunction by means of the GP energy functional in the fuzzy limit (137), identify those energy terms that arise due to the vortex and determine whether they lower or raise the total energy for given parameters of the system and the DM particles. We thereby assume that the Gaussian sphere with its vortex remains a viable approximate solution for SP. To this aim, the wavefunction ψ of the object in equilibrium is decomposed into an unperturbed and a vortex-part, as follows where "0"-indices indicate variables of the unperturbed system, ψ0 = |ψ0|e iS ′ 0 , and the vortex is included by means of the ansatz w = |w|e iS ′ 1 . Of course, the appearance of a vortex affects the density and the gravitational potential of the initially unperturbed (vortexfree) system. This raises the question how to identify the density and gravitational potential associated to the perturbation due to the vortex. The total density may be decomposed into and thanks to the linearity of the Laplace operator in the Poisson equation, we have The fact that the density of the unperturbed halo core is given by and hence ∆Φ1 = 4πGρ0(|w| 2 − 1) must be solved for the gravitational potential associated to the distortion of the density, brought about by the vortex. We apply the same method as in RS12 (derived in detail in their appendix B) to arrive at a convenient splitting of the energy contributions into vortex-free and vortex-carrying parts. According to that method, inserting the ansatz (138) into the energy functional (137) yields where is the vortex-free energy, while and are the energies due to the vortex. An FDM halo core with one central vortex is energetically favoured, as compared to an unperturbed, vortex-free one, if is negative, δE ′ < 0, i.e. the total energy of the system is reduced through the vortex. Our aim is to calculate δE ′ as a function of the parameters which define this model, and in order to do so, we apply our density models and a vortex wavefunction. For this halo model A 9 , we choose the Gaussian density (59)-(60) as the unperturbed background in (138), namely For the vortex in (138), we pick the ansatz in cylindrical coordinates (r, φ, z), with amplitude This is the same ansatz as in RS12; see also their figure 4, for d = 1. It corresponds to an axisymmetric, singly-quantized vortex along the axis of rotation (which we choose to be the z-axis) with vortex core radius s. The amplitude is dimensionless and the constant Cn will be given by a normalization condition 10 . The amplitude function (151) guarantees that the real and imaginary parts of the wavefunction tend to zero inside the vortex core, hence the density also vanishes. The form of the vortex ansatz (150)-(151) has some implications. First, it reflects the property that outside of the vortex, i.e. abruptly at r = s, the density is simply given by the unperturbed profile (149), and that there is a discontinuity. Moreover, each axisymmetric integration including the vortex wavefunction has to be split accordingly at r = s. Importantly, the above form of the vortex ansatz does not depend upon the regime we consider, i.e. it is valid in the fuzzy regime, as well as in the TF regime studied in RS12. The physics of the regime only enters, once we replace the characteristic vortex core radius s, e.g. with the healing length 11 . We will discuss this point in due course. Despite the fact that the global shape of halo model A is a sphere, the parametrization of the vortex demands cylindrical coordinates. As a result, the spherical domain over which the integration is performed is defined by the following intervals for the three cylindrical coordinates: where R is the radius of the Gaussian sphere, set by (66). First, the dimensionless constant Cn in Eq. (151) can be determined via the normalization condition of the unpertubed system, i.e. the vortex shall not change the overall mass of the system, Inserting (62) into (149) and using (151), we get 12 where Erf(x) is given in Eq. (84). It turns out that it will be useful in the course of calculating δE ′ to write C 2 n in the form of (157). The density of an FDM halo core with a singly-quantized vortex in its centre is then given by with C 2 n given by (156). As an illustration for FDM halo cores with vortex, we plot two-dimensional density profiles m|ψ| 2 = ρ/ρc for two choices of vortex core radii,s = s/σ = 0.8 ands = s/σ = 1.8 in Figure 3. They show that the vortex eats up the density in the very centre of the halo core, as expected from the vortex wavefunction. The white ring atr =s indicates a discontinuity of the overall density profile. This discontinuity can be understood by considering both one-sided limits, lim r→s ± ρ approaching r = s from above or below. According to Eq. (158), approaching ρ at r = s from below yields ρce −a(s 2 +z 2 ) C 2 n and from above yields ρce −a(s 2 +z 2 ) . These two expressions differ by C 2 n , which itself is a function of the vortex core radius s. However, since we will consider global energies only, this discontinuity poses no further problem.
Next, looking at Eq. (146), we see that we require to calculate the gravitational potential Φ1 associated with the perturbation of the density caused by the vortex ρ1.
The density perturbation due to the central singly-quantized vortex, is the source of the gravitational potential Φ1. Hence,

∆Φ1 = ∆Φ
12 Remember that the index n in Cn and Kn refers to "normalization" and is not to be confused with the polytropic index.
In fact, obtaining Φ1 is more complicated than was the determination of Φ0, the potential of the vortex-free Gaussian sphere in (90a)-(90c) with result in (93), because now we are looking for two different functions, Φ only valid outside the vortex (hence "(o)"). The general solution to Eq. (161) in cylindrical coordinates is where C1 and C2 are integration constants. We follow the same line of argument as in RS12 by imposing that the solution approaches a point-mass potential for large r at fixed z and for large z at fixed r. This yields where Msource is the mass of the assumed point-source. Basically, the idea is such that if we are far enough away from the source, i.e. the region in which the vortex-perturbation acts, its gravitational potential "feels" like a point-mass potential. However, how can Msource be determined? It is important to keep in mind, that the perturbed matter density inside the vortex region (159) is negative, which is comprehensible since the vortex removes, or rather redistributes the initial matter away from the vortex core region. Since this should be reflected by Φ where Mi is the mass inside the vortex core. In other words, the source of the outer-vortex potential associated with the perturbation of the density due to the vortex is set to be the negative vortex core mass. The mass inside the vortex core is given by and finally the outer-vortex potential is given by Finding an analytical inner-vortex solution for Φ (i) 1 posed a hard problem to our analysis. The corresponding partial differential equation (167) (or equivalently (169) below) seems to admit no closed-form solution as our attempts to find one has failed, or rather we found only a very approximate solution, based upon a multipole expansion, but its accuracy was determined to be too insufficient, so we will not reproduce it here. Therefore, we have solved the Poisson equation (167) numerically, and used that numerical solution in the forthcoming energy integrals. By multiplying both sides of Eq. (167) with and introducing the dimensionless variablẽ the Poisson equation in question can be written as with given in (157). Initially, Eq. (169) has to be solved forΦ However, since the density is an even function ofz, the potential will be too, and therefore we solved the system in the regioñ The boundary of that rectangular domain of integration in the (r,z)plane consists of four line segments connecting the four vertices (0, 0), (s, 0), (0,R), (s,R). There are two line segments connecting the vertices (0, 0) and (s, 0), and (0, 0) and (0,R) respectively, that lie within the vortex volume. Von Neumann boundary conditions, setting the normal derivative ofΦ (i) 1 to zero, were imposed on the differential equation along these line segments. Continuity requires that the solution of this partial differential equation matches with the analytical (closed-form) expression for the outer-vortex solution Φ along those two line segments. To arrive at Eq. (171), we have  (165) and (166) and inserted (156). The implementation of these boundary conditions returned a dimensionless, numerically interpolating function for our solution that we denote as Φs =Φs(r,z). A contour plot of this function can be seen in Figure 4. However, the calculation of the gravitational energy below requires the dimensional form, namely

Vortex energy
Now, the way is clear to proceed with the heart of the calculation, namely calculating δE ′ in Eq. (148) and determining whether it is negative or positive. We will consider each term in (146) and (147) separately.
The quantum-kinetic term, obviously requires to calculate the square of the absolute value of the gradient of the vortex wavefunction w, whereê φ andêr denote the azimuthal and radial unit vectors along the cylindrical coordinate directions φ and r, respectively. As a result, the absolute value of the complex vector field ∇w can be written as After splitting the axisymmetric integration domain according to (152), the quantum-kinetic term yields with the integral function in Mathematica). Let us briefly discuss this result. In the TF regime of strong SI, the quantum-kinetic energy due to the vortex always has a form ∼ ln (Z/s), where Z stands for a characteristic global length scale of the system. Likewise, the quantum-kinetic term in RS12, who incorporated a homogeneous Maclaurin spheroid for the halo (or halo core) with vortex, yielded a logarithmic term of that form, where Z was the length of the equatorial semi-axis of the Maclaurin spheroid, and s was the vortex core radius. In addition, Z/s ≫ 1 in the TF regime (compare also to Eq. (43)), and the quantum-kinetic energy due to the vortex is of leading order in the total energy. If we expand Eq. (176) for large R/s in a series representation, we would recover a similar logarithmic term, by means of the contribution of where γ is the Euler-Mascheroni constant. However, unlike in the TF regime, the quantum-kinetic term due to the vortex is not anymore of leading order in the energy in the fuzzy regime, because Z/s is not any longer much larger than one, an important feature of vortices in FDM that we will discuss in section 6. The second term of G ′ ρ 0 [w] in Eq. (146) is a gravitational potential energy term, which by itself does not include information on the vortex, and has already been calculated in subsection 4.2, given by Eq. (94). It is important to stress that the integration domain chosen for the calculation of this energy term is not exactly the same as the integration domain of all the other terms in this section. The reason for that lies in the splitting of the axisymmetric integration domain according to the vortex ansatz, where we integrate over a cylinder with height 2R and radius s and in addition over a domain, which we can visualize by imagining the result of shooting this cylinder through the centre of a sphere. With increasing s, the sum of the volume of this object and the volume of the cylinder differs a bit more from the integration volume of a sphere with radius R.
The fourth term of G ′ ρ 0 [w] contains the gravitational potential Φ1, associated with the distortion of the density by the vortex, and it is Now, we turn to the rotational energy in Eq. (147). First, we have already established via Eq. (77) that the unperturbed sphere shows no net velocity in the rotating frame, i.e. ∇S ′ 0 = 0. Finally, from (174) and Ω × r = Ωrê φ follows that where we have inserted the expressions for ρc and C 2 n , (62) and (156) respectively, in order to arrive at the last equality.
In conclusion, we have the energy difference between the total energy of the halo core with vortex and the total energy E ′ [ψ0] of the halo core without vortex, in the co-rotating frame, i.e. Eq. (148). In units of and using the characteristic mass introduced in (46), that energy difference reads as δE ′ ΩQM LQM = withΩ defined in Eq. (100). We can see that the energy difference in question is a function of several parameters, DM particle mass and halo model parameters, i.e.
Understanding the implications of this result requires some intuition, regarding the variables upon which δE ′ /(ΩQM LQM ) depends.
As we have already emphasized in section 1, in the fuzzy regime the vortex could in principle take up the whole halo, i.e. s R or in dimensionless notations R . Therefore, we will study δE ′ for the entire range 0 ≤s ≤R. In addition, plotting the energy difference due to the vortex for different particle masses provides even further insight, see Figure 5, left panel, where the contributions of the central singly-quantized vortex to the energy of the system is plotted as a function of the vortex core radius for different DM particle masses m for the lower bound on angular velocity, given in (133), i.e.Ω = Ω/ΩQM = 1.43564. It is evident that for this minimumΩ, FDM is not able to form a central vortex in the context of halo model A for any of the considered parameter values, since δE ′ /(ΩQM LQM ) > 0 everywhere. Before we proceed in adopting values forΩ, which correspond to spin parameters of interest, let us examine first the features of the vortex energy curves: The divergence fors → 0 is due to the quantum-kinetic energy term (176) and it is only weakly dependent on m. The m-dependent gravitational potential energy contribution, which does not include any information on the vortex in Eq. (94), merely fans out the individual curves; this term just adds an m-dependent constant value to the energy. The third term of G ′ ρ 0 [w], which combines the unperturbed gravitational potential with the total density ρ0|w| 2 does not show a strong dependence on the vortex core radius s in the given range. Its effect on the total vortex energy δE ′ /(ΩQM LQM ) amounts to a slight decrease with increasing s. The rough overall shape of the curves and their dependence on the vortex core radius s is to a great extent the result of the quantum-kinetic energy of the vortex and the gravitational potential energy due to the vortex potential Φ1. If we were to plot their sum, in units of (ΩQM LQM ), as a function ofs = s/σ for different particle masses m, and compare it to the left panel of Figure 5, it would reveal that these figures look quite the same; the disregarded terms only yield a slight horizontal stretching of the energy curves. Therefore, we established that, in the fuzzy regime, the vortex energy is strongly dominated by the quantum-kinetic and the gravitational potential energy due to the vortex. Now, the striking dependence of the gravitational potential energy generated by the vortex potential Φ1 on the vortex core radius itself, see Eq. (179), is shown in Figure 6, where the two integrals corresponding to the inner-and outer-vortex contributions are plotted separately. While the monotonic increase of the inner-vortex contribution (left panel), first integral in (179), for configurations with a larger and larger central vortex is easily comprehensible, the sdependence of the outer-vortex contribution (right panel), second integral in (179), is quite striking. In principle, this second integral should decrease with increasings since in that case the integration domain decreases, given a fixed halo core sizeR. However, the sdependence of the inner vortex mass (165) counteracts this general trend initially, yielding a local maximum.
Last but not least, there is the rotational energy term (180) which, in units of (ΩQM LQM ), yields the term −Ω in our energy difference (182). This term is responsible for a global shift of the energy curves in the left panel of Figure 5, downwards to lower energy values. This makes sense, because increasing angular velocities should make the vortex more favourable, in principle. Vortex formation would be shown to be energetically favoured, as soon as one of the curves would cross the abscissa. Therefore, we need now to consider angular velocities which correspond to spin parameters of interest.
But first, we revisit considerations of section 3 regarding the gravitational healing length ℓgrav. We have already mentioned that ℓgrav, as well as the corresponding healing length ℓ in the TF regime, can be regarded as the distance over which the wavefunction tends to its background value when subjected to a localized perturbation. For this reason, RS12 proceed in the course of their vortex energy analysis in the TF regime by replacing the vortex core radius s with ℓ, based upon the assumption that s is very close to the healing length ℓ, motivated by the vortex ansatz. In just the same spirit, we may proceed with the assumption that the vortex core radius in the fuzzy regime is well approximated by the gravitational healing length ℓgrav, given in Eq. (47). That is, i.e. for givenR = 2.576, the vortex core radius is now a function of the boson mass. Thus, replacings in (182), according to (185), constitutes the first step in order to arrive at the right panel of Figure 5. In subsection 4.2, we have calculated the relationship between the spin parameter andΩ, see Eq. (101). Once we employ that relationship, the angular velocity Ω corresponding to a given λ now depends on the particle mass m. We calculate the contributions of the central singlyquantized vortex to the energy of the system, δE ′ /(ΩQM LQM ), again as a function of m/mc by settings =lgrav, fixing the value of λ and calculatingΩ for every m/mc in the considered range. Finally, the right panel of Figure 5 shows the result of this procedure for λ-values (0.01, 0.05, 0.1, 0.15, 0.2). We can immediately see that δE ′ /(ΩQM LQM ) > 0 for all parameters considered. The plot also shows that vortex formation is even less energetically favoured for higher m/mc. Thus, while the necessary condition for vortex formation is better and better fulfilled, the higher m/mc (i.e. the necessary condition is fulfilled for an increasing range of λ), vortex formation is less energetically favoured for higher m/mc. Now, we checked whether a change of the definition of our gravitational healing length in equ.(45) affects our conclusion. The definition in (45) already picks the higher-end value, so we repeated our calculation for values of lgrav which are a factor of 0.5 and 0.1 smaller, respectively. As a result, the associated vortex core radii s are assigned smaller values, and we are effectively getting closer to the asymptotic energy range of small s with its quantum-kinetic energy, which diverges for small s. This implies even much larger (positive) values for the vortex energy δE ′ /(ΩQM LQM ), moving it orders of magnitude to higher values. Therefore, we can safely conclude that FDM halo cores as rotating Gaussian spheres are not subject to vortex formation for any value of allowed spin parameter λ.

Halo model B: irrotational Riemann-S ellipsoid transitions to Gaussian sphere with vortex
In this section, we will also employ an energy analysis in order to verify whether the improved modelling of the rotating equilibrium object would change the conclusion of the last section. In fact, our conclusion does not change: vortex formation is not favoured. To this aim, we will apply an energy argument, which is similar but not equivalent to the approach described in the previous section, however, it is also inspired by a similar consideration of RS12. Let us consider some unspecified dynamical process in the course of which the halo core may form one singly-quantized vortex in its centre. The actual occurrence of this vortex creation requires that the process in question transforms the initial state of the halo core into a final state with lower total energy compared to the initial state. As a result, we are going to look at two "snapshots" of that transformation, namely the initial and the final configuration, calculate their total energies and compare them for a given set of parameter ranges. The underlying halo core in the initial state will be modeled by an irrotational Riemann-S ellipsoid along the lines of subsection 4.3. In the context of model B, however, we will assume that this process transforms the halo core from a vortex-free configuration to a configuration with one central vortex: the vortex takes up all of the angular momentum of the system, once it has formed. Thereby, the initial ellipsoidal shape of the halo core disappears, it becomes spherical, allowing us to draw on the perturbed Gaussian sphere as a model for the second snapshot or final state. In the course of the following calculations, it is important to distinguish clearly between quantities and global properties of the initial Riemann-S ellipsoid, which will be denoted by the index "R", and the final vortex-carrying Gaussian sphere, denoted by the index "G".
We start the analysis by considering the total energy of the vortexfree initial configuration in the rest frame which can be written as where T and W are given by (112) and (105), respectively. On the other hand, KQ has no classical counterpart and is therefore absent in the studies of Chandrasekhar (1969) or LRS93. However, the internal energy of a Riemann-S ellipsoid is given by (104) and can be written as due to the ellipsoidal approximation, where ρc,R is the central density of the initial Riemann-S ellipsoid incorporating a polytropic density profile. We have already seen that the internal energy of a sphere arising from an (n = 2)-polytrope is related to the quantum-kinetic energy via Eq.(74), given the density profile (71) in section 3. Owing to the ellipsoidal approximation of LRS93, the total internal energy of a rotating polytrope is identical to that of a spherical one with the same central density. This implies that the quantum-kinetic energy can be written as with Kp given by (70) and n = 2. Before we can continue with the calculation of the total energy ER, the current analysis requires to establish relations between several quantities of the initial and final state of the halo core. While the underlying initial halo core configuration is modeled as an irrotational, (n = 2)-polytropic Riemann-S ellipsoid, whose respective energy terms LRS93 derive from initially spherical polytropes modified by rotation, we set the final halo core to be a sphere -a Gaussian sphere, not a polytropic sphere! However, LRS93 provide us with relation (123), i.e. the process we consider shall transform the system in such a way that it settles with the radius of the final Gaussian sphere, RG, which we associate with R0. In this sense, we heavily rely on the fact that the density profile of an FDM halo core can be approximated by both, the Gaussian profile and the (n = 2)-polytropic profile. Thus, we require and RG = R99,G is our cutoff radius (66). What about the central densities of the initial and final configurations? Given the central density ρc,s of the spherical polytrope with radius R0, we have the relation between the central and mean density of a polytrope, see appendix A. The analogue for the polytropic Riemann ellipsoid can be written as The halo core mass M is required to be conserved during this transition. Thus,dividing Eq. (191) by Eq. (190) yields ρR =ρsg(e1, e2) 6 and ρc,R = ρc,sg(e1, e2) 6 , where we have used relation (189). In the context of this model, the vortex-free Riemann-S ellipsoid shall transition to a Gaussian sphere with vortex. Hence, we set although, unlike the spherical polytrope, the Gaussian profile has no compact support, i.e. ρc,G is the central density of an infinitely extended system which we cut off at R99,G . Once more, we rely on the presumption that the two density models are equally appropriate and closely related. As a result, the total internal energy of the Riemann-S ellipsoid can be written as U = k1 2π 9 where we have used (70) and k1 in (B5) has to be evaluated for n = 2. Consequently, the quantum-kinetic energy of the FDM halo core, as a Riemann-S ellipsoid, in units of ΩQM,RLQM , is given by where we have used (119), (135), and (189), as well as Ωgrav,R ΩQM,R = √ 3 2 m mc,R = √ 3 2 m mc,G g(e1, e2) −1 . (197) By means of these relations, the gravitational potential energy (105) for n = 2, in units of ΩQM,RLQM , amounts to Finally, the rotational kinetic energy (112) reads as where κ2 can be found in appendix A, Eq.(A13), and the definition of the dimensionless factor h(e1) arises from setting fR = −2 and combining the relations (102) and (110):  (206) and logarithmically plotted) as a function of DM particle mass m/m c,G . The gravitational angular velocity is determined according to Eq. (117) for given polytropic index n = 2, f R = −2 and the axis ratios are listed in Table B1 for different λ. Right-hand-plot: Energy difference δE ′ RG /(Ω QM,R L QM ) as a function of spin-parameter λ for fixed total angular momentum L = L QM (light, dotted curve). For comparison, the particle mass m in units of m c,G at fixed L = L QM is shown again (dark, solid curve); see also Figure 2, but note the difference between linear and logarithmic scaling.
of the system. This implies i.e. the DM particle mass is now a function of the halo core geometry. In other words, given a shape of the halo core, requiring the amount of angular momentum to be enough in order to sustain just one vortex yields a condition on the particle mass. Meeting this condition implies that the particle mass is fixed for given λ, or equally for given axis ratios of the halo core. Hence, this condition can only be met at one single point when considering the energy difference as a function of particle mass (compare to the left panel of Figure 7), i.e. fixing the halo core geometry by λ immediately yields one value for δE ′ RG /(ΩQM,RLQM ). This can be seen in the right panel of Figure 7, where δE ′ RG /(ΩQM,RLQM ) and m/mc,G are each plotted as a function of λ. Although the energy difference decreases with increasing spin-parameter (for fixed L/LQM ), it remains much too high in order to make energetically favourable the formation of a vortex ever. Thus, we have shown that vortex formation is also not favoured in the context of halo model B.

CONCLUSIONS AND DISCUSSION
We have studied gravitationally bound halo structures made of ultralight DM bosons. These structures are referred to as SFDM haloes or BEC-DM haloes, whose bosons are described by a single scalar wavefunction. These haloes can be modelled, using the Gross-Pitaevskii-Poisson system of equations, which can be written in the form of quantum-mechanical fluid equations. The analysis in this work focused on the so-called fuzzy regime (fuzzy DM, or FDM), in which quantum pressure balances gravity. In order for the wave nature of SFDM to be potent on galactic scales, the de Broglie wavelength of bosons (evaluated with the virial velocity of the object) has to be of same order of magnitude than the global size of the system. In fact, ground-state solutions of the Gross-Pitaevskii-Poisson system (also called Schrödinger-Poisson system for FDM) have a size of roughly the de Broglie wavelength and constitute attractor solutions, so-called "solitons", of these equations. They were the subject of our investigation in this paper. These ground-states describe the "solitonic cores" of large haloes, as well as the entire halo, if the latter only consists of the core (e.g. appropriate hosts for the smallest galaxies).
The main objective of our analysis was the study of rotating FDM haloes and halo cores, and to study the question of vortex formation. The underlying equations of motion allow for vortex solutions; within vortices the density goes to zero and the velocity diverges. As such, they can change the smooth halo background, with dynamical and structural consequences. In general, vortices are manifestations of the quantized vorticity in superfluids and they are the building blocks of quantum turbulence, phenomena which have been also studied for BEC-DM haloes. As a quantum superfluid, BEC-DM is irrotational -vorticity-free, but, while starting from irrotational initial conditions, it can become unstable to vortex formation; outside of vortices, the rest of the system remains vorticity-free. This was demonstrated for the case of haloes and halo cores that form in BEC-DM with repulsive SI in the Thomas-Fermi (TF) regime, when angular momentum is present, as expected during large-scale structure formation, in Rindler-Daller & Shapiro (2012). In particular, the virialized halo cores in that case, in which gravity was balanced by SI pressure and rotation, were found to be unstable to vortex formation for a large range of boson parameters, mass m and SI coupling strength g. However, since this instability required a minimum SI strength, we suggested there that the small SI regime (and especially FDM without SI) would be found to be vortex-free inside halo cores. In Rindler-Daller & Shapiro (2012), we also introduced approximate, analytical equilibrium solutions for rotating BEC-DM halo cores, in particular we applied for the first time irrotational Riemann-S ellipsoids to BEC-DM haloes, for which an (n = 1)-polytropic density profile was used, appropriate for the TF regime. These irrotational Riemann-S ellipsoids are useful models to describe halo cores which do not form vortices.
In Rindler-Daller & Shapiro (2014), we further suggested that the polytropic SI pressure support that set the size of these halo cores in the TF regime would be supplemented on larger scales by the wave-motion-support generated by the wave nature of BEC-DM and its quantum pressure, during the virialization of haloes that assemble from infall and mergers, making it possible for haloes to be much larger than their polytropic cores in which only SI dominates. Later, BEC-DM without SI, i.e. FDM, with boson masses around m ∼ 10 −22 eV, has been studied with greater detail, including simulations that report all haloes have solitonic cores of the size of the de Broglie wavelength (as evaluated inside haloes), supported against gravity by quantum pressure. BEC-DM haloes also show a wave-supported envelope outside this solitonic cores, with a profile that resembles that in CDM haloes, but in which wave motions provide the random internal motions responsible for virial equilibrium to be dicussed shortly.
In this paper, we extended the analysis of Rindler-Daller & Shapiro (2012) and study the question of vortex formation in rotating FDM halo cores (without SI) in equilibrium, with the same analytical rigour than in this previous work. We stress that, in the fuzzy regime, all characteristic length scales are roughly comparable, including the size of perturbations of the system, like quantum vortices, which poses a hard problem for analytic studies. As in RS12, we associate the typical amount of angular momentum, expected from large-scale structure, with the spin parameter. We considered two different analytic models for rotating, solitonic cores, first a Gaussian profile cut off at a finite radius, we call it the Gaussian sphere, and second an irrotational Riemann-S ellipsoid with an (n = 2)-polytropic density profile. The latter equation-of-state was shown to be a viable approximation for FDM solitons. We thereby established that irrotational Riemann-S ellipsoids can be used in the fuzzy regime, as well.
We studied the necessary and sufficient conditions for the formation of one centrally located, singly-quantized vortex in FDM halo cores. The sufficient condition is studied by employing a detailed energy analysis. We found that, unlike the TF regime, for solitonic cores of FDM, vortex formation cannot be triggered by angular momentum as it is in the TF regime, for neither of our models that we considered. For vortex formation to be triggered by angular momentum, the specific angular momentum must first satisfy a necessary (minimum) condition, that it exceed the minimum value that gives each particle an angular momentum of . If this necessary condition is satisfied, then it is further required that vortex formation be energetically favoured, in order to establish that vortex formation will take place. In the TF regime, both conditions can be met for the typical amounts of specific angular momentum for cosmological haloes, for a large range of m and g. However, for FDM, we have shown in this paper that the necessary condition is generally not met for typical amounts of halo angular momentum. We have further shown that, even for angular momentum which is large enough to meet the necessary condition, vortex formation is nevertheless not energetically favoured. This is consistent with and can explain the fact that simulations of structure formation in the FDM model do not find vortices in the solitonic cores of FDM haloes. Now that we have described the implications of our results, let us finally compare them in more detail to some of the previous simulation works. The cosmological simulations of Schive, Chiueh & Broadhurst (2014a) first showed convincingly that in the centre of gravitationally bound FDM haloes, one finds coherent standing waves, i.e. stable solitonic cores. Schwabe, Niemeyer & Engels (2016) simulated the dynamics of these solitonic cores by investigating binary and multiple mergers of up to 13 such cores. Again, it is found that solitonic cores are embedded within bigger haloes, whose outer density profile declines like NFW density profiles. Furthermore, Schwabe, Niemeyer & Engels (2016) find that their emerging solitonic cores are rotating ellipsoids, if the system is initialized with non-zero total angular momentum. In fact, the respective volume rendered images and velocity fields of the cores strongly indicate that they resemble irrotational Riemann-S ellipsoids, as the authors point out. The same conclusion was drawn by Edwards et al. (2018) in their study of the dynamics of solitonic cores. While the previous work by Rindler-Daller & Shapiro (2012) had introduced Riemann-S ellipsoids, their study was limited to the TF regime. Now, the analysis in this paper has firmly shown that irrotational Riemann-S ellipoids can be used in the fuzzy regime, as well. Therefore, our work has shown that, in general, polytropic, irrotational Riemann-S ellipsoids provide useful analytical counterparts for the formed solitonic halo cores of BEC-DM halo formation simulations.
Moreover, vorticity is generated during structure formation (from vorticity-free initital conditions), but only outside of the solitonic cores, as found in Schwabe, Niemeyer & Engels (2016) and Mocz et al. (2017). The origin of this vorticity has not been wellstudied, but its absence from solitonic cores is consistent with the results of Rindler-Daller & Shapiro (2012), as well as with our new results of this paper. More precisely, Mocz et al. (2017) present a set of 100 numerical simulations in which a group consisting of 4 to 32 solitonic cores merge and form one final halo whose core is (like the final bound cores of Schwabe, Niemeyer & Engels (2016)) well-fitted by the density profile introduced by Schive et al. (2014b), see Eq.(58). However, in contrast to the simulations by Schive et al. (2014b) and Schwabe, Niemeyer & Engels (2016), the work of Mocz et al. (2017) specifically includes the study of quantum turbulence, found in the envelopes of their final haloes. Given the soliton fit (58), they find that the break between the soliton profile and the outer NFW-like profile within the final virialized BEC-DM halo occurs universally at ≈ 3.5rc, which approximately corresponds to the soliton radius. Regarding turbulence and vortices exhibited by the final halo, Mocz et al. (2017) conclude by analysing the energy power spectra E k , the radial energy density profiles and 2D slices of the wavefunction amplitude |ψ| of their 100 simulations. Granules and turbulence appear everywhere in the domain, except for the central solitonic core. The stable solitonic core remains free of substructure and turbulence. The radial energy density profiles show that the quantum-kinetic energy supports the structure up to 2.7rc. Beyond that radius, all three energy contributions become comparable, yielding a characteristic signature of turbulence, namely equipartition. The energy power spectra lack power for small k, show a mode which displays most of the turbulence, and finally follow a k −1.1 power law for large k. This resembles the spectrum of thermallydriven and hence isotropic turbulence of superfluids (Tsatsos et al. (2016)).
Furthermore, Mocz et al. (2017) show that the power spectra of their simulations peak at 2π/k peak ≈ 7.5rc which corresponds to a scale of twice the soliton radius. This explains why the filamentary distribution of the |ψ|-field (outside the soliton) show preferentially soliton-sized granules. Now, Mocz et al. (2017) explain the absence of vortices (and quantum turbulence, as a result) within solitonic cores by referring to their very equilibrium properties, as ground-state solutions of the fundamental equations. However, this assessment is not accurate. For example, the analysis in Rindler-Daller & Shapiro (2012) showed that vortices do arise in equilibrium halo cores in the TF regime, for a large parameter space, although these cores also correspond to ground-state solutions of the fundamental equations. The analysis in this paper suggests that vortices add energetic "penalty" to the equilibrium of the cores, hence are strongly unfavoured, while the addition of a positive, large-enough particle self-interaction overcomes this penalty and vortices can be favoured, as shown in Rindler-Daller & Shapiro (2012). Moreover, BEC-DM in the fuzzy regime requires m/mc 2, see Eq.(53) (corresponding to small m/mH in Rindler-Daller & Shapiro (2012)), and for the lower mass end of this ratio, even the necessary condition for vortex formation is not met, not even for high spin parameters. Thus, the FDM parameter space excludes vortices in halo cores, because either the necessary, or the sufficient condition fails to be fulfilled.
More recently, Hui et al. (2020) followed up on the question of vortex formation in FDM haloes. They use analytical arguments, as well as numerical studies of merging, isolated Gaussian "peaks" to model solitons. They were interested to model the outer regions of virialized BEC-DM haloes, that exhibit quantum turbulent dynamics due to vortex tangles. Indeed, similar to the papers discussed above, Hui et al. (2020) do also not find vortices within the central cores of their simulated haloes. We noticed that the "universal" properties of vortices, prominently listed in their abstract, have been found already in previous works (some of which we cited earlier, incl. our own work). Defects are investigated analytically in the absence of gravity, before gravitational effects of vortices are included in their numerical simulations. They argue that vortex gravity is not of great importance. However, since all characteristic length scales in FDM (unlike the TF regime) are of similar order (e.g. perturbations characterized by the vortex core radius s can be as large as the system size, in principle), the quantum-kinetic energy of the vortex need not be the leading-order term, as we have shown in our analysis. Our results suggest that the gravitational potential energy due to the vortex is not any less important than its quantum-kinetic energy, and should not be neglected. Moreover, Hui et al. (2020) make a point in their discussion that solutions exist with angular momentum that do not carry vortices, and that realistic haloes would be supported by velocity dispersion, rather than rotation. Of course, it is true that we do not expect rotationally supported FDM haloes or halo cores, either, a premise which we also made in this paper. Systems with angular momentum need not have vortices (long known for laboratory systems), and this was one key result in Rindler-Daller & Shapiro (2012), when irrotational (vortex-free) Riemann-S ellipsoids were shown to be viable approximate solutions for TF cores with spin parameters similar to CDM, for model parameters which would not favour vortex formation.
To re-iterate, we have shown in this paper, using analytical methods, that angular momentum by itself will never be sufficient to create vortices within the cores of FDM haloes, in perfect agreement with our earlier work Rindler-Daller & Shapiro (2012) and with simulations of FDM halo formation.
We wish to dedicate this work to the memory of Ernst Dorfi .