# Galaxy stability within a self-interacting dark matter halo

###### Abstract

This paper investigates spheroidal galaxies comprising a self-interacting dark matter halo (SIDM) plus de Vaucouleurs stellar distribution. These are coupled only via their shared gravitational field, which is computed consistently from the density profiles. Assuming conservation of mass, momentum and angular momentum, perturbation analyses reveal the galaxy’s response to radial disturbance. The modes depend on fundamental dark matter properties, the stellar mass, and the halo’s mass and radius. The coupling of stars and dark matter stabilises some haloes that would be unstable as one-fluid models. However the centrally densest haloes are unstable, causing radial flows of SIDM and stars (sometimes in opposite directions). Depending on the dark microphysics, some highly diffuse haloes are also unstable. Unstable galaxies might shed their outskirts or collapse. Observed elliptical galaxies appear to exist in the safe domain. Halo pulsations are possible. The innermost node of SIDM waves may occur within ten half-light radii. Induced stellar ripples may also occur at detectable radii if higher overtones are excited. If any SIDM exists, observational skotoseismology of galaxies could probe DM physics, measure the sizes of specific systems, and perhaps help explain peculiar objects (e.g. some shell galaxies, and the growth of red nuggets).

###### keywords:

dark matter — galaxies: elliptical and lenticular, cD — galaxies: kinematics and dynamics — galaxies: structure^{†}

^{†}pagerange: Galaxy stability within a self-interacting dark matter halo–C

^{†}

^{†}pubyear: 2011

## 1 Introduction

For decades it has been observed that visible matter in galaxies and clusters is insufficient to bind them gravitationally (e.g. Zwicky, 1937; Freeman, 1970; Roberts & Rots, 1973). The extra gravity of invisible “dark matter” (DM) is invoked, but its properties are uncertain, with many particle physics conjectures offering candidates. On Earth, DM particle detection experiments differ in their findings: some report annually modulated signals (e.g. Bernabei et al., 2008; Aalseth et al., 2011a, b; Angloher et al., 2012); others see nothing above their expected background (e.g. Angle et al., 2008; Aprile et al., 2010; Akerib et al., 2006; Ahmed et al., 2009). Reconciling them may require extra physics: inelastic or composite forms of dark matter, or extra internal degrees of freedom (e.g. Smith et al., 2001; Chang et al., 2009; Alves et al., 2010; Kaplan et al., 2010, 2011)

Large-scale cosmic filaments and voids emerge in -body computer simulations using gravity alone (e.g. Frenk, White & Davis, 1983; Melott et al., 1983; Springel, Frenk & White, 2006) but there is observational tension concerning smaller, self-bound systems. The predicted abundances of dwarf galaxies exceed those seen in many types of environments (e.g. Klypin et al., 1999; Moore et al., 1999; D’Onghia & Lake, 2004; Tikhonov & Klypin, 2009; Carlberg, Sullivan & Le Borgne, 2009; Zavala et al., 2009; Peebles & Nusser, 2010). The biggest simulated subhaloes appear too massive to match real satellite galaxies (Boylan-Kolchin, Bullock & Kaplinghat, 2011, 2012). Collisionless DM haloes develop singular central density cusps (e.g. Gurevich & Zybin, 1988; Dubinski & Carlberg, 1991; Navarro, Frenk & White, 1996). On the contrary, observations of diverse galaxy classes indicate or allow a dark density profile that is smooth and nearly uniform within a kpc-scale inner core, (e.g. Flores & Primack, 1994; Moore, 1994; Burkert, 1995; Salucci & Burkert, 2000; Kelson et al., 2002; Kleyna et al., 2003; Goerdt et al., 2006; Gentile et al., 2004; de Blok, 2005; Thomas et al., 2005; Kuzio de Naray et al., 2006; Gilmore et al., 2007; Weijmans et al., 2008; Oh et al., 2008; Nagino & Matsushita, 2009; Inoue, 2009; de Blok, 2010; Pu et al., 2010; Murphy, Gebhardt & Adams, 2011; Memola, Salucci & Babić, 2011; Walker & Peñarrubia, 2011; Jardel & Gebhardt, 2012; Agnello & Evans, 2012; Salucci et al., 2012; Lora et al., 2012). Galaxy cluster haloes should be centrally compressed by the steep density profile of cooling gas, but nontheless many examples are definitely or possibly cored (e.g. Sand, Treu & Ellis, 2002; Sand et al., 2008; Ettori et al., 2002; Halkola, Seitz & Pannella, 2006; Halkola et al., 2008; Voigt & Fabian, 2006; Rzepecki et al., 2007; Zitrin & Broadhurst, 2009; Newman et al., 2009, 2011; Richtler et al., 2011).

