Linear analysis of the gravitational beam–plasma instability

We investigate the well-known phenomenon of the beam–plasma instability in the gravitational sector when a fast population of particles interacts with the massive scalar mode of a Horndeski theory of gravity, resulting in linear growth of the latter amplitude. Following the approach used in the standard electromagnetic case, we start from the dielectric representation of the gravitational plasma, as introduced in a previous analysis of the Landau damping for the scalar Horndeski mode. We then set up the modified Vlasov–Einstein equation, using a Dirac delta function to describe the fast beam distribution. We thus provide an analytical expression for the dispersion relation, and we demonstrate the existence of a nonzero growth rate for the linear evolution of the Horndeski scalar mode. A numerical investigation is then performed with a trapezoidal beam distribution function, which confirms the analytical results and allows us to demonstrate how the growth rate decreases as the beam spread increases.


I. INTRODUCTION
One of the most intriguing features of plasma phenomenology is Landau damping [1,2], i.e. the decay of an electromagnetic wave amplitude even when it propagates through an ideal plasma [3].Such a peculiar property possessed by an electromagnetic plasma differently from any other medium, is due to the non-microscopic scale of the Debye length [4], which describes the plasma quasi-neutrality.From a phenomenological point of view, electromagnetic waves in a plasma acquire a longitudinal polarization (say the photon takes an effective non-zero mass [5][6][7][8]) and, when the phase velocity of the Langmuir modes is of the same order of the thermal velocity of the plasma constituents (ions and electrons), energy can be transferred from the former to the latter.Since the Landau damping has an ideal (reversible) character, it is also possible to demonstrate [9][10][11][12][13] that a fast population of particles can transfer its kinetic energy to a Langmuir mode, enhancing the wave amplitude up to a nonlinear saturation value.Indeed, the beam-plasma instability is conceptually a reversed Landau damping phenomenon, at least in its linear regime.In fact, since the interaction of an electromagnetic wave with an ideal plasma is a statistically reversible process, depending on the specific conditions of the system (mainly the resonance condition that the wave frequency is near the plasma one), the amplitude of the wave can be either suppressed (Landau damping) or enhanced (inverse Landau damping).For the beam-plasma physics of interaction, it is relevant the possibility to transfer energy from a fast beam of charged particles to the so-called Langmuir waves of the plasma, i.e. self-consistent electrostatic modes living in thermal plasma.As shown in [11], when the resonance condition is fulfilled the topology of the dispersion relation acquires a peculiar structure.Actually, in order to preserve the dielectric representation of the plasma linear response, the beam must have a tenuous number density with respect to the plasma one and the ratio of the two densities is a fundamental small parameter of the linear and non-linear dynamics.As a consequence of the smallness of this ratio, the dispersion relations in the (k, ω)-space (we are referring to a one-dimensional problem) are associated to the vanishing of either the plasma dielectric function, namely when we are dealing with a Langmuir mode, or when the wave phases are comoving with the beam, i.e. ω = kv B , being v B the beam speed.The overlap of these two conditions, naturally leads to the emergence of an unstable mode, corresponding to the amplification phenomenon of the self-consistent Langmuir wave present in the plasma.It is worth noting that, since the beam is fast, the Langmuir modes interested to this enhancing process have a phase velocity v ph ∼ v b , much greater than the thermal speed of the background plasma and, therefore, they do not suffer for the Landau damping, very efficient only near the point of inflection of the Maxwellian distribution function of the electrons.Thus, we see that the interaction of a fast beam with a thermal plasma is due, for what concerns its linear instability behavior, to the inverse Landau damping that the fast particle induces on the existing resonant Langmuir modes.The process is also characterized by a non-linear dynamics (see [12]), in which the wave amplitude saturates, but this part of the interaction is not discussed in the gravitational parallelism here proposed.Despite the fundamental differences between the electromagnetic and gravitational interaction it is possible to draw a parallelism between them, on the ground of the common results obtained for what concerns the phenomenology associated to direct and inverse Landau damping.Indeed, in the case of a gravitational plasma (see [14][15][16] for pio-neering treatments, [17] for a comprehensive review and section II for an in-depth discussion) the inertial forces can be considered to act as a neutralizing background; in other words, in the local inertial frame a set of massive particles feels the gravitational perturbation of an incoming wave only, i.e. the static background curvature generated by the whole medium can be gauged out with a convenient choice of the coordinates.In [18] it has been shown that, for the case of the dynamics predicted by a Horndeski theory [19][20][21][22][23], the associated scalar mode is subject to Landau damping, in an analogous way as in the standard electromagnetic sector.
This result opens an interesting scenario for studying the amplitude profile of the Horndeski scalar mode when it interacts with a gravitational medium (for a discussion of the standard gravity case, see [24][25][26][27][28][29]). Here, starting from the dielectric gravitational function derived in [18], we analyze the complementary question to the gravitational Landau damping, i.e. the so-called gravitational beam-plasma instability.Namely, we investigate the interaction of gravitational Langmuir modes, having a dispersion relation which annihilates the dielectric function, with a fast population of massive particles propagating through the medium.
We show that the beam-plasma instability takes place in the gravitational sector as in the electromagnetic case.A non-zero growth rate of the Langmuir waves is present as long as the beam velocity and the phase velocity of the gravitational scalar Horndeski mode are comparable.The analytical study is developed by describing the beam distribution via a Dirac delta function, say by dealing with a null temperature beam.As in the electromagnetic sector, the instability arises in the presence of a degenerate point of the dispersion relation of the plasma and beam components.The significance of this condition is to ensure that we are dealing with a Langmuir mode living in the gravitational medium, with a phase velocity matching the beam one, allowing for the resonance phenomenon at the ground of the energy transfer.
It must be remarked that our study demonstrates that the instability takes place also in the limit of vanishing scalar mode mass, so that a phenomenologically viable model for the beam-plasma instability is always attainable without violating the current constraints on the graviton mass [30][31][32][33][34][35].Moreover, as the beam velocity approaches the speed of light, the gravitational instability is suppressed, simply because the Jüttner medium distribution function has a natural population cutoff at that scale and no real Langmuir modes can emerge in that parameter region.
A numerical analysis is performed using a trapezoidal form of the beam distribution function, able to model a finite (non-zero) temperature, while still allowing for a semi-analytical treatment.This numerical study confirms the main results obtained in the cold beam scenario, with the additional important feature that the amplitude of the growth rate decreases as the width of the trapezoid increases, i.e. at larger beam temperatures.We conclude by observing that both the analytical and numerical treatments describe a profile of the growth rate peaked around the critical wavenumber of the degenerate point discussed above.
The paper is structured as follows: in section II we elucidate the physical setting we intend to work with, namely we introduce the concept of "gravitational plasma".In section III we briefly present the linearized field equations for tensor and scalar gravitational waves from Horndeski theories of gravity on a Minkowski background.In section IV we review the analysis of the interaction between gravitational radiation and a medium of massive particles in which collisions are neglected, outlining the possibility of Landau damping for the scalar mode.In section V we extend this previous result by studying the behavior of self-consistent scalar waves, here denoted as gravitational Langmuir modes, in the case in which a tenuous beam of particles, distributed as a Dirac delta function, is injected in the medium, evaluating explicit expressions for the degenerate wavenumber and the correspondent amplitude growth rate.In section VI we analyze, by making use of numerical methods, a more realistic setting for the beam-plasma scenario, namely we adopt for the beam distribution a trapezoidal shape, in order to investigate the dependence of the instability from the temperature of the beam.In addition to this, we confirm the analytic results obtained for the delta beam, particularly the behavior of the growth rate with respect to the mass of the scalar mode and to the velocity of the beam.Finally, in section VIII we comment the results obtained and draw the final conclusions.

