Constraining self-interacting dark matter with scaling laws of observed halo surface densities

The observed surface densities of dark matter halos are known to follow a simple scaling law, ranging from dwarf galaxies to galaxy clusters, with a weak dependence on their virial mass. Here we point out that this can not only be used to provide a method to determine the standard relation between halo mass and concentration, but also to use large samples of objects in order to place constraints on dark matter self-interactions that can be more robust than constraints derived from individual objects. We demonstrate our method by considering a sample of about 50 objects distributed across the whole halo mass range, and by modelling the effect of self-interactions in a way similar to what has been previously done in the literature. Using additional input from simulations then results in a constraint on the self-interaction cross section per unit dark matter mass of about $\sigma/m_\chi\lesssim 0.3$ cm$^2$/g. We expect that these constraints can be significantly improved in the future, and made more robust, by i) an improved modelling of the effect of self-interactions, both theoretical and by comparison with simulations, ii) taking into account a larger sample of objects and iii) by reducing the currently still relatively large uncertainties that we conservatively assign to the surface densities of individual objects. The latter can be achieved in particular by using kinematic observations to directly constrain the average halo mass inside a given radius, rather than fitting the data to a pre-selected profile and then reconstruct the mass. For a velocity-independent cross-section, our current result is formally already somewhat smaller than the range $0.5-5$ cm$^2$/g that has been invoked to explain potential inconsistencies between small-scale observations and expectations in the standard collisionless cold dark matter paradigm.


