Probing Self-interacting Dark Matter with Disk Galaxies in Cluster Environments

For colliding DM halos, the microscopic nature of the interactions between the DM particles determines the dominant macroscopic signatures that will be observable in such systems. The two main macroscopic effects that have been considered in the literature are the evaporation rate of the smaller halo and an effective drag that each halo experiences as it moves through the environment of the other halo. For short-range interactions with isotropic cross sections, immediate halo evaporation is the most dominant effect (e.g., Markevitch et al. 2004 and Kahlhoefer et al. 2014). As a large fraction of collisions have a high momentum transfer, particles may be expelled from their host's potential, leading to subhalo evaporation. The observation of surviving low-mass subhalos therefore puts a stringent upper bound on the cross section of isotropic scattering, as well as the fraction of collisions that result directly in the expulsion of DM particles from the subhalo.

On the other hand, if we consider interactions where the cross section is velocity-independent and anisotropic—that is, there are many more small-angle scatterings than there are large-angle scatterings—frequent interactions are possible without completely disrupting the subhalo. In this scenario, individual interactions are usually unable to expel particles from the subhalo. However, there is a nonzero cumulative evaporation rate resulting from multiple interactions. As shown in Kahlhoefer et al. (2014), frequent and anisotropic self-interactions also lead to an effective "drag" force, and the rate of deceleration due to the drag force is comparable to the rate of cumulative evaporation. For these interactions, therefore, the macroscopic effects of the drag force can be comparable to the effects coming from evaporation (Kummer et al. 2017).

To understand the origin of the drag, we consider a two-particle interaction in the center-of-mass (COM) frame. If the scattering angle in this frame is θ, then in the direction parallel to the relative velocity of the two particles, the change in the velocity of one particle is given by

$$\begin{eqnarray}&&\delta {v}_{\parallel }=-v(1-\cos \theta ),\end{eqnarray} \tag{ 1 }$$

where v is the initial velocity of the particle in the COM frame. For isotropic interactions, $\cos \theta$ is drawn from a uniform distribution between −1 and 1. For anisotropic interactions of the form that we are interested in, $\cos \theta$ is drawn from a distribution that is peaked near 1 and −1 and has a very low probability of collisions where $\cos \theta \sim 0$. In particular, we use the same form of the differential cross section for the anisotropic interactions as was used in (Kahlhoefer et al. 2014):

$$\begin{eqnarray}&&\displaystyle \frac{d\sigma }{d{\rm{\Omega }}}\propto \displaystyle \frac{1+{\cos }^{2}\theta }{1-{\cos }^{2}\theta }.\end{eqnarray} \tag{ 2 }$$

If collisions are frequent, and if individual interactions only change the initial velocity by a small amount, we can integrate over all possible interaction angles to obtain the average change in the velocity in the parallel direction of a particle passing through a sea of other particles. In a time interval dt, the number of interactions is given by

$$\begin{eqnarray}&&{dN}=\displaystyle \frac{\rho }{{m}_{\mathrm{dm}}}(v\,{dt})\displaystyle \frac{d\sigma }{d{\rm{\Omega }}}d{\rm{\Omega }},\end{eqnarray} \tag{ 3 }$$

where $d\sigma /d{\rm{\Omega }}$ is the differential cross section, ρ is the ambient density, and ${m}_{\mathrm{dm}}$ is the mass of the DM particle. Using this, we can write the total change in the parallel velocity as

$$\begin{eqnarray}&&{{dv}}_{\parallel }=-\displaystyle \frac{\rho \,{v}^{2}\,{dt}}{{m}_{\mathrm{dm}}}\int \displaystyle \frac{d\sigma }{d{\rm{\Omega }}}(1-\cos \theta )d{\rm{\Omega }}.\end{eqnarray} \tag{ 4 }$$

Since we assume that the DM particles are indistinguishable, we follow Kahlhoefer et al. (2014) and define the momentum transfer cross section σ_T as

$$\begin{eqnarray}&&{\sigma }_{T}=4\pi {\int }_{0}^{1}\displaystyle \frac{d\sigma }{d{\rm{\Omega }}}(1-\cos \theta )\,d(\cos \theta ).\end{eqnarray} \tag{ 5 }$$

Once again, this expression makes sense only when the overall direction and velocity change per particle in the time interval dt is small compared to the incoming velocity. For our choice of the differential cross section, this is a safe approximation to make for sufficiently small time steps. The integral runs over the scattering angle $0\leqslant \theta \leqslant \pi /2$, since, for indistinguishable particles, a scattering angle above π/2 is the same as an equivalent scattering angle below π/2 but with a relabeling of the two outgoing particles. Using the expression for σ_T, Equation (4) reduces to

$$\begin{eqnarray}&&{{dv}}_{\parallel }=-\displaystyle \frac{\rho \,{v}^{2}\,{dt}}{2{m}_{\mathrm{dm}}}{\sigma }_{T}.\end{eqnarray} \tag{ 6 }$$

Therefore, this change in the parallel velocity can be written as an effective drag deceleration due to the anisotropic self-interactions as a particle moves through a sea of other DM particles, with

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{\mathrm{drag}}}{{m}_{\mathrm{dm}}}=-\displaystyle \frac{1}{2}\left(\displaystyle \frac{{\sigma }_{T}}{{m}_{\mathrm{dm}}}\right)\rho {v}^{2}.\end{eqnarray} \tag{ 7 }$$