II. THE CONCEPT OF GRAVITATIONAL PLASMA
In this section we analyze the concept of a "gravitational plasma" and we summarize the basic analogies and discrepancies with respect to the electromagnetic case.What makes really different the gravitational interaction from the electromagnetic one is, apart from the non-linearity of its dynamics, the absence of charge of opposite sign, able to shield the gravitational field generated by the different type of sources.This fact could suggest that, in a neutral medium sensitive to the gravitational interaction, there is no chance to define an analogous quantity to the Debye length of an electromagnetic plasma.By other words, it seems impossible to recover the concept of a "neutralizing background" for a gravitational medium, like the role played by ions in a cold plasma.However, the Equivalence Principle offers an intriguing point of view to solve this puzzling question when describing a medium as a gravitational plasma, first explored in the pioneering work [14] and explicitly stated in [18].The point is that, in a local inertial frame, the background gravitational field can be actually canceled and therefore it emerges the idea that, at least in a small region of the space-time, a gravitational analogy of the neutralizing background is provided by the inertial forces.Hence, a gravitational medium in a local inertial frame is essentially in a "quasi-neutrality" condition and its selfconsistent gravitational fluctuations are not far from the physical character of the so-called Langmuir waves [2,36].Clearly, the analogous of the Debye length is here the spatial scale on which the tidal forces are not significantly appreciable in the considered system, i.e. the scale up to which it is possible to speak of an inertial frame of reference.Once fixed this theoretical framework, it is possible to study Vlasov equation coupled to the linear wave equation governing the dynamics of the gravitational radiation in the theory of gravity under consideration.When the analysis of gravitational waves from general relativity interacting with a neutral gas is addressed, however, the gravitational parallelism with ordinary plasma physics seems to fail, simply because no damping of the gravitational waves can take place in an ideal gravitational medium [24][25][26][27][28][29].Indeed, as clearly elucidated by a number of works in the literature, tensorial massless gravitational waves in a collisionless medium are characterized by a superluminal phase velocity at all frequencies, so that it is impossible for any massive particle to be resonant with a specific mode and exchange energy with it.It must be remarked that all the studies dedicated to the general relativity case analyze only the behavior of tensor polarizations, neglecting the possibility of the arising of scalar and vector effective modes within the medium, to which could be associated a different phenomenology (for the emergence of extra effective polarizations for gravitational waves in a molecular medium, see [37]).Thus, in [18] (see also [38] on the case of parity-violating theories of gravity), the gravitational version of the Landau damping has been searched and actually found in the context of Horndeski theories of gravity.Specifically, it has been found that the massive scalar mode, corresponding to an intrinsic longitudinal fluctuation mode [23], can be actually characterized by subluminal phase velocity within the gravitational medium and be resonant with specific modes.On the contrary, the tensorial part of the gravitational radiation can not be subject to such phenomenon, just as it occurs in the general relativity case.In this work our intention is to study the inverse Landau damping phenomenon through the beamplasma interaction, i.e. the possibility for a fast particle population to make the gravitational medium Langmuir modes unstable, so transferring their energy to the scalar mode of a Horndeski gravity model.Our findings demonstrate that the inverse Landau damping takes place as long as the particles composing the beam have velocities compatible with the allowed range of phase velocities for the Langmuir gravitational modes.The present result is of interest because the amplitude of the massive scalar mode in a propagating gravitational fluctuation is expected small and rather difficult to be detected with the present ground-based interferometers, as well as in the near future with space-based instruments.Thus the possibility to enhance the intensity of such a contribution via the interaction with a fast particle beam, could improve the detection chance of such non-Einsteinian modes, via their phenomenological impact on bounded gravitational systems across the Universe.Finally, we observe that the recent detection of a multi-messenger signal [39] has demonstrated that the velocity of propagation of a hypothetical massive mode has to be extremely close or coinciding with the speed of light.It is remarkable that the beam-plasma instability remains present even when the mass of the scalar mode tends to zero, according to the requirement that its propagation speed is very near the value of c.In this respect, we can claim that the new gravitational instability is still viable given the presentday detections of the gravity phenomenology.

