Merger of galactic cores made of ultralight bosonic dark matter

TL;DR

This paper investigates the dynamics of binary mergers of ultralight bosonic dark matter cores, revealing complex oscillatory behaviors and non-stationary final states through numerical simulations of the Gross-Pitaevskii-Poisson system.

Contribution

It provides a detailed analysis of the merger process, including gravitational cooling and oscillation characteristics, which were not previously comprehensively studied for ultralight bosonic dark matter cores.

Findings

Final configurations oscillate around a virialized state.

Oscillation amplitude and frequency depend linearly on impact parameter.

Mergers do not lead to a clear stationary state, even after relaxation.

Abstract

We study binary mergers of ultralight bosonic dark matter cores by solving the Gross-Pitaevskii-Poisson system of equations. The analysis centers on the dynamics of the relaxation process and the behavior of the configuration resulting from the merger, including the Gravitational Cooling with its corresponding emission of mass and angular momentum. The oscillations of density and size of the final configuration are characterized, indicating that for the equal mass case the dependency of the amplitude and frequency of these oscillations on the impact parameter of the pre-merger configuration is linear. The amplitude of these oscillations changes by a factor of two or more indicating the final configuration does not approach a clear stationary state even though it oscillates around a virialized state. For the unequal mass case, global quantities also indicate the final configuration oscillates around a virialized state, although the density does not show a dominant oscillation mode. The evolution of the angular momentum prior and post merger is analyzed in all cases.

Tables

Table 1. Table 1: Fitting parameters of ρ ​ ( x , 0 , 0 ) 𝜌 𝑥 0 0 \rho(x,0,0) for the case M ​ R = 1 𝑀 𝑅 1 MR=1 , p x ​ 0 = 0.1 subscript 𝑝 𝑥 0 0.1 p_{x0}=0.1 and three values of the impact parameter y 0 subscript 𝑦 0 y_{0} .

$y_0$	${\rho }_0$	$r_c$	$t$	comment
1	8.280	0.647	351	density at a maximum
	2.684	0.978	358	density at a minimum
5	6.862	0.776	348	density at a maximum
	2.566	0.982	355	density at a minimum
10	5.053	0.852	254	density at a maximum
	2.571	4.023	261	density at a minimum

Equations

$$i{\partial }_t\Psi$$

$${\nabla }^{2}V$$

$$M$$

$$K$$

$$W$$

$$L_z$$

$${\rho }_{soliton}(r)={\rho }_0{\left[1+0.091{\left(\frac{r}{r_c}\right)}^2\right]}^{-8}$$

Full text

I Introduction

One of the viable dark matter candidates currently under study is the ultralight spin-less boson Matos-Urena:2000 ; Ostriker:2016 , which is attractive because of some interesting properties consistent with observations. For instance when its mass is of order $m\sim 10^{-22}eV$ structures do not develop cusps due to the large de Broglie length Chen:2016 ; Schive:2015 ; Du:2016 ; Velmaat2018 whereas at large scale the behavior is consistent with that of CDM Schive:2014 ; Schive:2014hza . At the same time, this model is also consistent with the small structure abundance of the mass power spectrum Matos-Urena:2000 ; Ostriker:2016 ; Marsh-Ferreira:2010 ; Schive:2015 .

Local scale dynamics on the other hand, should indicate differences between CDM and ultralight bosoinc dark matter and impose constraints on the later. For instance, the relaxation process should be special, being the gravitational cooling process an option GuzmanUrena2003 ; GuzmanUrena2006 in which matter carries out kinetic energy leaving the structure under relaxation in a nearly virialized state, or other processes involving dynamical friction Mocz2017 , or damping ChavanisDamping could provide the relaxation mechanism. Also the collisions and interaction between structures can provide important restrictions to the model, for example the density resulting from head-on core mergers AvilezGuzman2019 that may result in the destruction of luminous matter clusters during the process for certain particular scenarios GonzalezGuzman2016 . Other restrictions, this time on the boson mass are found from the analysis of core oscillations that may or may not allow the formation of star clusters in galaxies Marsh2018 .

