© The Authors 2025

1 Introduction

Dark matter (DM) is essential for galaxy formation and evolution, making up about 84% of the matter of the Universe (Planck Collaboration 2020; Adame et al. 2025). It provides the gravitational framework for the appearance of galactic structure and helps explain discrepancies between observed rotation curves and the distribution of baryonic mass (Rubin & Ford 1970; Rubin et al. 1980), large-scale structures (Efstathiou et al. 1985), and colliding systems such as the Bullet Cluster (Markevitch et al. 2004). The DM halo and the stellar component of galaxies are dynamically coupled, meaning that their gravitational interactions influence each other. While the halo shapes the overall stellar structure of galaxies, the stellar component, in turn, induces gravitational responses within the halo, potentially leaving detectable signatures. Understanding these interactions improves our knowledge of DM, galactic dynamics, and the overall structure of the Universe.

One notable effect of this coupling between components is the gradual slowing down of the Galactic bar rotation, driven by the exchange of angular momentum between the stellar disc and the DM halo (e.g. Weinberg 1985; Hernquist & Weinberg 1992; Debattista & Sellwood 1998, 2000; Athanassoula 2003; Chiba & Schönrich 2022). This process is a result of the resonant nature of the system. Indeed, Tremaine & Weinberg (1984) found that the exchange of angular momentum in the halo is driven by particles that are nearly in resonance. The other classical example is the orbital coupling of satellite galaxies, governed by dynamical friction (Chandrasekhar 1943; Banik & van den Bosch 2022), whereby satellites lose orbital energy as they move through the DM halo, eventually spiralling inward (e.g., Tremaine et al. 1975; Weinberg 1986).

In recent years several studies have focussed on the impact of the interaction between the stellar and DM components on the precise structure of the halos. One of the most studied cases is the disturbance of the Milky Way (MW) halo caused by the infall of the LMC (e.g. Weinberg 1998; Gómez et al. 2015; Laporte et al. 2018a; Garavito-Camargo et al. 2019; Cunningham et al. 2020; Petersen & Peñarrubia 2021; Amarante et al. 2024). In addition, early theoretical works (Weinberg 1985; Hernquist & Weinberg 1989) predicted that stellar bars generate a quadrupole wake in the DM halo. Later N-body simulations confirmed the emergence of a DM overdensity in response to the bar potential (e.g. Weinberg & Katz 2002), though these studies primarily characterized the halo response as classical wakes. Petersen et al. (2016) demonstrated that DM particles can become trapped in the bar potential, along with stars, forming a persistent ‘shadow bar’ that contributes ∼ 10% of the total bar mass (see also Collier et al. 2019; Collier & Madigan 2021; Frosst et al. 2024; Marostica et al. 2024; Ash et al. 2024). Another example of interactions between components is the trapping of halo stars into bar-related resonances, proposed in Dillamore et al. (2023, 2024) as an explanation for the ‘chevrons’ observed in Gaia data (Belokurov et al. 2023).

The stellar spiral arms and the DM halo interact similarly. In general, the particles in the halo must respond to any perturbation in the potential. Therefore, a spiral perturbation in the potential will create a response in the halo. However, the nature and consequences of this response are surprisingly underexamined in the literature. Mark (1976) predicted a significant amplification of the stellar spiral density waves as a result of the loss of angular momentum to the halo. Later, Fuchs (2004) employed the shearing sheet model (Goldreich & Lynden-Bell 1965; Julian & Toomre 1966) to investigate spiral amplification in the presence of a live DM halo. He predicted that the maximum growth factor would increase significantly due to the coupling with the halo and pointed out that the transfer of momentum is exclusively mediated by halo particles whose orbits are in resonance with the spiral waves. This response can be described as a broader dynamical process: a ‘forced response’, whereby periodic perturbations (e.g. a bar, a satellite, and the spiral arms) systematically distort orbits, with the strongest effects occurring near resonances. In the same set-up, he also studied the response of the halo to a density wave in the disc, and reported the development of a wake. Subsequent studies found that halos with anisotropic velocity distribution increase the maximum growth factor of the spirals considerably (Fuchs & Athanassoula 2005). More recently, Sellwood (2021) tested these predictions using a complete N-body simulation and reported that a live DM halo had no significant effect on the growth rate of the spiral modes, although he showed hints of an increase in the final stellar spiral amplitude of about 20%.

In low-surface-brightness galaxies (LSBs), the disc-halo coupling plays a key role in the evolution of spiral arms, due to their low stellar densities and DM-dominated dynamics. Theoretical studies have predicted that halos with a small local vertical scale height near the disc or halos with a significant angular momentum can strongly amplify disc density waves through resonant interactions. Chiueh & Tseng (2000) proposed that such coupling mechanisms could not only excite spiral-arm formation but also influence their morphology, with tighter pitch angles associated with smaller halo scale heights. Recent simulations by Narayanan et al. (2024) confirm the importance of the geometry of the halo in driving spiral patterns in LSBs, showing that oblate halos can sustain long-lived global spirals.

While earlier studies have concentrated on the feedback effects of coupling on the stellar component, our work shifts focus to the morphology of the DM halo itself. In this work, we find clear signs of a spiral-arm-shaped structure in the DM halo in several MW-like simulations with cold DM. We start by characterizing the overdensities created in the DM halo by the stellar spiral arms in an isolated pure N-body simulation. We reproduce the response of the DM halo to the gravitational perturbation induced by the stellar spiral arms using a simple test particle simulation. This demonstrates that the relation between the stellar spiral arms and the DM substructure observed in the pure N-body simulation is dominated by a forced response to the modification of the potential. Finally, we observe similar interaction patterns in pure N-body simulations with a massive perturber, isolated hydrodynamical simulations, and cosmological simulations, proving the ubiquitous presence of DM spiral arms in the state-of-the-art simulations of MW-like systems.

In this work, we show that in any MW-like system with cold DM, fairly strong spiral arms will have an increase of about 10% in their total mass, due to the forced response of the DM halo. The existence of DM spirals implies that the local DM density varies non-axisymmetrically across the galactic disc and, in particular, close to the Sun. Thus, the torques in a MW-like spiral galaxy will have contributions from both the stellar and the DM components. In the future, the kinematic and density maps of the DM particles provided in our work might be used as test for the halo distribution of our Galaxy.

In Sect. 2, we briefly describe the simulations and methods used throughout the paper. In Sect. 3, we present our findings for the isolated pure N-body model. These results are then compared with a test particle model in Sect. 4, where we also perform some tests on the characteristics of the response. In Sect. 5, we test the presence of DM spiral arms in a set of different simulations. These results are discussed in Sect. 6. Finally, in Sect. 7 we provide a summary of our results and list our conclusions.

2 Simulations and methods

In this study, we analyse several simulations that present strong two-armed stellar spirals. The main analysis of this paper is based on a pure N-body model of a MW-like isolated galaxy, which we label Iso-NBody. This model is from our new suite of simulations run in the awarded project in the call from the Open Clouds for Research Environments (OCRE). The model is a slight modification of the L2 model in Laporte et al. (2018b), with three components (cold exponential stellar disc, stellar Hernquist (1990) bulge, and a Hernquist DM halo) but with no perturber. The resulting galaxy model displays a strong bar and prominent two-armed stellar spirals, with a nearly constant pattern speed of Ω_s ∼ 20 km s⁻¹ kpc⁻¹, reaching a density contrast of 20% at R ∼ 6 kpc (see Sect. 3.1).