To remedy (or defer) these DM crises, many astrophysicists invoke energetic “feedback” by starbursts or active galaxies: explosions and outflows violent enough to flatten the dark cusps and denude small subhaloes of baryons (e.g. Navarro, Eke & Frenk, 1996; Gnedin & Zhao, 2002; Mashchenko, Couchman & Wadsley, 2006; Peirani, Kay & Silk, 2008; Governato et al., 2010). These events would have happened early enough to be difficult to observe today. The details would involve radiative shocks in magnetised multiphase interstellar media that are thermally and dynamically unstable. The relevant scales are numerically unresolvable in present simulations. Results depend on phenomenological recipes for sub-grid processes, and the implementations are more varied and elaborate than any DM physics. For reasons of resolution alone, the proposition that feedback cures the DM crises cannot be decisively evaluated for at least another decade. Regardless, feedback may be energetically incapable of removing cusps whilst leaving realistic stellar populations in dwarf galaxies (Peñarrubia et al., 2012).

An alternative approach extends standard particle physics phenomenologies: taking the dark matter crises as signatures of dynamically important physics in the dark sector. In “SIDM” theories, dark matter is self-interacting yet electromagnetically neutral. Pressure support or thermal conduction could prevent cusps, while drag and ablation reduce subhaloes (e.g. Flores & Primack, 1994; Moore, 1994; Spergel & Steinhardt, 2000; Firmani et al., 2000; Vogelsberger, Zavala & Loeb, 2012; Rocha et al., 2012). SIDM interactions entail direct inter-particle scattering (with specific crossection ) possibly enhanced by short-range Yukawa-type forces, or a long-range “dark electromagnetism” (e.g. Ackerman et al., 2009; Loeb & Weiner, 2011). Collisionless plasma supported by “dark magnetism” is conceivable. Alternatively, the SIDM might consist of boson condensates of a scalar field (e.g. Peebles, 2000). Plausible early-universe evolution and cosmic structures emerge in cosmogonies with fluid-like SIDM (e.g. Moore et al., 2000; Stiele, Boeckel & Schaffner-Bielich, 2010, 2011) or scalar field dark matter (e.g. Matos & Arturo Ureña-López, 2001; Suárez & Matos, 2011). There were early fears that SIDM haloes would form sharper cusps via gravothermal collapse (e.g. Burkert, 2000; Kochanek & White, 2000). Later analyses put this beyond a Hubble time, and found two safe regimes that suppress thermal conductivity: weaker scattering (rarely per orbit); and stronger (fluid-like) cases with tiny mean free paths (Balberg, Shapiro & Inagaki, 2002; Ahn & Shapiro, 2005; Koda & Shapiro, 2011). We will refer to the and regimes as “weak-SIDM” and “strong-SIDM” respectively.

The ubiquity of dark cores supports SIDM, but there have been circumstantial counter-arguments. It was argued that observed halo ellipticities imply low collisionality (e.g Miralda-Escudé, 2002; Buote et al., 2002). Such limits cannot apply to systems that are elliptical due to tidal disturbance, rotation or non-radial pulsations. Newer SIDM simulations reject the spherical prediction and imply a looser limit (Peter et al., 2012). Other papers inferred tight limits on SIDM collisionality, from the growth of supermassive black holes (: Ostriker, 2000; Hennawi & Ostriker, 2002) and the survival of cluster members ( or , Gnedin & Ostriker, 2001). However those investigations were in a way self-fulfilling, because they assumed cuspy initial conditions, with a strong thermal inversion at the centre. This choice was inconsistent with cored profiles expected from SIDM.

Gravitational lensing models of merging galaxy clusters give conflicting signs about DM physics: sometimes the inferred DM seems collisionless, following galaxies and separating from gas (, e.g. Clowe et al., 2006; Randall et al., 2008). In other systems DM seems collisional, following the gas or separating from galaxies (e.g. Mahdavi et al., 2007; Jee et al., 2012; Williams & Saha, 2011). The derived limits depend on assumed projected geometries and histories of unique systems; these are each open to specific counter-explanations and rebuttals. For instance if the acclaimed “Bullet Cluster” were truly a direct-hit merger then ram pressure and shocks would quench star formation. This is not observed (Chung et al., 2009).

Saxton & Wu (2008) modelled galaxy clusters in which the halo could possess extra heat capacity (e.g. composite SIDM particles). Stationarity in the presence of gas inflows requires a central accretor exceeding some minimum mass. For the range of DM properties giving plausible cluster and core profiles, a part of the inner halo teeters towards local gravitational instability, implying likely collapse to form supermassive black holes of realistic sizes. Saxton & Ferreras (2010) modelled elliptical galaxies comprising a polytropic halo plus coterminous collisionless stars. Observed stellar kinematics of elliptical galaxies were fitted by the same class of DM that gives the most realistic clusters.

This paper is a further theoretical study of spheroidal galaxies. Their stability properties will imply constraints on the end-points of galaxy evolution. If galaxies do possess SIDM haloes then their oscillatory modes will open the future possibility of skotoseismology — acoustic diagnosis of dark matter haloes that pulsate. The paper’s organisation is as follows: §2 formulates the stationary model and perturbations; §3 describes the oscillatory solutions and instabilities; §4 discusses applications and interpretations; §5 concludes.

