Ultra-Light Dark Matter

Ultra-light dark matter (ULDM) is a class of dark matter models (DM) where DM is composed by bosons with masses ranging from $10^{-24}\, \mathrm{eV}<m<\mathrm{eV}$. These models have been receiving a lot of attention in the past few years given their interesting property of forming a Bose-Einstein condensate (BEC) or a superfluid on galactic scales. BEC and superfluidity are one of the most striking quantum mechanical phenomena manifest on macroscopic scales, and upon condensation, the particles behave as a single coherent state, described by the wavefunction of the condensate. The idea is that condensation takes place inside galaxies while outside DM behaves like a normal cold particle DM. This wave nature of DM on galactic scales that arise upon condensation can address some of the curiosities of the behaviour of DM on small scales while maintaining the successes of LCDM on large scales. There are many models in the literature that describe a DM component that condenses in galaxies. In this review, we are going to describe those models and classify them according to the different ways they achieve condensation. For that, we review the phenomena of BEC and superfluidity, and apply this knowledge to the DM in order to explain their construction and phenomenology. We describe the small scale challenges these models aim to solve and how ULDM alleviates them. These models present a rich phenomenology that is manifest in different astrophysical consequences. We review here the astrophysical and cosmological tests used to constrain those models, together with new and future observations that promise to test these models in different regimes. We finalize by showing some predictions that are a consequence of the wave nature of this component, like vortices and interference, that could represent a smoking gun in the search of these rich and interesting alternative class of DM. (Abridged)


Introduction and Motivation
An overwhelming amount of observational data provides clear and compelling evidence for the presence of dark matter (DM) on a wide range of scales. This component, which is believed to be responsible for the "missing" mass in our universe, is the main ingredient for all the structures we have in our universe. This is one of the oldest unsolved problems in cosmology, being traced back to the 1930's [1,2], and also one of the best measured ones. The evidence for dark matter first emerged from the study of the rotation curves of galaxies. From pioneering studies [3], it was already evident that the amount of matter necessary to fit the flat observed rotation curves did not match the theoretical curves predicted, assuming Newtonian mechanics and taking into account only the visible matter present in those galaxies. Dark matter was proposed as an additional (non-luminous) component to explain this discrepancy.
Nowadays, the evidence for dark matter comes from precise measurements on a wide range of scales. From sub-galactic and galactic scales, to clusters, going up to the large scale structure (LSS). On cosmological scales, the observed anisotropies of the Cosmic Microwave Background (CMB) [4], together with data from Type Ia Supernovae, determine the total energy density of matter with high precision. This together with the bounds on the abundance of the light chemical elements from Big Bang Nucleosynthesis, which constrains the amount of baryonic matter in the universe, strongly shows the need for a clustering component of non-baryonic 1 origin, that does not interact (strongly) with photons, and that dominates the matter content of the universe, accounting for approximately 85% of all matter. The same non-luminous and clustering component is necessary to explain the structures we see in our universe today, as is evident in observations of the large scale structure of our universe [5,6].
With all this evidence coming from precise astrophysical and cosmological observations, cosmologists have converged to a phenomenological model to describe our universe, the ΛCDM model. This model is currently the concordance model of cosmology and it accumulates a number of observational successes. It exhibits outstanding agreement with current cosmological observations [5], which is manifested in the parameters of this model being DE Ultra-light DM curious challenge is the regularity/diversity of rotation curves. One thing that is surprising about galaxies is that they are extremely diverse, but at the same time they are incredibly regular. This fact is manifest in several empirical scaling relations, such as the well known Baryonic Tully Fisher relation (BTFR) [9,10]. The BTFR shows the correlation between the total baryon mass (including stars and gas) of the galaxy with the the asymptotic rotation velocity in galaxies. The measured BTFR follows a scaling relation different from the one predicted by ΛCDM, and it holds for a range of 6 orders of magnitude in mass, with very small scatter. The significance of these discrepancies is disputed and addressing these challenges is an active field of research. Those challenges emerge on scales where baryonic physics is relevant and simulations including several baryonic effects have been perfected pointing in the direction that baryons could possibly explain some of these observations within ΛCDM.
As the physics of these baryonic processes is complex and as there is no final consensus about the status of theses discrepancies, an alternative explanation for these discrepancies on small scales could be that DM is not the usual CDM, but a more complex component that has different dynamics on small scales. With that, the dynamics on small scales can offer a chance to probe the properties of DM in the hope to help find the microphysics of this component. The simplest modifications that can address some of the small scales discrepancies is warm dark matter (WDM) [11], where DM has a small mass leading to a thermal velocity dispersion, small enough so it does not spoil the large scale prediction of CDM. Even with a small velocity, DM free streams out of potential wells and is enough to suppress the formation of small scale structures addressing some of the small scale problems. Another popular model inspired by those discrepancies is the self interacting DM (SIDM) [12], where the DM particles interact in a way to also suppress the formation of structures on small scales.
In the past few years, another class of alternative models has (re)emerged offering a different mechanism to explain the dynamics on small scales. In this class of models, DM forms a Bose-Einstein condensate (BEC) or a superfluid on galactic scales. The idea is that the wave nature of DM on astrophysical scales modifies the dynamics on small scales, while maintaining the successes of CDM on large scales. In these models DM is composed by ultra-light bosonic fields that are produced non-thermally in the universe. An example of such a non-thermal relic is the axion, and the connection between axions and the phenomenon of Bose-Einstein condensation on small scales attracted a lot of attention to these models. However, these models are much more general and offer a rich phenomenology.
There are many ways that such a macroscopic quantum phenomena can be modelled in galaxies. For this reason there are many different models in the literature that aim to address the small scale challenges using different realizations of these phenomena. We classify them in this review into three categories, according to the different quantum systems they describe. These three classes are the fuzzy dark matter, which forms a BEC on galactic scales; the self-interacting BEC DM, and superfluid DM, where these last two classes describe a superfluid that forms on galactic scales, but with distinct dynamics. Each of these categories have different properties which lead to different astrophysical consequences, that can be probed by current and future astrophysical experiments. For this reason, this class of models has regained interest in the community in the past few years, with new and exciting experimental effort to probe many aspects of the small astrophysical scales, opening the avenue to test these models and answer some questions about the nature of DM.
Motivation for this review and detailed plan: This review has the goal of giving a short overview of the ultra-light dark matter (ULDM) candidates. This class of models is described in the literature with many different names and different constructions. Sometimes the names of these models are used interchangeably, but one needs to be careful about the distinctions in those models where condensation (and the consequent phenomena) are reached from different methods. The current literature lacks a review that includes all these models, describing their differences and similarities. This review thus attempts to not only give a very general and brief overview of the motivation for proposing this class of models, but one of its goals is also to unify the language and classify the models according to the different quantum phenomenon they describe.
These classes of ULDM models are based on very well understood and well studied phenomena in condensed matter physics. Understanding their description and differences is particularly important in order to understand the difference between the classes of ULDM models that exist in the literature. For this reason, we include here a brief review of BEC and superfluidity, and the different descriptions of these phenomena.
With that, we aim to give a general picture, although not complete, of the state of the field to date. We hope this review can be a resource to researcher entering this exciting field.
The review is organized as follows. First, in Sec. 2, we start by describing the small scale challenges of ΛCDM, as a motivation to show the discrepancies these alternative models of DM aim to address. Next, in Sec. 3, we introduce the basic concepts of the quantum phenomena of BEC and superfluids, so we can, in later sections, apply those concepts for the case of DM. Following this we are ready to describe the ULDM models in Sec. 4. We start by describing the three classes that we propose to classify the models of ULDM based on the type of quantum system they describe. We then talk about the fuzzy DM and the self-interacting BEC DM models, showing the conditions for them to condense on galactic scales. We then focus on the fuzzy DM model, showing how and in which conditions the model attempts to solve the small scale challenges, and the interesting astrophysical consequences this class presents. We then talk about the superfluid DM model describing its condensation on galactic scales, the formation of the superfluid core and its observational consequences. We also discuss the stability of this construction, and its possible extension to cosmology.
Since there is no unique general literature this review is based on, but a series of articles referring to specific topics, the main references used are cited in the corresponding sections. The only exception is Section 2 that is based mainly in the following reviews [13,14,15]. The review will not address axions specifically, although axions are a subset of the ULDM models and share most of the gravitational consequences cited here. However, axions have a whole rich phenomenology in its interactions to baryons that will not be explored here, but that can be seen in the following reviews, e.g. [16,17].
Conventions: In the entire review natural units are used, where c = = 1, unless stated otherwise. The exception is section 3 where all the factors are present. With that, the reduced Planck mass is given by M 2 pl = 1/8πG, where G is the Newtonian gravitational constant. Unless stated otherwise, the metric signature used is (+, −, −, −), and Greek letters are indices going from µ, ν = 0, . . . , 3. Sometimes for simplicity we describe partial derivatives as ∂ µ = ∂/∂x µ . In the text gray boxes bring definitions necessary for the understanding of the topics in the section or following sections. Frames text and equations refer to important results or discussions that we would like to highlight.

Small Scale Challenges of Cold Dark Matter
In the concordance model of cosmology, dark matter is described by the CDM paradigm. This hydrodynamical description for DM is in very good agreement with observations from large scales. This can be seen in the power spectrum (P (k)), which probes the two-point correlation function of the matter density perturbations, represented here by the dimensionless power spectrum where k is the wavenumber of the fluctuation, shown in Figure 2. The large scales (around k 0.1 Mpc −1 ), measured by the CMB and LSS galaxy surveys, show a good compatibility with the CDM model. This agreement is also robust for the non-linear intermediary scales (k ∼ [10 −1 − 10] Mpc −1 ) with constrains from clusters, weak lensing and Ly-α forest. As we go to smaller and highly non-linear scales (k 10Mpc −1 equivalent to M 10 10 M ), these constraints are less strong, and might retain important information about possible deviations from the CDM paradigm. We can see on the right side of Figure 2 on galatic and subgalatic scales, different models of DM would behave very differently from the expected linear behaviour of CDM and this could be probed by the observations on those scales [18,19].
On small scales, the formation of structures is highly non-linear and the evolution of structures is studied using large numerical simulations. In the past few years, those simulations improved in size and precision, simulating the cosmological and small scales. But when compared to the observations of galaxies, a number of discrepancies emerged, revealing some curious behaviour on small scales. Given the enormous success of the concordance model, these discrepancies attract a lot of interest of the community. They might represent that we need to better take into account the astrophysical processes that happen inside those regions, which indeed have a complex dynamics. Or this might indicate that the CDM model is not good to describe the physics on small scales and the coarse grained CDM paradigm needs to be revised. An even a more radical approach would be to modify gravity on smaller scales.
In this section we present very concisely the theory of non-linear structure evolution. We show how the numerical predictions assuming the concordance model might be in tension with the current observations of galaxies. These tensions are seen in the counts and density of low-mass objects, and in the scaling relations that show the tight regularity that galaxies present. We highlight in this section some of the concepts that are going to be used in the ULDM section and that might not be too familiar for researchers from fields of dark matter phenomenology and cosmology. . The gray dotted horizontal line represents the limit from linear to non-linear regime, where ∆ ∼ 1. The power spectrum for ΛCDM and for WDM were generated using the Boltzmann code CLASS [20,21], and for the fuzzy DM using AxionCAMB [22,23] 3 .
Milky-Way (MW) galaxy: MW is a barred spiral galaxy and part of the Local Group of galaxies with mass ∼ 10 12 M . It has a stellar disk of approximately 30 kpc in diameter and 0.3 kpc thick, and v vir ∼ 100 km/s (virial velocity, defined below), with the halo of the MW being hundreds of kpc in size.
Dwarf galaxies: Dwarf galaxies are low luminosity, small size galaxies, with masses smaller than 10 9 M . Regarding their mass, they can be further divided into: Bright dwarfs (M ∼ 10 7−9 M ), classical dwarfs (M ∼ 10 5−7 M ), and ultra-faint dwarfs (M ∼ 10 2−5 M ). Regarding their characteristics, they can be divided into ellipticals, spheroidal and irregulars, that contain gas and star formation.
Dwarf Spheroidals (dSphs): Type of dwarf galaxy with a close to spheroidal shape, they have lowluminosity with a very small quantity of gas and dust, and no recent star formation. They present a large amount of DM and are usually the satellites.

Dark matter halos and substructures
A halo can be defined as a virialized spherical mass concentration of dark matter. Halos are formed by gravitational collapse of a non-linear overdense regions that stopped expanding to collapse into a sphere in virial equilibrium 5 . The 4 The masses are indicated in terms of the solar mass M which is equivalent to 2 × 10 30 kg in SI units. Distances are denoted in parsec (pc), where 1 parsec corresponds to 1 arcsecond of measured parallax, and it corresponds in SI units to 3.1 × 10 16 m. 5 Virial equilibrium means that it obeys the virial theorem E kin = −2 Epot and conservation of energy. So we can describe the system only in terms of the radius R and thee mass, M (or velocity V ) of the spherical mass concentration.
virialization of the halo happens through violent relaxation, where the DM particles scatter on small fluctuations of the gravitational field present in this distribution, taking a time t dyn , the dynamical time, to fully cross the sphere. Once this process is completed, at t coll , the dark matter halo has a radius approximately 1/6 of the radius of the region it collapsed from, and average density [24] ρ = (1 + δ vir )ρ(t coll ) , (2) whereρ is the mean density, and (1 + δ vir ) ≈ 178 Ω −0.6 m . Given this, the dark matter halo is defined as the spherical region where the density is approximately 200 times the critical density of the universe at a given redshift, with mass given by: where ρ cr = 3H 2 (z)/(8πG). The virial velocity is given by the circular velocity at the virial radius, We can see from these expressions that halos that form early in the evolution of the universe are less massive, while late-forming halos are more massive and larger.
This definition is not unique and depends on the choice of the virial overdensity parameter 6 , ∆, which above was taken to be ∆ = 200ρ cr /ρ. More generally, (3) can be written as M vir = (4π/3) R 3 vir ∆ρ. The values of ∆ can vary in the literature, with some common definitions being ∆ = 333 at z = 0 for a fiducial cosmology given by [4], which asymptotes to ∆ = 178 at high-z [25]; or a fixed ∆ = 200 at all redshifts, usually denoted by M 200m .
We identify the DM halos from numerical simulations and can extract from them the abundance of halos as a function of their mass for a given redshift. The individual halos can also be analyzed in those simulation and the radial mass profile can be determined. A surprising feature encountered in those simulations is that halos appear to have a universal density profile, averaged over spherical shells. Their functional form is characterized by the Navarro, Frenk and White (NFW) profile [26], ρ NFW (r) = ρ s (r/r s ) (1 + r/r s ) 2 → 1/r , for r r s 1/r 3 , for r r s where r s is the radius where the slope of the profile changes and ρ s = ρ(r s ). We can see that this profile diverges towards the center of the halo, presenting a cusp. The amplitude of the density profile can be written in terms of R 200 , as we can see from (3) where x = r/R 200 , and c := R 200 /r s is the concentration index and describes the shape of the distribution. With that, the NFW profile can be determined completely by R 200 (or M 200 or any other halo radius definition), and the parameter c. The shape of the concentration can be inferred from simulations, where c ∝ (M/M * ) −1/9 (1 + z) −1 . We can see that, early-forming halos have a smaller radius, and they are denser than the larger ones, given the higher concentration. The NFW profile can be generalized for a three-parameters profile that better fits the DM profile of halos for all ranges in mass [27,28,29].
Surface brightness of galaxies Our capacity of observing galaxies is limited by the brightness of the sky. They can only be observed if their surface brightness, which is the brightness per area since they are extended objects and not point sources, is higher than the sky surface brightness 7 where the area is the area of the survey (µ B = 23 mag/arcsec 2 ). This can limit our understanding of the distribution of galaxies, making us miss the fainter ones. Therefore, we have the following nomenclature for the galaxies with respect to their surface brightness.
Low surface brightness (LSB) galaxies: There is no formal definition for LSB galaxies, but in general they are disk galaxies that have surface-brightness smaller than µ B . They are believed to make the majority of the galaxies in our universe, and most of the LSB galaxies are dwarf galaxies. However, this is not necessarily the case, with LSB galaxies being galaxies in a broad range of masses, and very diverse morphologies. This low luminosity is likely associated to a small star formation rate in those galaxies. So those galaxies are believed to be DM dominated. Their rotation curves usually reach much smaller speeds than the ones from high surface brightness galaxies (see below), with a very slow rise before reaching the plateau region given their lower density, but broaden mass distribution.
High surface brightness (HSB) galaxies: They are usually defined as galaxies that are brighter than µ B . They are the usual galaxies we study. The rotation curves are known to reach high velocities with a steep rise, coming from the inner region that has a higher density of baryons with narrower mass distributions than LSBs, which is described by the Newtonian baryonic acceleration. This is followed by a Kleplerian fallout to the flat part of the rotation curve.

Discrepancies in comparison with observations
In this section we will show how some of the theoretical predictions from simulations of the small scales considering the ΛCDM model compare with respect to astrophysical observations. However, this comparison is not straightforward, since we indirectly probe the dark matter inferring it from the visible matter that traces the gravitational potential of galaxies and clusters. There are a few approaches to connect the information of galaxies and the dark matter halos like forward modelling, abundance matching and kinetic measurements, and each of those methods has its difficulties and limitations. The result of this comparison is a series of discrepancies that challenges the results of the simulations, and in some cases limitations in observations. I will present some of these challenges in this section. Some of those challenges might have complementary origin and solution, and are indeed connected, as we will discuss bellow.

Cusp-core
As we saw above, the expected density profile from colisionless simulations is the NFW profile which is cuspy towards the central region of the halo. Given the complex dynamics of baryonic matter in some galaxies, good laboratories to probe the halo structure are low surface brightness (LSB) galaxies and late-time dwarfs. Those systems are dominated by DM throughout their halo up until the central regions. Measuring the rotation curves of dwarf galaxies, refs. [31,30] found that those measurements preferred cored isothermal profiles. Many other measurements of the rotation curves of those systems [32,33,34,35,36,37,38,39,40,41,42,43] have confirmed this discrepancy, showing that a constant density core with a profile with a slope γ = 0 − 5 (considering the profile at small radius given by ρ ∼ 1/r γ ). The smallest values for this slope from dissipationless simulations are too large in comparison to the ones obtained by observations.
The recent measurement of nearby dwarf galaxies from the survey THINGS (HI Near Galaxy Survey) [44] and LITTLE THINGS [45] confirmed this discrepancy. Measuring the rotation curves from 7 and 26 nearby dwarfs, they found that the inner slope is much smaller than the NFW one (γ = −1), with γ = 0.29 ± 0.07 for the LITTLE THINGS survey, as we can see in Figure 3.
The situation is more complex for high surface brightness (HSB) objects given its complex inner density structure; or for galaxies with large mass, like spiral galaxies, where at small radii is dominated by baryonic matter; and for dwarf spheroidals that show a strong dependence on the bias in the modelling of the system. Even in the case of dwarf galaxies, it was pointed out that some systems present cuspy profiles, while others cored ones, presenting an unexpected diversity in the rotation curves [46]. Since different results were obtained by different techniques for the same system, this shows that determining the inner slope of galaxies is a hard task.
The origin for these discrepancies can come from the fact that the simulations take into account only DM, while the properties of galaxies are also influenced by the presence of baryons. The newest hydrodynamical simulations obtained by many independent groups have shown that baryonic feedback can in fact soften the inner cusps in the profile and generate core-like profiles like the ones observed for dwarf galaxies. The main effects are supernova feedback flattening and dynamical friction from baryonic clumps (for a more detailed list of these and other baryonic processes see [14]). These simulations show a threshold mass of M vir ∼ 10 10 M below which the simulation predict profiles that are cusped [47,48,49,50,51,52].
However, not all simulations agree with this result. Additionally, modelling those baryonic feedback effects is challenging, and introduce many new parameters and uncertainties in modelling assumptions. Finally, not all Fig. 3 In this plot we show the results from the THINGS survey [44]. The plot shows a comparison of the density versus the radius, normalized by R 0.3 and ρ 0.3 , of the theoretical parametrizations of the NFW potential (solid line) and the pseudo-isothermal (dashed line), with the simulated galaxies, represented in the plots by the legend DG1 and DG2, with the observational data. The observational data from 7 dwarf galaxies measured by THINGS is represented by the other points. We can see that the galaxies seem to follow a cored profile, while NFW predicts a cusp.
baryonic processes that might influence the formation and dynamics of galaxies were included in the simulations, and that might reveal to be important for the result. It is clear that the inclusion of baryonic effects is hinting in the right direction, but until consensus is achieved, alternatives need to be considered. As mentioned before, a modification of the properties of DM might in a simple way account for that, as we will show for the case of Bose Einstein condensate DM. An early solution to the cusp-core problem, and that explains the rotation curves with exquisite precision is a modification of the dynamics of gravity on small scales, the MOdified Newtonian Dynamics (MOND). This is also a solution for the regularity versus diversity challenge, and its main points and shortcomings will be presented at the end of this section.

Missing satellites
Structure formation is hierarchical in nature and it is expected that the DM halos are also populated by small subhalos. This is confirmed in ΛCDM simulations of Milky-Way sized halos, which show that the subhalo mass function diverges toward low masses, limited only by the numerical limit. Those simulations then predict several hundreds of subhalos with v max ∼ 10 − 30 km/s, that are large enough to host a galaxy (M peak ≤ 10 7 M ). On the other side, until 2005 only 12 MW classical satellites were known, with 15 more confirmed ultra-faint satellite galaxies until 2014, with the data from Sloan Digital Sky Survey (SDSS) [54]. To date, with the inclusion of Dark Energy Survey (DES) data, a few more ultra-faint candidates were discovered, with the known count of satellites of more than 50. However, the number of MW galaxies satellites is still much smaller than the number predicted from simulations. This is known as the missing satellites problem, and not only appears in the MW, but also in the Local Group.
DES and future observations are expected to discover more of those ultra-faint galaxies, which can alleviate this discrepancy, but there is still a debate if this will solve the problem. Another possibility is that low-mass subhalos are there, but we just cannot see them, since they have very low baryonic content. One can expect that for low mass subhalos, galaxy formation is suppressed since the photoionizing background heats the gas, reducing its cooling rate and inhibiting gas accretion for M vir ∼ 10 9 M [55,56,57,58,59]. Star formation is also suppressed since supernova-driven winds could strip the gas out of these halos [60]. Other mechanisms can also suppress the baryon content in the low-mass galaxies, see [14], like reionization suppression. So, the visible subhalos are only a set of the entire distribution of halos that contains the non-visible faint end. It was shown recently in hydrodynamical simulations [61,62,63] that apparently this mechanism can solve the difference in the number of predicted and observed satellite galaxies, thus solving the missing satellite problem. But the question remains if this process needs to be too finely tuned to solve the problem.

Too big to fail
The above mechanism that could solve the missing satellite problem leads to another challenge: the too big to fail problem. When we say that the visible subhalos of the MW are only a set representing the most massive subhalos in the total distribution of subhalos, to have agreement with ΛCDM simulations these visible MW subhalos need to correspond to the most massive subhalos predicted by the simulations. But, the most massive subhalos predicted by those simulations have central masses 8 (V max > 30km/s) that are too large to host the observed satellite galaxies [65,66], and the ones that have central mass like the expected by the MW (with 12 < V max < 25km/s) are not the most massive ones. So, the puzzle is why should the most massive subhalos, where the gravitational potential is the strongest and the striping gas mechanisms cited above are not important, be too big to fail to form stars and galaxies? This problem also appears in the galaxies in the Local Group and Local Volume [68,69], so it is not a specific property of the MW.
This higher central mass from the most massive sub-halos predicted by simulations in comparison to the MW dSphs, seems to be a more general feature that appears in simulations. This discrepancy might indicate a more delicate issue related to the internal structure of the sub-halos. In this way, the too big to fail problem is more than just the problem related to the missing satellites problem, as stated above.
Like for the other problems, it was proposed that some astrophysical processes driven by baryons could be important on those scales and solve the too big to fail problem. However, these solutions seem to only work for the MW and for very efficient feedback, like the supernova feedback that only solves the too big to fail problem if very efficient. This is an intense topic of debate and no consensus appears to have been reached. As these notes were being written, there has been claims that the the too big to fail problem has been solved [72].
As for the cusp core problem, different DM physics could solve those problems by having a mechanism that suppresses the formation of small scale subhalos, and that reduces the central densities of massive subhalos (or modifies the dynamics of the central regions). We are going to show how the models with Bose Einstein condensation address some of those problems.

Diversity vs. regularity: scaling relations
Although our universe came from very smooth initial conditions, nowadays the diversity of galaxies that we find in the universe is extraordinary. This incredible diversity of galaxies, though, presents a surprising regularity. This fact is manifest in several scaling relations, that are shown to hold very tightly for a diverse range of galaxies. These relations relate the dynamical and baryonic properties of galaxies, and hold even for DM dominated systems, and they are one of the most tantalizing aspects of galaxy phenomenology, representing the most pressing challenge for ΛCDM on small scales.
The most famous of those relations is the Baryonic Tully Fisher relation (BTFR) [9,10], which relates the total baryon mass (including stars and gas) of the galaxy to the asymptotic circular velocity in galaxies, V f (this is the velocity measured at the flat portion of rotation curves): where a 0 is the critical acceleration, a scale that appears in observations. Its value can be obtained from the data and given by a 0 ∼ 1.2 × 10 −8 cm/s. The BTFR expands the regime of validity of the Tully Fisher relation which relates the luminosity, instead of the total mass, to the circular velocity. Luminosity is a probe of the stellar mass, and in the BTFR the observed gas mass is also considered on top of the stellar mass. This extends the validity of the scaling relation by many decades in mass. This empirical scaling relation is shown to hold for large ranges of masses, 6 generations, with a very small scatter, compatible to the size of the error bars. The left panel of Figure 4 presents the BTFR. As we can see, the slope of the BTFR is different from the one predicted by ΛCDM, There is another general scaling relation that also displays the interesting behaviour of galaxies: the mass discrepancy acceleration relation (MDAR). This is more general since the BTFR can be obtained from the MDAR at large distances in the disk. The MDAR is a relation between the gravitational acceleration from baryons alone (g bar ), from the distribution of gas and stars in galaxies [74,75], and acceleration inferred from rotation curves (g obs = V 2 /r). As it can be seen in the right panel of Figure 4, this scaling relation shows a remarkably tight correlation between these quantities for very diverse and large number of galaxies This can be seen by comparing the interval determined by the solid red lines and the uncertainty in each individual point, represented by the red uncertainty bars on the top figure, with the dashed lines show that the data is compatible with negligible scatter.
This relation shows us that in regions of high acceleration, where g obs > a 0 and baryons dominate, one has g obs ∼ g bar . For low accelerations, in the central regions where it is expected to be DM dominated, this relation deviates from the unit line. This suggests a very curious behaviour: the baryon mass distribution dictates the  [15], which shows the relation between the baryonic total mass (M b ) and the asymptotic circular velocity. Dark and light blue points represent star and gas dominated stars, respectively. The dashed line represents the relation expected for ΛCDM, with slope equal to 3; while the dotted line which better fits the data, has slope 4. Right panel: Plot of the Radial Acceleration Relation from [74], for 153 SPARC galaxies. The fit to the data is given by the solid line while the dotted line is the unit line. The insert is a histogram of the residuals. The red uncertainty bars represent the uncertainty in each individual point. The lower panel shows the residuals, and the red uncertainty bar shows the mean uncertainty on individual points. The dashed lines represent the rms value in each bin and the solid red lines represent the observational uncertainties and variation between the stellar mass-to-light-ratio from galaxies.
behaviour of the rotation curve at all radii, even for the regions expected to have less baryons. And this behaviour holds even for galaxies that are DM dominated.
These empirical relations, coming directly from observations, show the surprising feature that in galaxies the dynamics is dictated by the baryon content, even when DM dominates. Even more unexpected these relations are very tight, showing very little spread, even if they come from very diverse types of galaxies. As pointed out in [13], in these correlations what dictates the dynamics is the baryon mass, which is the sum of gas and stars, and not only the stellar mass, which is the one that is expected to correlate more with the total feedback energy.
Recently it was shown by many groups that these relations can be explained within the ΛCDM paradigm [85,86,87,88,89] using the latest hydrodynamical simulations like EAGLE [95,96], APOSTLE [97], Illustris [94], ZOMG [90,91,92], and NIHAO [93]. Those simulations include several baryonic effects (like star formation, stellar evolution, metal enrichment, gas cooling/heating, galactic outflows and BH feedback, among others) to their ΛCDM simulations 9 . Those new large volume and very accurate simulations, like Illustris and EAGLE, have also been able to reproduce the features of the rotation curves of galaxies within ΛCDM. This large amount of progress in the simulation side is very encouraging.
However, some questions still remain. While the BTFR and the MDAR trends can indeed be reproduced by those simulations, it is pointed out by most of the authors that the scatter obtained in the scaling relations is larger than the one expected from data (some authors claim that this spread is correlated with the errors in the stellar feedback). The question remains though whether this is a matter of improving the feedback models and/or resolution of the simulation, or if given the stochastic nature of the feedback effects, they will ever be able to give such tight correlations. Another point that is important to be answered is about the importance of these baryonic feedbacks, since some of those authors claim that the results are not very sensitive to the feedback model, which is intriguing. The way those effects are introduced in the simulations is by parametrizing their effects, instead of introducing all of these feedback mechanisms from first principles, which is understandable given the complexity of each of those phenomena. This includes many new parameters to the simulations. And different simulations might use different parametrizations. Those simulations also still do not go all the way until dwarf galaxies 10 , which are DM dominated and where most of the tension is. In summary, this is a very challenging and exciting field and a lot of progress has been done on the simulation side with results that are very encouraging to explain the formation and dynamics of galaxies. But there are some uncertainties in those results and the simulations still do not fully reproduce the observations.