To complement the analysis, we also examine other simulations, including a test particle model (TP), a tidal N-body model with a Sagittarius-like perturber (Sgr-NBody), a high-resolution hydrodynamical simulation (Hydro), and a cosmological simulation from the Auriga project (Au27). These models provide additional insights and robustness into the formation and evolution of spiral arms under different physical conditions and simulation set-ups. Further details of the simulations are provided in Appendix A.

We analysed the response of the DM halo to stellar spiral arms by focussing on the region |Z| < 4 kpc and R < 30 kpc. The potential impact of the vertical cuts on our results is found to be small (see Appendix B.1). Key measurables include: radially normalized density (δρ_r), which is the fractional deviation of the local density (integrating over the vertical dimension) with respect to the average (over all azimuths) at the same radius (Eq. (B.2)), which clearly emphasized the non-axisymmetric variations; mean radial and azimuthal velocities (⟨V_R⟩, ⟨V_ϕ⟩); and vertical breathing motion (V_{Z, br}), which quantifies vertical velocity asymmetry between the upper and lower halves of the disc (Eq. (B.3)). V_{Z, br} traces the compression and expansion caused by the quadrupoles in the potential (spiral arms and bar). We used a Fourier decomposition to detect the spiral structure, with the m = 2 mode tracking the phase (ϕ₂) and amplitude (Σ₂) of the spirals. The spiral pattern speed (Ω_s) was derived from a linear regression of the phase over a ±0.1 Gyr window. Faceon images, smoothed with Gaussian kernels, were generated to aid interpretation. For simulations with perturbers, we applied a centro-symmetrization (rotate by 180^∘ and average with the original) to remove the dipole signatures coming from the response to the perturber, and highlight the quadrupole signature of the spiral arms. Further details about the methods used are provided in Appendix B.

Fig. 1 Spiral arms in the Iso-NBody simulation, at t = 7 Gyr. Left: surface density plot of the stellar component. Middle: radially normalized density (δρr) of the stellar component. The contours represent the regions where δρr = 0 (in black), and δρr = ±0.2 (in white). Right: δρr of the DM particles around the stellar disc (|Z| < 4 kpc). The contours are the same as the middle panel. We observe a clear spiral arm pattern in the DM halo, strongly correlated with the stellar spiral arms.

3 Dark spirals in a pure N-body simulation

3.1 Density

We start by studying the Iso-NBody model. In the stellar density we clearly observe strong two-armed spirals reaching to R ∼ 15 kpc with a peak δρ_r ∼ 0.2 at R ∼ 6 kpc (left and central panels of Fig. 1). Studying the temporal series (ochre lines in Fig. 2) of the amplitude (Σ₂/Σ₀) and phase (ϕ₂), we observe that these spirals appear at t ∼ 4.5 Gyr and grow with time with an approximately exponential profile to Σ₂/Σ₀ ∼ 20% at the end of the simulation. The pattern speed (Ω_s) of the stellar spiral arms (central panel of Fig. 2) is ill-defined before the appearance of the spiral, and stabilizes at Ω_s ∼ 20 km s⁻¹ kpc⁻¹ once the spiral emerges.

Studying the radial profile of the stellar spiral arms at a given time (t = 7 Gyr, solid ochre lines in Fig. 3), we confirm that the strength of the spiral arms peaks at R ∼ 6 kpc, and observe that their pattern speed is almost constant with the radius, going from the bar region up to the outer parts of the disc (solid ochre line in the third panel of Fig. 3). For discussion about the pattern speed for different models, we refer to previous results (e.g. Grand et al. 2012; Roca-Fàbrega et al. 2013; Antoja et al. 2022).

We now focus on the effect of the spirals on the DM halo (right panel in Fig. 1). We observe a DM spiral arm-shaped overdensity, nearly coincident with the position of the stellar spiral arms. The overdensity of the spiral arms in DM is δρ_r ∼ 0.02, one order of magnitude smaller than that of the stellar ones. In fact, the density contrast in DM spirals is roughly one order of magnitude lower than the density contrast of stellar spirals for all the models that we have studied (see Sect. 5).

We see a strong correlation between the two spirals in terms of temporal evolution (dashed grey line in Fig. 2, after t ∼ 4.5 Gyr). The relative amplitude of the DM spirals (top panel) evolves linearly with the stellar one, and is about 10% of the stellar relative amplitude at all times after 4.5 Gyr (as is indicated by the auxiliary vertical axis in the plot). The pattern speed of the spiral arms is also almost identical in the two components (stellar and DM) once the signal is strong enough (t ∼ 4.5 Gyr). Finally, in the bottom panel of Fig. 2 we show the phase lag between the DM spiral arms and the stellar spiral (long-dashed red line). We observe that initially there is no correlation: the dispersion is σ ∼ 45^∘, which is the expected dispersion for a uniform distribution. Once the stellar and DM spirals are formed, the DM spiral is consistently trailing the stellar spiral by ∼ 10^∘. This phase lag is a natural consequence of the delayed response of the DM halo to the evolving gravitational potential of the spiral arms, comparable to the mechanism at work in barred galaxies (e.g. Weinberg 1985; Athanassoula 2003).

Finally, we focus on the radial profile of the DM spirals at t = 7 Gyr (Fig. 3). In the top panel, we notice that the local density (Σ₀) of stars and DM is very similar in this model. The relative amplitude, Σ₂/Σ₀, of the DM spiral arms (second panel) varies with radius but remains consistently 10% of the stellar amplitude, as is shown by the secondary vertical axis on the right. We thus reveal a linear relationship between the stellar and DM amplitudes as a function of radius, matching the relationship observed in the temporal evolution (top panel of Fig. 2). This corresponds to adding about 10% to the total mass of the stellar spiral arm, since Σ₀ of both the stellar and DM components are of the same order. Finally, we confirm the consistent pattern speed of the stellar and DM spiral arms across all radii (third panel in Fig. 3), and the lagging of the DM spiral with respect to the stellar spiral at all radii where the signal is strong (bottom panel in Fig. 3). We see a small trend of the phase lag with radius that is a consequence of the local properties of the stellar spiral arms and DM distribution function (DF), which one could potentially predict using perturbation theory.

Fig. 2 Temporal evolution of the Iso-NBody simulation in the range R ∈ [5.5, 6.5] kpc. Top: relative amplitude of the Fourier mode (Σ2/Σ0), i.e. strength of the spiral arms for the stellar (solid ochre) and DM (dashed grey) components. We see that both spirals evolve in a correlated way, although the relative amplitude of the DM spirals is one order of magnitude smaller (notice the auxiliary axis, which refers to the amplitude of the DM spiral arms). Middle: pattern speed (Ωs) of the mode m = 2 structures in stars (solid ochre) and DM (dashed grey). Bottom: phase difference between the mode m = 2 of the stellar and DM structures (long-dashed red line). Shaded regions shows the dispersion of the phase lag. Once the spirals are strong enough (t ≳ 4.5 Gyr, vertical lines), we detect a constant pattern speed of Ωs ∼ 20 km s−1 kpc−1 for the stellar and DM spirals, and a constant phase lag of Δ ϕ2 ∼ 10∘, showing a strong correlation between the two spiral arms.

3.2 Velocity space

