© G. Aguilar-Argüello et al. 2022

1. Introduction

The current paradigm assumed for cosmic structure formation is the Λ cold dark matter (ΛCDM) model. This scenario encompasses collisionless cold dark matter (DM) and a cosmological constant (Planck Collaboration VI 2020). The dark halo structure of galaxies has been extensively studied in ΛCDM-related models, most frequently with rotation curves, disk dynamics, or galaxy scaling relations (Rodríguez-Puebla et al. 2016; Aquino-Ortíz et al. 2018). The dynamics of satellite galaxies are also critical tests (Wang et al. 2018), particularly for nearby galaxies (Cautun et al. 2020). Recent and upcoming surveys such as Gaia (Gaia Collaboration 2021), APOGEE (Majewski et al. 2017; Fernández-Trincado et al. 2020), DESI (Cooper 2021), WEAVE (Dalton et al. 2012), and LSST (Rich 2018) will produce an exquisite mapping of the Milky Way’s (MW) stellar halo. In addition, these surveys will open up an opportunity to track the subtle stellar response that is presumably paired with the DM halo response, constraining such processes as dynamical friction as well as the properties of halos, or even those of DM itself, as recently shown via observations by Conroy et al. (2021). Using high-resolution N-body simulations, Garavito-Camargo et al. (2021) investigated the halo response modes to satellite accretion following such pioneering studies as Weinberg (1989). In addition, Ogiya & Burkert (2016), Tamfal et al. (2021) used super high-resolution simulations to study the relative importance of global mode and local wake in dynamical friction.

The recent studies mentioned above (Ogiya & Burkert 2016; Garavito-Camargo et al. 2021; Tamfal et al. 2021) have revealed a rich response mode population in fully self-consistent calculations reaching particle numbers of around 10⁸ − 10⁹ based on large computational resources. However, exploring a large parameter space using N-body simulations, while possible, may prove challenging due to the large amount of computational resources required. Therefore, devising an alternative technique may be useful. In the field of numerical simulations, testing results using different techniques consistently contributes to their individual robustness. Additionally, a continuum medium approach to dynamical problems is always interesting thanks to shot-noise handling, which is critical when the main mechanism we aim to capture is low-amplitude overdensities, as touched upon in the discussion above.

In this work, we propose using a flexible alternative or complementary strategy based on a continuous medium description for matter density. Our code is based on recent Collisionless Boltzmann Equation (CBE) solver implementations aimed at studying cosmic neutrinos (Banerjee & Dalal 2016, hereafter BD16). In cosmological simulations with massive neutrinos, the shot-noise density fluctuations are comparable to the model’s actual initial inhomogeneities, triggering artificial power. Traditionally, this is alleviated thanks to the use of a large number of particles, but a continuous medium description such as BD16 is an interesting choice. The BD16 method solves the Boltzmann moment equations on a grid and uses tracer particles to estimate the higher-order velocity moments of the Boltzmann hierarchy (instead of the density). We illustrate the code’s performance by tracking isolated halo and satellite sinking simulations.

In Sect. 2, we describe the general aspects of the BPM code. Section 3 presents the simulations following the wake, the global density, and the kinematics response in sinking satellite experiments exploring the dependence of different host halo density profiles. Our results are presented in Sect. 4, including comparisons to N-body codes (Sect. 4.2). Finally, our discussion and conclusions are given in Sect. 5.

2. BPM code

Solving the CBE+Poisson equations directly allows us to address many important aspects of the dynamics of self-gravitating systems by maintaining the advantage of a continuous media description while avoiding particle shot-noise. However, such an approach is quite expensive due to the high dimensionality of the phase-space, however the world’s top supercomputers could make it feasible, as, for instance, illustrated in recent works by Yoshikawa et al. (2013, 2020). To reduce the cost compared to the direct approach, we propose, instead (following BD16) to use a moment hierarchy approach. Such a technique allows for a 3D description that is more manageable from the computational point of view because it reduces aspects of dimensionality, as compared to a full Vlasov one. We calculate the first two CBE moments or the continuity and Euler-like equations coupled with the Poisson equation on a fixed 3D Cartesian grid. The order of the moment expansion may be questioned, in particular, crossed terms are not included. However, in contrast to self-consistent field codes (e.g., Weinberg 1989), our approximation is not directly performed on density or gravitational potential. As a healthy measure, we present in Sect. 4.2 a quantitative comparison to full N-body codes such as a fast multipole method (FMM, Greengard & Rokhlin 1987; Dehnen 2000) and a particle mesh (PM, Hockney & Eastwood 1988) ones.