Introduction
On cosmological scales, dark matter (DM) is about five times as prevalent as ordinary matter [1], and known to be the main driver of structure formation. The paradigm of cold, collisionless dark matter (CDM), one of the main ingredients of the cosmological concordance model, has been remarkably successful in describing the observed distribution and properties of structures in the universe [2,3]. In apparent contradiction to this success, however, current observations do not actually constrain the DM self-interaction cross section to be smaller than that of the strong interaction between nucleons (for a recent review, see [4]), which is many orders of magnitude less stringent than corresponding bounds on DM interacting with standard model particles [5,6]. Self-interacting DM (SIDM) thus remains a fascinating option which, if directly confirmed observationally, would significantly reduce the number of possible DM candidates from particle physics. Such observations would, furthermore, offer a window into the particle properties of DM that may be impossible to access by other means -a fact which has created significant attention in recent years (see, e.g. Refs. [7][8][9]).
The typical phenomenological handle on SIDM models, and in fact the context in which the very idea of SIDM was proposed in the first place [10], are observables related to structure formation at galactic scales and below. In particular, it has been demonstrated [11][12][13][14][15][16][17][18] that SIDM could alleviate all of the potential small-scale problems of ΛCDM cosmology [19], most notably the "core-cusp" [20,21], "too-big-to-fail" [22,23], "diversity" [24,25] and (for late kinetic decoupling) "missing satellites" problems [26][27][28]. These solutions require a DM selfscattering cross section per unit DM mass of the order of σ/m χ ∼ 1 cm 2 /g, close to current exclusion limits. But even if all current small-scale discrepancies between ΛCDM observations and expectations are resolved, in the sense that they can be attributed to baryonic effects, the uncertainty in Milky Way mass, cosmic variance and observational uncertainties (see e.g. [25,[29][30][31]), SIDM remains an intriguing possibility -not the least given the absence of any undisputed positive results in searches for DM particles. Current constraints derive from a large number of observations on different scales, see again [4] for an overview, but their common feature is that they are typically obtained from individual objects. Given the modelling uncertainties of the respective objects, ranging from their individual formation history to their baryonic content and its potential interference with the effects of SIDM, this is not unproblematic.
Here, we introduce a new way to constrain DM self-interactions that instead relies on ensembles of many astrophysical objects, thereby reducing the systematic uncertainties related to individual objects. Concretely, we revisit the well-known observed scaling relation between surface density and halo mass [32][33][34][35], which is essentially understood in ΛCDM cosmology as a reflection of a similar relation between halo mass and concentration. We investigate how this relation is affected by the central cores observed in the DM distributions in various (dwarf) galaxies and (possibly) also galaxy clusters. Assuming that these cores can exclusively be explained in terms of SIDM, we derive an experimentally robust upper bound on the effect of DM self-interactions once we take into account additional constraints deriving from a direct comparison to SIDM simulations. The final translation of this bound to a constraint on the physical self-interaction cross section per unit mass, σ/m χ , is necessarily somewhat less robust as it involves less certain theoretical modelling of the effect of DM self-interactions. Even when taking this into account, we arrive at an upper bound on σ/m χ that is competitive compared to existing bounds in the literature. We note again that this is mainly the result of the relatively large number of objects that we include in the analysis, combined with crucial input from simulations about quantities that cannot be directly constrained observationally.
This article is organized as follows. We start in Section 2 by reflecting on how DM self-interactions would change the halo density profiles expected in ΛCDM cosmology. In Section 3, we then introduce the observed surface density scaling relation, as well as its ΛCDM explanation, and derive how the surface density is expected to change as a function of the DM self-interaction cross section. We finally derive constraints on the latter in Section 4.1, and discuss them, before concluding in Section 5. In two Appendices, we describe in more detail the halo objects that we include in our analysis (App. A) and briefly review how the halo age depends on its mass in ΛCDM cosmology (App. B).

Halo profiles for self-interacting dark matter
Numerical simulations of gravitational clustering in ΛCDM cosmology reveal that collisionless DM halos have a universal density profile that, at all redshifts and masses, roughly follows the form suggested by Navarro, Frenk and White (NFW) [36,37]: Here, the scale radius r s marks the position where the profile has a slope of d(log ρ)/d(log r) = −(3 + γ)/2, and hence the transition between the slopes encountered in the inner and the outer part of the halo. At redshift zero, CDM halos have mostly converged to a negative inner slope of γ ≈ 1, with some scatter, and this is what we will use in the following when referring to the NFW profile. More recent simulations sometimes tend to prefer an Einasto profile [38], which is slightly shallower in the very central parts of the halo, but this difference will not affect our discussion.
It should be stressed that the above result only holds for DM-only simulations. Including baryons can both lead to significantly steeper, even cuspier DM profiles -mostly due to a process known as adiabatic contraction [39][40][41] -and to the formation of much shallower profiles in the inner parts, often referred to as cores, due to feedback from star formation and supernovae [42]. Despite the enormous numerical challenges, there has recently been significant progress in including baryonic effects in hydrodynamical simulations of structure formation even at cosmological scales [43][44][45][46]. However, this comes at the price of implementing phenomenological prescriptions, rather than prescriptions based on first principles, that have to be tuned to match one class of observations in order to successfully "predict" another. In this sense, a complete understanding of all relevant scales and processes is still missing. We will try to avoid these difficulties by mostly focussing on systems where the effect of baryons is subdominant (but get back to this discussion in Section 4.3).
The main effect of DM self-interactions, as we will see, is to isotropize the DM phasespace distribution f . In other words, the claim is not only a Maxwellian velocity distribution, but that the velocity dispersion σ v is (approximately) constant at least in some region of the halo, namely in its central parts where the DM density and hence the collision rate is highest. This behaviour has been confirmed by numerical simulations [12,16,47] for cross sections in the ballpark that we are interested in, σ/m χ ∼ 1 cm 2 /g, showing results broadly in agreement with the expectations originally formulated along with the SIDM proposal [10]. It is worth noting that a constant σ v is only possible in the weakly interacting regime for the DM particles. The entropy increase in the inner parts then leads to a formation of an inner core of roughly constant density -again confirmed by numerical simulations. For much larger cross sections, on the other hand, DM would behave as a collisional gas, leading to an isothermal density profile with ρ ∝ r −2 (see, e.g., [48,49]). We can easily check that we are indeed in the weakly interacting regime, motivating the claim of a constant σ v , by looking at the mean free path λ of the DM particles, This is clearly larger than the sub-kpc cores reported in dwarf galaxies, for realistic core densities ρ core = O(1) M /pc 3 for dwarf galaxies [50] and ρ core = O(0.1) M /pc 3 for clusters [51]. For these cross sections, a DM particle thus typically only scatters at most a few times during the whole halo lifetime t age even though it may pass through the core region much more often. In the outer parts of the halo, on the other hand, scattering is so rare that the standard (NFW) DM profile should be unaffected. When estimating the region of equilibrium r < r SIDM , for which we have σ v ∼ const., we should thus refer to the average density inside this region, ρ χ SIDM , to compute the mean free path. We can do so by demanding that the average time between collisions should roughly equal the halo time, i.e. 1 Here, v χ = 4/ √ π σ v is the average relative velocity of the DM particles inside the core, and we have allowed for an unknown factor ξ 1 to account for the fact that, for scattering between two non-relativistic particles, the relaxation time (the time-scale needed to achieve thermal equilibrium) is indeed somewhat larger than the scattering rate (while for scattering between a DM particle and a relativistic particle of energy E m χ , the difference would be ξ ∼ m χ /E and hence much larger [52]).
In order to relate v χ to the circular velocity v c let us consider the Jeans equation, where Φ is the gravitational potential and it has been used that SIDM halos should be isotropic and spherically symmetric in their inner parts [10,12]. For the range of radii where σ v is constant, r < r SIDM , we can rewrite the Jeans equation in the simple form is the logarithmic slope of the DM profile at radius r, and M (r) denotes the halo mass inside r. This leads, finally, to the following implicit definition of r SIDM : (2.8) Here, α SIDM is the (negative) logarithmic slope of the SIDM profile at r = r SIDM ; as we will see below, its value falls into the range 1 α SIDM ≤ 2. For t age in the above expression, we will use the average halo age as a function of the virial mass adopted from Refs. [53][54][55], see Appendix B for further details.
Before continuing, let us comment on why it is so challenging to directly relate the scale r SIDM to the observed core size, r core ≤ r SIDM , which would provide a handle on the scattering cross section that is observationally easy to access. It has been argued [8,9,17,56,57] that this question can be fully resolved by considering the Jeans equation. In particular, when assuming σ v (r) ∼ const. and using Poisson's equation, as well as neglecting the contribution of baryons to the gravitational potential, Eq. (2.5) simplifies to It can be easily checked that the isothermal sphere, ρ(r) = σ 2 v /(2πGr 2 ), is a solution to this equation; as discussed above, this is the physical solution for cross sections much larger than what we are interested in here. There is, however, also another class of solutions, with ρ χ (0) = 0, which describe a roughly constant density at small radii, and a profile approaching the isothermal sphere solution for large radii. The phenomenological modified (or pseudo-) isothermal sphere,  Figure 1. The mass and density profiles of the dwarf spheroidal galaxy Carina. The shaded regions are reconstructed values from kinematic data measured within 1 kpc radius for this object [59]. The green, blue and red lines are examples of SIDM profiles adopting the method from Ref. [8], with the slight modification given in Eq. (2.4). for an increasing value of the self-interaction cross section as indicated. The dots represent the radius r SIDM in each of these cases, and the vertical dashed line is the half-light radius.
is a very good approximation to this class of solutions, at least for r r 0 and r r 0 . Further phenomenological profiles that are commonly used to describe cored density distributions, and fit observational data, are the cored NFW and Burkert [58] profiles, see Appendix A, which reproduce the expected scaling of the NFW profile rather than that of the isothermal sphere at large radii -where in fact we cannot expect a constant σ v anymore.
The cored solutions to Eq. (2.9) still have one free parameter, for a given σ v , which can equivalently be chosen as the central density or core size (for the Burkert parametrization, e.g., the core size is given by r 0 = σ v / √ 2πGρ 0 ). The issue is that the relation between the core size (or core density) and r SIDM very strongly depends on the assumed matching conditions to the asymptotic NFW profile -typically taken to be that self-interactions do not change the mass inside any radius larger than r SIDM , and that the density profile remains continuous at r = r SIDM -and hence on the part of the profile that should not be affected by DM self-interactions. We illustrate this in Fig. 1 for the case of the dwarf spheroidal galaxy Carina, where the grey area indicates the range of the mass profile consistent with dynamic observations of the stellar population [59]. We note that this is a somewhat extreme example, picked for the sake of our argument, where the whole range of kinematic observations is compatible with a cored density profile. Choosing one central density and core size consistent with these data, we plot the resulting profile for a large range of self-interaction cross sections σ/m χ . Clearly, dynamical observations of this object hardly constrain r SIDM , allowing r SIDM r c even if the core size r c is well measured (we obtain qualitatively the same result if we replace in Eq. (2.4) ρ χ → ρ χ , as in Ref. [8]). A much more stringent constraint on the scale of self-interactions instead results from requiring that the asymptotic behaviour of the NFW profile matches the average expectation of ΛCDM cosmology [8,9]. As already indicated, classical dwarf spheroidal galaxies like Carina are the most extreme objects in this respect. For other objects -like the dwarf galaxies considered in [8], which are not Milky Way satellites -the outer (NFW) part of the profile is typically much better constrained even without taking into account cosmological priors. Still, there is a remaining worry that the inferred (bound on the) cross section depends to some extent on the choice of matching conditions with the asymptotic NFW profile rather than only on the (not directly observable) physical size of the self-interaction region r r SIDM .
In conclusion, there is a considerable theory uncertainty in how to relate, from first principles, the observationally accessible core size of individual halos (as opposed to the theoretically well-defined, but observationally inaccessible scale of self-interactions) to the self-interaction cross section. Concretely, the method typically adopted so far may allow larger ratios between r SIDM and core size than what one would expect from the underlying physics. In fact, it is quite conceivable that the exact core size even depends to some degree on the formation history and hence on environmental effects. In this article we will instead fix a related quantity, namely the ratio between average densities inside the core and r SIDM , by directly comparing it to results taken from simulations. Technically speaking, from the point of view of the Jeans analysis, this amounts to adopting a different set of boundary conditions to solve Eq. (2.5). Furthermore, we will make use of observed scaling relations in ensembles of astrophysical objects which are more robust to astrophysical uncertainties than the analysis of individual objects.

Observations
Let us define the mean surface density, sometimes also referred to as Newtonian acceleration, of a halo as is the mass enclosed inside the radius r (we assume that this mass is dominated by the DM component). For a cored profile, this quantity is maximized at a radius close to the core radius. 2 For the NFW profile Σ(r) is approximately constant for r r s and reaches its maximum for r → 0, with Σ(0) = 1.15 Σ(0.1 r s ) = 2.62 Σ(r s ). We will denote this maximal value as Σ max . Refs. [32,33] argued that for galaxies (from dSph to ellipticals) Σ max is a constant, independent of the galaxy type. However, as Refs. [34,35] have demonstrated this conclusion has been based on a too small range of dynamical masses. When taking into account all type of DM halos, from dSph galaxies to galaxy clusters, one can see that Σ max (or quantities that can be related to it) increases with M 200 . As demonstrated in [35] this scaling can be explained within the secondary infall model. This result has later been confirmed with larger datasets (see e.g. [60][61][62]) with a scaling relation given by [62] when fitting a single power-law to systems (almost) exclusively composed of DM. Let us stress that we introduced here the maximal surface density only for the sake of simplicity. In practice, one would instead choose a radius which is small but where the enclosed mass and hence the surface density is observationally still well constrained. As long as this radius is directly proportional to the core radius, or the scale radius in the case of an NFW profile, this does not affect the scaling with M 200 (but can significantly reduce the observational scatter in this relation [34,35]).  For the purpose of constraining DM self-interactions as explained further down, we selected DM dominated objects over a large range of halo masses where we could find reasonable fits to both a cored profile and an NFW profile in the literature (for more details, see Appendix A). For an NFW halo, we have In Fig. 2, we show the surface densities of the objects in our sample, inferred from the fit to the NFW profile and evaluated at the scale radius. The errors of the data points we calculate by using uncorrelated 1σ errors on the NFW parameters quoted in the corresponding references (using instead directly the kinematic constraints on the parameter combination ρ s r s , as available for dwarf galaxies, would result in much smaller errors). We then fit a power-law to these data points, resulting in which is consistent with the scaling reported earlier. We indicate the best-fit power-law as a black dashed line in Fig. 2.

ΛCDM interpretation
To proceed, it is useful to introduce the halo concentration, where r 200 is the virial radius defined as the radius inside which the mean density ρ exceeds the critical density ρ c = 3H 2 0 /(8πG) by a factor of 200 (for other profiles than NFW, suitable generalizations of this definition of c exist [53]), i.e.
which allows us to exchange the parameters (ρ s , r s ) for (M 200 , c). For the surface density inside the scale radius, we thus find (3.11) From numerical simulations [63][64][65][66][67][68][69], but also observations [70][71][72][73][74][75][76][77][78], the concentration is not independent of the halo mass, but rather follows a simple scaling law. There is, however, a significant object-to-object scatter associated to this relation, and even the best-fit values for slope and normalization of this power-law differ in the literature. One of the more often used results is the one by Macciò et al. [65] (3.12) Using the mean value of this slope in Eq. (3.11) results in the slope d(log Σ max )/d(log M 200 ) ranging from roughly 0.13 at low masses to 0.25 at very large masses. In ΛCDM cosmology, the scaling of the surface density as given in Eq. (3.3) can thus consistently be interpreted as a reflection of the above concentration-mass relation [34,[60][61][62]. We note that, from Eqs. (3.11,3.12), we can infer not only the scaling of Σ with M 200 , but also the normalization of that relation in ΛCDM cosmology. In fact, we can turn the argument around, and provide an updated measurement of the concentration-mass relation by fitting Eq. (3.11) to the data points shown in Fig. 2. Assuming again a simple power-law for c(M 200 ), we find (3.13) In Fig. 2, we show the resulting surface density as solid red line, with the shaded red region indicating the uncertainty in the concentration-mass relation that we derived.

Cored halo profiles
Let us now discuss how the scaling of the mean surface density would change in the presence of cored profiles as produced by SIDM. For small core radii, much smaller than the scale radius r s of the NFW profile expected in the outer parts, the enclosed mass inside r s will obviously not be significantly affected. We are hence no longer interested in evaluating the surface density at r s , as before, but at a smaller radius. The optimal choice in this respect is the core radius itself. 3 While the surface density of NFW and cored profiles would differ even more at radii smaller then the core radius, in particular, the surface density in this regime is much less well constrained by observational data (and would anyway decrease with respect to its value at the core radius). We use this opportunity to stress again that the scale of efficient self-interactions, r SIDM , is essentially impossible to determine observationally, so we are also not interested in considering the surface density at this (somewhat larger) radius. For any cored profile parameterization, we could now in principle follow an approach in analogy to what we did for the NFW case, and analytically express Σ(r 0 ) in terms of the halo mass and the core scale parameter r 0 . This, however, is impractical because it would have to be done for each choice of profile parameterization separately. Besides, one would need independent prescriptions of how to relate r 0 to the scale of efficient self-interactions, r SIDM , for each of these cases, which complicates the analysis we are interested in here. In the spirit of keeping the discussion as model-independent and general as possible, we therefore consider instead the experimental core radius r c , which we define as the radius at which the best-fit NFW profile equals the central density of the best-fit cored profile, (3.14) Unlike the profile parameter r 0 that appears, e.g., in the Burkert profile and the modified isothermal sphere, the core radius r c defined in this way is relatively independent of which cored profile is used for the fit, and typically more robustly constrained observationally. The "observed" surface density at this core radius is then simply given by where we note that ρ cored c ≈ ρ cored (0) = ρ NFW (r c ) because the DM density inside r c is almost constant (we do not use this last approximation in our analysis). We plot Σ c, obs in Fig. 3, for the same objects that we used in Fig. 2, as a function of both virial mass and the circular velocity at the core radius. Concretely, we assumed that the  parameters of the NFW and cored profile follow a Gaussian distribution in each case, with standard deviation as quoted in Appendix A. Drawing parameters from this distribution, we then determined ρ c and r c for a large sample of profile realizations for each halo, to calculate v c (r c ) and Σ c, obs . The data points in Fig. 3 thus created provide the surface density corresponding to the best-fit profile parameters, with errors corresponding to one standard deviation in our sample. We note that this Monte Carlo approach results in conservative estimates for the errors ∆Σ c, obs , because it treats the two profile fits to identical halos as being independent. This will result in conservative limits when eventually using this dataset to constrain the DM self-interaction rate.
Let us now consider the theoretically expected surface density at the core radius, assuming that the existence of the core is the result of DM self-interactions. We start by recalling that for r > r SIDM we expect the standard NFW profile to be unaffected by DM self-interactions. This implies that the mass, and hence the average (surface) densities inside r SIDM remain the same: We can then use M NFW from Eq. (3.4) and the implicit definition of r SIDM given in Eq. (2.8) to relate the self-interaction scale to the scale radius, for a given NFW profile and self-interaction cross section: where η is a solution to the equation For illustration, we show in Fig. 4 how the resulting r SIDM the average density ρ SIDM (from Eq. (2.8)) scale with the self-interaction strength. In order to produce these curves, we take ξ = 1, α SIDM = 2 and implement an average halo ago t age following the prescription in Appendix B. The range of concentrations for each mass displayed in this figure very roughly corresponds to that given by the c − M 200 relation, see Eq. (3.13). For the displayed halo mass of 10 7 M , we allow for a larger scatter taking into account that dwarf satellites are not too well fitted by this relation, see Fig. 2 and the discussion in Refs. [34,35]. Once we have ρ SIDM , we obtain the theoretically expected surface density as In the last step we have introduced a phenomenological ratio κ of the average densities inside r SIDM and r c , respectively: This ratio, as argued in Section 2, is difficult to determine observationally or directly from first principles. On the other hand it is a quantity that turns out to be tightly constrained by simulations, and this is what we will make use of when determining limits on σ/m χ in Section 4. In the following, we will simply treat κ as a constant, i.e. independent of cross section and halo mass, which is consistent with the simulation data we have explicitly looked at. We note that the exact form of κ, and hence the robustness of our limits, can be improved by taking into account a larger sample of (new) simulations, which is beyond the scope of the present work.
Let us stress that ρ c here refers to the average DM density inside the same observationally determined core radius r c as defined in Eq. (3.14). If baryonic effects also contribute to the observed core, then the net effect of an increasing core size r c and a decreasing core density ρ c is a surface density Σ c, obs that is necessarily smaller than the theory expectation for a halo consisting exclusively of SIDM as stated in Eq. (3.19). This will allow us to place upper, but not lower limits on the self-interaction rate. We also note that the above definition of κ -via ρ c ≈ ρ(0) = ρ NFW (r c ) -implicitly fixes the ratio of r c and r SIDM (as a function of σ/m χ , κ and the NFW profile parameters). Using κ from Eq. (3.20), we can thus also calculate the logarithmic slope α SIDM of the density profile at r SIDM that appears in Eqs. (2.8) and (3.18). To do so, we numerically solve the Jeans equation (2.9) with ρ (0) = 0 and determine ρ χ at r = r SIDM , as a function of κ. As explained in Section 2, the cored solution to the Jeans equation with fixed central density ρ χ (0) depends only on one dimensionless parameter, r SIDM Gρ χ (0)/σ v ≡ r SIDM /r Jeans , so this mapping between α SIDM and κ must be unique.
Our formula (3.19) describes, as physically expected in the weakly interacting regime, a surface density that decreases with increasing cross section. This implies that it must break down once we leave this regime, and the core size instead should start to decrease again as the profile approaches the isothermal r −2 solution. We will not consider such large cross sections in our analysis. In the opposite limit of σ/m χ → 0, on the other hand, we recover the NFW expectation for the surface density at any given radius r c . This is an important property of our model when constraining the self-interaction strength σ/m χ in the next section.

Statistical treatment
In the previous section, we have discussed how DM self-interactions can affect the observed scaling of the surface density with halo mass. Here, we will derive constraints on the DM self-scattering cross section based on this prescription. As experimental input for our analysis, we consider the surface density at the experimental core radius r c for the objects described in Appendix A, as shown in Fig. 3. In order to determine limits on the "signal" strength, we then use the likelihood ratio test [80]. Here, we have introduced an effective cross sectionσ to reflect the degeneracy of the physical scattering cross section σ with the O(1) factor ξ as introduced in our basic SIDM ansatz of Eq. (2.4). We model the total likelihood as a product of normal distributions over each halo object ("data point") i, where f i is the value of the (logarithmic) surface density inferred from the kinematical anlaysis, σ i its variance, and µ i is the (logarithmic) surface density predicted by the model. The latter depends not only on the self-interaction strength S, but in principle also on a number of nuisance parameters {α k }. The contribution from DM self-interactions therefore enters with a single, non-negative degree of freedom, which implies that 95% CL upper limits on S are derived by increasing S from its best-fit value until −2 ln L has changed by 2.71, while refitting ("profiling over") all nuisance parameters {α k }. Concretely, we use Eq. (3.15) for the data points f i and Eq. (3.19) for the model prediction µ i . For the variance, we add observational and theory uncertainties in quadrature, σ 2 i = σ 2 i,obs +σ 2 i,theo . The errors in the "observed" surface density uncertainty, σ i,obs = ∆Σ c, obs , are determined as described in the previous Section and indicated in Fig. 3. In the "theory error" in this figure, we include two contributions: Here, the largest contribution, σ i,halo , is related to uncertainties in the observational determination of the halo parameters and picks up two contributions: i) from the variance in the core radius r c , which is determined in the same way as for σ i,obs and directly enters in the model prediction µ i via Eq. (3.19), and ii) from error propagation of the NFW halo parameters (ρ s and r s ) that enter in Eq. (3.18) for the average density ρ SIDM . For the halo age, we adopt a value of σ tage that corresponds to a rather generous factor of 1.2 [81] in the expected halo-tohalo scatter of the implemented t age (M 200 ) relation described in Appendix B. The value of κ in Eq. (3.19), finally, we extract from direct comparison with SIDM simulations [47], 4 finding κ = 3.9 with a scatter of σ 2 κ ≡ κ 2 − κ 2 = 1.4 2 . We vary κ freely within this range.

