Analyzing planar galactic halo distributions with fuzzy/cold dark matter models

We perform a numerical comparison between the fuzzy dark matter model and the cold dark matter model, focusing on formation of satellite galaxy planes around massive galaxies. Such galactic dynamics with controlled initial subhalo configurations are investigated using GADGET2 for the cold dark matter and PyUltraLight for the fuzzy dark matter, respectively. We demonstrate that satellite galaxies in the fuzzy dark matter side have a tendency to form more flattened and corotating satellite systems than in the cold dark matter side mainly due to the dissipation by the gravitational cooling effect of the fuzzy dark matter. Our simulations with the fuzzy dark matter typically show the minor-to-major axis ratio c/a of the satellite galaxy planes to be 0.21 ∼ 0.30; this well matches the current observed value for the Milky Way.


Introduction
Numerical simulations based on the cold dark matter (CDM) model successfully reproduce the observed large-scale structures of the universe. However, it encounters difficulties in explaining small-scale structures around galactic scales. For example, the CDM model predicts a cusped central halo density which is not observed and too many dwarf galaxies around massive host galaxies than the observed [1][2][3][4].
On the other hand, the satellite-plane problem is another serious small-scale challenge for the CDM model (see [5,6] for a review). Observations indicate that a substantial fraction of satellite galaxies around the Milky Way, Andromeda and Centaurus A co-orbit in flattened planar structures [7,8]. At least 10 highly flattened planes of dwarf galaxies have been discovered [9]. This is in sharp contrast with cosmological simulations with the CDM, which show much more random distributions and motions of satellite galaxies. The probability of finding such flattened and co-rotating satellite systems with the CDM is very low [5,7] although the statistic is controversial.
It is not easy to explain these satellite planes theoretically in the CDM model. Disks in astronomical objects such as a protoplanetary disk or a black hole accretion disk form when falling matter such as gas loses its energy by a dissipation (cooling) mechanism while keeping the angular momentum of the matter. For a given angular momentum, the lowest energy state is a disk orthogonal to the rotation axis [10]. Similarly, to form a flat disk-like structure small galactic halos should lose their gravitational energy by some dissipation mechanism when they fall towards heavy central galaxies. However, the typical collisionless CDM model lacks an efficient mechanism for this to happen. Due to the large scales (O(10 2 kpc)) baryonic processes could not strongly affect the satellite planes.
As an alternative to the CDM, there is a growing interest in the fuzzy dark matter (FDM) [11][12][13][14][15], also known as ultra-light axion or scalar field dark matter. In this model, dark matter (DM) particles satisfying the Schrödinger-Poisson equation have an ultra-light mass m 10 −22 eV/c 2 and a long de Broglie wavelength λ = 2π /mv kpc, where v is the velocity of the dark matter particles. (See [16,17] for more references.) Unlike CDM particles, condensed FDM particles behave as a coherent wave, and due to the quantum -1 -

JCAP12(2022)033
uncertainty principle the FDM has a typical length scale λ which suppresses the formation of DM structures smaller than a galaxy core and helps resolve the aforementioned small-scale issues [16]. While the CDM model has no efficient dissipation mechanism to make the observed satellite systems, a gravitational cooling effect in the FDM model provides a unique and efficient dissipation mechanism for them. The gravitational cooling [18][19][20] is a mechanism for relaxation by ejecting part of the FDM density and carrying out excessive kinetic energy and momentum. In this mechanism colliding FDM halos can have interfering wave profiles which contain high momentum modes. These modes can easily escape the gravitational potential of the halos carrying out energy.
The main aim of this paper is to show that halos in the FDM model have a tendency to form a more flattened and corotating satellite plane than in the standard collision-less CDM model. 1 To show this tendency we perform numerical simulations of formation of toy galactic systems using Gadget2 for the CDM and PyUltraLight for the FDM, respectively. Since it is challenging to simulate a full-fledged cosmological structure formation with the FDM due to the wave nature of the model, we restrict ourselves to studying the differences in structure formation at galactic scales with controlled initial conditions for seed satellite galaxies in the two models. (There are cosmological scale hydrodynamical simulations with FDM and baryon matter for structure formation [22].) In section 2 we describe our numerical simulation setup for formation of galactic systems in the FDM and the CDM models. We shall apply the numerical results to find the motion of satellite galaxies using the so called Mulguisin halo finder. Section 3 contains the results of our numerical study. We analyze and compare resulting satellite distributions of the two DM models. In section 4 we compare our results with the observed astronomical data. The last section is devoted to concluding remarks.