## 2 Formulation

This paper investigates spherical galaxies comprising self-interacting dark matter and collisionless stars in near-equilibrium configurations. The SIDM halo is polytropic, adiabatic, non-singular and radially finite (§2.1). The system is isolated: unconfined and unexcited by any external medium. Stellar orbits are isotropically distributed, in an observationally motivated density profile (§2.2). Oscillations and instabilities are calculated in response to radial perturbations (§2.4–§2.6). Mass displacements are non-uniform, and at any location the two constituents may be displaced with differing amplitudes and phases. The perturbed variables are non-singular at the origin, and the outer boundary of the halo is free (§2.7).

### 2.1 dark matter halo profile

For each mode in which ideal collisional gas particles are free to move, the energy per particle is where is Boltzmann’s constant, is temperature, is the particle mass and is a characteristic (one-dimensional) velocity dispersion. If SIDM gas has effective thermal degrees of freedom that are locally in equipartition (due to collisionality or collective effects), then the specific heat capacities at constant volume and pressure are and . The ratio of specific heats

(1) |

From the first law of thermodynamics, adiabatic processes imply a polytropic equation of state for the dark pressure:

(2) |

or (equivalently) the mass density is where is a pseudo-entropy, and is a generalised phase-space density. In adiabatic processes, and remain constant.

The micro-physical meanings of “” deserve further remarks. Point-like particles moving in three spatial dimensions have , due to translational kinetic energy. Composite particle substructure raises (e.g. rotational and vibrational modes). A diatomic gas has because two parts are moving in 3D space, at a fixed separation. Incompressible fluids are more constrained and . For a relativistic or radiation dominated plasma, . If particles were massive enough (short wavelength) to traverse hidden spatial dimensions, then extra translational modes would raise . An isothermal medium has .

Polytropic fluids also emerge in some non-classical theories of dark matter. Scalar field and boson condensate models usually imply , but possibly other values depending on an index in the interaction potential (e.g. Goodman, 2000; Peebles, 2000; Arbey, Lesgourgues & Salati, 2003; Böhmer & Harko, 2007; Harko, 2011a, b; Chavanis & Delfini, 2011; Robles & Matos, 2012). This paper is agnostic about which particle theory causes the value of , which should ultimately be constrained by observations.

The radial structure of a spherical halo is determined by its momentum equation,

(3) |

The stellar mass enclosed at radius is (see §2.2). The corresponding enclosed dark mass is obtained by numerical integration of

(4) |

In hydrostatic equilibrium, the condition (3) becomes

(5) |

Neglecting any compact central mass (e.g. black hole), a non-singular boundary condition applies at the origin (). The equilibrium dark halo is adiabatic, and its pseudo-entropy is spatially uniform. This may occur in nature if a galaxy halo is well mixed, or if is set by universal constants of the SIDM particle.

For a halo of finite mass truncates () at a finite outer radius . For other , any finite mass equilibrium model has infinite radius: it boils off into the void, or fails to condense from the cosmic background in the first place. The high limit () is a generalised Plummer (1911) model, with infinite radius and finite mass. Haloes may also seem radially infinite in some SIDM simulations with (e.g. Rocha et al., 2012). This may be due to ongoing supersonic infall of neighbouring objects, or because the scattering is weak enough to leave the outskirts collisionless. If the SIDM “strength” is increased (greater scattering , or long-range “dark magnetism” effects) then polytropic behaviour affects more of the halo. Beyond some threshold strength, a definite surface must occur. Rocha et al. (2012) note a temperature drop at large radii in their largest simulation. This is attributable to incipient truncation.

The value of determines the natural ratio of core to system radii: for greater the core is smaller (e.g. compare cases in Figure 1). If dark mass dominates, the circular orbit velocity profile () peaks at a radius which occurs at when and respectively. In these cases, the dark density index drops to at radii ; and index at . In other words, a strong-SIDM core with is radially smaller than a tenth of a core. Appendix A lists standard radii for other cases. Adding stars distorts these ratios but the general trend with persists. Figure 1 (right panel) shows how a stellar component can shrink and steepen the dark core. For comparison, the “scale radius” of simulated CDM haloes (Navarro, Frenk & White, 1996) is , and . For the Burkert (1995) cored profile with scale radius , the key radii are , , and . Simulated SIDM cores fitted by Rocha et al. (2012) () appear consistent with Burkert’s parameterisation. The ratios are similar for Burkert models and polytropes, but less similar for other .

These ideal core sizes are upper limits. In collisional SIDM, weaker delays the formation of full-sized cores. Current SIDM simulations focus on the case, which predicts cores larger than those observed, unless the scattering crossection is kept small (, e.g. Yoshida et al., 2000; Arabadjis, Bautz & Garmire, 2002; Katgert, Biviano & Mazure, 2004; Vogelsberger, Zavala & Loeb, 2012; Rocha et al., 2012). The alternative is to vary : the range of results in galaxy clusters with realistic –kpc cores (Saxton & Wu, 2008). In principle, haloes with could produce cores of realistic radius for arbitrarily large . Such values also provide small enough to be consistent with observations of large flat regions in disc galaxies’ rotation curves. The choice of is therefore an astrophysically preferable subdomain, but for the sake of generality this paper considers a wider interval, .

