On the stability of shocks with particle pressure

# On the stability of shocks with particle pressure

Stefano Finazzi111Currently at SISSA, Trieste. and Mario Vietri Scuola Normale Superiore, Pisa
###### Abstract

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, showing that they are stable. We then discuss the differences between our treatment and previous work.

shock waves – cosmic rays

## 1 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, , 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 () 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 levels222The situation is much less clear for UHECRs, because many of the proposed classes of sources are purely hypothetical, so that there are no constraints on their properties..

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 sub-shock), a convergent spectrum will be obtained, with a definite value for the parameter . 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 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 amplitude333In all of this paper, the term corrugation instability is taken to include also the spontaneous emission of sound waves by the shock.. 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, . 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 () 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 , and it will be allowed to be an arbitrary function of both , 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 on 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, , and to flow out of each of them, , be ordered like this: . 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 , with fluid coming from the left, and exiting from the right, so that all speeds are . 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 , which restricts our analysis to the regime , where is the perturbation wavelength perpendicular to the shock, and 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.

## 2 Perturbations in a homogeneous medium

We give here, for future reference, our basic equations. They are the conventional hydrodynamic equations:

 ∂ρ∂t+▽⋅(ρ→v)=0 (1)
 ∂∂t→v+→v⋅▽→v=−1ρ▽(P+Pc) (2)
 ∂s∂t+→v⋅▽s=0 (3)

which contain a term for the momentum exchange between the fluid and the non-thermal particles represented by the gradient of the particle pressure , plus the usual Boltzmann equation in Skilling (1975):

 ∂f∂t=▽⋅(D▽f)−→v⋅▽f+13▽⋅→vp∂f∂p. (4)

We assume to be a given function of and .

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

 δs∝exp(ıωt−ı→k⋅→r) (5)

so to obtain, from the equation of entropy conservation, eq. 3,

 (ω−ukx)δs=0→ω−ukx=0, (6)

where is the unperturbed fluid velocity in the x-direction. Perturbation of the mass conservation equation yields

 (ω−ukx)δρ−ρ→k⋅δ→v=0→→k⋅δ→v=0. (7)

Perturbation of the momentum equation yields

 (ukx−ω)δ→v+→kρ(δP+δPc)=0→δP=−δPc. (8)

Now, the perturbation of the Boltzmann equation yields

 ıωδf=−k2Dδf+ıukxδf−ı3→k⋅δ→vp∂f∂p (9)

which implies because, for entropy perturbations, and , as deduced above. Furthermore, since , necessarily , and thus . 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, :

 δs=δs∘eıωt−ıkyy−ıkxx ω−ukx=0 δP=δPc=δf=δ→v=0 δVV=γ−1γmδskB (10)

where the last equation applies to ideal fluids: is the average particle mass and Boltzmann’s constant.

We turn now to isentropic perturbations, which automatically satisfy eq. 3. Mass conservation implies

 (ω−ukx)δρ=ρ→k⋅δ→v (11)

while momentum conservation implies

 (ω−ukx)δ→v=→kρ(δP+δPc). (12)

Let us begin by assuming that the left-hand-side of eq.11 vanishes; the same is then true for . Multiplying eq. 12 by , we see that , unless . If the curl vanishes, then so does , 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 for perturbations with non-zero vorticity; we thus find that , and then, again from eq. 9, that , implying : vorticity perturbations do not couple to particles either. We have for them:

 δs=δρ=δP=δf=δPc=0 ω−ukx=0 →k⋅δ→v=0. (13)

We consider now those perturbations where does not vanish. We can show that here too there are two distinct classes of modes: in the first one, is not coupled to the fluid quantities, while in the second one it is through the term . 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

 ∂δf∂t=∇⋅(D∇δf)−u∂δf∂x, (14)

where we dropped the term 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 is the initial distribution function at time (possibly dependent upon ), the solution at later times for is

 δf(x,y,t,p)=∫dx0∫dy0ϕ∘(x0,y0,p)14πDte−(x−x0)24Dte−(y−y0)24Dt, (15)

and the solution for is just . At the same time, we must make sure that this solution does not ruffle the fluid: this obviously requires at all times. Now integrating eq. 14 over we find

 ∂δPc∂t=4π3∂2∂x2∫D(p)δfp3vdp+u∂δPc∂x, (16)

which clearly shows that, in order to have everywhere at all times we need everywhere at the initial time. If we now Fourier-expand the initial condition with respect to , the above conditions become

 ϕ∘(x,y,p)=∫dkx∫dkyg(p,→k)e−ıkxxe−ıkyy,∫g(p,→k)p3vdp=∫g(p,→k)D(p)p3vdp=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 , which does not have a Fourier transform with respect to .