What the small scales tell us
As we saw above, the small scales hold precious information that can help us understand astrophysical processes, or even the nature of DM. This is revealed by the challenges presented above, which show rich dynamics on galactic and sub-galactic scales. There are a number of ways that these discrepancies can be addressed. Within ΛCDM, this can done by including baryonic effects, which as we saw in the previous sections seem to address partially or completely some of those puzzles. Another proposal for solving some of the puzzles of galactic evolution is more radical and proposes a universe without DM that has a modified force law for small accelerations, the MOdified Newtonian Dynamics (MOND).
A third avenue is to modify the DM paradigm. Different models of dark matter can affect the formation of structures in distinct ways, both in the linear and in the non-linear regimes, as we can see in Figure 2. Therefore the small scales offer an opportunity to probe the microphyics of DM, beyond the hydrodynamical large scale CDM paradigm. The non-linear regime can be specially changed by modifications of this paradigm, with galaxies being a sensitive probe of DM. This could help find new properties of DM, that could help elucidate its nature.
MOND empirical law -Milgrom, in 1983 [73, 79, 80], motivated by the scaling relations and rotation curves of galaxies, made a remarkable observation about the mass discrepancy in galaxies. He observed that the mass discrepancy can be determined by the observed baryonic matter, and can be described by the simple empirical law, where a N,b = GM b (r)/r 2 is the Newtonian acceleration due to baryons. The scale a 0 appears naturally from observations, like we saw in the previous subsection, and its value can be fitted by the data 11 giving a 0 ∼ 1.2 × 10 −8 cm/s. This scale separates the regimes where the centripetal acceleration experienced by a particle is given purely in terms of the Newtonian (baryonic) acceleration at large acceleration, and at small acceleration, by the geometric mean of a N,b and a 0 . The relation works very well fitting the rotation curves of galaxies, both HSB and LSB galaxies. LSB galaxies (which were predicted by Milgrom), are DM dominated, or in the language of MOND, have low accelerations, given their low stellar surface density. It is also remarkably successful in explaining the empirical scaling relations (for a review see [81,82,83,15]). More importantly, this empirical relation reveals a very interesting and curious fact. It seems that the dynamics in galaxies is driven by the baryons, even for galaxies that are DM dominated. This seems to indicate a long range correlation between baryons on galaxies.
This fact made Milgrom think that a fifth force was responsible for this correlation, instead of DM, and that this relation could come from a modification of gravity at those scales. In order to try to get the empirical law (7) as a modified gravity theory, Bekenstein and Milgrom [84] described an effective theory for MOND, which we will call full MOND. This can be accomplished by having a scalar field coupled to gravity with effective Lagrangian which represents a scalar field with a non-canonical kinetic term that is conformally coupled to matter. This Lagrangian, for static and spherically symmetric source, results in a modified Poisson equation 10 To my knowledge. But as I said, it is a fast moving field.
where in the second equation the spherical symmetry was assumed with denoting derivative with respect to the radial coordinate. This theory describes that on top of the Newtonian force, there is a scalar field mediated force, which is given by the MONDian acceleration. This is a simplified version of their theory, since in their theory they have a way of making an interpolation between the different regimes. This theory also has a fractional power kinetic term, which might be problematic. We will discuss this in Section 4. Nowadays, with the precise observations on large scales, specially from the CMB anisotropies and lensing observations, any theory that does not have DM is not compelling. This full MOND theory cannot explain galaxy clusters, since it does not predict an isothermal profile. Many attempts were made to extend MOND, by including DM, to try to explain the observation on scales larger than galactic, or extending it to relativistic regimes (see reviews cited above). However, the empirical relation (7) is incredibly successful. That alone, without the assumptions of full MOND (no DM and modified gravity), even in the context of ΛCDM, is a powerful statement about how DM is distributed in galaxies: in regions where baryons dominate, the theory behaves like Newtonian theory, and in regions where the DM dominates, the DM mass is uniquely determined by the baryonic distribution, Given the shortcomings of the full MOND, but the great successes of the empirical law, instead of trying to obtain this theory from a fundamental Lorentz invariant theory, the idea is to obtain the MOND dynamics from a theory of DM. In this way, MOND dynamics emerges only at galactic scales while maintaining the CDM behaviour on large scales. This is achieved in the theory of DM Superfluid that will be presented in Section 4.
We discussed above some alternatives to CDM, like WDM and SIM, that can address some of the discrepancies on galactic scales. And this is the idea of the ULDM classes of models that we are going to present in this review. In general, the goal of these models is to have a DM component where the wave nature of DM, characterized by the formation of a Bose Einstein condensate on galactic scales, leads to the formation of a homogeneous core in the central part of the galaxy, leading to a suppression of the structure formation on sub-galactic scales, addressing some of the small scale discrepancies. DM forms a coherent state in the central core of galaxies, modifying the physics on these scales. More specifically, to describe the scaling relations one of the classes of models of ULDM, superfluid DM, has the goal of providing an effective modification of the dynamics on small scales, reproducing the MOND empirical law from this modified theory of DM. The idea to use modifications of DM is to offer a global explanation for the scaling relations, maintaining the tightness of those relations.
We have many different ways in which this can be implemented by having pure Bose Einstein condensates or superfluids with different types of interactions, and each of those can yield different modifications of the small scales and also different predictions. Those models rely on the quantum mechanical effect of Bose Einstein condensation and superfluidity, applied to a context very different than where they were discovered, in the realm of condensed matter physics. In the next section we present a small review of the physics of Bose Einstein condensate and superfluid as an introduction to the section where we describe dark matter as having those properties.

Bose-Einstein Condensation and Superfluidity
In this section we present a short review of Bose-Einstein condensation and superfluidity. We aim to give an introduction to the basic theory, properties and the methods used to describe those systems, so they can be applied for the case of DM in the next section. The different description of those systems and their limit of validity are very important to be able to understand the construction and validity of the DM models presented next, and why they present different phenomenologies and astrophysical consequences.
Bose-Einstein condensation is one of the most fascinating phenomenon of quantum mechanics. Since it was theorized in the year of 1920's, by Satyendra Nath Bose and Albert Einstein, its experimental realization opened the door for many advances in the physics of many-body systems, and even to the application of this phenomenon in other fields like cosmology. Its first experimental realization was done in 1995 by two independent groups using laser and magnetic cooling device to cool down rubidium atoms gas [98,99]. Nowadays, BECs are observed in helium, ultra-cold atomic gases, quasi-particles in solids, multi-component (mixtures) of BECs, among other systems.
Following on the works of Bose [100], which described the quantum statical properties of photons, Einstein extended this concept to a gas of non-interacting particles of integer spin, later called bosons as a tribute to Bose, that follows a Bose-Einstein statistics [101]. This Bose gas has the property that at low temperatures a large number of these bosons, described then as quantum oscillators, condense into the lowest momentum state, exhibiting long range coherence. This physical phenomenon initiated the idea of Bose-Einstein condensation (BEC).
A BEC is defined as a system where at very low temperatures a large fraction of the bosons of the system occupy the lowest energy state of a configuration. This macroscopic occupancy of the ground state is an inherently quantum mechanical phenomenon. Physically we can interpret it as a consequence of the wavelike nature of these particles at low temperatures, where the de Broglie wavelength of these bosons is larger than the inter-particle separation, and their wavepackets superpose and form a coherent macroscopic wavefunction describing the entire system. The BEC is then described by a single wavefunction of the system, linking to the long range coherence property of a condensate.
A few years after BEC was theorized by Einstein, another intriguing macroscopic quantum mechanical phenomenon was discovered: superfluidity. In 1937, Pyotr Kapitsa [102] and independently John Allen and Donald Misener [103], conducting experiments with helium-4 realized that after cooling down this liquid to a certain temperature, the fluid starts flowing without friction, even climbing the walls of the container where it was stored. Fluids that exhibit this behaviour, characterized by a zero viscosity, are called superfluids. Landau provided a phenomenological description of this effect which rendered him the Nobel prize in 1962. It was proposed by Fritz London [105], after the development of laser cooling techniques for atomic gases, that the properties of He 4 superfluid are related to BEC. This was not obvious given that the (textbook) description of BEC as an ideal non-interacting Bose gas, contrary to He 4 that is a strongly interacting fluid. This gave relevance to the, until then, only theoretical ideas of Einstein, and BEC became a rich topic of research. The relation between superfluids and BECs was confirmed years latter in ultracold atomic gases where almost the entire fluid at low temperatures is condensed and exhibits superfluidity.
It is very challenging to describe the strongly interacting helium system. A weakly interacting Bose gas was then proposed by Bogoliubov, as a modification of the non-interacting Bose gas model, in order to study the Bose-Einstein condensation and superfluidity. In this way, superfluids can be modelled by a Bose Einstein condensate that has self-interaction, and superfluidity is described as being achieved through interactions in a BEC. Notice that BEC can happen even in the absence of self-interaction, as seen above, since it is a statistical property of a gas of bosons in low temperature, but this system does not exhibit superfluidity. The weakly interacting theory is used to describe many superfluid systems at certain limits. This description, tough, evolved in the last few years in order to extend and generalize this framework to finite temperature systems, mixtures and even stronger interacting system corrections. New frameworks also emerged in order to describe different systems that cannot be modelled by the weakly-interacting theory. One of those ideas based on the hydrodynamical description is to write these systems as an effective field theory (EFT) in order to describe the system macroscopically using symmetry alone without the need of working its microscopic description. This EFT, depending on the symmetries of the system, can recover at some limits the weakly interacting superfluids, but also can be used to describe more general superfluids, superconductors and even systems like the unitary Fermi gas [106,107], which is a gas of fermions that interact through a strong 2-body coupling that is a superfluid in the ground state.
The theoretical description of those condensed matter systems, together with the experimental efforts, is a field of research that is in fast development. In this review we are going to describe the basic concepts on BEC and superfluid, and detail the different descriptions and properties of these systems. We start by describing the non-interacting ideal gas, where condensation was first conceptualized in order to present in more detail the definition and the conditions for condensation. We then start to describe superfluidity. We show first the definition of superfluidity as defined in Landau's theory of superfluidity. We then describe a more concrete model for a BEC where superfluidity is present, the weakly interacting Bose gas. This is the simplest example of superfluidity. We show how this model describes condensation and superfluidity. We then follow to show the field theory description of the superfluid, where the system is described as a system which undergoes spontaneous symmetry breaking caused by the condensate. This description brings advantages and makes clear the study of many features in BECs and superfluids. As a low energy description of the superfluid, we present the EFT of superfluids as another description of more general superfluid systems. We finalize describing what happens when we rotate a superfluid, showing the nucleation of vortices upon rotation.
Not linked to what we will discuss in the review, but an important fact. Nowadays, it is known that superfluidity is not necessarily linked to condensation. Recent investigations seems to point that there are states where you can have superfluidity for the majority or all the particles in the system, while only a small fraction is condensed. This happens for example for liquid helium below a certain temperature.

Non-interacting ideal gas
We start our discussion with the non-interacting Bose gas. The properties of this system are a consequence purely from the quantum statistic of indistinguishable bosons. We will see that in the grand canonical ensemble we can write the Bose distribution function in which we can see the conditions for Bose Einstein condensation.
We want to describe here a theoretical gas of many non-interacting bosons in a box. In a system with a large number of particles (N ), it is impractical to try to determine the state of each particle or even the collective many-body wavefunction that describes this system Ψ (r 1 , · · · , r N , t). In this sense, to describe a system with many particles that can occupy many different states, we represent the system using a statistical ensemble description. To describe the state of this collective system, one does not need to label the state of each particle, but to determine the number of particles in each state of the system. The ensemble that is specially convenient for this task of deriving the probability of microscopic states is the grand canonical ensemble. Since our system is composed by bosons, which are indistinguishable particles Ψ (r 1 , r 2 , · · · , r N , t) = Ψ (r 2 , r 1 , · · · , r N , t), called Bose symmetry, this ensemble is useful to describe the system where many particles can occupy the same state.
The GCE is a statistical ensemble that describes a system that is in contact in thermal and chemical equilibrium with a large reservoir, in a way that there is exchange of energy and particles with the reservoir. This exchange of particles with the reservoir makes the number of particles in the system to fluctuate, although the number of particle of the system plus reservoir is constant. As the system is in equilibrium the energy and particle number fluctuate around an average. This ensemble can be described by the following constants: the chemical potential (µ) and the temperature (T ), which hold for a system (of volume V ). If the GCE is applied to small systems, an additional condition is necessary: that the gas is diluted. In principle, the probability of finding the system in a state s with energy s and n s particles, or occupation number n s , is given by: where β = 1/(k B T ), with k B the Boltzmann constant. The chemical potential µ = (∂E/∂N ) S,V is the energy required to add one particle to the isolated system, fully determined by N , the total number of particles, and T . The chemical potential is defined to be negative (so no unphysical negative occupation occurs). The total energy of the system is given by E = s n s s . The normalization Z GC s is the grand canonical partition function. With the GC distribution function, we can then evaluate the average occupation number, where the sum converges for µ < . This is the Bose-Einstein distribution. This gives us the total number of particles in the system: We can separate the total number of particles into two contributions, where N 0 = 1/e β( 0−µ) − 1 is the number of particles with s = 0, which is the number of particles in the condensate, with 0 indicating the lowest energy of the single particle state. The number of particles that are not in the ground state, not in the condensate, also called the thermal component of the gas, N T = i =0 n i . We can replace the sum for an integral and, from the partition function, the thermal component is given by, where λ T = 2π 2 β. For a fixed temperature N T reaches a maximum when µ = 0 . So, N T is limited, meaning that in this limit there a finite number of particles not in the ground state. At this same point in this limit, N 0 can diverge showing that the number of particles in the ground state grows becoming macroscopically occupied.
This macroscopic occupation of the ground state is seen as a condensation and this phenomenon is called Bose Einstein condensation.
The critical temperature T c defines the temperature below which there is the formation of the BEC. We can define it as the temperature above which all the particles of the system are not going to be in the condensate: N T (T c , µ = 0 ) = N . The chemical potential can be zero at T c , and from (14), for the maximum µ = 0 = 0, we can get that where n is the total number density and ζ(3/2) ≈ 2.612 is the Riemann's zeta function. With that, for T < T c , we expect that most of the particles are going to be in the condensate, and the number of particles in the condensate is: Complementary, the number of particle in the thermal component is N T = N (T /T c ) 3/2 . From that expression we can see that the occupation becomes macroscopic towards small temperatures, as we can see in Figure 5. This indicates the formation of a BEC. This condition that a BEC can form for T < T c can be translated into the condition: nλ 3 dB 1, where λ dB = 2π /(mk B T ) is the thermal de Broglie wavelength that gives the coherence length of the gas. This condition indicates that the gas needs to be dilute in order for condensation to happen. This condition is also equivalent to having the de Broglie wavelength of the particles overlap and the system being described by a macroscopic wave-function. When T = 0, all the particles of the system will be in the ground state and the condensate is described by a single macroscopic wavefunction. As we see from Figure 5, at high temperatures, condensation is broken and the system behaves as a gas of individual massive particles.
de Broglie wavelength: associated wavelength of a massive particle given by λ dB = h/p = h/mv, where p and v are the momentum and velocity of the particle with mass m, respectively. For an ideal gas of temperature T in a volume V , we have the thermal de Broglie wavelength which determines the coherent length of the gas. The thermal de Broglie wavelength can be defined then, where the characteristic thermal momentum is When the thermal de Broglie wavelength is much smaller than the interparticle distance (d), the gas is classical. As the temperature becomes smaller, the de Broglie wavelength becomes comparable to the interparticle distance, the quantum nature of the particles becomes important, until when λ dB > d we have the formation of a condensate.
In summary, Bose Einstein condensation can happen for ideal gases. The condition for condensation is that the occupation number of the ground state is so large that becomes macroscopic, this happens when T < T c there is formation of BEC; and for T > T c there is no condensation. This can be translated in a condition for condensation: if nλ 3 dB 1, there is the formation of a BEC. For such large occupation, the quantum corrections are suppressed, and the BEC can be written as a good approximation in terms of classical theory. With that we can see that this very simple theoretical model already shows this intriguing macroscopic quantum phenomenon that would be confirmed experimentally many years later.

Condition for condensation of a non-interacting ideal gas
The condition for condensation of an ideal gas N bosons in thermal equilibrium with volume V and temperature T is: These conditions state that the temperature must be smaller than the critical temperature; or that we have a macroscopic occupation number of the ground state N 0 ; or that the de Broglie wavelength needs to be bigger than the mean space between particles in order to have quantum degeneracy. It is easy to see that these conditions are equivalent. In this figure we plot the number of particles in the ground state, normalized by the total number of particles, with respect to the temperature for the non-interacting Bose gas. We show schematically that for higher temperatures T ≥ Tc the system is in the normal state where the particles behave as free particles and occupy all energy levels. As the temperature is lowered, when T < Tc we have the formation of a condensate described by a macroscopic wavefunction. When T = 0, all particles of the system are in the ground state and we have a pure BEC, described by a single wavefunction.

Landau's superfluid model and criteria for superfluidity
A few years after BEC was theorized, another striking macroscopic quantum phenomena was observed, superfluidity. Landau constructed a phenomenological theory to explain the results of superfluidity in helium, which was observed to flow in thin capillaries. This phenomenological theory, however, is quite general to describe superfluids and gives general conditions for the appearance of superfluidity. This theory has the goal of explaining why in superfluids charge is transported without friction. As we described above, according to London's ideas, in order to have superfluidity one needs to have a BEC. The condensate has the role of transporting charge. So we consider a superfluid as the condensate that transports charge without losing energy. Dissipation of the condensate, which is equivalent to friction in the fluid, is caused by exciting particles out of the condensate. We have a supeffluid in the limit of no or low dissipation, and the superfluid is lost in the limit of high dissipation. We present now the conditions for that to happen.
Consider a superfluid moving through a capillary with velocity v s . The energy of elementary excitations is given by [111], in the rest frame of the capillary. The kinetic energy of the fluid is given by E kin , and p > 0 and p are the energy of the excitation and momentum in the frame of the fluid, and translated to the frame of the capillary. Dissipation happens when p + p · v s < 0. This can only be negative if its minimum, when p + pv s cos(θ) when θ = nπ for n integers, is smaller than zero: p − p v s < 0. With that we can determine the critical velocity: For v < v c , with v c = 0, the system transports charge without dissipation and the coherence of the BEC is maintained. This is the first criteria for superfluidity. The second necessary criteria is that v c cannot be zero, so we need to have a condensate that transports the charge. A non-interacting (pure) Bose gas like we saw in the previous section has v c = 0, so it cannot be a superfluid. A weakly interacting Bose gas has v c = 0 and it is a good representation of a superfluid. As we are going to see in the next section in the case of the weakly interacting BEC that because of the spontaneous breaking of the U(1) symmetry, a Goldstone mode appears, the phonon. This mode is gapless p=0 = 0 and it is an elementary excitation of the superfluid. Even for that mode, the critical velocity is not zero, so there is some cost for producing the gapless excitation. For this weakly interacting Bose gas, the critical velocity is the fluid sound speed, and Landau's criteria for superfluidity becomes: In summary, given that v c is nonzero, and that we have a condensate (by construction) in this system, if v < v c , we can say that there is superfluidity. This results, however, is only valid for zero temperature.
Landau also developed the theory for a superfluid at finite temperature, the two fluid model. At finite temperatures the fluid has two components: the superfluid component that flows without friction, and a normal fluid which describes the excitations. In this theory then there are two sounds speeds, for each degree of freedom. In the case of weakly interacting Bose gas, the first sound is c s associated with the oscillation in density, and the second sound is c s / √ 3 that corresponds speed of propagation of a temperature oscillations. This phenomenological theory is still an important topic of research as a condition for superfluidity. From simulations to experiments it is interesting to ask if the Landau criteria is fulfilled as a criteria for supefluidity. This criteria seems to be valid only in the regime of linear perturbations. This is the case since in Landau's theory the superfluid dissipates only into elementary excitations. However, we know that it is also possible to exist quantum vortices, topological defects present in rotating superfluids. This phenomenological theory does not take that into account, which can change the critical velocity of the superfluid to smaller values, reaching the dissipative regime of the superfluid before than if using Landau's critical velocity [110].

Landau criteria for superflduidity
Phenomenological conditions for superfluidity (at T = 0): 1. Existence of a condensate; 2. v c = 0 (Non-interacting Bose gas has v c = 0 -not a superfluid! Interaction is crucial for superfluidity.) 3. v < v c -system transports charge without dissipation and the coherence of the BEC is maintained. v c =velocity above which excitations can leave the condensate (v c = c s -interacting BEC) At finite temperatures, two fluid model: the superfluid component that flows without friction, and a normal fluid which describes the excitations.

Weakly interacting Bose gas -superfluid
We now turn to the discussion of interacting systems. Inspired by Landau's phenomenological theory, the weakly interacting Bose gas was proposed as the simplest system to study superfluidity, and as a realistic model to understand condensation. In this section we are going to model a superfluid by a Bose Einstein condensate that has self-interaction, and show that, although a BEC can be formed both in the case of the non-interacting and interacting Bose gas, the presence of interaction is crucial for superfluidity [111,112,113].
We present here a microscopic description of superfluid which arises upon condensation. The microscopic system where this happens is a weakly-interacting gas of bosonic particles. To describe this interacting gas, first we need to understand how to describe the excitations in this system.
One of the conditions for condensation is that the gas is dilute. However, even in a dilute gas the interaction can play an important role. The way we describe interactions or collisions in a Bose gas is somewhat different than in a classical fluid. Since now we describe it using their wavefunction, we need to have a interatomic potential V int to enable these collisions. In a dilute system at low temperatures, three-body interactions are suppressed, so we are going to describe this system with binary collisions. In such a system, the two-body collisions depend only on one parameter a, the s-wave (or coherent) scattering length [114], which is the zero energy limit of the scattering amplitude a = lim T →0 f scat . This is valid only for low energies when the other length scales of the problem λ dB , d a. In this limit, the elastic scattering cross-section becomes constant σ = 8πa 2 , and the two-body interatomic potential can be written as V int (r − r ) = (4πa 2 /m) δ(r − r ) ≡ g δ(r − r ), which is short-ranged and present only when the atoms interact. The s-wave scattering length a can be positive or negative depending on the system described, representing a repulsive or an attractive interaction. With this we can model the effective interaction Hamiltonian asĤ The dynamics of this many-body interacting system is given by the second-quantized N -body Hamiltonian, where the brackets is the single particle Hamiltonian, and V trap (r) is the trapping potential, an external potential applied to the system (that in the next section could be the gravitational potential). In the Heisenberg description, we can then write the Heisenberg equations of motion, with the brackets indicating the commutator. This is the Schrdinger equation for the Bose field operatorΨ (r, t). The Bose field operatorΨ † (r) andΨ (r) create and annihilate a particle at position r, and obeys the canonical commutation relations with only non-zero commutator given by [Ψ (r),Ψ † (r)] = δ(r − r ). The Bose field operator describes a continuum spectrum of single particle position eigenstates, and can be re-written in the single-particle basis asΨ where φ i is the states wavefunction, and the creation and annihilation operators,â † i andâ i , create and annihilate a particle from the state φ i . They obey the Bose commutation relations 12 , with only non-zero component given by [â i ,â † j ] = δ ij . Many-particle systems described by the Hamiltonian (20) are very difficult to be solved. With the exception of a few simple models, in order to find solutions to this problem and be able to study its properties we need to make simplifications. For that we use Bogoliubov's prescription or mean-field approximation. For the general case of a non-uniform gas, the mean-field approximation can be written, in the Heisenberg picture as, where ψ(r, t) ≡ Ψ (r, t) is classical field called the wavefunction of the condensate. The density of the condensate is fixed by: n 0 = |ψ(r, t)| 2 = n. Like we described for Landau's theory, δΨ (r, t) is a small perturbation of the system with δΨ (r, t) = 0 and describes depletion of the condensate. Effectively this approximation leads the many-body problem to be reduced to a single body problem by describing the averaging the effects of all other particles. Given that the interactions are weak, and that the gas is diluted, quantum fluctuations on the condensate are suppressed. The mean field approximation is valid for dilute system with n a 3 1. When this conditions is not met there are deviations of the mean-field approximation. We can treat these deviations in perturbation theory, where we invoke non-vanishing moments for the fluctuation operator, like the Hartree-Fock-Bogoluibov, which considers a non-zero δΨ , or the Hartree-Fock-Bogoluibov-Popov, for terms up to second order in the perturbation.
With this approximation, we can write the generalized Gross-Pitaevskii (GP) equation: The GP equation is a non-linear Schrdinger equation, with non-linearity arising from the third term on the right hand side of the equation. The GP equation describes the dynamics of the zero-temperature dilute weakly-interacting Bose system by allowing us to determine the shape of the single particle wave function, the condensate.
We can study the case of stationary solutions. The stationary solution can be taken as the solution that provides us with the condensate, the ground state wavefunction. The ground state is the lowest energy state of a quantum mechanical system, with the excited states being the states with higher energy than the ground state. The stationary states Ψ are the eigenfunctions of the Hamiltonian operator, with eigenvalues µ related to the energy of the system, in a way that for the wavefunction we have With that we can write the stationary solution as, 12 In terms of the creation and annihilation operator, the Hamiltonian of the many-body system is given by, where ij|V |km denotes the matrix element for the interaction, and H sp where the eigenvalue of the Hamiltonian µ is also called the chemical potential, and φ s is real field with dr ψ 2 = N 0 = N . The Gross-Pitaevskii equation becomes: In the Thomas-Fermi limit, which is the approximation where the interaction energy is bigger than the kinetic energy for a large number of particles, the kinetic energy can be neglected so we have µψ = (g n + V trap )ψ. As a solution of this equation, we get that in the Thomas-Fermi limit,