### 2.2 stellar mass profile

The assumed stationary distribution of the stars resembles the empirical Sérsic (1968) form, where the surface density at projected radius is

(6) |

and is the “Sérsic index.” The “effective radius” encircles half the starlight. Assuming that the stellar mass-to-light ratio is constant, there exist exact deprojections of (6) that express the 3D density profile (e.g. Mazure & Capelato, 2002), but these involve exotic special functions. We will instead use the approximate deprojection of Prugniel & Simien (1997), with density and mass profiles

(7) |

(8) |

where is the lower incomplete gamma function. The index determines and (Lima Neto, Gerbal & Márquez, 1999; Ciotti & Bertin, 1999; Márquez et al., 2000). Adopting to match the classic de Vaucouleurs (1948) profile, the values are and . For this standard choice, the total stellar mass is , and that within the half-light radius is .

The frequency scale of the stars will be abbreviated as

(9) |

Given the properties of large elliptical galaxy like M87, and (Murphy, Gebhardt & Adams, 2011). For a compact elliptical such as M32, and (Mateo, 1998; Choi, Guhathakurta & Johnston, 2002). For “red nuggets” at high redshift, and (e.g. Toft et al., 2012; van de Sande et al., 2011). If stars were the only component present, the keplerian period at would be , and the angular frequency .

It is convenient to write the density gradient index as

(10) |

noting that near the origin. In central regions a modified radial coordinate, , is useful to soften terms involving the density,

(11) |

as the factor cancels the singular part of in any ODEs.

The stellar velocity dispersions may in principle differ in the three orthogonal directions: radial (), polar () and azimuthal (). We will consider spherical models in which the transverse velocity dispersions are equal, . The radial and transverse velocity dispersions are related by defining an anisotropy profile,

(12) |

(e.g. Osipkov, 1979; Merritt, 1985; Mamon & Łokas, 2005). The radial “pressure” of stars is

(13) |

and its profile is calculated from a spherical Jeans equation for momentum, expressed in a frame co-moving with the mean radial velocity () of a mass “packet” of stars (Jeans, 1915; Binney & Tremaine, 1987),

(14) |

where is the gravitational field.

### 2.3 acoustic & orbital timescales

The acoustic crossing time of the dark matter from the origin to the outer surface is obtained by radial integration of

(15) |

For the stars, the analogous “acoustic” time is:

(16) |

Assuming stationarity, then substituting (10) and (14) simplifies (16): e.g. if then . From dimensional considerations, the crossing time should scale with the halo radius as . Appendix B shows examples of the exact relations. When calculating the natural modes of the galaxy, we should expect oscillations to scale in proportion to the angular frequencies and . Once the eigenfrequencies () of a galaxy are known, one can compare them to local keplerian frequencies,

(17) |

In regions where , perturbations should only affect stellar orbits adiabatically. Wherever , tidal shocking may occur when disturbances exceed the linear regime.

### 2.4 perturbative displacements and gravitation

For each of the variables describing local properties of the dark matter or stars, we define a first order perturbation. For instance for a density variable we have

(18) |

where the subscript “0” denotes the stationary solution; the “” part describes the response to perturbation; the function describes the perturbation in relative terms (e.g. ). Perturbed variables such as and are complex. The small constant is the scale of perturbation; Higher order terms are discarded (). The frequency of perturbation is complex () with the real part describing oscillations and the imaginary part describing the tendency to growth or decay.

As in Eddington (1918) we shall consider lagrangian perturbations. For a mass shell that sits at radius in equilibrium, we write its displacement and velocity under perturbation:

(19) |

(20) |

The dark matter and stars at any may be displaced differently, so we treat them in two slightly different lagrangian frames, and distinguish their coordinates carefully:

(21) | |||||

(22) | |||||

(23) | |||||

(24) |

As the perturbed variables are complex, the mass constituents may be displaced with different phases (i.e. it is possible that ). Time derivatives are straightforward, e.g. for density,

(25) |

Expressions for the perturbation of each radial gradient term depend on whether we consider the comoving frame of the stars or the dark matter. Let’s abbreviate “” and “” for “” and “” subscripts. Mass conservation requires that , and then in the frame we have the operators

(26) |

The gravitational field provides the only coupling between components. Consistent evaluation of its perturbations (in both frames) requires special care. The stellar and dark shells that rest at in the stationary solution will feel the same field strength there () but the perturbed versions of differ since the shells oscillate through displacements with different phase and amplitude. Define

(27) |

For newtonian gravity, we have with being the total mass contained within in equilibrium. Poisson’s equation applies at all orders of perturbation, which implies

(28) |

Combining (26) and (27) gives the perturbed part of as felt by the component. Denoting gives

which has the advantage of being nonsingular and at the origin.

