Novel Substructure and Superfluid Dark Matter

The recent observation of the distribution of accreted stars (SDSS-Gaia DR2) suggests that a non-trivial fraction of dark matter is contained within halo substructure. With this in mind, in this letter we construct novel solutions to the equations of motion governing condensate dark matter candidates, namely axion Bose-Einstein condensates and superfluids. These solutions are highly compressed along one axis and thus have a disk-like geometry. We discuss linear stability of these solutions, consider the astrophysical implications as a large-scale dark disk or small scale substructure, and find a characteristic signal in strong gravitational lensing. If observed, such substructure is a smoking gun signal of condensate models of dark matter. This indicates that dark matter substructure is a powerful new observable window for testing the nature of dark matter.

The recent observation of the distribution of accreted stars (SDSS-Gaia DR2) suggests that a non-trivial fraction of dark matter is contained within halo substructure. With this in mind, in this letter we construct novel solutions to the equations of motion governing condensate dark matter candidates, namely axion Bose-Einstein condensates and superfluids. These solutions are highly compressed along one axis and thus have a disk-like geometry. We discuss linear stability of these solutions, consider the astrophysical implications as a large-scale dark disk or small scale substructure, and find a characteristic signal in strong gravitational lensing. If observed, such substructure is a smoking gun signal of condensate models of dark matter. This indicates that dark matter substructure is a powerful new observable window for testing the nature of dark matter.

I. INTRODUCTION
Recent data [1] has revealed that a non-trivial fraction of dark matter in the local dark matter halo is in substructure [2,3], that is, in gravitationally bound clumps distinct from the halo. Meanwhile, the nature of dark matter remains elusive, and the strongest evidence for dark matter remains gravitational. It is tempting to speculate that dark matter substructure may be a powerful tool to discriminate between classes of dark matter models. Indeed, an intriguing feature of so-called 'Fuzzy Dark Matter' [4,5] is a resolution of the 'missing satellites problem' [6], though recent evidence indicates that there is no such problem in need of solving [7].
Particle dark matter predicts spherical sub-halos on all scales [8], while [9,10] has argued that a strongly interacting subcomponent can lead to a "dark disk" aligned with the visible disk [11]. Condensate dark matter scenarios, such as axion Bose-Einstein condensates [4,5,12], bosonic superfluids [13][14][15] and fermionic superfluids [16,17], also exhibit substructure, and to-date the known forms are spherical clumps [5,18,19] and vortices [20]. All three condensate systems, in the non-relativistic limit, are described by the same set of equations, exact solutions to which should exist as dark matter substructure. This substructure in principle provides a smoking gun signal of condensate dark matter scenarios.
However, it is a notoriously difficult to task to find bound state solutions without spherical symmetry, see e.g. footnote 21 of [5]. In this letter we study a novel form of substructure that can exist in condensate models. We look for and find disc-like structures, that can exist both as a large-scale dark disk or else as isolated substructures. This investigation reveals new solutions to old equations, and adds qualitatively new results to the existing mathematical literature.

II. THE EQUATIONS OF MOTION OF NON-RELATIVISTIC DARK MATTER
The coupled Gross-Pitaevskii (non-linear Schrödinger) and Poisson equations emerge in the non-relativistic limit of a classical scalar field theory interacting with gravity. This arises in an astrophysical context as cold dark matter at high number density, e.g. an ultralight scalar [4,5]. The non-relativistic limit is defined via the decomposition and the limit |ψ| mc 2 |ψ|/ [5]. The equations of motion are then given by [21], where ∇ 2 is the Laplacian in three spatial dimensions, and we set c, = 1. Here λ > 0 corresponds to an attractive interaction, and we will take λ > 0 throughout. We will consider a static solution with an energy E, such that ψ(x, t) =ψ(x)e −iEt . The equations of motion can be put in a dimensionless form via the redefinitions, where m pl is the reduced Planck mass. The characteristic distance scale R is thus set by the coupling λ and the mass m. The equations of motion are then given in terms of the hatted quantities by, where∇ denotes the Laplacian with respect to the hatted coordinates, ∇ 2 ≡ ∂ 2 ∂x 2 + ∂ 2 ∂ŷ 2 + ∂ 2 ∂ẑ 2 .

