Collective modes of a collisional anisotropic quark-gluon plasma

In this paper we consider the collective modes of a momentum-space anisotropic quark-gluon plasma taking into account the effect of collisions between the plasma constituents. Our analysis is carried out using a collisional kernel of Bhatnagar-Gross-Krook form and extends prior analyses in the literature by considering all possible angles of propagation of the gluonic modes relative to the momentum-anisotropy axis. We extract both the stable and unstable modes as a function of the collision rate and confirm prior findings that gluonic unstable modes can be eliminated from the spectrum if the collision rate is sufficiently large. In addition, we discuss the conditions necessary for the existence of unstable modes and present evidence that unstable mode growth rates are maximal for modes with momentum along the anisotropy direction. Finally, we demonstrate that when there is a finite collisional rate, gluonic unstable modes are absent from the spectrum at both small and large momentum anisotropy. These results pave the way for understanding the impact of collisions on a variety of non-equilibrium quark-gluon plasma observables.


I. INTRODUCTION
At sufficiently high temperatures nuclear matter undergoes a phase transition from hadronic bound states to liberated quark and gluon degrees of freedom. This phase, called the quark-gluon plasma (QGP), is expected to be weakly-coupled at very high temperatures due to the asymptotic freedom of quantum chromodynamics (QCD). Despite the fact that perturbative methods can be applied, in order to systematically understand the physics of a high-temperature QGP, one must resum a class of graphs called hard thermal loops (HTLs) [1][2][3][4][5]. The HTL dressing of the quark and gluon propagators introduces new mass scales that are proportional to the temperature. These temperature-dependent masses can be obtained from the static limit of the quark and gluon self-energies, which are non-trivial functions of the temperature and the four-momentum of the excitations [1,2]. Computation of these selfenergies allows one to determine the collective modes of the system, which in the case of gluons can be determined by finding the poles in the resummed propagators. In the case of an equilibrium QGP, one finds that there are two gluonic collective modes corresponding to transverse and longitudinal polarizations [6]. In addition, one finds that a non-trivial cut appears in the self-energies related to the physics of Landau damping of plasma oscillations.
Such resummations are necessary ingredients for the computation of both static and dynamical properties of the QGP. See Ref. [7] for a recent review of perturbative resummation in QCD. Despite the broad applicability of these results, they are limited in their assumption that collisions among plasma constituents are ignorable. In the high temperature limit this is due to the fact that the time scale between hard collisions (large momentum exchange), given by t hard T ∼ 1/α 2 s log α −1 s , is much longer than the time-scale for the soft momentum exchanges, given by t soft T ∼ 1/α s log α −1 s [8]. As the temperature is decreased and α s increases in magnitude it becomes necessary to account for both hard and soft collisions. In thermal equilibrium for a variety of observables, e.g., jet energy loss, the systematic inclusion of collisions can be accomplished using the effective kinetic theory of Arnold, Moore, and Yaffe [9][10][11].
In the case of a non-equilibrium QGP it is still possible to make studies of the collective modes of the system, however, the situation is complicated by the presence of unstable gluonic modes related to the chromo-Weibel instability [12][13][14][15][16]. For a recent review of plasma instabilities in the QGP, see Ref. [17]. In the collisionless limit such instabilities generically arise if there is any degree of momentum anisotropy in the underlying gluon one-particle distribution function, with the growth rates of the unstable modes depending on the temperature, degree of momentum anisotropy, and the angle at which the mode is propagating relative to the anisotropy axes. The existence of the chromo-Weibel instability makes it difficult to understand the late time dynamics of momentum-anisotropic plasmas. For example, in the case of the computation of the imaginary part of the heavy-quark potential, one finds that the chromo-Weibel instability results in a pinch singularity that causes the imaginary part of the integral defining the heavy-quark potential to diverge [18,19]. In Ref. [18] the analysis was performed assuming small anisotropy, while in Ref. [19] this finding was extended to all anisotropy strengths. The pinch singularities found were a direct result of the appearance of unstable modes in the spectrum and, in the static limit, there being modes with negative squared masses. A similar problem will occur in the perturbative calculation of many key quantities used in QGP phenomenology, e.g., jet energy loss and momentum broadening, if the perturbative calculations are carried to sufficiently high order in the gauge coupling.
In the case of the heavy-quark potential the issue is manifest already at leading order, so this makes it an ideal candidate for understanding the role plasma instabilities play. As discussed in Refs. [18][19][20][21] the real and imaginary parts of the heavy-quark potential can be obtained from three-momentum integrations involving the retarded and Feynman propagators, respectively. As a result, one needs to know the form of the resummed propagators and hence self-energies for arbitrary angle of propagation with respect to the anisotropy axes. In this paper, we study the gluonic collective modes of a collisional plasma subject to a momentum space anisotropy of spheroidal form, i.e., where T and L indicate momentum perpendicular and parallel to a given anisotropy axis denoted by a unit vector n, respectively. The degree of momentum-space anisotropy is quantified by an adjustable parameter ξ in the range −1 < ξ < ∞. In addition, λ is a scale which specifies the typical momentum of the particles in the plasma and should be understood as the temperature only in the thermal equilibrium limit [14,16]. For a spheroidal deformation, there is only one anisotropy axis and the gluon polarization tensor depends on the angle of propagation of the modes relative to n. In the collisionless limit, the collective modes of a QGP with a spheroidal momentum anisotropy were determined completely in Refs. [14,16]. In this limit, one finds two unstable modes whose growth depends on the strength of the anisotropy, ξ, and the angle of propagation θ of the mode relative to n.
In order to study the collective modes of the QGP beyond the collisionless limit, in Ref. [22] the authors computed the collective modes of an isotropic thermal QGP by making use of a collisional kernel of Bhatnagar-Gross-Krook (BGK) form. They found that a non-zero collision rate ν = 0 pushed the cut associated with Landau damping into the lower half plane, induced damping of the stable collective modes, and caused one of the two stable collective modes to terminate at finite momentum. 1 The work of Carrington et al in Ref. [22] was extended to the case of a momentum-anisotropic QGP in Ref. [8], wherein the authors considered a spheroidally deformed QGP, but only considered the special case that the momentum of the gluon was along the momentum anisotropy direction n. In this special case, it was possible to obtain analytic results for all four gluon polarization functions. The authors of [8] found that it was possible to eliminate the unstable modes from the spectrum if the collision rate was sufficiently large.
In this paper we extend the analysis of Ref. [8] to include all possible gluon propagation angles and general anisotropy strength ξ. 2 For general anisotropy strength and propagation angle it is not possible to obtain purely analytic results, however, the problem can be reduced to numerically evaluating a set of one-dimensional integrals. To prove this, we decompose the polarization tensor into a set of five independent structure functions. We then determine the full spectrum of stable and unstable modes as a function of both the collision rate ν and angle of progagation parameterized by t = cos θ with 0 ≤ t ≤ 1. Our findings confirm those found in the case t = 1 [8] and we find that for all ν > 0 the unstable growth rate is maximal at t = 1 (θ = 0 with respect to the anisotropy axis n). Our results set the stage for the computation of the real and imaginary parts of the heavy-quark potential in a collisional and momentum-space anisotropic QGP. Beyond this, our results can be used in the computation of a large set of inputs necessary for understanding the role of momentum-space anisotropies in the QGP at intermediate couplings, e.g., jet energy loss, shear viscosity at higher orders in α s , etc.
The structure of our paper is as follows. In Sec. II we review how one obtains the gluon polarization tensor from kinetic theory. In Sec. III we detail how to obtain the five structure functions necessary to describe a collisional and momentum-space anisotropic plasma. In this section we also discuss the static limit of the structure functions and their relation to the five non-trivial mass scales appearing in the problem. In Sec. IV we present our results for the dispersion relations of the (damped) stable modes in the system for various gluon propagation angles and collisional rates. In Sec. V we present our results for the dispersion relations of the unstable modes. Finally, in Sec. VI we present our conclusions and an outlook for the future. In App. A, we present some limiting cases in which it is possible to obtain analytic expressions for the five structure functions and static masses. Finally, in App. B we present results obtained when holding the energy density of the system fixed when varying the plasma anisotropy.