In a dynamical system, the appearance of an overdensity in the configuration space must leave a trace in the kinematics. In this section we study ⟨V_R⟩, ⟨V_ϕ⟩, and V_{Z, br} of the stellar and DM components in the Iso-NBody model.

We start by analysing the kinematic imprints in the stellar spirals. In the top row of Fig. 4, we show the kinematic maps of the stellar component, compared with the position of the bar and spiral arms (shown as contours in the figure). We observe three clear regimes: the bar, the inner disc, and the outer disc. In the third row of Fig. 4, we show the amplitude of the velocity quadrupoles at each radius (solid ochre line), and confirm that the in-plane velocity quadrupoles (⟨V_R⟩ and ⟨V_ϕ⟩) reach about 20 km s⁻¹ in the inner part of the disc, and 5 km s⁻¹ in the outer parts. As for the phase of the quadrupole (last row of Fig. 4), which we measure with respect to the phase of the stellar density structures, we observe a 90^∘ change in the phase of the quadrupole inside and outside the co-rotation¹ of the spirals (vertical red line). The kinematic signature of the spirals in our model is rather complex and depends on the radius and how the spirals connect with the bar at each time. However, since this is not the main focus of this article, we refer to the ongoing discussion on the topic (Siebert et al. 2012; Grand et al. 2015; Monari et al. 2016a,b; Antoja et al. 2016, 2022; Eilers et al. 2020).

In the DM particles (second row in Fig. 4), we also observe quadrupoles in the kinematics. In ⟨V_R⟩ and ⟨V_ϕ⟩, we observe that along all the stellar spiral overdensity (solid contours) we have positive values, reaching amplitudes of ∼ 2 km s⁻¹, much smaller than the stellar kinematic amplitudes. In the third row of Fig. 4 (dashed grey lines), we show the amplitudes of the kinematic signatures. For ⟨V_R⟩, the amplitude of DM behaves qualitatively similar to the stellar one, with a peak at R ∼ 4 kpc, and lower values at outer radii. However, the radial profile of ⟨V_ϕ⟩ shows a different behaviour when comparing the stellar and the DM components, with the latter peaking at R ∼ 8 kpc. In the last row of Fig. 4, we show the phase of the velocity quadrupoles in the DM halo (dashed grey lines). The distinct in-plane kinematic signatures between the stellar and DM spiral arms reflect their differing intrinsic dynamics: the stellar component, governed by rotationally dominated orbits, responds to the spiral potential through self-gravity and angular momentum exchange, while the DM halo, with its near-isotropic distribution, exhibits a kinematic response shaped primarily by the forced response rather than by self-gravity.

In the vertical response of the stellar component (V_{Z, br}, top right panel of Fig. 4), we also observe a clear quadrupole, but with a small amplitude of ∼ 2 km s⁻¹ in the inner parts of the disc, and <1 km s⁻¹ in the outer regions. It should be noted that V_{Z, br} is tracing the compression and expansion caused by the spiral potential. Outside co-rotation (R>10 kpc), the breathing motion is contracting (purple) in the leading part of the arm and expanding (orange) in the trailing. Inside co-rotation (R < 10 kpc), where the stars orbit faster than the spiral arms and ‘overtake’ them, there is a change in the phasing of the vertical force that leads to a change in the sign of the quadrupole (phase shift of 90^∘ in the bottom panels). These findings are compatible with the predictions for the vertical kinematics induced by grand design spirals (Debattista 2014; Faure et al. 2014).

The vertical response of the DM halo is shown in the third column of Fig. 4. In this case, the breathing motion contracts (purple) in the leading part of the arm and expands (orange) in the trailing part. This is similar to the stellar breathing motion only outside co-rotation. In a non-rotating halo, on average the particles encounter the spiral arm through the leading part at all radii, explaining the similarities between components of the vertical response outside co-rotation. Furthermore, in the third and fourth rows of the figure we observe that the amplitude (∼1 km s⁻¹) and phase of the DM mode are almost constant with radius. While the differences in the planar motions of the stars and DM particles hint that we are observing dynamically different components, the similarities in the vertical breathing motions are not unexpected as they naturally result from the response of a passing overdensity corresponding to the spiral arms.

Fig. 3 Radial profile of the Iso-NBody simulation at t = 7 Gyr. We show the Fourier amplitude and phase for the stellar (solid ochre) and DM (dashed grey) components. First: amplitude of the 0th Fourier mode (Σ0), in units of solar mass. We see that, close to the midplane, the surface density of stars and DM is comparable. Second: relative amplitude of the Fourier mode (Σ2/Σ0) at all radii. The relative profile of the spirals is similar, with the relative amplitude of the DM spirals being one order of magnitude smaller than the stellar one (notice the auxiliary axis). Third: pattern speed (Ωs) of the mode m = 2. Fourth: phase lag (Δϕ2, long-dashed red line) between the DM and stellar components.

Fig. 4 Kinematic analysis of the Iso-NBody simulation, at t = 7 Gyr. The contours represent the regions where δρr = 0 (in black) and δρr = ± 0.2 (in white). We show ⟨VR⟩ (left column), ⟨Vϕ⟩ (middle column), and VZ, br (right column). First row: kinematic maps of the stellar component. Second row: kinematic maps of the DM component. Third row: radial profile of the kinematic amplitude of the Fourier mode (A2) for the stellar (solid ochre) and DM (dashed grey) components. Fourth row: radial profile of the phase difference between each kinematic mode (ϕ2) and the stellar spiral in density (ϕ2(ρStars)). The co-rotation radii of the spirals (RC R ∼ 10 kpc) is shown as a vertical red strip.

Fig. 5 Response of the halo test particles to an analytical spiral potential (Cox & Gómez 2002), rotating at a constant pattern speed. From left to right, we show the radially normalized density (δρr), the mean radial and azimuthal velocities (⟨VR⟩,⟨Vϕ⟩), and the vertical breathing motion (VZ, br). In black and white, we show the 0 and 0.5 ρmax density contours of the potential, respectively. The response of the halo is compatible with the results shown in Figs. 1 and 4.

4 Response to a fixed potential

As we discussed in the introduction, we expect that the dominating dynamical process creating the DM spiral arms is the ‘forced response’ of the halo to the stellar spiral arms. To test this, we integrated halo orbits using a test particle set-up (TP; Appendix A.2), whereby halo orbits were integrated in an analytical galactic potential including spiral arms, isolating the forced response by removing self-gravity. By comparing the resulting density and kinematic signatures with those observed in the Iso-NBody model, we verified that the forced response remains identical even when self-gravity is absent.

In Fig. 5, we show the shape of the interaction signal in the four measurables we are studying. We see the clear correlation between DM and stars in density (first panel) and the slight lag of the DM spiral. For the kinematic modes, we see exactly the same patterns as for the Iso-NBody model. We observe clear quadrupoles in velocity, with positive ⟨V_R⟩ and ⟨V_ϕ⟩ on top of the spiral arms (second and third panels in Fig. 5), and amplitudes of ∼2 km s⁻¹, in the same range as the ones observed in the Iso-NBody model. In the vertical breathing map, we also observe the change of sign in V_{Z, br} on top of the spiral, with a breathing amplitude slightly smaller than the in-plane velocities. The signal approximately matches the one observed in Iso-NBody. Therefore, we conclude that the dominant dynamical process creating DM spiral arms in the Iso-NBody model is the forced response of the halo to the spiral arm potential.