III. GRAVITATIONAL WAVES FROM HORNDESKI THEORIES
As well known [19][20][21], any scalar-tensor gravitational theory with second order equation of motion can be derived from Horndeski action where the terms L i have the following expressions in which R is the Ricci scalar and G µν is the Einstein tensor.We have also introduced the compact notation A choice of the free functions K and G i corresponds to the selection of a definite second-order theory of gravity (for instance, general relativity is recovered for G 4 = 1 and keeping null all others functions).It must be remarked that the detection of the multi-messenger signal GW170817-GRB170817A [39] has set severe constraints on the parameter space of Horndeski action: specifically, in order to reproduce the observed physics the function G 4 must depend only on the scalar field ϕ and the G 5 must be constant.Moreover, a recent work [40] demonstrates that also the function G 3 should be independent from the kinetic term X, in order to satisfy Witten positive energy theorem [41].However, we point out that these limitations on the choice of the free functions characterizing Horndeski action do not affect the results obtained in the present work.We are interested in describing gravitational waves propagation on a Minkowski background, therefore we perform the usual splitting of the metric tensor where η µν is Minkowski spacetime metric, namely η µν = diag (−1, 1, 1, 1), and h µν is the metric perturbation satisfying |h µν | 1 in all spacetime points in which ( 4) is valid.We also operate an analogous decomposition on the scalar field by setting with φ 0 representing the background value and φ a perturbation, whose size is assumed to be |φ| |φ 0 |.By adding the matter contribution in the action (1), we write the linearized field equations for the metric and scalar field perturbations (for details of the derivation of the field equations in vacuum, see [22]) as The functions K and G i , together with their derivatives, are all evaluated at the point ϕ = φ 0 and the kinetic term is retained null, i.e.X = 0.The notation G µν indicates the linearized version, with respect to h µν , of the Einstein tensor, whereas T µν and T (1) are the perturbations in the stress-energy tensor and in its trace induced by the presence of gravitational waves.The scalar field perturbation φ obeys the Klein-Gordon equation (7) with mass given by and the positivity of this term must be ensured in order to prevent a tachyonic behavior in vacuum.The coupling constants are expressed as , where κ = 8πG is the Einstein gravitational constant 1 .The contributions from the metric and scalar field perturbations in equation ( 6) can be decoupled by defining the generalized trace-reversed tensor 1 In all sections we adopt natural units c = 1.
with h = η µν h µν and α = G4,ϕ(0) G4(0) .It results that the new metric perturbation is solution of an inhomogeneous d'Alembert equation which must be solved together with the additional constraints ∂ µh µν = 0, h0i = 0 and h = 0. Thus the tensor hµν contains only two degrees of freedom, carrying the usual plus and cross polarizations acting on the transverse plane with respect to the direction of propagation.