Results
In order to illustrate our approach, we compare in Fig. 5 the observed surface density (as shown in Fig. 3) with the SIDM predictions that we derived above, for various values of the effective self-interaction cross section. These plots clearly demonstrates that a too large self-interaction cross section would be inconsistent with the experimental data. In fact, we even see a slight preference for a non-zero value ofσ/m χ . Of course, the latter cannot be taken as an indication for SIDM, but is rather a reflection of the fact that we consider objects consistent with a cored profile -which leads to a smaller surface density than what would be expected for the NFW case.  We make these observations more quantitative by showing in Fig. 6 the full likelihood described in Section 4.1 as a function of S =σ/m χ . From this, we can read off an upper bound ofσ /m χ 0.12 cm 2 /g, (4.4) which roughly corresponds to a 95% C.L. limit as we have allowed κ to vary up to its maximal range within 2σ (allowing values up to κ ≤ 8.1, i.e. within the 3σ range of κ, the limit would relax toσ/m χ 0.14 cm 2 /g). We emphasize again that our method does not allow to put a lower bound on the cross section because there are also baryonic feedback processes, not modelled here, which could lead to a core and hence a reduced surface density. In order to understand how the above constraint depends on the type of objects that we include in our analysis, we include for comparison also the likelihood resulting from various subsets of our full halo sample. We will discuss this in more detail in Section 4.3 below. . The values of ξ obtained by applying our method to the simulated halos from Refs. [47] (green points) and [16] (red points). The black line is the best fit to the points.
Lastly, we would like to estimate the systematic bias of our model ansatz, i.e. the difference between the effective cross sectionσ and the physical cross section σ that we introduced in terms of the parameter ξ. To do so, we apply the same analysis as above, but now to 8 simulated halos with cross section 1 cm 2 /g in the mass range from 5 · 10 9 M to 2 · 10 14 M from Refs. [16,47]. In Fig. 7, we show the ratio of the reconstructed value of the cross section to the "true" cross section (i.e. the one used in the simulations). For this ratio, we find a best-fit value of ξ = (σ/m χ ) sim σ/m χ = 1.86 ± 0.32 . Ideally, we would of course use more simulated halos, for cross sections closer to the value of roughly 0.1 cm 2 /g that corresponds to our constraint, but such simulation results are currently not publicly available. The real theoretical uncertainty encoded in the factor ξ (as well as κ) may thus be somewhat larger, but a full investigation of this effect is beyond the scope of this work. Taking this caveat into account, we arrive at a limit of approximately σ/m χ 0.3 cm 2 /g (4.6) for the physical self-interaction cross section.