II. THE GLUON SELF-ENERGY FROM KINETIC THEORY
In the high temperature phase of QCD, the gauge coupling constant g becomes small enough and, therefore, there exists a hierarchy of scales and one can construct effective theories at various scales by eliminating the degrees of freedom at higher scales. For example, the collective excitations, which typically develop at a momentum scale ∼ gT , are well separated from the typical energy of hard partons and can be described by kinetic equations of the Boltzmann-Vlasov type. In the kinetic theory, the distribution of hard particles in the QCD plasma is described by the gauge covariant Wigner functions n i (p, X) where the superscript i denotes the species of the plasma constituents, i.e., quarks, antiquarks and gluons, which are treated as massless particles and assumed to satisfy the mass-shell constraint. Expanding the Wigner function around the color neutral background fields n i (p), namely, n i (p, X) = n i (p) + δn i (p, X), the linearized kinetic equations read where V = (1, v) with v ≡ p/p and θ q = θ g = 1, θq = −1. In the above equation, we have already taken into account the weak-coupling limit and neglected terms of subleading order in g. Thus, the theory becomes effectively Abelian. There is no coupling among different color channels in Eq. (2) and the field strength tensor is simply given by F µν = ∂ µ A ν − ∂ ν A µ since the term involving a gauge field commutator can be dropped consistently. In addition, the distribution function f i (p) and its color fluctuation δf i a (p, X) are related to the Wigner function via and where t a and T a are the SU(N c ) group generators in the fundamental and adjoint representations, respectively. On the right hand side of Eq. (2), we introduce a BGK-type collision term, which is given by with f i a (p, X) = f i (p) + δf i a (p, X). The above collision kernel ensures an instantaneously conserved number of particles N i a (X) = p f i a (p, X). Accordingly, the particle number in equilibrium is given by N i eq = p f i eq (p). In addition, the collision rate denoted by ν is inversely proportional to the equilibration time of the hot QCD plasma under collisions between the hard partons. In the remainder of this work, the collision rate ν is taken to be a free parameter which has no dependence on the momentum and particle species. Determining the magnitude of ν in the ultrarelativistic heavy-ion experiments is still an open question. Estimates from Ref. [8] place ν in the range ν ∼ 0.1 − 0.2 m D , however, it is possible that the collision rate is even larger than this. Herein, we present results for small ν and postpone the consideration of large ν to a forthcoming paper.
It is more convenient to solve the linearized kinetic equations for the fluctuation δf i a (p, X) in momentum space which leads to the result for the linearized induced current by each particle species. The total induced current J µ ind a (K) with K = (ω, k) is simply a sum of the contributions from quarks, antiquarks, and gluons. Furthermore, the gluon self-energy can be obtained by functional differentiation of J µ ind a (K) with respect to the gauge field A b ν (K) and has already been computed in Ref. [8]. Including the collision kernel as given in Eq. (5), the gluon self-energy Π µν ab (K) takes the following non-trivial form In the above equation, p ≡ d 3 p/(2π) 3 and v ′ = p ′ /p ′ . For convenience, we introduce the dimensionless collision rateν = ν/k and defineK = (ω,k) withω = ω/k andk = k/k. In addition, the function W(ω,ν) reads Inclusion of the collision kernel doesn't break the transversality of the gluons self-energy and one can show that K µ Π µν = K ν Π µν = 0. As the Lorentz components of the gluons self-energy are not all independent, it is sufficient to consider only the spatial components and the temporal components can be simply obtained by using the relations Π 0i (K) = k l Π li (K)/ω and Π 00 (K) = k i k j Π ij (K)/ω 2 . Since Eq. (6) is diagonal in color, color indices will be suppressed in the following.
In the above equation, f (p) = 2N c f g (p) + N f [f q (p) + fq(p)] and the distribution functions of the hard partons are completely arbitrary up to the requirement that the distribution function used results in convergent integrals for, e.g., the number density, energy density, gluon self-energy, etc. This sets constraints on the behavior of f (p) at small and large momenta. From here on, an anisotropic parton distribution function as specified in Eq. (1) will be adopted. It can be obtained by stretching or squeezing an isotropic distribution f iso along the direction of anisotropy n which is conventionally taken to be parallel to the beam-line direction. Due to early-time free-streaming of the QGP, it is expected that ξ > 0, corresponding to an oblate momentum-space distribution. As a result, an oblate distribution with ξ > 0 is the case most relevant to high-energy heavy-ion experiments and hence will be considered in the following. Plugging Eq. (1) into Eq. (6) and performing the radial part of the integrations, the gluon self-energy using the above distribution reads where It is important to point out that the above gluon self-energy is not symmetric in Lorentz indices due to the appearance of the BGK-type collision term, 3 namely Π ij (K) = Π ji (K). The same is also found in a magnetized plasma due to the lack of time-reversal symmetry [24]. Thus, a symmetric tensor basis as developed for an anisotropic system in the collisionless limit doesn't apply to the decomposition of Eq. (8). Here, we adopt a more general tensor basis and decompose the gluon self-energy into five structure functions as where the tensor basis reads In the above equation,ñ i ≡ A il n l , which is orthogonal tok i . As compared to the tensor basis used in Refs. [8,14], the only difference is that the structure functions forñ ikj andk iñj are not identical any more. The determination of the five structure functions can be carried out based on the following contractionŝ Clearly, if the gluon self-energy is symmetric, the last two contractions lead to the same result and δ = ρ. The explicit results for the five structure functions in the above decomposition will be presented in the next section. Given the gluon self-energy, one can further compute the resummed gluon propagator ∆ ij (K) by using the Dyson-Schwinger equation. In temporal axial gauge, the inverse propagator is written as Upon inversion, we find that The dispersion relations are obtained by finding the poles of the propagator ∆ ij (K), namely, Hence, we find three separate equations which determine the so-called α-mode, G + -mode and G − -mode, respectively.