IV. GRAVITATIONAL LANDAU DAMPING
In a recent work [18], the propagation of gravitational waves from Horndeski theories in an isotropic collisionless medium of massive particles was analyzed.In this section, we resume its key points and present the main result, namely the possibility of energy exchange between the background medium and gravitational waves through Landau damping for the scalar degree φ only.The medium is described through a distribution function f ( x, p, t), defined on the single-particle phase space and normalized in order to give the total number of particles when integrated over the entire domain and the density of particles when integrated on all momenta We choose to retain the contravariant components of the position x i and the covariant components of the momentum p i as generalized coordinates in phase space, in order to deal with simpler equations in the following.Therefore, the phase space volume element is d 3 xd 3 p ≡ dx 1 dx 2 dx 3 dp 1 dp 2 dp 3 .The stress-energy tensor of the set of particles is obtained from with g the determinant of the metric tensor g µν given in (4).The time evolution of the medium distribution function is ruled by Vlasov equation Let us consider a medium that, in the absence of gravitational waves, has reached thermal equilibrium and whose distribution function is an isotropic, homogeneous and time-independent configuration f 0 (p), where p ≡ δ ij p i p j is the momentum Euclidean modulus.At the initial time t = 0, we turn on the disturbance provided by gravitational radiation and the ensemble of particles is driven out of equilibrium by a force term where p 0 = m2 + g ij p i p j is the particle energy and m its mass, calculated from the geodesic equation for the metric (4).For continuity, the distribution function must obey f ( x, p, 0) = f 0 ( g ij ( x, 0)p i p j ), which at first order in the perturbations results in where f 0 (p) ≡ df0 dp .For any positive time, a perturbation δf arises in the distribution function, as response of the medium to the presence of gravitational perturbing fields.We demand that the size of this perturbation with respect to the equilibrium configuration be of the same magnitude of that of gravitational waves.Then, Vlasov equation ( 16) can be linearized with respect to the distribution function perturbation, resulting in The set of equations ( 7), ( 12) and ( 19) results closed once that the source terms in the wave equations are evaluated in terms of the distribution function perturbation, namely We search for plane wave solutions of this differential problem, choosing the z axis of our coordinate system to be coincident with the waves direction of propagation.Then, by performing a Fourier transform on the spatial coordinate z, labeled by a real Fourier parameter k, together with a Laplace transform on the time coordinate t, with complex Laplace parameter s, the solution of the linearized Vlasov equation ( 19) can be readily written as Here we denote with h(k) ij (0) and φ (k) (0) the initial value of the projections in the Fourier space of the unknown functions hij (z, t) and φ(z, t).The dependence of δf (k,s) ( p) from both h(k,s) ij and φ (k,s) could result in a coupling between scalar and tensor modes in the wave equations, but it turns out that the spurious contributions cancel out due to the symmetries of integrals (20) and (21).Hence, the wave equations for h(k,s) ij and φ (k,s) remain decoupled at linear order and we can calculate their solutions as where we have defined the complex dielectric functions describing the allowed modes for scalar and tensor gravitational waves throughout the medium.Specifically, by introducing the complex frequency ω ≡ is, with ω r and ω i its real and imaginary part, it can be shown that scalar and tensor perturbations within the medium are damped or growing waves of the form e −i(ωt−kz) = e ωi(k)t e −i(ωr(k)t−kz) , where the dispersion relation ω r (k) is the curve along which the real part of the dielectric function results null whereas the characteristic frequency ω i (k) is calculated from In order to obtain these quantities for scalar and tensor perturbations we first have to specify the equilibrium distribution function f 0 (p).For this purpose, we assume a Jüttner function with n the particles density, Θ the medium temperature in units of the Boltzmann constant k B and K l (•) the modified Bessel function of the second kind.The parameter ζ ≡ m Θ , namely the ratio between rest and thermal energy, quantifies the magnitude of relativistic effects, reproducing the ultrarelativistic limit for ζ → 0 and the Newtonian one for ζ → ∞.The explicit expression of the dielectric functions for scalar and tensor perturbation in the case of a Jüttner background is In plasma physics this kind of integral is usually addressed by assuming that the phase velocity of the signal v p ≡ ωr k is much greater than the mean thermal velocity of particles v T and expanding the denominator as a power series in the small parameter v T vp .Then, integrating term by term and truncating at a chosen order, a real quantity, identified with the real part of the dielectric function r , is obtained.From this function it is possible to extract the dispersion relation ω r (k).The imaginary part i is instead calculated by exploiting the residue theorem for the integration around the Landau pole (u − v p ) −1 on a semicircular path in the complex upper half-plane.In our case we can of course make the request that the phase velocity be much greater than the thermal velocity of particles.Indeed, for ζ 1 the mean velocity predicted by Jüttner distribution rapidly reaches values of order 10 −1 or less.Therefore, for a relativistic phase velocity value close to unity, the error associated to a truncated power series in v T vp can be made satisfactorily small.However, here the situation is a bit more involved with respect to the electromagnetic case, since we have to deal with a more elaborate form of the denominator contained in the integrals, where the phase velocity plays a role in determining the position of the pole.Specifically, we see that a necessary condition for the presence of a pole in the domain of integration is a subluminal phase velocity for the perturbations within the medium.Indeed, in this case, the dielectric functions will acquire a imaginary contribution coming from the poles at the points ).On the contrary, for a superluminal phase velocity the Landau pole lies outside the domain of integration, thus we have a purely real dielectric function and, consequently, a vanishing imaginary part of the frequency.In all the cases, the condition v p < 1 must be checked a posteriori from the analysis of the dispersion relation obtained, through formula (27), from the real part of the dielectric function calculated as a truncated power series.This being said, we report the dispersion relation for tensor modes where we have introduced the function γ(ζ) ≡ K1(ζ) K2(ζ) and the proper frequency for tensor excitations ω 2 h ≡ κ mn 6 = κ ρ 6 , with ρ = mn the mass density characterizing the medium.From this expression, it is immediate to show that the phase velocity for tensor waves results always greater than the speed of light.Therefore, the dielectric function is purely real and the propagation of tensor gravitational waves in the medium is featured by dispersion only, reproducing the results obtained in the context of general relativity [24][25][26][27][28].By applying the same expansion of the denominator in (30) we obtain the real part of the dielectric function for scalar waves, which we report here for future convenience, from which it is possible to calculate the connected dispersion relation, i.e.
Now it is easy to verify that a subluminal phase velocity is ensured for all wavenumbers as long as the inequality is satisfied.Here we introduced, in analogy with the tensor case, the proper frequency ω 2 φ ≡ ακ ρ 6 .Then, for subluminal phase velocities, the pole appearing in integral (30) falls inside the domain of integration and the connected dielectric function acquires a non-null imaginary part.By making use of ( 28) we calculate the characteristic frequency which can be easily shown to be negative for any value of the wavenumber.The threshold effect introduced by inequality (35) clearly separates, for a fixed value of M , media that are able to damp scalar waves from those which leave the wave amplitude unaffected.Nevertheless, the same inequality can be seen as a criterion that selects the maximum mass allowed for a Horndeski scalar wave in order to be damped by a medium of assigned mass density ρ, restricting the range of variation of the parameters involved in the definition of M given in (8).
As reported in Figure 1, the phase velocity converges to the speed of light in the short-wavelength limit k → ∞ for any value of the mass M .Note that here and in the following, the overbar denotes quantities normalized by a factor ω −1 φ , e.g.k ≡ k ω φ and equivalently for the others.Particularly, for M 2 < 6γω 2 φ the function v p (k) is bounded, monotonically increasing with k and exhibiting a global minimum reached for vanishing wavenumbers, namely When the mode mass is in the range 6γω , is reached for a massless scalar wave.Moreover, scalar massless radiation from Horndeski theories will be damped by any material medium, given that in the case M = 0 inequality (35) is satisfied by any ρ = 0. Also when we look at the behavior of the characteristic frequency ω i (k) at fixed ζ, we find that the maximum magnitude of the damping is expected for vanishing mass at any k, as it can be clearly observed from the curves depicted in Figure 2.For growing masses in the range M 2 < 6γω 2 φ we expect decreasing absolute values of ω i (k).The minimum phase velocity v min defined in (37) is a striking peculiarity of gravitational Landau damping and it is absent in the analysis of the standard electromagnetic phenomenon.In the subsequent section, in which we study the interaction of Langmuir scalar modes with a tenuous beam of massive particles, we will illustrate the consequences of this important feature.