Fluid description
We can decompose the complex macroscopic condensate wavefunction into, where, as we saw above , ψ is normalized to the total number of particles, |ψ(r, t)| = n(r, t) = ρ(r, t)/m, and θ(r, t) is the phase distribution. Inserting this into the GP equation we get two equations. We make the following redefinition v(r, t) ≡ m ∇θ(r, t) .
They are a representation of the GP equation into "hydrodynamical" equations since they have a similar form as the continuity equation and Euler equation for a perfect fluid. However, the second Madelung equation describes a fluid with a potential flow, given the definition of the velocity, with zero vorticity ∇ × v = 0. This corresponds to the main characteristic of the superfluid that it flows without friction, has irrotational flow. This equation also differs from the perfect fluid Euler equation by the presence of the quantum pressure term. The second Madelung equation, the Euler-like equation, reveals more interesting properties of the superfluid. This equation has two pressure terms, P int and P QP that are respectively the pressure term, and the quantum pressure: The pressure term comes from the self-interaction which gives a polytropic type of pressure. For the two-body interaction, which is the case we show here, j = 1. If we have a three-body interaction, for example, the pressure would have polytropic index j = 1/2, giving P int ∝ ρ 3 . The constant K depends on the interaction constant. The quantum pressure is defined in terms of the quantum potential Q = −( 2 /2m)∇ 2 √ n/ √ n. The quantum pressure comes from the quantum mechanical kinetic term that arises due to the uncertainty principle (see more details below), and is characteristic of this quantum system.
Quantum pressure: Quantum pressure (QP) is the name given to the term 13 where P ij, QP is the quantum pressure tensor. Together with the second term on the left hand side of equation (33) this term comes from the spatial part of the kinetic term. However, those two components are very different. The classical component describes the kinetic pressure due to the motion of the particle.
The quantum pressure comes from the quantum part of the kinetic term that arises due to the Heisenberg uncertainty principle. This is an additional force term that appear in the Mandelung equation due to the zero point motion of particles. This term is repulsive and counter-acts attraction from a potential or attractive interaction, supporting the system against collapse. With that, the system cannot have vanishing size. This term modifies the dispersion relations of the excitations of the condensate and it is important at small scales, for scales smaller than the healing length. This term is negligible for large scales, the Thomas-Fermi approximation, and for a uniform superfluid, since n = const..
Healing length: defined as the length for which P int = P PQ , given by, It is the length for which the interactions "heal" (coarse-grain) any density of phase perturbations in the condensate.

Condensate solution
Having established how we describe the weakly interacting system in the mean-field approximation, we can now describe what is the ground state, the condensate. We can solve analytically the GP equation for a few simple cases and obtain the condensate solution as for the cases with no interactions, in harmonic potential, and other simple systems [112]. We are going to present some interesting cases here.
In the case of a uniform gas, the wavefunction Φ 0 (r) = √ N φ 0 / √ V describes the condensate formed, while the remaining functions form a complete set of functions orthogonal to the condensate.
Solitons -Solitons are a localized solution of the one-dimensional GP equation with V trap = 0, which is integrable in this limit. They are also called solitary waves, and are solutions described by a permanent and localized wave, that maintains its shape and velocity upon collisions. They are obtained when the dispersion term, the kinetic spatial term, and the non-linear term from the interaction cancel out. If the interaction is repulsive (g > 0), a dark soliton is formed,which is given by the solution ψ d (x) = ψ 0 tanh(x/ √ 2 ξ), where ψ 0 = lim x→∞ ψ. A bright soliton is the solution for attractive interactions (g < 0), with ψ b (x) = ψ 0 e −iµt/ cosh( 2m|µ|/ x) −1 . They are called dark and bright since they represent a decrease and a concentration of condensate, respectively.
In reality, soliton-like objects are more interesting representing solutions in 3-dimensional condensates in a trapping potential of near to or non-integrable systems.

Collective excitations, dispersion relation and sound speed
The excitations in the superfluid which represent perturbations of the condensate, are an important part of this system. They represent sound waves, called the phonons that propagate through the condensate. We are going to show how they arise.
Here we are going to work in the case V trap = 0, for simplicity. The case where a trapping potential is present, like for example a condensate in the presence of a gravitational potential, is studied in Section 4.1. For a homogeneous condensate, to study the perturbations around the condensate, we perturb the classical wavefunction ψ(r, t), where we assumed that the motion is only in the x-directions without lost of generalization. The perturbations are small and δψ (i) indicate the i th order in perturbation. To linear order, we can re-write the linearized GP equation: ). We make an ansatz for the solution as travelling waves, δψ (1) The parameters A, B are determined by the initial conditions and the dispersion relation is given by, 13 One can notice that the form of the quantum potential coming from the Mandelung equations, and consequently the quantum pressure tensor, are very similar to the Bohm quantum potential [116,117]. Some authors point that Bohm rediscovered the quantum pressure in his new interpretation of quantum mechanics. However, some authors claim this equivalence is not so clear. Some authors also claim that since the quantum pressure tensor has non-diagonal components, and it is not isotropic, so it cannot be called pressure [118].
where the sound speed is defined as the term coming from the linear part of the dispersion relation, The sound speed appears because of the presence of the interaction, therefore since for a superfluid the presence of an interaction is crucial, we can say that a superfluid is a fluid that has a sound speed. With this we can easily see the definition and properties of a superfluid.

Superfluid:
In the presence of interactions, a sound speed is present which determines the behaviour of the excitations on large scales.
-For large wavelengths (small k), the higher k terms do not contribute and the dispersion relation is given by: which is the dispersion relation of a sound wave. The superfluid is characterized by excitations that propagate as waves, the phonons. Because of this property phonons can mediate a long range forces (F ∼ 1/r 2 ). This longrange force is the responsible for the effective dynamics of a superfluid: flowing without friction. We are going to see later (in the field theory description) that this gapless mode can also be viewed as the Nambu-Goldstone modes coming from the spontaneous symmetry breaking of the U(1) symmetry of the system caused by the formation of the BEC. This limit where the quantum pressure term can be ignored, or more specifically when P int P QP , is the Thomas-Fermi approximation we saw above, and can alternatively be defined by wavelengths bigger than the healing length. This is the case for a repulsive interaction (g > 0). The situation is different for an attractive interaction (g < 0), where w k is imaginary and the solutions is unstable given by exponentially growing or decaying functions. This means that it is not possible to form a stable condensate in these cases.
-For small wavelengths (large k), the quantum pressure term dominates, and the dispersion relation is given by ω k = k 2 /2m, which describes a free particle. In this limit, the system stops exhibiting superfluidity. In general, for intermediary frequencies, the full dispersion relation (38) does not propagate as a wave and shows two degrees of freedom: a gapless mode (the phonon), that propagates as a wave, and a massive mode, related to particle creation. We are going to discuss these again in the field theory section.
No interactions -BEC: In the limit where g → 0, the BEC stops exhibiting superfluidity. The phonon becomes gapped, the dispersion relation is given by: ω k = k 2 /2m, which is the dispersion relation of a free particle. Since ω 2 k > 0 in this case, the solution of the linearized GP equation (without a external potential) is a stable oscillatory solution.
We showed above that the superfluid can decay into collective elementary excitation, the long-wavelength soundwave quanta, the phonon. These are excitations with linear dispersion relation that behave as periodic fluctuations in density in the superfluid. When studying linear perturbations around a classical condensate background, only the phonon excitation is expected. However, a superfluid can also deplete into other excitations called rotons and maxons. The linear part of the dispersion relation is only a part of the dispersion relation at small momenta. For higher momenta, the dispersion relation presents a maximum and then a minimum. Near the maximum we have the maxons excitations, and near the minimum, for even higher momenta, we have the roton excitations. So, phonons, rotons and maxons are excitations described by different parts of the same dispersion relation. And lastly, for a rotating superfluid, there is also the vortex nucleation of the condensate, as we will see in Section 3.5.
U(1) symmetry group: The unitary group of degree n = 1 is the group associated with n × n unitary matrices (U * U = U U * = 1) under matrix multiplication, U = e iα , with α the parameter of the group. For n = 1, this is the unitary transformations complex numbers, which corresponds to phase rotations. The U (1) group is isomorphic 14 to the SO(2), the group of the 2 × 2 orthogonal rotations in R 2 with unit determinant. The U(1) symmetry can be global, when it acts in the same way in all points in space-time, or local, acts differently at each place in space and time.
Noether's theorem: One of the most important results theoretical physics, Noether's theorem definition can be found in basically every textbook in field theory and even classical mechanics. The theorem states that for each continuous symmetry of a action there is an associated conserved current j ν (x, t), ∂ ν j ν (x, t) = 0 . This has profound implications since the conservation laws that describe the dynamics of a systems are then associated with the symmetries present in the problem. We can illustrate this for a system with the global U(1) symmetry that, according to Noether's theorem has a conserved particle number density 15 . Consider a complex scalar field (Ψ (x, t)) theory with action: where V (x) is a potential. On top of being invariant under Lorentz transformations, this system has U(1) symmetry, with the action invariant under continuous rotations of the phase of the complex scalar field: Ψ → e iα Ψ (and c.c.). Since this is a global symmetry, α is independent of (x, t). With that, we can evaluate the Noether current and charge of this system: which are conserved: ∂ ν j ν = 0 and dQ/dt = 0. From j ν = (mn, j), the conserved charge can be written in terms of the local density of Q, Q = d 3 x ρ(x, t), and the conserved current implies a continuity equation for ρ. In this way, Q ∝ Ψ † Ψ , which implies the conservation of the norm and it is associated with number conservation.
Some comments are in order. When defining ψ as the condensate wavefunction, this quantity is actually a meanfield value of the wavefunction, the degree of freedom that defines the condensate. This description of averaging selecting the condensate is consistent with the theory of critical phenomena, like phase transitions. This Bose system can be seen as a system with spontaneous breaking of a symmetry of the description. In our case the U(1) symmetry which is the symmetry of the Hamiltonian. This is analogous to the the spontaneous symmetry breaking in a ferromagnet. The difference is that, since we have a Bose system, the idea of spontaneous symmetry breaking to the thermodynamical limit of a finite size Bose gas defines the number of particles of the system, and this is only consistent with the picture of having a condensate that can change the number of particles in the ground state, if the number of particles is conserved.
We can also understand this argument for the symmetry breaking of the system by analysing another approach to the weakly interacting system initially developed by Bogoliubov, Penrose and Onsager [121], and Beliaev [120]. In this approach the condensate wavefunction is identified, using the density matrixρ(r, r ) = Ψ (r )Ψ (r) , to the number density of particles n(r) =ρ(r, r) = i n i |φ(r)| 2 , where we are working in the stationary case. The formation of a BEC, which means that the ground state has a macroscopic occupation, leads to the factorization The field operator then can be factorized in the presence of a condensate intoΨ (r) = Φ 0 (r) + δΦ(r), with δΨ (r) = 0. This means when a condensate form, A state that has conserved particle number has Ψ = 0. So this condition above is seen as describing a symmetry breaking, more specifically Bose symmetry breaking, and the consequence is a system where particle number is not conserved in the condensate. In the absence of a condensate this goes back to the particle conserving condition: This suggests that condensation comes from a spontaneous symmetry breaking theory.
With that, the interacting condensate system above can be understood as a particle conserving system of bosons with U(1) symmetry, described by the classical field ψ, the wavefunction of the condensate, where the formation of a Bose Einstein condensate is a phase transition, coming from a spontaneous breaking of the symmetry (that can be seen as spontaneous coherence). This description of the condensate makes us see a parallel with the formalism used in field theory. 15 Isomorphism is a one on one mapping between the elements of groups preserving its group operations. In the case of U (1) and SO(2), the isomorphism takes a complex number of unit norm into a rotation, e iθ → cos θ − sin θ sin θ cos θ , where the rotation angle is the argument of z. 15 For a more general derivation of Noether's theorem in field theory for any continuous symmetry and the connection with the generators of the symmetry group see [?]