III. STRUCTURE FUNCTIONS AND MASS SCALES
To evaluate the structure functions, we need a master integral of the form In the above equation, n is an integer and n = 0, 1, 2, 3 are relevant in the following evaluations. The propagation angle θ parameterized by t = cos θ corresponds to the angle between k and n. In addition, we assume that n is parallel to the z-axis while k lies in the y-z plane. More generally, k y should be replaced by k ⊥ . The integration variable x is cosine of the polar angle between p and n. In general, ω is complex-valued. Then we can express the structure functions as the following 4 4 An overall factor m 2 D is omitted in these expressions. The same is done for the mass scales. where and z =ω + iν.
It is interesting to consider the static limit of the structure functions which defines the corresponding mass scales as the following, Taking the limitω → 0, the master integral in Eq. (19) reduces tõ Its imaginary part is an even function of x while the real part is odd in x. Therefore, the above integral can be written as where and When n is even, Eq. (29) is pure imaginary, while it is real-valued when n is odd. In the above equation, the sign function sign(k z ) or sign(t) has been dropped. In fact, one can only consider the region 0 ≤ t ≤ 1 because as we will see below, the mass scales m 2 α , m 2 β and m 2 γ are symmetric under t → −t. Although m 2 δ and m 2 ρ are odd functions of t, it is their product m 2 δ m 2 ρ that is relevant in our consideration, which is also even in t.
In terms of the integral defined in Eq. (29), one can obtain the five mass scales. We find The factor W(0,ν) in m 2 β (ξ, t,ν) is real-valued and can be expressed as 1 −ν arctan (1/ν). Therefore, the above five mass scales are all real-valued. For later use, we also define the following two mass scales with We point out that whenω is imaginary-valued, i.e.,ω = iΓ ≡ iΓ/k, the structure functions can be expressed in terms of the integrals defined in Eq. (29) which are more efficient than using Eq. (19) when performing the numerical evaluation. It is straightforward to show that I [n] (ξ,Γ, t,ν) =Ĩ [n] (ξ, t,Γ +ν), therefore, the five structure functions are all real-valued for imaginary-valuedω. In general, the above structure functions and the corresponding mass scales must be evaluated numerically. However, it is possible to obtain analytical results in some limiting cases, which we discuss in App. A.

