Interference of Dark Matter Solitons and Galactic Offsets

By performing numerical simulations, we discuss the collisional dynamics of stable solitary waves in the Schrodinger-Poisson equation. In the framework of a model in which part or all of dark matter is a Bose-Einstein condensate of ultralight axions, we show that these dynamics can naturally account for the relative displacement between dark and ordinary matter in the galactic cluster Abell 3827, whose recent observation is the first empirical evidence of dark matter interactions beyond gravity. The essential assumption is the existence of solitonic galactic cores in the kiloparsec scale. For this reason, we present simulations with a benchmark value of the axion mass $m_a = 2 \times 10^{-24}$ eV, which is somewhat lower than the one preferred for cosmological structure formation if the field is all of dark matter ($m_a \approx 10^{-22}$eV). We argue that future observations might bear out or falsify this coherent wave interpretation of dark matter offsets.


Introduction
The nature of dark matter is one of the most important open problems in fundamental physics. Projected experiments and astronomical observations are expected to shed new light on this question in the next decade [1].
In this context, the first evidence of dark matter (DM) non-gravitational self-interaction has been recently reported for the Abell 3827 cluster [2] (z ≈ 0.1), where a displacement of the stars with respect to the maximum density of its DM halo has been observed, for some of the merging galaxies [3]. Possible explanations for this offset within the ΛCDM model comprise casual allignment with other massive structures that might influence the results from gravitational lensing, astrophysical effects affecting the baryonic matter, tidal forces or simply wrong identification of lensed images [4]. Even if these causes cannot be fully excluded, meticulous observations and simulations have shown that any such interpretation is unlikely to explain the collected data [4,5]. This tension with collisionless dark matter models [5] suggests the necessity of considering other possibilities as, e.g. self-interacting dark matter, that yields a drag force slowing down the galactic DM distribution while leaving the standard model sector unaffected [3,4,5]. Nonetheless, requiring that the drag induces the offset implies a lower bound for the cross section that is in tension with upper bounds derived from other observations, as carefully discussed in [6]. Thus, the Abell 3827 cluster presents a challenging puzzle that opens up questions of crucial importance to understand the nature and dynamics of DM.
In this work, we address the problem of the measured offset using the scalar field dark matter (ψDM) model [7,8,9], which considers a Bose-Einstein condensate (BEC) of non-relativistic ultralight axions (ULAs) of mass m a subject to Newtonian gravity and that was introduced to solve difficulties of ΛCDM (e.g. missing satellites problem and cusp-core problem [10]), while maintaining the successful phenomenology of the model at cosmological scales [11,12]. Impressive numerical simulations [13] resolving largely different length scales have recently given support to this expectation. These extremely light scalar particles can arise in string theory constructions, e.g. [14] and other extensions of the standard model, e.g. [15]. Light scalars can also naturally appear as composites of hidden theories like the random UV field theory scenario [16].
We will show that the wave-like coherent nature of BECs severely affects the collisional dynamics of dark matter clumps, providing important effective forces even in the absence of explicit local interactions between the elementary dark matter constituents. We then discuss the possible relevance of this phenomenon to the puzzling observations described above.

Mathematical model
In the condensed scalar field scenario, the DM dynamics is governed by a Schrödinger-Poisson equation [17,18,19,20] for the mean-field wave-function ψ of the dark matter distribution: where |ψ| 2 is the particle number density, G the gravitational constant and t and x are time and position. For simplicity, we disregard cosmological evolution of the scale factor and the contribution of baryonic matter to the gravitational field, implicitly assuming that they do not play a prominent role in the processes studied below. Although a local interaction term λ|ψ| 2 ψ can be added to (1) [11,21,22], we will restrict ourselves to the simplest λ = 0 case [7,12,13] that, as we show below, is enough to describe the observed behaviour. Notice, however, that drag forces appear in similar mathematical models for optical systems with non-linear terms λ = 0, e.g. [23]. Equation (1) can be recast in terms of adimensional quantities: Following [24], the adimensional unit of length, time and mass correspond to: We have taken H 0 = 67.7km/(s Mpc) for Hubble's constant and Ω m0 = 0.31 for the matter fraction of energy today. Equation (3) yields localized, radially symmetric, self-trapped robust solutions which we will loosely call solitons. α is an arbitrary scaling constant, the propagation constant is β = 2.454α, the soliton mass is M sol = |ψ| 2 d 3 x = 3.883 √ α and its diameter (full width at half maximum) is d sol = 1.380/ √ α. f (.) and ϕ(.) are functions that can be computed numerically. In terms of dimensionful quantities, the mass and size of the solitons are related by: where M is the solar mass. In order to be reasonably self-contained, we give more details on these solutions and also discuss the numerical methods used for the computations in the appendix (section 7). Finally, let us remark that these stationary states have been independently discussed in several physical contexts: foundations of quantum mechanics [20,25], cold trapped atoms [26,27], QCDaxions [28] and ultralight DM [29,30]. This often overlooked formal coincidence indicates that studies concerning equation (1) can have deeply multidisciplinary implications.