The set of equations that we solve are listed below:

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\frac{\partial (\rho v^i)}{\partial x^i}& =0,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial (\rho v^i)}{\partial t}+\frac{\partial (\rho v^i{v}^j)}{\partial x^j}& =-\rho \frac{\partial \mathrm{\Psi }}{\partial x^i}-\frac{\partial (\rho {\sigma }^{\mathit{ij}})}{\partial x^j},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}^2\mathrm{\Psi }& =4\pi G\rho .\hfill \end{array}$$(3)

This set is equivalent to what is included in the cosmological case presented by BD16. A critical point when using a moment expansion approach is the closure condition used to truncate the expansion. Following BD16, we close the CBE hierarchy by using test particles to estimate moments of the distribution function, such as the velocity dispersion tensor and the bulk velocity. We use, at all times, the continuous density from the CBE moments solution to solve the Poisson equation (Eq. (3)). Additionally, we evolved in time the test particles using the leap-frog integrator (Kick-Drift-Kick, KDK) and the gravitational potential calculated from the Poisson equation. Then, we use the trajectories of the test particles to calculate the velocity dispersion tensor and the bulk velocity to close the CBE hierarchy. It is important to say that such kinematics estimation is the principal noise source; all the other quantities are computed at the continuous medium regime.

To solve the coupled moment equations (Eqs. (1) and (2)), following the BD16 implementation, we first use the directional splitting technique (Toro 2012). We then use a 1D piecewise linear advection solver, which is non-linear in space. In addition, we use the Superbee flux delimiter (Harten 1983), which has minimum diffusivity. The CBE moments are evolved in time by a total variation diminishing (TVD) Runge–Kutta on the order of two (RK2) integrator, which is known to preserve the TVD properties (Trac & Pen 2003). Therefore, our code (to which we refer to as BPM) is stable to artificial oscillations since it is TVD. Based on the previous discussion, the number of operations in BPM is similar to that of the PM code, exchanging the CIC cost for the advection part, which makes BPM more expensive than a PM.

Our implementation adds the possibility of following two separated species, either both of the species followed by the CBE moments (dubbed as CBE+CBE) or one species followed by CBE and the other by PM¹ (hereafter CBE+PM). In the CBE+PM case, the species are coupled only through the Poisson equation. This capability may be useful in attempting to follow a stellar component in addition to the DM density; however, in this paper, we are only considering the DM component. The CBE+PM version is also useful because it allows for savings in computational cost, as compared to the CBE+CBE version.

3. Simulations

To study the response of the underlying DM halo to a satellite accretion, we used spherical halos without any loss of generality. Since the typical particle speed is small compared with the speed of light, we only solve the continuity-like equation (1st CBE moment). This has been illustrated in BD16 in reference to their isolated Plummer halo case and cosmological CDM experiments. For the purposes of analysis, it is convenient to simulate the host halo and satellite as different species. As mentioned previously, BPM has two versions for simulations with two species: CBE+CBE and CBE+PM. We performed tests using both versions of the code; in particular, for the CBE+PM version, we evolved the host halo using CBE and the satellite using PM. We found that both CBE+CBE and CBE+PM gave similar results regarding the satellite sinking history; however, the CBE+PM is less computationally expensive. For this reason, the results that we present are based on the CBE+PM version of the code.

3.1. Models and initial conditions

The host halo configuration is either a Plummer (Plummer 1911) density profile or a cuspy NFW (Navarro et al. 1997) model, with spherical symmetry in both cases. For the satellite, we consider only a Plummer density profile. The satellite and the halo have the same mass resolution. The initial conditions were generated by sampling the Plummer distribution function and the NFW model using the publicly available MakeHalo code (McMillan & Dehnen 2007). The mass ratio between the satellite and host is 10%. We placed the satellite at a distance of 90 kpc with an orbital circularity of 0.8. Table 1 presents the parameters of the models. We generated the initial conditions for the BPM code directly using the analytical density solution from the equations in McMillan & Dehnen (2007). It is important to state that we did not pretend to simulate a specific galaxy model. We used three distinct halos with V_max close to 220 km s⁻¹ but with different density profiles; a Plummer model, an NFW with a similar circular velocity out of 30 kpc (NFW1), and a 20% less massive NFW halo with a higher concentration but with the similar maximum circular velocity (NFW2), as seen in Fig. 1. The satellite is represented by a Plummer model with a core radius of 3.6 kpc. Before starting the merger simulation, we evolved all the halos in isolation for 2.0 dynamical times, t_dyn.