We note that for isotropic cross sections where $d\sigma /d{\rm{\Omega }}$ is a constant, the momentum transfer cross section σ_T and the total cross section $\tilde{\sigma }$ are given by $\tilde{\sigma }=2{\sigma }_{T}$. Therefore, we can recast Equation (7) in terms of $\tilde{\sigma }$ so that our results can be compared directly to results for isotropic self-interactions. This yields

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{\mathrm{drag}}}{{m}_{\mathrm{dm}}}=-\displaystyle \frac{1}{4}\left(\displaystyle \frac{\tilde{\sigma }}{{m}_{\mathrm{dm}}}\right)\rho {v}^{2}.\end{eqnarray} \tag{ 8 }$$

We note that we assumed that the relative velocity between the two colliding halos, or the velocity of a subhalo falling into the host halo, is much larger than the velocity dispersion of the subhalo, and therefore the self-interactions between the subhalo and the particles of the host halo dominate over self-interactions between the particles of the subhalo itself.

We now describe the astrophysical system that is the focus of this study. Consider a spherically symmetric galactic subhalo that contains a disk galaxy seen edge-on by the observer. The COMs of both the subhalo and the disk initially coincide. The whole system moves along a trajectory r_sh with velocity v_sh through a host DM halo of density ρ_h (a galaxy cluster), assumed to be at rest. For typical orbital velocities and halo densities, Equation (8) leads to a nonnegligible contribution of SIDM to the motion of the subhalo, decreasing its infall acceleration and consequently distorting the galaxy disk. We present a schematic description of the system in Figure 1.

Zoom In Zoom Out Reset image size

To simplify our treatment of this problem, we make several approximations.

• 1.  
  The galaxy disk is made of stars only, and they are perfectly collisionless. We assume that the gravitational effect of the gaseous component of the disk is irrelevant to the dynamics.
• 2.  
  The galaxy disk experiences the gravitational attraction toward the DM subhalo in which it is contained, but the subhalo does not experience any attraction toward the stars. While this approximation is used for most of the simulations in this work, we do test its validity in Section 3.4.
• 3.  
  The DM subhalo is spherical and is not distorted while crossing the host. It is characterized by a static radial density profile, so we treat its trajectory semi-analytically. The same applies to the host halo. This approximation implies that any possible disruption of either halo during the merger or evaporation of the subhalo has a small effect on the stellar disk; the effective drag force is dominant.

With these approximations, plus Equation (8), the acceleration of the subhalo COM and a given star can be written as

$$\begin{eqnarray}&&{\ddot{{\boldsymbol{r}}}}_{\mathrm{sh}}=-G\displaystyle \frac{{M}_{{\rm{h}}}(\lt | {{\boldsymbol{r}}}_{\mathrm{sh}}| )}{{| {{\boldsymbol{r}}}_{\mathrm{sh}}| }^{3}}{{\boldsymbol{r}}}_{\mathrm{sh}}+\displaystyle \frac{1}{4}\left(\displaystyle \frac{\tilde{\sigma }}{{m}_{\mathrm{dm}}}\right){\rho }_{{\rm{h}}}{\dot{{\boldsymbol{r}}}}_{\mathrm{sh}}^{2},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}{\ddot{{\boldsymbol{r}}}}_{\star }^{i} & = & -G\displaystyle \frac{{M}_{{\rm{h}}}(\lt | {{\boldsymbol{r}}}_{\star }^{i}| )}{{| {{\boldsymbol{r}}}_{\star }^{i}| }^{3}}{{\boldsymbol{r}}}_{\star }-{{Gm}}_{\star }\displaystyle \sum _{i\ne j}\displaystyle \frac{({{\boldsymbol{r}}}_{\star }^{i}-{{\boldsymbol{r}}}_{\star }^{j})}{{| {{\boldsymbol{r}}}_{\star }^{i}-{{\boldsymbol{r}}}_{\star }^{j}| }^{3}}\\ & & -G\displaystyle \frac{{M}_{\mathrm{sh}}(\lt | {{\boldsymbol{r}}}_{\star }^{i}-{{\boldsymbol{r}}}_{\mathrm{sh}}| )}{{| {{\boldsymbol{r}}}_{\star }^{i}-{{\boldsymbol{r}}}_{\mathrm{sh}}| }^{3}}({{\boldsymbol{r}}}_{\star }^{i}-{{\boldsymbol{r}}}_{\mathrm{sh}}),\end{eqnarray} \tag{ 10 }$$

where ${{\boldsymbol{r}}}_{\star }^{i}$ is the 3D position of the ith star relative to the center of the host, ${M}_{{\rm{h}}(\mathrm{sh})}$ is the mass of the host (subhalo) enclosed in a given radius, and m_⋆ is the mass of a single star.

This system of equations is similar to that employed by Kahlhoefer et al. (2015), to which we add the mutual gravitational attraction between stars by introducing the summation term of Equation (10). In particular, disregarding that self-gravity term, direct numerical integration of Equations (9) and (10) was used by those authors for an approximate description of the offset created between a collisionless component and its subhalo's center.

Our choices for the halo profiles are shown in Figure 2. For both the subhalo and host, we use a Hernquist profile (Hernquist 1990) given by

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{M}_{\mathrm{dm}}}{2\pi }\displaystyle \frac{a}{r{(r+a)}^{3}},\end{eqnarray} \tag{ 11 }$$

where a is the scale factor and M_dm is the total mass of the halo. Following Springel et al. (2005), we determine a for each halo by matching their inner densities to an NFW profile (Navarro et al. 1996) with concentration c and virial radius r₂₀₀ using

$$\begin{eqnarray}&&a=\displaystyle \frac{{r}_{200}}{c}\sqrt{2[\mathrm{ln}(1+c)-c/(1+c)]}.\end{eqnarray} \tag{ 12 }$$

In our fiducial analysis, the host has a halo mass of ${10}^{15}\,{M}_{\odot }$ and concentration c = 5, while the subhalo has a mass of ${10}^{12}\,{M}_{\odot }$ and c = 8.