### 2.5 dark matter perturbations

The dark matter perturbations are considered adiabatic ( is spatially and temporally constant). The equation of state links the pressure and density perturbations:

(30) |

Mass conservation then yields

(31) | |||||

(32) |

The lagrangian equation for momentum (3) decomposes under (26) to give the hydrostatic condition,

(33) |

and a perturbation equation, writable in two forms:

(34) | |||||

(35) |

with the abbreviation

(36) |

The presence of the stars is felt through the gravitational perturbation obtained from (LABEL:eq.mass.perturbed).

### 2.6 stellar perturbations

For the collisionless stellar component, the perturbed mass conservation equation resembles the dark matter version, but without an equation of state to link density and pressure.

(37) | |||||

(38) |

Next, we define perturbed forms of the radial and transverse velocity dispersions of the stars,

(39) |

(40) |

where and are the perturbations of the radial and transverse velocity dispersions respectively. The stationary models in this paper assume isotropy ( everywhere) although orbital anisotropies occur transiently during perturbations. Nonetheless let’s retain terms in the general formulation. If the perturbations are radial then there is no torque. If the specific angular momentum of a mass element is conserved during a radial displacement then each local value of stays constant, eliminating a variable:

(41) |

Expressing the radial pressure perturbation as , provides a convenient relation,

(42) |

The Jeans equation (14) separates into stationary and perturbed parts,

(43) |

(44) | |||||

Given another Jeans equation for evolution of the radial stellar pressure (Appendix C),

(45) | |||||

and performing the perturbation expansion yields an equation in the stationary variables (which reduces to a tautology, ) plus an equation in perturbed variables. The latter, first order equation is:

(46) |

Substituting (37) into (46) leads to an algebraic relation,

(47) |

### 2.7 boundary conditions

At the origin, by symmetry, there cannot be any oscillatory displacement of mass, and therefore

(48) |

In relative terms the offset is free; any small nonzero value of is acceptable at the inner boundary. Regularity of the mass equations (32) and (38) then gives:

(49) |

(50) |

Combining (47) and (49) implies

(51) |

and using the central pressure from the stationary model. Analogously, the central pressure oscillation of dark matter is . Assuming that the stars and dark matter are comoving in the dense inner regions requires at the origin. This completes the set of inner boundary conditions.

As at the halo surface, (35) usually tends towards singularity. To avoid this, the natural modes must satisfy an outer boundary condition:

(52) |

Condition (52) is only satisfied at eigenfrequencies. The stellar mass profile lacks any finite outer boundary. Therefore there are no outer boundary conditions applicable to the perturbed variables.

### 2.8 numerical integration and solution