When , velocity perturbations are induced in the fluid by the non-vanishing particles’ pressure gradient, and eq. 9 can be solved for as

 δf=−δρρp3∂f∂pı(ω−ukx)ı(ω−ukx)+k2D. (18)

If we integrate this over , we find

 δPc=δρρ∫pMpm4π9dpvp4∂f∂p−ı(ω−ukx)ı(ω−ukx)+k2D, (19)

where we assumed the non-thermal particles to have a minimum () and a maximum () momentum. It is convenient to introduce a weighted diffusion coefficient defined as follows:

 11−z¯D(z)≡∫pMpm4π9dpvp4∂f∂p11−zD(p)∫pMpm4π9dpvp4∂f∂p, (20)

which shows to be a function of only.

The above can be simplified a bit by integration by parts. If the integral were to extend from to , we would have

 −∫∞04π9dpvp4∂f∂p=∫∞04π9dp(ddpvp4)f=4π3∫∞0dpfp3v(43+m23E2)≡γcPc (21)

where is the particle pressure in the unperturbed flow, and , the effective particle polytropic index, satisfies . Since however the integral extends only from to , we define

 −∫pMpm4π9dpvp4∂f∂p=γ′cPc. (22)

Thus eq. 19 can now be rewritten as:

 δPc=γ′cPcδρρ11−ık2ω−ukx¯D(ık2ω−ukx). (23)

We remark that both and 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 between eqs. 11 and 12, and then use (with the sound speed because we are considering isentropic perturbations) and the equation above to obtain

 Ω2=1+γ′cPcγP11−ık2ω−ukx¯D(ık2ω−ukx) (24)

where we have called

 Ω≡ω−ukxkcs (25)

the comoving eigenfrequency, in suitably scaled units.

Eq. 24, together with the definition of (eq. 20), is the sought-after dispersion relation for the coupled small perturbations in a homogeneous medium.

### 2.1 Properties of coupled modes

The case where is independent of 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 , and thus , is a constant, independent of . The eq. 24 then reduces to

 (Ω−ıkDcs)(Ω2−1)=Ωγ′cPcγP, (26)

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 . 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 . 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 , because one must solve simultaneously eqs. 24 and 20. We begin by considering the limit . In this case, and assuming a constant, we easily find, to leading order in , : 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 higher order terms in . We find:

 Ω=ık¯D(ık2/(ω−ukx))cs11+γ′cPcγP. (27)

Now we use the definition to simplify the above to

 1−z¯D(z)=−γ′cPcγP (28)

which can be now used with eqs. 20 and 22 to obtain

 γP=∫pMpm4π9dpvp4∂f∂p11−zD(p), (29)

where of course . As a function of real , the right-hand side above (where it exists) is easily seen to be a monotonically decreasing function of , vanishing for . In any realistic physical problem the integral must extend from a minimum () to a maximum momentum ; since is expected to be a monotonically increasing function of , we see that the integral above always exists for , and , and it diverges exactly at and . Thus the integral on the right-hand side of the equation above spans the whole range from to as varies between and , and the range to as varies between and . Somewhere in the range 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 , and the small perturbations were assumed to vary as , the result that implies that all modes are damped by diffusion, as physical intuition obviously suggests.

The discussion in the opposite limit, , is similar. If we assume a constant, we find the solution

 Ω2=1, (30)

without the correction to the sound speed due to the presence of the particles’ pressure: in the limit 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 lower order terms in . We obtain now:

 Ω=ık¯D(ık2/(ω−ukx))cs (31)

which can be rewritten as

 z¯D(z)=1. (32)

Comparing this with eq. 20, we see that the value of 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

 z=1D(pm). (33)

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.

For future reference, we give the expressions for all small quantities, in units of : the link between and is obtained from Eqs. 20 and 22, and the others follow easily. For later convenience, we use again rather than :

 δPc=δPc∘e(ıωt−ıkyy−ıkxx) δVV=−1γ′cδPcPc(1−ık¯DcsΩ)=−δPcγP(Ω2−1)≡−qδPcPc δ→v=−→kkcsΩδVV⎛⎜⎝1+γ′cPcVc2s11−ık¯DcsΩ⎞⎟⎠=−→kkΩcsδVV≡→zδVV δPP=−γδVV, (34)

where the definitions of and will come handy later on. In these equations, must be regarded as a known function of and , specified by eq. 24.