Zoom In Zoom Out Reset image size

The disk galaxy contained in the subhalo has an exponentially decaying radial profile and a squared hyperbolic secant vertical profile,

$$\begin{eqnarray}&&{\rho }_{\mathrm{disk}}(R,z)=\displaystyle \frac{{M}_{\mathrm{disk}}}{4\pi {z}_{0}{h}_{0}^{2}}\exp \left(-\displaystyle \frac{R}{{h}_{0}}\right){{\rm{sech}} }^{2}\left(\displaystyle \frac{z}{{z}_{0}}\right),\end{eqnarray} \tag{ 13 }$$

where h₀ is the scale length of the disk, z₀ is its scale height, and M_disk is its stellar mass. We choose ${M}_{\mathrm{disk}}=3\times {10}^{10}\,{M}_{\odot }$, ${h}_{0}\approx 3.5\,\mathrm{kpc}$, and a ratio of scale height to scale length of ${z}_{0}/{h}_{0}=0.1$ for the fiducial simulated galaxy disk. We also relate the velocity dispersions of the galaxy in the radial and vertical directions such that ${\sigma }_{z}=0.5{\sigma }_{R}$ and set ${\sigma }_{\phi }={\sigma }_{R}$ for the azimuthal direction $\hat{\phi }$. With these choices, the resulting disk is roughly compatible with a Milky Way–type galaxy (Binney & Merrifield 1998). The resulting galaxy is stable, and its Toomre Q parameter is greater than 1 at all relevant radii (Toomre 1964).

To set up initial conditions for the phase-space distribution of a galaxy disk, we use the publicly available code GalIC⁸ (Yurin & Springel 2014). GalIC iteratively minimizes a set of merit functions that measure the discrepancy between the target and realized spatial density field and velocity dispersions. It achieves this by using a Monte Carlo–style algorithm that forces the velocity distributions to have the correct second-order moments. We use $N={10}^{5}$ particles to populate the stellar disk and set the remaining parameters of GalIC to reproduce our fiducial scenario. The generated fiducial disk is close enough to an equilibrium state that it maintains its profile for a sufficiently long time.

To evolve the motion of the subhalo–disk system in time, we need to efficiently solve Equations (9) and (10). We use the publicly available code GIZMO⁹ (Hopkins 2015) in order to do so. GIZMO inherits some of its N-body algorithms from GADGET-2¹⁰ (Springel 2005) and GADGET-3 and solves gravity with an efficient tree method. It also allows for the use of external analytic gravitational potentials and forces, and we use that feature to reproduce the motion of the subhalo through the host. Both halos are introduced as analytic potentials and are not populated with particles, and the subhalo moves according to a prescribed trajectory while the host is fixed. We first integrate Equation (9) with the relevant set of initial conditions and a choice of $\tilde{\sigma }/{m}_{\mathrm{dm}}$ (0.5, $1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ or CDM), then feed to GIZMO a lookup table with this subhalo trajectory ${r}_{\mathrm{sh}}(t)$. We describe in Section 3.4 tests in which both halos and galaxy are treated as collections of particles in the simulation and evolve under self-gravity, without analytic shortcuts, to confirm that this method is valid.

As depicted in Figure 1, the subhalo-plus-disk system starts its infall at the virial radius of the host, roughly $1600\,\mathrm{kpc}$ away from its center. To simulate realistic cases, we use subhalo orbits with several choices of pericenters, defined as the closest approach distances from the center of the host. We set the pericenter by initially giving the subhalo–disk system a tangential velocity whose magnitude is a fraction of the circular velocity about the host. While we experimented with different values, we focus on pericenter distances of 200 and 300 kpc for our analysis, since these show appreciable SIDM effects and are not dominated by tidal distortions due to the host that could break our initial assumptions. The ratio of mean pericenter to cluster virial radius (apocenter) is about 1:6, as shown in simulations by Ghigna et al. (1998). We note that a more recent study by van den Bosch et al. (2017) found a mean pericenter that is larger by a factor of around two with respect to Ghigna et al. (1998). Our choices thus reflect a sizable fraction of orbits of subhalos inside clusters.

We also expect the orientation of the galaxy disk at its closest passage to the host to change the observable morphological distortions. This orientation is defined as the angle between the COM velocity of the galaxy disk and the normal to its plane. We study four different scenarios: orientation angles of 0°, 225, 45°, and 675 (0° is a face-on passing through the host, and 90° would be a passing along the diameter of the disk, which is not interesting for our analysis). Note that in all of these cases, the galaxy is still edge-on for the observer. In our simulations, each disk is evolved separately, so no simultaneous interaction with other galaxies is present. We describe the results for various choices of pericenters and orientation angles in Sections 3.1 and 3.2.

Consider a single star of mass m_⋆ in a nearly circular orbit around the disk center. To determine the disk warping, we are interested in the star's motion in the direction perpendicular to the galaxy plane, z(t). Under the approximation that the drag force is constant across the entire subhalo, and considering that the warp is small enough such that the disk potential is not severely disturbed, we can describe the motion of a star in the reference frame of the subhalo as

$$\begin{eqnarray}&&\ddot{z}=\left(\displaystyle \frac{{F}_{\mathrm{drag}}}{{m}_{\mathrm{dm}}}\right)-\displaystyle \frac{{F}_{\mathrm{sh}}\sin \theta }{{m}_{\star }}-\displaystyle \frac{\partial {{\rm{\Phi }}}_{\mathrm{eff}}(R,z)}{\partial z},\end{eqnarray} \tag{ 14 }$$