Locally, the dynamics of this dark matter model is ruled by the Gross-Pitaevskii-Poisson (GPP) system, that describes the evolution of a Bose-Einstein Condensate in the Gross-Pitaevskii mean field approximation, contained by the gravitational potential generated by itself. One point the various studies and approaches at local scale of the model have in common, is that this type of dark matter clumps into structures with a universal profile, either into an equilibrium configuration of the GPP system for isolated systems GuzmanUrena2006 ; BernalGuzman2006b , or composed of a core, sometimes called solitonic profile that matches the density profile of an equilibrium configuration Schive:2014 , and a surrounding cloud with a NFW profile obtained from simulations involving structure formation clustering Schive:2014 ; Schive:2014hza ; Schwabe:2016 ; Mocz2017 ; Niemeyer2016 ; ChavanisCore .

Among the common interactions between structures or cores, the merger of two of them is very important and is the subject of this paper. Configurations resulting from a merger with angular momentum, naturally inherit rotation from the original merging cores. Rotating structures within this dark matter model are interesting for various reasons. One of them is that rotation is an extra parameter for BEC dark matter halos that help fitting galactic rotation curves by keeping the boson mass unchanged BEC2014 ; BEC2015 , and will possibly help to reduce the dispersion of boson mass in rotation curve fitting of big catalogs Tula . In a similar context, ellipsoidal analytic solutions to the GPP with rotation have been associated with possible vortex solution RindlerShapiro2012 . And more recently, new exact solutions of the GPP system with rotation are also being constructed with the aim of studying this dark matter model at local scale DavidsonSchwetz2016 ; Hertzberg2018 .

The study of core mergers in orbit or during structure formation, within the context of ultralight bosonic dark matter is not new. In fact also multiple soliton mergers have also been studied Schwabe:2016 ; ParedesMichinel2016 ; Mocz2017 . Specially in Schwabe:2016 the mergers have been analyzed in detail, from the initial conditions to the properties of the final configuration. Among the most interesting results it was found that the final mass of the merger does not depend on the initial momentum of the orbiting objects and only depends on mass ratio, total initial mass and total energy of the system. Also in Niemeyer2016 the density of cores resulting from mergers is compared with the solitonic profile in the context of structure formation simulations.

The analysis in our paper is very similar to that in Schwabe:2016 , however some new results arise. Important differences are that we solve the GPP system without using the Madelung transformation, not for calculations nor for diagnostics of macroscopic quantities. We in fact confirm that the final mass of the merger does not depend on the initial angular momentum of the pre-merger configuration, however we find this result holds only for the equal mass case. We also find that the angular momentum of the final configuration depends on the initial conditions prior to merger, for both the equal and unequal mass cases.

On the other hand, in Velmaat2018 within the analysis of structure formation it is found that cores exhibit strong undamped oscillations. Our results are consistent with this evidence. From our analysis, we find that the configuration resulting from the merger of two cores exhibits a dynamical behavior, characterized by oscillations with considerable amplitude that depend on the parameters of the binary system. The final structure does not relax, however by fitting the density profile at different times we illustrate how the core radius and central density change in time.

The paper is written with the following structure. In Sec II we describe the method used to simulate the mergers. In Sec III we analyze the equal and unequal mass scenarios. In Sec IV we draw some conclusions.

II Evolution of the system

Like in the analyses of structure formation and binary mergers mentioned before, we assume the dynamics of the ultralight bosonic dark matter is ruled by the GPP system of equations. Likewise we assume the free field regime, where the self-interaction among bosons is neglected, the so called fuzzy dark matter regime. Finally, we solve the equations using numerical methods and initial conditions described below.

II.1 Numerical methods

We solve the time dependent GPP system of equations which in code units is written as

[TABLE]

that describes the evolution of the fuzzy dark matter. Here, $\Psi$ represents the wave function of the system and $|\Psi|^{2}$ is interpreted as the macroscopic density of the condensate and $V$ is the gravitational potential sourced by the condensate itself. We solve these equations for $\Psi$ in a cubic finite domain, with initial data for $\Psi$ consistent with the potential $V$. In this system, Poisson equation is a constraint that has to be solved on the fly as the bosonic gas density evolves.

We solve the Gross-Pitaevskii equation numerically in 3D using the method of lines for the evolution across spatial slices separated by intervals of time $\Delta t$. The spatial domain $D=[x_{min},x_{max}]\times[y_{min},y_{max}]\times[z_{min},z_{max}]$ is described with a Cartesian and uniformly discretized grid defined by $x_{i,j,k}=x_{min}+i\Delta x$, $y_{i,j,k}=y_{min}+j\Delta y$ and $z_{i,j,k}=z_{min}+k\Delta z$, for $i=0,...,N_{x}$, $j=0,...,N_{y}$, $k=0,...,N_{z}$, with an isotropic resolution $\Delta x=\Delta y=\Delta z=(x_{max}-x_{min})/N_{x}$.