Field theory description
Given the suggestion that we can describe condensation as spontaneous symmetry breaking process, we turn now to the description of this system using a field theory language. The methods from field theory are very appropriate to describe this type of system where spontaneous symmetry breaking is present. Given the description we had above for the superfluid and the properties of this system, we represent this system as a massive complex scalar field 16 with self-interactions with a global U (1) symmetry that is spontaneously broken by the presence of the ground state, the condensate, with superfluidity arising upon condensation. We are going to work here, in order to best illustrate this description, in the homogeneous case, where there is no trapping or external potential applied to the system. Given that, we describe this system by the Lagrangian density for a two-body interaction, where we consider g > 0 in order to get a stable condensate with long range coherence, as discussed above. As we saw above, this system has a U (1) global symmetry, which means that it is invariant under continuous rotations of the phase of the complex scalar field, where α is a constant 17 . The equation of motion is given by, In order to describe the condensate, like we did previously, we separate the condensate contribution to the perturbations, and assume the mean-field approximation, where ψ(x, t) is the background that gives the condensate field, and δΨ (x, t) are the excitations of the condensate that are considered small with respect to the background solution.
As we saw in (27), the stationary solution determines the ground state of the system. Using that, the condensate wavefunction can be written as where θ bg = µt is the phase of the ground state where µ is the chemical potential 18 . This is the background solution for the the equation of motion (47) as long as, When we have a spontaneous symmetry breaking, only the ground state, the state that broke the symmetry, is not invariant under U (1) anymore. This means that we still have a conserved current given by , as we saw above in the definition of SSB. It is easy to see using the equations of motion that this current is conserved. For the ground state (49), this current gives the number density: n = 2µv 2 . In the field theory description the conservation of the norm ( |Ψ | 2 ), related to the number density in the condensate and conservation of the number of particles in the condensate, seen in the quantum mechanical approach above is interpreted as a consequence of the global symmetry of the system which leads to a conserved Noether charge, Q = j 0 ∝ n ∝ |Ψ | 2 , conservation of the number density of field quanta. 16 We have used here a complex scalar field that has a continuous symmetry since, as we saw before according to Noether's theorem, this system has a conserved charge. The U (1) symmetry gives us a system with conserved particle number. We can see that using a complex scalar field is essential, since a real scalar field has instead a discrete Z 2 symmetry, which implies that there is no conserved charge. And this would not be a case that reproduces the known behaviour of a superfluid. 17 It is equivalent to re-write this, in a language of SO(2) symmetry, where the complex field can be written as where the system as an invariance under rotations: Potential of the weakly self-interacting system. Left panel: Potential of the symmetry restored phase, m 2 > 0 which has a minimum at ψs = 0. This ground state is symmetric, respects the symmetry of the system. Right panel: Symmetry breaking potential, for m 2 < 0. This phase has a degenerate minima not invariant under the symmetry of the system. It represents the condensate state.
Spontaneous symmetry breaking (SSB): SSB is the phenomenon where the state of the theory is not invariant (not symmetric) under the symmetry transformations (U ) of the Hamiltonian (or action) that describes the system. This stable state |ψ spontaneously broke the symmetry of the system. This mechanism offers an explanation for why there exist stable states, like a condensate, that do not respect the symmetries of a system. This allow for the existence of different symmetry related state U |ψ , of same energy but different phases defining a set of symmetry-broken states. To distinguish them, or their phases, we have the order parameter of the system, defined The order parameter can be used to identify if a symmetry was broken O = 0, where the system is said to have long range order, given by its two-point function (C(x) being proportional to a constant. For unbroken This is the background solution corresponding to the condensate. This ground state is responsible for spontaneously breaking the U (1) global symmetry. And we can see that explicitly. Since this solution represents the ground state, we can compute the energy functional for this system, where V ef f = m 2 Ψ * Ψ + gΨ * Ψ . The ground state is given by the stationary, minimum energy state. This can be found by finding the minimum of the energy (51), which amounts to finding the minimum of the potential energy. The set of solutions for the minimum of this potential energy is: The value of the minima v 0, ssb = ±m/ √ g are called vacuum expectation value and are the value that the field assumes at the ground sate, apart from a phase. As we can see from on the left panel of Figure 6, the symmetry restored phase when m 2 > 0 has a well defined minimum at v 0, s = 0, and this is the normal phase, with no condensation. This ground state ψ s = 0 is preserved under the symmetry of the system (U (1) symmetry), the rotations of the phase of the complex field. The symmetry breaking phase when m 2 < 0, there is a continuous set of minima with the ground state given by ψ ssb = v 0, ssb e iα , corresponding to all the possible phases in the circle α ∈ [0, 2π), as seen in the right panel of Figure 6. All of these classical backgrounds are not invariant under U (1) symmetry, the symmetry of the system, which means that φ i, ssb → φ j, ssb = v 0, ssb e iα = φ i, ssb . In this way we say that the symmetry is spontaneously broken by this condensate ground state. From the SSB we can see that the condensate has long-range order, with the field having the role of the order parameter of the system. 19 For a very complete and extensive review on SSB, see reference [123]. 19 An alternative way of writing this Lagrangian that leave explicitly the spontaneous breaking of the symmetry by the finite charge:

Excitations
We consider now fluctuations around the classical condensate configuration in order to study the spectrum of this system. Considering small fluctuations means perturbing each degree of freedom of the field around the condensate (48). This is equivalent, in the polar notation to: where ρ(x, t) can be interpreted as a perturbation in the radial direction, and π(x, t) a perturbation in the angular/phase direction. Plugging this into the Lagrangian, we have, With this expansion is already easy to see that ρ has a mass term (the term accompanying the ρ 2 term), so the perturbation in the radial direction is massive. The perturbation in the phase has no mass term and it is going to be massless. This massless excitation was already expected from the Goldstone theorem (see below), where a SSB leads to the appearance of a massless excitation, the Nambu-Goldstone (NG) boson π. In the context of a superfluid this gapless excitation is the phonon. For low energy theories, the phonon is the only degree of freedom that is excited in the theory, as the massive mode can be integrated out. Therefore, at low energies a superfluid is completely described by the phonon excitations.
Low energies here mean energies lower than the mass gap of the massive mode ρ. We are going to work on this limit here to obtain the spectrum of this NG boson. In this limit, the kinetic term of ρ can be neglected and the equation of motion becomes: Using this we can rewrite the on-shell Lagrangian in terms of X as, This Lagrangian depends only on the phonon π, as ρ was integrated out, and is invariant under shift symmetry, π → π + c, inherited from the U (1) symmetry of the complex scalar field. This is the effective Lagrangian at leading order in derivative expansion for the phonon. To obtain the dispersion relation we can expand this Lagrangian, using π c = µπ, where we used that n = 2µv 2 . Taking the Fourier transform of the field, we compute the equation of motion and can see that the dispersion relation of the phonon is given by, where c s is the sound speed of the superfluid. This expression shows us the general behaviour we already seen in the excitations calculated in the QM approach. For long wavelengths (small k), the higher order term in k is suppressed, and the dispersion relation is w 2 k, long = c s k 2 , which is the dispersion relation characteristic of a sound wave. Therefore, the phonon propagates as a wave with sound speed c s in the long wavelength regime.
The second term of the dispersion relation is characteristic of a massive particle. This means that for short wavelengths, the system behaves as a system of normal particles propagating, and not a superfluid.
We can re-write these expressions in the non-relativistic (NR) regime, which is the regime we are interested in comparing to the previous approach (and also the limit where we describe, in the next section, the behaviour of DM in galaxies). In this limit, g v 2 m 2 , which implies µ 2 ≈ m 2 . The dispersion relation in this case is given by This is equivalent to the usual way we introduce the chemical potential in the Hamiltonian: H − µN , where N = j 0 is the conserved charge. For this modified Lagrangian, the condensate would have trivial phase θ bg = 0. This is equivalent to what we did in the text where the Lagrangian had the canonical kinetic term, but the phase of the condensate had a time dependency with the chemical potential θ bg = µt.
which is equivalent to the one found in (38).
In the absence of interactions, we recover the ideal Bose gas from section 3.1, with background solution is given by Ψ ideal = ve imt , and dispersion relation of a massive particle w 2 k = (1/4m 2 )k 4 , showing again that although this system condenses in to a BEC, in the absence of interactions there is no superfluidity. It is easy to see from (58) that in the relativistic regime, g v 2 m 2 , c 2 s = 1/3. In this way we showed that the field theory description of the superfluid is very good to describe the general properties of the superfluid. In order to properly compare with the results obtained above in the QM approach, we are going to show that in the non-relativistic limit, the field theory description gives us the GP and Madelung equations obtained above.
Before doing that, one comment is in order. In (57), we showed that in the low energy regime we can re-write the microscopic theory of a superfluid as an effective theory only of the phonon with a non-canonical kinetic term. We did this here in the case of a weakly self interacting system with two body interaction, but in section 3.4 we are going to extend this idea to general superfluid systems and show the construction of the EFT of superfluids.
Goldstone theorem: The Goldstone theorem [124,125,126] or Nambu-Goldstone theorem, states when a spontaneous symmetry breaking occurs a mode with energy that vanishes as k → 0 is present in the spectrum of excitations of the system. This mode is called Nambu-Goldstone (NG) boson and is a massless particle, in the case of relativistic systems or collective excitations with zero energy gap for non-relativistic systems 20 . When a symmetry is spontaneously broken, the Noether theorem still applies so there is still conserved currents. The stable state responsible for the SSB is not invariant under this conserved charge, Q|ψ = 0 (or O = 0). This condition implies that there must be a state, the NG mode, with E k → 0 as k → 0, whose quanta is a massless boson. The NG boson still exists if the symmetry is not exact or broken by an external potential µ, but in this case the mode has a gap µ at k → 0. For ordinary NG bosons (type A), the number of NG bosons created is equal to the number of broken symmetry generators Q, n BG . The Goldstone theorem described here is valid for a system that is invariant under Lorentz transformation, with the appearance of n BG NB bosons. However, SSB is an important phenomena in many systems that are not Lorentz invariant, like the BEC or supefluid, with a number of bosons that will appear in the theory, n NG , called Nambu-Goldstone (NG) bosons, which might not be equal to the number of broken generators, like in the Lorentz invariant case. A generalization of the Goldstone theorem, which includes systems that do not have Lorentz symmetry, exists and can be found in Refs. [127,128]. In these works they classify and generalize the Nambu-Goldstone theorem for any symmetry, including non-relativistic systems invariant under Galilean symmetry, showing how to compute the number of NG bosons created by the breaking of such symmetry.
Recovering the other approaches: We want to show that we can recover the GP theory presented in the previous subsection, and emphasize that the field theory description is compatible to describe the superfluid. The field theory presented above is a fully relativistic theory, which means that the action is invariant under Lorentz transformations on top of the global U(1). However, GP description shown above is non-relativistic. Therefore, in order recover the GP and Mandelung equations we need to take the non-relativistic limit of the relativistic field theory above 21 .
Starting from our field theory for the weakly interacting bosons (45), we take the non-relativistic limit of the Lagrangian. For this we rewrite the fields as, and take the non-relativistic limit. With that, the Lagrangian can be written as: From this non-relativistic Lagrangian we can evaluate the equation of motion for the scalar field ψ, which is exactly the Gross-Piatevskii equation like shown above, in the absence of a trapping/external potential. From that, we can also derive the Mandelung hydrodynamical equations. If we make the following substitution: The vorticity of the superfluid is zero and the momentum density has a non-zero curl. Plugging this in the equations of motion we recover the Madelung equations (32) and (33) in the absence of an external potential. This shows again that we can recover the equations that describe the interacting BEC using the field theory approach.
In this section we showed how to describe a weakly interacting BEC. We showed the standard treatment of the theory, where the many-body quantum system is described by the GP equation. We showed that condensation can be thought as a spontaneous symmetry breaking process and showed that the system can be described in a equivalent way using the field theory approach. We specialized in both cases in the mean-field theory, which is valid for dilute systems, and simplifies the significantly the study of the system.
It is important to comment on the validity of this theory and the approximations made. The theory presented above is only valid for zero temperatures and in the mean-field approximation, that holds for na 3 1, where we can describe the condensate as a classical wavefunction and the limit where quantum corrections are not important. As we cited above, there are correction to the mean field and other approximations where one can study this model (see [111] for some examples). For finite temperature, one has to describe the superfluid as a two-fluid model.
In the cases we are going to study, we will extend a bit the validity of the zero temperature description, as a first approximation, since in galaxies the temperature is obviously not zero. However, since the occupation number will be very high in our problem, the classical description is safe.

Two fluid model
The description presented above for the superfluid is valid for a zero temperature dilute weakly interacting Bose gas. However, as already described in Landau's phenomenological theory, for finite temperature system the superfluid has to be described as a two-fluid model: a mixture of a superfluid and a normal fluid. As we saw previously in this description, for T < T c , the fluid is a mixture of normal and superfluid with most of the fluid in the superfluid phase, while for temperatures higher than T c , the coherence of the superfluid is broken and the fluid is mostly in the normal phase.
Connecting Landau's theory and the microscopic description of a superfluid by a weakly interacting theory, a two-fluid relativistic theory of superfluid was developed. This can be linked to the non-relativistic one described above, and a "hydrodynamical" and field theory approaches are developed. We are not going to describe this in details in this review, since this is not used to describe the current ULDM models in the literature. However, the two-fluid description should be used for finite temperature systems to describe a superfluid and we believe it is important to describe realistic DM superfluid models. For a review of the two-field model, see ref. [122]. In section 4, we present one work where the two-fluid formalism is used to describe the self-interacting BEC DM.
The weakly interacting Bose system studied in this section is the prototypical description of a BEC system that presents superfluidity. It is a microscopical description that shows the behaviour of the condensate and its excitations. But this description is limited and cannot be used to describe all the models of possible superfluids and realistic experimental systems. The interesting point is that this theory can be recast as a spontaneous symmetry breaking theory of the U(1) symmetry of the many-body system. BEC and superfluidity are a consequence solely of the spontaneous symmetry breaking, independent on the specific model chosen. Therefore, this hints us to describing these symmetry breaking systems using the effective field theory approach, where the Lagrangian of the system is described by symmetry alone. The hope is to be able to describe more complicated superfluid systems. We are going to describe this approach then in the next subsection.

Effective field theory of a superfluid
We described in the previous section the construction of a microscopic effective theory for the weakly interacting BEC that can be used to describe superfluids. This description is based on London's idea, that has its roots in the superfluid/superconductor hydrodynamics, that the BEC can be described by a theory with spontaneous symmetry breaking caused be the condensate, with superfluidity arising upon condensation and being described by the Goldstone boson, the phonon, at low energies. We saw that we can write an effective Lagrangian for the phonon that describes the behaviour of the superfluid, matching many observations, in the low energy and momentum regimes.
This procedure shows us that we can describe the hydrodynamical degrees of freedom of a theory by the Goldstone modes created by the SSB of a symmetry, the global U (1) symmetry in our case. This is more general than the simple weakly interacting two-body interaction case showed above. This is already the case of hydrodynamics that describes macroscopically the behaviour of low energy variables and interactions of system given a symmetry, coarse-graining over the smallest scales. This is the perfect playground for the use of effective field theories (EFT), and EFT techniques are very appropriate for this task. An EFT describes the low energy (long distances) behaviour of a system, without having to refer to its underlying microscopic theory, by parametrizing our ignorance of those short scales.
The idea is to use EFT methods in order to describe the dynamics of a superfluid. This was developed in refs. [129,130,131,132], where they develop the general formalism to describe fluid hydrodynamics without dissipation as a EFT. In this approach the Lagrangian that describes the system is constrained by symmetry alone. This approach is very powerful since not only can be applied to many different systems, nut describes the system without the need of its microscopic understanding. It is also very powerful since it is an expansion over momentum. At leading order in the expansion, we describe the low energy theory. But this description allows us to go beyond leading order in the long-wavelength expansion.
In our case, we want to construct the EFT that reproduces, in the long wavelength regime (low energies), the superfluid hydrodynamics, as presented above. This theory is a theory of the phonon, which is the only degree of freedom that is excited at low energies. This is the Goldstone mode produced by the SSB of the U (1) global symmetry by the ground state. We are going to work here in the non-relativistic regime, but one can see the references above for the relativistic case. Restricting to the non-relativistic regime does not imply any loss of generality of the argument, with the system only subjected to different symmetries than the relativistic case.
Inheriting the knowledge of a superfluid from previous sections, at low energies, the only dynamical degree of freedom that describes a superfluid is the phase of the condensate, the phonon. Therefore, in the non-relativistic regime, we need to construct the EFT of this phase π. The superfluid is described by a theory where the symmetry is spontaneously broken by the ground state, The theory that is described by this phase L eff (θ) is invariant under shift-symmetry which is inherited from the U (1) symmetry of the complex scalar field. EFT states that to construct the effective Lagrangian for the phonon, we have to write all the terms that are compatible with the symmetries of the problem. This system has shift symmetry and Galilean symmetry. For the shift symmetry in order for the Lagrangian to be invariant under this symmetry, only terms that are acted by a derivative can appear in the theory: This may contain terms with any power of the derivative of the field. However, this theory has more symmetries. In a generic space-time and adding a gauge field, which is a natural extension of the simple scalar field model, we require that this Lagrangian is invariant with respect to three-dimensional general coordinate transformations and gauge invariance. The most general Lagrangian L = L D t θ, g ij D i D j θ that is invariant under these symmetries, shift symmetry, gauge invariance and general coordinate invariance, is given by, where D t θ =θ + A 0 and D i θ = ∂θ − A i . In flat space, g ij = δ ij , the general coordinate invariance corresponds to a Galilean symmetry.
As it can be seen in [132], Galilean symmetry is not enough to constrain the NLO terms and one needs to consider the full general coordinate invariance. This is equivalent to considering an additional constraint for the theory, which is known from fluid hydrodynamics, where T 0i is the off-diagonal component of the energy-momentum tensor. As stated in [134,131] this additional constraint (with an analogous constraint in the case of relativistic systems [131]) states that only one degree of freedom carries all the current and momentum. Introducing a new symmetry, full general coordinate invariance, is equivalent to assuming this relation.
With the two symmetries of the system, in the absence of gauge fields A 0 = A i = 0, the Lagrangian that describes this system is given by: This is a Lagrangian that has a non-canonical kinetic term and it obeys (69). At the background θ = µt, and T = 0, this Lagrangian density is equal to the pressure P = P (µ). With that we can evaluate the particle number density, where indicates the derivative with respect to A 0 . For the condensate then the equilibrium number density at chemical potential µ is n(µ) = P (µ), where P (µ) is the thermodynamical pressure, defined up to a constant. Given the Lagrangian (70) with θ = µt + φ we can write the Lagrangian as: We can see from that the phonon speed of sound where ρ = mn is the mass density. One limit of this EFT Lagrangian is the quadratic Lagrangian for the weakly interacting BEC, shown in (57). We can see that by considering the special case where P (X) is written as a polynomial, P (X) ∝ (θ/m) n . Depending on the power chosen we will have a superlfuid with a different equation of state: These represent different systems. The case n = 2 is equivalent to the previous case, in section 3.3, where we had a superfluid with a two-body interaction described by the microscopic weakly interacting theory, where we obtained P (X) = X 2 for the low-energy (57). The case with n = 3/2 can correspond to the same theory as the previous case, a weakly interacting theory, but with a three-body interaction, with this effective Lagrangian obtained by integrating out the massive radial mode, like done in the previous section. These two equivalences shows us an interesting aspect of this EFT and from the hydrodynamics of superfluids: the interaction is linked to the equation of state of the superfluid, and this can be seen by a different choice of P (X) in the EFT. The case n = 3/2 can also represent another completely different superfluid system, like we will see in the case of the superfluid DM in the next section, where the this case does not come from a weakly interacting microscopic theory with three-body interaction. This is also the case for the other example we show here, the unitary Fermi gas, which cannot be described by a microscopic theory like we did in the previous section, and being described with this EFT if n = 5/2. One comment is in order here. Usually in quantum field theory having fractional exponents can be problematic, leading to caustics or superluminal propagation. However, in the case of the superfluid this is not a problem. Before reaching these regimes (like the formation of caustics), the superfluid coherence is broken, and the EFT description of the superfluid is no longer valid.
If one wants to add an external or trapping potential V ext , like for example if the gas is under the influence of a gravitational potential, this corresponds to making A 0 = V ext . With that, the Lagrangian is given by (68), with X =θ − (∂iθ) 2 2m − V ext . In the case of the a condensate in a gravitational potential this is given by where Φ is the gravitational potential. This is going to be the case studied in the next section for the DM superfluid. With that Lagrangian we are able to describe the theory as the other approaches we used to describe the BEC and superfluid theories.
Equivalence of the EFT description: For low energies, and in the non-relativistic case, the EFT of superfluids is equivalent to the microscopic description presented above. Considering P (X) ∝ (θ/m) n .
Superfluid hydrodynamics: From this formalism we can also describe the superfluid hydrodynamics. From (71) we re-write the field equation with respect to the number density to obtain the continuity equation: The gradient of the field π can define the velocity of the superfluid, v s = −∇θ/m = ∇π/m, we can derive the second equation of superfluid hydrodynamics:π Validity of the EFT: As we saw above, this theory is valid for low energies, or long-wavelength, and breaks for high energies. The Lagragian shown here is valid in leading order in derivative expansion. In this regime it reproduces the results from hydrodynamics of superfluids. But this framework also allows us to go beyond leading order in this momentum expansion, the next-to-leading order (NLO) Lagrangian. In [132], they show a prescription to take into considerations next-to-leading order terms. This can be done in this framework at arbitrarily order, only requiring that the NLO Lagrangian respects the symmetry of the system, and at the cost of introducing new free parameters. This might allow the study of those systems in a regime that is challenging for the microscopic perturbative description. However, the validity of this effective theory constructed here still needs to be checked as higher order terms in the Lagragian can only be neglected if the sound speed of the theory is not too small. In section 4.2.5 we describe the validity of the EFT for that concrete example of superfluid and show that the theory is valid for the parameters of the model. This theory is also only valid in the absence of dissipation. For a discussion of how to describe dissipative phenomenon in this EFT approach can be seen in [133].
We presented in this subsection the formalism of the EFT to describe the low-energy dynamics of the superfluid Goldstone mode, the phonon. This formalism is more general, tough, and translates into a EFT language the hydrodynamics of fluids at zero temperature and without dissipation, so it can be generalized to describe superfluids with different equation of state, supercondutivity, unitary Fermi gas, among other systems. This low energy EFT approach is very useful to describe the dynamics of various physical systems and writing the superfluid in this macroscopic effective Lagrangian offers us the chance to study the dynamics of this system without having to work out the, many times complicated, details coming from the microscopic short distances physics. With that this formalism allows us to study the behaviour of more complicated superfluids, with different equations of state that might come from different and more complicated interactions. This will be useful to describe the DM superfluid model, in Section 4.

Rotating superfluid -quantum vortices
When we rotate a normal fluid, the fluid rotates together with the recipient in a homogeneous way, like a rigid body, with vorticity ∇ × v n = ∇ × Ω × r = 0, where the normal fluid velocity v n of a rotational fluid is given by the cross product of the angular velocity and the position.
As we saw above, a superfluid, described by a weakly-interacting BEC, has irrotational flow: v s = ( /m)∇θ which gives ∇ × v = 0, where viscosity is absent. This is the defining property of a superfluid. In another language, this means that the circulation, around a closed contour C in a superfluid is given by: where dl is a length element on the path C and A is the area enclosed by this contour. When a superfluid is rotated, this property says that the superfluid would not rotate, but would remain stationary. So, how can we rotate a superfluid and maintain the irrotational flow? This is possible if the superfluid phase presents a singularity. A superfluid is a state where the system is described by one macroscopic wavefunction. In the presence of this singularity, this wavefunction is single-valued, ψ(θ) = ρ/m e iθ = ψ(θ + 2πn), which leads to the quantization of the circulation: This properties above describe a vortex. The way a superfluid rotates is inhomogeneous by forming quantized vortices [135,136,137]. From that we can see that the azimuthal velocity is of a irrotational fluid is given by v φ = n( /m)(1/r), where r is the distance to the center of the closed loop C and φ is the azimuthal angle. At the center we have the vortex core, as lim r→0 ψ → 0, of size equal to the healing length, where the density ρ vanishes and θ, the phase, rotates by 2π around the core. The flow in the center is given by the vortex line. The vorticity of the rotating superfluid is given by: where r i is the location of the N v vortices, and we considered that the vortex lines are in the z-direction. The vorticity is non-zero only at each vortex. So the flow is irrotational in most of the superfluid, except in the vortices. Given that the vortices are line singularities, they form a lattice of uniformly distributed vortices in the superfluid and carry the angular momentum of the rotation L z = n L qm , where L qm = N is the the minimum angular momentum necessary to have one quantized vortex and N = |ψ| 2 is the number of particles in the condensate. This configuration is energetically preferable (instead of for example concentric sheets around the superfluid). With that, the spatially averaged vorticity is given by: where n v = 2Ω/( /m) is the density of vortices, which is related to the angular velocity [138]. Although most of the superfluid has irrotational flow, the rotational flow is given by the vortices in a way that the whole superfluid then effectively flows as a normal fluid v s = v n , allowing the superfluid to rotate. Given that the superfluid has irrotational flow, if it is rotated by a small Ω, the superfluid will remain stationary. There is a critical angular velocity, Ω c , above which the nucleation of vortices in the superfluid occurs. The critical velocity is given by [139]: where R is the radius of the cylinder where the superfluid is contained. Through the dependence in the healing length we can see that the critical angular velocity depends on the interaction, for the case of a weakly interacting BEC that behaves as a superfluid. The case where the interaction is attractive, g < 0, it is known not to produce vortices. This is also the case where the BEC has a finite size coherence, producing only smaller soliton cores. This condition for the formation of cores can also be re-written as L L qm , where L angular momentum of the applied rotation to the system.
If there is the formation of vortices, the area density of vortices is given by n v = n/A BEC = mΩ/π . From that, we can calculate the size of the vortices created (πR 2 v = 1/n v ) which depends on the mass and the angular velocity of the fluid, R v = /mΩ. But what exactly is the structure of these vortices and how are they described? As we saw before, we defined the vortex as a singularity, where the wavefunction is zero at the vortex line, with a quantized circular irrotational flow around the vortex line. Given this, we represent the vortex as an object in cylindrical coordinates: ψ(r) = f (r, z)e inφ . To describe the density distribution, their structure and size, one has to solve the Gross-Pitaevskii equation for this object, which in the case of the weakly interacting Bose gas is, In the limit n = 0 we recover the standard Gross-Pitaevskii equation. So the term that contains n is called the centrifugal barrier and it is the kinetic energy term from the azymutal velocity of the vortex. Solving this equation gives the density profile of the vortex, from which we can determine the size, density and general structure of the vortex.
In the presence of a trapping or external potential, changes the density of the condensate which also changes the dynamics of the vortices. For details of the formation of vortices in the presence of a trapping potential see also ref. [140].
Since its discovery in Helium 4 superfluid, vortices have been one of the central topics in the research of supefluidity [141]. They have been experimentally observed in many systems like Helium 4 and Helium 3, superconductors and atomic BECs, including multi-component BECs. These new observational advances allowed us to vizualize and study those vortices, as it can be seen in left panel of Figure 7.

Summarizing
We saw in this section an introduction to two of the most interesting phenomena in quantum mechanics, the BEC and superfluids. BEC is a consequence of the statistics of bosons, and that a non-interacting Bose gas can achieve condensation on low temperatures. Superfluidity, following London's idea, arises upon condensation and is described as a fluid that flows without dissipation for velocities below a certain critical velocity, which is the Landau criteria for superfluidity.
We also saw that condensation and superfluidity are general phenomena that can happen in different systems. Here we showed that condensation can arise in the theoretical ideal gas, or in a weakly interacting system of bosons. BEC happens in many different systems that have very diverse dynamics. The same happens for superfluidity that, although it has some conditions for a superfluid to be formed, can emerge in many different systems.
The simplest description of a supefluid is a weakly interacting Bose gas. We showed here how this theory with an interaction describes the ground state, the condensate, in the mean field approximation, and how the excitations describe the phonons. This shows how these phonon excitations in the superfluid behave like sound waves with linear dispersion relation, with sound speed coming from the interaction. This shows that interaction is crucial for superfluidity. A non-interacting Bose gas achieves BEC but not a superfluidity, while an interacting BEC forms a superfluid upon Bose-Einstein condensation. We also showed that this theory has a hydrodynamical description and can be written in terms of fluid variables, where it is evident the main characteristics of superfluids, its non-vorticity. This theory can be described in the field theory language where a superfluid is described as a theory where spontaneous symmetry breaking happens because of the condensate ground state. This description is equivalent to the quantum mechanical one and allows us to evaluate the ground state and excitations, as well as the fluid equations, and is very appropriate to make these calculations.
This theory provides a microscopic description of the superfluid in the case of two-body interactions. We showed that we can also write the superfluid theory as a macroscopic theory at low energies that describes only the phonons. This theory can be obtained by symmetry alone and does not rely on details of the microphysics of the system. This description allows us to describe many different types of hydrodynamics fluids, including more complicated superfluids, superconductors, among other systems where describing the microscopic theory would be either very difficult or impossible. Those many descriptions above show us the richness of quantum mechanical effects that exist on macroscopic scales. All of those phenomena, from BEC and superfluidity, are very well studied in condensed matter physics and were experimentally confirmed and studied in many different systems, from helium-4, helium-3, to atomic BECs, unitary Fermi gases, and many more. These different descriptions are important to describe these different systems that have very different properties.
In the next section we will use this knowledge to describe the different models of DM using ultra-light fields. As we saw from this section, when considering ultra-light bosonic fields, BEC is expected. Now, we can describe different systems with these ultra-light fields, which depends on the dynamics these fields have. In one case we can choose a non-interacting trapped, gravitationally bounded, ultra-light field, which is a BEC. In the other we can describe this ultra-light field as being a superfluid, described by a weakly interacting two-body system. Or, we can even have a DM that is a different superfluid, with a more complicated dynamics, which can be described by the EFT of superfluids. In this sense we will explore the richness of these quantum mechanical phenomena and show that these behave as DM and explain some curious properties of DM on small scales, a consequence of the quantum mechanical nature of these fields.

Ultra-Light Dark Matter
Ultra light DM denotes a class of models where DM is composed by ultra-light bosons that form a Bose Einstein condensate or superfluid on galactic scales. These models were introduced as a new class of DM models with the goal of addressing the small scale challenges of ΛCDM. The idea is that those models behave like CDM on large scales, recovering the incredible observational successes of this description, while inside galaxies DM matter condenses into a BEC or superfluid. This happens due to their very small mass, which results in a macroscopic de Broglie wavelengths that leads to the formation of a condensate core in the inner regions of the halos of galaxies. Because of the presence of this condensate core in the halo, the density profile of galaxies is finite towards its the central regions of the halo, and the density cannot diverge and become cuspy. The profile outside the condensate can follow the expected from CDM, since outside the core DM behaves as a normal fluid.
The wave nature of this condensate core leads to the suppression of the formation of structures on those scales. The condensate is a stable regions where no clustering takes place, suppressing the growth of smallest structures like low-mass substrucutres and sub-halos. A smaller number of those substructures can address some of the small scale challenges of ΛCDM. For this condensate core to form on galactic scales, the mass of the ULDM particles needs to be small. We need to have masses that are small enough so that this particle can undergo BEC and produce cores that are relevant on galactic scales. We can determine the general conditions for this to happen in galaxy halos. From the conditions for condensation from Section 3.1, condensation arises when de Broglie wavelength of the particles start to overlap forming the single wavefunction of the condensate. This happens when the de Broglie wavelength is larger than the inter-particle distance of the condensate constituents, where, assuming an spherical halo, the interparticle distance is defined as the radius of a sphere with density ρ. This gives a bound on the mass of the DM particle. This condition is equivalent to requiring a macroscopic occupation number of the ground state for T < T c , where T c is the critical temperature that depends on the characteristics of the system. This is equivalent to requiring that nλ 3 dB 1, where nλ 3 dB determines approximately the occupancy number (N ).
Since we want the density and velocity of the dark matter halo, we are going to take this bound at virialization. Like described in Section 2, from standard spherical collapse, this is given by [178], where we derived these expressions assuming H 0 ∼ 70 km s −1 Mpc −1 and a halo mass of order of the MW. This gives the bound: Taking z vir ∼ 2 and for halos with mass of order of 10 12 M , we have an upper bound for the mass of the ultra-light DM particle in order to have a condensate core in the galaxy halo: m ∼ 2 eV. More massive halos, like the ones from clusters, will then be in the in the normal phase, where DM behaves like free particles.
We can also impose a lower bound on the mass. The size of the condensate core cannot be larger than the halo, since we want the condensate only on galactic scales and normal DM on larger scales reproducing CDM. The maximum case we can have for the formation of a condensate is where the de Broglie wavelength of the ULDM particle is of the order of the size of the halo. Therefore, there is a bound on how large the de Broglie wavelength can be λ dB < R 200 . Taking again z vir ∼ 2, we can see that for a spherical halo [145], this imposes a lower bound on the mass In this way, in order to have a model where DM forms a condensate core inside the halo of galaxies we need to have a DM particle with masses: Those masses are much smaller than the ones usually considered for DM candidates and cannot be produced thermally in the early universe. Therefore, ULDM is a non-thermal relic, having to be created by a non-thermal mechanism in order to be cold today and be able to form a condensate. There are many non-thermal production mechanism, as cited below, but since in this review we are being agnostic on the type of particle that consists our ULDM, we are not assuming any creation mechanism. One thing that is important to mention is that the calculation showed above uses the condition for condensation that comes from the non-interacting BEC, as we saw on Section 3.1. In this way this is only an approximation given that the DM in halos is not an ideal BEC, and only serve to gives us an order of magnitude estimate for the masses required for the ULDM model. For condensation to happen one also need to have thermal equilibrium. The condensate formed needs to be a stable coherent state, and for that one needs to study the perturbations, like we did in 3.3. We are going to study those conditions for each model in more detail in the next sections, and we are going to see that depending on the different type of BEC or superfluid used, there are other parameters that influence the formation of the condensate like the presence of an interaction. This leads to the possibility of having an condensate for difference ranges of masses, always respecting the bound above which is very general for this class of models. However, it is important to emphasize that to really determine if condensation can happen in the DM halo is a very complicated problem. To properly determine this one would need to determine if there is indeed a macroscopic occupation of the very complicated ground state of the system given by the DM halo. This is still an open question and an active research topic.
As we stated above, given these general conditions of the ULDM class of models, one can construct many different models of DM using all the different BEC and superfluid constructions available. Each of these different types of BECs and superfluids can lead to the formation of a condensate core inside the halo, each with a different condition for that, with distinct and rich phenomenology that can be probed by their different astrophysical consequences. We are going to see now the different types of ULDM models present in the literature and show that they fall into three global categories.

DM relics
Dark matter candidates can have distinct formation mechanisms, with the main ones being thermal or nonthermal relics produced in the early universe. Depending on this mechanism, different masses and couplings for these are allowed as the correct relic abundance of DM is obtained.
Thermal relics: This refers to the particles, including DM, that are produced from the hot and high density thermal bath of the photon-baryon plasma in the universe. Initially, in the early universe, the universe was in a state where it was hot and dense where particle and photons were very close to thermal and chemical equilibrium. This means that the time scale of the particle interactions (1/Γ ) in the plasma are much bigger than the expansion time of the universe H Γ = n σv , where σv is the thermally averaged cross Section [176]. As evolution follows and those quantities redshift, at some point H ∼ Γ , and the particle decouples from the thermal bath (at the temperature T fo ). This is a simplified description of the process called freeze-out or decoupling. Depending on the interaction rate of each particles, they decouple at different times. This process is described by the Boltzmann equation and it is how electrons, neutrons and neutrinos are formed. If T fo m, where m is the mass of the particle, the particle decouples as a non-relativistic particle, and it is called a cold relic; otherwise, if T fo m or T fo ∼ m, and we have hot and warm relics. We can assume that DM is a component that was in contact with the thermal bath and it is a particle produced through decoupling from the thermal bath like described above. For the thermal relics, the particle with smaller mass is hotter. So, from structure formation there is a lower bound on the mass of thermal DM of O(keV ). This considers that all DM is given by this thermal DM, but if this is not the only component for DM, this bound can be evaded. For cold relics, in the case of WIMPS, created through this mechanism, to obtain the observed DM abundance the averaged cross section needs to be of the order of the one characteristic of weak interactions σ DM−DM 10 −8 GeV. This is known in the literature as the WIMP miracle. Many models of DM are produced thermally like supersymmetric candidates, more complicated WIMP candidates or particle DM decays 22 .
Non-thermal relics: As we saw above, there is a limit for the mass of the DM particle that can be created thermally. The only way of having smaller mass candidates of DM is by having a non-thermal mechanism to produce those DM particles. There are many mechanisms that can produce non-thermal DM candidates that include decaying from topological defects, decaying from a massive parent particle, vacuum misalignment, among others [16,177]. Vacuum misalignment or vacuum displacement is one of the genesis mechanism for these ULDM and axions. The vacuum displacement mechanism [220] can be described, in a concise way in the following way. A massive scalar field in an FRW universe, when H > m ϕ , is overdamped and it behaves nearly as a constant. So, if we consider that initially this field was displaced from its minimum, ϕ = ϕ * , the field has a potential energy given by ϕ * . When H ∼ m ϕ , the field starts to evolve and begins to oscillate in its potential, and in turn redshifts like matter. The mass and the initial displacement fix the energy density of this misalignment field. If we consider that the ultra-light particles are created by this mechanism, it imposes a lower bound on the mass H(a eq ) ≈ 10 −28 eV, in order to start behaving like DM around equality. For more detail on this mechanism for axions, see [16].
Classification of ULDM models 22 Even the axion has a thermal production channel, if the axion is in contact with the thermal bath, with this axion being hot and contributing to a fraction of the effective number of neutrinos.
The idea of having DM condensation on small scales is not new and has been around for 30 years [146,147]. These models receive a lot of attention given their rich phenomenology on galactic scales while maintaining the successes of ΛCDM on large scales, and also because of their possible connection with axions, which are a strong candidate for DM. Although this connection is very attractive, these models are much more general and the light DM particles do not need to be axions, only need to be very light particles that can condense into a BEC on small scales.
Nowadays, there are many models that describe a DM component that condenses on galactic scales. These models receive many names in literature and invoke condensation in similar but distinct ways that have important implications for the observables of the models. We can classify these models basically in three categories (which somehow agrees with what was suggested in [148]). This distinction is not always presented in all the literature and sometimes the names of these models are used interchangeably. However, this classification is instrumental since it elucidates the physics responsible for the condensation, and it separates the different phenomenology those model present.
The first category is given by a gravitationally bounded scalar field model. In this model condensation under the influence of the gravitational potential is achieved in galaxies where the gravitational attraction is counteracted by the quantum pressure. This class of model can be called fuzzy DM (FDM), since this name is already very well established for these gravitationally bounded BECs. One of its main realizations, which coined the name fuzzy dark matter is presented in [163,164], where the DM is given by a light particle with m ∼ 10 −22 eV. The FDM model is the ULDM model that was studied the most in the literature both theoretically and numerically. With a particle with this mass, the FDM model is known to be able to solve some of the challenges from small scales presented above, and to be in agreement with large scale observations. Some interesting phenomenology also emerges from this model, as we will discuss in detail in the next subsections, that can be probed by current and future astrophysical observations. This model has also been called in the literature by wave DM, ψDM, among other names [146,147,165,166,167,168,169,170]. The second category is called self-interacting BEC (SIBEC) DM, but it also receives the names repulsive DM, scalar field DM, fluid dark matter, among others [150,151,152,153,154,155,156,157,145,158,159,160,161,162], in the literature. In these models, the BEC is described by a bosonic field with a (usually) 2-body self-interaction. The presence of the interaction makes this model present superfluidity upon condensation. This case is described by the interacting BEC presented in the previous section, which is the simplest example of a superfluid. The presence of the interaction controls the stability of the condensate and for this reason this model presents a different phenomenology depending not only on the mass of the particle, as for FDM, but given the strength and sign of the interaction. For a repulsive interaction, the condensate has a long range coherence and presents superfluidity. The 2-body case is characterized by having an equation of state (EoS), P ∼ n 2 , as we saw in the previous section. Higher order interactions, describe SIBEC with different equations of state.
The case of the QCD axion was studied both in the contexts of SIBEC and FDM [171,172,173] (although the axion presents a self-interaction, sometimes they use the QCD axion masses in the context of the FDM model). The validity of this description as a DM model was disputed in [174], since long-range coherence in the condensate can only be achieved in the presence of a repulsive interaction, and in the case of attractive interaction, like it is the case of the axion, a long-range coherence and superfluidity are not achieved. We will describe this in details in the next subsection.
The third category is called DM Superfluid [178,144,179,180,181]. This theory was proposed with the goal of reproducing the MOND empirical law on small scale, presenting a natural framework for the emergence of this theory. Different than in the case of SIBEC, in order to reproduce MOND on small scales it requires that the EoS is given by P ∼ n 3 , like what is expected by MOND, with a more intricate dynamics describing the small scales. To accomplish that, this model is described using the EFT of superfluids which allows us to describe superfluids with a more general dynamics and EoS.

ULDM categories
Classification is based on the different ways they achieve condensation.
Fuzzy DM (FDM) described by a ultra-light scalar field under the influence of gravitational potential. Forms a BEC on galactic scales.
Self-Interacting BEC described by a self-interacting scalar field with 2-body (or higher) interaction.
This classification divides the models according to the different condensate and superfluid construction they describe. FDM and SIBEC can be described by a microscopic theory for the non-interacting and interacting BEC, respectively. While the DM superfluid model is described by the EFT of superfluids that is represented by the macroscopic phonon Lagrangian, providing a superfluid with a different equation of state. All of them describe distinct system with different properties and dynamics in galaxies. Condensation will arise in each of these models for different values of the mass and of the other parameters of the model. Given these different mechanism to invoke condensation, each of those categories will have a very different mass for the ULDM particle. Those differences are of many order of magnitude, around the lower and upper mass values of the required range (88). For example, as we will see below, for the FDM model the mass is around 10 −22 eV, while for the DM superfluid m ≈ eV, with different condensate core sizes. This shows how important it is to understand the dynamics of each of those categories of ULDM model.
It is important to notice that the classification and the properties of each of this system in this review are only going to be about describing their gravitational consequences. We remain agnostic about the origin of this field. This scalar fields can be axion or axion-like particles, or any ultra-light field that can come from different non-thermal formation mechanisms. Here we are only interested in their description and gravitational consequences, not in any other interaction they might have with baryons. In this sense, this review presents a very general description that any ultra-light particle would present 23 . Now that we classified ULDM models, in this section we describe them in detail, showing their properties and the conditions for them to achieve condensation on galactic scales. These conditions depend on the characteristic of each model but also have to be enforced in order to address the small scale challenges, and to be compatible with the large scale cosmological and astrophysical observations. We also present some of the astrophysical consequences that each of those models predict, and discuss in the end some observations that are being studied that promise, in the near future, to allow us to test these model even further and using different and complimentary probes. To aid in this task it is important to simulate the behaviour of those model in different scales and regimes, and we describe some of the efforts being made in this front.

SIBEC and FDM
In this section we will describe the general features of the SIBEC and FDM models. These models can be described by a similar microscopic theory given by a scalar field that undergo spontaneous symmetry breaking by the condensate ground state, as we saw in Section 3.3. We can describe them as a field-theory of an interacting and trapped BEC, respectively.
We want to describe a class of the models cited above where a light scalar field is the DM component. This can be described as a very light scalar field φ, with a self-interaction that is minimally coupled to gravity given by the action, where g µν is the metric. If we wanted to specialize in the QCD axion, the potential in the Lagrangian above could be thought of coming from the QCD axion potential V (φ) = Λ 4 (1 − cos(φ/f a )) for small field values (φ f a ), where m = Λ 2 /f a are of order O(10 −5 )eV, and g = −Λ 4 /f 4 a < 0. In our case, this is a generic light scalar field. Since the particles we are considering are very light, the Compton wavelength m −1 is much smaller than the particle horizon at matter domination, when we are interested in DM being responsible for the formation of structures and the behavior in galaxies. In this limit, we can use Newtonian approximation and the non-relativistic dynamics for the field. Taking the non-relativistic limit, we re-write the field as: and take c → ∞ (if we recover the c's). The non-relativistic action then is given by where we are using a Friedmann Robertson Walker (FRW) metric given by: This Lagrangian gives the following equation of motion: If we consider time scales smaller than the expansion, we can ignore expansion and write the equation of our system as: which is a non-linear Schrodinger equation. The gravitational potential term can be re-written in the form −Gm 2 ψ d 3 x |ψ(x )| 2 /| x |. This non-linear equation is the Gross-Pitaesvkii equation described in the previous section. We can use this equation to analyse the properties of this system, analytically and numerically.
We can also rewrite the field theory above as a set of hydrodynamical-like equations, in this long wavelength limit. For that, if we identify (using the theory in the presence of expansion): The vorticity of the superfluid is zero and the momentum density has non-zero curl. The comoving equations of motion for ψ are:ρ These set of equations are the Mandelung equations, generalized for an expanding universe. The second term in the right hand side of equation (97) comes from the self-interaction term, where P int is the pressure from the interactions. The last term of the second equation is the quantum pressure. This is present even in the absence of interaction and it is going to be important for the effects and formation of the condensate for the FDM model. The quantum pressure, coming from the uncertainty principle, has the role of not allowing the FDM to cluster and collapse, which also makes the density of this collapsed region to have a finite value. In this way, this model has naturally a cored profile inside the condensate region, addressing the cusp-core problem. This form of the equations is useful for numerical simulations that can reveal some properties of the DM scalar field.

Stability and evolution of the condensate
We saw in the beginning of this section the general conditions for the ULDM to form a condensate in galactic scales. The mass of the ULDM particle needs to be very small, in the range 10 −25 eV m 2 eV for a MW-like galaxy. This means that this DM needs to be produced non-thermally in the early universe. For masses in this range we will have a high occupancy of the ground state and condensation. This very simplistic estimate based on the conditions for condensation for a non-interacting Bose-gas already gives a good estimate for the order of magnitude of masses for condensation to happen in the halo of a galaxy. However, our models of DM describe different types of condensate systems. In the case of the SIBEC and the FDM, both are described by the same microscopic theory, as shown above. But condensation arises very differently in those models. The SIBEC has the presence of self-interaction that will result in a pressure term that counter-acts or helps the action of the quantum pressure. While for the FDM model, a gravitational bounded BEC where interactions are absent, the condensation in the system will arise as a balance between the quantum pressure and the gravitational collapse. Therefore, these two systems have very different condensate states given by their background solutions. The conditions for these condensate state to exist depends on the behaviour of the perturbations of the system that shows the stability conditions for this system. This determines also the size of the condensate core that can be formed. This analysis was already seen in Section 3.3, when we analyzed the weakly interacting Bose-gas. Here the analysis has to be conducted in a similar manner, but now with the presence of a gravitational potential in order to understand the properties and conditions of the condensate formed in these models. This is what we investigate next. We are going to study now the evolution around the background solution that describes the condensate and the stability of this solution done by studying the spectrum of the perturbations around the condensate. The coherence size of the condensate is going to be determined by the stability conditions. Having determined the stable configurations that can lead to condensation, the ground state solution that describes this condensate is showed. We want to calculate this for the cases of SIBEC and for the FDM, since they are described by the same microscopic theory in the presence of self-interacting or without interaction and in the presence of a gravitation potential, respectively

SIBEC: self-interactions
We are first going to treat the SIBEC model, which is described by the presence of a self-interaction. We specialize here to the two-body interaction, since higher order interaction are usually suppressed in low-energy systems like the one we are interested. But it is easy to generalize this for higher order interactions. From what we saw in the previous section, the theory of a self-interacting condensate describes a superfluid in certain regimes. We are going to see here for which conditions this occurs in the SIBEC model.
We are going to work in the limit where we ignore gravity in order to investigate the effect of the self interaction in the model in a very similar way as described in Section 3.3. The Schrdinger equation that describes model is given by: We decompose the field into a homogeneous background solution, which represents the condensate, plus a perturbation part: ψ(x, t) = ψ c (t) + δΨ (x, t). The condensate part satisfies the Gross-Pitaevskii equation, that has a simple periodic solution ψ c (t) = ψ 0 e −iµct , where |ψ 0 | 2 = n 0 is the number density of particles that fixes the amplitude of ψ 0 , and µ c = gn 0 /8m 2 . The equation describing the evolution of the perturbation, making the field redefinition δΨ = ψ c δΨ , is Since δΨ is a complex scalar field, we can decompose the field into a real and a imaginary parts, Ψ = A + iB. We want to determine the dispersion relation of this system, so we write the equation of motion in Fourier space: where we call ζ k = k 2 /2m + gn 0 /4m 2 . The dispersion relation is given by: We an see that for ω 2 k > 0 we have an oscillatory solution: where Z is an arbitrary complex parameter. When ω 2 k < 0, the solution of the equation for δΨ k is given by exponentials, where γ k = (k/ √ 2m) √ −ζ k are the eigenvalues of the matrix, and c 1 and c 2 are constants given by the initial conditions. The regimes where the dispersion relation is positive or negative are separated by the modes with wavenumber, from where we can also determine the wavelength λ * = 2π/k * that divides these regimes. This wavelength is proportional to the healing length calculated in the previous sections. However, can see from this scale that what is actually going to determine if we have a stable oscillatory solution or a exponentially growing instability is the sign of the interaction. For a repulsive interaction, g > 0, the homogeneous configuration is always stable, and it is always going to be described by an oscillatory solution, either if λ is bigger or smaller than λ * . This means that for all wavelength we can have a stable solution that can describe a condensate. From the dispersion relation we can also see that, in the long wavelength regime ω k c s k , which is the dispersion relation of the phonon, that propagates as a wave, mediating long range correlation. This means that the SIBEC with a repulsive interaction is a superfluid, with the theory fully described by this propagating phonon. On the other hand, for very small wavelengths, large k, the term ω k = k 2 /2m dominates, which is the dispersion relation of a free massive particle, and the system stops exhibiting superfluidity. In the intermediary regime, where we should consider both terms of the dispersion relation, and the theory described by two degrees of freedom the phonon and the massive particle associated to particle creation away from the condensate. The scales that determine what we mean by long-wavelength regime is where λ λ * in a way that the linear term dominates the dispersion relation. Since λ * is proportional to the healing length, this condition for the long-wavelengths is equivalent to the condition in Section 3.3 that the healing length gives us the scale where quantum pressure (QP) can be neglected, which is what we are describing here. When we have an attractive interaction, given by g < 0 we have two regimes of stability. For k < k * (λ > λ * ) we have exponential growing solutions, which means that perturbations grow parametrically. Given this instability, the condensate cannot be formed on these scales. For k > k * (λ < λ * ), the solution oscillates and is stable, forming a condensate. This stable configuration, however, is different than in the case for repulsive interaction, forming a localized object, with maximum size given by λ * . This localized stable solutions are called soliton. Therefore it makes sense that this is the healing scale is the scale of the stability, since this is the scale below which the interaction "heals" perturbations of the condensate.
The case of the attractive interaction is not a superfluid even in the stable localized regions, the solitons. The only stable regions are for λ < λ * and this is the regime where the linear term is always subdominant in the the dispersion relation. The long-wavelength regions in this case are the regions where instability happens and there is no formation of a condensate.

Condensate solution:
Having studied the stability of the system and determined the regimes where we have stable and unstable solutions, we can now describe the background solution for each case. As we saw in the "condensate solution" part of Section 3.3, one possible solution for the weakly interacting BEC is given by the solitons. For repulsive interactions the condensate solution is a dark soliton, while for a attractive interaction one has a bright soliton.
What is important is that for the repulsive interaction we can create a condensate of any size desired. So the size of the condensate is not limited and it will depend only on the choice of mass and strength of the interaction. For the case of attractive interaction, we can only have localized stable solutions, giving a finite size bright soliton, with maximum size λ * = 2π 2m/|g| n 0 . Therefore we can see that we could only have a soliton, for an attractive interaction, that is relevant on galactic scales if the mass is very small but with a not very small coupling. An interplay between these two parameters need to occurs then if one wants galactic sized λ * .
This condition also shows us that the QCD axion is not a good candidate for a ULDM (at least not to represent most of the ULDM, but it can still represent a fraction of DM). For the QCD axion the interaction, given by g a = −Λ 4 /f 4 a is negative, and it is extremely small with g a ∼ −10 −48 , because Λ ∼ 0.1 GeV, for the typical QCD scale, and f a ∼ 10 11 GeV, for typical Peccei-Quinn scale. Given that the mass of the axion is approximately m ∼ 10 −5 eV, the soliton length, for n 0 ∼ n gal , is λ s ∼ 2.8 × 10 11 km ∼ 9 × 10 −6 kpc. This is much smaller than any galactic scales, which are of the order of tens to hundreds of kpc. The QCD axion, then, produces these small and localized clumps of axions. One thing that is being investigated by some researchers is the possibility of many of those small cores to be created in a halo, with the possibility of them clumping together to form a bigger condensate.
Occupancy number evolution: With those solutions in hand, we can understand how the evolution of the occupancy number for the condensate will behave for each mode. Determining the evolution of the occupation number is very important since having high occupancy number of the ground state of the system it what is actually the definition of a condensate. And this definition is independent of the system described, representing the best way of showing that condensation took place. This is given by: N = |ψ k | 2 /V , where ψ k is constructed from the exponential and oscillatory solutions described above, with a random phase. The average occupation number evolves as: for ω 2 k > 0 For g > 0, an repulsive interaction, γ k is imaginary, with γ k = iω k , so the occupation number oscillates and the oscillations are stable. The ratio N k (t) / N k (t i ) which has the largest value is obtained for modes that minimize ω k , which are the modes with k → 0. These are the longest wavelengths. Since the long wavelength dominates, this means that long range correlation is present, and we can have a long range condensate. For g < 0,the occupation number grows exponentially. The fastest growth is given by the modes k = k * that maximize γ k . So the modes k > k * or λ < λ * , where * denotes the characteristic scale where instability sets in, will dominate and the stable configuration of the system will be localized clumps. The size of these clumps will be given by the mass and interaction of the model.

FDM: only gravity
We are now going to describe a model without interaction, where the ultra-light particles are under the influence of the gravitational potential. We can think that gravitational potential has the same effect as an attractive interaction, in a way that quantum pressure has to counter-act the gravitational collapse. This gives a good picture to what to expect, from the knowledge obtained for the SIBEC. However, one needs to remember that the FDM model is described by a non-interacting theory, which means that it condenses into a BEC, but does not exhibit superfluidity in any regime.
The Schrdinger equation for DM in a gravitational potential, in the absence of interaction, is, which is coupled to the Poisson equation: where the average background density,ρ, was subtracted. Expanding the field as done previously ψ(x, t) = ψ c (t) + δΨ (x, t), the equation for the condensate is trivial and ψ c = ψ 0 = const.. For the fluctuations we can write the linearized systems of equations that govern the evolution of the perturbation: These can be combined into the equation: One can notice that the equation above is very similar to the equation we had for the interacting case (100) for an attractive interaction. With that, we expect that there is an instability for long wavelengths, and that the condensate stable solution is only given for a finite region, forming a localized core. To determine this lets take the Fourier transform of the fields. Like in the interacting case, the instability is divided by the regimes where the dispersion relation is smaller or bigger than zero. For the parameters of the FDM, we can determine the wavenumber that separates the regimes as:ω This scale is the Jeans scale and it separates the regimes where gravity dominates and collapse happens (k < k J ), and the regime where the quantum pressure dominates and the solution is stable and oscillates (k > k J ). In this regime we can have a condensate. The quantum pressure term counteracts the gravitational attraction and any attempt to localize the particle is accompanied by an increase in energy. So stability below the Jeans scale arises because of the uncertainty principle. We can estimate the size of the coherent condensate core. Rewriting the Jeans wavelength as: forming a condensate that is of the order of the scales of the halo of galaxies. As we can see, for a MW-like galaxy, this core formed is smaller than the halo. So we expect that in the outskirts of the halo the DM is not going to be condensed and is going to behave like normal matter, with the profile following the NFW profile. For the QCD axion, if we assume that we can have an axion without interaction, λ J ∼ 1.7 × 10 −7 kpc, which is a very small scale in comparison to galaxies. This stable bound system is called a Bose star.
Adding the expansion of the universe, from (91), we can see that the system of equations becomes Like before, we determine the linearized equations for the perturbations δΨ around the coherent homogeneous background that evolves as ψ c ∝ a −3/2 . The equation that describes the evolution of the perturbation then is where Ω a = mn 0 /ρ tot is the density parameter of our FDM particles. The real part of the perturbations obeys the following equation:Ä All the terms in this equation would also be present from a massive DM component like CDM, except for the last term, which is related to the quantum pressure. The Jeans length in this case is given by: This stability analysis shows us that if we want to construct a ULDM model with the desired feature of having a condensate in galactic scales, for the FDM model we need to have a very small mass, of the order of m ∼ 10 −22 eV. This mass is within the bounds we had obtained earlier for all the ULDM models. For the FDM model the only parameter that controls the size of the condensate is the mass. We need then to determine the mass of the FDM particle for which this model can address the small scale problems, and for which the model is still in accordance with the very precise cosmological observations. This is what we show in the next subsection. We are going to put bounds on this parameter according to the observations in galaxies and on large scales.
In the case of the SIBEC, we saw that what determines the size of the condensate is the mass and of the strength of the interactions. Like for the FDM, in order to construct a SIBEC model that addresses the small scale problems and that allowed by current observations, we need to put bounds on those parameters. This model was much less studied than the FDM model. There are a few works that present some bounds on those parameters (see references in the definition of this class), but with much less constraints than in the FDM case. This is also understandable since here we have an extra parameter in the model. A more complete analysis for the bounds of the SIBEC is going to be presented in a future publication and it is being considered by the authors while this review is being written.

Fuzzy Dark Matter
We studied above the evolution of a light scalar field behaving as DM in the presence of an interaction and under the influence of the gravitational potential. At late times, the dominant interaction is given by gravity, so it is natural to study the light scalar field in the absence of interaction. As we saw, in this regime the stability below the Jeans scale arises because of the uncertainty principle, and we have the formation of a condensate that has a maximum size given by the Jeans length. The core is going to be smaller than the halo, so in the outskirts are going to be described by normal particle DM with a NFW profile, like described in a schematic way in Figure 8. This size can be big if the mass of the field is very small, much smaller than the QCD axion mass.
We need now to determine for which masses this model can address the small scale problems of the ΛCDM, and determine if these values are compatible with current cosmological and astrophysical observations. After that we describe some astrophysical consequences and possible observables this model yields [164] 24 .

Astrophysical and cosmological bounds of the FDM model
We want to determine the mass of the FDM candidate. From the discussion above, we saw the there is a bound for the mass in order to condensate in galaxies 10 −25 m 4.7 eV for a typical galaxy MW-like galaxy. Later we saw that, for masses of order of the usual QCD axion mass, around m ∼ 10 −5 eV, the stable configurations are very localized and small, far from galactic scales. With that we can already see that m 10 −5 eV for the FDM model. Now, we are going to see other conditions that can bound the mass and show the mass range for the FDM particle that gives the best DM phenomenology.

-Cosmological constraints:
As we calculate above, the Jeans wavelength given by the FDM model is: where we assumed that ρ =ρ and the Ω m ∼ 0.35. The Jeans length it is known as the scales that determines the gravitational instability for the matter perturbations in the universe. Growth of structures occur for λ > λ J , and no structure formation happen for smaller wavelengths. For CDM, for example, since the sound speed is effectively zero, which leads to zero Jeans length, so structure for for all wavelengths. In the case of FDM, the Jeans length has a finite size, so structures on smaller scales are suppressed, in comparison to CDM. This scale is solely determined by the mass of the particle. Therefore, if we want to have a condensate that is relevant on galactic scales, in a way that suppresses the structures in the astrophysical scales in order to try to address the small scale challenges, we can impose that the Jeans length must be minimally on the kpc scales, therefore λ J > O(kpc). This implies that the mass of the FDM particle needs to be m 10 −20 eV, which already makes the range of possible masses for the FDM model much smaller. For this mass, the Jeans wavelength is λ J ≈ ×10 2 h −1/2 kpc. So modes with λ smaller than this Jeans length are going to form a stable condensate, and structure below those scales is going to be suppressed. This can be seen, for many different masses in the predictions for the matter power spectrum, on the left panel of Figure 10. The influence in the power spectrum on large scales is going to be more pronounced as the mass of the FDM particle is smaller, and this are the regions where we have very precise measurements of the power spectrum and that we know it should behave very close to ΛCDM. With that we can use current cosmological data to put bounds on the mass of the FDM model. In [23], the authors investigated that using a combination of CMB data from the Wilkinson Microwave Anisotropy Probe (WMAP), Planck satellite, and also from ground CMB experiments like the Atacama Cosmology Telescope, and South Pole Telescope, and galaxy clustering data from the WiggleZ. With that, as we can see in Figure 10, they concluded that given this suppression, the mass of the ultra-light field needs to be in order to this field to behave as dark matter and to reproduce the observations. A FDM model with mass in this range is going to be indistinguishable of CDM on large scales, which is what we wanted. If m 10 −32 eV the ultra-light field can behave as dark energy. We will discuss this in more details in the end of this section. This suppression of the power spectrum on small scales also suppresses the formation of galaxies. It is found in simulations that the number of sub-halos in FDM in comparison to CDM is reduced by a factor of ∼ (3M/M 1/2 ) 2.4 . This suppression of formation of small galaxies is larger in FDM at higher redshifts, in comparison to CDM. This opens up an important question about FDM being able to produce small scales structures at early times to be probed by Lyman-α forest.
-Halos: minimum size, maximum density and the cusp-core problem One limit to the size of the condensed core is that it needs to be smaller than the virial radius, λ dB < R = GM/v 2 , the same conditions we used to put a lower bound on m in order to have condensation in galaxies. Now, we want to use this condition to obtain the maximum size for the condensate core: R 1 GM m 2 , where M is the mass of the galaxy. We can write that in terms of the half-radius, radius where half of the mass of spherically symmetric system is contained: This bound on the radius is compatible with the half-light radii inferred from the densities of 36 Local group dwarf spheroidals [185] if the mass of the FDM particle is m ∼ 10 −22 eV. With the above condition, we can also compute the upper bound in the central density: If we compare this bound to the observations from 8 dwarf spheroidals, and we can see for the density to be comparable to the one measured for these dSph, the FDM mass needs to be m = 8 +5 −3 × 10 −23 eV for Draco and m = 6 +7 −2 × 10 −22 eV for Sextans [186]. For those masses, the distribution at the center of the galaxies seem to be cored, alleviating the cusp-core problem. So, taking a fiducial mass of m ∼ 10 −22 eV seems to be compatible and explain the observations from dSphs.
-Lower bound on the FDM halo masses and the missing satellites problem As we saw before, for the self-gravitating FDM systems, since gravity is attractive, we have coherence on a finite scale. The size of this core depends on the mass, being larger as the mass gets smaller. So the smallest radius to be produced in the FDM model are determined by the smallest mass allowed for the particle. Having a limit for the smallest cores that can be created has important consequences in the abundance of low mass halos, and it is going to be different in this model than what is given by ΛCDM. We can see that by calculating the smallest structures formed in the FDM model. This is given by when λ = λ J , where λ J represents the last scales that can suffer gravitational instability. With the Jeans length, we can calculate the Jeans mass: This is the minimum mass of substructure created in the FDM model. This is in contrast with CDM, where halos with mass below ∼ 10 8 M are highly created, with abundance dn(M h ) ∝ M −2 h dM h . In that sense, the missing satellites problem is addressed in the FDM model, since halos of smaller masses than M J are not created, if m = 10 −22 eV. If the mass of the FDM is smaller than that, sub-halos with smaller mass will form, and the FDM model is not going to address the missing satellite problem anymore. Therefore this shows that in order to solve this problem, the FDM has a preferred mas of around 10 −22 eV. We are going to see below that tidal disruption can also act suppressing small mass haloes, aiding FDM in solving the missing satellites problem. The too big to fail problem is also addressed by the FDM, since we have a mechanism to explain the fact that low-mass sub-halos are not formed, we do not need to invoke mechanisms that creates the too-big-to-fail problem.

-Maximum soliton mass
Since this self gravitating system produces localized cores, we can find the maximum mass for a FDM soliton in an analogous way as the Chandrasekar mass for self-gravitating fermion system. This is estimated making the calculation using general relativity [187,188,189] giving M max = 0.633 c/(Gm) = 8.46 × 10 11 M 10 −22 eV/m . If self-interaction was considered, this mass would have been smaller. This is the maximum mass of the soliton. The soliton is only part of the halo, localized in the central part of the halo.
If we combine all of those conditions, we an see that the preferred value for the mass of the ultra-light particle in the FDM model that is in accordance with CMB and LSS galagy measurements, and that addresses the cusp-core and the missing satellite problems is around m ≈ 10 −22 eV.
However, from recent measurements from Ly-α forest, which provides a data from intermediary scales, closer to the scales where these models start changing its clustering behaviour, a new lower bound for the mass of the FDM particle was found that presents a challenge for the FDM model.

Lyman-α constraints on FDM
Recent investigation of FDM models in light of Lyman-α forest finds new constraints on FDM model parameters [195]. It puts a bound in the mass of FDM in the case where more than 30% of the DM being composed by this scalar field of m 10 −21 . This value is larger than the value necessary for the FDM model to solve the small scale problems of ΛCDM, and puts some tension in the FDM model.
The Lyman-α forest is produced by the absorption of the light from quasar by clouds of neutral hydrogen localized at low-redshifts, in the line of sight between the quasars and us. For this reason, this is an important probe of the matter spectrum on small scales, on scales of order 0.5 Mpc/h λ 100 Mpc/h.
In this work, data from the XQ-100 survey was used which refers to 100 medium resolution spectra with emission redshift 3.5 < z < 4.5. This data is compared against a simulation of the FDM with different masses and abundance today. On non-linear scales, quantum pressure is added. The result is shown in the left panel of Figure 11. The right panel of this figure we also see the impact of the constraints obtained in cosmology and in the astrophysical implications. In cosmology, the constraints obtained give a bound in the value of the displaced field, assuming that the genesis mechanism for this light field is vacuum displacement. Combining this data with CMB data, they also derive bounds on inflation, more specifically on r the tensor to scalar ratio for an inflationary epoch in the presence of FDM. They also show how this bound impacts the resolution of the small scale problems presented by FDM. The cyan line indicates the bound where the missing satellites problem is solved by FDM. The constraint is very tight and it shows a tension with the Lyman-α measurements.
A possible caveat from this analysis is that they use hydrodynamical simulations, and they mostly neglect quantum pressure. However, quantum pressure can be very important and play a vital role in structure formation, which is what the analytical behaviour seems to show us [196]. Therefore, new and independent analysis need to be done in order to confirm if the intermediary to small scales hold more information about these models. Another probe that can help with that is the 21cm from neutral hydrogen, since it probes even smaller scales. We will briefly talk about that in Section 4.4.

Astrophysical consequences of FDM
Here we show some of the astrophysical consequences of the FDM model, using the parameters constrained in the previous section. These are usually consequences of the fact that DM forms a condensate core in the inner regions of galaxies. The hope is that these astrophysical consequences that help constraint the mass of the FDM model.
-Dynamical Friction: The change in dynamical friction is one of the most interesting consequence in BECs and superfuids. This emergent phenomenology can lead to consequent observations that might reveal the characteristics of those systems. It is interesting to see how dynamical friction behaves in the presence of a BEC core, as in the case of the FDM. The dynamical friction in a BEC is not expected to change as dramatically as in a superfluid, which has no friction, but some change is expected nevertheless. It is expected that the FDM changes this prediction because of three phenomena: (i) change in the rate of orbital decay because of the presence of the condensed core; (ii) since the FDM produces a homogeneous core, a mechanism similar to the "core stalling" observed in N-body simulations can take place and reduce or eliminate drag from dynamical friction; and (iii) the way dynamical friction is calculated must be modified by the presence of an object with large de Broglie wavelength, a quantum mechanical extension to the calculation of dynamical friction must be done. An interesting puzzle that can potentially be explained by a modified dynamical friction is the puzzle of the Fornax globular clusters. From the standard dynamical friction expected for CDM and baryons it is expected that globular clusters orbiting Fornax should have rapidly fallen towards its center to form a stellar nucleus. However, there is no signal of mergers and we detect 5 globular clusters orbiting Fornax.
In [164] only the last effect is described and simulated for different parameters the orbital decay times for Fornax in CDM and in the FDM cases. They found that in FDM the orbital decay time is longer, and four of the five decay times simulated are bigger than 10 Gyr or more, thus explaining the puzzle for why the globular cluster in Fornax survived. More simulations and observations are needed in order to confirm this claim, but the FDM model seems to address the dynamical friction puzzle. The ideal is to have the microscopic theory describing dissipation in the FDM model. Finally in [190], a detailed analytical and numerical study of dynamical friction in the FDM model was performed.
To describe the dynamical friction in this model, they describe the dissipation that a perturber moving in a condensate causes. They work this dissipation theory for point-sources (satellites), extended satellites and pointlike satellite in a FDM background with finite velocity dispersion. This analytical theory is then verified by their numerical-simulation that solves the Schrdinger-Poisson system in the presence of such perturber satellite, showing good agreement with the analytic methods. This framework is applied to the cases of the Fornax globular cluster, but also to the Sagitarius (Sgr) stream and the Large Magellanic Clouds (LMC) (we will see about those systems in more detail in Section 4.4). For the Fornax, they find that if the mass of the FDM model is m 10 −21 eV this model stops explaining the Fornax globular cluster merging times, which is in agreement with the mass bounds necessary for the FDM to solve the small scale problems. For Sgr and LMC it is found that the dynamical friction on those are described by the classical limit, described by the Chandrasekhar formula. More simulations need to be done to confirm this, and it is very important to understand the dynamical friction in this regime since these bodies have a strong influence in the MW, and are the target of many studies and observations. These results already show that this modified dynamical friction in the FDM model can maybe explain some interesting astrophysical observations which might be a good opportunity to measure the ULDM on galactic scales.
-Most massive halos -clusters: For distances larger than the de Broglie wavelength of the FDM, it is expected that DM behaves as standard CDM and that the halo enveloping the soliton has a NFW profile. This can be seen numerically for the mass range 10 9 M M vir 10 11 M [191,192,193,194], which gives an estimate for the mass of the central soliton. Extending this relation to larger halos, the FDM predicts that in the inner regions of clusters there will be a condensed core, a soliton, with mass: which is still below the maximum mass for the soliton calculated above. The corresponding half-mass radius is: So, the presence of this soliton with this mass and size would be a prediction of the FDM model. But a question that remains to be answered is the following: the presence of such solitons in the interior of clusters halos is compatible with merging cluster like the Bullet cluster or the anit-Bullet cluster?
In [164] they ask the question if solitons in the center of the galaxies have not been misinterpreted as super massive black holes. They compare the mass of the central dark region measured from Virgo and show that this is similar to the mass of the soliton core in a galaxy like Virgo for a mass of the FDM particle of m ∼ 10 −22 eV. However, since the observations from the Event Horizon Telescope of the black hole in the center of M87 were released, this hypothesis seem to be almost excluded, and it is indeed a super massive black hole that inhabits the center of this galaxy. We have to wait for more data to confirm this, but this is an exciting measurement that can also be used to test the FDM hypothesis.
Another interesting fact is that we know that the galaxies host a super massive black hole in their center. For this reason in [164] they investigate the possibility of a super massive black hole to be created in the center of a soliton. Apparently, the black holes do not grow for the fiducial mass, in a condensate core. Their creation only starts being significant for m 5 × 10 −22 eV, which is in tension with other bounds in the mass, like the one to solve the missing satellites problems.
Those observations show us that small scale observations of very varied phenomena can help us put constraints in the FDM model. We will see in Section 4.4 that other observations might also have the potential to constraint the ULDM models, and their usage for that is under investigation.

DM Superfluid
In this section we are going to describe the third category of models of ULDM, the DM superfluid. In previous sections we saw the small scale problems of ΛCDM and how MOND empirical law offered a very good fit to the rotation curves of galaxies and the scaling relations that emerge from the dynamics of galaxies, which might be challenging in the context of ΛCDM. However, as we saw there is no present framework that can explain MOND, given that the initial proposed theory, the full MOND, and its extensions present serious problems. We present here an alternative model to DM that has the goal of reconciling CDM and MOND: the DM superfluid. This model intends to not only solve the problems that the previous models attempt to address but also offer a mechanism to describe MOND on small scales. In this framework, DM behaves as standard CDM on large scales, while the MOND dynamics emerges on galactic scales. And this is possible through the physics of superfluidity.
The idea of the DM superfluid model, is that DM forms a superfluid on galactic scales, where superfluidity arises upon condensation. This superfluid core present in the inner regions of galaxies not only addresses the small scale challenges of ΛCDM in a similar way as the previous models, but the superfluid described in these regions behaves following a different dynamics which reproduces the MOND behaviour. This is possible given that this superfluid is described by a Lagrangian similar to thee one from the MOND theory which is allowed given the EFT description of the superfluid. In this way, the long-range correlation present in MOND is going to be given by the behaviour of the phonons which mediate long range forces. Outside the superfluid core, DM behaves as normal matter as in the previous models. This is the general vision how the DM superfluid model attempts to describe the behaviour of DM in galaxies.
In the following subsection we will construct the DM superfluid theory showing first in which conditions DM condensates on galactic scales, following that we will present the theory that describes this superfluid phase. With that in hand we can calculate the halo profile and rotation curves in order to compare with data and check the fit of the theory. We present how this model explain many astrophysical systems and possible predictions. After that we show the limits of validity of this description and its relativistic completion. We briefly describe how the cosmology works in this model. This model is constructed by using the fact that in a generic superfluid we can reproduce the exact action as in MOND. This is a very specific model and serves as a toy model for the understanding of theories of DM that present an emergent dynamics on small scales. It is important to point out that what is important is for the model to be able to explain the rotation curves and scaling relations described by galaxies. This translates into a theory that exhibits long-range correlation on small scales. We present here a very specific example of a model when this occurs, where the way the long-range correlation is obtained is by restricting the Lagrangian to behave like MOND on small scales, but this does not have to be the case. The search for a more general theory where this emerges, with a known microphysics is the final goal. Such a construction is currently being searched.  Fig. 12 Approximate calculation form [144] for the fraction of the particles in the condensate versus in the normal state, for a series of sub eV masses in accordance with our bounds. We assumed z vir = 0.

Conditions for DM condensation
Before describing how the DM superfluid behaves inside galaxies, we need to determine in which conditions DM condensates into a Bose Einstein condensate in galaxies. As we saw in the previous section, two conditions need to be met for condensation: first, we need that all the particles are in a single coherent quantum state, described by a single wavefunction of the condensate; a second condition is that the DM particles are in thermal equilibrium, in order to be described by a Bose distribution.
In this section we want to obtain a rough estimate of the bounds in the parameters of the model in order to obtain this condensed core in the inner parts of the galaxy. For that, for simplicity, we use the criteria for weakly interacting gases.

Condensate wavefunction and thermalization
We showed in the beginning of this section the condition on the mass for the ULDM particles to condense in galaxies, showing that they should be in the range 10 −25 eV m 2 eV. This is an approximate condition for the case of DM in galaxies, but it gives us an order of magnitude estimation for the mass of the DM superfluid particles. However, if we remember from Section 3, there is a second condition for condensation.
The second condition to form a condensate is that the particles are in thermal equilibrium. The condition to achieve thermal equilibrium is that the time scale of thermalization must be smaller than the time scale where dynamical processes happen in the halo, the dynamical time. If this condition is satisfied and thermal equilibrium is achieved, the condensate is coherent in the entire halo. The time scale of thermalization if given by the inverse of the self-interaction rate, and the condition for thermalization is given by: where Υ ∼ ρvir m (2π) 3 (4π/3)(mv) 3 is the Bose enhancement factor, which tells you if a boson is already in the state, the probability to another boson to be in that state will be enhanced by a factor of Υ . The dynamical time is taken here as the time it takes to a sphere of density ρ to collapse due to gravity. This condition gives a bound in the self-interaction cross section: where we assumed z vir = 2 and M = 10 12 M . If we want that our self-interaction satisfies the merging-cluster bound [198,199,200], which is ∼ 1cm 2 /g, this gives another bound in mass of the superfluid: m eV. From these conditions, we can obtain a few properties of the our DM superfluid condensate: -Critical temperature: With DM in thermal equilibrium, the temperature can be obtained by the equipartition theorem: k B T = m v 2 /2, which is valid for temperatures smaller than the critical temperature. Above that temperature, the condensate is broken. The critical temperature T c is associated with the "critical" velocity v c , than can be read when we saturate the bound (84): With that, the temperature in a halo of mass M is given by: -Condensate fraction: As we saw on Section 3 at T = 0 it is expected that almost all the particles are in the condensate. However, at finite but subcritical temperature, as seen in Landau's theory [104], it is expected that the fluid is going to be a mixture of superfluid and normal fluid, with the majority in the superfluid. Borrowing from the non-interacting BEC description, this can be estimated as: This formula is only valid for free-particles, and a particle with interaction and trapped in the gravitational potential has a different power than 3/2, but it serves as an estimate. We can see in the Figure 12 the condensate fraction with respect to the mass of the halo for different masses of the DM particle. From this we an see that for a MW-like galaxy for any mass smaller than V, the particles are in the condensate, while for higher masses of clusters, for example, this is not true. Therefore, the mass range from the bound (88) describes a condensate that condenses on galactic scales.
The above conditions were obtained assuming that the condensate will take the entire halo. However, as mentioned in Section 2.1, virialization occurs through violent relaxation, which is an out-of-equilibrium process. In this way, the DM superfluid cannot thermalize. What should happen is that first, the halo virialized and the profile is the expected NFW. After this process, DM particles start to enter thermal equilibrium in the inner, most central regions of the condensate, where the interaction is more pronounced. In this way, the halo would have an inner region (r < R T ) where DM is in a condensed state surrounded by the outer part of the halo (r > R T ) that follows the NFW profile [181]. Since in this model the goal is to be able to describe the rotation curves of galaxies, R T needs to be larger than the radius where the circular motion of stars and gas is observed. For r > R T , the density profile of the halo follows the NFW profile, ρ ∝ r −3 . So we can rewrite the density and velocity with respect to the virial quantities used above: which tells us that it is easier to reach thermal equilibrium in the center of the galaxies where the density is higher. This translates into a bound to the termalization radius: For a MW-like galaxy with M = 10 12 M if we can measure the the circular velocity up to approximately 60 kpc, this will translate into a bound for the mass: m 4.2 σ/m eV.

Superfluid dynamics
Since we have determined that DM condenses and forms a supefluidity in the central regions of the halo, we now need to describe the evolution of this superfluid inside this region. We need to determine the dynamics of the superfluid in order to be able to calculate the profile of the region of the halo comprising the superfluid and with the calculate the rotation curves of galaxies. In this section we will describe the effective field theory of superfluids and show how this is theory reproduces MOND at small scales. As we saw in the previous section, a superfluid at low-energies is described by the effective Lagrangian that is invariant under shift and Galilean symmetries: where Φ = −GM (r)/r is the external gravitational potential for a spherical symmetric static source. The thermodynamic pressure is given by P .
We want our theory to describe the MOND dynamics at the regions where it is superfluid. Given this general Lagrangian for the phonons (135), we want it to describe the MOND action (7). For this, we conjecture that our phonon action is given by: This fractional power might seem strange from the point of view of a quantum field theory of fundamental fields, leading to superluminal behaviour and caustics. However, as a theory for the phonons this is not problematic and it determines uniquely the equation of state of the superfluid. As we can see for the condensate, the background, where θ = µt, the pressure is given by the Lagrangian density: where, in the non-relativistic regime, ρ = mn and n = ∂P/∂µ is the number density of condensed particles. As expected from the result from MOND, this Lagrangian gives us an EoS for the superfluid P ∝ ρ 3 , which is what we wanted to reproduce MOND.
To evaluate the excitation spectrum, we write the action for the phonon excitations φ, that can be obtained by expanding (135) to quadratic order. Neglecting the gravitational potential: from where we can infer the sound speed of the phonon excitations: However, only those ingredients are not enough to reproduce a MOND-like force and a coupling between the phonons and the baryon density needs to be introduced: where α is a dimensionless coupling constant. Although necessary in order to obtain the MOND regime, this interaction Lagrangian breaks shift symmetry softly, only at the 1/M pl level. This term is here considered as a phenomenological term in order to reproduce MOND. In this way, this superfluid theory has 3 parameters: the mass m, the scale Λ and the coupling α. The present form of the Lagrangian to obtain MOND is not the only way of obtaining the MOND behavior in the context of the DM superfluid model. In [179] it was used higher order corrections to generate the non-relativistic MOND action, which is inspired in the symmetron mechanism. Using the same Lagrangian (136) as the leading order Lagrangian, higher order corrections involving gradients of the gravitational potential are added to effectively modify the gravitational force. This results in the spontaneous breaking of a discrete symmetry. The symmetry is broken for small accelerations leading to MONDian gravity, and is restored in the limit of large acceleration leading to Newtonian gravity. In this theory the shift symmetry of the entire system is maintained. A difference from the present mechanism, as we are going to see later, is that cosmologically all the DM is in the normal phase, behaving like CDM, and reproducing all the results from ΛCDM. Here we will describe the method of adding a photon-baryon coupling since this was studied in more detail in the literature.

Finite Temperature
The theory developed above is valid for a T = 0 superfluid. However, in reality, the DM in galaxies has a non-zero temperature. As we mentioned in Section 3, for finite temperatures, this Lagrangian needs to receive finite temperature corrections. In Landau's model the finite temperature superfluid consists is described as a two-fluid theory where a superfluid component and a normal component are present. Those components must interact with each other. Following this idea from Landau's theory, at lowest order in derivatives, we can write the general form of the EFT at finite temperatures and finite chemical potential as a function of three scalars [204]: where X = X(θ) was defined before with respect to the superfluid variables. The other new components are: B is defined with respect to the normal fluid three Lagrangian coordinates ϕ I (x, t); and Y represents the scalar product of the normal and superfluid velocities: where u µ is the unit 4-vector from ϕ I (x, t), and in the last equality of Y we have taken the non-relativistic limit, so v is the velocity vector of the normal fluid component. There are many ways to construct the finite temperature operators. Our restriction is that we want our finitetemperature theory to generate the expected MOND profile. To construct such a Lagrangian requires first-principle knowledge of the microphysics of the superlfuid. Since we still do not have a fundamental description of the DM superfluid model, we proceed empirically. We suggest the following finite-temperature Lagrangian for the model: where the finite temperature effects are parametrized by a dimensionless constant β. When β → 0, we recover the T = 0 result; we are using the fiducial value β = 2. We included the interaction term so we could represent the entire action of the model that we are going to use next.

Halo profile
With the Lagrangian of the theory, we can evaluate the halo profile in the superfluid region and, after matching with an outer NFW profile, calculate the rotation curves of galaxies. And this is what we are going to do in this section: estimate the halo profile. This will be done in steps. First, we estimate the DM halo profile taking into account only the density coming from (136). Next, we include the baryons, by calculating the profile for the full action including interaction. We are going to use here the finite-temperature effective action (143), since in the case of the T = 0 the perturbations around a static background configuration suffer from a ghost-like instability. Although phenomenological, it retains the features of the initial superfluid Lagrangian and can give a more realistic description of the system.

DM halo profile
We can now calculate the density profile of the condenstate, in the superfluid region, assuming that we have only dark matter and no baryons for simplicity. This is the halo profile given by the different equation of state that the superfluid has: P ∝ n 3 , given by equation (137). This analysis is almost the same for the zero-temperature and finite temperature cases, with accounts for the replacement:Λ →Λ = Λ √ β − 1. Assuming hydrostatic equilibrium, for a static and spherically symmetric halo, the pressure and acceleration are related by: By making a change of variables ρ(r) = ρ 0 Ξ and r = ξ ρ 0 /(32πGΛ 2 m 6 )

1/2
, where ρ(r = 0) = ρ 0 , this equation reduces to the Lane-Emden equation (with n = 1/2), Choosing boundary conditions Ξ(0) = 1 and Ξ (0) = 0, we can numerically solve this equation. We can see from the change of variables that the size of the condensate and the central density are given by [203], where ξ 1 is the numerical simulation vanishes. From the numerics ξ 1 ∼ 2.75 and Ξ (ξ 1 ) ∼ −0.5 gives the following halo radius and central density: With that, we can determine the chemical potential µ = ρ 2 / 8Λ 2 m 5 . For m ∼ eV and Λ ∼ meV, we obtain realistic core sizes, which are of sizes that cover a big part of the halo, as we wanted. For this reason, we choose the fiducial values: For these values, we have a cored density profile with a condensate core of radius 158 kpc for M DM = 10 12 M . The condensate does not make the entire halo, but we expect that this condensed core is surrounded by a NFW profile. The central density obtained is smaller than the expected from CDM simulations, which is preferred by observations. In this way the DM superfluid offers a simple resolution to the cusp-core and the "too big to fail" problems. We will see these results in more details in Section 4.2.4.

Including baryons
Now, we derive the condensate profile in the presence of baryons. We expect that there is this extra acceleration due to the interaction to baryons. This comes from the dynamics of the phonon excitation φ given the the Lagrangian (143). We are going to assume a static, spherically symmetric approximation: θ = µt + φ(r). The equation of motion for the phonon is given by, whereμ ≡ µ − mΦ. If we ignore the homogeneous curls term, in the limit where (∇φ) 2 2mμ the solution is where a b is the Newtonian acceleration due to baryons only. Then acceleration mediated by φ is, for a 0 = α 3 Λ 2 /M pl , which is exactly the acceleration expected in the deep MOND regime, as showed in Section 2.2.4. In the regime (∇φ) 2 2mμ, we recover the Newtonian acceleration given by the baryons. So, in this model, the total acceleration is given by a b , a φ , and also a DM the Newtonian acceleration from the DM halo itself (obtained in the previous section), since we have DM in this model (different than MOND).

Halo profile algorithm
Having developed the theory of the superfluid DM above, now we want to evaluate the density profile of the DM halo and the rotation curves, and compare it with the data to make a first proof of concept of the model. To evaluate the rotation curve, we need to determine the circular velocity with respect to the radius.
As discussed in our model the galaxy contain a superfluid core in the central region of the galaxy surrounded by a NFW profile envelope. So in order to calculate these quantities for the galaxy we first need to evaluate them inside the superfluid core, and then at R = R NFW match the density and the pressure obtained for the superfluid ρ SF and P SF , to the ones given by the full NFW profile.
For that, we need to evaluate these quantities in the superfluid phase. In order to obtain the halo density profile, we need to determine the total mass of the halo M (r). The rotation curve is the circular velocity with respect to the radius, given by a = v 2 circ (r)/r where a = ∂Φ/∂r. So we need to determine the gravitational potential Φ in order to calculate the rotation curve and, also to determine M (r). The Poisson equation in the superfluid region is given by: MOND for rotation curves (but particle DM for lensing) Bars and spiral structure in galaxies MOND Interacting Galaxies Dynamical friction Absent in superfluid core Tidal dwarf galaxies Newtonian when outside of superfluid core Spheroidal Systems Star clusters MOND with EFE inside galaxy host core -Newton outside of core Dwarf Spheroidals MOND with EFE inside galaxy host core -MOND+DM outside of core Clusters of Galaxies Mostly particle DM (for both dynamics and lensing) Ultra-diffuse galaxies MOND without EFE outside of cluster core Galaxy-galaxy lensing Driven by DM enveloppe =⇒ not MOND Gravitational wave observations

As in General Relativity
The baryon density is given by the observations, while the superfluid density we can obtain from our theory by differentiating our Lagrangian (143) with respect to Φ,: where we can see that ρ SF = ρ SF (Φ, φ). So, in order to solve the Poisson equation, we need the equation for φ, which is given by its equation of motion (150). The system of equations we need to solve is given by (153) and (150), which can be very intricate to solve. One approximation that can be done to simplify this is to assume that baryon distribution is spherically symmetric (which we know it is not true, but used as a simplification). With that, the system can be solved numerically. This is done in [181]. After having this, this solution needs to be matched to the NFW profile that describes the outskirts of the halo. With that, it is possible to evaluate the density profile and the rotation curves of galaxies.

Observational consequences
In this section we will describe the main observational consequences of the superfluid DM. A summary of all the effects already worked out can be seen in Table 1. Since it is going to be used a lot in this section, we remind Landau's conditions for superfluidity is that the fluid velocity (v s ) is smaller than the superfluid sound speed c s , v s < c s .
-Galaxy rotation curves: In [181] the rotation curves of IC 2574, a low surface brightness galaxy, and UGC 2953, a high surface brightness galaxy, were numerically calculated using the method developed above, as a proof of concept of the galactic dynamics that the DM superfluid is able to reproduce. Since for the theoretical predictions a spherical baryonic distribution was assumed to simplify the calculations, and since this is far from the actual distribution of baryons, in these calculations a hybrid method, mixing the results calculated with the spherical distribution, was implemented. In this, the acceleration is corrected for the actual distribution leading to: where a b,real is the acceleration computed from Poisson's equation for a non-spherical baryon distributions; a DM the Newtonian acceleration from the DM halo using spherical baryon distribution; and a phonon from (152) sourced by a b,real , but with Newtonian potential from the spherical case. Although calculated in the hybrid method, a phonon ∼ √ a 0 a N,b as expected in MOND regime.
The fiducial parameters used for this numerical calculation were m = 1eV, (σ/m) = 0.01 cm 2 /g, which are optimal for having a superfluid core that encompasses the baryonic disk of the galaxy, while still within the bounds to agree with cluster observations; Λm 3 = 0.05 meV × eV 3 ; and α = 5.7. The rotation curves can be seen in Figure 13.
1. LSB galaxy: As pointed out before since these type of galaxies are DM dominated, the rotation curves from LSB are expected to have a slow raise before reaching the plateau region. As we can see in the left panel of Figure 13, our model reproduces the observed rotation curve for IC 2574, represented by the orange points, very precisely for the parameters chosen. The size of the superfluid core obtained for this galaxy is R SF ∼ 40 kpc, which here is represented by the NFW radius where the profile is matched with a NFW profile and has a close value to R T . Relative to R 200 ∼ 57 kpc for this galaxy, the superfluid core is relatively large encompassing 58% of the total DM mass of the halo. 2. HSB galaxy: The rotation curve features of HSB galaxies are known to be hard to be reproduced. We saw that MOND empirical theory is successful in reproducing those features. It is interesting to see if the superfluid DM model is also able to reproduce it. The rotation curve for UGC 2953 is shown in the right panel of Figure 13, using the same conventions as for the LSB. The radius obtained for the superfluid core in this case is R SF ∼ 79 kpc, which is small in comparison to R 200 ∼ 245 kpc. Only 24% of the total mass of DM is in the superfluid core. The difference from the LSB results is the red curve, where the total DM mass is set by the ΛCDM abundance matching value of M = 65 M . For the red curve, we get a bigger superfluid radius, R SF = 93 kpc, which is still significantly smaller than R 200 = 446 kpc. The rotation curves seem to fit the data well, showing a smaller value but still compatible with observations for the velocity in the point where the curve turns to flat. Also, the superfluid DM show a slight rise in the end of the rotation curve, which is compatible to the data but not existent in MOND. In general, it seems that the superfluid model reproduces the rotation curves of LSB and HSB galaxies. Also, the BRTF relation is satisfied, as expected. Of course, this calculation shows a proof of concept and the rotation curves of many more galaxies with different characteristics need to be fitted, also to help determine the parameters of the theory, which were chosen here. However, it is expected that the behavior of the rotation curve is similar to the fits shown above for different types of galaxies inside the superfluid core, only changing the size of the core depending on the galaxy.
-Dynamical Friction: From the very definition of superfluidity, flow without friction, we can expect that in these models in the inner regions of galaxies, where superfluidity emerges upon condensation, dynamical friction to be absent. This might lead to interesting astrophysical consequences and help understand some puzzles with CDM [164,202], while testing the DM superfluid model. One example of an observation that can be explained by this characteristic of superfluids is the velocity of galactic bars in spiral galaxies, which are expected to have been slowed down by dynamical friction, but are measured to be nearly constant which is consistent with no dynamical friction. Another interesting puzzle directly linked with dynamical friction is the Fornax globular clusters, as we already mentioned above. In the presence of a superfluid in the halo, given the absence of dynamical friction, these globular cluster should not necessary have merged with Fornax. The effect expected for the case of a superfluid is more pronounced than in the FDM model, for example. This shows that these type of system can offer an opportunity to test these ULDM model.
However to use these observations to test these models we need to really understand what is the behaviour of dynamical friction in the DM halo. For that a microscopical description of a superfluid theory with dissipation is necessary, and that is what is shown next.
The simple picture that in a superfluid there is no friction is a simplification. To fully understand how friction and dissipation in the superfluid takes place, one needs to work out the superfluid theories presented in Section 3 in the presence of dissipation. For that, one needs to describe the motion of an impurity, a particle moving in the superfluid represented by a real scalar field χ, and estimate the rate of phonon emission, which describes the dissipation of the superfluid. This was discussed in detail in [207]. The regime of validity of such a theory in the presence of dissipation is discussed in [133], and it is an extension to the discussion presented in Section 3, where we assumed the limit without dissipation. The theory with dissipation was worked out for the simplest superfluid example, the interacting BEC described by the microscopic Lagrangian of a self-interacting complex field (45), with the presence of an impurity, We are going to work with the Lagrangian which actually describes the SIBEC model, but that is the simplest model to understand superfluidity. We are interested, as before, in the non-relativistic case and low energy regimes. To study the dissipation of phonons, we perturb this Lagragian and work with the linear theory. The process that we want to study is the dissipation of the superfluid radiating phonons caused by the motion of the impurity, with Φ, the Newtonian potential behaving as the mediator of this process. This can be described by the process (at first approximation) χ → χ + π, and the rate of this process can be computed.
As discussed in reference [133], it is important to reach the correct result to consider the higher order derivatives of the phonon effective action, like we did in (45) that gives rise to the higher order k 4 term in the dispersion relation, together with the higher order terms involving the impurity field. With that it is possible to calculate the energy dissipation, is the rate of the process described above, with 'in' indicating the initial values before the collision and 'f' indicating final values. The initial momentum of the impurity is p in χ = (M + M V 2 /2, M V) and with k * given by k 2 * /2M + ω k * = k * V . With that, we can determine the friction force in the system: This shows that the friction force is not discontinuous, given having friction or no friction in the case of having a superfluid or not, as suggested by Landau, but actually it varies monotonically with the velocity. In the limit where V is equal to the sound speed, then the friction force vanishes, as expected for a superfluid. If we include gravity in this system, we have an extra term coming from the coupling to gravity which modifies the dispersion relation for the superfluid as shown in Sections 3.3 and 4.1.1.1 where there is an additional tachyonic mass term from the gravitational contribution, given by m 2 g ≡ 4πGρ 0 . The presence of this term modifies the Jeans scale and Jeans instability occurs when k > k J = 2m 2 c 2 s (−1 + 1 + (m 2 g /m 2 c 4 s )). In this case, as shown in [207], the force evolves monotonically with the velocity. However, it never reaches zero friction for subsonic velocities, when there is the superfluid, because of the Jeans instability. This shows that the dynamical friction in a superfluid is more complex than the simple dichotomy of absence or presence of friction if there is or not a superfluid. This is an active field of research and might lead to interesting observational consequences for the DM superfluid and SIBEC, which is the model worked out here.
-Galaxy Clusters: In a simple way, following the analyzes in Section 4.2.1, clusters have large dispersion velocities, and at large distances, of order of R 200 , Υ is going to be small and thermal equilibrium cannot be achieved. The DM in clusters is in the normal phase. However, as we saw for galaxies, in the central regions the density increases and thermal equilibrium may be achieved. In clusters only a very small amount can be in the superfluid state, since observations exclude that clusters are largely in the superfluid regime. We can then see the bounds in our mass in order to have a small amount of superfluid component in clusters that is not in tension with data. We assume that R T /R 200 0.1, which gives, using the relations from Sections 4.2.1 and 2.1: We can now repeat the analysis of thermal equilibrium done in Section 4.2.1. However, for such a small R T in comparison to the cluster size, we use the full NFW profile for the halo. This yields a constraint in the mass of the DM particles: This combined with the constraints from thermalization in galaxies gives the allowed range for the DM mass: From the tightest constraints from approximately 30 merging systems [208], σ/m 0.5 cm 2 /g. This value is in accordance with the one from Section 4.2.1, and from the constrain above it gives a DM mass between 1.5 eV m 2.4 eV. For DM superfluid in this mass range, we have condensation inside galaxies and the condensation in the interior of cluster happens for very small radius, appearing not to be in conflict with what is expected from observations. This constraint can be made broader by assuming a more realistic and not constant cross section. A quantitative analysis via numerical simulations would be ideal to check this result.
-Galaxy mergers: The behaviour of merging galaxies is an interesting question, given the superfluid nature, the absence of friction, proposed for the inner core of galaxies. In the absence of friction, it is expected that the merger would make the galaxies pass through each other without interacting. But the existence of these superfluid phases in these merging systems is going to depend on the comparison between the infall velocity for the merging galaxy and the sound speed of the phonon, given by the Landau criteria.
v infall c s -In this regime, the halos are driven out of equilibrium, so coherence of the condensate is broken and the halo will be in the normal phase. The merging process will proceed as in ΛCDM, where mergers are fast due to dynamical friction. Thermal equilibrium and condensation will be achieved in the merged halo after some time.
v infall c s -In the case of subsonic velocities, the DM halo is in the superfluid phase, and the superfluid cores will pass through each other with almost no dissipation. In this case, dynamical friction is reduced taking a much longer time to the system to merge, and possible multiple encounters. In our case, the phonons have sound speed c s = 2µ/m, and for the fiducial values adopted (149), c s ∼ 220 km/s for a 10 12 M halo. This needs to be compared with the infall velocities of galaxies of a merging system to see how the merger dynamics proceeds.

-Merging Clusters:
The Bullet cluster is a system of two merging cluster that was very well investigate observationally. It represents one of the best evidences of the existence of DM (and against alternatives like MOND). This is seen by a segregation in the position of the mass peak (highest concentration of total matter) given by lensing that probes all the matter content, and the one from X-ray measurements, which measures the baryonic matter. This is consistent with the CDM picture, where the DM in the merging processes due to its negligible interaction passes through almost without interaction, while the baryons are slowed down. This poses a problem for theories that do not have DM, like full MOND. By construction, and as it was seen above, in this model clusters do not develop a condensed core of cluster size and have most of the DM in the normal phase. However, the galaxies inside the clusters have a condensed core, and the cluster can develop smaller sized cores in its inner regions. Therefore, when clusters merge, the presence of a core or not also depends if the merger is subsonic or supersonic, obeying or not Landau's criteria like we saw for the merging galaxies. The outcome of the merging depends on the infall velocities, which determines if most of the DM is in the superfluid phase or in the normal phase in each of the merging clusters. If the infall velocities are subsonic, the superfluid cores are present, and most of the DM will be in the superfluid phase. Any collision between a cluster where DM behaves like a superfluid will follow without dissipation, with the clusters pass through each other without friction. Now, the DM in the normal phase presents self-interactions. Therefore the collision of two clusters in the normal phase would be slowed down due to these interaction. In the case of the Bullet cluster [178], in order to be consistent with observations, at least the sub-cluster must be in the superfluid phase. As we can see, the sound speed of the phonon for the sub-cluster (M sub ∼ 10 14 M ) is, for our fiducial values, c s,sub ∼ 1400 km/s, while for the main cluster (M main ∼ 10 15 M ) is c s,main ∼ 3500 km/s. The relative velocity between the clusters is ∼ 2700 km/s [209,210]. If we take this to be the infall velocity, we can see that the sub-cluster is in the superfluid phase, while the main cluster is in the normal phase. With that, the clusters will merge without dissipation and pass through each without friction, as it is expected from observations. For the Abell 520 train wreck [211,212,213,214], another merging cluster system, the DM superfluid model predicts a subsonic merger, with two peaks representing the superfluid component, compatible with the lensing map, and a peak during the normal component, coming from the X-ray luminosity peak. This shows that the DM superfluid framework can accommodate not only the dynamics on galactic scales, but also explain clusters and its merger events.
-Gravitational Lensing: In the case of the full MOND theory, or its relativistic completion TeVes, because of the absence of DM to be able to explain the relativistic regime makes necessary the introduction of a complicated non-linear term between the scalar field of the theory and baryons, which should also couple to a time-like vector field in order to give the correct gravitational potential to be able to explain gravitational lensing.
In the case of the DM superfluid, since the theory has DM, we have the superfluid component described by the phonon scalar field, and we have the normal component which provides the time-like vector field u µ . The gravitational potential is then sourced by both dark matter and baryons, as expected.
As we have that the superfluid core resides in the inner part of the galaxy, surrounded by an NFW envelope, gravitational lensing will come primarily from this NFW outter part.
Recently, the DM superfluid model was studied in the context of strong lensing [215].
-Gravitational waves: In the superfluid DM, different than in MOND, the superfluid core is locates in the inner regions of the halo and the outskirts of the halo have a NFW profile. So the gravitational lensing signal comes from this outer part of the halo and it behaves like in the case of GR+CDM. This means that photons and gravitons propagate at the speed of light travelling along the same geodesics. This is in agreement with the recent constraints from the gravitational waves from neutron stars merger GW170817 [216], which rule out relativistic completions of full MOND [217]. The implications for the gravitational waves in the case where the phonon has a non-vanishing sound speed was considered in [218], together with its observational effects in future GW experiments. Specifying the microphysics of the DM superfluid particle can also yield other signatures in the produced gravitational waves, like chirality, as done in [219].

Validity of the EFT
In this section we are going to scrutinize the validity of the EFT construction used, verifying the regimes where this leading order EFT is valid, and the regimes where the theory obeys the Landau criteria for superfluidity.

Higher order derivatives
First, we need to check if in our regime and for the parameters of the model, it is valid to ignore higher order terms in the EFT. As we saw above, the EFT is constructed by including all the terms which are invariant under shift symmetry. We retained only the first order contributions, given that we are working in the low-energy limit. Higher order terms involve more than one derivative per field. Higher order contributions to the quadratic Lagrangian for the phonon (138), can contain terms of the form: where ∂ → ∂ t or c s ∇, and φ c = Λ 1/2 m 3.4 µ −1/4 φ is the canonical variable. The scale that controls these higher order terms is given by Λ s = Λm 3/2 µ 3/2 1/4 , which we call the strong coupling scale. So, higher order corrections can be neglected when: This is the general condition for ignoring the higher order corrections in a EFT given by the Lagragian we described here. Given this, we an easily see that the approximation of ignoring these terms breaks for small sound speeds. However, specializing to the parameters of the DM superfluid model described here, and using the profile obtained in (147), which determines µ, the strong coupling scale is given by the DM superfluid model: For the fiducial parameters, Λ s ∼ meV. So, higher derivatives are suppressed if r 0.2mm, which is clearly satisfied on astrophysical scales.

Criteria for condensate coherence
An important criteria to verify the validity of the superfluid description we are using is to check if our superfluid obeys the Landau criteria. As we saw in Section 3.2 the criteria for the system to transports charge without dissipation, leading to the coherence of the BEC to be maintained, is that the velocity of the superfluid is smaller than the critical velocity: where the critical velocity must be non-vanishing. In our case, in Section 4.2.1, we already evaluated the conditions for DM to be condensed in the center of galaxies. We can estimate v c by using the halo mass density ρ = (2m) 3/2 mΛ |X| ∼ 2m 2 Λ √ κ, where we assumed MOND regime in the last equality and κ = mμ, which gives us: where M b is the baryon's mass. The superfluid velocity id given by v s = ∂ r φ/m ∼ √ κ/m, which yields: With that, using (166) and assuming spherical symmetry, we an determine thee radius where superfluidity can occur: We can see that this condition is satisfied in the central regions of galaxies, and we have coherence of the condensate and superfluidity in those scales.

Solar system
At solar system scales the bounds on deviation from standard Newtonian gravity are very tight, and these measurements do not allow deviations from the Newtonian dynamics. Full MOND is in tension with these bounds. However, the DM superfluid scenario fits well into the Solar system bounds. We can see that using the coherence bound for the condensate (166). For that we need to evaluate v s and v c .
The superfluid velocity in the vicinity of the Sun (M b = 1 M ) is given by: where r is the distance to the Sun and AU is the astronomical unit, the average distance between the Earth and the Sun. The critical velocity of the Milky Way galaxy (for M b = 3 × 10 11 M ) evaluated at our solar system (r ∼ 8 kpc) is: We can see that the coherence bound v s v MW c is obeyed for distance much larger than the solar system scales: Therefore, on distances like the solar system, DM is in the normal phase since the condensate loses its coherence, and obeys standard Newtonian gravity.

Relativistic Completion
As we saw in Section 3.3, the description of a superfluid is given by a weakly self-interacting field theory with global U(1) symmetry. The symmetry is spontaneously broken by the superfluid ground state of a system at chemical potential µ. In the previous section, where we defined this field theory for superfluids, we added a 2-body selfinteraction, g − 3|Ψ | 4 . This gives an equation of state P ∝ n 2 . As we saw in the previous section, the pressure that describes the interaction in the Madelung equations has the form of the pressure of a barotropic fluid. For a threebody interaction, the equation of state is given by P ∝ n 3 . For the DM superfluid, in order to reproduce MOND, we wanted to have a theory that gave P ∝ n 3 . So one might think that the DM superfluid could be described by the microscopic theory of an interacting BEC with three-body interaction. However, we are going to show now that this is in fact not the case, since those theories give a Lagrangian with different signs.

3-body interaction
Lets consider now like before that the self interacting theory with U(1) symmetry that gives us the superfluid has a 3-body interaction, instead of a 2-body one. The relativistic action of this theory is given by: where g 3 > 0 for stability. Like before, this theory conserves particle number. Since we are interested in the nonrelativistic (NR) theory, replacing Ψ = ψe imt and taking the NR limit gives us: With that, we can calculate the equation of motion, which gives us the Schrödinger's equation: The condensate is described by the background solution, at zero temperature: ψ 0 = √ 2m ve iµt , where µ = λv 4 /2m. The excitations are given by: where ρ is the gapless mode and φ is the Goldstone boson associated with the broken U(1). At low energies, we substitute this into (174) and integrate out the gapless mode: which is the action to leading order in the derivative expansion, with X = µ +φ − (∇φ) 2 /2m. This is very promising since the theory with a 3-body interaction gives a low-energy Lagrangian with the same exponent as the one we need for MOND and for the effective Lagrangian P (X) for the EFT of superfluids. However, it has the opposite sign, given that g 3 > 0! As we saw before, the limit where g 3 < 0 is unstable, and we cannot have condensation on all scales, superfluidity and MOND.
So, the expected description as a theory with 3-body processes is does not work for the DM superfluid model, where we want to recover MOND behavior in galaxies. It was phenomenologically proposed in [144] a relativistic Lagrangian that is able to reproduce our expected Lagrangian (136) in the non-relativistic regime, which is given by: The scale Λ c was introduced in order for the theory to admit Ψ = 0 vacuum. It is easy to see that this action reduces, in the non-relativistic limit and when Λ c |Ψ | 2 , this action gives (136). The condition for MOND, given by Λ c can be rewritten as |X| Λ 4 c /(2mΛ 2 ), which corresponds to: where a φ is the acceleration from the phonon that can be obtained from the action and given by a φ = α(Λ/M pl )φ . According to observations, the deep MOND regime is very accurate for ∼ a 0 /10, which poses a bound for Λ c . This theory presented here is a phenomenological relativistic version of the DM superfluid. However, it would be interesting to have a relativistic complete microscopic theory for the superfluids.

Cosmology
After working out the galactic behavior of the DM superfluid model, we need to work out its cosmological behavior. One compelling feature of this model is that at the same time it describes the small scale behavior, it also recovers the large scale successes of CDM. In this section we show how DM superfluid behaves cosmologically.
The first question we would like to answer is if DM is in the superfluid or normal phase cosmologically. We saw in Section 4.2.1 that the critical temperature of the DM superfluid is given by (130), and T /T c today is around 10 −2 for massive galaxies (M ∼ 10 12 M ). Cosmologically, the temperature in units of T c is much colder. We can estimate that by known that ultra-light candidates for DM, like the DM superfluid, are non-thermal relics. They can be generated, for example, through a vacuum displacement mechanism (see below for a definition of this mechanism) like the axion. So, the particles are created when H i ∼ m, which corresponds to a temperature for the photon-baryon plasma: which is around the weak scale! With that, we can rewrite the condition for thermalization (84), given that the velocity and density redshifts as v ∝ a −1 and ρ ∝ a −3 . At matter-radiation equality, we can write this condition as: where ρ eq ∼ 0.4 eV 4 and by using that v eq = v i a i /a eq ∼ eV/ mM pl is much smaller than one, while v i ∼ 1 since it was created deep into the radiation era. Since T /T c = (v/v c ) 2 , we have that: So, cosmologically, all the DM is in the superfluid state. As we saw, the cosmological temperatures are many orders of magnitude different than the temperatures on galaxies. For the EFT built for the superfluid to be valid on such a scales (of time, temperature), the parameters of the EFT Λ and α need to evolve with the temperature. This dependence is estimated in [144] by making some phenomenological statements for the theory to match both regimes.

Simulating ULDM models
We described above the main characteristics of the ULDM models. We showed how we expect the small scale structures to be suppressed in this model by computing quantities in the linear limit, and showed how this model presents a core solution for a simplified model of the halo. However, to study the formation of structures at different scales, and the formation of galaxies which are highly non-linear processes, one needs to resort to simulations. Cosmological simulations have been one of the biggest tools for the understanding of the non-linear formation and evolution of structures and galaxies in the past few years, modelling diverse scales and physical process that are present in those processes. Therefore, to better understand how different the structures are going to be in the ULDM models, together with modelling some structures that are exclusively present in these constructions like the presence of cores, interference and vortices, we need to resort to cosmological simulations (for a summary of the current FDM simulation, see [197,221].
The traditional simulation methods present in the literature to study the formation of structures, like N-body simulations or hydrodynamical simulations, cannot be readily applied to the case of ULDM since they do not take into account the wave nature of these models. As we saw above it is this wave nature that can lead to important observational consequences present in these models of DM.
There are two approaches to simulate the ULDM models. One solves the Schrdinger-Poisson system, composed by the Gross-Pitaevskii equation for a given ULDM model coupled to the Poisson equation; and the other is given by solving the hydrodynamical equivalent of the GP equation, the Madelung equations. Each of those approaches have advantages and disadvantages, so they can be considered complimentary.

Schrdinger-Poisson
Hydrodynamical-Madelung equations In the absence of interaction, we have the Fuzzy DM model, and in the presence of interaction, the terms marked in blue, we have the SIBEC model. The interactions can be a two-body or three-body interactions, like showed respectively in the last two terms of the GP equation above, or even higher order body collisions represented in the ellipsis, if this is allowed in the system. These describe fluids with different equation of state. In the Madelung equations, the interactions are represented by the interaction pressure term P int , which has the form of a polytropic fluid. This polytropic fluid can describe different fluids with different equations of state, that arise depending on the type of interaction. So one can simulate either one of those models using these system of equations 25 .
For the simulation that solve the Schrdinger-Poisson equations, also called wave simulations, there are many groups that are attempting to solve this system using different methods [222,223,224,225,226]. In general this approach is very good to describe the small scales, being able to resolve the small structures and taking into account the wave nature of the condensate, as we can see in the left panel of Figure 14. With that, in this approach predicts the correct and expected structures on small scales. Being able to resolve the smaller scales, this simulation can resolve and show the presence of cores in the halos, the granular interference structure in the condensate or the presence of vortices. However, perturbation theory in this approach is only valid for scales smaller than the de Broglie wavelength of the model, being then hard to explain accurately structures on larger scales. This approach is also demanding numerically, since it requires a more finely resolution in order to resolve scales of the order of the de Broglie wavelength. This makes these simulations to be much smaller in size than the fluid ones, not describing cosmological scales or being able to span many decades in redshift. The simulations that solve the Madelung equations, the fluid simulations, have the advantage of being able to use the already written and well explored hydrodynamical codes available in the literature. They can be implemented by adapting those known codes to the case of FDM or SIBEC. The difference from a normal fluid simulation is the presence of the quantum pressure term, ∇P QP = −n∇Q, where Q = ( 2 /2m)(∇ 2 √ ρ)/ √ ρ. This term is singular when the density is zero, and the quantum pressure is not well defined in this regime. This restriction translate into those simulation not being able to resolve the smallest scales, coarse graining through the granular structure or any other substructure expected in these models. This leads to fluid simulations predicting a more pronounced gravitational collapse leading to an enhancement bigger than expected in the power spectrum at small scales [221].
The advantage of the fluid approach, though, is not only being implemented using the already mature hydrodynamical codes, but also being able to run much larger simulations than in the wave case, since it is less computationally expensive. With this method cosmological size simulations are possible. Many groups have been working on simulating the FDM using the fluid approach, with some variations in the implementations of the codes and the solvers [227,228,229,230].
Many research groups are attempting to perform those simulations so we can better understand the behaviour of the FDM model and reveal possible smoking gun signatures of this model. Those simulations are crucial so we can understand and better search for these signatures on observations. For this it would be interesting to have the small and large scales of the simulations resolved. Since the fluid simulations are good to describe the large scales and the wave simulations the small scales, some groups are exploring the possibility of having hybrid simulations where both methods are considered for the scales they work better [221]. More work in the presence of interactions, describing the SIBEC, would also be welcomed, since this would reveal more about the superfluid nature of this DM scenario and the possible consequences of having a superfluid core in the inner regions of galaxies.
Simulations of ULDM models are a fast moving and essential field to study these models, and the current advances are very exciting, which makes us look forward to new results in the near future.
SIBEC -As we mentioned in Section 3, based on the ideas by Landau, the correct description of a finite temperature superfluid is using the two fluid model, where there is a normal and a superfluid components. The description used above for the SIBEC and DM superfluid extend the limit of usage of the zero temperature description of the superfluids when describing DM in a halo, which should be treated with a finite temperature approach. An effective extension to include finite-temperatures was attempted for the DM superfluid case. When simulating the superfluid model, one should use the equations for a superfluid in the two-fluid formalism. This is what it was done in [231]. They want to study structure formation in the a model where DM forms a superfluid. For this reason, they solve the hydrodynamical equations that describe this superfluid. The finitetemperature hydrodynamical equations for the simplest superfluid, the weakly interacting Bose gas in a trapped potential, were worked out in [232,233]. With that, the authors set the trapping potential to be the gravitational potential and re-write the hydrodynamical equations in an expanding background. In this way, the model describes the finite-temperature hydrodynamical equations for the SIBEC model. The authors test the two-body and threebody interactions for comparison.
One challenge of this approach is that the hydrodynamical equations do not include the Landau criteria in them. To include this explicitly one needs to add dissipative terms. Since this is not completely worked out analytically, to do that one has to make some assumptions and assume a form for these terms. A very big step in including dissipation in this theory was worked by [207,133] and shown above in the "dynamical friction" subsection of the DM superfluid. This is a work in progress. Until this is worked out in detail, and avoiding adding unknown dissipative terms, the authors chose to do this numerically imposing Landau's condition at every position in the simulation.
They numerically integrate the hydrodynamical equations from redshift z = 100 until today. The authors found from the simulation that, in SIBEC model, the growth of structure proceeds less efficiently than in CDM, as expected for the ULDM models, although more efficiently than expected, with the suppression more pronounced on small scales and at high temperatures. They also study the role of the interaction strength and of the equation of state.
This numerical simulation presents some limitations giving some of the assumptions and limitations inherited from the hydrodynamical formalism. For example, in this approach we cannot see the complete dissipation of the superfluid. Also, correctly adding the dissipation physics would be a big improvement in this description. However, this simulation represents a very important step towards simulating more realistic superfluids, which already showed to lead to interesting observations consequences.

Astrophysical probes
We presented above some cosmological and astrophysical observations that help to constraint the ULDM models. There has been a huge advance in the observations of the small scales, with new windows of observations being opened that can help determine the nature of DM. We present in this section some of these new advancements that are still being tested and being developed, but that promise to help testing the ULDM models.

Local Milky Way observables and stellar streams
Our galaxy, the Milky Way (MW), is our closest source of information about DM and it is a very good laboratory for studying its behaviour on small scales. Here we present some observations from the MW that promise to help us study and test different models of DM. We are in a very special era for observations of the MW and Local group with data coming from many current and future experiments like Gaia, Large Synoptic Survey Telescope (LSST), Prime Focus Spectrograph (PFS), WFIRST, Maunakea Spectroscopic Explorer (MSE), among others. Using the incredible new data from these observations promises to be revolutionary in the studies of the MW, and hopefully for helping discover the nature of DM.

Stellar streams
Stellar streams are a stream of stars orbiting a galaxy which are remnants of a tidally disrupted globular clusters or dwarf galaxies, that was torn apart by a more massive system. These streams are usually long and thin, and wrap around the disrupting galaxy. Streams are good dynamical probes since they are initially cold and very sensitive to the gravitational potential [234,235,236,237]. This means that the streams when encounter substrutures present in the halo of the galaxies can be influenced by it, causing dynamical heating, which are changes in the velocities in the stream, but that are very hard to detect. These encounters also cause disturbances in the morphology of the stream, with the formation of gaps which are underdensities caused by the sub-halos encountered [237]. Only a part of the sub-halos hosts baryonic matter and can be observed directly, so stellar streams offer the opportunity to detect dark halos invisible by the traditional methods (see also discussion about lensing below). These gaps contain information about the substructure and can be caused by clumps that are even less massive than a what is expected of a DM sub-halo, showing how sensitive the streams are for detecting substrucutres [238,239,240]. Since the stream extend for large distances in the galaxy and outside the galactic plane, streams might contain detailed information about the gravitational potential and its variations of large part of the halo of a galaxy.
Since different models of DM predict a different amount of sub-halos, it is argued that the stellar streams can be used to test DM models. In the case of the ULDM models, which suppress the formation if substructures, having a much bigger size of minimal subhalo allowed to be created in the galaxy halo. Models like WDM and SIDM also have a modified abundance of sub-halos in comparison with CDM. So there is the hope that these different models would imprint very distinct signals in the streams given their different substructure distribution Streams can also have gaps coming from baryonic substructures, so one needs to be very careful in the analyses not to overestimate the presence of sub-halos.
Up to now, 22 MW stellar streams are known, being the Sagitarius stream one the most important 26 . One of the streams that has been used to determine the presence of sub-halos is GD-1. The MW stellar stream GD-1 was discovered using SDSS maps [242], originated from a globular cluster, and it is seen as a 63 • long structure in the North Galactic region. Gaps of scales of approximately 10 • were found in this stream using SDSS and these were associated with the encounter with sub-halos by [243]. This was confirmed by [244] analyzing the stellar density perturbations from accurate measurements of the morphology [245] of this streams using Gaia data combined with photometry from Pan-STARRS [246]. They found that the data indicates that these perturbations should come from sub-halos, and that their abundance and masses are compatible with the expected from CDM from simulations. The error bars of these abundances and expected masses are still large, but these measurements inaugurate and opened the avenue for searches of substructures using stellar streams.
Novel experiments like PFS and LSST promise to measure the streams from the MW and Andromeda in more detail and with higher precision, with the goal of not only detecting the signals of substructures but to understand in more detail the characteristics of the stream, its formation and the properties of its progenitor.
However, a recent publication [247], revises these conclusions and indicate that the features found in GD-1 can be explained by simple epicyclic motion in a smooth Galactic potential as shown by their N-body simulation, without the need of the presence of a subhalo. Therefore, they conclude that the measurement from GD-1 show no evidence for the influence of DM. This shows that obtaining information about DM substructure in the streams might be more complicated than expected, with some degeneracies with other effects, and that more modelling and understanding of the influences of the MW in the streams are necessary. In this way, streams are a promising probe for DM, and, although this analysis should be taken more carefully, the next few years are guaranteed to be very exciting in this fast moving field.
MWs gravitational potential: Very quickly we would like to point out that another observable that can help us test DM. Knowing the shape, mass and distribution of the halo in the MW can gives us clues on the DM model since different DM models predict distinct shapes for the halos. To understand the halo we need information not only in the inner regions of the halo, but on many different scales up to the virial radius. The measurement of the position and velocity of satellite galaxies and globular cluster can give information for the dynamics in a good range of distance from the center of the MW. The distribution of satellite galaxies is already used to put constraints in the ULDM models (and other models of DM), as we saw above.
Experiments like Gaia, can give us very accurate data for scales much smaller than the virial radius. PFS galaxy archaeology survey will also measure stars in the galactic disk, complementing complementing Gaia's survey. LSST is expected to provide information from stellar tracers on scales close to R vir , being able to measure many new satellites that are fainter and more distant than the known today, extending the determined halo mass function by three orders of magnitude.
This is linked to the study of streams discussed above, since streams given their long range in the halo, can give us information on the gravitational potential of the halo for of even larger scales.

Other local MW observables -vertical dynamics
Our galaxy holds more information that can be used to probe DM models. There are several other local MW observables that imprint information from the underlying DM that forms the halo.
One example of one observation that can bring information about the different dynamics imprinted by models like MOND and the DM superfluid is presented in [248,249]. MOND and the DM superfluid model are constructed to give a very good fit to the rotation curves of galaxies and to explain the scaling relations, like the MDAR, by having a different dynamics on small scales. On top of that, the DM superfluid model reproduces the expected CDM behaviour on large scales. However, in these papers the success of those models in explaining the dynamics of galaxies is challenged by observations from the MW dynamics.
The radial dynamics is what it is constrained by those models, since it is the information that enters in the rotation curves of galaxies and the scaling relations. Now, if one consider also the vertical velocities in the vicinity of the Sun, the authors show that then it remains a challenge for those models to explain the rotation curves and these vertical velocity dispersion data. MOND and the DM superfluid, in slightly different ways, modify the dynamics on small scales, where for low-acceleration regions, a different acceleration than the Newtonian emerges in the system. In the DM superfluid this is caused by the presence of the phonons, as we saw above. This modified acceleration changes the dynamics in the radial and in the vertical directions. And in these articles they constraint both the radial and the vertical velocities dispersion for these models in comparison to CDM.
To constraint the radial dynamics, the very precise data from Gaia was used giving data of the circular velocity between R = [5,18] kpc. For the vertical dynamics it was used data of K-dwarfs from the SEGUE sub-survey from SDSS (Sloan Digital Sky Survey) were it was inferred the velocity dispersions for three mono-abundance stellar populations. With this data a Bayesian likelihood analysis was conducted. The authors found a discrepancy between the vertical acceleration predicted by MOND and DM superfluid model in comparison to the one inferred from the data, giving values that are around 15% larger than the ones inferred from observations and larger than the ones predicted by CDM. And for this reason they claim those models are not preferred as the DM model. This work establishes an important new observable that should now be taken into account when constructing DM models, their impact not only in the radial dynamics, but also in the vertical dynamics.
However, as some authors have pointed out this result should be taken carefully. In [251] it is pointed out that this is a 2σ discrepancy and that maybe the data used to infer the vertical velocities implies a local DM density that is the double of the one inferred from the radial dynamics. This would give a halo that is not spherical, like assumed in the analyses of [248].
This brings an important point about the data for the vertical dynamics of the MW that we would like to highlight. The dynamics of the MW is very complex with a rich accretion history. Although the MW, with its thin disk, is considered a very stable galaxy with no major recent dynamical interactions, many observations indicate a major accretion event with Sagittarius (Sgr) dwarf spheroidal. This is seen in the streams that wraps around the MW, and it is expected that this event strongly influenced the dynamics of the MW. In [250] they show a study of the MW's major accretion events with Sgr and the Large Magellanic Clouds (LMC). They show that this encounter with Sagitarius produces oscillations in the MW disk. These vertical perturbations are an influence of many passages of Sgr, as it falls into the MW. The influence of the LMC also changes the dynamics, but much less than Sgr, with Sgr being the main influence in the dynamics of our galaxy in its recent evolution. This shows that the dynamics of the MW is complex and has many unknowns. These vertical oscillations are also found in recent observations both in the Galactic disk and around it, close and far away the solar neighbourhood (see [250] for a summary and discussion of those). Therefore, even if the explanation for these oscillations and other features in the Galactic disk is not solely the one above, this rich dynamics of the MW has to be considered when using these observations, specially to constrain models.
In this way there is the possibility that the vertical dynamics data used in [248,249] has the influence of these oscillations and cannot be used at face value in order to constraint the dynamics of DM. Or, maybe this complexity of features has to be considered as larger error bars in the inferred vertical dispersion. Therefore, vertical velocities are a new and important observable to test DM models, but it necessary to understand if this information can be disentangled from the complex dynamics from the accretion history of the MW. When this degeneracy is resolved this will be an important criteria for DM models to satisfy in the MW.
In [252] the authors also point out that the model of the DM superfluid used is too simplistic and that the DM superfluid interacts with baryons which might also change the vertical accelerations.

Substructure
Different models of DM predict different substructures, from its abundance to the minimum mass for the possible substructures formed, like the ULDM showed above. Therefore, the observables that probe substructures represent one of the most important tests for DM on small scales.
We already presented above a probe of substructures present in the halo, the stellar streams. This probe seems to be very sensitive to DM substructures and has the potential in the future to help test DM models.
Experiments like LSST will also make a more complete search for ultra-faint satellites, which is an important search for probing small substructures and that can be used for testing any DM models that suppress the presence of substructures below a certain scale.
However, there is another way of probing more generally and directly the presence of substructures, if they are luminous or not, which is gravitational lensing. Gravitational lensing is distortion of light from objects by the presence of a gravitational potential. One of the main observations used to search for substructures is strong gravitational lensing of quasar. Lensed quasar can present multiple images, arcs or even Einstein rings. The presence of substructures modifies the lensed images of quasar changing its morphology and flux ratios, in a way that the substructure can be mapped, as done in [254,255,256,257], including distortions from sub-halos in the line of sight [263]. There are many efforts to probe DM substructures from strong lensing using different frameworks [258,259,260,261,262], and machine learning techniques [264,265] (for a list of other machine learning approach see [264]).
As we saw above, generally all the ULDM models suppress the formation of sub-halos of a certain mass and size. Other types of substructures can also appear in these models, like vortices or even a substructures with different shapes like dark disks [253] in the case of a superfluid DM, which are not expected in CDM and are very distinct signatures of those alternative models. Strong gravitational lensing is an important technique to probe these substructures, and one of the techniques that is going to improve the most with the upcoming experiments.

21 cm cosmology
As discussed above, the ULDM models give a suppression of the matter power spectrum on small scales. Those scales can only be marginally probed by the known cosmological probes like CMB, LSS, cluster abundance, Ly-α forest, with those measurements putting lower bounds in the amount of suppression those models can have. Those measurements can only constraint the structures on scales k ≤ 10 Mpc −1 , not being able to probe the smaller scales. One new window of observation, the 21 cm line from neutral hydrogen (HI), promises to allow us to probe the matter power spectrum in much smaller scales, k ≥ 10 Mpc −1 , on scales even smaller than the Ly-α forest. This is possible since neutral hydrogen is only present in the universe, after reionization, inside dense clouds in damped Ly-α systems, which are small objects (k ∼ 10 2 Mpc −1 ). In this way the 21 cm HI signal is a biased tracer of the galaxies, and consequently of the underlying matter distribution in such small scales. Therefore, measuring the global 21 cm HI signal together with its fluctuations can gives information about the still largely unconstrained small scale matter power spectrum.
However, at those scales these system are dominated by astrophysical process, making it difficult to disentangle the behaviour of DM from those processes. To obtain cosmological information from these measurements is a challenge. Some recent studies [266] forecast that the matter power spectrum can be measured from the global 21 cm HI signal for an experiment with parameters close to EDGES [267], with a precision of O(10%) integrated over the scales k = (40 − 80) Mpc −1 , after imposing priors on the astrophysical effects like star formation rate and feedback amplitude. They also parametrize the effect of foregrounds, like the Galactic foreground, that plagues all 21 cm experiments and might represent a huge limitations for them if not well mitigated. Detecting the 21 cm HI fluctuations is a much harder task. Large interferometer experiments like HERA (Hydrogen Epoch of Reionization Array) [268], LOFAR (LOw-Frequency ARray) [269], LWA (Long Wavelength Array) [270], and SKA (Square-Kilometer Array) [271] have the goal of measuring this signal from the epoch of reionization (EoR) and also late times. For a HERA-like experiment it was found that the matter power spectrum can be probed with an accuracy of O(10%) integrated over the scales k = (40 − 60) Mpc −1 and k = (60 − 80) Mpc −1 . The measurements of the fluctuations probe the evolution of the matter density tomographically, carrying more information about the scale and redshift dependency of the HI signal, bringing more information on the astrophysical processes. This makes the parameters of these astrophysical processes to be better disentangled from the HI signal, allowing this constraint in the power spectrum be less dependent on the astrophysics in these regions.
Specifically for the case of ULDM models, forecasts using 21 cm HI signal were made in [272,273,274], and they specialize in the FDM model. Reference [272] studies how the recent EDGES measurement of the global 21 cm HI signal [267] can constraint the FDM model. The global signal is the average radio signal from 21 cm redshifted emission from z ∼ 15 − 20 in the case of EDGES. This measurement showed an absorption profile that had an amplitude two times bigger than the expected. This higher amplitude indicates that already at redshift z = 20, there was significant star formation, which leads to a also significant Lyα background. This fact shows that the smallest structures, and consequently the power spectrum on small scales, cannot be largely suppressed. This puts constraints in models of DM that have the feature of suppressing the small scale structures, like the FDM (or any ULDM model). This measurement alone is capable of putting a challenging constraint in the mass of the FDM particle: m ≥ 5 × 10 −21 . This bound is challenging for the FDM models, since, as we saw above, this value is not small enough for this model to address the small scale problems of ΛCDM. To obtain this constraint some assumptions on the star formation, and on the halo mass profile had to be made. Given the importance of this result for the FDM models, the bounds obtained from this data should be explored further, as well as the future data from 21 cm signal.
To better understand how the HI signal is affected by the FDM model, in reference [273] they study of the impact of FDM models in the 21 cm HI signal from the cosmic down and EoR analytically, together with some forecasts for future experiments. They use an analytic model which take into account the Ly-α coupling, X-ray heating and ionization to study the 21 cm in a ΛFDM cosmology. They find that suppression of structure from FDM models, which makes small sub-halos absent in this model, has the effect of postponing the formation of sources and the reionization of neutral hydrogen. This delay changes the global 21 cm signal showing a smaller absorption feature than expected from ΛCDM. The amount of suppression allowed considering the results from the EDGES experiment puts a lower bound on the mass of the FDM model, with m ≥ 6 × 10 −22 eV, which might already be considered challenging for the FDM to solve the small-scale problems, but marginally. They also show the potential of a SKA-like and LOFAR experiments to test these models in the future. This shows us that future experiments will be able to confirm the important bound of the FDM model imposed by the EDGES.
In [274] they analyze the impact of measurements of the 21 cm forest, which alternatively from the tomographic and power spectrum techniques to use 21 cm HI signal, proposes to use the 21 cm narrow absorption features from the intergalactic medium cause by high-z loud radio sources or collapsed objects, like minihalos. The 21 cm forest is expected to be measured by SKA. In this reference they show that the impact of this measurement can also constraint the mass of the FDM model, and that this is degenerate with the fraction of FDM that composes the DM.
The possibility of 21 cm HI signal to constrain alternative models to CDM, like WDM for example, was studied in many references [275,276,277,278,279]. They show how the signal changes for different DM models and also show how measurements like the one from EDGES can put constraints in the mass of WDM.
The study of the capabilities of 21 cm experiments to give us cosmological information is an active field of study and the references above are just some examples of those efforts. These studies give us hope that, in the near future, this new window of observation will allow us to probe the still unconstrained small scales, helping elucidate the nature of the DM component.

Black-hole superradiance
Until now we explored ways of probing DM by the gravitational effect they imprint in the structure and substructure of our universe. We explored in detail the astrophysical tests of DM which happen in environments where we are in the weak field regime, but that are dominated by baryonic effects and complex non-linear physics, but that can still give us hints of the nature of DM. We present here a very different way of probing specifically ultra-light fields in strong field environments, far from the linear cosmological scales.
Ultra-light particles can be largely produced around spinning BHs, a process called BH superradiance [281]. When a BH rotates faster than the angular phase velocity of an incoming wave, it amplifies the energy and angular momentum of the field in its vicinity. This superradiance effect [282,283,284,285] is a natural mechanism to create clouds of ultralight bosons around Kerr BHs (see [281] for a review). For ultra-light particles with Compton wavelengths of order or larger than the BHs gravitational radius, they will be efficiently produced by the superradiance, forming a large 'cloud' around the BH. This cloud is a condensate of ultra-light particles created through this instability carrying up to 10% of the BH's mass and angular momentum, diminishing the initial rotation of the BH [286,287,288,289].
In the non-relativistic limit, the eigenfunctions of the system are determined by a Schrdinger-like equation and the whole set up is sometimes referred to as a gravitational atom. Superradiance instability depends on the spin of the BH, the mass of the BH and the mass of the ultra-light particle created, where the modes are co-rotating with the BH. Depending on the mass of the BH, from stellar mass BHs to supermassive BHs spanning masses from a few to billions of solar masses, the ultra-light bosons produced through this mechanism can have masses from 10 −20 to 10 −10 eV [290,291]. Superradiance can also occur for ultra-light vector fields. In this case particles of mass 10 −14 eV to 10 −11 eV can be created by stellar mass BHs, and 10 −20 eV to 10 −17 eV for supermassive BHs [292].
This cloud emits GW, which allow us to probe its presence around BHs. This GW signal could possibly be detected by experiments like LIGO, when coming from stellar BH clouds, and LISA, from supermassive BH clouds [290,291,292]. The signature of this GW will depend if the cloud is made of real particles, which creates a non-axisymmetric cloud, or complex scalar particles, where the cloud is axisymmetric and the emission of GWs is suppressed. For the cloud made of real bosons given the non-axisymmetric configuration, it can emit GWs when the ultra-light bosons interact with gravitons, or when gravitons change levels, emitting monochromatic GWs [291,298]. The cloud can also collapse if there is an attractive interaction between those UL particles, emitting GWs in this process.
If the spinning BH and its respective cloud is in a binary BH system, the GW signal is modified due to the presence of the companion presenting a richer GW phenomenology [293,294,295]. The evolution and GW signature of the cloud is modified by the presence of a second BH where the waveform and the amplitude of the signal can be modified and even vanish given the resonant transitions between the growing and decaying modes of the cloud. Sharp features in the GW waveform appear, offering a window to probe the signal from these UL particles using GW experiments. This is an active field of research with the modelling of these effects, the study of back reaction and the observational signatures still ongoing topics of research. Detecting the GWs coming from these clouds would give us the opportunity to probe ultra-light particles in a range that is very interesting for DM. In this way this observational signature is very relevant to the ULDM models.
Recently, the possibility of producing part of the DM of the universe by the evaporation of primordial BHs (PBH) was studied in [297]. In the process of their evaporation, PBHs of mass around 0.1 g to 10 9 g which evaporate before BBN, decay into DM particles with masses around MeV − GeV. Although those masses are higher than the ones we are interested in this review, this mechanism serves as another illustration of a different way of creating DM particles from these strong gravity objects that is being studied in the literature.
As we saw here, those busy and complex environments like galaxies and BHs might still offer, despite their modelling complexity, new ways of probing de fundamental properties of DM. For a review on the same lines for probing fundamental physics in those environments see [299].

Probing the quantum nature of ULDM
We have shown above the astrophysical and cosmological consequences of DM being described by a BEC or a superfluid inside the halos of galaxies. From changing the behaviour of LSS on small scales, to changes related the halos having a minimum mass in those models, and the different profile the wave nature of ULDM leaves in the inner parts of halos, we have been investigating the gravitational consequences of the ULDM models. Some of those predictions might be degenerate with other DM models, like the WDM and SIDM, and even with possible consequences from unaccounted baryonic physics.
However, given the quantum nature of the ULDM, these models present some predictions that are a direct consequence of this wave nature and that are completely distinct from any other DM model. We discuss here vortices and interference patterns. These are distinct effects that can appear in galactic scales, and can also be connected in the formation of halos.
The detection of any of these effects would be a direct evidence of those models of ULDM and a smoking gun for the wave nature of DM. For this reason it is interesting to study how these effects arise in each of those classes of ULDM models and understand if they yield observable consequences that allow us to test them.

Vortices
Until now we neglected the fact that the halo might be rotating. As we saw in Section 3.5, rotating superfluids have an interesting new phenomenology, the appearance of a lattice of quantized vortices that allow for the rotation of this irrotational fluid. This is a purely quantum mechanical phenomena resulting in quantized vortices being produced 27 . This new phenomenology arising from rotating BEC halos might lead to observable astrophysical consequences that will represent a direct probe of the wave nature of these DM models. The hope is that measuring this unique signature from these models, will make it possible to distinguish this class of models from other alternative DM models which might present some signatures that are degenerate with the ULDM models.
Disk galaxies are one of the most common galaxies in our universe and those are rotationally supported systems. Therefore, the DM halos from those systems are also expected to be rotating. Galaxy halos acquire angular momentum in their formation via tidal torques coming from the neighbouring large-scale structures [300]. This angular momentum is conserved after those halos virialize, which leads to the rotating supported disks in galaxies. This process is still being fully understood, with N-body and hydro-simulations showing that the halos of galaxies are expected to have angular momentum. From CDM N-body simulations, the angular momentum obtained, represented by the dimensionless parameter λ = L|E| 1/2 /GM 5/2 where L is the angular momentum and E the total energy, are on the range λ ∈ [0.01, 0.1] [301]. Therefore, when considering realistic halos of DM, one needs to consider rotation, and if DM is made of ultra-light particles, this can lead to the formation of vortices.
Here we show the effects of the halos being rotating in the ULDM, an effect that is not as explored in the literature as expected, but that might present a decisive observational signature for these ULDM models. We saw in Section 3 above there are conditions for the formation of these vortices, depending on the angular momentum of the rotating halo. The presence of vortices in the DM halo can alter some properties of the halo, like the mass distribution or the presence of those substructures in parts of the halo that might lead to observable signatures. Detectability of these vortices or of the effects caused by their presence, and the observational technique used to probe them, will depend on their abundance and size, which should be studied solving the Schrdinger-Poisson system. Those present some theoretical solutions and estimates of the properties and formation of these vortices, but mainly vortices have to be studied numerically, with wave simulations like described above. We present here some of those that investigate this in the SIBEC model, and a rough estimation for the DM superfluid.
During the final preparation of this review a paper studying the possibility of the formation of vortices in gravitationally bounded BECs appeared [118]. In this paper they study the possibility of formation of a vortices in the FDM model and of other topological defects with different topologies. They also present results of simulations of this system for rotating halos and the possible observational signatures as consequence of the existence of such vortices.
Other topological defects like the creation of strings on the condensate in the DM halo in the SIBEC is investigated in [114].
Self-Interacting BEC -There is a small amount of studies in the literature that investigate the presence of vortices [145,303,302] (and recently [118]). Given that to investigate vortices one needs to study the Schrdinger-Poisson equations, the main model where those were studied is the simplest superfluid model given by the SIBEC model. But even in these simplified models it is still hard to find solutions to the GP equation.
The study in [302] discusses the possibility of the formation of a lattice of vortices in the halo of a galaxy with the same parameters as Andromeda. This aims to show a order of magnitude estimate for a known galaxy, the amount of vortices formed and put bounds in the mass of the SIBEC particle, and in its interaction. Assuming that a galaxy like M31 has M ∼ 10 12 M , average radius R ∼ 100 kpc, average density ρ avg ∼ 10 −23 kg/m 3 and angular velocity Ω ∼ 10 −16 rad/s. If vortices are created they cannot be larger than the size of the condensate core, R v ≤ R bec , the SIBEC must have a mass m 10 −24 eV, with radius R v ∼ 10 21 m ≈ 30 kpc. Therefore, in a halo like the halo of M31, there might be from 1 to 100 vortices. For the simplified case of having only one vortex in the halo, the GP equation can be solved in the Thomas-Fermi limit, ignoring the QP term, determining the structure of this cylindrical vortex core, and the critical velocity, showing that Andromeda galaxy could have formed a vortex given its angular velocity.
Reference [303], focus in a simplified study of the consequences for virialization and on the rotation curves if a lattice of vortices is present in the FDM model. Depending on the choice of mass and interaction, this vortices can lead to oscillations in the rotation curve of galaxies, which the authors claim resemble what is observed in rotation curve of spiral galaxies. Now, in reference [145], the authors try to study the vortex solutions in more detail, solving the GP for some assumed halo profiles that are more realistic than in the previous work (see also [175] for a review of vortices in SIBEC).
First, it is interesting to repeat something we presented in Section 3.5. The presence of a rotation in a condensate in a spherical halo is that without rotation, the BEC wavefunction is real and positive, however the angular velocity induces a superfluid current, making the wavefunction to be complex, ψ = |φ|e iS(t,r) , with a phase that gives a velocity flow, v = ∇S/m, and as we saw in Section 3.5, the fluid velocity is v = v − Ω × v. Therefore, a rotating BEC is a superfluid. As shown in Table 1 from [145], FDM models, where there is no interactions, do not form vortices. If the interaction is attractive, vortices are also not formed. For this reason we work in the case of SIBEC model with repulsive interaction.
The condition for the formation of vortices is Ω > Ω c or equivalently L > L qm . This imposes a bound in the mass of the SIBEC particles. They work this out and the solutions for the GP equations for two halo models. First, in halo model A, they assume a simplified halo model where the density and potential is given by a homogeneous Maclaurin spheroid, which gives a known form for the gravitational potential inside halos. This halo is not irrotational, having L L qm . With that simplified gravitational potential one can calculate the characteristics of a virialized rotating halo. This allows to put bounds in the mass and the interaction of the SIBEC particle, assuming the critical case where (L/L qm ) c and λ = 0.05 is the average of the dimensionless rotation parameter: m/m H ≥ (m/m H ) c ≈ 50 and g/g H ≥ (g/g H ) c ≈ 2550.
For a more realistic halo which allows to take into account the compressibility of the fluid, the halo is considered as (n = 1)-polytropic Riemann-S ellipsoids, halo model B, which is irrotational before forming halos L = L qm . This shows that if L < L qm , there is no formation of vortices and the rotating BEC can be described by halo model B. And for L L qm , there is the formation of, at least one vortex if the quantities are equal to the critical quantities or more in case they exceed these values, with the halo being described by halo model A. The case where L L qm , the halo is described by the halo model B again, with a single vortex in the center of the halo.
This study shows the conditions for the formation of the vortices and shows the complexity that can arise in the presence of rotations and of a vortex lattice.
These studies show important characteristics for the formation of vortices in the halo of galaxies in the SIBEC model. However, to fully study the presence and consequences of a rotating halo in those model one needs to perform wave simulations. One study that takes that in consideration is [118].
DM superfluid -As we saw above, to calculate the abundance and properties of the vortices, it is necessary to solve the equations of motion of the superfluid coupled to the Poisson equation. For the DM superfluid model, there is still no numerical study of the solution of the equations equivalent to the GP in the presence of gravity in order to understand the formation of vortices in this model. For this reason, here we present a dimensional analysis and order of magnitude estimation of the presence of vortices in the DM superlfuid context.
Vortices are formed when the angular velocity of the superfluid is larger than the critical velocity: Ω Sf Ω c , with the critical angular velocity given by (82). For R ∼ 100 kpc and m ∼ eV, we can see that Ω c ∼ 10 −41 s −1 (neglecting the logarithmic factor). This is much smaller than the rotation velocity of the halo: Ω ∼ λ √ Gρ ∼ 10 −18 s −1 , using a halo mass density ρ ∼ 10 −25 g/cm 3 . So, there will be the production of vortices in the halos of the galaxies in the DM superfluid model And this production seem to be very numerous.
We can estimate the number of vortices in the halo N v = Ω/Ω c ∼ 10 23 , with a core radius, given by the healing length, ξ = 1/(mc s ) ∼ mm (assumed a MW type galaxy and fiducial values). Although highly numbered, these vortices are small. It is still unclear if it is possible to detect those vortices via, for example, gravitational lensing or any other effect they might have in the galaxy.

Interference fringes
Another interesting effect that comes from the wave nature of DM superfluid model is interference. Since the condensate is described by a coherent wave function with the density if give by |ψ| 2 , in the formation of galaxies, interference patterns are expected to form.
As we discussed in the numerical simulations section, the interference patterns can only be seen in wave simulations, since when there is destructive interference, the quantum pressure term is not well defined and this is not present in the fluid simulations. We show in Figure 15, the result of two wave simulations from different groups that show the interference fringes that appear in the filaments of FDM structure. These patterns appear on scales of order of the de Broglie wavelength. General features like caustics, which are density singularities that usually appear in CDM and even WDM simulations, are regularized in FDM models due to the uncertainty principle, and do not appear in FDM simulations.  [304] showing the density distribution of FDM at redshift z = 0.1 at different scales. The code used describes the wave-like FDM using an adaptive-mesh-refinement (AMR) scheme. Right panel: Figure from [305] that shows the simulations of the slices of DM through a filament in the DM distribution at redshift z = 5.5. This simulation was made by modifying the magnetohydrodynamics code AREPO for FDM model. These figures clearly show the interference pattern that the FDM model imprints in the halos which is very distinct from CDM or even other DM models like WDM. These interference fringes that are present at the scales of the de Broglie wavelength. On top of that, in both stimulations we can see the soliton cores formed in the halos, also characteristic of the FDM model.
Another effect that could generate interference patterns is the collision of subsonic galaxies. Since those galaxies maintain the coherence of their condensate core described by the coherent wave function, the collision between these cores leaves an interference pattern [309]. If this interference fringes could be observed, this would be still another form of probing these models.
Interference patterns are going to be created also in the merging processes of dark/bright solitons [306,307]. Some authors suggest that this effect could be linked to the shells seen around elliptical galaxies [308].

ULDM as dark energy
Ultra-light fields can also behave as dark energy (DE), depending on their mass and on the different theory they are applied. We can see this in two different cases where the ultra-light field can be used to explain the acceleration of the universe.
Fuzzy DM: In the case of fuzzy DM, where we have a (non-interacting) ultra-light scalar field in FRW universe, the field can behave as dark matter, early dark energy or dark energy depending on the mass of this field. The behaviour of the field depends on how its mass is related to the Hubble parameter. At early times, the ultra-light field has m H. In this regime, the field is almost frozen and behaves as DE with w ≈ −1. As m ∼ 3H, the field starts to coherently oscillate around the potential minimum and to behave as DM, where the equations of state averaged over the oscillations approaches zero. Depending on the mass, the oscillating DM phase can happen at different times of the evolution of the universe. If the oscillating phase happens after radiation-matter domination, the ultra-light field behaves like DM. From [23], we see this happens for m 10 −27 eV. This bound comes from observations from the CMB and LSS galaxy clustering, and in this regime the ultra-light field can oscillates before the present day and the field redshifts as DM. More generally, for m 10 −32 eV we can still have the ultra-light field behaving as DM, since the field starts oscillating before the present time, but in this case the ultra-light field can only make a fraction of the DM. For m ∼ 10 −33 eV ≈ H 0 , the field behaves like a quintessence field, if in the presence of a potential, and it can be the responsible for the late time acceleration. Since the field is virtually frozen until the present time, or slowly-rolling the potential, with almost constant density, and it behaves very closely to a cosmological constant. With observations from CMB, we can see that in this case the ultra-light fields have a maximum bound on the energy density compatible with the expected amount of DE in the universe. For masses around 10 −27 eV m 10 −32 eV, the ultra-light field behaves like DE earlier than what it is expected, and can be thought as an early DE component.
Superfluid DM -Unified superfluid dark sector: There is another way to explain the late time acceleration using these ultra-light fields, where the acceleration is not given by this behaving like a quintessence field. This can be done in the context of the DM superfluid, where the dark energy behaviour is yet another manifestation of the same superfluid that emerges at cosmological scales at late times, as presented in [310,311]. In the previous case and in the case of quintessence, a ultra-light field is a component with a very small mass m ≈ 3H 0 , that dominates around the present times and drives the acceleration. Differently, in the framework we present here, the dark sector is composed only by DM described by a superfluid, without DE. The late acceleration emerges from the dynamics of this superfluid, and we have a unified model for the behaviour of DM and DE. We can see how this arises from the model described below. BECs and superfluids in the laboratory are usually made of atomic species. After the discovery of BEC and superfluidity, one big evolution in the study of these systems was to study mixtures of condensates where atoms that compose the fluid that condenses might be at different atomic configurations. This allowed researchers to explore the richness of the internal structure of the atoms that compose the superfluid, describes a more realistic system, where atomic transitions are allowed to happen, and also the different dynamics that appears in systems where more than one species of condensate and superfluid is present. In the model we present here, we assume that the superfluid is composed by two different species, which can be represented by the ground (Ψ 1 ) and excited (Ψ 2 ) state of the dark atom that composes the superfluid. These species interact through a Josephson interaction [314,316] 28 , which is a contact interaction between the components of the superfluid that has the simple form L int ∝ −(Ψ * 1 Ψ 2 + Ψ 1 Ψ * 2 )/|Ψ 1 ||Ψ 2 |. This interaction leads to an oscillatory potential for the low-energy Lagrangian of the phonons, L = P 1 (X 1 ) + P 2 (X 2 ) − (1 − 2Φ) V (θ 1 , θ 2 ) , with V (θ 1 , θ 2 ) = M 4 [1 + cos(θ 2 − θ 1 + ∆E t)] .
From the form of the interaction term, the oscillatory potential for the phonons is given by a cosine potential, where M is the explicit symmetry breaking scale coming from the interaction Lagrangian that breaks softly the shift symmetry of the phonon action, and that has to be of the order of M 4 = 2M 2 pl H 2 0 ≈ meV, in order to drive the late time acceleration. The parameter ∆E is the energy gap between the two species, between the ground and first excited state of the component of the superfluid. Considering the approximation ∆E m i , we can see that the superfluid has two distinguished behaviours: one degree of freedom that behaves like dust and one that evolves under the influence of the potential, like what is expected from a field that behaves like DE. The cosine potential is similar to the pNGB models of DE [312,313], and it is a special potential for explaining DE since it only softly breaks the shift-symmetry, and the flatness of this potential is still approximately protected against radiative corrections, which is one big problem in quintessence models of DE. The late time acceleration behaviour from this DM superfluid can be seen in the evolution of the Hubble parameter in the NR limit, From the figure above we can see that we have a decelerated evolution, following the behaviour of DM, followed by a period of accelerated expansion at present time. Therefore, this model behaves like what is expected by DE, even without the presence of a specific component responsible for the acceleration, and being a model of DM alone. At future times, this models deviates a lot from the predictions of ΛCDM, as the cosine potential becomes important.
Although the evolution in this model is very close to ΛCDM, this model presents distinct predictions. This can be seen by computing the perturbations in this model. One of those predictions is the growth factor, that in this model deviates from the ΛCDM one by around 10% at present times. Future galaxy surveys might be able to test this deviation.

Summary
In this review we studied an alternative class of DM model, the ultra-light DM. These models have been receiving a lot of attention in the literature nowadays given their interesting property of forming a BEC or superfluid on galactic scales. In this review we aimed to give a summary of the models of ULDM, suggesting for the first time a classification into three categories according to the ways they achieve condensation: the fuzzy DM, self-interacting BEC and the DM superfluid. Their different descriptions lead to different phenomenologies and observational effects that can be used to test these DM paradigms. We had also the goal to give a snapshot of the field as it stands at present. We highlight how important the observations on small scales are for helping to determine the nature of DM, and how the observations of the small scales, galaxies and our MW have been advancing very fast in the last few years. With current and new experiments like Gaia, LSST, PFS, HERA, SKA, just to cite a few, and the new and exciting probes like stellar streams, 21cm cosmology, MW observables, BH superradiance, the next few years promise to revolutionize the tests of DM on small scales. In this review we tried also to give a complete description of the theory behind these models, describing the striking phenomena of BEC and superfluidity. We stressed the difference between models, showing their different descriptions and phenomenological consequences. Since the understanding of this DM class of models require knowledge from several fields, including condensed matter physics, we believe that this theoretical summary of these constructions, is very important to better understand of all the features that these models of DM can present and even help in finding new observables for these models. It is also important for future progress, since there is still a lot of room for theoretical development of these DM models.
Maybe one of the biggest challenge for these models and testing them against observations is the need for numerical simulations. Those are necessary so we can understand how the formation of structures proceeds in these models. However, performing those simulations in a way that they resolve the small scales in order for us to see the interesting effects coming from the wave nature of these models, and also that they simulate the structure on cosmological scales is a challenge. Incredible advances in this field have been made in the past few years and there are many groups currently working to improve those simulations. Most of those advance, tough, are only for the FDM model. Simulations for the SIBEC are only a few and there is still no numerical framework to study the DM superfluid. Thus we should expect continue progress in this field in the future, which will lead to also a better understanding on how to probe those models with observations.
In summary, the study of the ultra-light DM is an active area of research and many challenges are still opened to be addressed theoretically, numerically and observationally. As this field becomes more and more popular we believe this progress will go even faster and we hope this review can help those entering or already in this field supporting their understanding of this fascinating new DM model.