Unified Superfluid Dark Sector

We present a novel theory of a unified dark sector, where late-time cosmic acceleration emerges from the dark matter superfluid framework. The system is described by a superfluid mixture consisting of two distinguishable states with a small energy gap, such as the ground state and an excited state of dark matter. Given their contact in the superfluid, interaction between those states can happen, converting one state into the other. This long range interaction within the superfluid couples the two superfluid phonon species through a Josephson/Rabi cosine potential. As a consequence of this potential, a new dynamics of late-time accelerated expansion emerges in this system, without the need of dark energy, coming from a universe containing only this two-state DM superfluid. The superfluid nature of the interacting species offers a number of advantages. First, because the phonon kinetic term is non-canonical, the decay constant appearing in the potential can be sub-Planckian, unlike the super-Planckian decay constants required in axion models of dark energy. Second, because the superfluid species are non-relativistic, their sound speeds remain suitably small throughout the evolution. We calculate the expansion history and growth of linear perturbations, and compare the results to $\Lambda$CDM cosmology. For the fiducial parameters studied here, the predicted expansion and growth function are close to those of $\Lambda$CDM, but the difference in the predicted growth rate is significant at late times. The present theory nicely complements the recent proposal of dark matter superfluidity to explain the empirical success of MOdified Newtonian Dynamics (MOND) on galactic scales, thus offering a unified framework for dark matter, dark energy, and MOND phenomenology.

Arguably, the most pressing challenge is the mass discrepancy acceleration relation (MDAR). The MDAR is a remarkably tight correlation between the dynamical gravitational acceleration inferred from rotation curves and the gravitational acceleration due to baryons only, as inferred from the distribution of stars and gas [21,22]. In plain DM parlance, the MDAR implies that by looking only at the baryon mass distribution, one can infer the DM density profile at every radius in the galaxy, even in galaxies where baryons are everywhere subdominant. At large distances in the disk, the MDAR implies the baryonic Tully-Fisher relation (BTFR) [23,24], which relates the total baryonic mass to the asymptotic/flat rotation velocity with remarkably small scatter [25]. In the central region, the MDAR implies the correlation between stellar and dynamical surface densities in disk galaxies [26].
The MDAR was of course predicted by Milgrom over thirty years ago [27]. Milgrom's law states that the total gravitational acceleration a is approximately a N in the regime a N a 0 , and approaches the geometric mean √ a N a 0 whenever a N a 0 , where a N is the usual Newtonian gravitational field generated from the observed distribution of baryonic matter alone, and the transition acceleration scale is a 0 ∼ 10 −8 cm/s 2 . Historically this empirical law was proposed as an alternative to DM known as MOdified Newtonian Dynamics (MOND). By now this possibility seems unlikely, given aforementioned large body of evidence for DM as a collisionless fluid on cosmological and cluster scales.
The model of DM superfluidity [28][29][30][31] proposes to successfully marry the phenomenological success of the ΛCDM model on cosmological scales with that of MOND on galactic scales. (See also [32][33][34][35][36][37][38].) Through the well-known physics of superfluidity, the collisionless fluid and "modified gravity" phenomena are two sides of the same coin, representing different phases of the underlying DM theory. The DM candidate in this case consists of axion-like particles with sufficiently strong self-interactions such that they thermalize in galaxies. With m ∼ eV, their de Broglie wavelengths overlap in the (cold and dense enough) central region of galaxies, resulting in Bose-Einstein condensation into a superfluid phase.
The superfluid nature of DM dramatically changes its macroscopic behavior in galaxies. Instead of behaving as individual collisionless particles, the relevant low energy excitations are phonons. It is crucial that the DM phonons must couple to baryons in a way that explicitly (but softly) breaks the global U (1) symmetry of the superfluid. As a result, the DM phonons mediate a long-range force that effectively modifies gravity. (See [39] for a recent proposal for obtaining the MDAR from particle interactions between DM and baryons.) For a particular choice of the superfluid equation of state, the resulting phonon effective Lagrangian is similar to the MOND scalar field theory [40]. Remarkably, this phonon effective theory is strikingly similar to that of the Unitary Fermi Gas (e.g., [41]), which has generated much excitement in the cold atom community in recent years. In recent work [42], the finite-temperature equation of state of DM superfluids was computed using a self-consistent mean-field approximation.