Fig. 1. Dark halos’ circular velocity curve. NFW1 (blue) has a similar mass as the Plummer (black) model and also the same circular velocity out of 30 kpc and NFW2 (red) is 20% less massive but more concentrated and with a similar Vmax.

3.2. Numerical experiments

We performed all calculations on a grid of 512³ cells and, unless explicitly stated, 256³ test particles. The time step was selected according to particles velocity dispersion or sound speed for continuous medium and stability analysis (BD16). We adopted 1 kpc as the spatial resolution for the BPM and PM codes, and 0.3 kpc for the Tree code. The first set of tests followed the stability of the isolated halo for each of the halo models with the BPM code. Figure 2 shows that BPM held the density profile of all the three halos for 4.0 t_dyn. Next, we performed the sinking satellite experiments. As expected, the satellite sinks due to dynamical friction. As expected, the cored Plummer profile shows the satellite core stalling (see Petts et al. 2015, and references therein), and the NFW1 and NFW2 models develop a core formation at the end of the run. More details of this process are presented in the following sections.

Fig. 2. Stability test for BPM. We present the density profile for each of the three halo models at different times up to 4.0 tdyn. The density profile inside 100 kpc barely changes, thus confirming that the BPM accurately calculates its evolution.

4. Results

4.1. Discussion of the results

It is critical to accurately find the host and satellite centers since a defect in such an estimation would ultimately trigger a false density response. We estimate such positions using the shrinking-sphere method (Power et al. 2003). The case for the Plummer halo was challenging. We tested our centroid estimation against the Rockstar halo finder (Behroozi et al. 2013), and the results of both estimations are compatible. Once the initial and the current snapshot share the same centroid, we proceed to calculate the overdensity maps. The BPM density is a direct code solution. The results are presented in Fig. 3, where we show a slice of the overdensity maps in the orbital plane with one cell thickness (1 kpc). From top to bottom: Plummer, NFW1, and NFW2 overdensity maps are presented for (left to right) t = 0.3, 1.0, and 2.0 t_dyn. At t = 0.3 t_dyn, the local wake behind the satellite is clearly seen in all models. It is also evident that there is an overdensity and an underdensity dipole (blue and green, respectively). The amplitude of such a dipole is larger for NFW2, followed by NFW1, and finally Plummer. For t = 1.0 t_dyn, the wake is still visible for all three models, but the relative amplitude regarding the global dipole is getting even for the NFW models. At t = 2.0 t_dyn, it is hard to see the wake, but the global mode is still clearly detected. Figure 4 shows a zoom in the overdensity map inside the Plummer core radius at t = 2.0 t_dyn; an over- and underdensity around the satellite is observed, showing that there are modes exchanging angular momentum besides the large scale one. A similar zoom in the NFW models did not show the same structure.

Fig. 3. Evolution of halo overdensity (ρt/ρ0 − 1) maps for the three host halos shown in Fig. 1. Purple colors are for negative overdensity values and red for positive. The upper row corresponds to the Plummer halo model, the middle row to NFW1 (same Vc as the Plummer model out of 30 kpc), and the bottom row to the smaller NFW2 mass model with higher concentration. Each column corresponds to 0.3 (left), 1.0 (center), and 2.0 (right) dynamical times (tdyn), respectively. The host center and the satellite positions are marked by a yellow cross and solid circle, respectively. The wake and the global response are evident in all panels. The Plummer cored model starts to deviate significantly once the core radius is reached by the satellite. There are small quantitative differences in the wake and dipole for each model. At 2.0 tdyn, the global anisotropy is smaller, but other modes are visible.

Fig. 4. Modes in over- and underdensity zoom map inside the core radius for the Plummer model. We can observe a small overdensity leading the satellite (cross) and also a trailing underdensity.