IV. DISPERSION RELATIONS FOR THE STABLE MODES
The solutions ω(k) of Eq. (15) can be obtained numerically, which determine the dispersion relations for the collective modes. The corresponding results for the stable modes are shown in Figs. 1-3 where all quantities are given in units of m g ≡ m D / √ 3. As found in previous studies, the three stable modes in the collisionless limit have poles at real-valued ω > k. However, with non-zero collision rate, these types of modes correspond to the complex-valued solutions ω(k) which are damped in time due to a negative imaginary part. 5 According to our results, the α-mode behaves very similarly to the G + -mode for large wave number k while for small k, it is more like the G − -mode. In particular, the α-mode is identical to either the G − -mode or the G + -mode at t = 1, thus, only two distinguishable modes exist, which were called α-mode and β-mode in Ref. [8]. It is interesting to consider the limit k → 0 where the values of ω in the dispersion relation can be obtained analytically. We find that all the structure functions are finite in this limit and α(k → 0) becomes angle-independent. Furthermore, although the other four structure functions still have a t-dependence, Ω 2 + and Ω 2 − as defined in Eq. (16) don't depend on t. Explicitly, we have As a result, the dispersion relations at vanishing wave number are determined by only the anisotropy parameter ξ and the collision rate ν. The corresponding values of ω α , ω + , and ω − in this limit are given by 5 Following the terminology in the previous works, we still call them stable modes.
where the imaginary part of ω equals to −ν/2 independent of the anisotropy parameter ξ. We note that the above complex solutions for the dispersion relations in Eq. (18) exist when the collision rate is not very large so that the square roots in the above equation are real. For the isotropic case, it requires ν/m D < 2/ √ 3. With increasing ξ, the upper bound of ν/m D will be reduced.
It is also worth noting that the G − -mode exhibits some new features in the collisional cases. It could become space-like (Re(ω) < k) at finite collision rate where a strong damping rate is also found. In addition, the G − -mode terminates at certain wave number k where the magnitude of its imaginary part approaches the collision rate. As mentioned in footnote 1, the reason for the termination of the stable G − mode is that the stable mode solution moves from the physical Reimann sheet to a higher Reimann sheet when it passes through the cut that exists between −1 < Re(ω)/k < 1 and Im(ω) = −ν. This is most easy to see in the case t = 1 where analytic forms for all structure functions can be obtained (see App. A and Eqs. (39) and (40) of Ref. [8]). In Ref. [8] the authors understood the termination of this mode by tracking the movement of the solution from the physical Reimann sheet to higher Reimann sheets associated with the logarithm appearing in the structure functions (see Eq. (41) of Ref. [8]). We refer the reader to Figs. 5 and 6 of Ref. [8] and the associated discussion of this figure. When t = 1, the cut remains in the same location, however, one must carefully analyze the discontinuities in the integrand of the master integral I [n] defined in Eq. (19) to reach the same conclusion. We note that the termination of this mode also occurs in the isotropic case with a finite collision rate [22].
According to Figs. 2 and 3, at a given anisotropy ξ, Im(ω − ) drops to −ν more quickly for larger t. In general, determining the value of k at which the G − -mode terminates involves a complicated numerical evaluation. In the special case t = 0 and t = 1, we show in Fig. 4 how the values of k ⋆ /m g as well as the ratio Re(ω ⋆ − )/k ⋆ at the termination point change with the collision rate ν. Notice that for clarity, the values of k and Re(ω − ) at the termination points are denoted by k ⋆ and Re(ω ⋆ − ), respectively. Our results show that Re(ω ⋆ − )/k ⋆ deviates from the light cone (which corresponds to the line Re(ω ⋆ − )/k ⋆ = 1 in this figure) more significantly when k is parallel to n. For a fixed anisotropy, the deviation becomes larger with increasing ν. However, as ξ increases, our numerical results indicate that the deviation can become either larger or smaller depending on the propagation angle. Finally, we note that as a function of ν, k ⋆ /m g can have a non-monotonic behavior depending on the degree of momentum-space anisotropy. In order to see the parameter dependence of the dispersion relations, we consider the α-mode as an example and show the corresponding t-, ν-and ξ-dependence in Fig. 5. As one can see from this figure, the real part of the dispersion relations is almost independent on the propagation angle, but becomes sensitive to the degree of anisotropy. On the other hand, the collision rate has a visible influence only when ν is very large. The dependence of the imaginary part in the dispersion relations on these parameters is more significant. Notice that similar behaviors are observed in the stable G + -and G − -modes. Although we only present results at a set of specific values of t, ν, and ξ, we internally considered many cases which are not listed here but support the same conclusions.