A. A unified description of superfluid DM and DE
A natural question is whether late-time cosmic acceleration can also emerge in the DM superfluid framework, as yet another manifestation of the same underlying substance as DM and MOND. In this paper we will argue that this is indeed possible. We propose a model where the DM consists of a mixture of two superfluid species, arising from two distinguishable states of DM separated by a small energy gap. Although we will not commit to a particular microphysical realization, a natural possibility is that DM is made of dark atoms [43,44] with hyperfine splitting of the ground state. In this case, the two-state superfluid framework can be thought of as a refinement of the previous model that now takes into account the atomic structure of DM.
Since DM particles in these different states are in contact with each other in the superfluid, it is natural to assume that they can scatter via contact interactions. Such interactions lead to a conversion between the different states (or species) of DM particles of the superfluid. This interaction is known as a Rabi coupling or an internal Josephson interaction e.g., [45,46]. The presence of the two states and their interaction has consequences for the collective behavior of the superfluid, altering some of its properties, such as changing the ground state of this system. The study of such systems and its properties is a very active field in condensed matter physics, being studied theoretically [47] (see [48][49][50] as examples of recent developments) and experimentally (e.g., [51,52]). We are interested in understanding the implications for cosmology of DM being described by a two-component interacting superfluid.
A superfluid is described in terms of its collective variable, namely the condensate wavefunction describing the entire superfluid. At low energy, the relevant degree of freedom is the phase of this wavefunction, which is the Goldstone mode for the spontaneous U (1) breaking. Its excitations are phonons (sound waves). In a two-state superfluid mixture, each superfluid specie has different phonon excitations, θ 1 and θ 2 , describing two different sound waves propagating through the system. The contact interaction between particles in the two different states results in an oscillatory potential that depends on the difference of the phonon phases: where ∆E is the small energy gap between the two distinguishable states of DM. (We work in natural units throughout the paper, unless stated otherwise.) This potential explicitly breaks the individual global U (1) symmetries down to a diagonal U (1) subgroup. Physically, this encodes the fact that particles of one species can be converted into the other, but the total number of particles is conserved. Such interaction potentials are ubiquitous in systems with multiple superfluid/superconducting species. Through this potential energy, a multi-state interacting DM component can give rise to new emergent dynamics coming from the collective behavior of this superfluid -a dynamical phase of cosmic acceleration. This phase will occur at the present time provided that M 4 is of order ∼ meV 4 . The dimensionless symmetry breaking parameter f , related to the physical decay constant, must be large enough so that the accelerating phase is approximately de Sitter. Despite the similarity with axions and pseudo-Nambu Goldstone Boson (pNGB) models [53,54], there are important differences. Firstly, the kinetic terms for the angular fields θ i are non-canonical. They are described by non-linear P (X) functions, which encode the superfluid equation of state. Because of this difference in kinetic term, we will see that the DM phase is obtained without oscillations of the field around the bottom of the potential, thus avoiding potential instabilities [55]. Furthermore, unlike the super-Planckian decay constants required in pNGB models of DE, our model can accommodate sub-Planckian values. Thirdly, the microphysical origin of the cosine potential is of course different in the superfluid context. The potential depends on the difference of the phases of the superfluid, which are low-energy collective excitations superfluid, and it arises from a long range interaction within the entire superfluid. This is different than a self-interaction potential for a fundamental scalar field, as in the case of an axion.
This type of superfluid system is realized in Nature in known effects such the Josephson tunneling [45,46]. It is observed in laboratory experiments in models like Iron-based superconductors [56], MgB2 [57], high-T c cuprate superconductors or the XY model [58]. This gives us a good motivation for searching for an analogue of these models in the universe. The other advantage is that this model is minimal, only requiring the presence of this DM superfluid, providing a dynamical explanation for the cosmic acceleration without the presence of DE. This modification of the dynamics is obtained without modifying the underlying fundamental gravity theory of our universe. So, we can say that this model unifies the DM and DE paradigms (and MOND), since we have an universe that has only DM, and where DE or the accelerated expansion is an emergent characteristic when considering a two-state superfluid.
The idea of unified DM and DE is not new, and is commonly known in the literature as Unified Dark Matter (UDM), Unified Dark Energy (UDE), or quartessence [59]. It comprises many models like the generalized Chaplygin gas [60][61][62] and k-essence [63,64] (see [65] for a extensive list of models). These models usually have a stress-tensor that describes the observational eras, mimicking the background evolution from the matter-dominated and DE-dominated eras. Another recent example is [66], where cosmic acceleration arises from suitable interactions between DM and baryons, without sources of negative pressure or new degrees of freedom beyond DM and ordinary matter.
However it has proven challenging to build a successful model that gives the desired background evolution together with a realistic and well-defined growth of perturbations, for the following reasons: • In Chaplygin gas models, which postulates an exotic equation of state P = −ρ −α for the dark fluid, the adiabatic sound speed becomes relativistic at late times, resulting either in large oscillations or late-time growth in the matter power spectrum [62,65]. This can be avoided if α assumes unnaturally small values, less than 10 −5 , in which case the model is indistinguishable from ΛCDM [62].
• In ghost condensation [67] (and the closely related k-essence examples [63,64]), the background energy density at the ghost condensate point mimics a cosmological constant and acts as DE, while small deviations in the scalar field away from the ghost condensate point behave as a pressureless fluid. Thanks to the spontaneous breaking of Lorentz invariance, linear perturbations have nearly vanishing sound speed, resulting in well-behaved perturbations in the linear regime. When these perturbations grow non-linear, however, the pressureless fluid develops caustics, whose resolution requires a UV completion for the ghost condensate [68]. Any attempt to describe DM as a fundamental scalar field, such as DBI models, faces the same issue [69][70][71][72].
Our model naturally avoids both issues. Superfluid phonons are described by a non-linear kinetic term of the P (X) form. The superfluid is assumed non-relativistic, hence its sound speed is suitably small throughout the evolution, which is given completely by the dynamics of the system, from the matter-dominated to the DE-dominated phase. At this level the situation is similar to ghost condensation (or, more precisely, tilted ghost condensation [73,74]). A key difference, as mentioned earlier, is that the phonon field is not a fundamental scalar field. When perturbations become non-linear, the fluid caustics are naturally resolved -the large gradients result in a local breakdown of superfluidity, the P (X) effective description is no longer valid, and the physics is instead described by individual DM particles in the normal phase. This is the mundane "UV completion" of the P (X) phonon effective theory. (Once the halo virializes and achieves a quasi-static state, the superfluid phase is restored in the central region where DM thermalization and condensation occur.)

B. Outline
This work is organized as follows. Section II contains a review of the superfluid DM model. We begin by illustrating in Sec. III our mechanism in the simplest example of two weakly-coupled BEC DM superfluids. In Sec. III C, in particular, we use the results of the cosmological analysis to place constraints on the theory parameters, to ensure that the fluids' sound speeds are sufficiently small, and that the DE component approximates a cosmological constant, consistent with observations. This will place a lower bound on the decay constant appearing in the potential. In contrast to canonical scalar field models of DE, where the decay constant must be super-Planckian, we will see that the aforementioned bound can be satisfied with sub-Planckian decay constants in our model. We then generalize the analysis in Sec. IV and consider two coupled superfluids with arbitrary equations of state.
In Sec. V we turn our attention to the background cosmological evolution. Remarkably we will be able to derive a closed equation for the Hubble parameter that holds for any superfluid equation of state. We numerically integrate this background equation in Sec. V C for a set of parameter values that respect the various constraints derived in earlier Sections. In Sec. VI we study the growth of linear perturbations in the model and derive an equation for the evolution of the total matter perturbation. This equation is numerically integrated in Sec. VII, where we compare the growth function and growth rate for our model to the ΛCDM predictions. In a nutshell the growth function is close to ΛCDM, but the difference in growth rate is substantial at late times. We briefly summarize our results and mention future lines of research in Sec. VIII.