where the first term is given by Equation (8) and comes from the change in reference frames, and the effective potential of the axisymmetric stellar disk is ${{\rm{\Phi }}}_{\mathrm{eff}}(R,z)={\rm{\Phi }}(R,z)+{L}_{z}^{2}/2{R}^{2}$, where R is the radial distance from the star to the center of the disk along its plane and L_z is its conserved angular momentum. The term ${F}_{\mathrm{sh}}\sin \theta$ is the projection of the subhalo gravity onto the direction perpendicular to the plane of the galaxy, with θ being the angle between the galaxy plane and the vector that goes from the galaxy center to the displaced star. One could also add to Equation (14) a term that corresponds to the tidal force caused by the host, but that component is at least an order of magnitude weaker than the gravity of the subhalo at the scale radius of the disk.

Recalling the approximation that a star does not wander too far from the plane of the disk, we have $\sin \theta \approx z/R$. Using the spherical symmetry of the subhalo and the epicycle approximation (Binney & Tremaine 2008) to Taylor expand the potential up to second order in z-derivatives, we can write

$$\begin{eqnarray}&&\ddot{z}=\displaystyle \frac{1}{4}\left(\displaystyle \frac{\tilde{\sigma }}{{m}_{\mathrm{dm}}}\right){\rho }_{{\rm{h}}}{v}_{\mathrm{sh}}^{2}-\displaystyle \frac{{{GM}}_{\mathrm{sh}}(R)}{{R}^{2}}\left(\displaystyle \frac{z}{R}\right)-\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}z,\end{eqnarray} \tag{ 15 }$$

where the second derivative of the potential is evaluated at the equilibrium position z = 0 on the plane of the galaxy. Describing the subhalo by a Hernquist profile, the cumulative mass within the radius R is given by (Hernquist 1990)

$$\begin{eqnarray}&&{M}_{\mathrm{sh}}(R)={M}_{\mathrm{dm}}\displaystyle \frac{{R}^{2}}{{(R+a)}^{2}},\end{eqnarray} \tag{ 16 }$$

where a is the scale radius of the halo, equivalent to that of Equation (11). To treat the potential, we use a fully analytic expression that approximates the 3D potential of an axisymmetric disk with characteristic radius h₀ and characteristic height z₀, given by Miyamoto & Nagai (1975):

$$\begin{eqnarray}&&{\rm{\Phi }}(R,z)=-\displaystyle \frac{{{GM}}_{\mathrm{disk}}}{{[{R}^{2}+{({h}_{0}+\sqrt{{z}^{2}+{z}_{0}^{2}})}^{2}]}^{1/2}}.\end{eqnarray} \tag{ 17 }$$

A final approximation is that most of the warp occurs in a relatively short time span, becoming maximal near the pericenter of the subhalo trajectory by the host. We verify with the simulations that this is especially true for the mild warps with which we are concerned in this derivation. With that in mind, the drag force acts as a "kick" on the star under consideration. Using Equations (16) and (17) on Equation (15) and keeping the dominant orders in R, the motion of a point on the disk is described by

$$\begin{eqnarray}&&\ddot{z}(t={t}_{\mathrm{imp}})={\left[\displaystyle \frac{1}{4}\left(\displaystyle \frac{\tilde{\sigma }}{{m}_{\mathrm{dm}}}\right){\rho }_{{\rm{h}}}{v}_{\mathrm{sh}}^{2}\right]}_{t={t}_{\mathrm{imp}}}-{\omega }^{2}z,\end{eqnarray} \tag{ 18 }$$

where ${t}_{\mathrm{imp}}$ is the time at which the subhalo is at pericenter and

$$\begin{eqnarray}&&{\omega }^{2}=\displaystyle \frac{{GM}_{{\rm{d}}{\rm{m}}}}{R{(R+a)}^{2}}+\displaystyle \frac{\left(1+{\textstyle \tfrac{{h}_{0}}{{z}_{0}}}\right){R}^{2}}{{[{({h}_{0}+{z}_{0})}^{2}+{R}^{2}]}^{5/2}}{GM}_{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{k}}},\end{eqnarray} \tag{ 19 }$$

such that the system resembles a driven harmonic oscillator with characteristic frequency ω.

Equation (18) suggests two things: that the gravitational pull of the subhalo and disk act together as a restoring force, opposing the warp, and that this distortion is oscillatory. Especially at larger R, the squared frequency given by Equation (19) goes as ${R}^{-3}$, so the restoring force is smaller at large radii and particles closer to the disk edge are less tightly bound gravitationally. The galaxy disk thus gets warped due to differential offsets along the disk, its longitudinal section becoming "U-shaped," and the warp shape may oscillate around the initially thin configuration. Note that the edge-on disk initially bends forward toward the direction of motion of the galaxy, as depicted in Figure 1. We indeed find the warping effect in the simulations described in the following sections, not only for the mildly warped disks but also for the intense warps that perturb the galaxy potential beyond our analytic approximations.

We can use Equation (18) to obtain a theoretical estimate of the maximum amplitude of the warp. Consider the fact that very shortly after the drag-induced "kick," the galaxy reaches its new distorted equilibrium position, and that the kick takes a finite but small amount of time. In that case, we can find the equilibrium displacement z by setting $\ddot{z}=0$ on Equation (18). The maximal displacement obtained this way is

$$\begin{eqnarray}&&{z}_{\max }(R)=\displaystyle \frac{{\rho }_{{\rm{h}}}{v}_{\mathrm{sh}}^{2}}{4{\omega }^{2}}\left(\displaystyle \frac{\tilde{\sigma }}{{m}_{\mathrm{dm}}}\right),\end{eqnarray} \tag{ 20 }$$