We emphasize that no fine-tuning of the model parameters was required; we tested other combinations of spiral and MW-like model parameters and all produced a similar response. Thus, we conclude that the precise simulation parameters have a second-order effect on the overall response of the halo, making our conclusions robust to these variations.

4.1 Temporal response to the potential

The flexible TP set-up allows us to test the properties of the response that creates DM spiral arms, such as its temporal evolution. To do so, we conducted a different test in the same set-up as the TP simulation but in which, instead of a smooth rise in the spiral potential amplitude, we switched it on and off every 1 Gyr (Fig. 6). At each snapshot, we computed the strength (Σ₂/Σ₀) of the DM spiral arm (grey line in Fig. 6). With this set-up, we can study the temporal delay (or the impedance, here referring to the resistance of the system to rapid change) of the response. We observe that the response of the DM spirals to the sudden appearance or disappearance of a spiral potential is very fast (∼20 Myr, as is shown in the lower panel of Fig. 6). This timescale closely matches the local dynamical timescale given by the inverse epicyclic frequency (1/κ ∼ 24 Myr) at R = 8 kpc in the TP model. This correspondence supports our interpretation that the response is governed by local orbital dynamics.

We also checked whether the presence of a long-lived spiral potential can align some halo orbits with the disc plane, increasing the presence of DM particles with disc-like orbits. These disc-aligned DM particles would be more sensitive to ‘traditional’ spiral modes. To test this alignment, we extended the TP model to t = 6 Gyr. The fraction of DM particles with disc-like orbits increased by less than 1%, indicating that this alignment does not play a significant role.

Fig. 6 Strength (Σ2/Σ0, grey line) of the DM spiral arm at R = 8 kpc for the toy model. In this set-up, the potential amplitude is switched on and off every 1 Gyr. The long-dashed red line represents the shape of the potential (its height is arbitrary, for illustration purposes). Top: whole evolution. Bottom: detail of the strength of the potential around the three increasing steps of the potential (t = 1, 3, 5 Gyr). We observe that the potential takes ∼ 20 Myr to react to a change in the potential (impedance).

4.2 Response dependence on the orbital structure of the halo

As we reviewed in the introduction, the interactions between the components of the Galaxy have a resonant nature (Tremaine & Weinberg 1984). Therefore, we expect the halo particles that are closer to co-rotation with the stellar spiral to be more affected by it. In Petersen et al. (2016), the efficiency of torque transfer is explained using a formal framework based on the torque expressions derived in Weinberg (1985). Specifically, the LBK torque (Lynden-Bell & Kalnajs 1972) and non-linear torque (Eqs. (29) and (41) in Weinberg 1985) depend on rotation matrices, R(β), where β is the instantaneous angle between the angular momentum vectors (L_z) of the orbits and the disc plane. These expressions show that the torque transfer efficiency varies with β, with low-β orbits being more sensitive to the torque than high-β orbits.

In Fig. 7, we test this difference in response by splitting the halo particles according to cosβ in the TP and Iso-NBody simulations. At each snapshot, we evaluate the strength of the DM spiral arms for each sub-population. We clearly observe that the disc-like orbits in the halo (cosβ ∼ 1) are about twice as sensitive to the stellar spirals as orbits that are counter-rotating or perpendicular to the disc. It is worth noting, however, that even the DM particles that are counter-rotating with respect to the disc (dark lines in Fig. 7) are influenced to some extent by the stellar spiral arms.

There is evidence that the DM halo of the MW is triaxial (e.g. Han et al. 2022), a common feature in MW-like galaxies, as has been shown by cosmological simulations (e.g. Prada et al. 2019; Dillamore et al. 2022; Han et al. 2023). This could influence the strength of the DM spiral arms. To investigate the effect that rotating and/or flattened halos have on the response of the DM halo to the stellar spiral arms, we tentatively tried employing a novel AGAMA module. This tool enables the construction of self-consistent, rotating, and flattened DM halos by fitting a DoublePowerLaw DF to a target density profile². We then performed the same orbit integration as in the TP model for each DF-potential pair generated. Flattening the halos unavoidably introduced changes in the DF to maintain the self-consistency, which made comparing the response of the flattened halos nontrivial. On the other hand, prograde rotation clearly amplified the spiral arm response, while retrograde motion diminished it.

Fig. 7 Strength (Σ2/Σ0) of the DM spiral arms at R = 6 kpc for the TP (left) and Iso-NBody (right) models. The DM halo particles are split according to β, the instantaneous angle between the Lz of the orbit and the disc.

5 Ubiquity of the dark spiral arms

To complement the analysis, we explored the presence of dark spiral arms in other simulations: a pure N-body with a Sgr-like perturber (Sgr-NBody), an isolated hydrodynamical (Hydro), and a cosmological (Au27), all of which are described in Appendix A. In Fig. 8, we show the results for these models. In the Sgr-NBody model, we have enough resolution to study the velocity maps. However, the number of particles and the lower strength of the stellar spiral in the models Hydro and Au27 makes it very difficult to observe the kinematic signature of the DM spirals (Table A.1).

The Sgr-NBody model reveals a clear two-armed spiral pattern after the second and third passages of the perturber. These stellar spiral arms have a tidal origin, and thus their dynamical characteristics, such as pattern speed and kinematic signatures, are different compared to the Iso-NBody (e.g. Antoja et al. 2022). However, the observed response in the DM halo remains consistent with that seen in the Iso-NBody and TP models. We observe a clear DM spiral arm aligned with the stellar spiral arms and a similar kinematic signature to that of the Iso-NBody and TP models, with positive ⟨V_R⟩ and ⟨V_ϕ⟩ values on top of the spiral overdensity. In addition, we have studied the Sgr-NBody model with the same Fourier techniques as for the Iso-NBody, and obtained very similar results in terms of strength of the DM spirals and lag in phase. Thus, we conclude that the dynamical signatures of the forced response in the DM halo is independent of the dynamical origin of the stellar spiral arms.

To show the DM spiral arms in the Hydrodynamic simulation (Hydro), we selected an early snapshot (t = 5 Gyr) and one near the end of the simulation (t = 12 Gyr), which are the ones with a stronger two-armed spiral signature. In both cases, we observe signatures of the interaction between the dominant stellar spirals and the DM halo in the density space (bottom left groups of panels in Fig. 8). As in previous models, we observe the slight phase shift and the difference of about one order of magnitude in strength.

Finally, we show the z = 0.023 and z = 0 snapshots of an Auriga halo, 27 (Au27, Fig. 8, bottom right quadrant). Both snapshots show a small bar with two grand design spirals, with a lower strength than the previously studied models (δρ_r ∼ 0.4). Once again, we can identify the signatures of interaction in the density space quite clearly, shown by the spiral arm shaped overdensities in the halo, located on top of the stellar overdensity, with a slight lag in phase.

Fig. 8 Dark matter spiral arms in other simulations. Sgr-NBody (upper two groups of panels): stellar density (first column), DM density (second column), ⟨VR⟩ (third column), and ⟨Vϕ⟩ (fourth column) after the second and third passages of a Sgr-like perturber (first and second rows, respectively). Hydro (bottom left group): stellar and DM density for two snapshots with strong spiral arms. Au27 (bottom right group): stellar and DM density for the last two snapshots in the simulation, z = 0 and z = 0.023.

6 Discussion

6.1 Formation mechanism