II. REVIEW OF SUPERFLUID DM
In order for DM to be in a superfluid phase in the central regions of galaxies, the particles must be light enough that their de Broglie wavelengths overlap. This is a necessary condition for Bose-Einstein condensation, assuming weak coupling. In other words, the de Broglie wavelength λ dB ∼ 1 mv must be larger than the mean interparticle separation ∼ (m/ρ) 1/3 . This translates to an upper bound on the mass: See [31] for a more refined version of this bound in the superfluid DM context using explicit density and velocity profiles. The general lesson is that superfluidity requires m to be less than order eV. As a second requirement to form a superfluid, the DM particles should interact sufficiently strongly (σ/m 0.1 cm 2 /g) to achieve thermalization. Thus the particles are axion-like in the sense of being light and produced out of equilibrium, but because of the need for strong self-interactions they cannot be QCD axions.
As shown in [31], the resulting density profile consists of a superfluid core, within which DM particles interact sufficiently frequently to thermalize and form a Bose Einstein condensate, surrounded by an envelope describing approximately collisionless particles following an NFW profile. The size of the superfluid region depends both on the model parameters as well as the total mass of the halo.
In order for phonons to explain rotation curve properties, the superfluid core should at least encompass the entire range over which rotation curves are measured, e.g., out to ∼ 60 kpc for a Milky Way-like galaxy. Meanwhile, it is well known that galaxy clusters pose a problem for MOND. Therefore, in galaxy clusters, the superfluid region should be small enough that most of the DM is in the normal phase. In [31] it was shown that these requirements can be simultaneously satisfied for a range of parameters. For instance, with mass m = 1 eV, cross-section σ m = 0.5 cm 2 g (satisfying the bound from merging clusters [75][76][77][78]), and other fiducial parameter values reviewed below, the superfluid cores make up a modest fraction of the total DM mass in galaxies, ranging between ∼ 20 % for a M ∼ 10 12 M highsurface brightness galaxy to ∼ 50 % in a M ∼ 10 10 M low-surface brightness galaxy. Therefore most of the mass is in the approximately collisionless NFW envelope, which has two advantages observationally: 1) Halos are triaxial near the virial radius, exactly as in ΛCDM, and consistent with observations [79]; 2) Most of the gravitational lensing signal comes from the NFW envelope, thus no phonon-photon non-minimal coupling is required to reproduce galaxy-galaxy lensing statistics. Consequently both photons and gravitons propagate at the speed of light, consistent with GW170817 [80].
Within the superfluid region, the physics is described by a phase of a spontaneously broken global U (1) symmetry. At low energy the relevant degrees of freedom are phonons, which are excitations of the Goldstone boson θ for the broken symmetry. The effective theory of phonons must be invariant under the shift symmetry, θ → θ + c, which nonlinearly realizes the U (1) symmetry, and Galilean symmetry, appropriate for a non-relativistic superfluid. Therefore, its most general form at leading order in derivatives and zero temperature must be a "P (X)" theory [81,82] L phonons = P (X) ; where Φ is the gravitational potential. The equation of state of the superfluid is encoded in the form of P (X). Indeed, at finite chemical potential θ = µt, and with the gravitational potential set to zero, the Lagrangian reduces to P (µ), which defines the grand canonical equation of state. In turn, the number density of particles in the condensate is n = P ,X (µ), with mass density ρ = mn. These relations can be combined to deduce the pressure as a function of density, P (ρ). The phonon sound speed is given by c 2 s = P ,X ρ ,X = 1 m P ,X P ,XX . Superfluids are often described by a polytropic equation of state, P (X) ∼ X n , corresponding to P (ρ) ∼ ρ n n−1 . For instance, a standard, weakly-coupled superfluid corresponds to n = 2: This is the case studied in the context of BEC DM [83][84][85][86][87][88][89][90][91][92][93][94]. As another example, the Unitary Fermi Gas, describing a gas of ultra-cold fermionic atoms tuned at unitary, corresponds to n = 5/2 [95]: where c 0 is a dimensionless constant.
The DM superfluid considered in [28,29] corresponds to n = 3/2. More precisely, to ensure that the action is well-defined for X > 0 (time-like profile) and X < 0 (space-like profile), the actual superfluid action is The square-root form also ensures that the Hamiltonian is bounded from below [96]. Modulo the square root and the fact that θ is a non-relativistic field, (6) closely resembles the Bekenstein-Milgrom action to describe MONDian dynamics [40]. Importantly, unlike Bekenstein-Milgrom the additional long-range force is not fundamental. Instead it is an emergent phenomenon due to collective excitations of DM.
To mediate a force between baryons, the DM phonons must couple to the baryon density as where α is a constant. (The reduced Planck scale M Pl is related to the gravitational constant by M 2 Pl = 1 8πGN .) This coupling explicitly breaks the U (1) shift symmetry, albeit softly given the M Pl suppression, making θ a pseudo-Goldstone boson. Physically, this operator implies that DM particle number conservation is violated in the presence of baryons. Its condensed matter origin is unclear and begs for a compelling interpretation. (Alternatively, the MOND-like behavior can come from next-to-leading order terms in the derivative expansion for phonons [30].) In the regime where phonon gradients dominate, such that X − ( ∇θ) 2 2m , the phonon-mediated acceleration matches the deep-MOND expression where a b is the Newtonian gravitational acceleration due to baryons only. The critical acceleration a 0 is related to the theory parameters as The total force experienced by baryons is the sum of the phonon-mediated force and the Newtonian gravitational acceleration due to baryons and the DM condensate itself. Using a combination of analytical and numerical calculations, [31] showed that explicit rotation curves of both high and low surface brightness galaxies could be reproduced in the superfluid model, with excellent agreement with observations. Interestingly, in contrast with MOND where the rotation curve becomes flat, [31] found a slight rise in the asymptotic velocity of massive galaxies. This is due the gravitational pull of the superfluid DM mass itself, which can be as large as ∼ 30% of the total force. To compensate for this effect, [31] found a best-fit value for the critical acceleration of which is somewhat lower than the MOND best-fit value of 1.2 × 10 −8 cm/s 2 . For the record, the corresponding parameter values were Λ = 0.05 meV and α = 5.7.
We have glossed over a subtlety of the superfluid framework, namely that perturbations around the static "MON-Dian" background are unstable (ghost-like). As argued in [28,29], this instability can naturally be cured by finitetemperature corrections, which are expected given the non-zero temperature (i.e., velocity dispersion) of DM particles in galactic halos. The temperature dependence in the effective theory is required independently to obtain an acceptable background cosmology and linear perturbation growth. To simplify the present analysis, we will ignore the issue of temperature dependence and focus our attention on the zero-temperature P (X) Lagrangian.