with ω² given by Equation (19). Using physical parameters similar to our fiducial system on Equation (20) for the disk and subhalo dimensions and mass, we estimate the magnitude of the distortions in units of kpc in Table 1 for the pericenters of interest. More explicitly, we use ${h}_{0}=3.5\,\mathrm{kpc}$, ${z}_{0}=0.1\,{h}_{0}$, ${M}_{\mathrm{dm}}={10}^{12}\,{M}_{\odot }$, ${M}_{\mathrm{disk}}=3\times {10}^{10}\,{M}_{\odot }$, and a = 30 kpc. We use the local DM density of the host at pericenter, and the subhalo velocity is around 3200 km s⁻¹.

Table 1.  Analytic Prediction of the Magnitude of the Displacement z_max as a Function of Disk Radius, in Units of kpc, for the Fiducial Galaxy Scenario when $\tilde{\sigma }/{m}_{\mathrm{dm}}=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$

	$R={h}_{0}$	$R=2{h}_{0}$	$R=3{h}_{0}$
200 kpc	0.48	1.05	2.39
300 kpc	0.19	0.42	0.95

Note. These results match the simulations with reasonable accuracy.

Download table as:  ASCIITypeset image

While Equation (20) itself suggests a way of fitting the shape of the distortion, we choose a slightly different method in Section 3.1, which also accounts for asymmetric warps. An added effect, not accounted for in the previous derivation, is an "S-shaped" warp resembling an integral sign. The presence of an S-shaped warp is a standard dynamical feature that arises quite generally from tidal interactions with close neighbors (Binney 1992; Reshetnikov & Combes 1998; Shang et al. 1998).

Previous works have focused on the overall COM offset between the DM halos and their luminous components but overlooked the differential offset from an extended body such as a disk galaxy. In particular, while an offset between the centroids of light and total matter in cluster mergers must reach tens of kpc to be measurable, we find that displacements smaller than 1 kpc between the subhalo and the stellar disk can lead to clear morphological distortions, as explored in Section 3.1.

However, perturbations that bend galaxy disks tend to damp out quickly, in a few rotation times (Binney & Tremaine 2008), so the warps are not expected to be permanent features. Our simulations show that, after the warping phase, the stellar disk also does not revert back to a thin plane. Rather, in its final and relaxed state, sufficiently far from the densest parts of the host such that the drag force is again negligible, the stellar orbits spend significant time away from the galactic plane. We interpret this final state as an enhancement of the thickness of the initially thin disk. In Section 3.2, we quantify how much thicker the disk contained in the SIDM halo gets with respect to the baseline CDM case.

As expected from our analytical arguments, we indeed find that galaxy disks in our simulations warp forward, oscillate and warp backward, and then finally relax in a puffier state than at the beginning. To quantify the intensity of the warps, we fit a third-order polynomial to the 2D projection of the edge-on stellar disk, $z(R)={{aR}}^{3}+{{bR}}^{2}+{cR}+d$, where z is the coordinate perpendicular to the plane of the disk and R is the radial distance along the disk's plane from its center. Such a polynomial is able to capture the initial and final states of the disks, when it should be well fit by a straight line, as well as the U and S warps, which will force a fit to large b and a coefficients, respectively. This choice of metric is simple and convenient for our present purposes, but we note that others are available and have been used in the literature (e.g., Vikram et al. 2013). To fit the polynomial when there is some orientation angle, we first rotate the galaxy so that it aligns with the vertical axis. Since the warps oscillate and then vanish, the fitting coefficients reproduce this behavior by getting increased in magnitude and then eventually decreased back to nullity. While the prediction from Equation (20) assumed a symmetric warp and employed the approximation of a constant disk potential, which breaks down for strong warps, we do check that the predicted values in Table 1 are compatible with the simulations.

We take the magnitude of coefficient b as a simple proxy of the U-warp intensity. To quantify some degree of theoretical uncertainty, we also obtain the error bars for b from the covariance matrix of the polynomial fitting. For a distant galaxy, one would hardly observe stars significantly further than a few scale radii from the disk center, due to the rapidly decaying flux caused by the lower stellar density. To avoid overestimating the observable warp intensity, we restrict the polynomial fit only to particles that lie within a sphere of radius equal to 3h₀ from the center of the subhalo, where h₀ is the scale radius of the galaxy disk.

For a better understanding of the distortions, Figures 3 and 4 show simulation snapshots of the evolution of the warps with time. With all stellar disks starting from the same initial conditions, we show in Figure 3 a snapshot of the $\tilde{\sigma }/{m}_{\mathrm{dm}}\,=0.5\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ and $\tilde{\sigma }/{m}_{\mathrm{dm}}=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, 300 kpc–0° simulation, as well as the CDM case. The galaxies in the SIDM halos exhibit prominent U warps. Notably, the stellar disk that inhabits the CDM subhalo displays a modest S warp at the same time frame. It is clear from the figure that the galaxy in the standard CDM subhalo should not display a significant b coefficient as given by the polynomial fit. Notice that the $\tilde{\sigma }/{m}_{\mathrm{dm}}=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ panel in Figure 3 provides a check for the analytical prediction of Table 1. There is overall good agreement between the values in the 300 kpc row of that table and what is presented in Figure 3. For $R=3{h}_{0}\gtrsim 10\,\mathrm{kpc}$, the predicted displacement of $0.95\,\mathrm{kpc}$ is in concordance with the visual displacement of around $1\,\mathrm{kpc}$ in the figure.

Zoom In Zoom Out Reset image size

Once the subhalo's gravity dominates over the drag again, the warp oscillates backward, as represented by Figure 4 by the snapshots at a time step around 200 Myr after those of Figure 3. Similarly, the S-shaped warp on the galaxy belonging to the CDM subhalo changes orientation.

Zoom In Zoom Out Reset image size