III. SOLUTIONS AWAY FROM SPHERICAL SYMMETRY
It is conventional to solve the above equations assuming spherical symmetry of the solution [18,19] [22]; see also [23][24][25][26][27] for the case λ = 0. In our work, we are interested in disk-like solutions, which instead have an axial symmetry, i.e. invariance with respect to rotations about a fixed axis. As previously noted, it is a difficult to task to find bound state solutions to (5) without spherical symmetry. This difficulty is due in large part to the non-linearity of (5) and to the fact that axially symmetric bound state solutions must have (at least) two independent variables. However, symmetry-reduced phenomena are known to exist in closely related systems (e.g. [28]), and thus it is not unreasonable to propose that such stable or metastable solutions exist in three dimensions with gravity.
To this end, we note that a simple class of disk-like geometries in R 3 can be described in terms of a squeezed radial coordinate with real-valued D sq > 2. For large D sq , functions of the form ψ(x, y, z) = ψ(r sq ) which decay at infinity are disk-like in R 3 . Interestingly, the Laplacian acting on a wavefunction ψ(r sq ) takes a simple form when expanded in powers of z/D sq : where we drop the "sq" subscript for notational convenience. When D is a positive integer, one recognizes the leading order term of (7) as the radially symmetric Laplacian in D+1 spatial dimensions. Our present work allows for D to be an arbitrary positive real number so that it may be interpreted as interpolating between dimensions corresponding to integer values of D.
It follows that disk-like solutions may be described in the region z/D 1 as solutions to the radially symmetric equation (7) in very large spatial dimensions. Inserting (7) into (5), and truncating at lowest order in z/D, we have upon dropping the hats and restricting to real functions ψ(r). Equation (8) resembles a core differential equation obtained through the traditional spatial dynamics method of far-field/core decomposition (see e.g. [29]). In the context of physics, this expansion is similar to conventional methods of electrodynamics, e.g. in Fresnel Diffraction [30], and in a more modern context, a similar expansion is appears in the famed Anti-de Sitter/Conformal Field Theory correspondence [31], where an Anti-de Sitter space emerges in the near-horizon limit of stack of coincident three-branes. We expect that solutions ψ(r) of (8) which rapidly decay to 0, i.e. are localized within the region z D, well approximate solutions to the full equation (5). In what follows we will numerically solve (8).

IV. DISK SOLUTIONS
We use the numerical fixed point algorithm described in the appendix to find solutions of (8). Throughout this work we have fixed V (0) = −1, but remark that solutions appear to exist for any choice of V (0) < 0, and therefore we leave an exploration on the dependence of the choice of V (0) to a subsequent investigation.
Our results indicate the existence of solutions of (8) for arbitrarily large D > 0 for which ψ(r) monotonically decreases and approaches 0 as r → ∞. Moreover, we have identified for D a family of solutions which can be parametrized by their value at r = 0, simply written ψ(0), that are well fit by a Gaussian  We note that the emergence of a family of solutions is not entirely surprising since in the absence of the selfinteraction term ψ 3 one may show that solutions exist for any value of ψ(0) > 0 [25]. Hence, when ψ(0) is taken sufficiently small the self-interaction ψ 3 is sub-dominant to the other terms in the Schrödinger equation, and therefor acts as a slight perturbation. Moreover, based on the approximate solution (9), it appears that finite-energy spherically symmetric solutions exist in arbitrarily large number of spatial dimensions. These solutions are characterized by a distance scale in physical, unhatted, coordinates, which defines the 'core radius', or planar extent, of the disk solutions. This can take a broad range of values, for example, for D = 10 5 , m = 10 −10 eV, and λ = 10 −45 , this evaluates to 32 kpc. The disk thickness, i.e. the extent in the z-direction, is independent of D, and is therefore suppressed relative to the core radius (10) by a factor of 1/ √ D. Thus, as anticipated, for large D we find thin disk solutions. Moreover, the super-exponential decay in z is consistent with the assumption of localization of ψ to the region z/D 1, providing an a posteori justification of the use of the Laplacian (7).
The total mass of these solutions is given by, where again x is the physical spatial coordinate, where we have used the analytic (9). Most strikingly, the bound ψ(0) ≤ 1 provides an upper bound on the mass (11) for each fixed D. Therefore, for a fixed M DD > 0, we are able to obtain a lower bound on the squeezing factor D, which in terms of the conserved particle number N ≡ M DD /m, is simply D > 2π 3/2 N √ λ(m/m pl ). The existence of this lower bound hints at the stability of these solutions. Namely, given a perturbation which leaves unchanged the boundary condition ψ(0), the conservation of N forbids D from dynamically relaxing to a small value, and therefore prevents a disk from relaxing to a spherical solution, the latter of which has D = 2 in our convention.
This can be further probed by extremizing the total energy H, with respect to the particle number N [18,19]. At fixed ψ(0) this leads to an expression for the various contributions to H in terms of N and D, which exhibits a stable minimum at particle number N ∝ D, as expected.
Finally, we turn to a numerical investigation of the linear stability of our solutions. In this preliminary investigation we only examine stability with respect to perturbations that are also axially symmetric, i.e. are functions of just r as well. The linear stability is carried out numerically using the methods outlined in the appendix. We find that these solutions have spectrum entirely contained on the imaginary axis, indicating linear stability. A plot of the spectrum for a solution with D = 10 3 is provided in Figure 2 and we remark that for all tested values of D we have obtained linear stability as well. Of course this work only provides linear stability with respect to a limited class of perturbations, but does provide valuable information pertaining to these disk-like structures. We leave a full stability analysis to a follow-up analysis, and note that absolute stability is not required for such objects to be present in the galaxy, which is out of equilibrium [2,3].