III. A SIMPLE EXAMPLE: WEAKLY-COUPLED NON-DEGENERATE SUPERFLUIDS
In the previous section we showed that a single species superfluid can describe the DM of our universe describing the MOND and particle-like dynamics of DM on small and large scales. In this section, we are going to generalize the model to a model where DM consists of two non-relativistic superfluid species, described in terms of two distinct phonon excitations. For instance, these could represent two distinguishable states of DM with slightly different energies, ∆E m, such as a ground state and an excited state. As shown in the Appendix, our two-component superfluid can be described in the mean-field approximation by two complex scalar fields, To set the stage, it is instructive to consider the simplest example of two weakly-coupled BEC DM non-relativistic superfluids (4), with mass m 1 and m 2 , respectively: where θ i is the phonon excitation for the i th field, coming from Θ i = m i t + θ i . In the Appendix we review how the above theory derives as the non-relativistic limit of a Lorentz-invariant theory of two complex scalar fields with quartic interactions. Thus far the two species are non-interacting, and the above theory enjoys a U (1) × U (1) global symmetry, describing particle number conservation of each species separately. The conserved number densities are The global U (1) symmetries act non-linearly on the Goldstones as shift symmetries θ i → θ i + c i . We introduce now the Rabi coupling. In the mean-field approximation, it is given by an interaction term that breaks the symmetry group to a residual global symmetry U (1) × U (1) → U (1): As shown in the Appendix, in the non-relativistic regime this interaction between the different states results in an oscillatory potential that couples the two species together: where ∆E ≡ m 2 − m 1 is the energy difference between the two states. 1 In other words, one can either think of the potential (14) as being added directly in the non-relativistic theory of the phonon fields θ 1 and θ 2 , or as resulting from the interaction term (13) in the mean-field description. In addition to the parameters from the DM superfluid model, Λ i 's and m i 's, the potential introduces two parameters: i) The scale M , which sets the amplitude of the potential.
To generate late-time cosmic acceleration, this will be set to the standard value M ∼ meV; ii) The dimensionless symmetry breaking parameter f . As mentioned in the Introduction, a cosine interaction is ubiquitous in systems with multiple superfluid/superconducting species, interacting through Josephson or Rabi couplings. Consistent with the non-relativistic approximation, we assume the mass splitting is small compared to the mass, ∆E m i . Moreover, to simplify the analysis we will assume that the mass splitting is large compared toθ i . Since the latter is ∼ Hθ i cosmologically, this amounts to requiring ∆E H init , where H init is the initial value of Hubble for the period of interest. In this case, the potential can be approximated as To summarize, in the non-degenerate case of interest our assumption is that the mass splitting satisfies This parametric window exists as long as H init m i , which is in any case required by the non-relativistic approximation.
The potential explicitly breaks the individual U (1) symmetries down to the diagonal U (1) subgroup that shifts the Goldstones by the same constant: θ i → θ i + c. The charge densities n i are no longer separately conserved, but the total density, is conserved. This represents the total number density of DM particles in the condensate. Our effective description is valid as long as the number density of particles in each condensate is positive, that is, as long as X i ≥ 0.
Adding the potential to the action, we obtain Note that the coupling of V to the gravitational potential follows in the weak-field limit by integrating out the spatial scalar potential Ψ and ignoring M Pl -suppressed non-local terms. This action should be supplemented by the gravitational action L grav = −M 2 Pl ( ∇Φ) 2 . For later purposes it will be convenient to record the sound speed of perturbations. Expanding the kinetic term in (19) around a time-dependent backgroundθ i (t) to quadratic order in field perturbations ϑ i = θ i −θ i (t), we obtain where the sound speed of each component is We can already see that for each of the superfluid DM species, the sound speed is always small given the non-relativistic approximation X i m i , which tells us that these components cluster like the usual dust component. It is important to notice that this is a dynamical statement and it does not depend on the parameters of the model.

A. Diagonalization
We can gain further intuition by performing the field redefinition This diagonalizes the kinetic term, L = 1 2ξ 2 + 1 2χ 2 + . . ., but the spatial gradients are still mixed. For our purposes it will suffice to work at leading order in ∆E mi 1. The relevant terms in this approximation are where we have defined The last gradient term in (23), proportional to ∆E m ( ∇ξ) 2 , appears naively suppressed, but as we will see it makes an unsuppressed contribution to the hydrodynamical equations and hence must be kept.
In this form, we recognize a Goldstone boson ξ associated with total particle number conservation, and a pseudo-Goldstone boson χ with periodic potential Since both fields have canonical kinetic terms, we can immediately identify the physical decay constant in the cosine potential as We will see that, unlike pNGB models of DE, f χ is not forced to assume super-Planckian values in order for the potential to drive approximate de Sitter expansion. Meanwhile, expanding V to quadratic order in χ, the effective mass for the pseudo-Goldstone boson is readily identified:

B. Hydrodynamical equations
Since our theory describes two interacting superfluids, it is instructive to write down their equations of motion in terms of fluid variables. The continuity and Euler's equations are first-order equations, hence to derive them we must work in the Hamiltonian description.
The conjugate phonon momenta follow immediately from (23): The Hamiltonian density H =ξΠ ξ +χΠ χ − L is then Hamilton's equations of motion arė where in the last step we have assumed, as in (15), that Moreover, here and henceforth, we have defined In other words, V (∆E t) ≡ dV (∆E t) d(∆E t) . To cast the above equations as fluid equations, we must identify the hydrodynamical variables. The fluid densities ρ ξ and ρ χ can be read off from the coefficient of Φ: Meanwhile, the fluid velocities can be identified by taking the spatial gradient of the first two equations in (30). By the Equivalence Principle, the result should be interpreted as˙ u ξ = − ∇Φ + . . .,˙ u χ = − ∇Φ + . . . This allows us to identify: In terms of these fluid variables, and the relative velocity u ≡ u ξ − u χ , the equations of motion (30) becomė The first pair of equations are Euler's equations; the second pair are continuity equations. Later on, in Sec. VI, we will use the above hydrodynamical equations to derive the evolution of density perturbations.