Simulation setup
In this section, we would like to describe our simulation setup for the above two DM models. As was mentioned already, we are not trying to do a full-fledged cosmological simulation, which is certainly beyond scope of our current study. Rather, we setup specific initial configurations for the late-time formation of galactic satellite systems, focusing on the purely gravitational dynamics in a fixed flat background. Our primary aim is to see the resulting differences between the FDM and the CDM dynamics.

FDM subhalo seeds
Our FDM dynamics are governed by the Schrödinger-Poisson equations where m (= 10 −22 eV/c 2 ) is the mass of ultralight scalar particles, G for the Newton constant, and M tot denotes the total mass of the galactic system of interest. The wave function ψ(x, t) is a complex function of spacetime normalized as d 3 x |ψ| 2 = 1, the corresponding FDM JCAP12(2022)033 mass density then given by ρ = M tot |ψ| 2 , and V for the gravitational potential produced by the density ρ. The m value is the fiducial mass suggested to solve the small scale problems such as the missing satellite problem. If we use a much larger value for m we expect the wave nature of FDM becomes less relevant and the difference of FDM from CDM becomes less prominent.
As an initial subhalo configuration, we shall use the ground state solution (see e.g. [16]) of the Schrödinger-Poisson system with a mass M s , which will be called a 'soliton' below for simplicity [23]. Let us denote this soliton wave function centered at x = 0 by h(x; M ), whose explicit functional form is known numerically. 2 It is also known that the half mass radius r 1/2 of the soliton is given by r 1/2 3.3541 × 10 8 kpc M /M s [16]. This soliton configuration is spherically symmetric, stable under a small perturbation, and may work as a seed for latetime satellite galaxies in our simulation. Although we are not explicitly presenting it in this note, a moving soliton solution with a constant velocity is rather well known in terms of the above wave function h(x; M s ) [25]. Using this, one may obtain an initial wave function [25] describing an initial condition for a soliton with an initial velocity v 0 centered at an initial position x 0 . The phase ϕ 0 (∈ [0, 2π)) will be assigned randomly for each soliton below. In this note, we shall take the initial subhalo mass to be M s 3.142 × 10 8 M leading to r 1/2 1.068 kpc. This choice, which is slightly larger than the typical satellite mass of nearby galaxies, was made because, through interactions with the rest, the seeds generically lose their mass in the formation of their satellite system.
For our FDM simulation, we shall use the python pseudo-spectral solver, PyUltraLight, developed in [25], which is publicly available. In this package, the whole system is placed in a box of size L with a periodic boundary condition x i ∼ x i + L. We shall take the box size to be L 204.0 kpc and the number of the grid points (N g ) in each direction is set to 900. Hence our spatial resolution becomes ∆x = L/N g 0.2267 kpc. We choose the simulation time-step, ∆t 0.7861 Myr, which is slightly smaller than the default value of the python code, ∆t d 0.8342 Myr. See [20] for a detailed explanation of a similar choice. In PyUltraLight, the full initial wave function is prepared by a simple superposition of the initial soliton wave functions describing initial subhalo seeds with given velocities. To avoid any significant overlap contributions, we separate initial positions such that one may ignore their overlap contributions. With this initial configuration, the Schrödinger-Poisson equations are solved with the above specified boundary condition, which produces time-series sets of wavefunction maps leading to an evaluation of relevant physical quantities such as density maps.

CDM subhalo seeds
We now turn to our simulation setup in the CDM side. Here we shall basically perform a gravitational N -body simulation using Gadget2 [26] which is also publicly available and we would like to imitate our initial FDM halo configurations as much as possible. Specifically, we set the mass of the initial subhalo to be M s as before and use the same sets of initial positions and velocities as the FDM counterparts. However, a stable solitonic configuration is not available in the CDM side. We instead prepare an initial CDM subhalo based on the 2 We take this to be real without introducing any extra phase at this stage. In [24], an approximate fit function of this soliton configuration was introduced as h fit = 0.51411/r 3/2 1/2 /(1 + 0.19191r 2 /r 2 1/2 ) 4 where r 1/2 is the half mass radius. This form will not be used in our simulation below.