We discretize the equations with second order accurate finite difference stencils for spatial derivatives. For the sake of accuracy in the region of the merger, we use fixed mesh refinement based on the Berger-Oliger algorithm BergerOliger , with concentric refinement boxes. The resolution factor between successive refinement levels is one half. Considering that for the stability of the evolution, time and space resolution are limited by the condition $C=\Delta t/\Delta x^{2}<0.25/\sqrt{3}$, we choose the value of $C$ to be that corresponding to the most refined level.

We solve Poisson equation for $V$ with a Multigrid algorithm with subcycles that use the Successive Over Relaxation method. This equation is solved at initial time and during the evolution. Due to its computational cost, the integration of Poisson equation represents the major bottleneck of the code during the simulations.

Since we want to avoid reflections of matter from the boundary of the numerical domain, and because the Gravitational Cooling depends on the emission of matter that carries kinetic energy with it, we implement a sponge consisting of the addition of an imaginary potential such that $V\rightarrow V+V_{im}$, acting as a sink of particles following the recipe in (GuzmanUrena2004, ). We make sure that the transition region of the sponge lies exclusively in the coarsest refinement level.

II.2 Initial conditions

We assume the colliding objects are equilibrium configurations, which are spherical stationary solutions, constructed by assuming a harmonic time dependence of the wave function $\Psi=e^{-i\omega t}\psi(r)$, where $\omega$ is the eigenvalue of the Sturm-Liouville problem resulting from the spatial and time symmetries of $\Psi$ as described in GuzmanUrena2004 .

The initial wave function for the collision of two configurations is the superposition of the wave functions of two of these equilibrium configurations with different masses and linear momentum. For the superposition we use the method in AvilezGuzman2019 , specifically, we do not solve the Sturm-Liouville problem for two equilibrium configurations with different masses. Instead, we exploit the scale invariance of the GPP system of equations GuzmanUrena2004 , namely that by scaling physical quantities as $t=\lambda^{2}\hat{t}$, $x=\lambda\hat{x}$, $\Psi=\hat{\Psi}/\lambda^{2}$, $V=\hat{V}/\lambda^{2}$, where $x$ represents any of the spatial coordinates and $\lambda$ is a number, the GPP system (1) remains unchanged. Thus, the solution of the GPP system for a given configuration means one can construct all other equilibrium configurations using this scaling property. In fact, a consequence of this scaling is that density and mass also scale as $\hat{\rho}=\lambda^{4}\rho$, $\hat{M}=\lambda M$ which are important physical parameters of a scaled configuration used below.

In practice the typical equilibrium configuration is that with the central value of the wave function $\psi(r=0)=1$ that we will call $\psi(r)_{1}$ and has mass we call $M_{1}$. We choose one of the two configurations that will collide, to be precisely this standard configuration.

The second configuration that will collide, is a configuration constructed with the scaling relations above, represented by the wave function $\psi(r)_{\lambda}=\lambda^{2}\psi(r)_{1}$, with mass $M_{\lambda}=\lambda M_{1}$. Notice that the scaling parameter happens to be the mass ratio $\lambda=M_{\lambda}/M_{1}=\lambda=MR$ between the first and the second configurations used for the collision. In the analysis we consider the convention $0<\lambda<1$ in all cases, so that $M_{\lambda}<M_{1}$ always.

We then interpolate and superpose the two configurations in the numerical domain $D$. In order to maintain the system evolving within the numerical domain, we set the center of mass of the configuration at the coordinate origin. We parametrize the initial conditions by fixing the coordinates of the lighter configuration with mass $M_{\lambda}$ at $(x_{0},y_{0},0)$ with $x_{0},y_{0}>0$. Then, in order for the center of mass to lie at the origin, the center of the heavy configuration with mass $M_{1}$ must be centered at coordinates $(-\lambda x_{0},-\lambda y_{0},0)$. In this set up $y_{0}$ will play the role of impact parameter prior to merger.

