{'fraction_non_alphanumeric': 0.03518204305251753, 'fraction_numerical': 0.010896007522947032, 'mean_word_length': 4.366758091058695, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the stability of shocks with particle pressure

We perform a linear stability analysis for corrugations of a Newtonian shock, with particle pressure included, for an arbitrary diffusion coefficient. We study first the dispersion relation for homogeneous media, showing that, besides the conventional pressure waves and entropy/vorticity disturbances, two new perturbation modes exist, dominated by the particles' pressure and damped by diffusion. We show that, due to particle diffusion into the upstream region, the fluid will be perturbed also upstream: we treat these perturbation in the short wavelength (WKBJ) regime. We then show how to construct a corrugational mode for the shock itself, one, that is, where the shock executes free oscillations (possibly damped or growing) and sheds perturbations away from itself: this global mode requires the new modes. Then, using the perturbed Rankine-Hugoniot conditions, we show that this leads to the determination of the corrugational eigenfrequency. We solve numerically the equations for the eigenfrequency in the WKBJ regime for the models of Amato and Blasi (2005), showing that they are stable. We then discuss the differences between our treatment and previous work.


Introduction
One of the major uncertainties still surrounding particle acceleration around shocks is the level at which this saturates: i.e., the total energy content in non-thermal particles, per unit volume, R, in units of the incoming fluid kinetic energy density. This saturation level plays a very important role, for instance, in discussions of the origin of cosmic rays as observed at Earth: it has to be rather large ( 0.1) to allow SNe to provide the observed flux. Also, discussions on the origin of UHECRs are influenced by similar considerations: the idea that UHECRs are accelerated by GRBs (Vietri 1995;Waxman 1995) draws some measure of support by the coincidence between the energy release rate of GRBs in γ-ray photons and that required to account for observations at Earth (Waxman 2002;Vietri et al. 2003).
It is certainly possible that this level of saturation is somehow connected with the very uncertain particle injection mechanism, and that saturation occurs, in many astrophysically important situations, at very low levels. Yet, given the arguments concerning SNe and GRBs mentioned above, it appears that at least occasionally saturation must occur at rather large levels. This is especially so in the case of Galactic cosmic rays: if in fact we were to discard SNe as sources of cosmic rays, all other possible classes of sources, being both less numerous and less energetic, would force us to require saturation at even larger levels 2 .
A distinct possibility is that the saturation level is determined not by the injection mechanisms, but by modifications which the particles' pressure induces in the shock properties. After all, in the test particle limit the particles' spectrum is ultraviolet divergent (though marginally so, of course), and one may hope that, since the particles' back-reaction on the fluid makes acceleration less efficient (by reducing the velocity jump around the gas subshock), a convergent spectrum will be obtained, with a definite value for the parameter R. After Malkov's seminal papers (Malkov 1997;Malkov et al. 2000;Malkov and O'C. Drury 2001), fully self-consistent solutions with particles' pressure properly included have been obtained by Amato and Blasi (2005). In these models, it appears that particles can carry away an arbitrarily large fraction of the incoming energy flux, for sufficiently large Mach number at upstream infinity. Even more intriguing is the fact that these solutions have particles' spectra which are more, not less, ultraviolet divergent than those in the test particle limit. In fact, this divergence is arrested only by arbitrarily limiting the largest individual energy to some fiducial value (Amato and Blasi 2005).
These solutions then seem to beckon for a stability analysis, in the hope that they are found to be unstable once R exceeds some critical value. The instability we envision is of course the corrugational instability for shocks, whereby ripples on the shock surface become of larger and larger amplitude 3 . If this instability were to exist, we might easily imagine that the shock is substituted by a region of strong fluid turbulence, where particles moving a few Larmor radii perceive only a small velocity difference at the two ends of their free wandering, and thus acceleration to high energies is made somewhat less efficient.
However, it is well-known (Landau and Lifshitz 1987) that shocks in polytropic fluids are stable against this kind of perturbations, and against the spontaneous emission of sound waves as well. The hope for the existence of an instability is based on a rather more subtle argument. When the shock surface flaps, it sheds in the downstream region pressure waves as well as entropy and vorticity perturbations. We will show later that the last two do not couple to particle perturbations, but pressure waves do, thus generating small perturbations in the particles' distribution function, δf . These particles will however return upstream by means of diffusion, and here generate, by their sheer perturbed pressure, some more fluid perturbations. Thus, not all energy shed by the shock is lost forever toward upstream infinity, as is the case in the purely hydrodynamical case, but some fraction of it returns to the shock to generate more flapping. In fact, since pressure waves are strongly damped by diffusion, and nearly all particles return to the upstream section since the shock is Newtonian, we may conjecture that most energy shed by the shock flapping makes it back to feed more flapping, even though account must be taken of diffusion and phase mixing. Is it possible that this sets up some strong reinforcement, making the whole system unstable? This is the question we address in this paper.
This question has of course been studied before, under somewhat different conditions. The diffusion coefficient of non-thermal particles by the fluid (≡ D) has been taken (by other authors) as independent of the particle impulse: this is conventionally referred to as the two-fluid approximation, because the Boltzmann equation for the non-thermal particles can then easily be recast into an equation for their pressure, thusly erasing all references to the underlying distribution function. Besides making the two-fluid approximation, Mond and O'C. Drury (1998) also neglected diffusion altogether. A more complex analysis has been presented by Toptygin (1999), who included energy transport and particle injection at the shock, but still in the two-fluid approximation. In the following, we shall make no such approximations: we will consider a finite diffusion coefficient D, and it will be allowed to be an arbitrary function of both p, the particle momentum, and of ρ, the local fluid density, so as to at least mimic the increase in magnetic field strength due to flux freezing. The arbitrary nature of the dependence of D on p automatically prevents the use of the two-fluid approximation, and forces us to use the full Boltzmann-Skilling equation.
A word of caution is in order about some assumptions. We will neglect all energy in the form of magnetic field and Alfvén waves; of course this is necessary because the zero-th order solutions of Blasi and collaborators do not include these effects, but in our case this neglect requires one extra assumption, i.e. that the time scales for energy to accumulate into any of these energy sinks, T in , and to flow out of each of them, T out , be ordered like this: T out T in .
If this inequality were severely violated, the magnetic field or Alfvén waves might acquire a significant fraction of the total energy, despite their negligibility in the zero-th order solution, and make our treatment completely irrelevant.
In the following, we shall follow closely the treatment of the shock corrugational instability given by Landau and Lifshitz (1987). In particular, we shall consider a steady-state shock in its own frame, located at z = 0, with fluid coming from the left, and exiting from the right, so that all speeds are > 0. Exactly like Landau and Lifshitz, we shall consider perturbations generated by the shock flapping, so that there can be no incoming waves, from either upstream or downstream infinity. In Section 2, we shall consider perturbations in a homogeneous medium; this is perhaps a tad boring, but it contains a number of new results which are of the utmost importance later on. In Section 3 we briefly discuss perturbations upstream, i.e. where the flow is inhomogeneous; we show here that we can easily obtain the perturbations in the WKBJ limit k y → ∞, which restricts our analysis to the regime λ L, where λ is the perturbation wavelength perpendicular to the shock, and L is the typical size of the region of inhomogeneity upstream. In Section 4 we discuss the boundary conditions on the perturbed particle distribution function, and derive how it relates to the amplitude of the modes to which particle perturbations couple. We present in Section 5 the perturbed Rankine-Hugoniot conditions and in Section 6 what fixes the global corrugation mode eigenfrequency. In Section 7 we present our numerical computations for the stability of the exact solutions of the zero-th order problem by Amato and Blasi (2005). In Section 8, we will compare our results with other works in the literature, and briefly summarize our work.

Perturbations in a homogeneous medium
We give here, for future reference, our basic equations. They are the conventional hydrodynamic equations: which contain a term for the momentum exchange between the fluid and the non-thermal particles represented by the gradient of the particle pressure P c , plus the usual Boltzmann equation in Skilling (1975): We assume D = D(p, ρ) to be a given function of ρ and p.
We consider small-amplitude perturbations around a homogeneous solution where the particles are supposed to exert a non-negligible pressure. First, we consider entropy perturbations. Perturbations can be taken in the form so to obtain, from the equation of entropy conservation, eq. 3, where u is the unperturbed fluid velocity in the x-direction. Perturbation of the mass conservation equation yields Perturbation of the momentum equation yields Now, the perturbation of the Boltzmann equation yields which implies δf = 0 because, for entropy perturbations, ω = uk x and k · δ v = 0, as deduced above. Furthermore, since δf = 0, necessarily δP c = 0, and thus δP = 0. It follows that entropy perturbations and particle perturbations are completely decoupled: in fact, entropy perturbations are just advected by the zero-th order flow.
In summary, entropy perturbations have the following characteristics, where we define, for future convenience, V ≡ 1/ρ: where the last equation applies to ideal fluids: m is the average particle mass and k B Boltzmann's constant.
We turn now to isentropic perturbations, which automatically satisfy eq. 3. Mass conservation implies while momentum conservation implies Let us begin by assuming that the left-hand-side of eq.11 vanishes; the same is then true for ▽·δ v. Multiplying eq. 12 by k∧, we see that ▽∧δ v = 0, unless ω = uk x . If the curl vanishes, then so does δ v, because any vector with vanishing divergence and curl is a constant, which can always be set to zero by a suitable choice of reference system. So we must have ω = uk x for perturbations with non-zero vorticity; we thus find that δP + δP c = 0, and then, again from eq. 9, that δf = 0, implying δP = δP c = 0: vorticity perturbations do not couple to particles either. We have for them: We consider now those perturbations where δf does not vanish. We can show that here too there are two distinct classes of modes: in the first one, δf is not coupled to the fluid quantities, while in the second one it is through the term k · δ v = 0. The first class of modes, which we call d-mode, cannot be obtained directly from eq. 9, for reasons to be made clear shortly. We use instead the correct form where we dropped the term ∝ k · δ v in keeping with our desire to find a solution totally uncoupled from the fluid. The solution of this equation is known from elementary courses: if φ • (x, y, p) is the initial distribution function at time t = 0 (possibly dependent upon p), the solution at later times for u = 0 is and the solution for u = 0 is just δf (x − ut, y, p). At the same time, we must make sure that this solution does not ruffle the fluid: this obviously requires δP c = 0 at all times. Now integrating eq. 14 over 4πp 3 vdp/3 we find which clearly shows that, in order to have δP c = 0 everywhere at all times we need δP c = Dδf p 3 v dp = 0 everywhere at the initial time. If we now Fourier-expand the initial condition φ • with respect to x, y, the above conditions become φ • (x, y, p) = dk x dk y g(p, k)e −ıkxx e −ıkyy , g(p, k)p 3 v dp = g(p, k)D(p)p 3 v dp = 0 .
(17) This completes the derivation of this purely damped d-mode, which will not perturb the fluid. Though it may look at this stage like a mathematical oddity, this mode plays a key role in the matching of boundary conditions between the upstream and downstream regions. It is also worth remarking why it cannot be derived from its Fourier-analyzed counterpart: the solution in question contains a term ∝ e −1/t , which does not have a Fourier transform with respect to t.
When δP c = 0, velocity perturbations δ v = 0 are induced in the fluid by the nonvanishing particles' pressure gradient, and eq. 9 can be solved for δf as If we integrate this over 4πvp 3 dp/3, we find where we assumed the non-thermal particles to have a minimum (p m ) and a maximum (p M ) momentum. It is convenient to introduce a weighted diffusion coefficientD defined as follows: which showsD to be a function of z only.
The above can be simplified a bit by integration by parts. If the integral were to extend from 0 to +∞, we would have − ∞ 0 4π 9 dp vp 4 ∂f ∂p = ∞ 0 4π 9 dp ( d dp where P c is the particle pressure in the unperturbed flow, and γ c , the effective particle polytropic index, satisfies 4/3 ≤ γ c ≤ 5/3. Since however the integral extends only from p m to p M , we define Thus eq. 19 can now be rewritten as: We remark that both γ ′ c andD are obtained by weighing the zero-th order solution, so that they can be immediately computed as soon as this solution is available. We can now eliminate k · δ v between eqs. 11 and 12, and then use δP = c 2 s δρ (with c s the sound speed because we are considering isentropic perturbations) and the equation above to obtain the comoving eigenfrequency, in suitably scaled units.
Eq. 24, together with the definition ofD (eq. 20), is the sought-after dispersion relation for the coupled small perturbations in a homogeneous medium.

Properties of coupled modes
The case where D is independent of p has been derived before (Ptuskin 1981) and coincides with the above equation. Surprisingly, the existence of this mode was not noticed by Toptygin (1999), even though a careful treatment of his equations yields precisely the same dispersion relation as above; this neglect of this mode has important consequences to be discussed later on.
It is best to begin our discussion with the case when D, and thusD, is a constant, independent of p. The eq. 24 then reduces to which is a simple polynomial equation of the third order. In this case two modes reduce to pressure waves, as is most easily seen in the test-particle regime P c = 0. There is however a third solution which is only slightly more mysterious: in the test particle regime these modes represent a local over-density of particles dissipated by diffusion. When the test particle regime does not apply, a particle contribution to the sound speed is introduced by the term ∝ P c . This new mode (which we call the third mode) is the equivalent of the d-mode when however the conditions 17 are not satisfied: the basic idea is still the same, the particle overdensity is dissipated, but since the particle pressure does not vanish, the fluid is consequently ruffled. Notice also that there are two third modes, traveling in opposite directions.
The situation is slightly more complex when D = D(p), because one must solve simultaneously eqs. 24 and 20. We begin by considering the limit k → 0. In this case, and assuming Ω → a constant, we easily find, to leading order in k, Ω 2 = 1 + γ ′ c P c /(γP ): as it must, the dispersion relation allows pure pressure waves, with the particle pressure providing a correction to the (pure gas) sound speed, because for large perturbation wavelengths particles are entrained by the perturbation. We remark that, in this limit, pressure waves are supersonic, in the sense that they are faster than pressure waves propagating in a pure gas of the same thermodynamical state, a result already noticed by Toptygin (1999).
Assuming Ω → a constant, we lost a solution, so we search for the third mode solution in the form Ω = αk+ higher order terms in k. We find: Now we use the definition z ≡ ık 2 /(ω − uk x ) = ık/(Ωc s ) to simplify the above to which can be now used with eqs. 20 and 22 to obtain where of course ∂f /∂p < 0. As a function of real z, the right-hand side above (where it exists) is easily seen to be a monotonically decreasing function of z, vanishing for z → ±∞. In any realistic physical problem the integral must extend from a minimum (p m ) to a maximum momentum p M ; since D(p) is expected to be a monotonically increasing function of p, we see that the integral above always exists for z < 1/D(p M ) ≡ 1/D M , and z > 1/D(p m ) ≡ 1/D m , and it diverges exactly at z = 1/D M and z = 1/D m . Thus the integral on the right-hand side of the equation above spans the whole range from 0 to −∞ as z varies between −∞ and +1/D(p m ), and the range +∞ to 0 as z varies between 1/D m and +∞. Somewhere in the range 1/D m < z < +∞ there is the one and only solution of the above equation. An illustration of the integral on the right-hand side of the previous equation is shown in Fig.  1, for a specific distribution function from Amato and Blasi (2005): the qualitative features of this plot are generic to all distribution functions we have tried.
Furthermore, since z = ık 2 /(ω − uk x ), and the small perturbations were assumed to vary as e ı(ω−ukx)t , the result that z > 0 implies that all modes are damped by diffusion, as physical intuition obviously suggests.
The discussion in the opposite limit, k → +∞, is similar. If we assume Ω → a constant, we find the solution without the correction to the sound speed due to the presence of the particles' pressure: in the limit k → +∞ particles escape by free streaming, and do not contribute to the sound speed.
Again, we lost a solution, so we now search for the third mode as Ω = αk+ lower order terms in k. We obtain now: which can be rewritten as zD(z) = 1 .
Comparing this with eq. 20, we see that the value of z we are searching for is the one that makes the integral on the numerator of the right-hand side of eq. 20 diverge. Following the previous discussion, we see that this is This is always positive, so that the solution is always damped by diffusion. This result has a simple physical explanation: when a small overdensity of particles is generated locally, the time-scale for damping of this overdensity is dictated by diffusion of the slowest particles.
Again for illustrative purposes, the real and imaginary parts of sonic and third modes are displayed in Fig. 2 for a specific distribution function from Amato and Blasi (2005). Again, the qualitative features are generic to all models we investigated.
In summary, we have seen that the introduction of particles modifies the modes of a homogeneous medium by adding two new modes, one coupled and one uncoupled to the fluid, both strongly damped by diffusion. -Real and imaginary parts for the three solutions of eq. 24 as a function of perturbation wave-number k, for the same distribution function as in Fig.1 and γ ′ c P c /(γP ) ≈ 158.629 (upstream pressure at 0 − ). The pressure waves are seen to be damped with the same rate (and opposite oscillation frequencies), while the third mode is pure imaginary, corresponding to a pure damping.
For future reference, we give the expressions for all small quantities, in units of δP c : the link between δP c and δρ is obtained from Eqs. 20 and 22, and the others follow easily. For later convenience, we use again V ≡ 1/ρ rather than ρ: where the definitions of z and q will come handy later on. In these equations, k x must be regarded as a known function of ω and k y , specified by eq. 24.

Initial value problem
As a preparation for later work, we discuss the initial value problem. This has some interest because perturbations in δf belonging to the various modes are not mutually orthogonal, and thus it may appear that initial conditions, especially when given only in terms of δf , cannot be decomposed into mutually independent modes. To fix ideas, let us consider a homogeneous zero-th order solution where all fluid quantities are unperturbed, but a small perturbation δf 0 (x, y, p) at time t = 0 is given: what are the amplitudes of the four modes that will be excited (two pressure waves, the third mode, and the d-mode)? First of course we Fourier-analyze δf 0 , calling a(p) its amplitude. We must then have where we have introduced some notation that will be useful in the following: a summation convention over i is understood, the A i 's are the amplitudes of the pressure waves and third modes for i = 1, 2, 4, respectively; g(p) is the amplitude of the d-mode, as in eq. 17, and the δf i 's are Comparing this with eq. 18, it becomes clear that the mode amplitudes A i 's are simply δρ i /ρ. Here the quantities ω i , k xi , k 2 i are supposed to be linked by the appropriate branch of the dispersion relation. For g(p) to represent a proper d-mode, we know that it must satisfy the two constraints in eq. 17; thus we derive two conditions on the mode amplitudes: The last condition can be obtained by satisfying the requirement that the perturbations to both density and velocity vanish at the initial time. With reference to eq. 34, we see that the two conditions, the vanishing of the density and of the velocity, imply which are two more linear relations which, together with eqs. 37 and 38, determine the A i 's. This simple example illustrates the importance of the d-mode, a kind of elephants' graveyard because, once particles join it, the perturbations they generate can only be dissipated away without ruffling the fluid ever again.

Perturbations upstream
It is useful to remark that, upstream, entropy perturbations are not allowed for arbitrarily dishomogeneous flows, while the same is not true for vorticity perturbations.
In fact, the perturbation of the entropy equation is: where we could neglect the term δv x ds/dx because the fluid is isentropic in the unperturbed state. From the above, we see that entropy perturbations are advected by the background flow all the way from upstream infinity to the shock; we cannot however accept this, since in our problem all perturbations must have as a source the shock flapping. Perturbations riding all the way from upstream infinity belong to perturbations in the boundary conditions, and are thus irrelevant. Thus the perturbed fluid will be assumed adiabatic, from now on.
The same argument does not apply to vorticity perturbations when the fluid is stratified, because they do couple to particle perturbations: in fact we easily obtain from eq. 2 that the equation for the vorticity, η ≡ ∇ ∧ v, is: In the absence of particles, this equation tells us that vorticity is exactly (i.e., not just to zero or first order) Lie-advected by the flow because, for adiabatic flows, ∇ρ ∧ ∇P = 0. But in the presence of particles and of spatial gradients it is easily seen that the particle pressure P c acts as source of vorticity perturbations.
In order to make progress, we compute the spatial dependence of the various modes upstream in the WKBJ approximation: i.e., in the limit k y → ∞, and in particular k −1 y ≪ L, the size of the upstream precursor. In other words, we take all perturbed physical quantities to be of the form and all derivatives in the x direction (i.e., perpendicular to the shock) are considered small when compared with terms proportional to k y . This is a perturbation analysis in which the transverse wavenumber k y is assumed large and as a consequence the longitudinal wavenumber ≈ W X is also large. The presence of a non-constant amplitude Q X (x) is equivalent to keeping the first two terms in an asymptotic expansion in the small parameter (k y L) −1 . This is often called the physical optics approximation (Bender and Orszag 1978).
This analysis is quite standard, but the amusing thing is that we don't even need to carry it through. In fact, we shall show later that the stability analysis requires knowledge of the physical quantities immediately before the shock, while knowledge of the perturbations run with x further from the shock is immaterial. We see from the above that all physical quantities close to the shock satisfy This is precisely the same form that holds in the homogeneous medium, so that eqs. 34 and 17 still hold. Also, the space-time dependence of the d-mode for the half-plane x > 0 is derived in Appendix A, assuming spatial homogeneity of the background solution; and this, for the argument above, applies in the WKBJ limit also to the upstream fluid. Let us also remark that, in this spatially homogeneous limit, vorticity perturbations do not couple to particles: in other words, spatial gradients are negligibly small, and thus vorticity perturbations, which can only be Lie-advected from upstream infinity, must vanish identically.
This is the result we need: since we are using a WKBJ approximation, we can treat the upstream fluid as if it were homogeneous, with the values for the physical quantities taken to be those immediately before the shock: we call this the Homogeneous Approximation: it clearly breaks down when the WKBJ analysis does, which occurs for k −1 y ≈ L, the size of the upstream precursor.

Boundary conditions for the particle distribution function
Because of diffusion, particles are not restricted to the downstream part of the flow: there will be a δf = 0 also upstream, so that we need to consider appropriate conditions to match δf in the two regions. The first condition is obviously the continuity across the shock, It is also well-known that the spatial gradient of δf needs to satisfy a boundary condition at the shock: this is derived by integrating eq. 4 on an infinitesimal interval straddling the shock. The unit vector normal to the surface of the flapping shock isn = (1, ık y ζ) where ζ is the shock corrugation amplitude. So, we obtain: Writing its first-order linearization and using eq. 44 we find: We note that (∂f /∂y) 2 = (∂f /∂y) 1 = 0 because at zero-th order the medium is uniform in coordinates parallel to the shock surface. The result is: Eqs. 44 and 47 are the appropriate boundary conditions for our problem.
We now show how to satisfy the boundary conditions, eqs. 44 and 47 at the shock. Using the notation introduced in Sect. 2.2, we know that δf + on the downstream side of the shock (the factor e ıωt−ıkyy will be omitted for simplicity in this subsection) satisfies Please notice that both downstream and upstream the summation is over 2 modes (a pressure wave and a third mode). In fact, we know on the one hand that particles perturbations are not coupled to entropy and vorticity perturbations, while, on the other hand, we know that pressure and third modes, for given values of k y and ω, exist for opposite values of k x (see eq. 24), and thus at least one pressure mode and one third mode will diverge exponentially toward infinity: this is unacceptable because we are studying flow instabilities, not variations in the boundary conditions at infinity. This discussion will be completed in Section 6.
On the upstream side we have analogously The continuity of δf at the shock allows us to derive a relationship between the g's: Using into eq. 47 we find which we regard as an equation for g u (p), whose solution can be written as where the functions w's are derived in the Appendix B.
We can now impose the conditions 17 on g u , obtaining: C p 3 vw C dp + A iu p 3 vw iu dp + A id p 3 vw id dp + (δv 2x − δv 1x ) p 3 vw 0 (p)dp = 0 , (55) And now we can impose the very same conditions on g d , eq. 50, to obtain: This set of four linear equations in five unknowns (C and the A i 's) can be solved in terms of one of them, say A 1d , the pressure wave downstream.
We thus see that the conditions at the shock, plus knowledge of the modes, allows us to determine the amplitude of all modes (except the vorticity and entropy modes downstream) in terms of the amplitude of the pressure wave downstream. Remembering eq. 34, we now see that all fluid quantities at the shock upstream are determined in terms of the amplitude of this very same wave; as for the two d-modes, their explicit space-and time-dependence is not needed, but is given for completeness' sake in Appendix A. Now, before discussing how to fix ω, the shock eigenfrequency, we briefly summarize how to perturb the standard Rankine-Hugoniot conditions.

The fluid conditions at the shock
We follow closely the discussion in Landau and Lifshitz (1987, Section 90), which requires only a small adaptation to our problem. In their consideration of the corrugational instability, in fact, there were no perturbations on the upstream side of the shock: deviations from the unperturbed state were generated by the corrugation of the shock surface (and its motion), and led to non-zero perturbations only downstream. In our problem, instead, the particles crossing the shock manage to generate new perturbations on the upstream side, leading to a slight modification to the Rankine-Hugoniot conditions. From now on, we indicate with the subscript 1 the quantities on the upstream side of the shock, and with the subscript 2 those on the downstream side. Also, to make contact with Landau and Lifshitz' work easier, we use as a variable V ≡ 1/ρ rather than the density ρ directly.
We consider a small-amplitude corrugation of the shock surface, away from the x = 0 plane, by a small displacement of the form: with respect to which the unit vectors parallelt and normaln to the surface have components in the xy plane:t = (−ık y ζ, 1) ;n = (1, ık y ζ) while the surface speed in the direction normal to the surface, with the respect to the reference frame of the unperturbed shock, is: All quantities are, of course, accurate to first order only.
The first two Rankine-Hugoniot conditions to be perturbed involve the fluid speed, and they are the continuity of the fluid speed parallel to the shock surface, and the discontinuity of the perpendicular component in terms of perturbed pressure and density (eq. 85.7 of Landau and Lifshitz 1987). We have: whose first-order linearizations are: The next equation to be perturbed is the shock adiabat, which is convenient because it is independent of all speeds involved. For a polytropic gas like the one we are considering, the unperturbed Hugoniot adiabat is given by (Landau and Lifshitz 1987, eq. 89.1): where we have defined h ≡ 4γ ((γ + 1) + (γ − 1)P 2 /P 1 )((γ − 1)P 1 /P 2 + (γ + 1)) (66) while the ratio P 2 /P 1 can be expressed as with M 1 the Mach number of the upstream fluid.
We need one more equation: we can use the equation relating the mass flux to the discontinuities in pressure and density: eq. 85.6 of Landau and Lifshitz (1987) recites: where j is the mass flux. When the shock surface is perturbed, the above becomes The perturbation to first order of the above yields 2δv Lastly, the condition that δf be continuous across the shock has already been imposed, see eq. 44.

The equation for the perturbation eigenfrequency
We discuss here first how this problem differs from the one without particles, keeping in mind that we are interested in the stability of the shock flapping, so that all modes present must have this flapping as their source: we cannot allow modes to come in from spatial infinity, because this amounts to perturbation of the boundary conditions, not of the shock geometry.
When we can neglect the particles' pressure, we know that there can be no perturbations upstream, since they are either generated at upstream infinity (in which case we would not be treating the case of shock instability) or, if generated at the shock, they cannot propagate away from it fast enough (the shock is supersonic). In the presence of particles the situation changes because particles can diffuse back to the upstream region, so that a perturbation δf generated downstream can return to upstream and perturb the fluid quantities. The situation is even more remarkable when one notices that, in this way, there can be a generation of pressure waves (even though they are just sonic in a supersonic medium) in the upstream region. The reason is shown in eq. 2: the gradient in particles' pressure is a source of perturbations, and since particles scatter a finite (i.e., non-zero) distance from the shock, there is no obvious reason why even sonic perturbations should not be generated. The impossibility of having sonic perturbations in a supersonic medium arises only when the point of generation of the perturbations is the shock itself, not a finite distance from it.
We must also remark that the presence of particles shuffling between downstream and upstream, and viceversa, has important consequences also for the kind of waves present downstream. In fact, while in the absence of particles the only waves present can be those shed by the shock, i.e. entropy, vorticity and pressure disturbances propagating to downstream infinity, when particles are included in the picture we find, by complete analogy with the argument above for the upstream region, that they can seed the third mode and the d-mode.
We have seen in Section 4 that the amplitude of all modes (and thus of all physical quantities) upstream, and of the third mode downstream, can be expressed in terms of a single free parameter, the amplitude of pressure waves downstream. Thus all amplitudes are fixed, except for three modes downstream (entropy, vorticity, pressure) and the shock displacement ζ; but we also have the four perturbed Rankine-Hugoniot conditions, which can be expressed in terms of these four quantities. The vanishing of determinant, once substitution of amplitudes for the upstream modes, and for the third mode downstream by means of the relations in Section 4 has been made, thus fixes the eigenfrequency.
A word of caution is in order: we have tacitly assumed that the pressure wave propa-gating away from the shock in the downstream region is also the wave with the physically acceptable (i.e., decreasing) behaviour at downstream infinity. This is not certain: in fact, even in the standard case with no particles, a proper mode exists when the mode explodes exponentially in time and the pressure wave propagating away from the shock dies down at downstream infinity (Landau and Lifshitz 1987). This condition can only be checked a posteriori, i.e. after having found numerically the value of ω.
We now explicitly derive the set of equations whose determinant must vanish, which is what fixes ω for a given k y .
In the four Rankine-Hugoniot equations there appear the total perturbations of specific volume, pressure and velocity. They can be written as sum of the respective perturbation for each mode. Downstream we have: where the subscripts s, v, p and t label quantities for entropy, vorticity, pressure waves and third mode. Upstream there are neither entropy nor vorticity perturbations: We must use these expressions in the perturbed Rankine-Hugoniot. Furthermore, we can write one of the two components of δ v 2v in terms of the other one, through k 2v ·δ v 2v = 0, where k 2vy ≡ k y and k 2vx = ω/v 2 (see eqs. 13). Then we can simplify our system by eliminating the shock displacement ζ, the remaining component of the velocity perturbation due to vorticity, and δV 1s . One equation remains, linking only pressure waves and third modes' perturbations: +2 ω 2 k y P 1 (P 1 − P 2 ) P 2 v 2 V 1 (V 1 − V 2 ) (δv 2py + δv 2ty ) = 0 .
Now, eqs. 34 hold both for pressure waves and third modes. We can use them to write the above equation in terms of volume perturbations. We remark that each mode has its own z and its own k x but, for given k y and ω, they are fixed by the dispersion relation, eq. 24. We obtain: where the sum is over the pressure mode and the third mode. Now, recalling the definition of A i , we see that A i = −δV i /V (see also the discussion following eq. 36). Then we have the four equations 55, 56, 57, 58, plus eq. 78, for a total of five linear homogeneous equations in five unknowns: four A i and C. The system has non-trivial solution only if its determinant vanishes. This condition determines the eigenfrequency of the system.

Results
In this section we apply our stability analysis to two shock solutions: a single solution (Fig. 3, upper panel) by Amato and Blasi (2005) and a multiple solution (Fig. 3, lower panel) by Amato et al. (2008).
As we stated above, instability takes place when perturbations exist and grow exponentially with time. Furthermore, they must decay away from the shock. These conditions may be summarized as follows: where the subscript i indicates that these conditions must hold for all waves.
For illustrative purposes, first we consider the solution of D'yakov's equation (see Landau and Lifshitz 1987, eq. 90.10): which fixes the shock eigenfrequency in the linear regime, with We apply D'yakov's analysis to a test particle solution for a strong shock (upstream Mach number M 1 → ∞) in a polytropic fluid with index γ = 5/3. From eqs. 89.6 and 89.9 in Landau and Lifshitz (1987) we obtain, respectively, the compression factor R = V 1 /V 2 = 4 and the downstream Mach number 1/ √ 5. The downstream sound speed is c s2 = u 2 /M 2 = u 1 /RM 2 = u 1 √ 5/4. In such system sound waves can propagate away from the shock only downstream. The dispersion relation for such perturbations is straightforward: where k x is the x-component of the sound wave. To write down the eigenfrequency equation for corrugations of the shock we need to calculate h L . This task is straightforward for a strong shock because the upstream pressure vanishes. As a consequence, the derivative (∂V 2 /∂P 2 ) V 1 ,P 1 vanishes too (see eq. 64) and h L = 0. Eq. 80 becomes: Eqs. 82 and 83 form a system to be solved with respect to ω and k x for a given value of k y . Since these equations are third-degree homogeneous in ω, k y , k x , if (ω, k x ) is a solution for a given k y , then (λω, λk x ) is a solution for λk y . Thus the problem is completely solved once all the solutions for a given k y are found. Let be k y = 1. Also the upstream fluid speed u 1 can be set to 1 by redefining the ratio between the units of measure of frequency and wave number. We calculated all the solutions of the above equations, for these values of u 1 and k y , and display them in Tab. 1.
The first four solutions must be discarded because they correspond to waves propagating from downstream infinity to the shock. The fifth solution must be discarded too because it diverges exponentially at downstream infinity. More intriguing is the last solution. It seems to satisfy all the requirements in order to be a real physical solution and it actually satisfies conditions 79. Is it a real instability? No, because it represents a pressure perturbation advected by the fluid since ω − uk x = 0. From eq. 90.5 in Landau and Lifshitz (1987) we see that such a perturbation should have δP = 0. This is clearly absurd: we assumed k x to be the wave vector of a pressure perturbation in order to obtain eq. 80 but we have now found a solution with vanishing pressure perturbation. Note that the fifth solution is affected by this problem too.
Summarizing, we solved the equation for the shock eigenfrequency in the test particle regime. We found six solutions, but none has physical meaning, since they correspond to sound waves either propagating from downstream infinity to the shock surface or with vanishing pressure perturbation.
Below we apply our theory to two particular shock structures (one with multiple solution) in order to search for an instability. We do not want to carry out a systematic analysis of a set of solutions but we shall just show how our machinery works. For the sake of completeness in Fig. 3 we report the particle distributions f we used in our calculations. We adopted a system of units of measure so that the following three quantities equal 1: the speed of light, the proton mass and the numerical constant of the Bohm diffusion coefficient, i.e., D p = v(p)p. Below everything will be expressed in these units.
We proceeded in the following way. We calculated the determinant of the system of five eqs. 55, 56, 57, 58 and 80 and we set it to 0, obtaining an equation with 5 unknowns: ω, k xpu , k xtu , k xpd , k xtd , which are respectively the frequency and the x-components of the wave vectors of the upstream pressure wave, the upstream third mode, the downstream pressure wave and the downstream third mode. This forms a system of five equations together with four dispersions relations as in eq. 24, each one linking ω and k y with their respective k x . This is exactly what we did above when we solved eqs. 82 and 83. In Fig. 4 and in Figs. 5, 6, 7, we illustrate the solution of our system of equations as a function of k y , respectively for the single and the multiple solution of the shock structure. Absolute values of real parts (left panels) and imaginary parts (right panels) of ω (upper panels), k xpu , k xtu , k xpd , k xtd (lower panels) are plotted. Signs of these quantities are reported in Tab. 2. These solutions must  be all discarded because they have some waves diverging at infinity, just like the first four solutions in Tab. 1. We used as starting values in our searching algorithm each one of the six corresponding limit solution in the linear case. All the solutions we have found cannot be accepted for the same reasons discussed above for the linear regime.

Discussion
In essence, our method is exact, except for the short-wavelength (WKBJ) approximation necessary to treat analytically perturbations in the inhomogeneous upstream precursor. One may however wonder where our method differs from previous work (Mond and O'C. Drury 1998;Toptygin 1999), which has reported the existence of corrugational instabilities. Mond and O'C. Drury (1998) have reported the existence of both genuine corrugational instabilities, and of spontaneous emission of sound waves, for some (not all) of their models. There is of course a number of differences between this paper and theirs: we do not use the two-fluid approximations, and are interested in small-wavelength perturbations, contrary to them; we also notice the existence of perturbations in the upstream fluid, which they do not discuss. Toptygin (1999) included a number of novelties in his treatment, but he too did not notice that particles would diffuse upstream, so that he neglects the third and the d-modes altogether. Despite this, he does have perturbations upstream (pressure waves!), because he remarks that, for sufficiently long-wavelength perturbations, these become supersonic with respect to the fluid alone. This occurs because particles and fluid are tightly coupled in long wavelength perturbations by diffusion, which traps non-thermal particles, so that the restoring pressure is the sum of particles' and fluid's pressures, which is larger than the pure fluid speed. This is of course correct, and exists in our computations as well.
Contrary to these authors, we have found that shocks with non-vanishing particle pressure are stable to corrugational instabilities, even in the region of parameter space where multiple solutions are possible. We believe that the reason for this discrepancy lies in our inclusion of diffusion, and the abandonement of the two-fluid approach, as we now argue. We stated in the Introduction that the only possibility for the shock destabilization (polytropic shocks without particles are well-known to be corrugationally stable) lies in setting up a loop whereby perturbations shed by the shock in the downstream region return to the upstream region via particle diffusion and excite more perturbations; thus perturbation energy is not lost to downstream infinity, but returns to the shock to create more havoc. However, the very same mechanism that brings particles back upstream, also causes perturbations' damp- ening: we have seen that diffusion leads to damping both pressure waves and the so-called third modes. Also, we have established that a part of the perturbation occurs in the d-mode (where d stands for damping), where particles are perturbed in phase-space, but the total perturbation to the pressure vanishes. Thus, some fraction of the perturbation goes into a totally useless form (the d-mode), the new mode (the third mode) is always damped, and even pressure waves acquire a damping which is altogether neglected in the two-fluid approximation. In the end, while diffusion brings particles (and perturbations) back to the shock, the simultaneous damping is so strong to make the excitation of an instability ineffective.
We should remark that damping of pressure waves is weakest for the largest wavelengths, a limit which is unaccessible to us because of our WKBJ approximation, but is exactly the limit investigated by Mond and O'C. Drury (1998). While we have not found a large wavelength beyond which the shock becomes unstable, we cannot exclude that a proper treatment of the perturbations in the space-dependent precursor may yield a transition to the unstable regime.
A further comment is in order: of all possible geometries, the planar one is probably the least likely to display instability. Consider in fact a spherically symmetric explosion like a SuperNova. In this case, the perturbations (except of course for entropy and vorticity) generated downstream do not escape to infinity, as they do in the planar case, but return to the shock because the downstream region is finite and because they are deflected by a spatially-dependent refraction index; once they reach the shock, they may generate further perturbations. The situation is even more promising when the shock is due to an accretion flow, or is stalling: in fact, except for the presence of particles, this is exactly the scenario proposed (see Laming 2007 for an analytic approach and discussion) for the generation of asymmetries in proto-neutron stars: in this case, the mechanism for the instability of a stalled accretion shock is the reflection by the hard star surface into outgoing pressure waves of advected entropy perturbations, which return to the shock to generate more mischief. In the problem with particles, there is the extra complication due to diffusion, to be overcome to generate instability. Given the relative complexity of even a pure fluid analysis (Laming 2007), it is likely that this problem will require a numerical approach.
±∞, depending on whether we are considering downstream or upstream regions, respectively. The solution can be obtained by separation of variables, obtaining: δf = α(p, k x )e νt e kxx e ±ıkyy , subject to the constraint ν = D(k 2 x − k 2 y ) .
The sign of k x is the one that allows the solution to remain finite at infinity. Solutions belonging to different values of k x and p can obviously be superposed, but we know that there is another boundary condition (eqs. 48, 49 and following discussion) to be satisfied: at the shock, x = 0, δf = g(p)e ıωt e −ikyy (A4) which obviously gives ν = ıω , Dk 2 x = ıω + Dk 2 y , α(p) = g(p) .
The function g(p) was chosen so that p 3 vg(p) dp = D(p)p 3 vg(p) dp = 0 .
We have already seen (see the discussion leading to eq. 17) that these are the conditions for the vanishing of δP c and its first time derivative at the initial time, and thus at all times. The same property is of course acquired by α(p), so that the mode we just found is surely a d-mode for the upstream region, in the Homogeneous Approximation.
When we assume u = 0, the above formulae remain correct except for the substitution ν → ν + uk x .

B. Explicit derivation of functions w
We need to find the solution to − DA iu δf iu (ık xiu + k xd ) + DA id δf id (ık xid + k xd ) .
This is an equation for g u (p):