To calculate the stationary profiles,
trial values are set for the pseudo-entropy constant
and the central velocity dispersion of dark matter
().
The system of coupled differential equations
(4)–(5)
is integrated numerically as an intial value problem in radius
(using the rk4imp and rk8pd integrators from the
Gnu Scientific Library^{1}^{1}1http://www.gnu.org/software/gsl/).
In practice it is necessary to change the independent variable
(the integration coordinate) partway through the integral.
Near the origin a softened radial coordinate
(11)
improves accuracy in the stellar cusp.
In the outskirts it is more practical to use
or
as the independent variable,
and integrate to a standard limit at
or .
(The chosen limit doesn’t change any results,
as long as it is consistent across all runs.)
At the outer boundary (radius )
the total dark mass is recorded
().

At a higher level of the program, an amoeba search routine reiterates this integration for different trial values of , seeking to minimise the difference between the obtained and target values chosen by the user. Usually these targets are set to collect a sequence of solutions of different but equal , i.e. describing galaxies of identical mass composition but different compactness.

Once the values are known for a desired stationary model, the stellar pressure profile (43) is integrated inwards from infinity and tabulated finely. Trial values are chosen for the complex frequency of perturbation, . The stationary model’s ODEs are integrated again simultaneously with the coupled ODEs for the perturbed variables: (LABEL:eq.mass.perturbed), (31), (32), (34), (35), (37), (38), (43) and (44). When the halo surface is reached () the outer boundary condition (52) is tested there, to assess the closeness of to an eigenvalue. Using a complex test score,

(53) |

a root-finding routine refines the trial values until . To seek eigenfrequencies on the real and imaginary axes of , a bisection root-finder is sufficient. To find eigenfrequencies elsewhere in the complex plane, a two-dimensional amoeba searcher is unleashed on the landscape. Due to symmetries of the ODEs in their factors, the eigenfrequencies can be: purely imaginary (§3.1); purely real (§3.2); or occuring in complex conjugate pairs (§3.3).

If numerical integration traverses many orders of magnitude in then unfortunately there are instances in which one or more of the variables overflows or underflows the computer variables. In practice, numerical integration towards the outer boundary produces diverging , if is poorly chosen. This must be tackled through a practical numerical trick. One can pick the worst offending variable () and assign its magnitude to a new dependent variable, . The code tracks the logarithmic profile of this quantity, by numerically integrating an extra ODE simultaneously with the other ODEs:

(54) |

All the variables can then be expressed in a locally normalised form, . If each corresponding ODE () is abbreviated as a function then

(55) | |||||

The simplification depends on the linearity of the perturbations. Evolving (54) for in conjunction with the ODEs of the modified form (55) keeps all the variables normalised to moderate values that do not overflow the floating-point representations. Keeping the running values small also avoids driving the numerical integrators into stiff or unreliable domains. Near the outer boundary, where tends to become large, the running values of and can be kept on the order of unity or smaller. This prevents numerical noise in the sign of a score that is based on (53). For presentational purposes, the profiles can be rescaled by afterwards in post-processing.

For brevity and computational ease, we choose physical units such that , along with the stellar half-light radius and the stellar density there . Results can be scaled into metric or astronomical units by inserting fiducial values of these constants for a known galaxy.

Families of models will be labelled by their total stellar mass fractions, . This paper presents galaxy models for a variety of choices of . A fiducial “cosmic” baryon fraction is defined as in Saxton & Wu (2008) and Saxton & Ferreras (2010). For poorly understood reasons, observed densities of baryons (e.g. Persic & Salucci, 1992; Bell et al., 2003; Xue et al., 2008; Anderson & Bregman, 2010; McGaugh et al., 2010) and dark matter (Karachentsev, 2012) differ from cosmic mean estimates (except in rich galaxy clusters). There may be an unseen ambient sea of both materials, free from self-bound objects. Matching the abundances of observed galaxies to simulated CDM halo distributions implies baryon fractions depending on galaxy mass (e.g. Guo et al., 2010; Papastergis et al., 2012). This may mean that galaxies lost most of their primordial gas endowments, or it may indicate that real halo populations differ from CDM. Cosmological initial conditions are beyond the scope of the present analyses, so is for now a free parameter. Nonetheless §3.1 shows that spherical SIDM galaxies have stability limits that may influence the emergent distribution of .

In collisionless CDM, the outer halo density drops gradually without a well-defined boundary. These simulated objects are conventionally circumscribed by the “virial radius,” containing some multiple of the cosmic critical density, e.g. . If the natural density and radius scaling of SIDM haloes is not too different from CDM then we can estimate upper limits to the halo radii of represenative galaxies, . For the example of M87, ; for M32, ; for a red nugget, ; and for the Fornax dSph (data: Walker et al., 2009). It is however unclear whether this rule-of-thumb applies to SIDM, which warrants new cosmic structure formation calculations.

## 3 Results

### 3.1 stability properties

Classical arguments based on secular conditions suggested that polytropes should be unstable to runaway expansion or collapse (e.g. Ritter, 1878; Emden, 1907; Chandrasekhar, 1939, pp 51–53). This tendency manifests itself (in the present formulation) as a conjugate pair of modes with imaginary values. These “collapse modes” or “growing modes” involve no oscillatory motion, only exponential growth or decay. In practice such modes are found by searching along the imaginary axis of the Argand plane.

Next, a polytropic halo can be tested in a static stellar background (acting as a fixed potential, with ). The stability properties can be compared to the generic system where both the stars and dark matter are enabled to move radially. Figure 2 illustrates modes of a particular model with , comparing the results if the stars are fixed or free. With stars held fixed, a collapse mode appears. With the stars enabled to move, the collapse mode vanishes from the model shown. This stabilisation is typical of galaxies where the halo extends farther than .

The stabilisation can be explained qualitatively as follows. SIDM fluid carries pressure waves, and modes are collapsible because pressure may not rise strongly enough to counteract compressive perturbations. Stars however tend to persist on their orbits unless tidally disturbed, and do not conduct true sound waves. They deafen the propagation of overdensities, and exert a local damping and stabilising effect. For many galaxy models this suffices to confer global stability. However, as we shall see below, in extreme cases the stabilisation is insufficient (§3.1.1) or a new type of instability emerges (§3.1.2; §3.1.4).

The rest of this paper assumes that the dark matter and stars are gravitationally coupled and both respond to radial perturbations. Exhaustive numerical surveys have located the collapse modes as functions of: the equation of state parameter (), the radial extent () of the halo, and the global richness of stars (). Figure 3 illustrates the collapse frequencies as functions of and (colours), for models where dark matter dominates (small , in respective panels). Figure 4 compares the modes when DM is present in the “cosmic” abundance, and a case that is DM-poor. Realistic haloes could have radii of tens or hundreds of kpc, perhaps a range depending on the galaxy. Ultra-dense and ultra-diffuse cases are also plotted to clarify some theoretically notable modes. Two distinct types of modes occur in particular domains.

#### 3.1.1 slow growing modes

If is large (a very diffuse halo) and then an instability occurs. This is a reemergence of the instability of classic one-component polytropes with . In these cases the stars are unable to suppress instability because the stellar density is too low in the outer halo. This happens over a larger domain of if is lower (dark matter dominates stars). As is increased, the mode turns on suddenly: for below this threshold, a galaxy is safe; for slightly larger the mode operates on its shortest timescale; for very large the eigenfrequency diminishes (). In a survey of configurations with in the range , the frequency scale of these “slow” growth modes is lower than the effective stellar frequency scale, . In other words, the -folding duration of the growth / collapse would last several tens of orbits for most stars. The slow mode straddles the accoustic frequency scale: always for haloes of the largest radii considered, but more usually. The ratios are lower for than , and lowest for .

The perturbed mass flux density through a sphere of given radius is , so the quantity is a useful expression of the radial flows of stars and dark matter. Figure 5 depicts the radial flow structure during “slow mode” evolution. Amplitudes for stars tend to be highest in the centre and decline radially. The dark flow has its minimum in the centre. The stars tend to move in the same direction everywhere. Dark matter has alternating ingoing and outgoing layers, separated by nodes at rest. In the interior, dark matter moves the same way as the stars. Then there is a node at some radius around to (farther out for larger ). Outside this node, the dark flow is in antiphase with the stars. Some models have further nodes and flow reversals at larger radii in the outer halo (e.g. two far outer nodes in the example in Figure 5). The stability analyses are symmetric under time reversal, so it is impossible to say whether a particular sign implies collapse or expulsion of matter. The slow mode could imply contraction of all stars and the inner dark matter while an external dark layer escapes. It could also mean radial expansion of stars and inner halo, while dark outskirts converge inwards.

#### 3.1.2 fast growing modes

If is sufficiently small then one or two growing modes appear. They occur for all choices. The affected models have hot and dense haloes, with much of the dark mass residing inside the half-light radius. These “fast” modes occur at high frequencies: usually , (and often for larger cases). Where there are two “fast” eigenvalues, the uppermost is quick compared to the acoustic crossing of the halo (, in many models by a few dex). The lower eigenvalue is faster than if , but can be slower than for .

Figure 6 depicts the mass flows due to the lower- fast-mode for illustrative models with halo radius and DM enriched over “cosmic” composition. In some ways these profiles resemble the waveforms of slow modes: the stellar amplitude is high in the interior; the dark matter amplitude is peaked farther out. The dark matter has one node near the outer boundary, well outside the half-light radius of stars. Except in the outskirts, the dark mass is flowing in the same direction as the stars. For the corresponding higher- fast modes (Figure 7) there is again just one node in the dark matter, but both constituents show a dramatic rise in amplitude towards the outer boundary. The rise is steeper for larger and larger . These fast modes act mainly in the diffuse outskirts: rapidly driving the stars and dark matter in opposite directions. The core is essentially unaffected.

Since the fast instability occurs for all , its main causes must be traits of the stellar profile. We might suspect the stellar cusp and its coldness ( at the origin) however tests show fast modes afflicting non-cuspy models with a Plummer (1911) profile of . If the exponential stellar fringe is a factor, the fast modes might differ in tidally truncated models (generalising King, 1966) which will be analysed in a future paper. In the present model, the index of the stellar gravitational field () curves significantly near but settles to the asymptote farther out. If the halo surface occurs at then the oscillatory restoring force varies in form during the cycle. Fast instabilities may result from the inadequacy of restoring forces for haloes that truncate in this region.

#### 3.1.3 safe configurations

Fast modes and slow modes are absent throughout a special intermediate range of . In these cases the interaction between stars and the halo suppresses the collapse mode (e.g. in the “coupled” case of Figure 2). These haloes with are stabilised by the mingling of collisionless stars, contrary to naïve textbook expectations. The three radial domains — defined by the “fast” instability, stable zone, and “slow” instability — imply limits on the compactness of the stable endpoints of galaxy evolution. The fast growing modes mean that models are unstable unless exceeds some lower limit (which is a slowly varing function of and ). The slow mode only afflicts cases, and it applies an upper limit on . This upper limit is a sensitive function of . For models with that are poor in stars (), the slow- and fast-mode domains leave only a narrow band of safe . This safe band becomes wider as increases. A numerical survey indicates that the slow mode (upper limit) vanishes for sufficiently high . (It persists till for , for , and for .) From this point till , galaxies with arbitrarily wide outer radii are stable. Figure 8 illustrates the stability limits on the central dark density () as a function of . Figure 9 shows the corresponding limits on the dark mass fraction within the half-light radius (). The latter quantity can be measured through gravitational lensing and measurements of kinematic tracers. Significantly, it appears that a SIDM halo around an isolated de Vaucouleurs type galaxy is only stable if the dark mass within is less than of the total, or less depending on and .

Mass modelling of observed elliptical galaxies often indicates values in the range of a few tens of percent (Loewenstein & White, 1999; Ferreras, Saha & Williams, 2005; Thomas et al., 2005, 2007; Cappellari et al., 2006; Bolton et al., 2008; Tortora et al., 2009; Weijmans et al., 2009; Saxton & Ferreras, 2010; Grillo, 2010; Memola, Salucci & Babić, 2011; Murphy, Gebhardt & Adams, 2011; Norris et al., 2012; Tortora et al., 2012; Wegner et al., 2012). If these galaxies are stable within polytropic haloes, we could broadly infer independently of . Observational estimates of the stellar mass are consistent with this range. No large nearby galaxy is yet known to challenge the density limit of the fast instability. A few exceptional galaxies appear DM-poor at large radii (Méndez et al., 2001; Romanowsky et al., 2003; Salinas et al., 2012). If their impoverishment is not an illusion due to orbital anisotropies, then we might explain them in terms of runaway halo loss via either of the instabilities (§4.4).

#### 3.1.4 fast modes in DM-poor systems

Entering the domain where stellar mass is globally dominant (Figure 4) a new fast mode appears. Beyond some threshold () an island of the fast modes begins to affect large- haloes with . This implies that DM-poor systems have a maximum stable radius, just like the DM-rich systems that are subject to slow modes. The values of these fast modes are slowly varying functions of (horizontal curves in bottom panel, Figure 4) which implies they arise from the stars rather than the DM crossing time.

As is increased further, the fast mode affects haloes with all possible values, and then there are no stable configurations with . For example when there are no stable models with or . The fast modes that afflict these DM-poor models are found to have steeply rising waveforms (mainly acting at the halo surface) resembling those in Figure 7. Results were not sought beyond , because the system of equations becomes stiff in ways that make the root-finding slow and arduous.

### 3.2 neutrally stable oscillations

As in single-fluid asteroseismology, (e.g. Eddington, 1918; Cox, 1980; Pijpers, 2006; Aerts, Christensen-Dalsgaard & Kurtz, 2010) each polytropic halo exhibits an infinite sequence of radial pulsation modes at discrete real values of . These radial modes are neutrally stable: neither growing nor decaying. This is to be expected, as the present analysis is purely adiabatic: omitting any dissipation or mixing processes (e.g. shocks in the dark matter; violent relaxation of the stars; chaotic migration of stars).

#### 3.2.1 flow profiles of standing waves

Figure 10 illustrates examples of eigenfunction waveforms expressed in terms of oscillating mass flux density, () for the dark matter (coloured curves) and stars (black curves). For the stars, flow amplitudes are higher in central regions and decrease with radius. For the dark matter, the peak amplitudes occur at intermediate radii. Such profiles are generic to the radial pulsations of all models with the values inspected so far. Near the centre, the dark matter and stars oscillate in phase, but in outer regions they can oscillate in anti-phase.

For each constituent there are nodes at specific radii, where the matter remains locally at rest, and where there is a local reversal in the flow direction. The number of nodes of the dark matter oscillations is a harmonic number, . Values of uniquely label the eigenmodes. As increments, additional nodes appear in the outer fringe of the halo (e.g. between and for the modes in Figure 10). Most of the nodes appear far enough out that the keplerian orbital frequency is locally less than the pulsation frequency: (marked with arrows in Figure 10). If wave amplitudes grew to the nonlinear regime, some tidal shocking and orbital migration of stars could occur. Such dissipative effects scale as which is of non-linear order , while the radial displacements are of order . Inwards/outwards pulsational displacements are an immediate effect, but orbital heating (and net outward expansion) could become significant gradually as more cycles elapse.

Depending on and the galaxy model, a number of nodes can also occur in the stellar profile (). For a galaxy with , and , the node numbers are and so on. Figure 11 shows the dark and stellar node positions for the first thirty oscillatory modes of this model, and for equivalent models where and . It appears that never decreases as increases, but there is no explicit link between these numbers. The sequence of and their radii provide a fingerprint of a particular galaxy, and must be obtained numerically.

Figure 11 also shows that the dark nodes are more crowded near the outer boundary () in halo models with larger . This is a consequence of slower wave propagation in the fringes, since small-cored (large ) models have proportionally lower profiles in the outskirts. Stellar nodes are more sublty affected by . For models with small , the stellar nodes are more numerous and reach smaller radii for given . Haloes with large tend to have fewer stellar nodes, and they are distributed farther out. Varying the stellar mass fraction has slight effects on the positions of the nodes in the outskirts. Across a wide domain (but with fixed , and ), the stellar nodes do not shift significantly. The outer nodes of the dark matter appear more affected: they occur at radii a few percent farther out if is very small.

For waves that have at least one DM node (), the radius of the innermost node varies more widely than any other node: from to in the domain of , and . Its location varies only sublty with harmonic number . Figure 12 depicts the detailed distribution of nodes in the inner, astronomically observable regions () of galaxies with and where the halo radius is and the DM content is rich or “cosmic” (). All else being equal, the innermost dark node tends to occur at smaller radii if the galaxy is DM-poor (large ), and farther out if the galaxy is DM-rich (small ). For DM-rich haloes, a larger radius shifts the inner dark node to larger radii. For some DM-poorer haloes, larger shifts the inner dark node inwards. Comparison of the cases in Figure 12 also shows how the lower- models have more stellar nodes at radii small enough that ripples might be observable in the profiles of starlight and stellar velocities.