V. DISPERSION RELATIONS FOR THE UNSTABLE MODES
We can study the unstable modes by considering an imaginary-valued energy ω = iΓ. Accordingly, the dispersion relations in Eq. (15) become Since Ω 2 + (iΓ + ) is numerically found to be positive for Γ + > 0, there is only one unstable G-mode. Some examples can be found in Fig. 6. Notice that in the special case where t → 1, one can prove that there is only one unstable  mode [8] because the unstable α-mode becomes identical to the unstable G-mode.
The existence of the unstable modes depends on the specific values of ξ, t andν = ν/m D . As shown in Fig. 6, with fixed ξ and t, the unstable modes can only survive in certain region ofν. Thus, it is interesting to find the conditions for k max (ξ,ν, t) → 0 which indicates the absence of an unstable mode. Here, k max corresponds to a non-zero value of the wave number k at which the unstable mode terminates, i.e., Γ → 0. By setting Γ = 0 in the dispersion relation, k max (ξ,ν, t) for the unstable α-mode is determined by the following equation Besides a trivial solution k max = 0, we are actually interested in the non-zero solution k max (ξ,ν, t), which can be obtained by numerically solving the above equation. In Fig. 8, we show how the dimensionlessk max = k max /m D depends on its variables. On the other hand, Eq. (41) can be rewritten as  with x ≡ ν/k max . This equation is well defined because m 2 α (ξ, t, x) is numerically found to be negative. In the collisionless limit the value of k max can be simply determined by k max (ν → 0) = −m 2 α (ξ, t, x → 0), where the corresponding mass scales have been analytically computed in Ref. [20]. Notice that m 2 α (ξ, t, x → 0) is negative and vanishes in the small and large ξ limit, therefore, the unstable α-mode exists for any non-zero and finite anisotropy in the collisionless limit. 6 On the other hand, according to Fig. 8, this is no longer true when the collision kernel is included.
In general, one can expect that the instability of an anisotropic system will be eventually eliminated by continuously increasing the collision rate. As a result, the survival of the unstable α-mode requiresν <ν(x → ∞) which is the collision rate at vanishing k max . Therefore, it becomes very important to analyze the behavior ofν(x → ∞). To do so, we Taylor expand the mass scale m 2 α (ξ, t, x) for infinitely large x and in the leading-order approximation, the result readsν where It can be shown that the function f (ξ) is positive and vanishes as ξ → 0 or ξ → ∞. In addition, the maximum of f (ξ) is numerically found to be ∼ 0.0562 which is located at ξ ≈ 6.40. Given the above results, we can obtain the conditions for the existence of the unstable α-mode based on Eq. (43) which are: • with fixed ξ and t, the collision rate should satisfyν < t f (ξ). The largest possible value for the collision rate is given by f (ξ) max ≈ 0.237.
• with fixed ξ andν, the propagation angle should satisfy t >ν/ f (ξ). For small enoughν, the propagation direction of the unstable α-mode becomes (almost) perpendicular to the direction of anisotropy which indicates t → 0.
• with fixedν and t, the anisotropy parameter should satisfy ξ min < ξ < ξ max where ξ min and ξ max denote the two solutions of Eq. (43).
Theν/t-dependence of ξ min and ξ max has to be determined numerically, which is given in Fig. 9. In addition, the two solutions become identical atν/t ≈ 0.237. It is possible to obtain analytical results for ξ min and ξ max for extremely small and large anisotropies, respectively. Explicitly, we have These two equations are good approximations to the exact results provided that the ratioν/t is small. They also confirm that the unstable α-mode exists for any non-zero and finite anisotropy in the collisionless limitν → 0 when t = 0. The above discussions provide a simple way to locate the zeros of k max (ξ,ν, t) as shown in Fig. 8 which actually correspond to the critical condition for the existence of the unstable α-mode. Thus, given a set of parameters ξ, t and ν, one can easily determine whether or not the unstable α-mode exists.
It is also important to study the maximal growth rate Γ max in the dispersion relations of the unstable modes. As shown in Fig. 10, Γ max depends on the collision rate, propagation angle, and the degree of anisotropy, which exhibits very similar behaviors when compared with the dependence on k max . In practice, it is of particular interest to find the maximum possible value for Γ max . According to our results, at a given ξ and t, the maximum of Γ max appears in the collisionless limit. In other words, the collision effect reduces Γ max . In addition, for fixed ξ andν, we find that Γ max increases with increasing t and thus, confirms that the growth rate of the filamentation instability is largest when the wave vector k is parallel to the direction of the anisotropy, i.e., t = 1. Finally, we note that, as can be seen in the rightmost panel of Fig. 10, Γ max is a non-monotonic function of the anisotropy parameter ξ and the peak position needs to be determined numerically.
Next, we discuss the unstable G-mode which exhibits some new and complicated features as compared to the unstable α-mode. We start with the dispersion relation in the limit Γ → 0 which can be written as Since the mass scale m 2 + (ξ, t, ν/k max ) is numerically found to be always positive, one can only consider the second term in the above equation. Following the discussion of the unstable α-mode, a similar expression for the collision rate ν can be obtained as Unlike the negative mass scale m 2 α (ξ, t, x), depending on the specific values of the variables, m 2 − (ξ, t, x) can be positive which leads to the absence of the unstable G-mode as well as an ill-definedν(x) in Eq. (47). Therefore, the above equation applies only for m 2 − (ξ, t, x) ≤ 0.  Following a similar discussion on the unstable α-mode, we can expand m 2 − (ξ, t, x) for large x. When keeping the leading order contribution, Eq. (47) becomesν where g(ξ, t) = g 2 (ξ)(2t 2 − 1) + g 1 (ξ) with Numerically, we find that m 2 − (ξ, t, x) is positive if it approaches to zero from above when x → ∞. Consequently, only positive g(ξ, t) is relevant for studying the unstable mode since g(ξ, t) and m 2 − (ξ, t, x → ∞) have opposite signs. Based on Eq. (48), the following conclusions can be drawn.
• with fixed ξ and t, the collision rate should satisfyν < g(ξ, t). Because both g 1 (ξ) and g 2 (ξ) are positive, the function g(ξ, t) increases with increasing t at a given ξ. Thus, the largest possible value for the collision rate is given by f (ξ) max ≈ 0.237 which coincides with our finding for the unstable α-mode.
• with fixed ξ andν, the propagation angle should satisfy t > t 0 (ξ,ν) where t 0 denotes to the solution of Eq. (48) and is given by In the above equation, the values of ξ andν cannot be arbitrary and a prerequisiteν < g(ξ, 1) must be met. Furthermore, at a given anisotropy ξ, t 0 (ξ,ν) increases with increasingν and the minimum of t 0 (ξ,ν) is found to be 1/ √ 3 ≈ 0.577 whenν ≪ g(ξ, 1) and ξ is infinitely large. 7 Accordingly, the propagation of an unstable G-mode is forbidden when θ is larger than arccos(1/ √ 3) ∼ 0.3 π.
• with fixedν and t, the anisotropy parameter should satisfy ξ min < ξ < ξ max where ξ min and ξ max denote the two solutions of Eq. (48).
As shown in Figs. 8 and 11, the most obvious difference between the two kinds of the unstable modes lies in thẽ ν-dependence of k max . In general, with fixed ξ and t, k max for the unstable G-mode is a non-monotonic function of ν and the peak position moves toν = 0 very quickly when increasing t. In the limit t → 1, the unstable G-mode becomes identical to the unstable α-mode, and thus k max decreases with increasingν. In particular, the function k max (ν) can even pass through the origin. According to Eq. (47), this happens due to the vanishing m 2 − (ξ, t, x) at non-zero x where bothν and k max are infinitely small but their ratio x is finite.
It is worth mentioning that the peak values of k max correspond to the minima of m 2 − (ξ, t, x). According to Fig. 7, the minima appear at some non-zero x for m 2 − (ξ, t, x), nevertheless, m 2 α (ξ, t, x) becomes smallest as x → 0. This explains the differentν-dependences of k max as observed in our results. Another interesting finding is that the unstable G-mode survives in a wider region of t in the collisional case which can be directly seen from Fig. 7. For example, given ξ = 10 and t = 0.68, m 2 − (ξ, t, x → 0) is positive and thus, no unstable mode exists in the collisionless limit where x → 0. However, a negative m 2 − (ξ, t, x) shows up at some finite x which indicates the existence of the unstable modes. As a result, one can expect that for certain propagation angles, the unstable G-mode, which is absent in the collisionless limit, would become unstable at finite collision rate. It can be also understood by considering the sign of m 2 − (ξ, t, x) in the small and large x limits. We find that at a fixed anisotropy ξ, a negative mass scale exists above a critical value of t, which are denoted as t c (ξ, x → 0) and t c (ξ, x → ∞) for m 2 − (ξ, t, x → 0) and m 2 − (ξ, t, x → ∞), respectively. According to Fig. 12, the critical value of t in the collisionless limit is always larger than that for infinitely large x. Thus, for any given value of t which is between the dashed and solid curves in this figure, there are unstable G-modes only for non-zero collision rate. Furthermore, the minimum possible value of t for finding an unstable G-mode when x → 0 is numerically found to be ∼ 0.689 where ξ ≈ 5.64, while it can be lowered to 1/ √ 3 ≈ 0.577 as ξ → ∞ if the collision term is included. On the other hand, as a function of t or ξ, k max shows a similar behavior for the two different kinds of unstable modes. However, recall that the existence of a non-zero lower bound t 0 (ξ,ν) for finding an unstable G-mode differs from the unstable α-mode where t can take any value between 0 and 1 provided thatν/t < 0.237 can be satisfied. Given the fact that the values of k max as t → 1 are the same for both unstable α-mode and G-mode, a larger increasing rate of k max can be expected for the latter. In addition, unlike the unstable α-mode, ξ min and ξ max of Eq. (48) don't simply depend on the ratioν/t. In Fig. 9, bothν-and t-dependence of these two solutions are presented. With fixed t, the two solutions become identical whenν equals the maximum of g(ξ, t). On the other hand, if the collision ratẽ ν is fixed, the equality of the two appears when the propagation angle t in Eq. (48) becomes smallest when varying ξ. Analytical results can be also obtained for ξ min and ξ max for extremely small and large anisotropies, respectively, which read These two equations are good approximations to the exact results whenν is small and t is close to 1. They also indicate that for t > 1/ √ 2 ≈ 0.707, the unstable G-mode can exist for any non-zero and finite ξ in the caseν → 0. In Fig. 13, we show the behavior the maximal growth rate Γ max for the unstable G-mode. Comparing with Fig. 10, a non-monotonic dependence onν is observed here, which is in accordance with theν-dependence of k max . Finally, as a function of t and ξ, Γ max behaves similarly for the two different kinds of unstable modes.