Numerical simulations
In ψDM, galactic dark matter distributions consist of a core which can be identified with a soliton surrounded by a background also governed by Eq. (3) and evolving in time and space with uncorrelated phases [13,24,31].
In this work, we propose that the offset of Abell 3827 [3] can come from the repulsion between coherent DM clumps (the solitonic cores) in phase opposition, without any extra local interactions. We show by numerical simulations that destructive interference can provide a large effective force acting on the cores. This repulsion between robust wave lumps is well known in soliton systems, from nonlinear optics [32,33] to atomic physics [34,35,36], where the mathematical description of the phenomena is similar to the theory of coherent DM waves.
In ref. [3], observations of DM concentrations with mass of the order of 10 11 M surrounding stellar distributions separated by around 10 kpc were presented. This value of 10 11 M does not correspond to a galactic mass, but to clumps within the cluster that we will identify with solitonic cores. Most of the mass is in the halo, which behaves incoherently and therefore does not feel interferential forces. We will come back to this point in section 4. Taking the aforementioned values for M sol and d sol in equation (8), we find m a c 2 ≈ 2 × 10 −24 eV. We will fix this benchmark value for the simulations below. In section 5, we provide a discussion on previous observational constraints on m a and on their relevance to the phenomenon described here.
First, we have analyzed the collision of two DM solitons by numerically integrating (3) with the initial condition: where 2|x 0 | is the initial separation, 2v the initial relative velocity, ∆φ the relative phase and α is related to normalization as described in section 2 (adimensional units). Previous studies of this sort with ∆φ = 0 can be found in [22,37]. We use a split-step pseudo-spectral algorithm, known as beam propagation method [38,39] (see the appendix for technical details). It is worth quoting other powerful numerical methods that have been recently developed for the Schrödinger equation with nonlocal terms [40,41]. As expected, see e.g. [42] for a discussion in nonlinear optics with a particular nonlinear potential, the outcome largely depends on the relative phase and speed. In the case of phase opposition, destructive interference creates a void region between the solitons which can induce a bounce. For phase coincidence, the solitons merge into a single matter lump (which for large initial velocities eventually splits again). Interference fringes appear for large velocities [22,37].
We must underline that in this work, for the first time to our knowledge, the effect of coherent DM waves on luminous matter has been calculated, by adding to our simulations test particles following classical trajectories in the gravitational field generated by the DM wave. These particles, initially located at the soliton centers, are a toy representation of the stars and can be shifted from the DM density peaks in a collision, as we show in figure 1. In figure 2, we plot the comparison between the trajectories of the point particle and the DM projected mass maximum. In order to check the limitations of this particle model, we have made use of the well known fact that Schrödinger equation can be cast into a hydrodynamic form through the Madelung transformation [43]. This allows us to develop a fluid toy model in which luminous matter is described as a spatially extended cloud (see the appendix). As it can be seen in the inset of figure 2, both models display a good qualitative agreement.
Even if the collision in phase opposition is the simplest case, luminous vs. DM shifts can happen in more general situations. Figure 3 shows an example with four galaxies. Initial conditions are four solitons of mass M sol = 0.72M each, located at the vertices of a square of diagonal 40 kpc and initial velocities of 100 km/s towards the center, with phases 0, π/2, π and 3π/2, respectively. When the solitons approach each other, their phase gradients induce a rotation of the DM cloud, with ordinary matter lagging behind. It is worth mentioning that stationary rotating solutions of Eq. (1) have been discussed in [44,45].
Interference also plays an important role in asymmetric collisions if the phase difference between the lumps remains a well-defined quantity during the process. As an estimation, take |β 1 − β 2 |∆t col 1 where β 1 and β 2 are the propagation constants of each soliton and ∆t col is the duration of the collision. In terms of the soliton mass: This yields: In section 4, a collision involving four solitons will be considered. It is easy to check that the condition (11) holds in that case. The relative phase between two solitons is a decisive factor for their collisional dynamics, and, in the ψDM model, it is important for galactic mergers. This phase is ultimately determined by initial conditions at galactic formation. Moreover, the relative phase for any pair solitons changes in time since the propagation constant β depends on the mass, Eq. (10). In theory, given precise initial conditions, the cosmological evolution of the axion field can be computed deterministically in the semiclassical description of Eq. (1) [13]. In practice, uncertainties in the initial conditions and the evolution imply that the relative phase for a particular collision can actually be considered as random.