V. ASTROPHYSICAL IMPLICATIONS
We have found self-gravitating solutions to the coupled Schrödinger and Poisson equations. These solutions can exist as a dark disk that coincides with the visible disk of the Milky Way, or else as isolated substructure. These two possibilities have distinct observational signatures, which we now consider. To simplify the discussion, we will fix ψ(0) = 1 in what follows.

A. Dark Disk Universe
If a component of particle dark matter has nonnegligible self-interactions, and a mechanism for dissipating energy, it can collapse into a disk aligned with the visible matter [9,10]. The properties of the disk, e.g. the thickness, can be estimated by accounting for the particle physics processes at play (e.g. Compton cooling), but precise estimates can only inferred from yet to be performed N-body simulations. However, if the dark matter particle in question is light, then the correct description following the collapse is the non-relativistic limit utilized here, and hence must be described by one of the solutions found in this work [32]. Interestingly, this possibility is already constrained by data, and in particular, by recent data from the Gaia telescope [1].
The Gaia data provides a loose upper bound on the density of dark matter in a thin dark disk [33,34]. Bounds are typically formulated in terms of the surface density Σ DD , related to the density inside the disk ρ by Σ DD = h DD ρ, where h DD is the disk thickness. The current bound on the surface density is given by Σ DD 5M /pc 2 [33,34], which of course only applies provided the core radius stretches out to kpc scales, R c > kpc in equation (10). This translates to a bound on the self-interaction and the mass, independent of D, m eV For a given mass m, the Gaia constraint thus provides a lower bound on λ. For example, for m = 10 −5 eV and with D fixed to give a disk of radius R c = kpc, the Gaia constraint is a lower-bound on the self-interaction λ > 10 −46 . Interestingly, saturating the bound (12) can correspond to a very small fraction of the total dark matter in the disk. The mass fraction is given by which follows from equation (11) and the NFW profile M N F W ∼ (4π/3)R 3 s ρ 2 0 , which for the Milky Way gives M N F W 3.03 × 10 68 GeV. For example, with λ = 10 −46 and D = 10 4 , the disc is a negligible fraction of the total matter. One could instead consider D = 10 24 , which gives a percent level fraction of dark matter in the disk.
This also has implications for fuzzy dark matter [4,5], namely, the bound (12), which follows from explicit solutions found here, allows for astrophysically interesting dark disks, i.e. with a thickness > 10pc and Σ DD 10M /pc 2 , even for the requisite small mass range m ∼ 10 −22 − 10 −21 eV [4,5]. Demanding Σ DD 10M /pc 2 translates to the bound m/eV 10 12 λ 3/8 . For m = 10 −21 eV this is saturated for λ 6 × 10 −89 , corresponding to a disk of thickness 240pc. Thus astrophysically interesting dark disks can exist in the fuzzy dark matter scenario, provided that λ has an extremely small value.