To further study the response of the host halo to the sinking satellite, we calculated the angular power spectrum (APS) of each halo. To that end, we first projected the halo’s density map (previously recentered on the initial halo centroid) onto a sphere with a 100 kpc radius. We then calculated the spherical harmonics coefficients (using the python package HEALPY, Górski et al. 2005; Zonca et al. 2019) c_ℓ, m of the projected density and summed over all m values, for a given m: C_ℓ := Σ_m|c_ℓ, m|². We refer to C_ℓ as a “mode”. Figure 5 shows the APS, normalized to the initial one, for each host halo and for the three different times. At t = 0.3 t_dyn (solid lines), the cuspy models are similar, however, the compact NFW2 one (red) seems less responsive, and the modes at ℓ = 3² and ℓ = 5 are similar in amplitude. In contrast, for the NFW1 halo (blue), the ℓ = 3 mode is dominant over the others. The Plummer model (black) seems the most responsive of the three halos based on the APS amplitude. The differences are more dramatic for the following times. At t = 1.0 (dashed lines) and 2.0 t_dyn (dotted lines), the APS amplitude of NFW2 grows considerably less than the other models, and the shape becomes smoother, suggesting that the compact and smaller NFW2 halo is less responsive than the NFW1 and the Plummer halos. The latter seems the most responsive because its APS amplitude grows considerably. The APS behavior is in agreement with the overdensity maps. There is a peak at ℓ = 3 and ℓ = 5 for all models, and particularly the Plummer one additionally captures peaks at ℓ = 7, ℓ = 9, and possibly more. At the end of the simulation (t = 2.0 t_dyn), both the NFW1 and the Plummer APS are heavily dominated by the ℓ = 3 mode, while the NFW2 one is less bumpy.

Fig. 5. Evolution of halo density angular modes for the three halo models calculated with the BPM code. The peak at ℓ = 3 (see footnote 2 2 for consistency with previous works) has the largest amplitude for Plummer (black) and NFW1 (blue) models at all times. However, the NFW2 (red) model shows an increasingly flatter spectrum. This is consistent with maps shown in Fig. 3, where the dipole has a comparable amplitude with the wake.

It is also interesting to see the evolution of the bulk velocity map in the orbital plane. Figure 6 shows the maps of the X (upper figure) and Y (lower figure) bulk velocity components for the Plummer, NFW1, and NFW2 halos (top to bottom, respectively), and from t = 0.3, 1.0, and 2.0 t_dyn (left to right, respectively). The left and central column panels show a negative or positive (blue or red) velocity gradient centered on the satellite illustrating the local wake formation. There is also a global velocity gradient (blue or red). At the rightmost panels, in agreement with the overdensity map, there is no clear signal of the local wake. We may ask if the velocity gradient indicates rotation. In Fig. 7, we show a zoom into the core region of the bulk velocity components maps (at t = 2.0 t_dyn), now including the Z component. It is evident that the rotation and all angular momentum transfer occur inside the orbital plane. All the density and velocity maps show dependence on the halo velocity profile. Although the results are encouraging, we need to test the robustness of numerical code dependence.

Fig. 6. Velocity maps. The map shows the X (upper) and Y (lower) bulk velocity components (positive red, negative blue). As in Fig. 3, the columns correspond to, from left to right, t = 0.3, 1.0, and 2.0 tdyn, respectively. The rows from top to bottom refer to the Plummer, NFW1, and NFW2 halos, respectively. Satellite and host centers are pointed by a yellow circle and a cross, respectively. The halo response behind the satellite is clearly seen in the first two columns. The map for NFW2 starts to deviate from the other two models at t = 1.0, suggesting that the momentum redistribution is different. The differences between the middle and right panels, for all the models, suggest that angular momentum transfers from the satellite to the background halo.

Fig. 7. Central kinematics at the end of the simulation for the three simulations models. Each row corresponds to a halo model, from top to bottom: Plummer, NFW1, and NFW2. The columns from left to right show the X, Y, and Z bulk velocity components, respectively. The first evident feature is that all the angular momentum exchange happens in the orbital plane (X − Y); therefore, the Z velocity component shows no pattern. Nearly cylindrical rotation is evident in the simulations’ final stage.

4.2. Comparison with N-body codes

Besides the convenient properties of BPM handling the shot-noise effects, it is also important to test the reliability of BPM in comparison with standard N-body codes. For this reason, we additionally ran some of our numerical experiments with our own PM code and also with a Tree code (Gyrfalcon, Dehnen 2000). Figure 8 presents the sinking satellite history for the three halo models and the three codes. In general, the three codes capture the same sinking satellite history with differences of around 10%. After 1.5 t_dyn, the sinking stalls for the Plummer model, as well as with all codes, as reported in previous studies (Petts et al. 2015). For the NFW models, the sinking rate also decreases; however, it is close to the resolution.

Fig. 8. Sinking satellite history for different models and codes. Solid, dotted, and dashed lines correspond to BPM, Tree, and PM codes, respectively. The three codes are consistent with each other up to t = 1.3 tdyn when small differences start to appear. The blue, red, and black colors correspond to NFW1, NFW2, and Plummer models simulations, respectively. The core stalling of the sinking satellite is evident in the Plummer model simulation.