C. Phenomenological constraints
We derive various constraints on the weakly-coupled model to ensure consistency with what is known observationally about DM and DE. The model involves four parameters m 1 , m 2 , Λ 1 and Λ 2 , to describe the superfluid species, as well as two additional parameters M and f describing their interaction. As before we assume m ∼ eV in order for DM to form a superfluid inside galaxies, and M ∼ meV in order for the potential to drive cosmic acceleration at the present time.
We now derive a constraint on the Λ i 's, which later on will be used to constrain the decay constant f χ . To simplify the analysis, we assume that the parameters of the two superfluids are all of the same order. This is already true of the masses, m 1 m 2 m, given the non-relativistic assumption ∆E m i . But we further assume where Λ was defined in (24), as well as The assumption (35) is fairly natural if the two superfluids originate from different states of DM. It follows from (21) that the sound speeds are nearly the same, and given by The number density is conserved, as mentioned earlier, which in a cosmological context implies that ρ = mn ∼ 1/a 3 . Hence the sound speeds both redshift as where the subscript"eq" denotes matter-radiation equality, and ρ eq 0.4 eV 4 is the matter density at that time.
An important constraint is that the fluids' sound speeds must be sufficiently small to satisfy observational constraints from the CMB and the large-scale structure. The sound speed is important to determine the scale, the Jeans length, above which perturbations can grow through gravitational instability. Perturbations smaller than the Jeans length do not grow, and instead oscillate. In order to be consistent with the large scale structures we observe in our universe, the sound speed of our components must be small, smaller than 10 −6 at matter-radiation equality, so the Jeans length is sufficiently small that perturbations on scales of cosmological interest can grow unimpeded.
To our knowledge, the constraints placed on the DM sound speed usually assume c s = const., e.g., [97,98], with the result c 2 s < ∼ 10 −6 . As such, this bound does not readily apply to our sound speeds, which rapidly decrease in time as c 2 s i ∼ a −3 . Nevertheless, as a highly conservative bound we will impose, at matter-radiation equality, This is conservative in the sense that, once (39) is satisfied, then c 2 s 10 −6 at subsequent times. Using ρ eq 0.4 eV 4 , this implies a lower bound on the superfluid scale: Interestingly, Λ can be of the same order as m. Since c 2 s ∼ a −3 , by the present time the sound speed has dropped to < ∼ 10 −15 . Thus the sound speeds remain suitably small throughout the epochs of matter domination and cosmic acceleration. This ensures that the growth of density perturbations proceeds uniformly on sub-horizon scales and avoids the undesirable features in the matter power spectrum seen in Chaplygin models [62].

IV. GENERAL SUPERFLUIDS
The generalization of the non-relativistic action (19) for arbitrary superfluids is straightforward: We could of course consider derivative interactions between X 1 and X 2 , but we focus on (41) for concreteness. 2 We will keep the P i (X i )'s general, though our primary interest lies in the cases P (X) ∼ X 2 , corresponding to BEC DM, and P (X) ∼ X |X|, corresponding to the DM superfluid theory of [28,29]. A simplification in the latter case is that we ignore baryons in our analysis, hence the interaction term (7) is irrelevant. The potential V will be kept general as well, though we will be primarily interested in the cosine potential (15). As before, we will assume V (θ 1 − θ 2 + ∆E t) V (∆E t) in the equations of motion. The conserved charge is the total number density of DM particles in the superfluid state: n P 1 ,X1 + P 2 ,X2 .
As before, our effective description is valid as long as the number density of particles in each condensate is positive, In the non-relativistic approximation, the energy density of the superfluids is given by their rest mass energy plus the potential energy Meanwhile, the pressure is given as usual by the Lagrangian density The adiabatic sound speed of each species, governing the growth of perturbations, can be obtained once again by expanding (41) to quadratic order in field perturbations ϑ i = θ i −θ i (t). The result is with the sound speed given by This agrees with (21) in the particular case P i (X i ) = To ensure stability, c 2 s i should be positive definite. In light of (43), this requires As we saw in Sec. III, the sound speed is always positive and small in the non-relativistic regime for the case where the kinetic term was quadratic. This is also automatically satisfied in the case where the non-standard kinetic term can be described by a power-law, P (X) ∝ (X/m) n , as seen in Sec. II. The sound speed of each component in this case is given by c 2 s i ∝ (X i /m i ) n which is positive and much smaller than unity in the non-relativistic regime X i m i . This smallness of the sound speed is a dynamical statement in our theory, not depending on the parameters of the model, and is important so the perturbations of this model produce the large scales structures we observe in our universe.
As before, it is helpful to gain intuition on the parameter f appearing in the cosine by translating it to a physical mass scale analogous to an axion decay constant. Similarly to the field redefinition (22), the pseudo-Goldstone χ is now given by the difference of perturbations: The physical decay constant in the cosine potential is readily identified as The mass of χ is once again given by (27).

V. BACKGROUND EVOLUTION FOR GENERAL SUPERFLUIDS
We now turn to the study of the background cosmology. Remarkably we can derive a universal equation for the Hubble parameter, valid for any superfluid equation of state. To see this, note that, using the total energy density (44), the Friedmann equation for a spatially-flat universe is where we have defined with m = 1 2 (m 1 + m 2 ). On the other hand, using (44) and (45), the second Friedmann equation (the "Ḣ" equation) isḢ where we have used the non-relativistic approximation m i P i ,Xi P i . Combining these two equations, we obtain As claimed, this is a closed equation for H(t), which holds for any choice of P i (X i ). This equation is similar to what one finds in ghost condensation [67]. Equation (54) is all we need to solve for the background cosmology. For completeness it is also instructive to consider the phonon equations of motion: where we have used the fact that V only depends on the difference θ 2 − θ 1 . As before, V denotes differentiation with respect to its argument -see (31). The sum of these two equations implies d dt This is nothing but the conservation of the total number density of DM particles (42). Meanwhile, the difference of the two equations (55) impliesρ This is similar to the equation for a canonical relativistic scalar, except that the potential is explicitly time-dependent. 3 We now show that the background evolution of this model can describe the evolution of our universe from a matterdominated phase to an accelerated phase at late times. We will solve analytically the background evolution only for each of the phases separately. We will come back in Sec. III C and find exact numerical solutions to the background equations in the simplest of weakly-coupled superfluids.
In the following we show the solutions for the dust dominated epoch and the late-time accelerating phase.

A. Dust-dominated phase
Let us consider the dust-dominated epoch, well before the onset of cosmic acceleration, when the potential energy is negligible (i.e., when H 2 M 4 /M 2 Pl ): The solution is H = 2 3t , describing a matter-dominated universe. In this regime (55) reduce to d dt a 3 P i ,Xi 0, telling us that each species is separately conserved: Not surprisingly, the total phonon energy density (44) describes a pressureless fluid: The above pressureless behavior was derived to leading order in Xi mi 1. Of course, more precisely our fluids do have a small pressure, and correspondingly a small adiabatic sound speed given by (47). We must make sure that the c 2 s i 's are sufficiently small to satisfy observational constraints from the CMB and the large-scale structure. In Sec. III C we will derive detailed constraints from observations on the theory parameters for a fiducial example, P (X) ∼ X 2 .
Incidentally, the scalar equations of motion (55), ignoring the potential, boil down tȯ With c 2 s i > 0, we see that X i decreases in time during the dust-dominated phase, hence the non-relativistic approximation X i m i becomes increasingly accurate.