### 2.2 Initial value problem

As a preparation for later work, we discuss the initial value problem. This has some interest because perturbations in belonging to the various modes are not mutually orthogonal, and thus it may appear that initial conditions, especially when given only in terms of , 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 at time 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 , calling its amplitude. We must then have

 a(p)=Aiδfi+g(p), (35)

where we have introduced some notation that will be useful in the following: a summation convention over is understood, the ’s are the amplitudes of the pressure waves and third modes for , respectively; is the amplitude of the d-mode, as in eq. 17, and the ’s are

 δfi=−13p∂f∂pı(ωi−ukxi)ı(ωi−ukxi)+k2iD(p). (36)

Comparing this with eq. 18, it becomes clear that the mode amplitudes ’s are simply . Here the quantities are supposed to be linked by the appropriate branch of the dispersion relation. For 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:

 ∫p3va(p)dp=Ai∫p3vδfidp (37)
 ∫D(p)p3va(p)dp=Ai∫D(p)p3vδfidp. (38)

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. 2.1, we see that the two conditions, the vanishing of the density and of the velocity, imply

 ∑iAi=0,ΩiAi=0, (39)

which are two more linear relations which, together with eqs. 37 and 38, determine the ’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.

## 3 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:

 ∂δs∂t+u(x)∂δs∂x=0 (40)

where we could neglect the term 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, , is:

 (∂∂t+→v⋅∇)→ηρ=(→ηρ⋅∇)→v+1ρ3∇ρ∧∇(P+Pc). (41)

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, . But in the presence of particles and of spatial gradients it is easily seen that the particle pressure 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 , and in particular , the size of the upstream precursor. In other words, we take all perturbed physical quantities to be of the form

 δX≈QX(x)eıωt+ı∫0xWX(x′)dx′−ıkyy, (42)

and all derivatives in the direction (i.e., perpendicular to the shock) are considered small when compared with terms proportional to . This is a perturbation analysis in which the transverse wavenumber is assumed large and as a consequence the longitudinal wavenumber is also large. The presence of a non-constant amplitude is equivalent to keeping the first two terms in an asymptotic expansion in the small parameter . 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 further from the shock is immaterial. We see from the above that all physical quantities close to the shock satisfy

 δX≈QX(0)eıωt−ıWX(0)x−ıkyx. (43)

This is precisely the same form that holds in the homogeneous medium, so that eqs. 2.1 and 17 still hold. Also, the space-time dependence of the d-mode for the half-plane 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 , the size of the upstream precursor.

## 4 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 also upstream, so that we need to consider appropriate conditions to match in the two regions. The first condition is obviously the continuity across the shock,

 δf1=δf2. (44)

It is also well-known that the spatial gradient of 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 is where is the shock corrugation amplitude. So, we obtain:

 ^n⋅(D∇f)2−^n⋅(D∇f)1+13p∂f∂p^n⋅(→v2−→v1)=0. (45)

Writing its first-order linearization and using eq. 44 we find:

 D(∂δf∂x|2−∂δf∂x|1) + ıkyζD(∂f∂y|2−∂f∂y|1) (46) + 13p∂f∂p(δv2x−δv1x)+13p∂(δf)∂p(u2−u1)=0,