JCAP12(2022)033
NFW profile [2] whose mass density function is given by where R s is the length scale of the profile. In particular we generate N p (= 2000) CDM particles with a mass M s /2000, whose initial positions are randomly allocated following the NFW profile with a cut-off scale set by 3R s . The two parameters, ρ 0 and R s are fixed by the conditions where r 1/2 will be further adjusted in the following way. We set these 2000 particles to move only in angular directions, selected randomly, and their speed to (1 − h ) GM (r)/r with a judicious choice of h = 1/30 where M (r) denotes the total mass of particles within a radius r. Now performing a separate N -body simulation, we let the subhalo system evolve for 1.0 Gyr such that the subhalo particle system is relaxed into a dynamical equilibrium. We adjust the initial parameter r 1/2 such that the resulting half mass radius agrees with that of the FDM soliton (r 1/2 ). In fact we found that r 1/2 ≈ r 1/2 with the above choice of h . We then use the resulting particle positions and velocities as our actual initial subhalo configuration. This relaxed configuration slightly differs from the original NFW profile we began with. However we view that this does not matter due to the following two reasons; first the NFW profile would not be that precise configuration that follows from the CDM simulations [27,28]. Second we tried other models such as the Plummer and the above soliton and found that the resulting differences are rather negligible. Gadget2 is a cosmological N -body simulation code based on the so-called TreeSPH code, by which one may compute the nonlinear regime of gravitational dynamics and hydrodynamics. In the simulation, with N particles, we first assign their initial positions and velocities, and evaluate their gravitational time evolution by computing the gravitational force on each particle using a TreePM algorithm (see [26] for the details).
In our CDM simulation, we again put the N -body system inside a box with a size L box together with the spatially periodic boundary condition. In the simulation with Gadget2, we also set the particle mesh grid (PMG) scale to L box /128 1.594 kpc. This choice may offer a computational efficiency without much loss of the accuracy in the calculation of forces based on the purely Tree algorithm. We set the softening length (SL) scale to value acc 2r 1/2 / N p . This requirement (SL > acc ) is needed to prevent any strong discreteness effects where acc plays a role of setting the upper limit of acceleration in two-body encounters by the mean-field value of a subhalo [29]. We have also tested more smaller scale SL = 2b 4r 1/2 /N p , which ensures preventing large-angle deflections in two-body encounters [29]. We have tried various such choices and found that there are no significant differences in their performances.
We make our resulting datasets of CDM density fields with the same resolution as those of the FDM side for an unbiased comparison of the two sides. In this analysis, we use the Pynbody package [30] which is an open source code especially suitable for astrophysical Nbody and smoothed-particle hydrodynamics (SPH) problems. We down-size the datasets of JCAP12(2022)033 three dimensional density fields from original 900 cells along each side of the simulation box to 450 cells, which is mainly for our computational convenience while keeping the accuracy needed for an identification of the satellite galaxies. These satellites, in the present case, as small clusters of DM distribution will be identified with Mulguisin halo finder, whose details are described below.