Discussion
In the previous section we have derived an upper bound on the effective self-scattering cross section, and stated the result in Eq. (4.4). Let us stress that this bound is not exclusively driven by the kinematic data, which allowed us to construct the surface densities shown in Fig. 3, but also by what we assume about the presently unmodelled relation between the average densities inside the core radius and the radius of efficient self-interactions, as parameterized by κ defined in Eq. (3.20). For the range of cross sections that we are interested in here, we find in fact that the bound onσ/m χ scales roughly linearly with the maximal value of κ allowed in our analysis. While data from simulations already strongly constrain the (average) value of κ to be much larger than what we have adopted in our analysis, more simulation data is thus needed in order to make this bound more robust. We note that the translation to a bound on the physical cross section adds, as stated in arriving at Eq. (4.6), another uncertainty to the present analysis -though one should stress that this uncertainty is shared by most analyses constraining SIDM (which is rarely spelled out explicitly). In other words, this does not affect the conclusion that our method to compare theory with observations is more robust than those based on individual objects. In any case, we clearly expect this to be an O(1) effect, at most, something which we have checked explicitly by comparison to a limited set of simulations (see Fig. 7 and, for a similar test, Ref. [8]), and which we think is captured in the approximate limit stated in Eq. (4.4). Still, we caution that in order to make the limit on the physical cross section more robust would require a further refinement of the analytical model, i.e. the formula (2.4), which has been used in a very similar way previously in the literature. This, in turn, calls for a more systematic study of the properties of simulated halos, which is beyond the scope of this article.
Another potential worry might be that our bound is mostly driven by a small subset of objects -which one then could argue may suffer from large systematic uncertainties, similar to bounds derived from individual objects. The dashed and dotted lines in Fig. 6 essentially demonstrate that this is not the case: both when excluding clusters (dashed line) and when excluding (all other) objects that are baryon-dominated in their central parts (dotted line) from our analysis, the limit does not change by more than what can be explained by the reduced sample size. Let us point out that these two model classes are indeed the main suspects when looking for such an effect: • Observations of colliding clusters lead to the strongest currently existing constraints on SIDM, with σ/m χ 0.47 cm 2 /g, and it has been argued that their small cores (if any) lead to even stronger bounds [8]. While not undisputed, this has triggered much phenomenological interest in velocity-dependent self-interactions in order to evade cluster bounds and at the same time allow for σ/m χ ∼ 1 cm 2 /g at (dwarf) galaxy scales, where the typical DM velocities are up to one order of magnitude smaller [82][83][84][85] (and more recently [8,11,15]). Our (cluster) bounds, hence, do not show a significant velocity dependence.
• Large baryon densities in the inner halo parts can lead to contracted profiles, instead of cores, counteracting the effect of core-formation due to self-interactions [17]. Such an effect could "hide" large SIDM cross sections, and hence potentially spoil our claim of deriving upper bounds on SIDM irrespective of the role of baryons. While we on purpose did not include any objects that are close to being as baryon-dominated as the corresponding examples in Ref. [17], a remaining worry might be that a similar effect could be seen in the small number of objects in our sample where the impact of baryons on the gravitational potential in the inner parts of the halo is similar to that of DM. As the dotted line shows, this is not a concern as removing those objects from our analysis does not significantly weaken the bound onσ/m χ .
We also checked, independently, the constraining power of dSphs alone (dash-dotted line in Fig. 6). This much weaker bound (σ/m χ 1.4 cm 2 /g) is mostly driven by Sculptor and Draco, leading respectively toσ/m χ ≤ 3.5 cm 2 /g and σ/m χ ≤ 4.4 cm 2 /g, while Carina alone formally allows for a cross section ofσ/m χ ∼ 50 cm 2 /g. While this appears broadly consistent with what was found in Ref. [9], we stress that one cannot easily compare these results as we do not impose a cosmological prior on the c − M 200 relation (or the distribution of circular velocities). From the discussion in Section 2, in fact, we would expect that kinematic data from classical dSphs alone would lead to even much weaker constraints (see for example Fig. 1). The solution to this apparent paradox is that, as already stressed above, also in our analysis it is not the kinematic data alone that set the constraints. For example, allowing κ to be as large as 20 -which is far larger than anything found in simulations -would imply that our dSphs bound relaxes fromσ/m χ 1.4 cm 2 /g toσ/m χ 4.2 cm 2 /g. In general, we find again that the bound scales linearly with the maximally allowed value for κ.
In this article, we have presented a new method to test SIDM and, as a proof of concept, derived stringent constraints already from a relatively small sample of objects. Prospects to derive more stringent limits, simply by increasing the sample size, are thus obviously promising. This would be particularly interesting for dwarf-scale field halos, for which currently no fits to cored profiles exist An even more promising way to significantly improve our limits is to reduce the errors on the observational parameters that enter our analysis, i.e. the experimental core radius and the average density inside this radius. This can be achieved if (the product of) these quantities is directly constrained in the kinematic analysis, rather than taking the detour via first fitting a density profile to the data. Recalling that the surface density is also very useful for understanding basic scaling laws of ΛCDM cosmology, viz. the concentration-mass relation, we use the opportunity for a general "plea" to observers: for many applications, it is more advantageous to present measurements of the average (surface) density of DM halos than fits to given profile parameterizations.