V. BEAM-PLASMA INSTABILITY FOR A MONOCHROMATIC TENUOUS BEAM
In this section, we analyze the selfconsistent scalar perturbations in a material medium which can be described as the superposition of two distinct populations of particles of the same mass m: an isotropic component at thermal equilibrium, which we refer to as "gravitational plasma", and a suprathermal anisotropic distribution, describing a beam of particles injected in the medium from an external source.As in the previous section, in the absence of gravitational waves the phase space properties of this gravitational beam-plasma system is depicted by a static and homogeneous total distribution function in which f P 0 represents the plasma contribution and f B 0 the beam component.It must be remarked that, with respect to the case presented before, the total background distribution function is not endowed with isotropy.The presence of the beam introduces a preferred direction.Of course, the background distribution function is normalized such that the integration on the entire phase space must return the total number of particles where the subscripts P and B stand again for plasma and beam contributions, respectively.As previously done, we assume that for any negative time no gravitational waves are present.For t ≥ 0 we perturb the medium with a scalar mode φ, in terms of which the external force dpi dt entering Vlasov equation has the form Following the steps already described we linearize the Vlasov equation ( 16) with respect to the distribution perturbation δf , obtaining and we readily obtain, in the Fourier-Laplace space, plane wave solutions traveling along the z axis, namely for the distribution perturbation, and for the scalar field.Then, we see that the presence of the beam alters the scalar field solution only through a modification in the dielectric response of the medium.Indeed, in the denominator of ( 43) together with the plasma function already defined, namely we have the appearance of an additional term, taking into account the beam contribution to the dispersion relation, i.e.
The plasma population will be again described through a Jüttner distribution as given in (29), where in this case we indicate the density of particles with the symbol n P .For what concerns the beam distribution we start by considering the simplest case, namely a monochromatic beam of the form in which p B and n B are the beam momentum and density of particles, respectively.Under this hypothesis, the term B can be exactly integrated and results in where the complex frequency ω and the plasma characteristic frequency ω φ were previously defined, and we have introduced the beam velocity and the densities ratio η ≡ n B n P .For what concerns the plasma term P , we first proceed by ignoring any imaginary contribution from the Landau pole, identifying this object with the real part of the dielectric function for the scalar mode obtained from the truncated series, as given in (33).Hence, the allowed modes for plasma disturbances in interaction with the beam are found from the condition P + B = 0, which in our case can be cast combining equations (33) and (47) in the form Solving this equation with respect to ω would return the exact dispersion relation, but it turns out that this task implies the treatment of a complete sixth degree polynomial.The problem can be simplified by making the assumption to deal with a tenuous beam, namely considering the ratio of the particles densities much smaller than unity, i.e. η 1.In this case the right-hand side of ( 48) is almost vanishing and we can obtain approximate solutions by investigating the cases in which the left-hand side results exactly null.We can enumerate three separate scenarios: 1.A first class of solutions is generated by the points (k, ω) which satisfy P (k, ω) = 0 and ω = v B k.These points correspond to the plasma dispersion relation reported in (34).
2. Secondly, we have the points (k, ω) pertaining to the beam dispersion relation ω = v B k which do not lay on the curve describing the plasma dispersion relation, i.e. points that satisfy P (k, v B k) = 0.
3. The third possibility is that there exist a number of degenerate points (k 0 , ω 0 ) which satisfy both dispersion relations, namely P (k 0 , ω 0 ) = 0 and ω 0 = v B k 0 .We find that this condition is fulfilled by a single wavenumber k 0 , which is determined by the following expression In order to ensure the reality of k 0 the following inequality must hold and this translates into a lower bound for the beam velocity, namely As anticipated, the minimum phase velocity v min plays an important role in the gravitational beamplasma interaction, namely it regulates the range of beam velocities that allow for the existence of a degenerate point of the two dispersion relations.
Now we focus our analysis on the properties of the dispersion relation in the proximity of the degenerate point.
The physical relevance of this family of solutions is given by the coexistence of both the plasma and the beam coupling to the oscillation mode.We proceed by expanding (48) up to third order in terms of the small quantities δω = ω − ω 0 and δk = k − k 0 .Thus, neglecting quartic terms, we obtain the following cubic equation where we defined By considering the form of the plasma dielectric function (33) and the expression of the degenerate point (k 0 , v B k 0 ), with k 0 given by ( 49), the parameters just defined have the following explicit expressions: For v B > v min and M 2 < 6γω 2 φ it turns out that the quantities P and C are negative, whereas Q is positive.It is found that equation ( 52) possesses a single real and a couple of complex conjugates solutions for any allowed value of the quantities involved in the definition of the parameters P , Q and C. Thus scalar radiation with wavenumber k ≈ k 0 and frequency ω ≈ ω 0 will be affected by both a shift in the phase, i.e. dispersion, and a modification in the amplitude, due to the imaginary parts of the complex solutions.In particular it is expected the arising of an instability region for wavenumbers around k 0 , caused by the presence of a dispersion relation with positive imaginary part.We denote this particular branch of the dispersion relation as δω (δk).We report the expression of its imaginary part in correspondence of the degenerate point, i.e. for δk = 0 Im (δω (0)) It must be remarked that this expression is exactly coincident with the analogous formula derived in the electromagnetic case in [11], except for the supplementary factor which is a genuine outcome of this gravitational version of the beam-plasma instability.In particular, the presence of this extra term implies that the amount of energy conveyed from the beam to the Langmuir scalar modes goes to zero as the beam velocity approaches the speed of light.For v B in the allowed range [v min , 1), it can be easily shown that G is a monotonic decreasing function and attains its maximum value when v B = v min , namely 1 − v 2 min 5 6 .Thus, given that the minimum value for v min is reached for vanishing mass of the scalar mode, we claim that the maximum instability is predicted in the case M = 0 and for a beam velocity equal to 3 − 1 2 , with a numerical value Im (δω (0)) It must be remarked that, despite the formal resemblance between formulas (55) and ( 57) with the analogous ones in the electromagnetic sector as given in [11], here is the gravitational coupling constant κ that appears in the definition of the proper frequency ω φ , affecting also the values of the degenerate wavenumber and frequency, namely k 0 and ω 0 respectively.Thus it is easy to figure out that, given the much smaller value of κ compared to the elec-tromagnetic coupling constant, in the gravitational version of the inverse Landau damping phenomenon we will deal with much smaller effects with respect to the electromagnetic analogue.In all cases, we postpone the discussion about the magnitude of the predicted effects to section VII, where we will provide quantitative estimates of the amplification phenomenon in a concrete physical scenario.In the subsequent section we will proceed to analyze the gravitational beam-plasma system employing numerical techniques, in order to enlarge the range of applicability of our predictions.