The angular momentum is added through the imprint of linear momentum to the configurations along the $x$ direction only. For this we parametrize the momentum with the $x-$component of the heavy configuration with mass $M_{1}$ that we set to $p_{x0}$. Then again, in order to keep the center of mass approximately at the coordinate origin, the momentum of the light configuration must be $p_{x0}/\lambda$. The momentum is applied to each of the configurations by redefining $\psi_{1}\rightarrow e^{ip_{x0}x}\psi_{1}$ and $\psi_{\lambda}\rightarrow e^{-ip_{x0}x}\psi_{\lambda}$. Finally the wave function of the binary system at initial time is $\Psi=\psi_{1}+\psi_{\lambda}$ and the scheme in Fig. 1 illustrates the initial conditions.

II.3 Diagnostics

We monitor the dynamics of the system by evaluating some macroscopic variables. These include the mass $M$, kinetic energy $K$, gravitational energy $W$ and the $z$ component of the angular momentum $L_{z}$. These quantities are

[TABLE]

where the integrals are calculated using the second order accurate trapezoidal rule. A first important quantity is the total energy $E=K+W$, whose sign determines when a system is bounded ($E<0$) or unbounded ($E>0$). A second one is $Q=2K+W$ which should be zero for a virialized system and allows one to determine when a system is near, tends to or is far from a virialized state.

III Analysis

III.1 Parameter space

There is a garden variety of possible configurations that can be explored. However three parameters influence the behavior of the configuration resulting from the interaction between the two cores, namely, mass ratio $MR$, momentum $p_{x0}$ and parameter impact. These three parameters determine wide ranges of angular momentum and total energy values at initial time. It would be ideal to have the possibility of exploring a wide range of this parameter space. Nevertheless, due to the expensive computational cost of simulations, we restrict the exploration to illustrate the influence of some parameters using specific values.

First, we set two possible values of the mass ratio $MR=\lambda=0.5,1$, which are the equal mass scenario and the two to one mass ratio case, which will illustrate well the behavior of the system in unequal mass encounters.

Second, we consider various values of the impact parameter $y_{0}$. The radius of the configuration with mass $M_{1}$ is $r_{95}\sim 3.93$ in code units GuzmanUrena2004 , whereas that of mass $M_{1/2}$ is twice as big. Thus we study the range of values $y_{0}=1,2,...,10$ that accounts for scenarios ranging from nearly head-on to a separation various times bigger than the size of the structures.

Third, we distinguish between merger and unbounded scenarios. In the first scenario the two configurations end up together and form a final configuration. In the second scenario, either the configurations flyby each other or behave as solitons. The momentum $p_{x0}$ is useful to generate the two scenarios, because it sets the amount on kinetic energy $K$ of the two configurations together. With a low value of this parameter the gravitational energy $W$ dominates, implying that $E<0$, otherwise a high momentum contributes to $K$ that can contribute importantly to the energy to be positive and produce unbounded configurations. The threshold value for the head-on scenario is found to be $p_{x0}\sim 0.7$ BernalGuzman2006 which serves as a guide to avoid non-merging cases.

We empirically found a range of values of $p_{x0}$ for which at initial time the total energy is negative for the two values of $MR$ and all the values of $y_{0}$. Values in the range $p_{x0}\in[0,0.3]$ produce configurations with negative energy. In what follows we use the case $p_{x0}=0.1$ to illustrate the generic properties of mergers.

The values of these physical parameters suggest the numerical parameters to be used. The first parameter is the location of the lighter configuration at $(x_{0},y_{0},0)$ with $x_{0}=10$ in all cases. We use this value because the interference at the origin between $\psi_{1}$ and $\psi_{\lambda}$, $\langle\psi_{1},\psi_{\lambda}\rangle$ is less than $10^{-8}$. We consider the domain to be the box $D=[-40,40]^{3}$ and cover it with two refinement levels, and maximum resolution $\Delta x=0.1r_{95}$ in the inner box, which covers the region where the dynamics is more important $D_{h}=[-20,20]^{3}$.

III.2 Global quantities

In a merger scenario the two cores collide and form a single final configuration whose density profile can eventually be fitted with a simple function that can be further used to understand and analyze the physics of different processes.

We study now this scenario using $p_{x0}=0.1$ for the two values of $MR$ and all the values of the impact parameter $y_{0}$. The system of equations (1) is solved numerically for the initial conditions described above and we show the evolution of some of the scalars defined in Sec. II.C in Fig. 2, for the ten values of $y_{0}=1,...,10$.

The energy $E$ is shown normalized with the absolute value of its initial value. Notice that the total energy becomes more negative than at initial time, which indicates that the gravitational energy plays a more important role with time. The energy is also lost in a bigger proportion for smaller impact parameter $y_{0}$, and for $MR=1$ than for $MR=0.5$.