Our results show that the DM halo responds to the perturbation caused by the stellar spiral arm potential. This is consistent with the predictions of the first-order perturbation theory (Mark 1976), which states that the response of the halo should be proportional to the mass of the spiral arms. The novelty is that with state-of-the-art simulations we can directly measure this interaction (Sects. 3 and 5), and we can go a step further in studying its nature. We find a DM spiral amplitude of ∼ 10% of the stellar spiral amplitude, a consistent phase lag (trailing with respect to the stellar spirals by ∼ 10^∘), and a characteristic kinematic signature in the halo. We also show that the DM halo responds to the presence of a spiral density perturbation within ∼ 20 Myr, a time compatible to the expected timescale, 1/κ, and that DM particles with disc-like orbits are more likely to be affected by the spiral potential.

Throughout this work, we have shown evidence that the DM spiral arms are formed through the forced response of the DM halo to the stellar spiral arms. First, the amplitude of the response of the halo is linear with respect to the amplitude of the stellar spiral arms in all the studied simulations. Second, with the TP (Sect. 4) we are able to reproduce the DM halo response to the spirals in both density and kinematic space. Since this model excludes self-gravity, it prevents the emergence of selfexcited dynamical modes in the halo, meaning that we see a pure response to the perturbation. In addition, using simple test particle set-ups, we show the fast response of the halo and the absence of significant subsequent alignment of the orbits with the disc. Therefore, we conclude that the signal observed in our simulations is predominantly due to the forced response of the DM halo to the stellar spiral arms.

6.2 Subdominant mechanisms

While the forced response dominates, secondary mechanisms may also contribute to the DM spiral arms formation. Here, we examine two of them: self-excited halo modes and stellar bar imprints in the DM halo.

In an isolated, non-cosmological DM halo, the formation of DM spiral-arm structures as self-excited gravitational modes is unlikely. This is mainly due to two factors: the lack of coherent rotation within the halo and its high velocity dispersion. These prevent the alignment necessary to form spiral arms, making coherent patterns unsustainable.

For the second mechanism, we note that here we have used a barred spiral galaxy as our fiducial model (Iso-NBody). Therefore, there is a possibility that the presence of a bar is contributing to the appearance of DM spirals. To test whether DM spiral arms could appear with ‘only’ a bar, we modified the test particle set-up and changed the spiral arms potential (Cox & Gómez 2002) for a Portail et al. (2017) bar (Appendix C). With this set-up, we are able to easily reproduce the shadow bar signature mentioned in the introduction. The bar produces a spiral signature in density in the DM halo. However, the strength is one order of magnitude lower than the one observed when there is a spiral potential. In the kinematic space, there is a spiral signature in ⟨V_ϕ⟩, with a smaller amplitude and a much larger pitch angle than the spirals in Iso-NBody. In ⟨V_R⟩ and V_{Z, br}, there is no signature. In addition, we do not see evidence of DM particles following invariant manifolds (e.g. Romero-Gómez et al. 2006; Voglis et al. 2006) in the DM halo. Therefore, from this set-up we do not expect a significant direct contribution from the stellar bar to the appearance of the DM spiral arms of interest. Nonetheless, bars have been proposed to be able to induce spiral arms in the stellar disc (Dobbs & Baba 2014; Sellwood & Masters 2022), which would generate DM spiral arms, again through a forced response.

6.3 Implications for the spiral arms morphology

We have shown that, in MW-like halos, the relative amplitude of the DM spiral arms is consistently one order of magnitude smaller than the stellar spiral arms (Figs. 1, 2, and 8). For instance, for a stellar spiral with a relative over-density of 20% (e.g. Widmark & Naik 2024), we predict a DM spiral arm with a relative over-density of 2%. This is translated into an increase of the total mass of the spiral arms of about 10% due to the overdensity in the DM halo, very similar to the predictions of the contribution of the “shadow bar” (Petersen et al. 2016). In the planar kinematic space, we observe a similar relation: while the signature of a stellar spiral arm is ≈10 km s⁻¹, the measured kinematic signature of the dark spiral is ≈1 km s⁻¹ (Figs. 4 and 8). Given these amplitudes, we expect the internal dynamics of the stellar spiral arms in the MW to dominate over the coupling with the DM halo.

However, these interaction effects are resonant, meaning that certain configurations amplify the DM spiral arm response. In Sect. 4.2, we show that DM particles on orbits aligned with the stellar disc are twice as sensitive to stellar spiral arms perturbations. We argue that this is consistent with analytical expectations: particles co-rotating with a spiral arm experience its gravitational influence over longer timescales. We also found that prograde rotating halos enhance the DM spiral response. While we have attempted to link the DM spiral response to the flattening of the halo, the interplay of the flattening with the other DM halo properties complicates the dynamical interpretation. A future statistical analysis of DM spiral amplitudes in cosmological simulations – correlating them with stellar spiral and halo properties – would be an interesting follow-up project.

7 Summary and conclusions

In this work, we present the first examination of the morphological and kinematic imprint of stellar spiral arms on the DM halo. We show that this interaction produces spiral-shaped overdensities in the DM halo. The appearance of a spiral overdensity in the halo is expected by a pure gravitational response. Our analysis confirms that the dominant mechanism forming these DM spiral arms is the forced response of the halo to the passage of the stellar spirals. Halo modes and bar forcing are not expected to form spirals with a density and kinematic signature consistent with the ones observed in the simulations. Our main findings and conclusions are the following:

• Whenever we observe stellar spiral arms (with Σ₂/Σ₀ ≳ 5%), we see an increase in the total mass of the spiral arms of about 10% due to the added mass from the DM halo response. The exact value is likely to depend on the properties of the spirals (e.g. pitch angle and pattern speed) and halo (e.g. rotation and anisotropy), and the region of the halo;

• The phase of the kinematic signature of the spiral arm in the DM also correlates with the stellar spirals phase. However, the DM spiral arms show a distinctive quadrupole kinematic signature. On top of the spiral overdensity, we consistently observe positive ⟨V_R⟩ and ⟨V_ϕ⟩, and a change in the sign of the vertical breathing motion;

• The amplitude of the planar kinematic signature in the DM halo is ∼10% that of the stellar component. However, the amplitude of the vertical breathing motion is similar in the stellar and DM components;

• The interaction between the stellar spiral arms and the DM halo is a resonant effect. The DM particles rotating with disc-like orbits are more susceptible to the torque produced by the spiral potential (Fig. 7). However, all the orbital families of DM particles are somehow affected by the interaction and produce spiral arm-shaped overdensities;

• We show the presence of DM spiral arms in a broad range of different state-of-the-art simulations, indicating that DM spiral arms are a common feature in simulated MW-like galaxies.

With this work, we contribute to the understanding of a surprisingly overlooked dynamical phenomenon: the interaction between spiral arms and the DM halo. The expected signal in density and kinematics is faint and convolved with the stellar spirals. However, this interaction occurs close to the Sun, offering a potential probe for exploring the nature of DM. If we were able to measure both the dynamical and baryonic mass distribution across a large region the MW disc, the difference between those would reveal the spiral substructure of the DM halo, rather than the expected first-order axisymmetric (or triaxial) distribution. Infrared data modelled by Drimmel & Spergel (2001) suggest a two-armed logarithmic spiral with a relative strength of 10%, yet different samples and techniques produce results that can vary from 10 to 20% (e.g. Eilers et al. 2020; Khanna et al. 2024). As for the dynamical mass, Widmark & Naik (2024), using the vertical Jeans equation and Gaia DR3 data, measured the distributions in the disc and inferred a local relative over-density of roughly 20%. The uncertainty in the current measurements of both the baryonic and dynamical mass makes the detection of DM spiral arms in the MW challenging but this could improve in the future. Once we reach that point, detailed predictions of the dark contributions will allow us to use the spiral arms as new tools to explore the nature of DM.