B. Late-time accelerating phase
The dust-dominated phase ends once the potential energy becomes a significant contribution to the energy density at the present time. This occurs when H ∼ M 2 /M Pl . In order for the universe to experience an approximate de Sitter phase, the potential energy should be approximately constant on a Hubble time, that is, For the cosine potential (15) we have Thus, for typical values of ∆E t f the condition (62) will be satisfied if where in the last step we have used (16). (Recall that H init is the value of Hubble at the onset of the evolution, such as matter-radiation equality.) Using (50) this translates to a bound on the decay constant: In particular, for the simplest example of two weakly-coupled superfluids discussed in Sec. III, with f χ given by (26), this gives In this case we can derive further constraints using the bound (40) on Λ. Assuming Λ 1 Λ 2 Λ, as we did in (35), as well as m ∼ eV, the above inequality becomes In other words, since H 0 ∼ meV 2 /M Pl , we obtain If we demand that f χ < M Pl , then this places a bound on H init : This corresponds to an initial (baryon/radiation) temperature of < ∼ GeV. To summarize, by choosing f χ < M Pl , which is desirable, our approximations are self-consistent from an initial temperature of ∼ GeV all the way to the present time. The resulting sound speed is consistent with observational bounds, and the phase of late-time cosmic acceleration is nearly de Sitter. Thus, unlike the super-Planckian decay constants required in pNGB models of DE, our model can accommodate sub-Planckian values. The difference traces back to the different kinetic structure in our case.

C. Numerical solution
To verify the above analytical arguments, we numerically solve the background cosmology of our model, for general superfluids. This amounts to integrating the evolution equation (54) for the Hubble parameter. Since only the potential V (∆E t) appears in that equation, we only need to specify the potential parameters M , ∆E and f , subject to the constraints derived above: • We set the scale of the potential to with H 0 ∼ 10 −33 eV in order to derive cosmic acceleration today.
• For concreteness we integrate the equations from matter-radiation equality onwards, hence for our solution.
• We set the energy gap, which is a characteristic of the microphysics of our superfluid, at Note that this lies in the desired window (16) for all times of interest.
• We choose which therefore satisfies (64). Figure 1 shows the solution (54) for the Hubble parameter for our model (blue solid curve) with the above parameters, as well as a ΛCDM model (black dashed curve) with the same matter density at early times and asymptotic Hubble constant H 0 . We can see that our model describes a phase of dust domination followed by a phase of accelerated expansion. The evolution is indistinguishable from ΛCDM until the present time (t 0 ∼ H −1 0 ). Because our interaction potential is dynamical, the two expansion histories deviate from each other in the future, as the effects of the cosine potential come into play. The characteristic time at which the two histories deviate is set by the period of the cosine, i.e., ∼ 2f H0 ∆E = 10 H −1 0 . Had we continued the evolution further in time, one would see the Hubble parameter oscillate due to the cosine. However it is worth remembering that our superfluid effective theory is only valid in the regime whenḢ ≤ 0, i.e., whenever the Null Energy Condition is satisfied.

VI. GROWTH OF DENSITY INHOMOGENEITIES
We have thus far showed that our model can describe the expansion history of our universe, with a period of decelerated evolution followed by a period of accelerated expansion, as in the ΛCDM model. A viable alternative to the ΛCDM model must not only reproduce the evolution of the background, but it should be able to describe the growth of density perturbations that leads to the structures we observe in our universe. In this section we turn to the analysis of density perturbations.
For simplicity we will focus on the interacting BEC DM superfluids of Sec. III. Our starting point is the set of hydrodynamical equations given by (34). Because these were derived in the weak-field approximation, they can be applied to the cosmological context in the freely-falling coordinate system of the Friedmann-Robertson-Walker (FRW) metric: The proper distance is related to the comoving distance via = a(t) x. The FRW background corresponds to This coordinate system is valid on sub-Hubble scales, H 1. To make the assumed coordinate system explicit, let us rewrite the fluid equations (34) as where (∂/∂t) denotes time-differentiation at fixed . The gravitational potential is determined as usual by Poisson's equation: We now transform to the expanding coordinate system (see, e.g., [99]), with = a(t) x. The fluid densities are The fluid velocities can be decomposed as The first terms account for the Hubble flow, whereas v ξ , v χ are the peculiar velocities. The Hubble flow of course drops out of the relative velocity: Finally, the gravitational potential Φ can be split into a background piece and an inhomogeneous term: Making these substitutions, we obtain the hydrodynamical equations in the expanding coordinate system: A. Non-linear evolution of inhomogeneities Each fluid density can be decomposed into a background piece and an inhomogeneous term: Note that δρ ξ and δρ χ are not assumed small at this stage. The background densities obey the equationṡ This confirms, in particular, thatρ ξ describes dust and redshifts as 1/a 3 . Meanwhile, the evolution of ρ χ is influenced by the potential. It will be convenient to define perturbations relative to the total background density,ρ =ρ ξ +ρ χ , whereρ satisfiesρ Note that the total density perturbations is given by Substituting the decomposition (83) making use of the background equations (84), the fluid equations becomė together with Poisson's equation, The above set of equations are non-linear, Newtonian hydrodynamical equations in an expanding universe. They can be solved numerically to describe the formation of non-linear structures in the Newtonian regime.

B. Linear perturbations
In the linear regime, where the δ's and v's are all small, (88) simplify tȯ In particular, by adding the first and third equations we obtain an equation for δ = δ ξ + δ χ , The above equations, together with Poisson's equation (89), can be combined as usual to obtain a second-order equation for the density perturbation: where we have introduced the sound speed parameters: In deriving the first of (92), we have made use of ∆E V V ∼ ∆E f H, which follows from (67). By differentiating the δ equation once more, we could eliminate v using the second equation. In what follows we will instead argue that the spatial gradient terms can be ignored to obtain a single ordinary differential equation for the density perturbation.