Initial halo configuration
In order to prepare our initial configuration for seed subhalos, we use the Plummer model [31] whose density profile is described by where r 0 denotes the half light radius. This model is chosen to simply provide enough randomness for our initial halo configuration. We shall set r 0 10.423 kpc, which seems an appropriate choice compared to the size of our simulation box. Including the velocity dependence, the distribution function becomes Our initial halo configuration will be consisting of 200 seed subhalos. In each set of this initial configuration, we begin by randomly generating subhalo positions and velocities following the above distribution function. In order to prepare an out-of-equilibrium configuration, we shall then multiply an overall factor q 1 (< 1) to the above generated velocity v P for each seed subhalo. For the smaller q 1 (and also for q 2 introduced right below), the initial halo configuration becomes the more tightly bound gravitationally and, of course, the more subhalos are then infalling to the central region. Now in order to provide an overall angular momentum for our initial halo configuration, we proceed as follows. With the above prepared initial set for total 200 subhalos, we randomly select Z subhalos and add an angular velocity q 2 v + = q 2 v rot (r) sin θφ (q 2 ≤ 1) to each of them where v rot (r) is the rotation velocity GMtot(r) r with M tot (r) denoting the sum of masses within the radius r from the center of the box. To the remaining (200 − Z) subhalos, we add for the Z or the remaining (200 − Z) subhalos, respectively. Finally we translate, rotate and boost the above halo system such that it has a vanishing center of mass position and velocity involving only a z component of total angular momentum. The parameters q 2 and Z are introduced to obtain an appropriate spin of our halo. Of course, the random choice of the Z subhalos may provide an extra randomness for our initial subhalo configuration. In order to consider a well-bound system, one has to restrict the parameter regime to q 2 1 + q 2 2 1 where q 2 1 + q 2 2 is roughly characterizing the kinetic energy contribution of the initial halo system. In measuring the spin of our halo system, we shall adopt the dimensionless spin parameter [32] where J tot denotes the magnitude of total angular momentum and E tot the total energy of the system. From the CDM-based Millennium simulation [33], the average value of the spin parameter for its field halos is estimated as λ avg = 0.0422 [34].
In this note, we shall consider five cases of q 1 = 0.25, 0.3, 0.35, 0.4, 0.45. The parameter q 2 will be fixed to be q 2 = 0.75. Then, for each choice of q 1 and q 2 , we shall adjust Z to match the initial spin parameter of our initial halo configuration with the above λ avg value approximately within 5% errors. In this manner, one finds Z = 90, 90, 90, 91, 91 for q 1 = 0.25, 0.3, 0.35, 0.4, 0.45 respectively. Finally we repeat our simulation 30 times for each choice of (q 1 , q 2 , Z) where initial subhalo positions and part of velocities ( v P ) are fully randomly generated for each simulation. Below we shall refer the resulting initial condition specified by (q 1 , q 2 , Z) to "initial condition q 1 " since, with q 1 specified, q 2 and Z are fully fixed in the above choices of initial conditions. In our CDM simulations below, these choices of initial parameters will lead to a typical CDM flattening of satellite distributions derived from full cosmological simulations under ΛCDM (see section 4 for the details). One may try larger values of q 1 such as q 1 = 0.8 or even larger. The resulting initial configurations are no longer gravitationally bound; with q 1 = 0.8, one indeed finds that part of subhalos go out of the simulation box and also that the late-time relaxation of the system is hardly achieved due to the relatively large kinetic energy. Therefore such choices will be inappropriate as our proper initial setup.
We set our total simulation time to be 1887.0 Myr, which corresponds to a total of 2400 time steps (in ∆t). However, we find that some of the subhalos after initial collisions may go out of the box through one side and reappear through its opposite side, which is basically due to the periodic boundary condition imposed in this note. To avoid such instances with enough margin, we shall take the validity limit of our simulation time to be 1226.4 Myr.

Halo reconstruction
Given the initial halo configuration, the simulations using PyUltraLight for FDM and Gad-get2 for CDM are performed. For each time step, we save both FDM and CDM density maps for further halo reconstruction and data analysis. Galaxies and halos are the gravitationally bound objects and they appear as clusters in the density map data. A cluster can be defined as a large collection of points or cells in a small area or a localized volume. Cluster finding has been a popular topic in computational science. Many different clustering algorithms are available which have been developed and maintained by many different groups. As each clustering algorithm has its own strength and weakness, one should carefully choose the right algorithm for their own research purpose.
In this study, we have adopted the Mulguisin clustering algorithm, which has been used in particle physics for finding jet structures from the collection of particle tracks [35]. The algorithm was modified to find clusters out of input cells in 3D space. It first finds the most massive cell in the input data and name it a seed. Then it looks for the next massive cell which will be used as a test cell. If the test cell is close enough to the seed, we then attach the cell to the seed and these two cells become a group. A group can therefore be regarded as a list of cells. If the distance between the test cell and the seed is longer than a certain cut length, then it becomes a new seed. In this study we use 1.4 kpc as out distance cut length. Now the same procedure is repeated for the next massive cell. The cell will be attached to a group when the minimum distance between the test cell and the group is shorter than the cut length. For this we calculate all distances between the test cell and the cells in the group and choose the minimum value. Once again, if the test cell does not join to any -6 -JCAP12(2022)033 existing groups then it becomes a new seed. The algorithm keeps doing this simple process until there are no test cells left. At this stage, cells are all converted to groups. Some groups are big, i.e. have many cells, and some are not. There are also single-cell groups that are isolated from other groups.
At this stage a group is called a cluster when its mass is larger than a certain mass cut. Clusters are then further tested to see whether they form a super-cluster. A super-cluster can therefore be regarded as a pack of clusters that are attached with each other. The minimum distance between two clusters is calculated from the whole combinations of cells from each cluster and it is used to decide whether the two clusters are to be merged or not. We use the same distance cut of 1.4 kpc for the merging process.
All these procedures are implemented in our Mulguisin halo finder software. Figure 2 shows the results from cluster finding by the Mulguisin algorithm. The clusters that were found by the halo finder are marked with black circles.