Hopefully, though, continuous improvements in the precision of the kinematic measurements of the stars and in the completeness of the samples used, aided by our ever-improving understanding of the baryonic dynamics, will allow us to confirm the predictions of this work. Once we reach that point, we shall be able to use the spiral arms as new tools to explore the nature of DM.

Acknowledgements

We thank the anonymous referee for her/his helpful comments. The authors thank Eugene Vasiliev for his support on the usage of AGAMA. We also thank Jason A. S. Hunt, the EXP collaboration, especially the Bars & Spirals working group, and the Gaia UB team for the valuable discussions that contributed to this paper. The OCRE simulations were run in virtual machines in the cloud provided by the OCRE awarded project Galactic Research in Cloud Services (Galactic RainCloudS). OCRE receives funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 824079. The M1_c_b (Hydro) simulation was run at the DiRAC Shared Memory Processing system at the University of Cambridge, operated by the COSMOS Project at the Department of Applied Mathematics and Theoretical Physics on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/J005673/1, STFC capital grant ST/H008586/1 and STFC DiRAC Operations grant vST/K00333X/1. DiRAC is part of the National E-Infrastructure. We acknowledge the grants PID2021-125451NA-I00, PID2021-122842OB-C21 and CNS2022-135232 funded by MICIU/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, by the “European Union” and by the “European Union Next Generation EU/PRTR”, and the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant CEX2019-000918-M. MB acknowledges funding from the University of Barcelona’s official doctoral program for the development of a R+D+i project under the PREDOCS-UB grant. TA acknowledges the grant RYC2018-025968-I funded by MCIN/AEI/10.13039/501100011033 and by “ESF Investing in your future”. JA is supported by the National Natural Science Foundation of China under grant No. 12233001, by the National Key R&D Program of China under grant No. 2024YFA1611602, by a Shanghai Natural Science Research Grant (24ZR1491200), by the “111” project of the Ministry of Education under grant No. B20019, and by the China Manned Space Project with No. CMS-CSST-2021-B03. JA thanks the sponsorship from Yangyang Development Fund. J.A. & C.L. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 852839). RG acknowledges financial support from an STFC Ernest Rutherford Fellowship (ST/W003643/1). OJA acknowledges funding from “Swedish National Space Agency 2023-00154 David Hobbs The GaiaNIR Mission” and “Swedish National Space Agency 2023-00137 David Hobbs The Extended Gaia Mission”. MSP is supported by a UKRI Stephen Hawking Fellowship. SRF acknowledges support from a Spanish postdoctoral fellowship, under grant number 2017-T2/TIC-5592.

Appendix A Simulations

In this section we describe the different simulation set-ups and the main characteristics of each model, summarized in Table A.1.

A.1 Isolated disc – Iso-NBody

The main analysis of this paper is based on a pure N-body model of a MW-like isolated galaxy. For the DM halo, a Hernquist (1990) profile with a scale length of R_h = 52 kpc and a mass of M_h=10¹² M_⊙ is used. The initial stellar disc is exponential, with scale length R_d = 3.5 kpc, scale height h_d = 0.53 kpc, and total mass M_d = 6 × 10¹⁰ M_⊙. Finally, a stellar bulge with M_b = 10¹⁰ M_⊙ and a Hernquist profile with a scale length of 0.7 kpc is included. The ratio σ_R/σ_z = 1.5 is fixed, so that the Toomre Q parameter is ∼ 1 in the disc (i.e. R ∈ [3, 10] kpc). The initial conditions were generated using GalIC (Yurin & Springel 2014), and evolved using the AREPO code (Weinberger et al. 2020) for 7.12 Gyr.

A.2 Test particle (TP)

To compare the N-body response with a controlled set-up, we used a suite of test particle simulations. We constructed a self-consistent model of a galaxy in equilibrium using AGAMA (Vasiliev 2019). We used a McMillan (2017) potential, and modelled the DM halo DFs as DoublePowerLaw (Binney 2014; Posti et al. 2015, implemented in AGAMA)³. Finally, we iterate the model using the self-consistent-modeling (Binney 2014; Piffl et al. 2015) module of AGAMA to obtain a self-consistent DF of the halo.

Once the potential and DF are set up, we include a Cox & Gómez (2002) spiral arm model (with a null total mass), $${\mathrm{\Phi }}_{\mathrm{s}}(R,\varphi ,z)=-4\pi G{\mathrm{\Sigma }}_0(t){\mathrm{e}}^{-R/R_{\mathrm{s}}}\sum _{n=1}^3\frac{C_n}{K_n{D}_n}\cos n\gamma {[\cosh \left(\frac{K_{n}z}{{\beta }_n}\right)]}^{-{\beta }_n},$$(A.1)

where Σ₀ is the central surface density, R_s is the scale radius, and C_n are the amplitudes of the harmonic terms. The functional parameters are: $$\begin{array}{rl}& K_n=\frac{nN}{R\sin \alpha },\\ & {\beta }_n=K_n{h}_{\mathrm{s}}(1+0.4K_n{h}_{\mathrm{s}}),\\ & \gamma =N[\varphi -\frac{\ln \left(R/R_{\mathrm{s}}\right)}{\tan \alpha }-{\mathrm{\Omega }}_{\mathrm{p}}t-{\varphi }_0],\\ & D_n=\frac{1}{1+0.3K_n{h}_{\mathrm{s}}}+K_n{h}_{\mathrm{s}},\end{array}$$(A.2)

where N is the number of arms, h_s is the scale height, α is the pitch angle, and ϕ₀ is the initial phase.

Following Dehnen (2000), the spiral amplitude is initialized at Σ₀(t) = 0 for t = 0, then increases smoothly between t = 0.5 Gyr and t = 2 Gyr as $${\mathrm{\Sigma }}_0(t)={\mathrm{\Sigma }}_{0,f}[\frac{3}{16}{\xi }^5-\frac{5}{8}{\xi }^3+\frac{15}{16}\xi +\frac{1}{2}],\xi =2\frac{t}{1.5}-1,$$(A.3)

and it remains at a fixed Σ_{0, f} for t > 2 Gyr until the end of the simulation at t = 6 Gyr.

In this run, we select parameters that visually approximated the shape and amplitude of the stellar spirals in the Iso-NBody model: N = 2, α = π/2, R_s = 4 kpc, C_n = [8/(3 π), 1/2, 8/(15π)], h_s = 1 kpc, Σ_{0, f} = 2.7 × 10⁸ M_⊙ kpc⁻², ϕ₀ = 0, and Ω_s = 20 km s⁻¹ kpc⁻¹. We integrate N_halo = 1.7 × 10⁸ particles for 6 Gyr, sampling every 30 Myr.

For completeness, we rerun the same model using a Portail et al. (2017) bar instead of spiral arm potential. Details of the model and results are presented in Appendix C. In brief, we do see signatures of interaction between the DM halo and the bar, but the signature in configuration and kinematic space is very different to the observed for the spiral arms.