Finally, at a later stage, around 1 Gyr after the first warp and pericenter passage, the bending modes have decayed, and the remaining effect is an enhancement of the thickness of the stellar disk. Such disk heating is also a feature of gravitational interactions in a standard CDM scenario. Our simulations suggest, however, that a population of disks under SIDM should be thicker than their CDM counterparts due to the extra heating caused by the initial warping. For a quantitative analysis, we determine a metric for the disk thickness and apply it to simulation snapshots in the following section, where the thickness enhancement is made clear in Figure 8.

In Figure 5, we show snapshots of the simulations for both choices of pericenter and all orientation angles at the maximal warping instant according to the polynomial fit for the case $\tilde{\sigma }/{m}_{\mathrm{dm}}=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$. All panels show galaxies that inhabit SIDM subhalos, and the disk initial conditions are the same. The units on each axis correspond to the separation of the edge-on disks from the centers of their respective subhalos. The general trend is that the 0° orientation angle (face-on onto the host) leads to more prominent, symmetric warps, while other inclination angles are somewhat asymmetric around the center of the subhalo. Also, a shorter pericenter leads to stronger distortions as it probes denser parts of the host halo.

Zoom In Zoom Out Reset image size

Figure 5 also shows that, for inclinations that are not exactly face-on (orientation angle differing from zero), the stellar distribution can become significantly skewed within the disk plane. This skewness is another potential signature of DM self-interactions, in addition to the U-shaped warps and enhanced disk thickness. Our simulations suggest that SIDM-induced skewness can persist within satellite disks for times similar to the durations of the U-shaped warps. We defer a detailed study of SIDM-induced skewness to future work and instead focus on quantifying warp signatures and disk thickness.

In Figure 6, we show the evolution of the quadratic coefficient of the polynomial fit as a function of time for the 0° orientation angle in our two pericenters of interest. Negative values of b correspond to a parabolic U shape in the forward direction of motion, while positive b values bend the disk backward. In both panels, we see oscillations of the coefficient that eventually vanish as the disk aligns again with the subhalo center. The error bands are multiplied by a factor of 10 on those plots for better visualization.

Zoom In Zoom Out Reset image size

We define the threshold for the detection of the warps to be $| b| \gt 0.003$, which is roughly an order of magnitude larger than the values this same coefficient reaches in the baseline CDM case due to fluctuations. For a very symmetric warp with $z(R)\approx {{bR}}^{2}$, this threshold value of b corresponds to a perpendicular displacement of approximately 0.03 and 0.3 kpc as measured around the disk's scale radius and three times the scale radius, respectively. After the first backward warp (the first positive bump of the quadratic coefficient b in Figure 6), the galaxy disk is significantly thicker than how it was initialized, looking "puffy" rather than warped. Conservatively restricting ourselves to only one forward plus one backward warp, we find from the curves in Figure 6 that the period during which the warp could be detected in the least-affected 0° orientation case ($\tilde{\sigma }/{m}_{\mathrm{dm}}=0.5\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, 300 kpc) is around 0.1 Gyr, while the most affected disk ($\tilde{\sigma }/{m}_{\mathrm{dm}}\,=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, 200 kpc) is detectably warped for at least 0.4 Gyr. The survival time of the U-shaped warp ultimately translates into a potentially observable sample size, which we further discuss in Section 4.

An important difference then becomes clear: for SIDM cross sections within the range that we have explored ($0.5\lt \tilde{\sigma }/{m}_{\mathrm{dm}}\lt 1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$), U-shaped stellar disks are expected to be quite general on galaxies that make a fast passing through a dense DM environment, such as a galaxy cluster. In standard CDM, such stellar disk warps should be rare when compared to their S-shaped counterparts. We expect that further work to measure the intensity of warps and how frequently they are found in observations can place tight constraints on $\tilde{\sigma }/{m}_{\mathrm{dm}}$.

To complement Figures 5 and 6, we present Figure 7, comparing different orientation angles for the case $\tilde{\sigma }/{m}_{\mathrm{dm}}\,=0.5\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ with a 300 kpc pericenter. We again show error bands that are multiplied by a factor of 10 for visualization. The trend is similar to that shown in Figure 5, with steeper angles displaying a less intense—but still detectable—warp signal when compared to the symmetric 0° case.

Zoom In Zoom Out Reset image size

Once the transient warps decay and the galaxy reaches a new equilibrium configuration, we find that the resulting stellar disks are thicker in the SIDM subhalos when compared to their standard CDM counterparts. This happens because stars in SIDM galaxies get scattered to outer orbits following the warped phase. To quantify the thickness enhancement, we use simulation snapshots at around 3 Gyr, when the disk is around its furthest distance from the host. We verify that the thickness is fairly constant around that time frame. On the 2D projection of the disk as seen by the observer, we first select only the particles that lie inside a box of radial range $R=[-3{h}_{0},3{h}_{0}]$ and height range $z=[-5{z}_{0},5{z}_{0}]$, where h₀ ≈ 3.5 kpc is the initial scale length of the galaxy disk and ${z}_{0}=0.1{h}_{0}$ is the initial scale height. In the cases with some nonzero orientation angle, we first rotate the galaxies by the known angle to align them with the vertical axis and then select particles. This cut is supposed to represent only particles that could be detected as part of the stellar disk in actual observations in the same way we select particles before fitting the third-order polynomial.