VI. CONCLUSIONS
In this paper, we have studied the collective modes of an anisotropic QCD plasma with a BGK collision kernel describing the equilibration of the plasma. This is an extension of a previous work where only a single propagation direction of the collective modes was considered. By solving the linearized kinetic equations of the Boltzmann-Vlasov type, we rederived the gluon self-energy. It turned out that the BGK collision kernel could break the symmetry of the gluon self-energy, namely, Π µν (K) = Π νµ (K) if the plasma has an anisotropic hard particle distribution in momentum space, ξ = 0. As a result, one needs a more general tensor basis to decompose the gluon self-energy for such an anisotropic system, which leads to five structure functions in the decomposition as shown in Eq. (10). This corrects the previously obtained result [8] where the gluon self-energy was taken to be symmetric in the Lorentz indices. It should be noted that in Ref. [8] the wave vector k was assumed to be directed along the anisotropy direction n, i.e., t = 1 and in this special case, the gluon self-energy obtained therein coincides with Eq. (10) where only two structure functions α and β are relevant. Consequently, the resulting collective modes with this specified propagation direction are indeed correct, however, for t = 1, five independent structure functions are needed.
We further determined the resummed gluon propagator, based on which we analyzed the collective modes by finding the corresponding poles. Dispersion relations for both the stable and unstable modes were presented. With the BGK collision kernel, the solutions ω(k) become complex valued for the three stable modes. It was found that Im(ω) equals to −ν/2 as k → 0 and for the G − -mode, it eventually drops to −ν at some finite k where this mode terminates. The value of k at the termination points turned out to be sensitive to the propagation angle, however, they also had a nontrivial dependence on the collision rate ν and the anisotropy parameter ξ. In addition, the k-dependence of Re(ω) seemed very similar as that shown in Ref. [8] which was actually consistent with our finding that Re(ω) had a very weak dependence on the propagation angle. Here, the most notable behavior was the spacelike dispersion in the G − -mode at large k. Our numerical results suggested that the ratio Re(ω − )/k at the termination point deviated from the light cone more significantly if the propagation is parallel to the anisotropy direction n. It was also found that such a deviation could increase with increasing ν at a given anisotropy ξ.
Our results also confirmed that when propagating along the direction of anisotropy, the corresponding modes would become most unstable even in the collisional case. In other words, the largest possible value of the maximal growth rate Γ max appeared at t = 1. On the other hand, different from what can be expected intuitively, namely that increasing the collision rate slows down the rate of growth of the unstable modes, we found that the maximal growth rate was not always a monotonically decreasing function of the collision rate. Exceptions were found in the unstable G-mode in the small ν region. Especially, for certain propagation angles, the unstable G-mode which is absent in the collisionless limit can be activated by introducing a small collision rate. We also carried out a systematic discussion of the parameter dependences of two key quantities associated with the unstable modes, namely k max and Γ max which were introduced in Sec. V. The former could be used to determine whether or not the unstable modes exist, while the latter described the degree of the instability. In general, no unstable mode could survive when the collision rate ν > 0.237 m D and the elimination of the unstable modes occurs at both small and large anisotropies, depending on the values of ν and t. Furthermore, for any non-zero and finite ξ, the unstable α-mode could exist for any propagation angle provided that the collision rate is small enough. On the contrary, large angle propagation of the unstable G-mode was prohibited regardless of the value of ν consider herein.
We would also like to mention that in equilibrium, with a specified temperature, other thermodynamic properties of the system such as the number density or energy density are uniquely determined. However, this is no longer true in an anisotropic system when the hard momentum scale λ in the distribution function is specified because both the number density and energy density will change with the anisotropy parameter ξ. Thus, it is also possible to study the collective modes in an anisotropic QCD plasma by holding either the number density or energy density constant. For example, as one changes the anisotropy parameter ξ, the energy density can be held constant which equals to that in the thermal equilibrium as long as we adjust the hard momentum scale as λ → λ(f 0 (ξ)/2) −1/4 . Here, f 0 (ξ) is given by Eq. (25). On the other hand, such an adjustment only leads to a simple modification on the Debye mass defined in Eq. (9), namely, m D → m D (f 0 (ξ)/2) −1/4 . As a result, it is trivial to compute the corresponding structure functions as well as the mass scales because the only change is an extra overall factor (f 0 (ξ)/2) −1/2 in Eqs. (20)- (24) and (31)-(35). With fixed energy density, the analysis of the collective modes doesn't involve anything new technically and the resulting collective modes have no qualitative difference as compared with those present in Secs. IV and V where we hold the hard momentum scale fixed while varying the anisotropy parameter ξ. Some examples of results obtained with fixed energy density are given in App. B.
As is well known, the appearance of the unstable modes is the most distinctive feature of an anisotropic QCD plasma. On the other hand, including a BGK collision kernel further complicates the gluonic collective modes. Our preliminary results indicate that there are some previously undiscovered behaviors which could emerge at large collision rate. For example, for a given anisotropy ξ, the complex solutions for ω as given in Eq. (39) cannot exist when the collision rate is large enough. Instead, pure imaginary solutions appear. It is certainly important to investigate these new modes and their possible physical implications. Another interesting path forward for future work is to assess the effects of the collision kernel on some other physical observables, such as the in-medium properties of the heavy quarkonia. These studies are currently in progress.
With the collision kernel, when ξ is small, up to linear order in ξ, the five structure functions read α(ξ,ω, t,ν) =ω In this appendix, we present some example results obtained when one requires that the energy density is held fixed as ξ is varied. Our results are shown in Figs. 14 and 15. As can be seen from these figures, the qualitative behavior of the stable and unstable modes is the same as in the unnormalized case presented in the main body of the text. FIG. 14. Parameter dependence of the dispersion relations of the stable α-modes where we fix the energy density in the anisotropic QCD plasma. Left: the t-dependence at fixed ξ and ν. Middle: the ν-dependence at fixed ξ and t. Right: the ξ-dependence at fixed ν and t. In these plots, mg = mD/ √ 3.