A.3 Tidal pure N-body (Sgr-NBody)

To test the interaction between the spiral arms and the DM halo in a perturbed galaxy, we also studied a pure N-body model from the same OCRE suite as Iso-NBody with a Sagittarius (Sgr)-like dwarf spheroidal galaxy perturber. Minor mergers like these can excite spiral arms and bar formation (e.g. Toomre & Toomre 1972; Barnes & Hernquist 1992; Dobbs & Baba 2014; Sellwood & Masters 2022; Jiménez-Arranz et al. 2024), and are thought to be causing some of the substructure observed in the Gaia data (e.g. Antoja et al. 2018; Binney & Schönrich 2018; Laporte et al. 2019).

The MW model is very similar to Iso-NBody, with a Sgr-like satellite and a slightly hotter disc. The orbit of the perturber and the MW properties of the simulation are almost identical to the L2 model in Laporte et al. (2018b), while the Sgr-like perturber is modified to have a less concentrated DM halo, resulting in an increased mass loss history, more compatible to current observations (Laporte et al. in prep.).

From this simulation, we select two snapshots, at t = 4.94, and t = 6.06 Gyr, located slightly after the second and third pericentres of the perturber. These are the moments where the simulation contains stronger spirals (see Grion Filho et al. 2021 and Antoja et al. 2022 for a detailed description of the spirals in L2).

A.4 Isolated hydrodynamical (Hydro)

We complement the set-up of pure N-body models with a self-consistent, high-resolution simulation with star formation. The model we use is a hydrodynamic MW-like model, initially introduced in Fiteni et al. (2021), as model M1_c_b. This model develops through the cooling of a hot gas corona inside a DM halo. From this model, we select two snapshots with strong spiral \operatorname{arms}(t = 5, and t = 12 Gyr).

Roškar et al. (2012) investigate the general evolution of the spiral structure using a lower-resolution version of the same model. They identify different coexisting pattern speeds, some of which are comparable to those obtained in the model Iso-NBody (Ω_s ∼ 20 km s⁻¹ kpc⁻¹). This is confirmed in Debattista et al. (2025), which uses the high resolution model M1_c_b. Finally, Ghosh et al. (2022) used one of the snapshots we use (t = 12 Gyr) to study the breathing motions induced by its prominent spiral arms.

A.5 Cosmological (Au27)

Finally, we complement our analysis with a simulation from the Auriga project, in order to test the presence of DM spiral arms in cosmological state-of-art simulations. The Auriga project (Grand et al. 2017) consists of 30 zoom-in cosmological magnetohydrodynamical high-resolution simulations of MW-mass dark halos, that have been proved to form realistic spirals (Grand et al. 2021). Using the latest public data release of the augmented Auriga Project (Grand et al. 2024), we select the Original/3 version of the Au27 halo. From the Au27 halo, we select the two last snapshots (z = 0.023, and z = 0), which present a MW-like galaxy model with strong spiral arms.

Appendix B Methods

In this section we describe the methodological details of each step.

B.1 DM halo selection

We are interested in the response of the DM halo to the presence of stellar spiral arms. Throughout this paper, when referring to the DM halo, we will be applying the following cuts: $$|Z|<4\mathrm{k}\mathrm{p}\mathrm{c},R<30\mathrm{k}\mathrm{p}\mathrm{c}.$$(B.1)

These cuts are applied to increase the signal of the faint interaction we are studying. We tested modifications of the vertical cut over the range |Z| ∈ (0.2, 20) kpc. The overall morphology of the DM spiral arms remains very similar throughout these tests. Notably, when selecting DM particles close to the disk (i.e., |Z| ≲ 1 kpc), the DM spiral signal can increase by up to 40%. However, this increase does not alter our conclusions, as the amplitude remains approximately one order of magnitude lower than that of the stellar spiral. Additionally, a narrower vertical range reduces the number of particles, thereby degrading the signal quality. For these reasons, we adopted a compromise by maintaining the cut at |Z| < 4 kpc.

B.2 Measurables

Radially normalized density (δρ_r): To highlight the non-axisymmetric features, we normalize the particle count in each bin by the mean density computed at that radius, where the average is taken over all azimuthal angles (integrating over the vertical dimension). This is defined as: $$\delta {\rho }_r(x,y)=N(x,y)/\overline{N}(R)-1$$(B.2)

with $R=\sqrt{x^2+y^2}$ and N̄(R) representing the mean particle count at radius R calculated by averaging over all azimuthal angles. This normalization was used only in the density maps. The same result is achieved with the 0-th order mode of the Fourier analysis.

Mean radial velocity ⟨V_R⟩ and mean azimuthal velocity ⟨V_ϕ⟩: To study the overall dynamical properties of the samples, we focus on the mean velocities. In the case of the ⟨V_ϕ⟩, in the planar maps we subtract the mean V_ϕ at each radius (⟨V_ϕ⟩_r) for visualization and direct comparison with the other quantities. This subtraction is not done in the Fourier transform, since, again, V_c is captured by the 0-th order of the Fourier analysis.

Vertical breathing motion (V_{Z, br}): To study the vertical response to the spiral arms, we compute the breathing motion of the disc, which is the vertical asymmetry of the vertical velocity: $$V_{Z,\mathrm{b}\mathrm{r}}(x,y)=\frac{1}{2}(⟨V_Z(x,y,z_{+})⟩-⟨V_Z(x,y,z_{-})⟩)$$(B.3)

where z₊ = z ∈ [0, 4] kpc, and z₋ = z ∈ [−4, 0] kpc (Debattista 2014; Widrow et al. 2014). The vertical breathing motion will be positive in the regions where the disc is expanding vertically, and negative where the disc is contracting.

B.3 Fourier decomposition

We also use a Fourier analysis at each radius and snapshot as our quantitative measure to detect spirals. We split each snapshot in radial bins of 1 kpc with 100 bins in azimuth. For each bin in radius, we compute the 1D Fourier transform of the azimuthal profile of ρ, ⟨V_R⟩, ⟨V_ϕ⟩, and V_{Z, br}. Since we are studying systems with strong, two-armed spirals, the m = 2 mode is the dominant one and we use it to track the phase (ϕ₂) and amplitude (Σ₂) of the spiral. When studying the density, we refer to the relative amplitude of the m = 2 mode (Σ₂/Σ₀). For the kinematic modes, we refer directly to the amplitude of the 2nd mode (A₂). To compute the errors of the Fourier amplitude and phase, we use a window of ± 0.1 Gyr (20 snapshots) around the given time and compute the standard deviation of the value in the window.

We compute the pattern speed of the stellar and DM spiral arms in the simulation Iso-NBody. To compute the pattern speed (Ω_s) of the m = 2 mode at a given time, we select the ϕ₂ values in a ± 0.1 Gyr window. We unwrap the phase values to avoid the 2π discontinuity and compute the linear regression of t−ϕ₂ in the window. The slope of the regression is the pattern speed at the center of the window, once it is transformed to have km s⁻¹ kpc⁻¹ units. This method can be unstable if more than one spiral is present. However, in the simulation where we compute the pattern speed (Iso-NBody) there is only one pair of strong spiral arms and this method provides robust and accurate results.

B.4 Face-on images

We include face-on representations of the studied quantities to help the physical interpretation of the Fourier results. We smooth the images using a Gaussian kernel of width 0.3 kpc for the stellar component and 0.6 kpc for the DM component.