In addition to testing the consistency of the sinking satellite history across codes, we also tracked the momentum evolution in experiments with the different codes. Figure 9 shows the linear (upper panel) and angular (lower panel) momentum during the 2.0 t_dyn. Indeed, the BPM (solid line) case has lower momentum conservation, likely because it only solves the continuity equation. The rate of change of momentum conservation is significant for BPM but does not grow systematically over time. Hence, the main dynamical effects are consistent with traditional full N-body codes, at least for the number of simulated dynamical times, as we can see in Fig. 8.

Fig. 9. Momentum conservation. Upper panel: linear momentum evolution for the three halos with two different codes: BPM (solid lines) and full N-body Tree code (dotted lines). Lower panel: angular momentum conservation for the three models with the corresponding codes.

Finally, we may consider whether BPM does indeed hold an advantage in tracing the overdensity signal in the simulations. Figure 10 shows the overdensity maps (from top to bottom) built based on: BPM, PM, and Tree codes (with two different cloud-in-cell, CIC, resolutions: 1.0 and 0.3 kpc, respectively). For the N-body simulations, we used the CIC interpolation to build the density and overdensity maps, while for BPM, we used the numerical solution of the continuity+Poisson equations. Density, in the BPM code, inherits some noise from the velocity estimation; however, it is less noisy than the N-body counterparts at a similar grid resolution. The fourth row shows the analysis of the Tree simulation using a CIC grid with equal spatial resolution as the simulation (i.e., 0.3 kpc, N_g = 1720). In this case, the satellite or the global mode is hardly seen. However, the structures start to appear when a lower resolution (1.0 kpc) CIC grid is used (third row). The above outlines that the analysis of traditional N-body simulations requires special treatment, for example, looking for convergence when smoothing the grid or using adaptive smoothing procedures. The PM case is similar to the Tree smoothed case. Finally, the BPM case directly shows the density and the overdensity fields out of the continuity equation with no smoothing required. All three codes show similar features (see first three rows); however, BPM makes it easy to identify small changes in the density field without any special analysis of the simulations.

Fig. 10. Overdensity maps noise effects with different codes and smoothing. The top row corresponds to BPM, where the overdensity is calculated directly with the solution of Continuity+Poisson equations. The second row shows the PM overdensity after applying CIC to the particle distribution. The two lowest rows correspond to the Tree code after applying CIC, using a grid resolution of 1.0 (5123 cells) and 0.3 (17203 cells) kpc, respectively. The highest resolution hardly shows the wake and global response. However, the smoothed case shows a host halo response close to the PM case, demonstrating a dependence on the grid smoothing.

Figure 11 shows the evolution of the APS estimated in simulations of the NFW1 halo using three different codes: BPM, PM, and Tree (solid, dashed, and dotted lines, respectively) with a spatial resolution of ∼1.0 kpc in the density field. The behavior of the APS is similar between codes: the ℓ = 3 mode has the largest amplitude, while the other odd modes decrease in amplitude as ℓ increases. This behavior remains similar for three snapshots. The BPM and PM codes capture the dominant ℓ = 3 mode with the same amplitude, and the Tree code captures the mode with a lower one. For the ℓ = 5 mode, the codes differ in amplitude, from highest to lowest: BPM, PM, and Tree, respectively. The relative amplitude between codes of the ℓ = 7 mode varies over time. For example, at t = 0.3, BPM has a higher amplitude, but PM and Tree have a lower similar one, while at t = 2.0, BPM and Tree have equal amplitude, and PM has a lower one.

Fig. 11. Evolution of density angular modes for the NFW1 halo with different codes. The peak at ℓ = 3 has the largest amplitude at any time. It has the same amplitude for BPM (solid line) and PM (dashed line), whereas it has a lower amplitude for the Tree code (dotted line). The amplitude of the peaks at ℓ = 5 and ℓ = 7 differs between codes.