Comparison with observations
We now show that, starting with separate solitons, the wave dynamics of equation (1) can generate the gross features of the Abell 3827 cluster [3]: there are two DM blobs, one comprising galaxy N.1, for which dark matter and stars are separated; and the other one comprising galaxies N.2-N. 4. In  Figure 4 shows the result of a simulation. Initially, four separate solitons with masses 0.72, 0.95, 1.28 and 1.1 times 10 11 M are considered. Solitons 1 and 3 are heading soliton 2 with relative velocities of 220 and 180 km/s. Soliton 4 has an initial velocity of 900 km/s in the transverse direction, in order to agree with the redshift measured in [3]. After evolution, we find offsets similar to those displayed in [3]. Matching these qualitative features as in figure 4 obviously requires an appropriate choice of initial conditions but we remark that no special fine tuning is needed. Obviously, the real conditions are far more complicated. Apart from the coherent solitonic core, DM of field galaxies includes a non-coherent halo with an approximate Navarro-Frenk-White profile, see [24] for a detailed discussion. When the cluster is formed, most of its matter will be in an incoherent state with, at most, coherent lumps around the initial galactic cores (namely, around the stellar distributions). However, the presence of a large incoherent background does not necessarily change the qualitative features of the dynamics. Clearly, the effect of soliton-cluster halo interferences averages out to zero and can be neglected. Moreover, since the background density varies only mildly within the cluster, the gravitational forces it generates will not be dominant. A natural concern is whether incoherent matter might be attracted by the larger densities at the solitons, leading to smaller and more massive lumps. This is avoided if the kinetic energy of the incoherent wave is enough to impede its absorption, as in [24] for single galaxies. In fact, we have checked by numerical simulation that an incoherent background does not severely affect the process of figure 4.

Discussion on the axion mass
In this section, we briefly review the observational constraints on the axion mass m a (see e.g. [31] and references therein) and their relation to our dark wave interpretation of the offsets. The essential hypothesis for our modeling is the existence of kiloparsec scale coherent cores. We remark that the dark matter density distributions for R 5 − 6kpc of galaxies like the Milky Way are subject to large uncertainties [46], [47] and therefore the assumption is neither confirmed nor excluded by direct inference of the galactic profiles. It is natural to assume that we are in a scenario in which the cusp-core problem is solved solely by ψDM. By studying the Fornax dwarf galaxy in this context, the authors of ref. [13] found a best fit of m a c 2 = (8.1 +1.6 −1.7 ) × 10 −23 eV. A related analysis of Fornax and Sculptor in ref. [31] yielded a one-sided constraint m a c 2 < 1.1 × 10 −22 eV. See also [9] and references therein.
On the other hand, a lower bound m a c 2 > 10 −24 eV comes from requiring that ψDM is indistinguishable from ΛCDM for the probes studied in [48]. This is a conservative lower bound, derived only from linear constraints on the cosmic microwave background. More stringent but also more model dependent lower bounds were derived from nonlinear probes in [49,50,51], see also [9], [52] and references therein. Let us quote the result of [50], where it is found that data from the Hubble Ultra-Deep Field exclude axions with m a 10 −23 eV contributing more than half of DM. Signals from pulsars might soon give new information on the existence of ultralight axions and their mass [53].
The benchmark value of m a = 2 × 10 −24 eV that we have fixed in the simulations shown in the sections 3 and 4 comes from requiring solitonic cores with radius of the order of few kiloparsecs for masses of the order of 10 11 M . We allow ourselves to use this value of m a since it complies with the conservative lower bound of [48]. However, we envisage two possibilities in which these large cores could be present for larger values of m a : First, ULAs could be just a fraction of dark matter, relaxing to some extent the aforementioned stringent mass constraints, see e.g. [50]. Moreover, the total mass constituting each solitonic core would be smaller (for a fixed total dark matter mass), leading to larger radii by virtue of Eq. (8). The mechanism introduced in this paper can only cause a displacement of the ULA fraction of dark matter from the stars, but that can anyway render an offset for the DM center of mass.
Second, if a repulsive term λ|ψ| 2 ψ, as first introduced in [11,21], is added to Eq. (1), the soliton radius would be larger than the one given in Eq. (8), see [27] for a detailed discussion.
If future observations and/or analysis indicate that ψDM can only be realized with subkiloparsec cores for Milky Way-class galaxies, it would then be unlikely that soliton interactions play any role for providing relevant offsets within clusters like Abell 3827. Nevertheless, the analysis in this work would still play a role for the interaction between the solitonic cores. Understanding whether it might yield observational consequences is left for the future.