Numerical results
In this section we shall present the outcome of our numerical study. In figure 1 we depict the time evolution of DM density during the formation of sample galactic halos (either in the CDM or the FDM model) as a time series of projected DM density maps in the unit of 10 10 M /kpc 2 . To make these projected maps, the mass density fields of each halo system are integrated over the range [−L box /2, L box /2] along y-direction (the left two columns) or z-direction (the right two columns) of the simulation box. Also for these maps, we used the samples from the simulations for the q 1 = 0.35 initial condition. Note that both the CDM and the FDM halo density maps shown in the figure are from simulations with the same initial positions and velocities of subhalos.
The CDM and the FDM halo system look quite similar at the early stages of the evolution process. But after the collapse of initial subhalo systems as illustrated in the second row of figure 1, their difference becomes clear as will be further described below: in the CDM side, there appear far more surviving subhalos, especially around the central region, than the FDM counterpart, and so is more substructure on the smallest scales. On the length scales of the whole parent halo, however, the tentacle-like structure found at late times is more prominent in the FDM than the CDM side.
Among clusters found by the Mulguisin algorithm, we shall identify objects with mass greater than 2 × 10 7 M , 3 as satellites, which were all marked with small black circles on the projected maps in figure 2. We find that the CDM system has roughly eight times more satellites than the FDM side. While there are many satellites in the central region of the CDM halos, the FDM halos involve almost none in that region. One may also see that the FDM satellite distributions appear more flattened than the CDM counterparts, even before analyzing any details. One can also see the outgoing waves from the gravitational cooling in the FDM side.
To be more concrete, we shall focus on the following three choices of initial conditions, q 1 = 0.25, 0.35, 0.45, particularly for figures 3-6 and their discussions.
First let us introduce the ratio of the semi-minor to the semi-major axis (denoted as c/a) of an ellipsoid derived from a satellite distribution [36][37][38] in the following manner. Its JCAP12(2022)033 principal axes are defined using the mass tensor of satellites where x (k) refers to the position of the k-th satellite and N s denotes the total number of satellites in each halo system. Note here that the origin of the position is defined by the center of mass in each halo system where its density peak is also located approximately. In fact, with our choice of initial halo configuration, this center of mass was set to agree with the center of the simulation box, and shall be simply called the center below. The eigenvalues of the tensor (ordered by λ 1 ≥ λ 2 ≥ λ 3 ) define the corresponding principal semi-axes respectively by a = √ λ 1 , b = √ λ 2 , and c = √ λ 3 . In figure 3, we depict the time evolution of c/a , i.e. the average of c/a for the 30 simulation samples prescribed in section 2.3, from the initial time t = 0 to 1226.4 Myr. In each plot, the black dash-dotted line is for the c/a of satellites system in the CDM halo model, the red solid for that of the FDM, the green dashed for their difference, and the blue dotted for the ratio of the latter over the former. The error bars in each plot represent the estimation for the standard deviation of the 30 simulation samples. It is clear that the FDM satellite profiles lead to a smaller c/a than their CDM counterparts. In addition, for both the FDM and the CDM, one finds that c/a gets larger as q 1 becomes larger. This is because, as q 1 becomes larger, the average strength of each random velocity v P in our initial halo configuration increases.
The r.m.s. height δ is the r.m.s. of vertical distances of satellites from their best-fitplane [38] that is determined by minimizing with respect ton wheren is a unit normal vector to the would-be best-fit-plane. In other words the best-fit-plane is obtained with all satellites in the whole simulation box. In fact, -9 -JCAP12(2022)033  the corresponding unit normal vectorn determined in this manner almost agrees with theê 3 eigenvector given above. One further introduces the r.m.s. height δ(< ρ) calculated only with the satellites within the horizontal distance ρ from the center, in which the best-fit-plane and unit normal vectorn are in advance fixed by the minimization procedure of δ in (3.2). In figure 4, we depict the r.m.s. height δ(< ρ) of satellites as a function of the horizontal distance ρ from the center at t = 1179.2 Myr, where we have taken ρ ≥ 10 kpc because we would like to avoid the centroid residing in the region within 10 kpc. We may confirm that the FDM systems have thinner satellite planes than the CDM counterparts. In particular, one finds that, for 10 ρ 30 kpc, the FDM r.m.s. height grows while the CDM counterpart does not change appreciably for the same range. This happens because, in the CDM side, for the above range of ρ already many satellites are found with a significantly build-up value of the -10 - An alternative to the ratio c/a is the r-dependent cumulative distribution function (CDF) f 1 (| cos θ|; r), 4 of a satellite system where one includes only satellites within the radius r from the center and cos θ is given byê 3 ·r with r being the position of a satellite andê 3 denoting the unit eigenvector corresponding to the eigenvalue λ 3 . We shall use the Kolmogorov-Smirnov probability P ks by comparing CDF samples of the FDM/CDM satellite systems where each CDF denoted by f 30 (| cos θ|) is drawn from the angular positions of satellites in all 30 samples for each choice of initial condition. 5 With the CDF's as well as P ks , we may deduce the main underlying factors that make the FDM halo system have a more flattened satellite system than its CDM counterpart.
In figure 5 we illustrate the CDF f 30 (| cos θ|) for the three initial conditions of our main interest. For each initial condition, this CDF with a subscript 30 is made by merging all the data points in each set of the 30 samples as explained in section 2.3. The black dash-dotted line is for the CDM, the red solid for the FDM, and the blue dashed corresponds to the reference spherically symmetric distribution. These plots clearly show that, for both the FDM and CDM models, the CDF's are quite distinguished from the spherically symmetric distribution. They also imply that the FDM model leads to more anisotropic satellite systems than their CDM counterpart. Figure 6 shows the r-dependent Kolmogorov-Smirnov probability P ks (< r) where one includes those satellites within the radius r from the center in order to build f 30 (| cos θ|; r). From figure 4, one may easily guess that the Kolmogorov-Smirnov statistic D 12 may decrease as a function of r especially for the region of r near the centroid. However, the number of satellites increases as r grows. Indeed the latter growth 4 Slightly abusing notation, we shall denote the full CDF including the whole satellites simply by f1(| cos θ|). 5 In general this quantity P ks represents the probability that any two sets of samples were drawn from the same cumulative probability distribution, where the probability is determined by the so-called Kolmogorov CDF whose argument is given by with Ni (i = 1, 2) denoting the size of the i-th sample and D12 the Kolmogorov-Smirnov statistic characterizing the maximal difference of two CDF's.  As mentioned in Introduction, all the satellites systems in our FDM halo models result in more flattened planes than those in the CDM counterparts with respect to all measures in this work. This trend may also be seen from the Kolmogorov-Smirnov probability P ks given in the table 1 for various initial conditions. Roughly speaking, the more similar the two distributions come to be, the larger is the corresponding Kolmogorov-Smirnov probability P ks . As q 1 becomes smaller, the initial collapse gets stronger, more subhalos are colliding into the central region, and the resulting satellite system becomes more flattened. Notice that these collisions mostly occur at an early stage of the evolution processes. The relatively -12 -