We may wonder whether the PM code is already capable of capturing the response modes. For that reason, we performed a test with BPM and PM codes using N_p = 128³ particles and a grid of $N_{\text{g}}^3={512}^3$ cells (1 kpc resolution). Figure 12 presents the results for the Plummer halo at t = 1.0 t_dyn. The top and middle panels show the overdensity field (built similarly to Fig. 10) for the BPM and PM simulations, respectively. The bottom panel shows the corresponding normalized APS (green lines). We clearly see a significant amount of damping in the PM case (green dashed line); even the relative satellite-halo-center distance is shorter at the scale of 15%, suggesting that unresolved modes in the PM code slow down the satellite sinking in the BPM. Increasing the particle number to 256³ is not enough (magenta line in Fig. 12), and possibly 512³ may approach BPM; however, the memory cost increment would be a factor of 64 larger for PM in comparison with BPM, considering that the memory cost of both codes is scaled as $N_{\text{p}}+N_{\text{g}}^3$.

Fig. 12. Noise damped modes for the Plummer model. Upper panels: overdensity maps at the orbital plane using BPM (top) and PM (middle). The wake and the overall map are affected by noise in the PM test. Bottom panel shows the APS (green lines), where the ℓ = 3 peak is severely damped for the PM code (dashed line). Both codes use 5123 cells, 1283 particles, and a spatial resolution of 1 kpc. Additionally, we include a PM case using 2563 particles (magenta dashed line); the APS amplitude is greater than the previous test. The dependence on particle number suggests that there is a shot-noise effect in the PM APS.

Regarding the execution time, we made a comparison using a single core. Table 2 summarizes the setup of the simulations and their performance. It is relevant to state that we did not make an exhaustive comparison of code performance. We found that the PM code is ∼40% faster than the BPM one (for a similar simulation setup). Tree and BPM run at similar times, both codes using the same number of particles. However, BPM easily captures subtle fluctuations in the density field (see Fig. 10). Furthermore, from Figs. 11 and 12, we can see that the amplitude difference of the mode of order ℓ = 3 between BPM and PM is lower for the simulations with 256³ particles (Fig. 11) than in the case with 128³, suggesting that there is a shot-noise effect. If we increase the number of particles to reduce such an effect in the N-body simulations, then the memory cost (and possibly the execution time) would be higher than the BPM simulation (with 256³ particles).

5. Discussion and conclusions

The study of triggered modes in N-body simulations of sinking satellites, either to assess their signature in the stellar halo or to study their role in dynamical friction, has been a challenging task requiring large particle numbers and computational resources. Using a new hybrid non-linear CBE+Poisson solver, we successfully detected the overdensity and kinematics response in sinking satellite experiments. We explored the effect on the host halo response related to its density profile. We also discussed the properties of the modes during the core stalling phase.

The upper panel in Fig. 13 shows the evolution of the density profile for the host+satellite system assuming the three host halo models (initial correspond to dashed lines and final to solid ones). It clearly shows the formation of a core in the central NFW1 (blue) that turns close to the Plummer (black) density profile. The NFW2 model (red) develops a smaller core. The behavior is not surprising due to the mass ratio between the satellite and the enclosed mass (Goerdt et al. 2010). The new code, BPM, presented in this paper correctly captures such a process.

Fig. 13. Core formation. Upper panel: initial (dashed) and the final (solid) density profiles for the three halo models. Both NFW models clearly show that a core was developed at 20 (NFW1) and 10 (NFW2) kpc, respectively. Middle panel: initial circular velocity (dash-dotted line) and the initial 3D velocity dispersion (solid line with stars). Lowest panel: same kinematic profiles at the end of the simulation.

The evolution described above is consistent with the stalled sinking shown in Fig. 8. After 1.5 t_dyn, it is possible to observe that the satellite falling holds. Such core stalling phase has been discussed by several studies, such as Petts et al. (2015) and the references therein. However, there is no consensus on the explanation behind this effect. In our case, the stalling is observed at around 10 kpc for Plummer and around 2–1 kpc for NFW2 and NFW1.

Figure 4 shows the overdensity modes within the Plummer core. A similar figure in the NFW models shows no structure, possibly because the core radius is smaller and hardly resolved. We may speculate that the central modes in the Plummer core compete with the global mode to weaken dynamical friction. The low amplitude of the density modes in Fig. 4 and the high V_rms/V_c value in Fig. 13 (see also Petts et al. 2015) may contribute to the weakening of dynamical friction (Inoue 2011), in agreement with recent linear CBE perturbation study (Kaur & Stone 2021, appeared after our submission).

