{'fraction_non_alphanumeric': 0.037170109427116856, 'fraction_numerical': 0.011410849084206214, 'mean_word_length': 4.5901864573110895, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Long range dark matter self-interactions and plasma instabilities

So far, the observed effects of dark matter are compatible with it having purely gravitational interactions. However, in many models, dark matter has additional interactions with itself, with the Standard Model, and/or with additional hidden sector states. In this paper, we discuss models in which dark matter interacts through a light vector mediator, giving rise to a long-ranged force between dark matter particles. Coherent scattering effects through this force can lead to the exponential growth of small perturbations, in analogy to electromagnetic plasma instabilities. These instabilities can be significant at couplings many orders of magnitude below those for which the usual particle-by-particle constraints on dark matter self-interactions apply. While this possibility has been noted in the literature, we provide the first systematic study of such instabilities, including the case where the mediator has finite mass. The latter is relevant for models of kinetically mixed `dark photon' mediators, which represent an important target for proposed dark matter detection experiments. Our analyses are of the growth of small perturbations, so do not immediately provide observational constraints on dark matter models - however, they do motivate further study of large regions of parameter space.


Introduction
There is overwhelming evidence that most of the universe's matter density consists of states other than Standard Model (SM) particles. All of the evidence for this 'dark matter' (DM) comes from its gravitational effects on astrophysical objects, and current observations are compatible with gravity being the only way in which DM interacts, either with itself or with the SM. Nevertheless, in many models of theoretical and experimental interest, DM must have additional interactions. Examples include WIMP DM, whose interactions with the SM radiation bath are responsible for its cosmological abundance, and the QCD axion, which has a self-interaction potential arising from its coupling to QCD. In all of these cases, DM interactions must be weak enough that its behaviour looks effectively non-collisional, on the galactic scales that we have observational data for.
In this paper, we will focus mainly on DM self-interactions. For models such as WIMP DM, the DM number density is low (the occupation number of particle modes in halos is 1), and the interactions between DM particles are short-ranged compared to the interparticle distance. Consequently, scattering is dominated by two-particle events. The standard DM self-interaction bounds from systems such as the Bullet Cluster [1] are based on this scenario. For DM masses in the GeV range, these constrain the interaction crosssection at around the barn level [2], comparable to SM nuclear cross-sections. The usual WIMP dark matter candidates, which have significantly heavier masses and weak-scale interactions, are safe from such constraints. In other models, DM self-interactions can be long-ranged. For the simplest examples, DM particles interact with a light mediator, whose mass is small compared to the inverse spacing between DM particles. 1 Quantitatively, where 10 GeV cm −3 corresponds to the highest DM halo density for which we have good evidence (unless otherwise noted, we use natural units with = c = 1). If the mass m of the mediating particle satisfies m n 1/3 DM , then coherent interactions between multiple DM particles can take place. These 'many-many' interactions are coherently enhanced, and can have much larger rates than standard 2 → 2 scatterings. Of course, they are also softer (transferring less energy and momentum), and their behaviour may be much more complicated.
There are various theoretical motivations for long-range hidden sector mediators. Examples include hidden sector counterparts of the photon (e.g. in scenarios such as Twin Higgs models [4,5]), light scalars [6][7][8][9][10], or other forms of 'modified gravity' (e.g. [11]). Such models are also important as an experimental target for low-threshold direct detection experiments, as we discuss in Section 5.
The effects of coherent, many-many scattering processes will depend on the form of the DM-mediator interaction. For a spin-0 mediator φ, the simplest kinds of interactions (e.g. φχχ, for a DM fermion χ) will be universally attractive or repulsive, depending on the sign of the coupling. This is basically the scenario of 'modified gravity' in the DM sector, a number of models for which have been explored in the literature (e.g. [6][7][8][9][10]).
For a spin-1 mediator A µ , the simplest forms of interactions, e.g. A µχ γ µ χ, are analogous to electromagnetism. 2 If the cosmological DM abundance is net neutral, with equal numbers of positive and negative charges, then (in the absence of large-scale charge separation) there will not be coherent forces between bulk matter, unlike in 'modified gravity' scenarios. In cosmological terms, inflation naturally sets any charge asymmetry to zero, and any processes regenerating it later must break the associated U (1). However, even for DM that is bulk neutral, there can still be significant 'plasma effects', driven by coherent many-many scatterings. Analogous coherent EM effects are extremely important in SM astrophysics, driving processes such as shocks, jets and other collective phenomena [12]. This paper will investigate scenarios in which DM interacts through a light vector mediator. As we will see, coherent effects can be important even with very small DM-mediator couplings; in some cases, many orders of magnitude smaller than those constrained by 2 → 2 scattering processes. This raises the possibility of much stronger constraints on such DM models, or equivalently, observational signatures at much weaker couplings. While the possibility of these effects has been noted in the literature [13], we provide a systematic analysis. In addition to considering massless vectors, we also consider mediators with a small but finite mass, as well as those with couplings to SM matter as well as to DM. These more complicated models are important phenomenological targets, in particular for proposed experiments aiming to detect light DM candidates.

Parametrics of plasma instabilities
Despite the simple underlying physics, plasma dynamics -the behaviour of charged matter interacting via a light vector mediator -is an immensely complicated subject. The great variety of different behaviours of SM matter in astrophysical systems, much of which is driven by the combination of EM and gravity, illustrates this richness [12]. Accordingly, this paper will not attempt to categorise all of the different kinds of behaviour that could occur in the presence of long-range DM interactions. Instead, we will focus on identifying the smallest couplings for which coherent effects would be important in astrophysical settings. This should help guide future investigations of such models, and identify regions of parameter space where we might expect deviations from collisionless behaviour.
In many circumstances, the first effects to become important, as we increase the coupling from zero, are exponentially-growing instabilities. For a non-relativistic plasma, the fastest-growing of these instabilities are 'electrostatic' ones, in which small charge separation perturbations grow, by extracting kinetic energy from a non-Maxwellian velocity distribution [14]. The simplest case is a 'two-stream' instability, in which two uniform 'streams' of plasma pass through each other at high velocity (relative to their individual velocity dispersions). Perturbations, in the form of longitudinal charge-density waves, are exponentially amplified, with growth rate ∼ ω p , where ω 2 p = q 2 n/m q is the 'plasma frequency', q the charge of individual particles, m q their mass, and n their number density [15].
Taking our 'plasma' to be composed of dark matter particles, interacting through a hidden-sector vector mediator, this two-stream scenario approximates a collision between the DM halos of two galaxy clusters, such as the Bullet Cluster [1,16]. In that system, the central densities of the colliding clusters' DM halos are modelled to be 0.1 GeV cm −3 (see Section 3). This DM density gives a plasma frequency of ω p = g 2 ρ χ m 2 χ 4 × 10 −7 yr −1 ρ χ 0.1 GeV cm −3 GeV m χ g 10 −17 (1.3) where m χ is the DM mass, ρ χ is the DM density, and g is the DM-mediator coupling. The scale radius of the DM halos is ∼ 100 kpc, and their relative velocity is ∼ 3000 km s −1 , giving a crossing time of ∼ 100 kpc/(3000 km s −1 ) ∼ 3 × 10 7 yr. Comparing ω −1 p to this crossing time, we see that the exponential growth of perturbations could occur even for g ≪ 1.
This estimate shows that coherent effects may become important at DM-mediator couplings orders of magnitude smaller than those at which 2 → 2 scatterings have significant effects. In the remainder of this paper, we analyse the growth of small perturbations more generally and quantitatively. In particular: • We review the linear theory of instabilities in a homogeneous plasma (Section 2).
Standard electromagnetic theory can be applied directly to the case of DM interacting via a massless hidden-sector mediator, simply by modifying the charges and masses of the particles.
• We extend this theory to the case of a massive vector mediator (Section 2). We find that mediator masses m ω p /v th suppress the growth of perturbations, where v th is the velocity dispersion of the plasma. For velocity distributions without multiple peaks (such as an anisotropic Maxwellian distribution), only magnetic instabilities exist, and these are suppressed for m O(ω p ).
• We identify regions in the m χ , m A , g χ parameter space (where m χ is the DM mass, m A is the mediator mass, and g χ is the DM-mediator coupling) where we expect small perturbations to grow in astrophysical systems (Section 3). At small m A , the potential constraints on g χ are many orders of magnitude stronger than those from 2 → 2 scattering (this is illustrated in Figures 6 and 7). We also comment on the potential for improved constraints when the DM is so light (m eV) that it has occupation number 1 in halos (Section 6).
• We study models in which DM interacts non-gravitationally with both itself and with SM matter, through a 'dark photon' vector mediator that kinetically mixes with the SM photon (Section 5). This is one of the simplest such possibilities, and represents an important target for low-threshold direct detection experiments. In large parts of the theoretically and experimentally interesting parameter space, we find that DM-DM plasma instabilities would be expected to grow exponentially in astrophysical systems ( Figure 10). This motivates further investigations to determine whether such models are astrophysically viable. We also discuss why such effects may not be as important for the case of millicharged DM (Section 4).

Previous literature
To our knowledge, the first paper to raise the possibility of strong constraints from coherent DM-DM scattering was [13]. This considered DM interacting through a massless, purely dark-sector vector mediator, and focussed mainly on the constraints and signatures arising from 2 → 2 scattering. However, it also commented on the possibility of a magnetic Weibel instability (see section 2.3), and noted that such phenomena may provide significantly stronger constraints than 2 → 2 scattering, while leaving further investigations to future work.
In unpublished work, [17] identified the electrostatic two-stream instability as being faster-growing than the Weibel instability (again, for the case of a massless hidden-sector mediator), and estimated the parameter space where constraints from Bullet Cluster type systems might exist. In addition, they conducted preliminary simulations of such cluster collision scenarios.
[18] also identified the electrostatic two-stream instability as the fastest-growing one, and argued that since this instability would make DM effectively collisional, 'models where all of DM is charged are ruled out by observations of cluster collisions unless the charge is extremely small' (though they did not attempt to quantify this). Their main focus was a scenario in only a sub-component of DM is self-interacting through a massless hidden-sector mediator. In [19], cluster collisions with such a subcomponent were simulated, assuming that it behaves like a fluid with some effective viscosity.
Coherent interactions involving DM have been investigated extensively in the context of millicharged DM, where DM is taken to have a (very small) charge under Standard Model EM. As discussed in Section 4, the most important interactions are generally 'onemany', with individual DM particles scattering off the large-scale EM fields created by SM charges and currents. While some of this literature (e.g. [20]) claims that similar bounds apply to models with massive mediators, such as kinetically mixed dark photons, this is only true when the mediator mass is extremely small (see Section 5). There have also been investigations of whether 'many-many' collective instabilities can be important for millicharged DM [21], but futher work would be needed to establish whether this could occur for couplings not excluded by other effects.
DM coupled to a light vector mediator is considered in many other papers (e.g. [22][23][24][25]), but in most cases, they do not consider plasma effects (or if they do, only insofar as they affect Coulomb scattering through Debye screening -see Section 5.1). [26] does consider plasma effects in a model with a massive hidden-sector mediator, but focusses on the properties of the DM plasma near Earth (in the context of direct detection experiments), rather than on systems such as cluster collisions.

Plasma instabilities
In this section, we will discuss the linear theory of plasma instabilities. Astrophysically, the DM distributions that arise from collisionless evolution will often be very far from thermal equilibrium (for example, the basically-bimodal velocity distributions in cluster collisions). Given strong enough interactions, such a distribution will relax towards equilibrium. 3 For some distributions, the initial stages of this relaxation are driven by exponentially-growing instabilities, which grow by converting particle kinetic energy into field energy.

Instabilities in a uniform plasma
We will begin by considering perturbations of a spatially uniform plasma, in a freely falling frame, with no large-scale currents or charge separation in the initial state. In real astrophysical systems, spatial inhomogeneities and gravitational forces are, of course, often important. However, if the spatial scale of a perturbation is small, these approximations can be good, and as we will see, the spatial scales associated with instabilities are often small enough. If the initial state does not feature large-scale currents or charge separation, then (for large enough DM number densities that particle number effects are not important), these must have arisen through interactions becoming important at earlier times. Consequently, to analyse the earliest point at which large-scale vector fields arise, we can start from approximately zero-field conditions.
In terms of the distribution functions f s (v, x, t) = dNs d 3 vd 3 x of the particle species s in the plasma, the collisionless Boltzmann equation 4 (also referred to as the Vlasov equation) is [29] ∂ is the acceleration due to the Lorentz force, with q s and m s the charge and mass of species s. This expresses that the phase space density does not change along particle trajectories. The charge density, ρ = s q s d 3 vf s , and the current density, J = s q s d 3 vvf s , give the source term J ν in the Proca equations (∂ µ ∂ µ − m 2 )A ν = J ν (where m is the mass of the vector mediator, and we use the + − −− signature). Together, the Vlasov and Proca equations determine the behaviour of the plasma distribution function.
For a particular configuration of charged particles, the distribution function will be a set of delta-functions corresponding to the positions and momenta of these particles, and solving the Vlasov-Proca equations will correspond to solving the interacting N -body system. To simply things, we need to perform some kind of averaging. One approach is to coarse-grain in position and velocity space, so that there are many charged particles in each phase-space 'bin' -another is to track the evolution of an ensemble-averaged distribution function, where we average over e.g. different possible starting configurations of the particles [29]. In either case, for the averaged solution to be a good approximation to the particular solution, we need the particles in the latter to be sufficiently dense in phase space. We will return to this condition in Section 2.2.
If our unperturbed state is a spatially uniform distribution function f s , with no vector field background, then the linear-order equations for a small perturbation δf s , A µ are We will work with non-relativistic velocity distributions, so will take γ 1. If we decompose the perturbations into Fourier modes, so δf s (t, x) = δf s e −i(ωt−k·x) etc, then we can express δf s in terms of A µ , giving a linear equation for A µ ; ( where Π νµ is the leading order response function (the 'self-energy', in the language of thermal field theory [30]). For a given k, there will be some particular ω for which this equation can be satisfied. These (k, w) pairs give the dispersion relation for perturbations.
To take the simplest possible example, a cold, uniform plasma has f s = n s δ 3 (v). We can evaluate equation 2.3 via integration by parts; if we assume that f decays fast enough at large velocities, then where ω s ≡ q 2 s n s /m s . For m = 0, this gives the usual dispersion relations for transverse and longitudinal excitations in a cold plasma [31]. Some more realistic examples, such as Maxwellian velocity distributions, can also be treated analytically (see appendix A).
One of the simplest examples that gives rise to exponentially growing perturbations is the 'two-stream' instability. If we consider a plasma consisting of positive and negatively charged species with the same mass (as we expect for symmetric DM models), this has ), corresponding to two cold 'streams' of plasma passing through each other. For longitudinal perturbations, i.e. those with k A v 0 , we have where ω p = q 2 n p /m is the plasma frequency of each of the streams, with the charge of the two species being ±q and their mass being m. If we start with some localised perturbation at an initial time, we can decompose it into spatial Fourier modes with real k. Consequently, there are exponentially growing perturbations if there are solutions of the dispersion relation with Im ω > 0 for real k. As we discuss in Section 2.2, such solutions The plasma is assumed to consist of positive and negatively charged species with the same mass, and having the same velocity distributions. We take ω 2 p = q 2 n/m q , where n is the number density in each stream, ±q is the charge of each species, and m q is their mass. exist for m ≤ ω p /v 0 , and the fastest growth rate is ω i 0.5ω s (attained for k ω s /v 0 , and m k). We will discuss the physical understanding of this instability (and others) in detail below, but to start with, we can understand the parametrics of the growth rate simply by considering the energetics of the instability. The vector field's energy density is u A = 1 2 (E 2 + B 2 ) (or larger for a massive mediator -see equation 2.11). The force exerted on a charged particle is F = q(E + v × B), so over a time δt, the magnitude of the change in the particle's velocity is at most δv ∼ qδt m (E + v × B). For a collection of particles with charges q i and velocities v i , the rate of energy transfer from particle KE to field energy is If we suppose that our perturbation is growing exponentially, and we take δt to be the e-folding time, then we must have δP ∼ P and δU A ∼ U A , so δP δt ∼ U A . For a perturbation where u A is dominated by the electric field, such as the longitudinal modes in equation 2.7, we have |δP δt| q 2 m N (δt) 2 E 2 , where N = nV is the number of particles in the volume V . Since |δP δt| needs to be u A V ∼ 1 2 E 2 V for the perturbation to grow, we have For a magnetic-field-dominated perturbation, the maximum growth rate is suppressed by  Figure 2. Dispersion relation for longitudinal oscillations in a cold, uniform plasma, with a massive mediator. In the left-hand panel, the mediator mass is m = 0.1ω p . Longitudinal plasmons, which for a massive mediator have ω = ω p , are mostly unaffected (apart from a small avoided crossing). At small ω, there is also the φ const Goldstone mode (which would be pure-gauge for a massless mediator), and at high ω, there are relativistic φ excitations. In the right-hand panel, we have m = 2ω p . Here, longitudinal plasmons are strongly affected for k m, linking with the Goldstone dispersion relation at smaller k. At higher ω, the dispersion relation is approximately that for a free massive particle. some power of v/c 1 (for a non-relativistic plasma). 5 Consequently, the ω i ∼ ω p growth rate obtained from the two-stream instability is as fast as perturbations can grow, starting from a background without larger fields present (equivalently, it is the fastest rate at which the overall field energy can increase).

Electrostatic instabilities
As we will see, unstable modes in a non-relativistic plasma generally have |k| |ω|; parametrically, for the velocity distribution of the plasma to be significant, we need |kv| ω. Consequently, instabilities in which the electric field perturbations are much larger than the magnetic field perturbations are dominantly electrostatic, with |φ| |A| (since B = ∇ × A and E = −∂ t A − ∇φ). In this case, we can simply the Vlasov-Proca equations to a one-dimensional integral. Reinstating c for clarity, we have so taking E −ikφ,  so the behaviour of electrostatic perturbations with a particular wavevector k is equivalent to longitudinal perturbations with the appropriate projected velocity distribution (see appendix B). For a massless mediator, longitudinal perturbations have been extensively analysed in the plasma theory literature [15]. In the symmetric, two-stream scenario, a basic physical picture of the instability can be obtained (following [15]) by considering the longitudinal oscillations of each of the streams individually. For a stream with velocity v in the lab frame, oscillations with ω = ω p in the stream's rest frame have ω = k · v ± ω p in the lab frame 6 , and the 'slow' oscillation (with smaller frequency) has less kinetic energy then the unperturbed beam, ∆KE < 0. If we take k such that |k · v| = ω p , then the slow oscillations for both streams are stationary in the lab frame. If the two oscillations are in phase, they reinforce each other, extracting KE from the streams to increase the E field energy, and 6 In the case of symmetric DM, each stream will be composed of equal numbers of positive and negative charges, with both species having the same velocity distribution. If each has number density n± then the plasma frequency for the stream is ω 2 p = 2q 2 n±/m. In the case of an SM plasma, we might instead have streams of electrons on top of a uniform, slowly-moving background of ions. The important quantity is simply the self-energy that the distribution function gives rise to, which determines the dispersion relation.
giving rise to an instability. This is the relevant picture for k · v ω p , where the peak instability rate is obtained. As Figure 1 shows, there are also growing perturbations with small k, for which ω i 0 as k 0. These can be thought of as arising from out-of-phase slow oscillations in each stream, such that (in the lab frame) the charge density perturbations almost cancel each other out.
A similar picture applies in the case of a massive mediator, with the difference that, for m ω p and k m, the dispersion relation for longitudinal oscillations in a stationary plasma is modified (see Figure 2). Consequently, the peak instability rate is modified for m ω p /v.
Even without solving equation 2.10 directly, we can determine whether an instability exists using the Penrose criteria [32] (see appendix B). If the two streams have Maxwellian (projected) velocity distributions, then for a massless mediator, an instability exists if v 0 1.3σ, where σ is the velocity dispersion for each stream, and 2v 0 is the closing velocity. For a massive mediator, even for well-separated stream velocity distributions (i.e. closing velocity velocity dispersions), we need m 1.3ω p /v 0 in order for an instability to exist. Figure 1 illustrates a particular example of two Maxwellian streams moving through each other (with parameters roughly analogous to a cluster collision scenario), taking the closing velocity to be 2×10 −3 c, and the velocity dispersion for each stream to be σ v = 10 −4 c. The height of the surface corresponds to the fastest instability rate for each k, which can be found from equation 2.10 (in fact, the figure uses the full Vlasov-Proca equation, rather than the longitudinal approximation, though the results are almost identical). As we go to more oblique k, the projected velocity distribution has a smaller closing velocity, so past the v p 1.3σ threshold, where 2v p is the projected closing velocity, the instability does not exist (though magnetic instabilities may still be present, as we discuss in Section 2.3).
For a massive mediator, if the closing velocity is significantly greater than the velocity dispersions, then going to oblique k can be helpful, since it decreases v p , and so increases ω p /v p . This is illustrated in Figure 3, which shows the effects on increasing m on the instability rate as a function of k. For m ω p /σ, the instability is suppressed entirely.
This behaviour matches what we would expect from basic energetic considerations. The energy density associated with the A field is where the second equality is for an electrostatic perturbation. Consequently, we expect the mediator mass to become significant when m k, as we observed above. As mentioned in Section 2.1, for the collisionless Boltzmann equation to be a good approximation to the evolution of a particular plasma configuration, the particle density in phase space needs to be sufficiently high. If the largest wavenumber of interest is k, then having many particles per ∼ k −3 volume means that we have a sensible spatial coarsegraining. Within each k −3 , we ideally want enough particles to resolve the important features of the velocity distribution. For the two-stream instability, the relevant condition is n χ (ω χ /σ) 3 . If we impose the 2 → 2 Coulomb scattering bound (equation 5.5), and also that g ≤ 1/(4π) (when m χ is large enough that the Coulomb bound is ineffective), then Consequently, the dark matter number density is always high enough that, in a spatially uniform plasma, there would be a good coarse-grained description. We discuss how well the approximation of spatial uniformity applies to astrophysical systems in Section 3. The other way in which the collisionless Boltzmann equation can break down is if collisions are important. For hard collisions, the obvious estimate is that, if the collision rate for a particle is comparable to the growth time for the instability, then the instability will be suppressed. This can be confirmed by explicit calculations [15]. For Coulomb collisions, the rate of soft collisions is enhanced by a Coulomb logarithm factor (Section 5.1), but this alters things by a factor of O(10) at most [13,24]. As we noted above, since only a few hard collisions per DM particle during the lifetime of universe are enough to significantly affect e.g. halo shapes [13,24], neglecting collisions will be a good approximation for couplings that are not already excluded by 2 → 2 scattering calculations.

Beam-plasma
If we consider asymmetric streams, as illustrated in Figure 4, the behaviour of instabilities can be somewhat different. In the plasma literature, these are known as 'beam-plasma' instabilities, with the 'beam' being the less dense stream and the 'plasma' the more dense one [15]. For us, they can serve as a model for astrophysical situations such as a denser, lower-velocity-dispersion subhalo (the 'plasma') moving through a larger DM halo (the 'beam'), as discussed in Section 3.
Setting up our problem (for a massless mediator, our calculations will follow those in [15]), we will take the denser stream to be at rest, with plasma frequency ω p , and the less dense stream to have relative velocity v 0 , with plasma frequency ω 1 . If we ignore the velocity dispersion of the streams, the dispersion relation for a massless mediator is If we assume the density ratio is small, Consequently, for ω ω p , the modes are stable, while for ω ω p , we have ω i ∼ √ Rku. The maximum instability rate is obtained for ω ω p , for which where n is an integer. Consequently, the maximum growth rate is Neglecting the velocity dispersions of the streams will be a good approximation if the phase velocity of the excitation is well outside the velocity distributions of both streams. In the plasma literature, this is referred to as the 'cold-beam' case. Since ω/k − v 0 = δ/k, this will be true if the fractional velocity dispersion of the beam is σ v 0 max δ kv 0 ∼ R 1/3 , and the plasma velocity dispersion is v 0 . The physical picture for the fastest-growing instability is that a longitudinal oscillation in the denser stream, which has frequency ω p , interacts with a ∆KE < 0 forced oscillation in the 'beam', extracting kinetic energy. For kv 0 ω p , we have a similar picture as for the two-stream instability, with partially-cancelling perturbations in the beam and the plasma. As we take R → 1, we move smoothly towards the two-stream case.
In many circumstances, the beam's velocity dispersion will be large enough that the 'cold-beam' approximation does not hold. For example, in the subhalo case, the velocity dispersion of the larger halo is generally comparable to the orbital velocity of the subhalo through it, so σ ∼ v 0 . The behaviour in this 'hot-beam' scenario is rather different.
Since we need to take the beam's velocity distribution into account, the dispersion relation for longitudinal oscillations becomes (again, for a massless mediator) (we will see that as long as the plasma's velocity dispersion is u, treating it as cold will not affect the most unstable modes). If we write ω/k = u r + iu i , then we can write the integral over v as an integral over a contour slightly below the real axis, which for u i sufficiently small is approximately For the imaginary part of equation 2.16 to be zero, we need For R 1, the real part of the RHS of equation 2.16 is dominated by the first term ω 2 p /ω 2 , so we need ω r ω p . Consequently, For a generic beam velocity distribution, with velocity dispersion ∼ σ, the maximum value For this to be valid, we need f (u r + iu i ) f (u r ); for a simple beam velocity distribution such as a Gaussian, this requires u i σ. From equation 2.20, this holds for for R σ 3 /v 3 0 . This is simply the opposite of the cold-beam condition derived above, showing that the hot-beam and cold-beam regimes are complementary. The maximum instability rates also have the same parametric value at R ∼ σ 3 /v 3 0 , consistent with this picture. The different behaviour in the hot-beam regime corresponds to a somewhat different physical picture of the instability, compared to the cold-beam case. In particular, the approximation of the counter integral used above is very similar to the usual derivation of Landau damping [14,29]. For the latter, the physical picture is that beam particles with velocities very slightly smaller the phase velocity ω/k of a plasma oscillation are accelerated by it, while those with velocities very slightly larger are slowed down. If the gradient of the velocity distribution is negative, then there are more slightly-slower particles than slightlyfaster particles, and energy is on average transferred from the plasma oscillation to particle KE. If the gradient is positive, then the opposite occurs; energy is transferred from particle KE to the plasma oscillation, resulting in an instability. This 'inverse Landau damping' picture breaks down when the instability growth rate becomes too fast, so that 'slightlyslower' and 'slightly-faster' encompass too wide a range in the velocity distribution.
The top-right panel of Figure 4 shows the instability behaviour for an example velocity distribution with R = 10, and σ = v (i.e. in the hot-beam regime). While the ω r values of the unstable modes are not shown, these are close to ω p . As the figure illustrates, we need k · v ω p for there to be an instability, corresponding to the phase velocity of the perturbation being in the increasing part of the v z distribution.
For a massive mediator, the same basic physical picture applies, once the modified dispersion relation of longitudinal oscillations is taken into account. Equation 2.16 is modified to From the real part of this equation, we need to have ω ω p (1 + m 2 /k 2 ) −1/2 . This corresponds to the lower branch of the dispersion relation (for the dense plasma) in the right-hand panel of Figure 2, with the usual relation ω ω p being modified for m k. Parametrically, we expect instabilities to be suppressed entirely for m ωp σp . For longitudinal oscillations, we can confirm this using the Penrose criteria; from appendix B, instabilities exist for where v c is the point at which f v = 0, balanced between the falling f from the plasma velocity distribution and the rising f from the beam. The second equality applies if both the beam and plasma have Maxwellian velocity distributions. The basic physical picture is that, for given k, increasing m pushes ω r smaller, until u r = ω r /k is no longer on the rising part of the velocity distribution (see Figure 12). This behaviour is illustrated in Figure 4, which shows the effects of increasing m on the instability rate as a function of k. Since the velocity distribution of the dense component is isotropic, and the dense component's velocity dispersion sets the largest k for an unstable mode, it does not help to go to oblique k.

Magnetic instabilities
The electrostatic instabilities discussed above relied on velocity distributions with sufficiently well-separated peaks. As well as these, there can also be 'magnetic' instabilities, in which the magnetic field energy dominates the electric field [33,34]. In situations where both magnetic and electrostatic instabilities exist, the growth rate of the former is generally suppressed relative to the latter by ∼ v s /c, where v s is the stream velocity. However, magnetic instabilities exist for a wider range of velocity distributions, including those which are very close to Maxwellian.
The simplest examples of magnetic instabilities can be seen in the (cold) two-stream velocity distribution [34]. If, instead of looking at k v 0 perturbations, we take k ⊥ v 0 , then the response function is (taking the streams to be in the z direction, and k ∝x)  [34]. These currents do work on the field via the E · J term (since E is in the vertical direction), enhancing the field energy and leading to an instability.
For A ⊥ẑ, we have the usual dispersion relations for plasma oscillations (as we would expect, since non-relativistic motion in theẑ direction has not effect on those). For A ẑ, the dispersion relation becomes (where we have reinstated c for clarity) which has solution Expanding in small v 0 /c, we have a solution so there are growing modes, with max ω i ω p v 0 /c (attained at k max(ω p , m)). The physical picture behind this instability is illustrated in Figure 5. A transverse magnetic field bends the trajectories of particles in the streams, creating current sheets in the stream direction. These currents do work on the electric field, growing the instability.
More generally, if the velocity distribution is separable in the x and z directions, so (assuming symmetric positive and negative species) with dv where σ 2 z = dv z v 2 z g z (v z ) and ω 2 p = 2q 2 ± n ± /m ± , and G is defined as in appendix B. We also have Π zx , Π xz ∝ dv z v z g z (v z ), so to simply the analysis, we can make Π diagonal by boosting in the z direction so thatv z = dv z v z g z (v z ) = 0.
For the case of a Maxwellian velocity distribution g x , we have (for u with positive imaginary part) is the plasma dispersion function [29] (inserting this into equation 2.28 reproduces the dispersion relation given in the original paper of Weibel [33]). 7 From the dispersion relation we can find that for a Maxwellian v x distribution, an instability exists for k ≤ k max , where If the RHS is negative, then there is no instability. In many astrophysical situations, the velocity anisotropy will be small, σ z /σ x 1. In this case, we can expand the plasma dispersion function in small ω i , and obtain (for a massless mediator) For a mediator appreciably lighter than the threshold mass ω p , the same ω i expression will hold approximately.
Comparing these results to the electrostatic instabilities we discussed previously, the mass range in which instabilities existed for a two-stream velocity distribution was m ω/v th . Consequently, in a scenario such as a cluster collision, the magnetic instability is present for a smaller range of mediator masses (though it can be present for other velocity distributions where electrostatic instabilities do not exist). When both instabilities are present, the parametric growth rate ω i ∼ ω p v/c, as compared to ω i ∼ ω p for an electrostatic instability, means that magnetic instabilities grow slower in non-relativistic plasmas. This is illustrated in Figure 1, for which the computation includes instabilities of all types. While for some k (e.g. k x) the fastest-growing instability will be magnetic, the closing velocity of ∼ 10 −3 c means that these growth rates are too small to be visible on the plot.

Astrophysical systems
As discussed in the Introduction, our goal in this paper is not to fully calculate the consequences of DM self-interactions. Instead, our aim is to identify the parameter space 7 To understand the parametric behaviour of the response function, it can also be useful to consider simpler situations, e.g. a top-hat velocity distribution, for which Figure 6. Estimated parameter space for which small perturbations would experience significant growth in astrophysical systems, for a vector mediator of mass m A , coupling to DM particles of mass m χ with coupling g. The shaded parameter space above the broken red lines, labelled with mediator mass m A , is the space in which perturbations would grow for that mediator mass. The region above the solid red line is the region in which perturbations would grow for sufficiently light mediators. The astrophysical situations we consider are cluster collisions such as the Bullet Cluster, and dense DM subhalos (such as dwarf galaxies) moving through a larger halo (see Section 3). The 'Coulomb' region shows the parameter space in which 2 → 2 Coulomb collisions would have a significant impact in DM halos (see section 5.1). Note that, for DM masses m χ eV, the DM occupation number in Galactic halos must be large, which may result in additional coherent scattering effects (Section 6).
in which self-interactions are likely to be significant, so that future work can investigate this space more thoroughly, and in particular, identify whether there are observational consequences.
To map out this parameter space, we want to identify the astrophysical systems where, as we increase the mediator coupling from zero, the first significant effects arise. In this section, we identify different classes of candidate systems -cluster collisions, subhalos within larger halos, and halos with anisotropic velocity distributions -and discuss their properties.
Cluster collisions, as observed in systems such as the Bullet Cluster, are standard laboratories for constraining DM self-interactions [2]. The ionised gas in each cluster interacts during the collision, through SM plasma processes. In constrast, the post-collision DM distribution, as mapped out by gravitational lensing, is compatible with no non-gravitational interactions [1,2].
To estimate when DM-DM plasma instabilities would arise during such collisions, we need some model for the DM distributions in the colliding clusters. [16] estimates these distributions, for the Bullet Cluster collision, using gravitational lensing information. As-suming cored DM density profiles 8 , their best-fit parameters have core radii ∼ 100 kpc, and central densities ∼ few × 0.1 GeV cm −3 , with inferred halo velocity dispersions of σ ∼ 800km s −1 . Simulations of the gas dynamics in the collision suggest that the relative velocity of the halos is ∼ 3000 km s −1 [2]. These parameters give a central plasma frequency of compared to a crossing time of ∼ 100 kpc/(3000 km s −1 ) ∼ 3 × 10 7 yr. If we take a conservative threshold of O(100) e-folding times, we would expect plasma instabilities to grow significantly for if the vector mediator is sufficiently light. The latter condition should be met for The corresponding region in m χ , g parameter space is illustrated in Figure 6 (for various m A ), and in m A , g parameter space (for fixed m χ ) in Figure 7. Compared to constraints from 2 → 2 Coulomb scatterings (discussed in Section 5.1), plasma instabilities could be important at far smaller couplings, for light enough mediators. From equations 3.2 and 3.3, the threshold for the mediator mass to alter the parameter space in which perturbations grow is m few × 10 −26 eV, basically independent of m χ . For a two-stream instability, the wavevector for the fastest-growing mode (with a massless mediator) has k ∼ ωp v 0 ∼ few × ω i v 0 . For the threshold coupling we considered, which has ω i ∼ O(100)/t cross ∼ O(100)v/R cross , this implies that k ∼ O(100)/R cross . Consequently, the wavelength of the growing modes is small compared to the scale of the DM halo. If the DM halo is smooth, then our approximation of a spatially-uniform, freelyfalling plasma should be a reasonable one (though more careful analysis, and potentially simulations, would be necessary to confirm this). For a massive mediator, the k of the fastest-growing modes is even higher. On the other hand, if the DM halo has significant substructure, e.g. a dense central spike, or a large number of 'miniclusters', this may have important effects.
For more massive mediators, we are limited by the m ω χ /v th condition on the existence of an instability. To maximise the mass range we can cover, we want DM systems with large n χ and small v th . Of the DM systems we have observational evidence for, the most extreme values of both n χ and v th are obtained in small subhalos of galactic halos, such as DM-dominated dwarf galaxies. Even assuming cored dark matter profiles, these can have central DM densities of 10 GeV cm −3 and velocity dispersion ∼ 10 km s −1 [35][36][37] (it should be noted that inferences about the DM halos of dwarf galaxies are controversial, 8 profiles with central cusps, such as NFW, would lead to higher central DM densities -to be conservative, we consider cored profiles (gravitational lensing data does not have the spatial resolution to distinguish these possibilities). Figure 7. Estimated parameter space for which small perturbations would experience significant exponential growth in astrophysical systems, for a vector mediator of mass m A , coupling to DM particles of mass m χ with coupling g. In this plot, we fix m χ = MeV, and vary m A (the inset plot shows m χ = 10 GeV). The dashed curve corresponds to the parameter space in which instabilities would grow in a Bullet Cluster-like collision, while the dotted curve corresponds to the parameter space in which instabilities would grow for a dwarf-galaxy-like subhalo moving through a galactic halo, as discussed in Section 3. The dot-dashed grey line shows the estimated instability region for magnetic instabilities, in halos with anisotropic velocity distributions. The horizontal parts of these curves are set by the growth rate of the instability being fast compared to the crossing time in the system, while the sloping parts are set by the mediator mass being small enough for instabilities to exist. The 'Coulomb' region shows the parameter space in which 2 → 2 Coulomb collisions would have a significant impact in DM halos (see section 5.1). The inset illustrates how these regions move to higher couplings as m χ increases; the threshold g values for the plasma instabilities increase ∝ m χ , while the Coulomb threshold value increases approximately as ∝ m with some studies claiming that Milky Way dwarf galaxies may contain no dark matter at all [38,39]). These parameters correspond to a plasma frequency of As illustrated in Figure 7, this means that for large enough g, instabilities persist up to larger m than in the cluster collisions case. The low density and high velocity dispersion of the galactic halo in which the subhalo moves means that we are in the (hot) beam-plasma scenario, so the instability rate is suppressed by the density contrast between the halos. Taking a standard model for the Milky Way DM halo (e.g. [40]), the density at radii 10 kpc is For typical dwarf galaxies at tens to hundreds of kpc from the galactic centre, this gives a MW halo density of ∼ 10 −2 − 10 −3 GeV cm −3 . Consequently, the maximum instability rate is In the case of cluster collisions, the available timescale for perturbation growth is set by the crossing time of the DM halos (as noted above, further investigation would be needed to understand the effect of spatial inhomogeneities in these halos). For the case of a subhalo moving through a larger halo, the crossing time (e.g. ∼ kpc/10 −3 ∼ 4 × 10 6 yr) is significantly smaller than the timescale over which the subhalo's surroundings will significantly change (e.g. ∼ 100 kpc/10 −3 ∼ 4 × 10 8 yr). Consequently, it is plausible that perturbations can build up in the subhalo over timescales significantly longer than the crossing time. This is compatible with the fact that the group velocity of growing electrostatic perturbations, in the hot-beam case, is suppressed relative to the stream velocity, since the real part of the frequency is set by the dispersion relation for plasma oscillations . Here, we leave a proper analysis to future work, and conservatively restrict ourselves to perturbations that grow substantially during a single crossing time. With this assumption, equation 3.7 shows that for very light mediators, cluster collisions will be more prone to instabilities at smaller couplings.
The unstable modes in a beam-plasma setup have k ∼ ω p v 0 , and ω i ∼ Rω p . At the threshold we took, ω i ∼ O(100)v 0 /R cross , so k ∼ R −1 O(100)/R cross , and the separation between the mode's wavelength and the scale of the halo is larger than in the two-stream case. For more massive mediators, the threshold coupling is higher, and the fastest-growing modes have larger k (if a significant instability is present).
Cluster collisions and dense subhalos are the most obvious examples of bimodal DM velocity distributions, of the kind that lead to electrostatic instabilities. For standard DM halos, the (virialised) velocity distribution is expected to be roughly Maxwellian, so electrostatic instabilities should not be present. However, as we reviewed in Section 2.3, for light mediators, magnetic instabilities exist even for close-to-Maxwellian velocity distributions. In fact, numerical simulations indicate that DM halos may, outside their central regions, have quite anisotropic velocity distributions [41][42][43][44][45][46], with the radial velocity dispersion being ∼ 1.5 times larger than the transverse ones. This could lead to Weibel-type instabilities.
As per Section 2.3, the growth rates of magnetic instabilities are velocity-suppressed, and the maximum mediator mass for which they exist is ∼ ω p (for O(1) velocity anisotropies) rather than ∼ ω p /v th . Consequently, a basic estimate of the parameter space in which they would arise within galactic halos gives a subset of the electrostatic estimate, as illustrated in Figure 7. In addition, the evidence for anisotropic DM distributions is somewhat indirect [47,48], compared to the evidence for bimodal velocity distributions in subhalos and cluster collisions. Nevertheless, the coupling threshold estimates for cluster collisions and for anisotropic halos are quite close at small mediator masses, indicating that a more detailed investigation is worthwhile. Bound structures in the universe most likely came about through a process of hierarchial structure formation, in which smaller structures formed first, and then merged to form successively larger structures [49]. During this process, the dark matter velocity distribution was often far from Maxwellian, allowing instabilities to occur. In some sense, the cluster mergers that we observe are the late-time part of this process. Generically, we might expect these later events, which occur over longer timescales, to allow more time for instabilities to grow, and so to be the first that become important as we increase the coupling from zero. This is why we have focussed on such events (along with the fact that we have observations of them). These arguments are, of course, schematic, and would need to be checked more carefully.
Going further backwards in time, we can ask whether DM-DM instabilities could be important during the radiation era. We do not expect the DM velocity distribution then to have been significantly bimodal, but it could have been slightly anisotropic, raising the possibility of magnetic instabilities. However, even being optimistic, it does not seem likely that these would be competitive with late-time systems, in terms of possessing instabilities at small couplings. Since H ∝ z 2 then, while ω p ∝ n 1/2 DM ∝ z 3/2 , the increased instability rate possible in denser environments is more than offset by the shorter instability time available, and the most favourable ratios occur at later times. Performing basic estimates at around the time of matter-radiation inequality, it seems likely that the late-universe systems we discussed above would allow instabilities at smaller couplings. There is, of course, a separate question of whether SM-DM couplings could lead to interesting earlyuniverse behaviour -in particular, the pressure in the baryon-photon fluid makes it behave differently from DM, giving rise to a relative DM-baryon bulk velocity [50]. We leave such investigations to future work.
As we have emphasised, the instability regions in Figure 7 do not represent observational constraints. Nevertheless, it is interesting to consider how one might be able to obtain constraints, at least in the close-to-threshold regime (that is, assuming that the relevant process is the first time in the evolution of the system that collisions become important, so that we can take the initial state to be the same as for collisionless CDM). For cluster collisions, a potential issue that the gravitational lensing measurements are compatible with O(1) effects on the DM halos, e.g. ∼ 30% mass loss from the smaller halo [51]. While our calculations have been purely at the level of small perturbations, which can be treated using a linear approximation, understanding the eventual consequences once perturbations become large would require a more complete analysis. It should be noted that the behaviour of SM matter in such collisions, which does undergo significant momentum exchange due to plasma effects [52], and DM simulations such as [19], suggest that the effects may be large enough to be visible. However, establishing this carefully would require significant extra work.
For the subhalo case, even if only a small fraction of the stream kinetic energy and/or momentum is transferred to the subhalo during one crossing, then the fact that the subhalo's lifetime of billions of years is much longer than the crossing time (of only a few million years) means that there could be large cumulative effects. Momentum transfer would lead to an effective drag force, causing the subhalo's orbit to decay, while kinetic energy transfer would heat up the subhalo, expanding or eventually evaporating it. Consequently, it may be the case that effects on subhalos provide more robust signals of DM plasma instabilities (though again, further work would be required to understand these).
It is possible that, if structure formation does result in 'cuspy' DM profiles with significantly higher densities than we have considered here, then plasma instabilities could be important at even smaller couplings. Relatedly, it is possible that the first situations in which self-interactions are important could occur earlier in time than our example scenarios, and so change the initial conditions that we have been assuming. If one is attempting to identify the observational signatures of a model, then such issues would be need to be investigated carefully. Figure 7 represents a conservative estimate of the parameter space in which such investigations may be warranted (in the sense that the true space may be larger).
An important related point is that, for couplings significantly larger than the threshold value where effects first become significant, the DM's behaviour may be very complicated. In extreme situations, it is conceivable that for large enough couplings, DM may look effectively collisionless -for example, if early cooling processes are sufficiently strong that most DM is in isolated, dense, overall-neutral objects. Something similar happens in the SM sector, where some fraction of matter is in the form of stars, which behave as collisionless particles. While this possibility seems somewhat far-fetched for the simple model of symmetric DM we have been considering, in more complicated hidden sectors, DM particles may form small neutral bound states. These 'mirror atoms', in analogy to electron-nucleon bound states in the SM, occur in a variety of models (e.g. [5,[53][54][55][56]), and if the ionised fraction is small enough, would not display the plasma behaviour we have been investigating. Of course, they would also have different behaviour in other settings, e.g. different scattering properties with the SM in direct detection experiments, so they do not represent a way to open up the same parameter space.

Millicharged DM
So far, we have considered entirely hidden-sector dynamics, with DM-SM interactions occurring only through gravity. However, in some models, DM interacts directly with SM states. An important example is 'millicharged' DM, which has a (small) charge under SM electromagnetism. In this case, the interaction dynamics are literally those of a usual EM plasma, with the addition of a different particle species. Millicharged particles can arise in a number of ways [57,58] (it should be noted that models incorporating an extra hidden-sector photon have different behaviour, due to this extra mediator, as we discuss in Section 5).
Coherent scattering effects for millicharged DM, in the form of DM particles scattering off large-scale coherent EM fields, have been investigated extensively in the literature [20,[59][60][61][62][63]. These can be viewed as 'one-many' scatterings, with a single DM particle scattering off a field sourced by many SM particles, as opposed to the 'many-many' scatterings we have been considering. For example, [20] analyses scattering of millicharged DM with coherent magnetic fields in the ISM, finding that if q DM /m DM 10 −13 GeV −1 , the galactic disk would have been spun down by interactions with the DM. [61] analyses the interactions of DM with magnetic fields in galaxy clusters, and finds a similar coupling limit, above which the DM profile would be strongly affected. While the plasma dynamics of these scenarios are complicated, and such limits should be viewed as plausible estimates rather than established constraints, they illustrate that, even for rather small couplings, one-many effects may be important.
An obvious question is whether, for even smaller couplings, 'many-many' scatterings of the kinds we have been considering could be important. For example, from Figure 7, for a massless hidden-sector mediator, and m χ = GeV, the electrostatic instability may be important in cluster collisions for g few × 10 −16 . In a purely-DM environment, such effects would be of interest. However, in most settings, the SM matter cannot be neglected. In particular, since the DM charge must be ≪ 1 (for m χ 10 13 GeV), the relatively large electron charge (e 0.3) means that SM matter will usually dominate the plasma dynamics, even when the DM number density is higher. As we discuss below, this can suppress plasma instabilities (or makes them dominantly SM-sector, so that DM scatters against SM-sourced fields as per the 'one-many' scenario considered above).
To take an example, we can consider a beam-plasma scenario in which the high-velocity 'beam' is DM-dominated, and the denser 'plasma' is SM-dominated. This might arise when e.g. a bound structure with DM and SM matter is passing through a DM halo. From section 2.2.1, the electrostatic instability rate is ω i ∼ Rω p ∼ q 2 nχ mχωp , where q is the DM millicharge (the hot-beam limit is appropriate, as R is very small). Consequently, the instability rate is ∝ q 2 , rather than ∝ g as in the purely hidden-sector case. 9 Physically, the SM matter makes the plasma 'stiffer', suppressing the instability. The instability rate being suppressed by two powers of the small DM charge means that one-many scattering probes, for which the momentum transfer rate is ∝ q, are generally more powerful.
On the other hand, there may be other situations in which collective DM effects are important. An obvious example would would be systems in which the SM density is small enough. Another possibility is that coherent SM fields provide a starting configuration for which collective instabilities involving the DM can arise. [21] investigates DM instabilities in supernova shock waves, finding that if the SM plasma is strongly magnetised, collective instabilities involving the DM can be important. In order for this to happen, the effect of the initial SM magnetic field on the DM trajectories must be significant, and further investigation would be required to understand whether such instabilities can occur in parameter space regions which are not already constrained by one-many scattering processes.

Kinetically mixed mediator
Another important class of models are those in which the hidden-sector mediator interacts with the SM sector. Of these, the simplest have the new vector interacting with a conserved SM current; interactions with a non-conserved current lead to (energy/vector mass) 2enhanced rates for production of the vector's longitudinal mode [64], resulting in more complicated behaviour, and strong constraints on light vectors [65][66][67]. The only conserved currents in the SM are B − L (if neutrinos are Dirac) and EM. Couplings to B − L result in a long-range fifth force between bulk matter, so we restrict ourselves to an EM coupling.
If χ does not couple directly to the massless SM photon (i.e. it does not have a millicharge), then the Lagrangian for a 'dark photon' mediator is This is equivalent to the standard 'kinetic mixing' interaction through the field redefinition A = A + A , giving (to leading order in 1) (note the corrections simply modify the normalisations of these terms, rather than introducing new ones). Such an interaction can arise through integrating out heavier matter that couples to both A andÂ [68]. Another useful field redefinition is to define the 'active'/'sterile' basis, withÃ = A+ A the 'active' field, andÃ = A − A the 'sterile' field (again, we assume that 1, and ignore normalisation changes). Then, If m is small compared to the appropriate scales of interest (e.g. the SM plasma mass, as we will discuss), then the m 2 mass-mixing term is small, and SM matter interacts only withÃ. This illustrates why constraints on a dark photon from SM processes decouple as m 0 (these constraints are shown in Figure 8). Kinetically mixed mediators are of considerable theoretical and phenomenological interest [69]. The kinetic mixing interaction can arise from arbitrarily high scales [68], and is generically present whenever there is charged matter than interacts with the hidden sector and SM. Phenomenologically, the effectively medium-dependent coupling of the dark photon to SM matter enables models to access parameter space that would otherwise be constrained [70]. This latter point will be important for models targeted by low-threshold DM direct detection experiments, as we discuss in section 5.

Existing constraints
Some papers in the literature (e.g. [20]) claim that the 'one-many' bounds from DM scattering off coherent SM fields, as discussed in section 4, apply in the case of a 'sufficiently Figure 8. Constraints on a dark photon of mass m A , possessing a kinetic mixing with the SM photon. These constraints do not assume anything about the coupling of DM to the dark photon, relying only on SM processes (astrophysical constraints are plotted in blue, and laboratory constraints in red). The constraints come from Jupiter's magnetic field [71], CMB observations [72], laboratory tests of Coulomb's law [73,74], the CROWS microwave cavity experiment [75], and stellar cooling bounds [76]. The red dots correspond to the dark photon masses illustrated in Figure 10.
light' massive dark photon. However, unless the DM interactions directly with the massless SM photon, the threshold for 'sufficiently light' is very small. If the coupling is entirely through the kinetic mixing term, then we can see from equation 5.1 that χ interacts only with the massive A field, which decays exponentially at distances m −1 from sources. Consequently, if an SM magnetic field is sourced by currents on a scale L m −1 , then the corresponding B field is suppressed by ∼ (mL) −1 compared to B (since it is dominantly sourced by currents within a range ∼ m −1 ). For the galactic magnetic fields considered in [20], which have coherence length ∼ kpc, this suppression is ∼ 10 −11 10 −15 eV m , where the nominal m = 10 −15 eV value is small enough to alleviate almost all constraints on from SM processes (see Figure 8). This illustrates that 'small' masses by particle physics standards can nevertheless be high enough to very significantly alleviate millicharge-type constraints.
The 2 → 2 DM-DM scattering constraints for a dark photon mediator can be estimated from modified versions of the Coulomb scattering calculations in [13,24]. Since we are most interested in these constraints for m ω p /v th , where plasma instabilities are suppressed, the mediator's Compton wavelength will be smaller than the Debye length, λ D ∼ v th /ω p . This means that the maximum impact parameter for Coulomb collisions is set by m −1 , rather than λ −1 D , so the Coulomb logarithm is reduced compared to the massless mediator case. 10 There are a number of disagreements between [13] and [24] regarding the details of the Coulomb scattering calculation; we will not attempt to resolve these issues, but will use the (more conservative) expressions from [24], with a modified Coulomb logarithm. This gives a relaxation time, through Coulomb collisions, of To avoid making DM too collisional, we need τ iso few Gyr in galaxies [24], which gives a bound where we have taken a typical value of the Coulomb logarithm (Λ C will depend on g, m χ and m A , but these effects will not be large). As mentioned in Section 3, it would require further analysis to determine whether all couplings above this bound are constrained, since strong enough self-interactions could lead to very complicated behaviour. However, in the absence of bound state formation (which, in the symmetric DM model we are considering, would lead to annihilation), it seems likely that frequent collisions would not be compatible with observations.

Plasma instabilities
In a situation with no charged SM matter, the growth of plasma instabilties would be exactly as in the massive-mediator case treated in Section 2. To understand how the presence of SM matter affects this, we can analyse the coupled equations of motion for Fourier modes of the vector fields, Here, Π is the self-energy contribution from charged SM matter, while Π χ is from the DM (for electrostatic oscillations, we can reduce this six-dimensional system to a twodimensional one, so that the equations below become scalar equations). For a null eigenvector (A, A ), and writing K 2 ≡ ω 2 − k 2 , we have (K 2 − Π)A = ΠA ; using this in the A equation, we have where K 2 ≡ ω 2 − k 2 . In the DM-only case, we have ( λD (in the case of a massless mediator), if n −1/3 < λD. This does not agree with the results of standard scattering integrals (e.g. [77]), or analyses of electron-positron plasmas [78]. Conversely, [13] does not take into account collective effects at all, and takes the maximum impact parameter to be set by the Galactic radius. For a massless mediator, the Coulomb logarithm will be set by λD, between these two results.  Figure 9. Longitudinal instability rate for the two-stream velocity distribution from Figure 1, for DM interacting via a light mediator that is kinetically mixed with the SM photon (Section 5). The SM plasma distribution is taken to be cold (for simplicity), with a plasma frequency 100 times that of the DM plasma, and a collision rate of ∼ 10ω p (where ω p is the plasma frequency of the DM plasma). The different curves show the instability rate for different kinetic mixings. As these illustrate, the effect of the SM plasma is small, despite its high density and large collision rate.
If the SM charge density is appreciable, and the DM-mediator coupling is g 1, then the SM plasma self-energy may be much larger than the DM contribution Π χ . If our solution is a small perturbation of the DM-only case, so K 2 A (m 2 + Π χ )A , then Π being much larger than m 2 + Π χ allows us to expand the (K 2 − Π) −1 term in equation 5.8 by treating K 2 as small, giving The physical intuition behind this is that the SM matter is 'stiff', due to its strong coupling to the photon, and so decouples from the hidden sector dynamics.
In the opposite limit, where Π can be treated as small compared to K 2 in equation 5.8, we obtain Consequently, if Π is small compared to m 2 + Π χ , then we again obtain behaviour close to the DM-only case. Thus, in both of these scenarios, there are modes whose behaviour is close to the DM-only case. Even when the SM and DM plasma effects are comparable, we can solve for the behaviour of the combined system (equation 5.6). The upshot of this analysis is that, in most scenarios, DM-DM plasma instabilities are important in the same parts of parameter space as they would be if the mediator was not kinetically mixed. Figure 9 demonstrates this in a particular situation, by showing the instability rate for a two-stream DM velocity distribution, in the presence of a denser, collisional SM plasma. This illustrates that, even for fairly large , the main effect is to slightly rescale the singlemediator instability behaviour, as per equation 5.9. 10 -8 Figure 10. Constraints on the effective charge q eff ≡ g of a dark matter particle χ with mass m χ , coupling to a dark photon of mass m A with coupling g, where the dark photon has a kinetic mixing with the SM photon. The dark blue shaded region corresponds to existing constraints from stellar cooling [79], from SN1987A [80], and from the Xenon10 experiment [81,82]. The 'freeze-in' curve corresponds to the q eff required for early-universe production from the SM radiation bath to produce the correct DM abundance [83,84]. The SENSEI [85,86] and GaAs [87,88] lines correspond to projected sensitivites for future dark matter direct detection experiments. The light blue shaded regions correspond to parameter space in which 2 → 2 Coulomb scattering would have a significant effect in galactic halos (Section 5.1). The orange regions correspond to parameter space in which dark-sector plasma instabilities would undergo significant exponential growth in astrophysical systems, as discussed in Sections 3 and 5. The different panels correspond to different dark photon masses; for each dark photon mass, the largest kinetic mixing parameter consistent with SM constraints is assumed (illustrated in Figure 8), allowing the smallest g for a given q eff . We can see that, while the smaller g permitted by smaller m A relaxes the 2 → 2 scattering constraints, the lighter mediator mass promotes plasma instabilities.

Dark matter parameter space
For a direct detection experiment, if m ω p in the target material (or m L −1 , where L is the scale of the conductive shielding around the experiment, etc), then the interactions of a DM particle with the target can be treated as occurring via the 'active'Ã state (equation 5.3). Consequently, the DM scatters as if it had an effective charge q eff = g. Processes such as the production ofχχ pairs in stars, and in the early universe, also scale with g, so the effective DM parameter space, for given m A , is the (m χ , q eff ) plane, illustrated in the plots of Figure 10.
For a given scattering cross section (set by q eff ), we can decrease the required g by increasing . A major difference between a kinetically mixed mediator and others (such as B − L mediator, or a scalar), is that we can increase the SM coupling all the way up to O(1) by taking m A small enough; as discussed above, the constraints on from SM processes decouple as m 0. These constraints are shown in Figure 8; for m A 10 −16 eV, ∼ O(1) is allowed. Consequently, the 2 → 2 DM-DM scattering constraints discussed in Section 5.1, which are only weakly dependent on m A , can be significantly relaxed, to the point where they do not constrain parameter space. This is illustrated quantitatively in Figure 10, where the different panels shows the effect of decreasing m A and increasing . As noted in Section 5.1, while the m A required are small on particle physics scales, they are still significantly larger than the inverse length scales of astrophysical magnetic fields, so such models evade millicharge-type constraints. For other light mediators, strong constraints on the SM coupling from e.g. stellar cooling, or fifth force tests, significantly reduce the maximum allowed scattering cross-section in laboratory experiments. This generally results in near-term experiments having no realistic prospects of detecting such DM models, at least for DM masses 10 MeV [70]. However, since we must go to rather small m A to open up this parameter space, this raises the possibility that coherent DM-DM scattering processes become important. Figure 10 also shows the parameter space regions in which plasma instabilities will start growing in astrophysical systems. Unlike the 2 → 2 scattering bounds, the enhancement of instabilities due to decreased m A compensates for the increased that is possible, and much of the low-mass target range is within the plasma instability region, even for ∼ O(1). In particular, the 'freeze-in' cross section [83,84], for which thermal production of χ from the SM radiation bath is sufficient to produce the whole DM abundance, is within this regime up to O(10) MeV. Conversely, when m A is large enough that plasma instabilities would be unimportant, then is constrained to be small, and the Coulomb scattering bounds constrain much of the target parameter space for proposed light DM detection experiments.
As we discussed in section 3, these plasma instability regions should not be viewed as observationally constrained, since it is logically possible that the growth of perturbations does not lead to observable effects. However, our calculations do motivate more careful investigations into whether these regions of parameter space are astrophysically viable. Of course, even if it turns out that such particles cannot be all of the dark matter, it is possible that O(10%) of the DM mass density is strongly collisional [70]. To predict the signals of such subcomponent models at DM detection experiments, the effects of this collisionality on the DM distribution at Earth would need to be taken into account.

Light bosonic DM
In our discussions so far, we have assumed that dark matter can be treated as point particles, with some interparticle separation ∼ n −1/3 . However, when the de Broglie wavelength of of the DM becomes comparable to n −1/3 , i.e. the occupation number becomes 1, this picture breaks down. For DM in the vicinity of the Earth, this occurs at m DM eV. While most experiments searching for DM scattering do not cover such small masses, some types of experiments [89] and astrophysical observations [90] may have sensitivity there, and in any case, understanding the behaviour is of theoretical interest.
In the m DM eV regime, we can treat the (necessarily bosonic) DM as a large occupation number classical field. If, before structure formation, the particle and antiparticle fields have some independent perturbations 11 , then (ignoring non-gravitational interactions) the virialization process will result in them becoming effectively independent Gaussian random fields in halos, with characteristic wavelength λ B ∼ 1/(m DM v). If we consider a particular spatial point, then the amplitude of the particle field will generically be O(1) different from that of the antiparticle field. Supposing that the amplitude of the particle field is larger, then there is a net positive charge density over a spatial scale ∼ λ B . In charge density terms, the particle and antiparticle fields result in an array of charged 'blobs', each of typical size ∼ λ B , and typical charge ∼ ±n DM λ 3 B . If λ B is significantly smaller than the wavelength of perturbations that would grow in a uniform plasma, then the smooth approximations that we used above will be valid (c.f. the discussion around equation 2.12). If not, then a modified treatment would be required. For a two-stream velocity distribution, and a light mediator, the wavenumber of the most unstable perturbations is k ∼ ω p /v, so at the threshold coupling from equation 3.2, and taking the nominal cluster collision parameters from Section 3, Bearing in mind that the DM mass must be 10 −21 eV to be compatible with small scale structure observations [91], we see that for a light enough mediator, significantly abovethreshold couplings would be required to affect the plasma instability behaviour. For a heavier mediator, the perturbation wavenumber near the coupling threshold is k ∼ m A , so we need m DM m A /v there for the smooth plasma description to be accurate. For a heavy enough mediator and small enough coupling, plasma instabilities will be suppressed. However, 'collisional' effects on the plasma evolution will still exist, and will no longer be well modelled as 2 → 2 Coulomb collisions. Instead, 'blob-blob' collisions will occur coherently (this is in direct analogy to the picture of gravitational relaxation for light bosonic DM [91][92][93]). Treating scattering off a blob via standard scattering theory [94], the scattering rate (for a single DM particle) will be parametrically set by where N B ∼ n DM λ 3 B is the number of particles in a blob, n B ∼ λ −3 B is the number density of blobs, and log Λ C indicates the appropriate Coulomb logarithm. This form will be valid if 11 for example, as would arise from inflationary isocurvature fluctuations, or from small-scale fluctuations arising from a phase transition production mechanism. We leave the question of how large such fluctuations need to be for future work. m A m DM v, so that the mediator is long-range compared to the blob size. 12 Equation 6.2 has the same density and velocity dependence as for gravitational relaxation of light bosonic DM, but differs from the 2 → 2 Coulomb collisions rate (equation 5.4), which may result in stronger constraints than the latter at small m DM .
As an example, we can consider parameters motivated by the experiment proposed in [89]. If we take m DM = 10 −4 eV, so that the de Broglie wavelength is not larger than the detector, and m A = 10 −7 eV, giving a mediator range on the same scale, then for g 5 × 10 −14 , astrophysical plasma instabilities are suppressed (according to our estimates in Section 3). If we naively applied the 2 → 2 Coulomb constraints from Section 5.1, these would give g 6 × 10 −13 . However, if we take the estimate of the coherent scattering rate from equation 6.2, and set Γ ∼ 1/Gyr as a rough constraint, this would give Consequently, the coherence-enhanced scattering rate is more significant than plasma instabilities in this regime, and potentially more constraining than a naive application of the 2 → 2 Coulomb scattering bounds.
Despite these rates, it is unclear whether coherent scattering processes of these kinds will provide observational constraints. One issue is that, if the scatterings drive the particle and antiparticle fields towards identical spatial profiles, then we will no longer have charged blobs. Consequently, it is not obvious whether the scattering process will saturate before it has observable effects on halo shapes.
For light enough mediators (m A m DM v 2 ), it would also be necessary to check whether coherent emission of the mediator is ever important. While particle-antiparticle annihilation to the mediator is possible when 2m DM ≥ m A , the rate of this process is ∝ g 4 , so would likely be unimportant (there will be a rate ∝ g 2 from 'off-shell' parts of the DM field's wavefunction, due to it being bound in a gravitational potential, but these will be extremely small in a typical halo).
These caveats illustrate that more investigation would be needed to understand the observational consequences of light bosonic DM interacting via a light vector mediator. However, our estimates in this section do suggest that coherent effects may be important in new regions of parameter space. 12 If the mediator is short-range compared to the de Broglie wavelength, then the behaviour will be as per a contact interaction, giving Γ ∼ 1 16π (the constant factor here should be not taken seriously). Similar parametric behaviour would occur for e.g. a scalar of mass m interacting via a quartic interaction λϕ 4 , which would give Γ ∼ 1 16π where we have used the usual scaling of the quartic coupling λ ∼ m 2 /f 2 in terms of the symmetry breaking scale f , for a potential arising from high-scale breaking of a shift symmetry. This matches the parametric relaxation timescale from [95,96].

Conclusions
As we have emphasised throughout, the goal of our paper has not been to constrain models of self-interacting DM, but to map out parameter space regions in which long-range coherent self-interactions may be important. The obvious follow-on would be to more carefully investigate the physics of these models, and to determine whether they are observationally constrained. One of the most important models is the kinetically mixed mediator scenario, for which it would be very interesting to determine whether the models being targeted by proposed DM direct detection experiments [85,87,89,[97][98][99] are astrophysically viable. Even for purely hidden-sector interactions, such investigations may be quite complicated. In principle, one would like to start from the very early universe (e.g. from when the DM was produced) and track the DM's evolution until today. Without this, one would not know the initial DM distribution in scenarios such as cluster mergers, so unless one could demonstrate that all observationally-acceptable initial conditions behave in a certain way, relating the model to observations would be difficult. Of course, in practice, it would hopefully be possible to make significant approximations, e.g. that self-interactions have negligible effects before a certain point, or conversely, that they effectively thermalise the distribution. SM plasma behaviour may provide some guide, at least in the regime of a light mediator and strong coupling.
As well as understanding the constraints that could be obtained from plasma instabilities, it would be useful to understand the observational signatures that could arise, in models with a coupling that is just large enough. In particular, understanding the relationships between observational signatures in different astrophysical systems would be important. For DM that is not super-heavy, such threshold-level couplings would need to be extremely small -however, there are potential theoretical reasons why such couplings might arise [100].
There are many dark matter models which feature long-range vector mediators, beyond the simple examples we have considered. The most obvious examples are models with additional species, having different charges and/or masses [5,23] (such as the proton and electron in the SM). If most of the present-day dark matter is ionized (i.e. consists of isolated charged particles), then analogues of the plasma calculations we have carried out here will apply, though the presence of different-mass species will change the behaviour in potentially-interesting ways (c.f. ion effects in SM plasmas). On the other hand, if dark matter is dominantly in the form of overall-neutral bound states ('dark atoms'), then the interactions between these will be effectively short-ranged, and will not give rise to plasma instabilities.
As mentioned in Section 3, it may be observationally viable for a ∼ O(10%) subcomponent of DM to be strongly self-interacting [70]. At the threshold of viability, there may be interesting observational effects of long-range interactions in such a subcomponent. In section 1.2, we briefly reviewed [18,19,101], which simulated the dynamics of a DM subcomponent interacting via a massless hidden-sector mediator. Other questions, including the DM distribution at Earth of a subcomponent [26], would also be interesting to consider, and relevant to the direct detection phenomenology of such a subcomponent.  Figure 11. Normalised longitudinal self-energy for a two-Gaussian-stream velocity distribution, evaluated at zero phase velocity (see equation B.4). As discussed in appendix B, this function is related to the existence of unstable longitudinal modes. that mode will be unstable to growth. For a massless mediator, it can be shown that the existence of such a growing mode is equivalent to there existing some real u for which G(u) is real and ≥ 0 (effectively, if there are growing modes, then there must be a marginally stable mode). Since to have G(u) ≥ 0, we need g (u) = 0. For the example of a symmetric 'two-stream' velocity distribution, where each stream is Maxwellian with velocity distribution σ, and the closing velocity is 2v, we have where D is the Dawson integral [102]. This function is is plotted in Figure 11. For a massless mediator, an unstable mode exists if G(0 + i0 + ) ≥ 0, which is true for v 1.3σ. For a massive mediator, the dispersion relation is By the same logic as in the massless case, an instability exists if there is some real u for which G(u)−c 2 m 2 /ω 2 p is real and non-negative. From Figure 11, we can see that this will be true only for some bounded range of v/σ, if at all. If m 2 c 2 /ω 2 p ≥ max G(0 + i0 + ) 0.3/σ 2 , then there are no solutions, and consequently no instability. In terms of the notation in Section 2.2, where ω p labelled the plasma frequency for each stream, rather than for both streams, we need m 2 c 2 /(2ω 2 p ) ≥ G(0 + i0 + ). Less symmetric velocity distributions are also of interest. An example is the case of a dense, low-velocity-dispersion stream passing through a lower-density, higher-velocitydispersion stream -the 'beam-plasma' scenario discussed in Section 2.2.1. This kind of  Figure 12. Velocity distribution for a 'beam-plasma' scenario, in which the beam velocity dispersion is σ s = 1, the beam velocity is v s = 1, the plasma velocity dispersion is σ p = 0.1, and the density ratio is R = 0.1. The blue curve shows the velocity distribution of the 'plasma' stream, the orange curve that of the 'beam', and the dotted curve the combined velocity distribution.
velocity distribution is illustrated in Figure 12. If both of the streams have Maxwellian velocity distributions, and we take the denser stream to be at rest, so (where we have assumed that v s v, which is true if σ s σ p ). If the beam velocity distribution is broad, so v s ∼ σ s , then we have u ∼ Cσ p , where C is at most logarithmically large. The integral for G(u + i0 + ) is then dominated by the 'plasma' peak near the origin, which (via integration by parts) is Consequently, for a massive mediator, an instability exists if for some C ∼ O(10).