C. Ignoring spatial gradients
The spatial gradient terms in (92) are all proportional to the sound speed parameters c 2 s, ξ and c 2 s, χ . We now argue that these parameters are small, hence the spatial gradient terms can be neglected on linear scales. For c 2 s, ξ , this immediately follows from the analysis of Sec. III C, in particular (39) and (40), since For c 2 s, χ , we shall assume for simplicity that Λ 1 Λ 2 Λ, as we did in (35). Moreover, from (84) we estimate that Unlike c 2 s, ξ , which redshifts as 1/a 3 , we see that c 2 s, χ grows as 1/H and therefore assumes its largest value at late times. Using the inequalities f Hinit H0 (see (64)) and ∆E H init (see (16)), the above implies It is easy to choose parameters that satisfy all of our constraints and yield c 2 s, χ 1. Ignoring spatial gradient terms, the first of (92) reduces Furthermore, the second of (92) reduces to˙ v + H v 0, which implies v ∼ 1/a. Thus the peculiar velocity difference redshifts away.

VII. NUMERICAL EVOLUTION OF DENSITY PERTURBATIONS
The second-order equation (97) for governing the linear growth of density perturbations is straightforward to integrate numerically. For this purpose we assume the same parameter values as used in Sec. V C. We also set m = 1 eV, and This satisfies the bound (40) and ensures that the sound speeds are small, thereby justifying our neglecting the spatial gradients. Note that because we can drop the spatial gradients, the resulting equation of motion (97) is independent of Λ. Furthermore, with f = 10 ∆E 2H0 10 23 (see (73)), the decay constant (26) evaluates to This is M Pl , as claimed earlier.
A. Background solution for the matter densities The first step consists of solving (84) for the background densitiesρ ξ andρ χ . This requires specifying initial conditions, i.e., the value of the densities at matter-radiation equality where we begin our integration. Since our superfluid is a mixture of particles in their ground and excited states, the relative fraction of the two populations depends on how the energy gap ∆E compares to the DM temperature at matter-radiation equality. This in turn depends on the production mechanism of our DM particles.
Assuming that the production mechanism is an axion-like vacuum displacement, then by definition this will take place when H ∼ m = eV, corresponding to redshift z production ∼ 10 16 . Although the DM particles are initially relativistic, T prod ∼ eV, soon thereafter they become non-relativistic, and their temperature subsequently redshifts as 1/a 2 . Hence by matter-radiation equality (z eq ∼ 10 3 ), the DM temperature has dropped to T eq ∼ 10 −26 eV .
Thus, with our chosen value of ∆E = 5 × 10 −11 eV, we have T eq ∆E. Hence, to a good approximation, all particles are in the ground state at that time. That is, the matter density is entirely in θ 1 .
Equations (101) and (102) give the initial conditions for the evolution of the energy densities. Note from (101) that ρ χ is initially negative, however this is an artifact of the definition of the variables ξ and χ in (22). It does not imply any violation of energy conditions, since, as seen in Secs. III and IV, the density of the phonon variables θ 1 and θ 2 is always positive.  Figure 2 shows the evolution, from matter-radiation equality onwards, of the fractional density parameters Ω x = ρx 3M 2 Pl H 2 for our different components: Ω ξ (red), Ω χ (blue), the total matter density Ω ξ +Ω χ (yellow), and the potential energy Ω V driving late-time accelerated expansion (dotted gray). Since we assume a spatially-flat universe, the Ω's add up to unity at all times, Ω total = 1 (dashed black). 4 We can see that, initially, during the matter-dominated period, ξ dominates the density budget of the universe, being responsible for almost all the matter density, while the contribution of χ remains small. As the matter density redshifts, the potential becomes increasingly important. Around the present time (H 0 t ∼ 1), the interaction potential starts to dominate, resulting in a period of accelerated expansion. Note that the model does not solve the coincidence problem, because the onset of cosmic acceleration is set by the choice of parameters.
The evolution can also be understood in the original (θ 1 , θ 2 ) variables. Cosmologically, the phonon fields are rotating along the symmetry-breaking direction with cosine modulations. In the early universe, the fields rotate rapidly, being oblivious to the small cosine modulations, and their kinetic energy dominates. Because of their non-canonical kinetic terms, this results in a DM-dominated universe. As the fields slow down due to Hubble expansion, they eventually feel the influence of the oscillatory potential, which triggers late-time acceleration. Meanwhile, on small, non-linear scales, the phonon gradient terms dominate again, resulting in DM behavior. This offers a compelling, unified picture for the dark sector of cosmology.
Coming back to Fig. 2, for future times (H 0 t > ∼ 1) we see deviations from the ΛCDM evolution, since our potential is dynamical and oscillatory. The density of ξ decays from equality until the future due to the expected redshifting of a dust component, but also because of conversion from the ground state θ 1 to the excited state θ 2 . Given this conversion, the density of χ starts to grow until it dominates the matter energy density, meaning that almost all the particles of the superfluid are in the excited state. As expected, the potential oscillates, becoming subdominant with respect to the DM component, around 7.5 H 0 t 9.5, indicating a new phase of dust-like decelerated evolution, and then dominating the energy density again (for H 0 t > ∼ 9.5) giving a new phase of accelerated expansion. Thus the future behavior of our model is quite distinct from that predicted by ΛCDM, showing an oscillation between epochs of decelerated, matter-dominated expansion and accelerated, potential-dominated expansion. We do not evolve this model past H 0 t 10, since after that time the null energy condition is violated and the effective field description breaks down. . We see that our predicted growth is somewhat suppressed at low redshift compared to ΛCDM. Right: Fractional difference between our model and ΛCDM with the same matter density at early times and asymptotic Hubble constant H0. The difference is less than 2.5% at z = 0 for our chosen parameters.

B. Growth factor
We are finally in a position to solve (97) for the evolution of linear density perturbations for our superfluid mixture. From the form of these equations, notice that the interaction potential affects the growth in two ways: i) through its impact on the background expansion, i.e., its contribution to H(t); and ii) its explicit appearance in the equation for δ. Deep in the matter-dominated epoch, however, the potential is negligible, therefore δ evolve as in CDM.
Let us say a word about our initial conditions, specified at matter-radiation equality. To match the observed primordial amplitude, we set δ eq = 10 −5 at the initial time. Moreover, since δ evolves as a CDM perturbation in the matter-dominated era, we set its initial time derivative to the growing mode behavior δ ∼ a. In other words, we seṫ δ eq = H eq δ eq . With these initial conditions and the parameters mentioned at the outset, we numerically integrate (97). Figure 3 (left panel) shows the rescaled growth function, defined as 5 for our model (solid blue curve) and ΛCDM model (dashed grey curve). Our model predicts a slightly suppressed growth compared to ΛCDM. As seen in the right panel of Fig. 3, the difference is quite small, less than 2.5% at z = 0, for the parameters considered here. The difference with ΛCDM is more pronounced with the growth rate, defined as This quantity is interesting since various probes of structure growth, such as redshift-space distortions, are sensitive to this quantity. In the left panel of Fig. 4 we plot the evolution of the growth rate of our model (solid blue line), and, for comparison, the expected growth rate for the ΛCDM model (gray dashed line) with the same asymptotic Hubble parameter H 0 . Our model deviates substantially ( > ∼ 30%) from the ΛCDM model, and predicts a steeper suppression of density perturbation growth. We should caution that the quantitative difference depends of course on the parameter values chosen here. In a future paper we plan to perform a more systematic study of model predictions and comparison to data, along the lines of [100,101].