Conclusions
We have discussed the phenomenon of soliton interactions based on wave interference, which is relevant for any model of BEC dark matter relying on a Schrödinger-Poisson equation [7,8,9,28,54]. Large effective forces can be induced during collisions, in analogy with well known experiments in laboratory BECs and nonlinear optical systems.
If an ultralight scalar represents a significant fraction of DM, it is plausible that interference between dark waves can have observational consequences for galactic mergers and, in particular, it can explain the consequential results of [3]. The simplest setting for generating offsets is that of head-on collisions in phase opposition (see figures 1 and 2), but we stress that they appear in more general situations (figures 3 and 4). The paramount hypothesis is the existence of coherent cores with radii of the order of few kiloparsecs. There are important qualitative differences with other models of DM: the force acting on DM is between the solitonic cores and not between a galaxy and the cluster halo. Moreover, the outcome depends on the value of the relative phase at the moment of the collision, which in a realistic situation could be taken as random. These two points can be tested if other similar mergers are observed with the level of detail achieved by [3] and can potentially reconcile the offset in [3] with the lack thereof in other systems [55,56], which are in tension in models with particle-like interactions [6]. It is worth emphasizing that the dark matter shift discussed here is different from the one observed in other systems like the Bullet cluster [55], where the offset is between dark matter and gas, not stars, it is observed after the halos have traversed each other and is perfectly consistent with collisionless dark matter. A natural question is whether interference effects could then spoil the standard description of those systems. That is unlikely because the effective forces only act on the solitonic cores and are therefore confined to the kiloparsec scale or less. We expect the corrections to average out to zero in larger collisions like the Bullet cluster, with a size of a few megaparsecs. Direct numerical confirmation of this assertion is left for future work.
In the present work, we have considered an extremely simplified description of galactic dynamics which is enough to understand the gross features that can be expected from a wavelike behavior of DM. It would of course be desirable to incorporate these features in more detailed simulations as, for instance, those reported in [47] or [56].
Thus, we expect that through future theoretical progress and astrophysical observations, the scientific community will be able to discern the present scenario from models with explicit DM self-interactions or other logical possibilities. In a broader perspective, it is worth emphasizing that continuously improving observational evidence increasingly calls for precise descriptions of nonlinear phenomena. For instance [57] has studied the consequences of nonlinear evolution in the formation of cosmic voids. It is of great interest to understand whether alternatives to ΛCDM lead to differences that might be experimentally tested. Whether the discussion of the present contribution or the ψDM model in general may have implications in this respect is a compelling question for the future.
Finally, it is worth pointing out that, due to the formal coincidence of the governing equations, refined control of trapped atoms [26] and optical media [58] introducing gravity-like interactions might allow for laboratory analogue simulators of galactic-scale phenomena.

Appendix: Technical details on numerical methods
In this appendix we describe a number of technical issues related to the numerical treatment of the Schrödinger-Poisson equation. We will use the form of the equation in terms of dimensionless quantities (3).

Stationary solution and initial conditions
Consider an ansatz with radial symmetry ψ(t, x) = e iβt f (r), Φ(t, x) = ϕ(r) where we have defined r = |x|. Equation (3) is reduced to: Moreover, β can be reabsorbed asφ(r) = ϕ(r) + β. In order to find the soliton solution, we fix f (r = 0) = 1 and perform a standard shooting method by varyingφ(r = 0). Regularity at r = 0 ensures that there are no more free parameters andφ(r = 0) is fixed by requiring that f (r) vanishes as r → ∞. Since lim r→∞ ϕ(r) = 0, the value of β is read from the large r behavior ofφ(r). We find that β = 2.454, M sol = 4π ∞ 0 r 2 f (r) 2 dr = 3.883 and the full width at half maximum of the density is fwhm≡ d sol =1.380. The solution is depicted in figure 5. Notice that ϕ(r) has been rescaled in order to refer both functions to the same axis (ϕ(r = 0) ≈ −4.76). Due to the symmetry of equation (3), a whole family of solutions can be found by scaling. If we take |ψ|(r = 0) = α for any α > 0, then there is a solitary wave solution with M sol = 3.883 √ α, β = 2.454α and fwhm= 1.380/ √ α. Once a stationary solution is found and the functions f (r), ϕ(r) are explicitly known, Galilean invariance of the Schrödinger-Poisson equation allows us to write down the general boosted solution representing a soliton moving with constant velocity v: This solution is used to define the initial conditions for the simulations, in which we consider initially separated solitons. For instance, the soliton collision of figure 1 has initial condition: (14) where (−x 0 , v), (x 0 ,−v) are the (adimensional) initial position and velocity of each soliton, ∆φ their relative phase and 3.883 √ α the adimensional soliton mass. The test particles modeling the stars are initially placed at the center of the solitons with the same initial velocity. It is worth noticing that, due to Galilean invariance of equation (3), the same process can be thought of as, for instance, one soliton initially placed at −2x 0 moving with velocity 2v toward a static soliton with center at x = 0. However, it is necessary to change the initial phase accordingly, namely: Similarly, the initial condition used for the simulation leading to figure 4 is, explicitly:

Dynamical evolution
In order to compute temporal evolution in (3), we have used a split-step Fourier method, widely used in the integration of the nonlinear Schrödinger equation because of its stability and precision, see for instance references [38] and [39]. Schematically ψ(t + dt, where F is the three-dimensional Fourier transform. The preservation of the norm |ψ| 2 d 3 x is automatic. At each step, we need to compute Φ(t, x) from Poisson equation. We do so using a finite difference scheme in a N x ×N y ×N z spatial grid. The discrete Poisson equation can be written as a linear problem A · B = C where A is the corresponding heptadiagonal N x N y N z × N x N y N z sparse matrix, B is a N x N y N z × 1 matrix with the value of the function Φ at each point of the grid and C includes the source and the boundary terms. This algebraic problem can be efficiently solved by an iterative symmlq method. We introduce boundary conditions as if all the mass were concentrated at the center of the computational grid (we compute in a reference frame in which the center of mass coincides with the center of the grid and the total momentum vanishes). The error introduced by these boundary conditions becomes negligible as the size of the computational box becomes much larger than the size of the region of interest. As an additional cross-check, we have also used a second method for solving Poisson equation, namely that of directly using Fourier transformation to deal with the laplacian. That implies periodic boundary conditions for Φ, which also approach the physical boundary conditions as the box is made larger. In order to compute the classical trajectories for the point particles representing the standard model matter, we use Heun's algorithm (a second order Runge-Kutta scheme), making use of the gravitational potential computed at each time step. In order to test the precision of the algorithm, we first track the solution corresponding to a soliton moving with constant velocity, see figure 6, where the evolution of the modulus of the wave function at the center of the soliton is plotted. Numerical errors introduce two types of oscillations around the theoretical constant value. The short period oscillation is due to the velocity whereas a fluctuation with a a longer period appears also for v = 0. In any case, the figure shows that the method is stable and that fluctuations can be kept small.
In order to further check the accuracy and validity of the computational methods, we have performed a series of auxiliary numerical simulations. In particular, we have calculated the effect on the dark matter soliton motion of different numerical schemes for integrating Eq. (3). Results for the example corresponding to figure 1 (∆φ = π) are shown in fig. 7, where we plot the evolution of the DM soliton clouds and the ordinary matter and the corresponding offset. The maximum density position of the dark matter distribution is found by quadratic interpolation around the maximum value in the discrete grid. As it can be seen in the pictures, the differences between the methods are negligible for the evolution times that we have used in the paper.

A fluid toy model for ordinary matter
In this paper, we have used the simplest possible toy model for the stars, considering just test particles. More realistic descriptions would take into account that luminous matter is not pointlike and that it sources the gravitational potential. However, these effects should not affect the qualitative conclusions presented above, since the offsets are generated by the interference force acting on the dark matter solitonic cores while not affecting the ordinary matter. As a first test of this assertion, we have repeated the simulation for the collision in phase opposition (figs. 1 and 2) by considering the, arguably, simplest fluid model for ordinary matter: each galaxy is modeled by an independent Schrödinger equation coupled to Eq. (3). Thus, i ∂ t g i = − 1 2γ ∇ 2 g i + γΦg i , where g i correspond to the luminous matter distribution of each galaxy and γ is a parameter controlling its size (γ = 2 was taken for the figure). The gravitational potential Φ is still determined by Eq. (3), being the ordinary matter distributions g i considered as test fields. It is well known that nonlinear Schrödinger equations can be recast as fluid equations using a Madelung transformation (see, e.g. [43]). Results are displayed in figure 8. showing a qualitative behavior similar to the point particle model. For the simulation, we used a system of two mutually incoherent nonlinear Schrödinger equations, coupled to the system of Eq (3) given in the paper. The trajectory of the center of gravity of the luminous matter has been tracked and compared with the evolution of the "point galaxies" in the inset of fig. 2, with good qualitative agreement. Red lines correspond to dark matter projected density contours.