The mass is normalized with the mass of the standard equilibrium configuration $M_{1}$, therefore for $MR=1$ the total mass is initially $M=2$, whereas for $MR=0.5$ the initial mass is $M=1.5$. Notice that the mass decreases because matter is ejected and eventually captured by the numerical sponge. The combination of these two observations indicates that the mass lost during the process carries kinetic energy with it, exemplifying the Gravitational Cooling process SeidelSuen1991 ; GuzmanUrena2006 .

Notice also that for $MR=1$, the total mass is higher at initial time, but is also lost in a bigger percentage compared to the case of $MR=0.5$. It can also be seen that the bigger the impact parameter $y_{0}$, the smaller the mass ejected during the process. For $MR=1$ the final mass converges to the same value independently of $y_{0}$, or equivalently to the initial angular momentum of the pre-merger configuration as discovered in Schwabe:2016 . Nevertheless for $MR=0.5$ this is not the case, at least within the time window of our simulations.

Another interesting result is that the matter also carries angular momentum with it. In the bottom panel of Fig. 2, the proportion of angular momentum during the merger is shown.

The evolution of angular momentum shows an interesting behavior. For $MR=0.5$ the amount of $L_{z}$ released is between $\sim 20$% for $y_{0}=1$ and $\sim 40$% for $y_{0}=10$. In this sense, the simulations indicate that the merger process can produce final configurations with a wide range of values of angular momentum that could give origin to rotating galactic cores. However, for $MR=1$ the loss of angular momentum radiated away is of $\sim 65\%$ for $y_{0}=10$ and even turns negative for $y_{0}=1,2,3$ in the time window of the simulations and could hold also for other values in a bigger time domain. This result is interesting and the reason for the change of sign is that small values of $y_{0}$ correspond to nearly head-on situations. Since $y_{0}$ is the value of the impact parameter only at the center of each configuration, part of the matter ejected should be that initially located farther from the $x-$axis which carries angular momentum with it when it abandons the domain. This turn in the direction of rotation could be an interesting sign that eventually may provide restrictions to the model or predictions.

III.3 Equal mass case

The evolution of a specific simulation is shown in Fig. 3 for the equal mass case. The final configuration remains centered at the coordinate origin, rotates and has an ellipsoidal density profile. Animations of this and cases with various other parameter values are available in the supplemental material suppl .

In order to learn more about the dynamical behavior of the final configuration, we track the value of the central density and $Q=2K+W$ as functions of time that are shown in Fig. 4 for the two extreme values of the impact parameter $y_{0}=1,~{}10$. It can be seen that the quantity $Q$ oscillates around zero with a decreasing amplitude as expected for the Gravitational Cooling AvilezGuzman2019 .

The central density oscillates changing values by factors between two and three in the nearly head-on case $y_{0}=1$ and smaller oscillations for $y_{0}=10$. Fig. 4 suggests that the amplitude of the oscillations and the central value of the density depend on the impact parameter. In order to find a dependency on $y_{0}$ we calculated the average density $\rho_{avg}$ and its standard deviation to have a measure of the amplitude variation around the average $\rho_{dev}$, for $t>200$. The results are shown in Fig. 5, which suggest that both, the central density and oscillation amplitude depend on $y_{0}$ linearly. Finally, calculating a Fourier Transform within the same time domain, we obtain the peak frequency associated to the dominant density oscillation mode, which also depends on the impact parameter as shown in the third panel of Fig. 5. Knowing that the final mass is the same for all values of $y_{0}$, the oscillation frequency is genuinely different for different values of $y_{0}$.

In structure formation simulations Schive:2014 ; Schive:2014hza ; Schwabe:2016 ; Niemeyer2016 ; Mocz2017 the density distributions resulting from the interaction of two or more configurations are associated to density profiles with a solitonic core and a tail, however it is not quite specified whether these are final, relaxed configurations or not. As far as we can tell, the oscillations shown in Fig. 4 do not correspond to a relaxed structure. Even though the density profile can be fitted with the core profile

[TABLE]

as indicated in Schive:2014 ; Mocz2017 , where $\rho_{0}$ is the central density and $r_{c}$ is a core radius.