VIII. CONCLUSION
In this paper we presented a novel unified framework for the dark sector, based on the idea of superfluid DM. The model relies on two non-relativistic superfluid species, which physically represent the ground state and excited state of DM. The two superfluid species interact through a Josephson/Rabi cosine potential, and this potential energy is responsible for driving late-time cosmic acceleration. The present theory nicely complements the recent proposal of DM superfluidity to explain the empirical success of MONDian phenomenology on galactic scales, thus offering a unified framework for DM, DE, and MOND. (A finite-temperature superfluid has also been shown recently to degravitate a large cosmological constant [102].) Unlike pNGB models of DE, which also rely on a softly-breaking periodic potential, the non-canonical nature of the phonon kinetic term allows for sub-Planckian decay constants. Furthermore, thanks to the non-relativistic nature of the superfluids, their sound speeds remain small throughout the evolution and results in acceptable growth of density perturbations depending only on the expansion rate of the universe. This is unlike other earlier attempts of a unified dark sector, notably the Chaplygin gas models, where the adiabatic sound speed becomes relativistic at late times, resulting either in large oscillations or late-time growth in the matter power spectrum.
The interaction potential assumed in this work is of the general Josephson/Rabi form, which is ubiquitous in systems with multiple superfluid/superconducting species. Mixtures of BECs and superfluids with such interactions have already been created in the laboratory. This motivates us to study the implications of such potential for cosmology. Our study further emphasizes the potentially rich phenomenology of the superfluid framework, in particular the use of collective modes, in generating novel cosmological dynamics.
Leveraging the knowledge from condensed matter systems, many avenues offer themselves along which the present model could be further developed. An interesting future development for this theory could be to describe the Rabi interaction coefficient as a quantity that depends on temperature (or time) or even in the intra-particle distance, trying to give a dynamical explanation for the parameter M that is liked to the onset of acceleration expansion. Since the theory was thus far described in the mean-field approximation at zero temperature, and the conversion is purely quantum, it would be interesting to study corrections to the mean-field approximation including finite-temperature effects. A first step in this direction was taken recently in [42] in this case of a single component DM superfluid.
Our preliminary study of the predicted expansion history and linear growth of perturbations was done assuming a fiducial choice of parameter values, for illustrative purposes. The small, yet significant deviations from ΛCDM (particularly for the growth rate) are tantalizing. They motivate us to explore more systematically the model predictions over a broader range of parameters. In a future publication we plan to perform such a systematic analysis and compare the results to data.
Although the numerical analysis focused primarily on the background expansion history and linear growth of perturbations, we derived analytically the general hydrodynamical equations governing the non-linear evolution of perturbations -see (88). We plan to solve these equations to study structure formation, for instance within the spherical collapse model. Along these lines, it would be interesting to study local effects coming from higher-temperature regions, such as galaxy clusters, where one might expect different relative populations for the species than the cosmo-logical populations.
Appendix A: Derivation of the non-relativistic superfluid action from a Lorentz invariant theory In this Appendix we show how the non-relativistic theory (19) for two weakly-coupled superfluids can be derived from a Lorentz-invariant action. Our starting point is the relativistic theory of two complex scalar fields with quartic interactions: where the form of the dimensionless couplings λ i = 2m 4 i /Λ 4 i was chosen for later convenience. Decomposing the fields in polar variables, the action can be written as Thus far the two species are non-interacting, and the above theory enjoys a U (1)×U (1) global symmetry, describing particle number conservation of each species separately. The conserved currents are Below we will turn on a weak interacting potential that will break this symmetry down to a diagonal U (1) subgroup. By definition, a superfluid is a phase of spontaneously broken U (1) symmetry, at finite charge density n i = j i 0 ∼ Θ i = 0. Integrating out ρ i to leading order in derivatives, we obtain For consistency, the state of interest must have − (∂Θ i ) 2 ≥ m 2 i . The corresponding charge densities are Substituting (A5) back into (A3) gives the action for the phonon fields Θ 1 and Θ 2 , to leading order in the derivative expansion: The global U (1) symmetries act non-linearly on the Goldstones as shift symmetries Θ i → Θ i + c i . At face value (A7) is identical to the action for two ghost condensates [67]. However there are two important differences. The first difference pertains to the underlying description. Ghost condensation describes a fundamental scalar field, which spontaneously breaks Lorentz invariance. The superfluid effective description instead describes a collective degree of freedom, i.e., sound waves. Secondly, the ghost condensate effective theory allows for field configurations withΘ 2 < m 2 , resulting in violations of the Null Energy Condition. In the superfluid description, however, such a regime corresponds to negative particle number density, as can be seen from (A6), and is therefore forbidden.
We introduce now the Rabi coupling, that in the mean field approximation is given by an interaction term that breaks the symmetry group to a residual global symmetry U (1) × U (1) → U (1): In the non-relativistic regime, this becomes as oscillatory potential for the phonon difference: This potential explicitly breaks the individual U (1) symmetries down to the diagonal U (1) subgroup that shifts the Goldstones by the same constant: Θ i → Θ i + c. The charge densities n i are no longer separately conserved, but the total density, is conserved. The Lagrangian is the sum of (A7) and (A9). At this point we take the non-relativistic limit, writing the phonons as withθ i m i . Similarly the metric takes the weak-field, Newtonian form g 00 −(1 + 2Φ), with Φ being the gravitational potential.
To leading order inθ i , the densities n i for each species, given by (A6), reduce to their non-relativistic expressions (12). Furthermore, the Lagrangian reduces to This agrees with the non-relativistic action (19) in the main text. Clearly, by generalizing the potential in (A1) for Ψ 1 and Ψ 2 , we can obtain different P (X i ) superfluid kinetic terms, and arrive at (41) through similar steps.