B. Halo Substructure
We now consider the possibility that these disk solutions could exist as isolated substructures, analogous to traditional particle dark matter sub-halos. These are in principle subject to the constraints traditionally applied to primordial black hole and MACHO dark matter [35][36][37][38], which generally constrain massive compact objects to be a sub-percent fraction of the mass of the halo. For concreteness, here we will consider the effect of a disc that is 1% the mass of the halo. The characteristic size of such a disc is, This can range in size from much less then a pc to vastly more. Thus a great diversity of disk-like subhalos is possible. An interesting probe of such substructures is strong gravitational lensing. Assuming the orientation with respect to the line of sight is random, gravitational lensing by disk-like substructure will mimic the effect of highellipticity sub-halos. Depending on the orientation, a disk can mimic the effect of line-like defects (vortices) or spherical sub-halos.
Using the software package PyAutoLens [39,40], we simulate the lensing of a galaxy at z = 1.0 by a spherical halo at z = 0.5 and Einstein radius of 1.5 arcseconds (mass M halo ∼ 5 × 10 11 M ), with and without a disk of mass 1%M halo and with orientation orthogonal to the line of sight. Figure 3 shows the fractional change in the lensing image introduced by the presence of the disc. The change ranges from a few percent to a ten-fold increase in the brightness, which indicates that gravitational lensing will be a powerful tool in distinguishing traditional particle dark matter from condensate models.
Finally, we consider the possibility that these disklike substructures could seed satellite galaxies, e.g. white dwarf galaxies and globular clusters. Typical scales for dwarf galaxies are a radius of ∼ 100pc and a mass < 10 10 M [41,42]. For example, the Hercules dwarf galaxy has a half-light radius of 350 pc and a mass of 7 × 10 6 M [41]. A disk solution of this mass and radius exists for m and λ satisfying m/eV = 1.2×10 12 λ 3/8 , with a range of masses and radii available by changing D and ψ(0). Thus these disk solutions can indeed provide satellite galaxies, however any estimate of the number density requires a more detailed study of the formation of such structures.
FIG. 3. Fractional change in brightness to the lensing from a spherical halo introduced by the presence of a disk of mass 1% that of the halo, cutting across the left side of the halo. The fractional change is defined as (with-without)/without.

VI. CONCLUSION
Observational evidence has shown our galaxy is out of equilibrium and has dark substructure [2,3]. With this in mind, we have argued that disk-like solutions should exist to the equations of motion describing a gravitating Bose-Einstein condensate, as emerge in models of dark matter involving ultra-light scalars and superfluids. This analysis indicates that dark matter substructure indeed provides a new opportunity to test the nature of dark matter.
We have not studied the evolution during or impact on structure formation, e.g. the direct seeding of baryonic disks or the disruption of dark disk substructure by tidal forces, the latter of which may result in interesting stellar streams. These are time-dependent problems, which will be studied in future work. Another interesting avenue, and potential smoking gun substructure of superfluid dark matter, is vortices [20]. These will generically be present inside the dark disk solutions if the disk has a net angular momentum. We leave this, and a complete mathematical analysis of the solutions presented here, to upcoming work. D r ∂φ ∂r whereφ is the complex conjugate of φ. We seek values of µ ∈ C for which nontrivial (ψ, w) can be found to satisfy (B2). System (B2) is in fact an eigenvalue problem, but the term D r ∂ ∂r can be difficult to deal with numerically. Therefore, we use the identity to introduce Φ(r) = r D 2 φ(r) and W (r) = r D 2 w(r) and transform (B2) to where denotes differentiation with respect to r. Separating Φ into real and imaginary parts by Φ = A + iB gives so that along with the boundary conditions W (0) = 0 and W (r) → 0 as r → ∞, we may solve for W to find This reduces reduces (B5) to solving which is now a proper eigenvalue problem for µ with eigenfunction [A, B] T . Implementing the above numerically can be obtained by considering some r * 1 to restrict r ∈ [0, r * ] and discretize the interval [0, r * ] into equally spaced points {r n } N n=0 so that 0 = r 0 < r 1 < · · · < r N −1 < r N = r * , where N 1 is chosen appropriately large. The boundary conditions on W are implemented numerically by the Dirichlet conditions W (0) = 0 and W (r * ) = 0, and using standard finite difference approximations we have that these boundary conditions give that the numerical discretization of the linear operator L is invertible. From here we again use the finite difference approximation of L to write (B6) as a matrix eigenvalue equation with A and B as vectors each of length N .