Conclusions
Self-interacting dark matter (SIDM) has been the subject of increasing interest in the last few years, both because it may provide a solution to the long-standing small-scale problems of the cosmological concordance model, and because it opens up interesting avenues for modelbuilding that involve DM particles which could not be detected by traditional means. In this article we have introduced a new method to constrain such DM self-interactions by re-visiting the surface density of astronomical objects ranging from dwarf galaxies to galaxy clusters. The main advantage of our method is that it is based on ensembles of objects, implying that the resulting constraints are rather robust and less sensitive to the often poorly understood astrophysical properties of individual objects. The other crucial input to our analysis is how the average density inside the region of efficient self-interactions is related to the core density; we obtain this ratio from a direct comparison to simulations.
We illustrated our method by selecting a sample of around 50 objects, as described in App. A, where both cored and NFW profiles have been fitted to the available kinematic data. We inferred the surface density at the experimental core radius for these objects, and compared it to the surface density expected for an SIDM scenario with given cross section, cf. Eq. (3.19) and Fig. 5. This allowed us to construct a total likelihood, containing all halo objects, as a function of the self-scattering cross section (Section 4.1) and derive an upper limit of approximately ∼ 0.3 cm 2 /g on the cross section per unit mass, σ/m χ (Section 4.2). There are two main uncertainties entering this result: i) the limit on the effective cross section stated in Eq. (4.4) is driven to a large extent by the ratio of average densities as introduced in Eq. (3.20), which we currently constrain only from a relatively small number of simulations; ii) translating this to a bound on the physical cross section involves an inevitable uncertainty in our theoretical prescription of the effect of SIDM, which we currently capture in the effective parameter ξ in (the commonly adopted) Eq. (2.4). While we believe that our quoted final bound does encapsulate this uncertainty, a more systematic investigation of the properties of simulated SIDM halos -e.g. by considering larger samples than what is shown in Fig. 7 -has the clear potential of resulting in even stronger and more robust bounds. As a byproduct of our analysis, we also revisited the scaling relation for the surface density evaluated further away from the center, at the scale radius of the NFW profile, and derived an updated version of the standard c-M 200 relation of ΛCDM cosmology, c.f. Eq. (3.13) and Fig. 2.
We note that for a velocity-independent self-interaction strength our upper bound is formally below, but still relatively close to, the ∼ 1 cm 2 /g that have been reported as a requirement to fully address the ΛCDM small-scale problems without invoking baryonic explanations. On the other hand, our bound can currently not (yet) constrain this idea in scenarios where the self-interaction cross section is velocity-dependent and drops significantly below 1 cm 2 /g for masses larger than those of dwarf galaxies. While still subject to some uncertainties, the reason for the relatively strong bound we report here is a combination of the large number of objects that we include in our analysis and the fact that the parameter κ introduced in Eq. (3.20) is so well constrained by simulations. As we have discussed in more detail in Section 4.3, furthermore, prospects for future significant improvements of these bounds are very good. This may finally settle the question of whether DM has astrophysically relevant self-interactions, thereby providing yet another example of how useful observational scaling laws are in astronomy.