VI. NUMERICAL ANALYSIS
In this section we aim to extend our previous analysis, made assuming a monochromatic beam distribution, to more general and realistic scenarios.The description of the beam particles in terms of Dirac delta functions ( 46) is ideal.Even supposing a perfectly monochromatic mechanism for the generation of the supra-thermal component, these particles can subsequently interact with the surrounding medium leading to a spreading of the distribution function; alternatively, the dependence of the generating mechanism to some external phenomenon can lead to the same result.In order to investigate this aspect, we perform a new analysis of the total dispersion relation P + B = 0 in the presence of a more realistic beam distribution function f B 0 .However, even for reasonable assumptions like Gaussian or Cauchy distributions, the integral in (45) does not have an explicit expression and requires a fully numerical approach.To overcome this, in the following we will use a toy model which can be integrated explicitly, but which still retains some freedom in order to gain insights on the warm (non δ) beam scenario.For this reason, let us define the following trapezoidal distribution function: where H is the Heaviside step function, p B represents the mean beam momentum, and the free parameters B and b indicate the major and minor bases of the trapezoid, respectively.It is worth noting that choosing delta functions on the p 1 and p 2 momenta does not affect the obtained numerical results, since the imaginary contribution from the Landau pole comes from the beam-aligned direction only.
We adopt the ratio B/p B as a qualitative measure of the beam temperature, i.e. the beam spreading in momenta space, with the value 0 as a limiting case corresponding to a perfectly cold beam.In addition, the ratio b/B regulates the shape of the distribution and it can be tuned in order to analyze different scenarios, with the limiting cases of b/B = 1 and 0 corresponding to a boxshaped distribution and a triangular one, respectively.The choice of the parameters describing the trapezoid is made in order to obtain the best resemblance with a standard Gaussian distribution, according to the following criteria: i ) the distribution must be centered around the mean of the target Gaussian µ, for symmetry reasons, ii ) the distribution must have the same maximum value as the target Gaussian, and iii ) the quadratic difference between the distributions, integrated over the whole domain, must be minimal.By consequence of these requests, we obtain the correspondence of the trapezoid to a Gaussian distribution with parameters µ and σ by setting Concerning the background plasma, we still assume a Jüttner distribution, so that (33) holds, while the beam dielectric function can be derived from (45) performing suitable substitutions.In particular, it is possible to explicitly integrate over the transverse plane by writing the integral in terms of the relativistic velocities u i = p i /p 0 .Then, the beam dielectric function results in: where we defined The integral has the explicit solution and the integration boundaries are given in terms of the distribution parameters by We emphasize that the calculation of the dispersion relation from the equation P + B = 0 is turned from an integral problem into the solution of a transcendental equation involving logarithms, which can be solved rather easily through numerical methods.We also remark that in this numerical analysis we consider the frequency ω as a complex quantity ab initio.Therefore, when treating the total dielectric function we determine the curves ω r,i (k) by requiring that both its real and imaginary parts vanish simultaneously.We are interested in characterizing the branch of the dispersion relation with positive imaginary part (still assuming ω i ω r ), i.In the following figures, the barred quantities correspond to variables normalized to ω −1 φ , as stated in section IV.The curves depicted in Figure 4 are obtained as numerical solutions of the system of equations Re P + B = 0 and Im P + B = 0.More specifically, for each fixed value of the wavenumber k we search for the only couple of real numbers ωr ( k) and ωi ( k) (with ωi > 0) which guarantees the vanishing of the complex total dielectric function.The real part of ω (top panel of Figure 4) turns out to be mostly linear and this is consistent with the tenuous beam scenario here considered.Indeed, for η 1 the presence of the beam alters the dielectric response of the medium only as a small perturbation.As a result, the real part of the frequency, describing dispersion of the wavelengths within the medium, is essentially similar to the one described by formula (34), derived with analytical treatments.Concerning the imaginary part, it can be noted from the bottom panel of Figure 4 that the instability exists only for wavenumbers larger than a minimum value, and it reaches a maximum in close proximity of the critical wavenumber k0 calculated from (49) according to the chosen parameters, equal to 3.01.
In Figure 5, we show how the growth rate results non-null in a narrow range of phase velocities around v B = 0.8.In particular, the instability is active for velocities smaller than v B , corresponding to the positive slope of the beam distribution function, as is expected from the standard electromagnetic theory.In fact, the inverse Landau damping process depends on the sign of the derivative of the distribution function taken at the considered velocity.Thus, positive values of ωi can only be obtained for velocities smaller than v B , where the beam distribution function has a positive slope, as can be inferred from Figure 3. Now we perform three sets of computations with a running on different parameters: i ) the width of the distribution B/p B , i.e. the beam temperature, ii ) the mass of the gravitational wave M , iii ) the beam mean momentum p B .In Figure 6 we display the imaginary part FIG. 6. Imaginary part of ω( k) for different temperature values as reported in the legend.The blue dots are the numerical solution of ( 48) corresponding to a perfectly cold beam scenario.
of the frequency obtained with the same method as the previous figure, for increasing values of the ration B/p B , which is a quantitative measure of the beam temperature, i.e. its spreading in the velocity space.As shows in the figure, the increase of the beam temperature determines a decrease of the maximum allowed growth rate of the instability.This result is coherent with the standard electromagnetic analysis in [11], according to the decrease of the slope of the distribution function.Indeed, as previously stated, the magnitude of the effect is correlated to the derivative of the beam distribution function calculated at the considered phase velocity.Additionally we also remark that, with the increase of the beam temperature, the wavenumber at which the maximum growth rate is reached decreases and departs from the critical value k0 calculated for a perfectly cold beam through the analytical treatment.Furthermore, we stress that in the finite width case a maximum wavenumber k arises for the existence of the instability, and the allowed range shrinks with increasing temperature, as can be noticed from the B/p B = 5 • 10 −2 case depicted in figure .Now we investigate the role of the scalar mode mass M in determining the magnitude of the amplification effect, by calculating the growth rate ωi at increasing values of the mass, with the same technique previously outlined.The obtained values are plotted in Figure 7 as a function of the phase velocity v p in order to emphasize once again (see Figure 5) the localization of the effect in velocity space.We see that the instability growth rate increases as the mass of scalar mode decreases.In this respect, it is worth noting that this behavior is consistent with the analogous feature of the Landau damping case presented in Figure 2.This result has a relevant physical meaning because it allows the survival of the instability within the severe constraints obtained on the mass of the scalar mode from the observations of gravitational waves propagation [42].We remark that in the case of vanishing mass we expect the possibility of energy exchange between the scalar mode and any material medium, irrespectively of the proper frequency ω φ ∝ √ ρ, given that the condition (35) is identically satisfied in the massless scenario.To conclude our analysis we are interested in studying the dependence of the instability on the beam velocity, which in the finite temperature case is just the mean velocity calculated from the beam distribution function.As shown in formula (56) we expect the growth rate to decrease for beam velocities approaching the speed of light.This feature has a simple physical explanation, since the background particle distribution is nearly vanishing close to v = 1 and thus the beam interacts with a Langmuir mode supported by a small number of plasma particles.In order to test this analytical prediction we perform three runs of simulations at different temperatures.In each run we calculate ωi ( k) for a definite value of the beam velocity, localizing the maximum of the curve together with the corresponding phase velocity.Then varying the parameter v B we repeat the calculation.The data obtained are displayed in Figure 8, with the values of the beam velocity used in the calculation reported in caption.We observe that the numerical analysis follows the expected behavior, regardless of the beam temperature.Still, the aforementioned dependence of the instability from the beam temperature is confirmed, i.e. warmer beams attain smaller values of the growth rate.In conclusion we point out that all the tests performed through numerical techniques have confirmed the results outlined in the analytical section, while extending our predictions to the case of a finite beam temperature.Although the beam distribution is taken as a trapezoid for numerical convenience, this choice preserves the essential features allowing for the investigation of a more realistic scenario with respect to the perfectly cold case.However, despite the good accuracy offered by our simplified model in characterizing the linear approximation of the beam-plasma system, we remark that an eventual analysis of the instability evolution throughout its quasi-linear regime would require an even more realistic description of the beam contribution, leading to a fully numerical treatment, as argued at the beginning of this section.