JCAP12(2022)033
stronger flattening in the FDM side is closely tied to the gravitational cooling effect, which is especially effective around the central region of the halos. It should be commented that P ks should be correlated with the parameter q 1 but, in table 1, it does not increase monotonically nor does it show any clear trend with q 1 . This is simply because P ks concerns about the difference between the two CDF's and the flattening trends of FDM and CDM may differ from each other.
We now turn to the corotation ratio η = N + s Ns where N ± s denotes the number of satellites involving a positive/negativeφ-component respectively. With the q 1 = 0.35 initial condition, one begins with η 0 = 0.4922 ± 0.0197 both for the FDM and the CDM halo. At t = 1179.2 Myr, we find the FDM corotation ratio η FDM = 0.9240 ± 0.0916, which is bigger than the CDM value η CDM = 0.5888 ± 0.0396. The resulting FDM ratio appears to be more consistent with observational data, but this requires further studies.
Finally let us comment upon satellite orbital poles, and in particular their standard deviation [39] defined by wheren avg is the normalized mean vector of orbital poles given byn avg = n avg / n avg with n avg = (1/N s ) Ns k=1n (k) . Note that this definition is somewhat different from the usual one in literature where the value from 7 or 8 best aligned angular momenta of the most massive satellites is taken.
We find a narrower distribution of FDM satellite orbital poles with