In the simulations with perturbers (Sgr-NBody and Au27), the motions in the DM halo (⟨V_R⟩, ⟨V_ϕ⟩, V_{Z, br}) are dominated by the presence of a dipole in velocity (Weinberg 1994; Heggie et al. 2020; Johnson et al. 2023). Since we focus on the study of the features with quadrupole signatures (two armed spiral arms), in Sect. 5 we show the centro-symmetrized component of the velocity maps of Sgr-NBody (average each pixel with its antipodal). This is $$\widetilde{A_{i,j}}=\frac{1}{2}(A_{i,j}+A_{n-i,n-j})\text{for all} i,j\in 0,\dots ,n-1$$(B.4)

where A is the 2D (n × n) grid of the studied quantity in Cartesian coordinates, i, j are the indices in the grid, and à its centrosymmetrized component. In practice, we are blindly removing all the odd modes of the image, highlighting the features we are interested in.

Appendix C Test particle simulation with bar

In this section we present the results of a test particle simulation with a purely barred potential to test if we can create DM spiral arms in this set-up. We start with the exact same set-up as the one described in Sect. A.2, constructing a self-consistent model of a galaxy in equilibrium using AGAMA. In that set-up, we then include a Portail et al. (2017) bar, rotating at Ω_b = 20 km s⁻¹ kpc⁻¹, and growing with the same law as Eq. A.3.

Fig. C.1 Response of the halo test particles to an analytical bar potential (Portail et al. 2017), rotating at a constant pattern speed. In black, we show the 0.1 and 0.5 ρmax density contours of the bar potential.

The results are shown in Fig. C.1. In δρ_r (first panel), we observe the clear presence of the shadow bar, and faint hints of the presence of spiral arms in the outer parts of the disc, with a strength below 0.1%. In ⟨V_ϕ⟩ (third panel in Fig. C.1), there is a clear spiral kinematic signature that could resemble the one observed in the pure N-body models. However, the amplitude of this kinematic signature is maximum in the inner part, and significantly smaller than the observed in the other simulations. Finally, we observe hints of a quadrupole structure in ⟨V_R⟩, very different to the signature observed in the other models, and we do not observe any signal in V_{Z, br}.

In conclusion, we are able to produce substructure in the DM halo when having only a bar, as expected from the results in Dillamore et al. (2023, 2024). In addition, this substructure is shaped as spiral arms. However, the signature in density and kinematics does not resemble the spiral signature we are studying.

References

All Tables

All Figures

Fig. 1 Spiral arms in the Iso-NBody simulation, at t = 7 Gyr. Left: surface density plot of the stellar component. Middle: radially normalized density (δρr) of the stellar component. The contours represent the regions where δρr = 0 (in black), and δρr = ±0.2 (in white). Right: δρr of the DM particles around the stellar disc (|Z| < 4 kpc). The contours are the same as the middle panel. We observe a clear spiral arm pattern in the DM halo, strongly correlated with the stellar spiral arms.

In the text

Fig. 2 Temporal evolution of the Iso-NBody simulation in the range R ∈ [5.5, 6.5] kpc. Top: relative amplitude of the Fourier mode (Σ2/Σ0), i.e. strength of the spiral arms for the stellar (solid ochre) and DM (dashed grey) components. We see that both spirals evolve in a correlated way, although the relative amplitude of the DM spirals is one order of magnitude smaller (notice the auxiliary axis, which refers to the amplitude of the DM spiral arms). Middle: pattern speed (Ωs) of the mode m = 2 structures in stars (solid ochre) and DM (dashed grey). Bottom: phase difference between the mode m = 2 of the stellar and DM structures (long-dashed red line). Shaded regions shows the dispersion of the phase lag. Once the spirals are strong enough (t ≳ 4.5 Gyr, vertical lines), we detect a constant pattern speed of Ωs ∼ 20 km s−1 kpc−1 for the stellar and DM spirals, and a constant phase lag of Δ ϕ2 ∼ 10∘, showing a strong correlation between the two spiral arms.

In the text

Fig. 3 Radial profile of the Iso-NBody simulation at t = 7 Gyr. We show the Fourier amplitude and phase for the stellar (solid ochre) and DM (dashed grey) components. First: amplitude of the 0th Fourier mode (Σ0), in units of solar mass. We see that, close to the midplane, the surface density of stars and DM is comparable. Second: relative amplitude of the Fourier mode (Σ2/Σ0) at all radii. The relative profile of the spirals is similar, with the relative amplitude of the DM spirals being one order of magnitude smaller than the stellar one (notice the auxiliary axis). Third: pattern speed (Ωs) of the mode m = 2. Fourth: phase lag (Δϕ2, long-dashed red line) between the DM and stellar components.

In the text

Fig. 4 Kinematic analysis of the Iso-NBody simulation, at t = 7 Gyr. The contours represent the regions where δρr = 0 (in black) and δρr = ± 0.2 (in white). We show ⟨VR⟩ (left column), ⟨Vϕ⟩ (middle column), and VZ, br (right column). First row: kinematic maps of the stellar component. Second row: kinematic maps of the DM component. Third row: radial profile of the kinematic amplitude of the Fourier mode (A2) for the stellar (solid ochre) and DM (dashed grey) components. Fourth row: radial profile of the phase difference between each kinematic mode (ϕ2) and the stellar spiral in density (ϕ2(ρStars)). The co-rotation radii of the spirals (RC R ∼ 10 kpc) is shown as a vertical red strip.

In the text

Fig. 5 Response of the halo test particles to an analytical spiral potential (Cox & Gómez 2002), rotating at a constant pattern speed. From left to right, we show the radially normalized density (δρr), the mean radial and azimuthal velocities (⟨VR⟩,⟨Vϕ⟩), and the vertical breathing motion (VZ, br). In black and white, we show the 0 and 0.5 ρmax density contours of the potential, respectively. The response of the halo is compatible with the results shown in Figs. 1 and 4.

In the text

Fig. 6 Strength (Σ2/Σ0, grey line) of the DM spiral arm at R = 8 kpc for the toy model. In this set-up, the potential amplitude is switched on and off every 1 Gyr. The long-dashed red line represents the shape of the potential (its height is arbitrary, for illustration purposes). Top: whole evolution. Bottom: detail of the strength of the potential around the three increasing steps of the potential (t = 1, 3, 5 Gyr). We observe that the potential takes ∼ 20 Myr to react to a change in the potential (impedance).

In the text

	Fig. 7 Strength (Σ2/Σ0) of the DM spiral arms at R = 6 kpc for the TP (left) and Iso-NBody (right) models. The DM halo particles are split according to β, the instantaneous angle between the Lz of the orbit and the disc.
In the text

Fig. 8 Dark matter spiral arms in other simulations. Sgr-NBody (upper two groups of panels): stellar density (first column), DM density (second column), ⟨VR⟩ (third column), and ⟨Vϕ⟩ (fourth column) after the second and third passages of a Sgr-like perturber (first and second rows, respectively). Hydro (bottom left group): stellar and DM density for two snapshots with strong spiral arms. Au27 (bottom right group): stellar and DM density for the last two snapshots in the simulation, z = 0 and z = 0.023.

In the text

	Fig. C.1 Response of the halo test particles to an analytical bar potential (Portail et al. 2017), rotating at a constant pattern speed. In black, we show the 0.1 and 0.5 ρmax density contours of the bar potential.
In the text