VII. QUANTITATIVE ESTIMATES
In this section we aim to give some estimates of the predicted effect magnitude of scalar waves amplification due to the gravitational beam-plasma interaction.As explained in the previous sections, when we look at the imaginary part of the frequency ω i as a function of the scalar mode mass M , its maximum value is obtained in the massless limit M → 0. Therefore we will always consider a vanishing mass for the scalar mode, given also the fact that this scenario is strongly favored by recent observations.Moreover, as shown in Figure 6, warmer beams correspond to weaker effects.For this reason we will assume a perfectly cold beam, well described by a Delta distribution, by making use of the analytical formulas ( 49) and (55).Now a simple consideration holds: from a direct inspection of (49), which we report here in the case M = 0, namely the order of magnitude of the critical wavenumber k 0 is roughly given by the medium proper frequency ω φ ∝ √ ρ, being ρ the medium mass density.An appreciable deviation of k 0 from ω φ is predicted only in the case of ultra-relativistic beam v B 1.For this sort of beams, however, the size of ω i is strongly suppressed by the factor 1 − v 2 B 7 6 , as clearly displayed in (55).Hence, if we want to look at a value of the critical wavenumber in principle measurable with available technologies, like interferometers and pulsar-timing arrays, avoiding the case of ultra-relativistic beams for which ω i 0, we have to consider material media characterized by the highest values of mass density.Unfortunately, by considering the typical dark matter densities in a system like the Milky Way galaxy, the proper frequency ω φ results far too small.To give an example, by taking the average value of dark matter density measured in the Solar system as reported in [43], i.e. ρ = 0.36 GeV cm −3 , we obtain a proper frequency ω φ ≈ 10 −15 Hz, which is something like 5 orders of magnitude smaller than the smallest frequency measurable with pulsar-timing arrays [44].However, one of the most robust hypothesis on the dark matter distribution in the proximity of the galactic center is that the latter is well described by a highly peaked, cuspy profile [45][46][47].In addition to this, it has to be considered that a Kerr supermassive black hole is capable of hugely enhancing the dark matter density in a region of size L = 50 r S around itself, being r S the black hole Schwarzschild radius.Indeed, as shown in [48], a rotating black hole compatible with the one probably hosted in the center of our Galaxy [49] is able to raise the dark matter density in its proximity up to a value of about ρ = 10 18 GeV cm −3 .In such an extremely dense environment we calculate a proper frequency ω φ = 7 • 10 −7 Hz.Then, by assuming ζ 1 so that γ ≈ 1 and a beam velocity v B = 0.66 we obtain a critical wavenumber k 0 = 1.4 • 10 −6 Hz.We remark that the frequency associated to k 0 , i.e. ω 0 = v B k 0 = 0.9•10 −6 Hz, falls well inside the range of frequencies to which pulsar-timing arrays are sensitive [50,51].Now plugging these numbers into formula (55) yields an imaginary part of the frequency ω i = η 1 3 3.6 • 10 −7 Hz.We can give a rough estimate of the total time of interaction between the beam and the Langmuir modes by dividing the length of the system L by the phase velocity relative to the critical wavenumber k 0 , which is clearly equal to the beam velocity v B .By taking the black hole mass as 4 • 10 6 solar masses we calculate a total time of interaction ∆T = 3 • 10 3 s.Then, computing the product between ω i and ∆T we obtain an amplification factor Γ = η 1 3 10 −3 .Defining A in and A f in = e Γ A in as the amplitude of the scalar waves before and after the interaction with the gravitational plasma, we calculate a relative amplification Now let us spend a few words on the nature of the beam.As shown in [52], Kerr black holes are able to generate highly collimated dark matter beams thanks to the Penrose process.The density contrast with respect to the en-vironment on scales smaller than the parsec is expected to be roughly η = 10 −3 .Then, we obtain a relative amplification This value of relative amplification seems rather small to be measured with a pulsar-timing arrays detection, given the typical strains inside the sensitivity curves of these instruments [44].Nonetheless, it has to be remarked that we provided our estimate by considering the scenario offered by the supermassive black hole hosted in the Milky Way galactic center.An analogous calculation could be carried out by looking at the case of a supermassive black hole with mass of billions or tens of billions solar masses, lying in the typical mass interval of black holes located in Active Galactic Nuclei [53].For such massive central objects, the predicted amplification factor would result even larger given the typical density values expected in their vicinity.