A Datasets
For our analysis, we select objects where the density profile is reasonably well fit to both a cored profile and an NFW profile at large radii, and where the two resulting masses inside the best measured distance are consistent with each other. The cored profiles that we take into account are: • Cored NFW profile: The presence of baryons can significantly modify the "isothermal" solutions of the Jeans equation that feature σ v ≈ const. for r < r SIDM , up to the point where instead of cored profiles one finds cuspy profiles [17]. These extreme cases are obviously not among the objects that we consider, given that we demand a good fit to a cored profile. In order to keep the discussion of core radius vs. self-interaction radius as model-independent as possible, however, we would still like to make sure that we only select objects that are sufficiently DM dominated for the baryons not to affect our analysis, e.g. via poorly known mass-to-light ratios for stars. Concretely, we demand that at the largest measured radius the mass of the DM is at least 4 times larger than the mass of the baryons. This implies that we keep some objects that have a larger baryon content at the smallest observationally accessible distances; in Section 4.3, we discuss explicitly the (small) impact of those objects on our final results.
The following is a complete list of objects that we selected for our analysis, followed by a brief description of the characteristics of each object class: • Dwarf Spheroidal galaxies, [59,86]   Dwarf Spheroidal galaxies Dwarf spheroidal galaxies (dSphs) are satellites of the Milky Way and nearby galaxies. They are the objects with the highest known mass-to-light ratio and DM dominated even in the central parts, thus very suitable to study DM properties. There are so-called classical dSphs with a relatively large amount of stars, such that their DM profile can be reconstructed sufficiently well, see Fig. 8 for velocity dispersion profiles and Fig. 9 with the reconstructed mass profiles. To minimize the uncertainty in our analysis, we will only consider classical dSphs. The observable quantities in dSphs are the velocity dispersion along the line of sight σ los and the half-light radius r h . The main uncertainty on the mass distribution comes from the fact that one instead needs the full 3D velocity dispersion σ to calculate the mass inside   some radius r. One of the best mass estimates is obtained for the mass inside the half-light radius [93][94][95], as clearly illustrated in Fig. 9. To measure the slope of the central part of the DM profile it was proposed to generalize this idea to 2 or 3 distinct star populations, which led to the (still being disputed) claim of evidence for cores [96].
Spiral and irregular galaxies In spiral galaxies, the measurable quantities are rotation curves for stars and neutral hydrogen. As these objects are dominated by baryons in their central parts, the main uncertainty on the DM profile in this region comes from modelling the baryon mass. As a result, the DM profile can often be fitted (almost) equally well with NFW and cored profiles, see e.g. the rotation curves displayed in Ref. [89]. The best measurement of the DM profile derives instead from the flat part of the rotation curve, where DM dominates over baryons. Low surface brightness galaxies (LSBs) are diffuse galaxies with a surface brightness that, when viewed from Earth, is at least one magnitude lower than the ambient night sky. Most LSBs are dwarf galaxies, and most of their baryonic matter is in the form of neutral hydrogen gas, rather than stars. They appear to have over 95% of their mass in DM, and are DM dominated even in the central parts. As in the case of spiral galaxies, the observables are rotation curves v(r) of stars and neutral hydrogen. Reconstructed DM profiles appear to be more consistent with cored profiles [97].
Galaxy clusters The exceptional class of objects where we allow a larger baryon contribution in the inner parts are clusters of galaxies. This is because all clusters are baryondominated in the center by the existence of the brightest cluster galaxy and gas. On the other hand, it was claimed that some clusters show evidence for cores in the DM profile [51]. Such an evidence for cores in clusters is highly disputed, see e.g. the discussion in Ref. [98]. For the current work we simply take the cored fits of the DM profiles from Ref. [51], without any additional selection of the objects (but demonstrate in Section 4.3 that these objects only have a very minor impact on our results).

B Halo age
In this Appendix, we briefly describe how we implement the average halo age that enters the definition of the scale of self-interactions in Eq. (2.4). The "formation redshift" (the half-mass formation time) of a halo with present mass M is implicitly given in terms of the critical overdensity, δ crit = (ρ − ρ c )/ρ c , by [53] δ crit (z f ) = δ These expressions assume a flat universe, Ω m (z) = 1 − Ω Λ (z), and we take the current value of the matter density Ω m = Ω 0 m = ρ 0 m /ρ 0 c as measured by Planck [1]. The second step in Eq. (B.1) involves a fitting parameter from the accretion histories, f = 0.068, and the linear rms fluctuation in spheres of mass M . The latter is given by [54]  (B.6) Solving Eq. (B.1) for the formation redshift, one finally gets the halo age as [55] t age = t 0 − t(z f ), (B.7) where t 0 is the current age of the Universe and We plot this function in Fig. 10.