The issue is that the density is oscillating with considerable amplitude as seen in Fig. 4. Nevertheless, the fitting was performed on the density profile when the central density is at a local maximum and at a local consecutive minimum. The fitting parameters results appear in Table 1 for the projection of $\rho$ along the $x-$axis. Notice that the central density changes by a factor between two and three from a minimum to a maximum, whereas the core radius changes by nearly 50%. An example of how the density profile changes in time is illustrated in Fig. 6, where the projection of the density along $x$ is shown at two specific times for $y_{0}=10$ at a minimum ($t=254$) and at a maximum ($t=261$).

In order to have an idea of the physical time scale of these oscillations, we use the recipe in Mocz2017 . Considering a boson mass value $2.5\times 10^{-22}$eV and that core radius of the final configuration is converging to $r_{c}=1$kpc, using the range of frequencies $\nu\in(0.084,0.1)$ from Fig. 5, the period of the density oscillations is in the range $T\sim 0.76-0.91$Gyr. If the core radius is considered to be $r_{c}=0.25$kpc the period is within the range $T\sim 47-57$Myr.

III.4 Unequal mass case

As described before, the case we look at in detail corresponds to $MR=0.5$. The time dependence of $M$, $Q$ and $L_{z}$ appears in Fig. 2. General properties are very similar to those of the equal mass case. The loss of mass and angular momentum is smaller when the impact parameter is bigger.

Snapshots of the unequal mass $MR=0.5$ merger with $y_{0}=7$ are presented in Fig. 7. The resulting high density region wobbles around the origin due to the asymmetric distribution of matter and at some point evolves toward the coordinate origin. Animations for other values of the parameters are also shown in the supplemental material suppl .

What is different from the equal mass case is the relaxation process. The evolution of $Q=2K+W$ and the central value of the density are shown in Fig. 8 for the two extreme values of the impact parameter $y_{0}=1,10$. The value of $Q$ oscillates around zero with amplitude an order of magnitude smaller than in the equal mass case. For density on the other hand, since the configuration is wobbling around the coordinate origin, instead of tracking the central value of the density we track its maximum value $\rho_{max}$. The result in the Figure is generic behavior for the unequal mass cases with values of $MR$ between 0.5 and 1 we experimented with. The highly dynamical behavior is due to the fact that the small configuration with mass $M_{\lambda}$ approaches with a higher velocity and the distribution is much less symmetric than for $MR=1$. This explains a quick ejection of kinetic energy so that $Q$ acquires small values.

The density does not show any clear sign of relaxation or a particular dominant mode during the time window used in the simulations. This is perhaps a major obstacle when the density is fitted with a space-dependent density fitting function. Unlike the equal mass case, where the average of the density is a good estimate of the asymptotic value, here the expected value of the central density is uncertain. Nevertheless, Fig. 8 indicates that the central density of the final configuration depends on the impact parameter $y_{0}$.

IV Conclusions and discussion

We have presented the merger process of ultralight bosonic dark matter cores, with detailed illustrations of the equal mass case $MR=1$ and a representative unequal mass case $MR=0.5$.

In the equal mass case it was found that the final configuration oscillates with amplitudes that depend on the parameters of the binary prior to merger, namely, the mass ratio of the two initial cores, linear momentum and impact parameter. The resulting final configuration was fitted with a solitonic density profile at different times during the relaxation process. It was found that the density may change by factors of nearly three whereas the core radius can change by nearly 50% percent, and that the amplitude and frequency of the oscillations can be linearly related to the impact parameter of the merger.

In the unequal mass case, due to the size of the initial configurations, the interference becomes important in the symmetry of the final high density zone, which wobbles around the center of mass before it settles toward a nearly fixed location. The density in this case oscillates, however with an irregular superposition of modes, although with values of $Q$ indicating that globally the system evolves around a virialized state.

In both scenarios, it calls the attention the fact that the density seems far from a stationary state in cosmological time scales. The reason is that the amplitude of oscillations of the configuration resulting from a merger is not small, and perhaps it would be useful to consider time averages in such fittings.

In order to determine observational restrictions of this dark matter model, it seems unavoidable to systematically analyze the effect of the dynamics of a configuration resulting from a merger on the luminous matter that can be involved in the process. For example their survival questioned for specific scenarios of the head-on case in GonzalezGuzman2016 or restrictions from the existence of star clusters near galactic cores Marsh2018 .

Acknowledgments

This research is supported by Grants CIC-UMSNH-4.9, CIC-UMSNH-4.23, CONACyT 258726 (Fondo Sectorial de Investigación para la Educación). The runs were carried out in the computer farm funded by CONACyT 106466 and the Big Mamma cluster at the IFM.