We note that because at zero-th order the medium is uniform in coordinates parallel to the shock surface. The result is:

 (47)

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 on the downstream side of the shock (the factor will be omitted for simplicity in this subsection) satisfies

 δf+=Aidδfid+gd. (48)

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 and , exist for opposite values of (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

 δf−=Aiuδfiu+gu, (49)

The continuity of at the shock allows us to derive a relationship between the ’s:

 gd=gu+Aiuδfiu−Aidδfid. (50)

Using

 ∂δf+∂x = −ıAidkxidδfid+kxdgd=−ıAidkxidδfid+kxdgu+Aiukxdδfiu−Aidkxdδfid (51) ∂δf−∂x = −ıAiukxiuδfiu+kxugu (52)

into eq. 47 we find

 D(kxd−kxu)gu+u2−u13∂gu∂lnp= − u2−u13Aiu∂δfiu∂lnp−13∂f∂lnp(δv2x−δv1x) (53) − DAiuδfiu(ıkxiu+kxd)+DAidδfid(ıkxid+kxd)

which we regard as an equation for , whose solution can be written as

 gu(p)=CwC(p)+Aiuwiu(p)+Aidwid(p)+(δv2x−δv1x)w0(p), (54)

where the functions ’s are derived in the Appendix B.

We can now impose the conditions 17 on , obtaining:

 C∫p3vwCdp+Aiu∫p3vwiudp+Aid∫p3vwiddp+(δv2x−δv1x)∫p3vw0(p)dp=0, (55)
 C∫p3vwCDdp+Aiu∫p3vwiuDdp + Aid∫p3vwidDdp (56) + (δv2x−δv1x)∫p3vDw0(p)dp=0,

And now we can impose the very same conditions on , eq. 50, to obtain:

 Aid∫p3vδfiddp−Aiu∫p3vδfiudp=0 (57)
 Aid∫p3vδfidDdp−Aiu∫p3vδfiuDdp=0 (58)

This set of four linear equations in five unknowns ( and the ’s) can be solved in terms of one of them, say , 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. 2.1, 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.

## 5 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 rather than the density directly.

We consider a small-amplitude corrugation of the shock surface, away from the plane, by a small displacement of the form:

 ζ=ζ∘eıωt−ıkyy (59)

with respect to which the unit vectors parallel and normal to the surface have components in the plane:

 ^t=(−ıkyζ,1);^n=(1,ıkyζ) (60)

while the surface speed in the direction normal to the surface, with the respect to the reference frame of the unperturbed shock, is:

 →q⋅^n=ıωζ. (61)

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:

 (→v1+δ→v1)⋅^t=(→v2+δ→v2)⋅^t (→v1+δ→v1)⋅^n−(→v2+δ→v2)⋅^n=√(P2+δP2−P1−δP1)(V1+δV1−V2−δV2) (62)

whose first-order linearizations are:

 δv2y−δv1y=ıkyζ(v2−v1) δv2x−δv1x=12(v2−v1)(δP2−δP1P2−P1−δV2−δV1V1−V2). (63)

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):

 V2V1=(γ+1)P1+(γ−1)P2(γ−1)P1+(γ+1)P2 (64)

whose perturbation yields

 δV2V2=δV1V1+h(δP1P1−δP2P2) (65)

where we have defined

 h≡4γ((γ+1)+(γ−1)P2/P1)((γ−1)P1/P2+(γ+1)) (66)

while the ratio can be expressed as

 P2P1=2γM21−(γ−1)γ+1 (67)

with 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:

 j2=P2−P1V1−V2 (68)

where is the mass flux. When the shock surface is perturbed, the above becomes

 ((→v1+δ→v1)⋅^n−→q⋅^n)2(V1+δV1)2=P2+δP2−P1−δP1V1+δV1−V2−δV2. (69)

The perturbation to first order of the above yields

 2δv1xv1−2ıωζv1−2δV1V1=δP2−δP1P2−P1−δV1−δV2V1−V2 (70)

Lastly, the condition that be continuous across the shock has already been imposed, see eq. 44.

## 6 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 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 propagating 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 .

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:

 δV2 = δV2s+δV2p+δV2t, (71) δ→v2 = δ→v2v+δ→v2p+δ→v2t, (72) δP2 = δP2p+δP2t, (73)

where the subscripts , , and label quantities for entropy, vorticity, pressure waves and third mode. Upstream there are neither entropy nor vorticity perturbations:

 δV1 = δV1p+δV1t, (74) δ→v1 = δ→v1p+δ→v1t, (75) δP1 = δP1p+δP1t. (76)

We must use these expressions in the perturbed Rankine–Hugoniot. Furthermore, we can write one of the two components of in terms of the other one, through , where and (see eqs. 2). Then we can simplify our system by eliminating the shock displacement , the remaining component of the velocity perturbation due to vorticity, and . One equation remains, linking only pressure waves and third modes’ perturbations:

 ωP2(v1−v2)V1{V1P1(ω2−ky2v1v2) −[ω2((1+h)P1−hP2)+ky2((h−1)P1−hP2)v1v2]V2}(δP1p+δP1t) −ωP1(v1−v2)V1{V1P2(ω2−ky2v1v2) −[ω2(hP1−(h−1)P2)+ky2(hP1−(1+h)P2)v1v2]V2}(δP2p+δP2t) −2ω2kyP1(P1−P2)P2v2V1(V1−V2)(δv1py+δv1ty) +2ω3P1(P1−P2)P2V1(V1−V2)(δv2px+δv2tx) +2ω2kyP1(P1−P2)P2v2V1(V1−V2)(δv2py+δv2ty)=0. (77)

Now, eqs. 2.1 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 and its own but, for given and , they are fixed by the dispersion relation, eq. 24. We obtain:

 −hγ(P1−P2)(v1−v2)(