Our results suggest that the response of the host halo is sensitive to its density profile (see Fig. 3). The smaller and compact NFW2 halo shows the strongest amplitude dipole, which is also observed in the flatter APS (Fig. 5) compared to the two massive halos, NFW1 and Plummer. All of them show the local wake with a good signal to noise ratio. The evolution of the amplitude and modes of the APS during the simulation (Fig. 5) indicates that NFW2 is the least responsive halo, as compared to the Plummer one. This behavior can be understood by looking at the middle panel in Fig. 13, where we show the initial circular velocity (dash-dotted lines) and 3D rms velocity (solid line with stars). It is clear that the NFW2 model (red) has had a higher velocity dispersion since the beginning, possibly blurring the orbital resonances. The lower panel shows the same (V_c, V_rms) curves at the end of the simulation; the behavior of V_rms is correlated with the core stalling radius. The final stage where core stalling is present is even more different; however, a common feature is that the local wake is not observed. We summarize our conclusions as follows.

1. Our code, named BPM, is able to capture the response modes (local wake, global dipole) in the density and velocity fields with a comparable computational cost to traditional N-body simulations. Since BPM is based on resolving a moment hierarchy of the CBE+Poisson equations on a 3D Eulerian grid and using test particles kinematics to close the Boltzmann hierarchy, BPM is less sensitive to shot-noise than the standard N-body approach. Therefore, BPM facilitates the detection of subtle fluctuations in the density and velocity fields without special treatment (e.g., higher number of particles, post-processing analysis) that would otherwise be required to achieve the same level of accuracy in a standard N-body simulation.

2. Both the spatial distributions of over- and underdensities and the normalized angular power spectrum depend on the halo density profile, particularly the global dipole strength. The dependence of the overdensity map is more evident when the satellite approaches the core region, and it is related to the central V_rms/V_c ratio.

3. Most of the angular momentum transfer occurs in the orbital plane. At advanced stages of the simulation, we found a central cylindrical rotation (Fig. 7). In the case of the presence of a stellar component, this may be a sign of past massive accretion events.

4. During the satellite core stalling stage, there is no clear signature of the local wake; however, low-density contrast modes are detected (Fig. 4). An interesting follow-up project would consist of studying their contribution to the slower sinking rate.

5. Based on the leftmost panels of Fig. 3, we conclude that the recent detection of the galactic stellar halo response (Conroy et al. 2021) may contribute to additional constraints on the mass distribution of the dark matter halo in the Milky Way. However, the situation is degenerated both in mass and density profile.

6. The capability of BPM for avoiding shot-noise makes it a convenient choice for problems that require resolving low-amplitude modes, with a moderate computational cost in memory. Situations such as the modes discussed in galaxy dynamical problems or low surface brightness structures in galaxies are natural choices. In addition, BPM allows us to study the response modes even in the case of multiple concurrent perturbers.

7. The right-hand-side term (source) of the Euler-like equation might consider additional dark matter physics such as self-interaction or other effects (equation of state) that produce a possible change in the halo density response (sound speed). An exploration of such an interesting avenue in density and kinematics will be the subject of future studies.

Acknowledgments

GA acknowledges useful exchange with A. Banerjee. We thank S. Roca-Fábrega and H. Velázquez for valuable discussions. GA and AT thank support from a CONACyT PhD fellowship. GA, AT, and OV acknowledge support from DGAPA-UNAM grants IN112518, IG101222, and AG101620. The authors thank for the facilities of cluster computers:Atocatl/LAMOD-UNAM and Miztli/DGTIC-UNAM. LAMOD is a collaborative project of IA, ICN, and IQ Institutes at UNAM.

References

All Tables

All Figures

	Fig. 1. Dark halos’ circular velocity curve. NFW1 (blue) has a similar mass as the Plummer (black) model and also the same circular velocity out of 30 kpc and NFW2 (red) is 20% less massive but more concentrated and with a similar Vmax.
In the text

	Fig. 2. Stability test for BPM. We present the density profile for each of the three halo models at different times up to 4.0 tdyn. The density profile inside 100 kpc barely changes, thus confirming that the BPM accurately calculates its evolution.
In the text

Fig. 3. Evolution of halo overdensity (ρt/ρ0 − 1) maps for the three host halos shown in Fig. 1. Purple colors are for negative overdensity values and red for positive. The upper row corresponds to the Plummer halo model, the middle row to NFW1 (same Vc as the Plummer model out of 30 kpc), and the bottom row to the smaller NFW2 mass model with higher concentration. Each column corresponds to 0.3 (left), 1.0 (center), and 2.0 (right) dynamical times (tdyn), respectively. The host center and the satellite positions are marked by a yellow cross and solid circle, respectively. The wake and the global response are evident in all panels. The Plummer cored model starts to deviate significantly once the core radius is reached by the satellite. There are small quantitative differences in the wake and dipole for each model. At 2.0 tdyn, the global anisotropy is smaller, but other modes are visible.