Comparison with observational data
For a thorough understanding of the galaxy formation problem, one has to perform a fullfledged cosmological simulation from the very early universe as was mentioned before. Since our results in this note are based on the simplified toy galaxy models of typical size O(10 2 ) kpc while adopting specific initial conditions, it is hard to compare our numerical results directly to observational data from real galaxies. Nonetheless, it is interesting to see that our toy galaxy models in the FDM side much resemble the observed at least qualitatively as discussed right below.
Our simulations in the FDM side typically show c/a = 0.21 ∼ 0.30 and sph 27 • (with the q 1 = 0.35 initial condition), which well agree with the above observed values for the Milky Way. Our simulations in the CDM side, on the contrary, exhibit typically c/a = 0.41 ∼ 0.53 and sph 49 • (with the q 1 = 0.35 initial condition), which are larger than the observed ones. These c/a values for the CDM are consistent with those of the planes of satellite galaxies in the EAGLE simulation with the CDM [36]. Thus one may say that the FDM 6 As previously mentioned, the definition of ∆ sph in the references is different from ours. On the other hand, the r.m.s. heights for the both DM models are around δ 10 kpc, which appears smaller than the observed. This discrepancy may be attributed to the relatively smaller size of the simulated centroid (carrying a typical mass M tot 5×10 10 M ) than that of the observed central galaxies possessing a heavier typical mass M tot = O(10 12 ) M . Assuming a similar average DM density for all galaxies, the spatial size of the heavier galaxies may become 2.7 [∼ (10 12 /5 × 10 10 ) 1/3 ] times larger than the simulated value. Therefore, we expect δ 27 kpc for the heavier galaxies, which is indeed comparable to the observed. To confirm this we need a simulation with a larger box, which is beyond the scope of this paper.
We interpret our numerical results in the following way. Energy loss with a conserved angular momentum during a collapse usually leads to a formation of a highly flattened astronomical structure. In both dark matter models, during the collapse satellites gravitationally exchange their energy and momentum with others and experience a kind of an averaging process. Those satellites getting a relatively large angular momentum during the process may escape the central region and easily survive, while satellites with a small angular momentum tend to fall into the central region.
In the FDM side satellites which get small angular momentum and fall into the central region can be absorbed into the centroid by losing their energy via the gravitational cooling effect. Apparently, in our simulations a large portion of the falling satellites in the FDM side are tidally disrupted and lose their cores and disperse completely. As a result surviving satellites in the FDM side have a tendency to form corotating planar structures orthogonal to the total angular momentum vector. Note that this effect is clearly seen in figure 1.
On the other hand, a large portion of similar falling satellites in the CDM side may survive the crossing through the centroid in the same situation, mainly due to the particle nature of the CDM model. Thus, we expect satellites in the CDM model without an efficient energy loss mechanism redistribute themselves more spherically than in the FDM model as seen in figure 2.
In our scenario we also expect that the width of a satellite plane in the FDM side could be roughly related to the size of the centroid, because the spatial size of the centroid provides a typical length scale for the satellite absorption process near the centroid which is responsible for the energy loss and the formation of the satellite planes as previously mentioned.
Another interesting observation is that, even with the same initial conditions, there are less FDM satellites surviving compared to the CDM side due to the same effects making satellite planes. This hints yet another route in the FDM model to solve the missing satellite problem of the CDM model, because the typical solution in the FDM model to the problem is usually attributed to the suppression of the initial matter spectrum [11] compared to the CDM model.
For more conclusive results we need more realistic simulations for larger galaxies, which is beyond the scope of this work. Since our simulations are DM-only simulations, they do not provide information of inclination angles of the satellite planes with respect to the galactic disks.

Conclusions
In this note, we have numerically shown that galactic satellite systems in the FDM model are more flattened and corotating than their CDM counterpart. This is basically due the -14 -

JCAP12(2022)033
gravitational cooling effect in the FDM side especially near the central part of galaxies. This energy loss mechanism is a unique feature of the FDM which is absent in other alternatives of the CDM such as warm dark matter. Our toy galaxies in the FDM model seem to reproduce the observational features such as axis ratios of the satellite planes.
Our work implies that the FDM model could be a way to solve the problem of the satellite-galaxy planes as well as other small scale tensions. Thus the galactic satellite planes could serve as a good test bed for dark matter model discrimination. For more conclusive results, cosmological simulations with baryons starting from more realistic initial conditions are certainly required, which we leave for our future studies.