We then make histograms of the height of these selected particles with respect to the disk plane, which in the edge-on view of the observer is simply given by the horizontal coordinate z. As a primary metric for the disk thickness, we use the standard deviation of the resulting distribution of particles ${{\rm{\Delta }}}_{\mathrm{disk}}\equiv \sqrt{\langle {z}^{2}\rangle -{\langle z\rangle }^{2}}$. While the initial disk's transverse profile is, by construction, very well fit by a sech² function, according to Equation (13), from which a scale height can be obtained, we make this choice because the final, thicker disks are not well fit by that expression. As the thickness differences between CDM and SIDM play a more significant role at around z ≳ 3z₀, we also fit an exponential profile proportional to ${e}^{-z/\tau }$ as a secondary metric, where τ is analogous to a new disk scale height. This fit is restricted within the range $[2{z}_{0},4{z}_{0}]$, where the profiles are well described by a single exponential length scale τ. For a comparison with a control disk, we also measure both thickness metrics in a simulated system that does not go through the host halo but rather evolves in isolation in a CDM subhalo. After 3 Gyr, the vertical profile of this isolated disk remains consistent with the initial conditions due to GalIC's success in finding a quasi-equilibrium configuration. We find ${{\rm{\Delta }}}_{\mathrm{disk}}^{\mathrm{isol}.}\,=0.24\,$ and ${\tau }^{\mathrm{isol}.}=0.20\,$ kpc in this control system.

Figure 8 represents the obtained histograms at a time of 3 Gyr for different cross sections in the 200 and 300 kpc pericenter cases, with orientation angle 0°. The vertical bands in both panels show the values of the thickness metric, and the uncertainty represented by their width corresponds to an averaging over ∼1 Gyr around the turnaround time. For the 200 kpc case, we obtain ${{\rm{\Delta }}}_{\mathrm{disk}}(\mathrm{CDM})=0.32\,\pm 0.006\,\mathrm{kpc}$, ${{\rm{\Delta }}}_{\mathrm{disk}}(\tilde{\sigma }/{m}_{\mathrm{dm}}=0.5\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1})=0.39\pm 0.008\,\mathrm{kpc}$, and ${{\rm{\Delta }}}_{\mathrm{disk}}(\tilde{\sigma }/{m}_{\mathrm{dm}}=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1})=0.46\pm 0.010\,\mathrm{kpc}$. Notice that the CDM disk also gets thicker due to tidal interactions alone, which is a potential systematic uncertainty. In both panels, a larger cross section leads to a thicker disk, but the discrepancy is more prominent for smaller pericenters, which probe a higher number of SIDM interactions, a trend similar to that of the U-shaped warps. Dashed lines correspond to the secondary metric, the exponential fit, which captures mostly the difference in the histogram tails. For the 200 kpc case, we obtain $\tau (\mathrm{CDM})=0.34\,\mathrm{kpc}$, $\tau (\tilde{\sigma }/{m}_{\mathrm{dm}}=0.5\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1})\,=0.49\,\mathrm{kpc}$, and $\tau (\tilde{\sigma }/{m}_{\mathrm{dm}}=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1})=0.64\,\mathrm{kpc}$. Both pericenters show a slightly stronger SIDM signal in the τ metric than in Δ_disk.

We note that including particles away from the limits $| z| \lt 5{z}_{0}$ in the calculation of Δ_disk increases the enhancement effect. That is mostly due to the fact that the disks in SIDM halos actually have a considerable number of particles spread out to much further heights, while the CDM histogram decays very quickly and has little contribution at $z\gt 5{z}_{0}$, as can be seen from Figure 8. This implies that our selection of particles is conservative for observational purposes.

Zoom In Zoom Out Reset image size

The survival time for the enhanced thickness is considerably longer than that of the U warps. We find that during at least 1 Gyr, the histograms presented in Figure 8 remain very stable despite the orbital motion of particles. In fact, due to the collisionless nature of stars, it is not expected that the stellar disk will completely revert back to a thin plane, since that would require some energy to be radiated away.

We thus conclude that this is another interesting difference between morphologies in SIDM and CDM in cluster environments. While in both scenarios, a disk that has made a fast passage through a dense DM environment should become thicker, the currently allowed SIDM cross sections can lead to an enhancement of that effect. A consequence of this conclusion is that field galaxies that have never been inside cluster environments should be relatively thinner, both in CDM and SIDM, than those inside clusters.

A caveat is worth noting: one should expect the cluster environment to be more disturbing to disks than our single-galaxy simulations can capture. This could lead to a statistically increased thickness of edge-on disks even in the absence of SIDM, for instance, due to tidal interactions with many potential close neighbors other than the host itself.

As expected from our theoretical estimates in Section 2.4, our results in the previous sections are somewhat sensitive to the initial choice of the fiducial system—the host and subhalo masses and their concentration, the initial galaxy thickness, etc. Here we consider variations of these physical parameters. In the next section, we also vary the "unphysical" parameters of the simulation, for instance, the number of particles and the gravitational-softening scale. We also test, to some extent, the initial simplifying approximations made explicit in Section 2.2. To accomplish these tests, we simply run a new suite of simulations, changing, whenever possible, only one relevant parameter at a time.

The physical characteristics of the system, like halo profiles and galaxy length scales, are ultimately a set of nuisance parameters that one would marginalize over, in some sense, when trying to infer the SIDM cross section $\tilde{\sigma }/{m}_{\mathrm{dm}}$ from actual observations. That is due to the fact that the results shown in Sections 3.1 and 3.2 are expected to have a dependence on these parameters.

When we vary the host and subhalo profiles, the resulting warp strength follows the trend expected from the analytical derivation in Section 2.4. Increasing the concentration of the subhalo, at fixed mass, from the fiducial choice c = 8 to c = 12 yields a more cuspy halo that consequently produces a stronger restoring force on the distorted galaxy disk. The disk warp is expected to be less intense. Indeed, the resulting warp curve for the 200 kpc, $\tilde{\sigma }/{m}_{\mathrm{dm}}=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ case then reaches a minimum coefficient $b=-0.024\,{\mathrm{kpc}}^{-1}$, approximately 30% less intense than the fiducial case (see Figure 6). The same trend is expected when the subhalo mass is increased from the fiducial ${M}_{{\rm{s}}{\rm{h}}}={10}^{12}\,$ to $5\times {10}^{12}\,{M}_{\odot }$, and we confirm that the most intense warp coefficient is 50% smaller than the fiducial scenario. Conversely, for a less massive subhalo with ${M}_{\mathrm{sh}}=5\times {10}^{11}\,{M}_{\odot }$, the maximum warp is 50% larger than the fiducial case in Figure 6. In these tests, the galaxy and host halo properties were kept the same as the fiducial choice.