VIII. CONCLUDING REMARKS
We analyzed the gravitational version of the wellknown process of the beam-plasma instability, starting from the results obtained in [18] regarding the Landau damping effect in the same context.The gravitational mode which is taken as subject of the interaction with fast flowing particles is the scalar mode of a Horndeski theory of gravity.This linear perturbation was thought here as the gravitational Langmuir wave living in the background plasma, modeled via its dielectric function, derived in [18].
Following the standard methodology introduced in [11], we pursued an analytical treatment by assuming a null-temperature beam morphology, well described by a Dirac delta function, searching for those points in the plane (k, ω) which simultaneously satisfy the two conditions ω = kv B and P (ω, k) = 0, where P is the gravitational dielectric function.This double constraint provided a critical value of the wavenumber k 0 , which had a crucial role in the characterization of the instability.Then, the obtained results were validated and extended via a numerical treatment employing a trapezoidal beam velocity distribution.
The main merit of our study is outlining the presence of a beam-plasma instability in the gravitational sector, described by a non-zero growth rate of the scalar Horndeski mode, i.e. a region of wavenumbers around k 0 for which ω i (k) > 0 strictly holds.Particularly, the maximum of the instability is located in the proximity of the critical wavenumber k 0 discussed above for sufficiently cold beam distributions.We observed a shift from this value when the beam temperature is increased, i.e for larger spread of the velocity distribution.
We also argued the existence of a threshold for the phase velocity of the wave in the plasma in the long wave-length limit, corresponding to c/ √ 3 in the case of vanishing mass of the scalar mode.Such a feature is extremely important from a phenomenological point of view, since the limits on the mode mass recently derived from gravitational wave observations as well as other methods [30][31][32][33][34][35], do not rule out the possibility of an efficient energy transfer from the fast particles to the longitudinal scalar perturbation.Moreover, we found that when the phase velocity approaches the speed of light, the growth rate tends to zero, since the plasma particle population is rarefied, due to the natural cutoff at superluminal velocities given by the Jüttner distribution used to model the background configuration.
The relevance of our analysis, in the perspective of future accurate detection of the linear gravitational modes, is that, if a Horndeski longitudinal scalar mode is really present in nature, then it could be searched with much more probability of success in those conditions in which it is enhanced by the interaction with a fast particle population.More specifically, within bounded gravitational systems, here dubbed gravitational plasmas, Langmuir modes exist and they are not significantly suppressed by the Landau damping effect, as discussed in [18].These modes can be unstable with respect to the interaction with a fast massive particle beam of external origin, leading to the growth of the amplitude up to a saturation value, fixed by the nonlinear dynamics of the system.In this respect, it remains open the important issue of discussing the nonlinear beam-plasma interaction for a Horndeski mode, extending to the gravitational sector the results in [12,13].

2 φ < M 2 < 9γω 2 φ
FIG. 1. Phase velocity as a function of the normalized wavenumber k for mass M of the scalar mode in the ranges shown in the legend.

FIG. 2 .
FIG. 2. Damping rate as a function of the wavenumber for values of the mass of the scalar mode shown in the legend.

FIG. 3 .
FIG. 3. (Top) The beam distribution function, marginalized to its p3 component, in momentum space, superimposed with a matching Gaussian distribution (dashed line).(Bottom) The single particle distribution function in velocity space F T OT 0 , as the sum of plasma (blue) and beam (red) distributions (here we choose η = 0.1 for illustrative purpose).
e. the one corresponding to inverse gravitational Landau damping, resulting in the growth of the wave during the initial time evolution, until nonlinear effects become important.Hence, in what follows we will focus on this root of the total dispersion relation only.As for the free parameters, we start from the following baseline scenario: m = ω φ , M 2 = γω 2 φ , ζ = 100, p B = 4/3 ω φ , B/p B = 10 −4 and b/B = 1/9.With this choice for the parameters values we calculate a beam velocity v B = 0.8.Additionally, all computations are run with a fixed value η = 10 −5 , which guarantees the validity of the weak damping condition.

FIG. 4 .
FIG. 4. Real (top panel) and imaginary (bottom panel) parts of the quantity ω( k), with parameters corresponding to the baseline scenario.

FIG. 7 .
FIG. 7. Growth rate ωi as a function of vp, for values of the mass of the scalar mode shown in the legend.