In the text

	Fig. 4. Modes in over- and underdensity zoom map inside the core radius for the Plummer model. We can observe a small overdensity leading the satellite (cross) and also a trailing underdensity.
In the text

Fig. 5. Evolution of halo density angular modes for the three halo models calculated with the BPM code. The peak at ℓ = 3 (see footnote 2 2 for consistency with previous works) has the largest amplitude for Plummer (black) and NFW1 (blue) models at all times. However, the NFW2 (red) model shows an increasingly flatter spectrum. This is consistent with maps shown in Fig. 3, where the dipole has a comparable amplitude with the wake.

In the text

Fig. 6. Velocity maps. The map shows the X (upper) and Y (lower) bulk velocity components (positive red, negative blue). As in Fig. 3, the columns correspond to, from left to right, t = 0.3, 1.0, and 2.0 tdyn, respectively. The rows from top to bottom refer to the Plummer, NFW1, and NFW2 halos, respectively. Satellite and host centers are pointed by a yellow circle and a cross, respectively. The halo response behind the satellite is clearly seen in the first two columns. The map for NFW2 starts to deviate from the other two models at t = 1.0, suggesting that the momentum redistribution is different. The differences between the middle and right panels, for all the models, suggest that angular momentum transfers from the satellite to the background halo.

In the text

Fig. 7. Central kinematics at the end of the simulation for the three simulations models. Each row corresponds to a halo model, from top to bottom: Plummer, NFW1, and NFW2. The columns from left to right show the X, Y, and Z bulk velocity components, respectively. The first evident feature is that all the angular momentum exchange happens in the orbital plane (X − Y); therefore, the Z velocity component shows no pattern. Nearly cylindrical rotation is evident in the simulations’ final stage.

In the text

Fig. 8. Sinking satellite history for different models and codes. Solid, dotted, and dashed lines correspond to BPM, Tree, and PM codes, respectively. The three codes are consistent with each other up to t = 1.3 tdyn when small differences start to appear. The blue, red, and black colors correspond to NFW1, NFW2, and Plummer models simulations, respectively. The core stalling of the sinking satellite is evident in the Plummer model simulation.

In the text

	Fig. 9. Momentum conservation. Upper panel: linear momentum evolution for the three halos with two different codes: BPM (solid lines) and full N-body Tree code (dotted lines). Lower panel: angular momentum conservation for the three models with the corresponding codes.
In the text

Fig. 10. Overdensity maps noise effects with different codes and smoothing. The top row corresponds to BPM, where the overdensity is calculated directly with the solution of Continuity+Poisson equations. The second row shows the PM overdensity after applying CIC to the particle distribution. The two lowest rows correspond to the Tree code after applying CIC, using a grid resolution of 1.0 (5123 cells) and 0.3 (17203 cells) kpc, respectively. The highest resolution hardly shows the wake and global response. However, the smoothed case shows a host halo response close to the PM case, demonstrating a dependence on the grid smoothing.

In the text

	Fig. 11. Evolution of density angular modes for the NFW1 halo with different codes. The peak at ℓ = 3 has the largest amplitude at any time. It has the same amplitude for BPM (solid line) and PM (dashed line), whereas it has a lower amplitude for the Tree code (dotted line). The amplitude of the peaks at ℓ = 5 and ℓ = 7 differs between codes.
In the text

Fig. 12. Noise damped modes for the Plummer model. Upper panels: overdensity maps at the orbital plane using BPM (top) and PM (middle). The wake and the overall map are affected by noise in the PM test. Bottom panel shows the APS (green lines), where the ℓ = 3 peak is severely damped for the PM code (dashed line). Both codes use 5123 cells, 1283 particles, and a spatial resolution of 1 kpc. Additionally, we include a PM case using 2563 particles (magenta dashed line); the APS amplitude is greater than the previous test. The dependence on particle number suggests that there is a shot-noise effect in the PM APS.

In the text

Fig. 13. Core formation. Upper panel: initial (dashed) and the final (solid) density profiles for the three halo models. Both NFW models clearly show that a core was developed at 20 (NFW1) and 10 (NFW2) kpc, respectively. Middle panel: initial circular velocity (dash-dotted line) and the initial 3D velocity dispersion (solid line with stars). Lowest panel: same kinematic profiles at the end of the simulation.

In the text