We run another set of tests modifying the host profile parameters. We simultaneously change the initial infall distance of the disk galaxy to match the virial radius of the new host. The other galaxy parameters remain fixed at their fiducial values, and again we look at the case of 200 kpc, $\tilde{\sigma }/{m}_{\mathrm{dm}}=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ with a 0° orientation angle. Decreasing the host mass from 10¹⁵ to $5\times {10}^{14}\,{M}_{\odot }$ at fixed concentration reduces the maximum warp by around 60%, resulting in b = −0.014 kpc⁻¹ at maximal warping. Also, further reducing it to ${10}^{14}\,{M}_{\odot }$ significantly reduces the maximum warp intensity to b = −0.001 (kpc)⁻¹, below our defined detection limit of $| b| \gt 0.003\,{(\mathrm{kpc})}^{-1}$. This again follows the trend described in Section 2.4: reducing the mass at fixed concentration makes the host less dense at its center, leading to a weaker SIDM drag force from Equation (8) and, consequently, subtler distortions.

We also modify the disk's structural parameters. Reducing the stellar mass from $3\times {10}^{10}\,$ to ${10}^{10}\,{M}_{\odot }$ and changing the initial scale height from ${z}_{0}=0.1{h}_{0}$ to ${z}_{0}=0.05{h}_{0}$ and ${z}_{0}=0.2{h}_{0}$ did not produce changes larger than 10% on the warp intensity.

In all of the tests above regarding the physical parameters of the system, we also look at the effect on the disk thickening. With respect to our fiducial case, all variations caused only marginal changes on Δ_disk below ±10%, to the final thickness of the stellar disk in the $\tilde{\sigma }/{m}_{\mathrm{dm}}=1.0\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ SIDM subhalo (see Figure 8).

We proceed to test the numerical "unphysical" simulation parameters and our initial approximations. The (constant) softening length used in our simulation runs was 50 pc, considerably shorter than the length scales of the effects we have studied. We performed several runs with softening as small as 25 pc and as large as 200 pc and found no significant difference in the obtained results. We were focusing especially on the evolution of the disk thickness when this parameter changed, since spurious two-body scatterings can lead to heating of the galaxy disk followed by an increase in its thickness.

To make sure our results are not an artifact of integration accuracy, we also experiment with GIZMO's tree construction frequency and the maximum allowed size of the time steps. We find no significant difference in our measurements when the tree updating is made faster by an order of magnitude, and the same for when the maximum time step allowed is made lower by one and two orders of magnitude. We also increase and decrease the number of disk particles by a factor of 2 and find subpercent differences in the warp intensity of thickness measurements. We thus conclude that our results are fully converged.

Finally, we test approximations (1) and (3) in Section 2.2 by implementing simulations with both the subhalo and host described by dynamic gravitating N-body particles. We run these N-body simulations by modifying GADGET-2 (Springel 2005) to include SIDM with an anisotropic scattering cross section (A. Banerjee et al. 2018, in preparation). The disk was again generated using GalIC. The simulations used 2 ×10⁷ particles for the host halo of mass ${10}^{15}\,{M}_{\odot }$, 2 × 10⁴ for the subhalo with mass ${10}^{12}\,{M}_{\odot }$, and 10⁴ in the $3\times {10}^{10}\,{M}_{\odot }$ disk. In general, anisotropic self-interactions will give rise to both evaporation of the subhalo and drag. To isolate the effect of drag on the galaxy disk system, we run a drag-only implementation where each simulation particle is considered to be representative of an ensemble of microscopic particles and the net drag force can be evaluated using Equation (7). In this case, we find that strong warps are produced in the galaxy disk even when the outer regions of the subhalo are largely distorted by tidal interactions. This is consistent with our assumption that the warp is sensitive to the subhalo profile within the scales of the disk radius, which does not change significantly due to tidal interactions. The drag-only case by itself is not energy-conserving. We run another set of simulations where we implement the full physical picture with drag and evaporation. We treat interactions between simulation particles as if they represent actual microscopic particles—i.e., we use the total cross section σ to decide if two neighboring particles interact or not and then choose the scattering angle θ using the probability distribution for θ from the form of $d\sigma /d{\rm{\Omega }}$ (A. Banerjee et al. 2018, in preparation). Since these interactions change the velocities of the interacting particles in both the direction parallel to the relative velocity of the particles and the direction perpendicular to it, this method naturally incorporates the effects of both drag and subhalo evaporation arising from self-interactions.

Since evaporation affects the entire subhalo profile, we find that for the cross section used in this work (Equation (2)), the warping can be suppressed in our fiducial 300 kpc pericenter case. Conversely, the disk thickening is enhanced due to mass loss from the subhalo. However, when performing a full simulation for the idealized scenario of a completely radial trajectory for the subhalo, we find measurable disk warps, even though their magnitude is smaller than the drag-only case. It must be noted that the amount of evaporation and drag can be affected differently, as both depend on the nature of $d\sigma /d{\rm{\Omega }}$ (Kummer et al. 2017). These simulations are currently in the process of development, and it will be of interest to do a detailed analysis on how different cross sections affect the warping of galaxy disks.