Collective modes in anisotropic plasmas

We study collective modes in anisotropic plasmas of quarks and gluons using a quasi-particle picture and a hard loop approximation. We use a general class of anisotropic distribution functions, and we consider chirally asymmetric systems. We introduce a complete tensor basis to decompose the gluon polarization tensor into a set of nine scalar functions. We derive and solve the corresponding dispersion equations. Imaginary modes are particularly important because of their potential influence on plasma dynamics. We explore in detail their dependence on the chiral chemical potential and the parameters that characterise the anisotropy of the system. We show that our generalized distributions produce dispersion relations that are much richer in structure than those obtained with a simple one parameter deformation of an isotropic distribution. In addition, the size and domain of the imaginary solutions are enhanced, relative to those obtained with a one parameter deformation. Finally, we show that the influence of even a very small chiral chemical potential is significantly magnified when anisotropy is present.


I. INTRODUCTION
We consider a relativistic plasma of charged particles that is not necessarily in equilibrium, but is in the regime where a quasi-particle description is valid. We use an effective phasespace description in terms of distribution functions that are not assumed thermal or isotropic.
A convenient parametrization for the distribution function of a spheroidially anisotropic system was introduced in [1,2]. The idea is to introduce one parameter that effectively stretches or squeezes the spherically symmetric isotropic distribution in one direction. This parametrization is interesting in the context of heavy ion collisions where one expects the parton distribution that is produced by the collision to be squeezed along the beam axis. Dispersion relations are obtained from the poles of the dressed propagators of the partonic quasi-particles. One finds that the spectrum of collective excitations is much richer than the spectrum of a system in thermal equilibrium. It is especially interesting that there are imaginary modes, which are associated with instabilities and could play an important role in the thermalization of the plasma [3][4][5][6][7][8]. Anisotropic plasmas have been studied intensively in the context of quark-gluon plasmas, which thermalize earlier than expected, for reasons that are as yet unclear. Collective modes have been calculated for both gluons [1,2,9] and quarks [10]. An ellipsoidal generalization was developed in [11,12]. These distributions have been used to study various properties of quark-gluon plasma (QGP) including transport coefficients [13][14][15][16], heavy quark bound states [17][18][19][20], and photon and dilepton production [10,[21][22][23][24].
Isotropic plasmas also support imaginary collective modes, if there is a chiral imbalance [25]. Chiral plasmas with spheroidal anisotropy were studied in Ref. [26]. Chiral systems are of interest in a wide variety of different situations in particle physics, nuclear physics, condensed matter physics, and cosmology. One important example is the chiral magnetic effect (CME) which is the production of a parity violating current in a plasma, in the presence of a magnetic field and an asymmetry between left and right handed fermions [27].
A chirally asymmetric plasma can be described using a chiral chemical potential defined as µ 5 ≡ (µ R − µ L )/2, where µ L and µ R are the chemical potentials of the left and right handed fermions. The CME is particularly interesting in the context of heavy-ion collisions because it has been argued that the quark gluon plasma that is produced in non-central collisions may contain regions where µ 5 is locally finite [28][29][30][31][32]. Instabilities in chiral plasmas have also been studied in electroweak theory at large lepton chemical potential [33] and in the context of the early universe at T µ 5 [34].
Collective modes can be studied using either kinetic theory, or effective field theories. The equivalence of the kinetic theory method and a hard thermal loop (HTL) effective theory approach, for isotropic and chirally symmetric systems, was shown in [35][36][37]. For anisotropic chirally symmetric systems, the connection between kinetic theory and a generalization of the HTL effective theory which is called a hard loop (HL) effective theory has also been established [38,39]. Isotropic chiral kinetic theories have been developed in [40][41][42][43][44], and applied to a variety of physical problems, see for example [45][46][47].
In this paper we use anisotropic distributions that are more general than those used in previous works, and we study the spectrum of collective modes in anisotropic plasmas with nonzero chiral chemical potential. The method can be applied generally to either a QED plasma of ultrarelativistic electrons and positrons or a QGP. We use natural units where = c = 1.
We use capital letters for four-vectors so that, for example, P 2 = p 2 0 − p · p = p 2 0 − p 2 . We define the unit vectorp = p/| p| and we will also usep 0 = p 0 /p. Real solutions to dispersion equations will be denoted ω( p) and imaginary solutions are written iγ( p).

II. ISOTROPIC FORMALISM
In vacuum the photon polarization tensor can be written in terms of one scalar function using a transverse projection operator as Π µν (p 0 , p) = (P 2 g µν − P µ P ν ) Π(p 0 , p). At finite temperature the rest frame of the heat bath, which we define with the four-vector n µ = (1, 0, 0, 0), breaks Lorentz invariance. Time-like axial gauge (TAG) is particularly useful at finite temperature, because the gauge condition is imposed in the heat bath rest frame. In TAG we only have to consider the spatial components of the propagator, for which the Dyson equation has the form In an equilibrated chirally symmetric plasma, the polarization tensor is symmetric and has two independent components, one transverse and the other longitudinal with respect to the three-vector p. An isotropic system with finite chiral chemical potential was first considered in [25], and we review the formalism presented in that paper below. In the presence of finite µ 5 the polarization tensor has an additional asymmetric component and can be written The propagator, obtained from inverting the Dyson equation, is The dispersion relations are obtained from the dispersion equations which give the poles of the retarded propagator. From equation (3) we have that the dispersion equations for an isotropic chiral plasma are Throughout this paper we refer always to the retarded polarization tensor, and we will not include any subscripts to indicate this.
The distributions for right/left handed particle/anti-particles are written n R (k),n R (k), n L (k) andn L (k), where, for example, n R (k) = n(k) µ=µ R .
We calculate the 1-loop photon polarization tensor in the HTL approximation (see Appendix A for details). The contribution from right handed fermions is and the corresponding expression for left handed fermions is obtained from the transformation µ R → µ L . We use P · V = p 0 − p · v and k/ √ k 2 + m 2 ≈k ≡ v since massless fermions are consistent with the hard loop approximation. It is straightforward to modify the HTL integrals to the case of a QCD plasma by changing the QED coupling constant to the QCD one, and using [48] n(k) +n(k) → N f 2 (n q (k) +n q (k)) + N c n g (k) where the subscripts q and g refer to quark (or anti-quark) and gluon distribution functions.
Equations (8,9) can be rewritten in a form that is sometimes more useful by integrating by parts. A straightforward calculation gives From the definitions in equations (A9, A10) one finds The functions Π T (p 0 , p), Π L (p 0 , p) and Π A (p 0 , p) in equations (2,3) are calculated by applying the projection operators P ij T = (δ ij −p ipj )/2, P ij L =p ipj , and P ij A = −i ijmpm /2 to Π ij (p 0 , p). The resulting expressions for the three scalar components of the polarization tensor are where we have used µ 5 = (µ R − µ L )/2 and defined the Debye mass parameter The transverse and longitudinal components of the polarization tensor are the familiar HTL results, and in the chirally symmetric limit µ R = µ L gives µ 5 = 0 and therefore Π A (p 0 , p) = 0. Throughout the rest of this paper we use the shorthand notationμ 5 = g 2 µ 5 /(2π 2 ). In all of our numerical calculations we set the Debye mass to one, which is equivalent to defining all dimensionful quantities in units of the Debye mass.
The dispersion relations are obtained by solving equations (4)(5)(6) where the functions Π L (p 0 , p), Π T (p 0 , p) and Π A (p 0 , p) are given in equation (14). There are pure real and pure imaginary solutions that appear in positive and negative pairs. Equation (4) has two real solutions which are just the usual longitudinal HTL modes, and are written ±ω L (p). Equation (5) also has two real solutions, which we denote ±ω + (p). Equation (6)

A. Distribution function
Next we consider the situation that the momentum distribution is not isotropic. A simple way to account for momentum anisotropy is to replace the distribution function in equations (8,9) using n(k) → n( k) = C ξ n kH ξ ( v) where v = k/k, and similarly forn(k). The subscript ξ indicates dependence on a set of anisotropy parameters that can be used to construct a distribution that is deformed relative to the isotropic one. The factor C ξ is a normalization and can be defined in different ways depending on the calculation being done.
Our choice is explained at the end of this section.
We can construct a completely general expression for the function H ξ ( v) as a sum of terms that are products of anisotropy parameters and dot products of the vector v with two perpendicular unit vectors, which we will writen 3 andn 1 . We are primarily interested in distributions of partons produced in heavy-ion collisions. In this context, we will taken 3 to define the beam axis andn 1 to give the direction of transverse anisotropy. We restrict to functions that satisfy the condition H ξ ( v) = H ξ (− v) and use an expression of the form   When ξ 0 = 1 and ξ i =0 = 0, we have H 2 ξ ( v) = 1, and the distribution is isotropic. For an arbitrary choice of the anisotropy parameters ξ i , the isotropic distribution is expanded in The values of the anisotropy parameters must be chosen so that H 2 ξ ( v) is positive for all orientations of the vector v, which is equivalent to the requirement that H ξ ( v), and therefore the argument of the distribution function, is real and positive. In section III C we will show that our expression (16) is equivalent to including spherical harmonics up to second order.
Equation (16) includes the parametrizations used in several previous calculations. In the original work of Ref. [1] the authors used ξ 0 = 1 and ξ 9 ∈ (−1, ∞) with all other parameters set to zero. The chosen value of ξ 9 allows one to consider either a slightly prolate momentum distribution (−1 < ξ 9 < 0), which is elongated along the beam axis, or an arbitrarily oblate distribution (ξ 9 > 0), which is squeezed in the direction of the beam axis. In the calculations of Ref. [11,12], where ellipsoidal asymmetry was included for the first time, the authors used ξ 0 = 1 and various values of ξ 2 and ξ 9 with all other parameters zero. In Ref. [9] the authors considered a one parameter deformation using ξ 0 = 1 + σ and ξ 9 = −σ with σ ∈ (−1, ∞). The resulting distribution is slightly oblate for −1 < σ < 0 and prolate to any degree for σ > 0.
Our definition (16) can be extended by including terms with higher order products of the scalar products (n 3 · v) and (n 3 · v). The only restriction we impose, besides positivity, is the condition . This means there are no terms where the sum of the exponents is odd, like for example (n 3 · v)(n 1 · v) 2 . The reason we require that H ξ ( v) is even under the transformation v → − v is as follows. We said above that when general distributions are considered the polarization tensor is obtained from equations (8,9) by making the replacement n(k) → n( k) = C ξ n kH ξ ( v) . In fact, there is an extra contribution to the anisotropic version of equation (8) that integrates to zero for any distribution that is even under the transformation (17) would produce a nonzero contribution to the polarization tensor that dominates over the HL terms, and is not present in the result obtained from semiclassical kinetic theory (which is equivalent to the HL expression).
We will calculate the anisotropic polarization tensor using equations (11,12) with the replacement n(k) → C ξ n (kH ξ ( v)) and similarly forn(k). We definek = kH ξ ( v) and perform a straightforward change of variable so that the integral over k is written as the product of an integral over the new variablek and a factored integral involving the angular variables. Equations (11,12) become where we have defined and used (m 2 D ) R = g 2 (T 2 /3 + µ 2 R /π 2 ). The full expression for Π ij is obtained using equation (13), as in the isotropic case.
We close this section by explaining how we define the function of the anisotropy parameters which we call C ξ , which can be thought of as a normalization of the anisotropic distribution function. This factor is determined by enforcing the condition from which we obtain This choice of normalization gives that when the chemical potentials are zero the mass parameter in equation (15) is independent of the anisotropy parameters. At zero chemical potential the entire spectrum of collective modes depends on this one scale (since we have assumed massless partons -see under equation (7)), and it is therefore natural to adopt a normalization that leaves the Debye mass invariant under a change of the anisotropy parameters. We also note that using equations (18,20,22) it is easy to see that the even part of the polarization tensor is unchanged when the complete set of anisotropy parameter are scaled uniformly: ξ i → Λξ i where Λ is any positive constant. This means that when we compare the spectrum of collective modes obtained from different choices of anisotropy parameters, we are comparing the effects of a specific deformation, and not the influence of an overall change of scale.
The anisotropic dispersion relations are the solutions of the dispersion equation, D = 0.

C. Coordinate system
To make further process we must define a coordinate sytem. The vector along the beam axis can be chosen asn 3 = (0, 0, 1) and the vector that defines the direction of asymmetry in the transverse plane can then be chosen asn 1 = (1, 0, 0). The external momentum vector then has the general form p = p(sin θ cos φ, sin θ sin φ, cos θ). We note that using these coordinates it is easy to see that the parametrization in equation (16) is equivalent to an expansion in spherical harmonics Y lm (θ, φ) with l ∈ (0, 4) and m ∈ (−l, l).
We can define coordinates that are more convenient for calculational purposes by performing a rotation [12]. We rotate counter-clockwise about the z-axis with angle θ z = π/2 − φ and then counter-clockwise about the x-axis with θ x = θ. The result is that the basis vectors becomen 1 = (sin(φ), cos(θ) cos(φ), sin(θ) cos(φ)) n 3 = (0, − sin(θ), cos(θ)) Any scalar or pseudo-scalar quantity will give the same result when calculated with either the original or the rotated coordinate system. This can be easily be seen from the fact that the scalar products (p ·n 1 ), (p ·n 3 ) and (n 1 ·n 3 ), and the pseudo-scalar ijkni Another useful thing about the coordinates in equation (29) is that the isotropic limit has a simple form. It is easy to show using equation (29) that δ ij −p ipj = P ij 1 + P ij 3 ,p ipj = P ij 2 and ijm p m = P ij 8 . In the isotropic limit there are only three different non-zero components in equation (26) which are π 1 = π 3 = Π T , π 2 = Π L and π 8 = Π A . Using these results one can see that equations (27,28) reduce to (3).
We will use the notation x = cos(θ) and x = cos(θ ).

D. Anisotropy parameters
The anisotropy parameters ξ 9 and ξ 2 produce uniform stretching or squeezing of the isotropic distribution along the beam axis (in the case of ξ 9 ), and along the vectorn 1 in the transverse plane (for ξ 2 ). We have introduced additional anisotropy parameters that can be used to produce distribution functions that are deformed relative to the isotropic distribution in more general ways, and produce dispersion relations with more structure. As discussed in section III A, we consider only distributions that satisfy the condition f ( k) = f (− k). This restricts the form of the terms in equation (16). Each term is chosen so that e 1 + e 3 is even, where the exponent of the factor ( n 1 · v) is denoted e 1 and the exponent of the factor ( n 3 · v) is called e 3 . There are two different possibilities: e 1 and e 3 individually even, or e 1 and e 3 individually odd.
The elliptically anisotropic distribution considered in [12] depends only on the two param- To illustrate more directly the roles of the different anisotropy parameters we look at plots of C ξ n kH ξ ( v) with β = 1 and µ = 0, in the (p x , p y ) and (p x , p z ) planes. In figure 3 we show the results for several choices of anisotropy parameters. The first panel shows, for reference, the isotropic distribution in the (p x , p y ) plane, the figure for the (p x , p z ) plane is identical. In the second panel we show a prolate distribution in the (p x , p z ) plane, and one sees clearly that the distribution is stretched along the beam axis. The third panel shows the squeezing of an oblate distribution along the beam axis. In the fourth and fifth panels we show an elliptically oblate distribution, which is asymmetric in both planes. In the sixth panel we show an example of the kind of structure that can be obtained with a specific choice of the extra anisotropy parameters we have introduced.

IV. DISPERSIONS RELATIONS
A. Analytic structure The isotropic and anisotropic dressing functions have the same analytic structure, which can be determined from equations (11,12), using that equations (29,30) give P · V = p(p 0 − x ) wherep 0 = p 0 /p. Whenp 0 is pure real the epsilon prescription is needed to define the integrands atp 0 = x . Shifting the integration variable v → − v one shows that π * i (p 0 ) = π i (−p 0 ), where the subscript i indicates any component of the polarization tensor. This means that whenp 0 is real, the real part of each dressing function is even inp 0 and the imaginary part is odd. We also note that all dressing functions are pure real forp 0 > 1, since the imaginary part comes from the discontinuity atp 0 = x and |x | < 1. Whenp 0 is not pure real the epsilon prescription is not needed, and we have π * i (p 0 ) = π i (p * 0 ). When p 0 is pure imaginary π i (p * 0 ) = π i (p 0 ), which can be shown by shifting the integration variable, and therefore the dressing functions are pure real and even inp 0 . The information above can be summarized as Re(p 0 ) > 1 and Im(p 0 ) = 0 → π i (p 0 ) real and even p 0 imaginary → π i (p 0 ) real and even Re(p 0 ) < 1 and Im(p 0 ) = 0 →    Re π i (p 0 ) even Im π i (p 0 ) odd .

B. Method
The integrand for the polarization tensor is given by equations (18,19,20). The nine scalar components π 1 -π 9 are obtained by contracting with the appropriate projection operators, which are easily obtained from equations (25,26). These calculations are straightforward but tedious, and they are therefore done with Mathematica. The resulting integrals are calculated numerically using Gauss-Legendre quadrature. For real valuedp 0 , the real part of the dressing functions has a single integrable singularity at x =p 0 that is easily handled with a numerical principal part prescription, and the imaginary part can be obtained from the residue of the pole. For imaginary valuedp 0 there are no difficulties. After the dressing functions are obtained, the dispersion equation is solved with a standard binary search algorithm.

D. Imaginary solutions
We are primarily interested in imaginary solutions, because of their potential to influence plasma dynamics. We will compare the imaginary modes produced by different choices of the anisotropy parameters, and explore their dependence on these parameters as a function of wave vector. We have previously used the notation γ( p) to indicate the dependence of the dispersion relation on the wave vector, but in this section we will use the more explicit notation γ(p, x, φ) where x = cos(θ). For different sets of anisotropy parameters we compare the maximum value of the imaginary solution, the position of this maximum, and the integral of the solution over its domain, which we write γ ≡ ∞ 0 dp 1 −1 dx 2π 0 dφ γ(p, x, φ).
We define three variables that characterise the degree of oblateness, the transverse eccen-tricity, and the azimuthal asymmetry of a given distribution: The angle brackets denote averaging over the momentum phase space, for example k 2 x = d 3 kk 2 x C ξ n(kH ξ ( v)).
In table I  For a given distribution, the number of imaginary solutions is 0 or 1 or 2, depending on the magnitude and direction of the wave vector. For each mode there is a critical vector p crit (x, φ, ξ i ) such that the imaginary mode appears for p less than this critical value. Analytic results for the critical wave vectors can be obtained in the limit of weak anisotropy using a Nyquist analysis. The method is explained in detail in Appendix B, where we show how to identify the angular variables that produce the largest critical wave vectors for a given set of anisotropy parameters. To understand the structure of the solutions obtained from a given distribution in more detail, we show several figures below.
In figure 5 we compare the imaginary solutions obtained from the distributions 3 and 6 for a specific choice of angles. The figure illustrates the concept of the critical wave vector.
For example, in the right panel of figure 5, the critical wave vectors for modes 1 and 2 are, respectively, 2.33 and 0.62. It is true in general that when the domain of the solution increases, the maximum value of the solution also increases, as seen in figure 5.   The distributions 6-9 involve anisotropy parameters that correspond to higher spherical harmonics, and have a much richer structure. In figures 8 and 9 we show some contour plots of the largest imaginary solution obtained from these distributions withμ 5 = 0, for different choices of x = cos(θ) and φ. give s 2 = 0, which corresponds physically to symmetry across the reaction plane. If any of these anisotropy parameters are non-zero, the value of s 2 will also be non-zero. Physically this effect could be produced by a collision at non-zero impact parameter of two ions of different size. Distribution 6 and 7, for which we should find s 2 = 0, give s 2 ≈ 9 × 10 −10 and s 2 ≈ 7 × 10 −8 , which can be taken as a measure of the numerical error in our computation.

V. CONCLUSIONS
We have studied plasmons in anisotropic plasmas using a hard loop approach within the quasi-particle regime. We have used more general distribution functions than have been previously considered, with a technique that could be easily generalized to account for the anisotropies that are relevant in any given physical situation. The distribution function we have used depends on 9 anisotropy parameters and is normalized so that at zero chemical potential the Debye mass parameter is the same for all distributions. In an ultra-relativistic plasma at zero chemical potential this parameter is the only dimensionful scale. It is therefore meaningful to compare the solutions obtained from the dispersion equations for different choices of anisotropy parameters. Imaginary solutions to the dispersion equation are important because they are associated with plasma instabilities, which have the potential to strongly influence plasma dynamics.
The important effect of imaginary collective modes on physical quantities in anisotropic plasmas, like transport coefficients, production rates, and bound state properties, is well documented in the literature. We have shown that the more general distribution functions we have introduced can significantly increase imaginary solutions of the dispersion equation in terms of both the magnitude of the mode, and the domain of wave vectors over which it exists. The distribution that is most commonly used in the literature depends on only one anisotropy parameter, which describes the oblateness of a distribution that is squeezed in one direction. An example is our distribution 3. In comparison, our distribution 7, which has approximately the same value of the parameter ξ that characterises its oblateness (see table I), produces a largest imaginary mode that is almost twice as big, and the integral of the solution over the momentum phase space is approximately 10 times larger. We have also shown that the effect of the chiral chemical potential is greatly enhanced by anisotropy. This is seen clearly from the results in table I which show that, relative to the isotropic result, the imaginary modes are much greater when even very moderate anisotropy is introduced. In this Appendix we calculate the one-loop retarded polarization tensor at finite temperature and chemical potential with the real-time formulation of finite temperature field theory, using the Keldysh representation [50][51][52].
The electron propagator is a 2×2 matrix of the form where the retarded, advanced and symmetric propagators are given by The self-energy has the form and the vertex function is a 2×2×2 tensor which can be written The contribution to the retarded self-energy from the one-loop diagram is The sum over Keldysh indices {i, i , j, j } ∈ {r, a} is easily done because G aa = 0 and a vertex function with an odd number of a indices vanishes. The result is where we have used the short hand notation r(K) to represent r(k 0 , k), and similarily for other propagator components. From this point on we will suppress the subscript that indicates the retarded component of the polarization tensor. In addition we consider only the relativistic limit T m and therefore we drop the electron mass in the integrand in equation (A8). We also introduce the short-hand notation dK ≡ d 4 k/(2π) 4 .
We will consider a chirally asymmetric plasma which is characterised by a difference in the chemical potentials of the right and left handed fermions, which are denoted µ R and µ L . We The integrals in equations (A9) can be rewritten by performing a shift of the 4-vector K → K + P on the term in square brackets that does not have a factor f (K). All terms in the resulting expression have a factor f (K) and, using equation (A4), this factor is proportional to δ(K 2 ) in the relativistic limit. The integral over k 0 can be done using this delta function. We work in time-like axial gauge and calculate only the spatial components of the polarization tensor. After doing the k 0 integral we obtain the results in equations (8,9). and define the function We consider the contour integral where the contour is a positively (counter-clockwise) oriented closed loop, which is chosen so that F (p 0 ) is analytic inside the loop except at isolated points. The integral is equal to the sum of the residues. The residue of F (p 0 ) at a zero of f (p 0 ) of order l is l, and the residue of F (p 0 ) at a pole of f (p 0 ) of order l is −l. This gives where n Z and n P are the numbers of zeros and poles of f (p 0 ) inside the contour C, taking into account the fact that each zero and pole of order l is counted l times. The purpose of the Nyquist analysis is to determine n Z .
The anisotropic dressing functions are obtained numerically from the integrals in equations (18,19). Using P · V = p(p 0 − x ) (see equations (29,30)) one finds that the dressing functions have a cut along the real axis for |p 0 | < 1, and the dispersion equation therefore has a cut for |p 0 | < p. This means that the contour C can be chosen as in Fig. 11. The integrals along the lines connecting the circular contour C ∞ to C cut always compensate each other and therefore the contour integral (B4) equals The contribution from the big circle is calculated by writing p 0 = |p 0 |e iφ and taking |p 0 | → ∞ which gives The integral along the cut can be calculated using that F (p 0 ) is the logarithmic derivative of f (p 0 ) which gives where p 0start is the (arbitrarily chosen) starting point of the contour which encloses the cut, and p 0end is the end point. The quantity n W is the number of times that the curve in the plane of complex f travels counter-clockwise around the point f = 0, and is called a winding number. Combining the results (B6, B7), we have that equation (B5) gives mode, so that n (2) ∞ = 0, and for all of the other dispersion equations it is straightforward to see that n ∞ = 2. tions When the chirality parameter µ 5 is non-zero we consider equations (B12, B13) instead of equations (B9, B11). For weak anisotropy the winding number is determined from the sign of the functions f (4) (0, p) = − (p + p 5crit (ξ i , x, φ)) (B18) The results for p 2 1crit , p 2 3crit and p 5hat are given in equations (B20, B21, B22).
In the isotropic case p 2 1crit (ξ i , x, φ) = p 2 3crit (ξ i , x, φ) = 0, the right side of equations (B16, B17) cannot be positive, and therefore the winding number is zero and both dispersion equations have two solutions. Remembering that in the isotropic limit π 1 = π 3 = Π T , we see that these are the HTL transverse solutions (±ω T (p)). When non-zero anisotropy parameters are introduced, and when p is small enough, there can be choices of the angles and anisotropy parameters for which p 2 1crit (ξ i , x, φ) > p 2 > 0 and/or p 2 3crit (ξ i , x, φ) > p 2 > 0. In this case, the winding number is two, which means that extra solutions to equations (B9, B11) exist.
From equation (B22) we have that in the isotropic limit p 5crit (ξ i , x, φ) =μ 5 , which means that n (4) W = 0 and n (5) W = 2 if p <μ 5 and zero otherwise. Thus we have that for µ 5 = 0 an additional pair of solutions, which turn out to be imaginary, exist even in the isotropic limit (see figure 1). The domain of these modes can be larger for certain orientations of the wave vector if the system is anisotropic.
The expressions in equations (B20, B21, B22) are only valid in the weakly anisotropic limit, however, we can use them to obtain information about the orientations of the wave vector p for which large imaginary solutions could be found with specific sets of anisotropy parameters. To obtain more insight, we can plot the coefficient of each anisotropy parameter in equations (B20 -B22). We show the results for three of the coefficients in figure 13. The first panel shows the dependence of the critical wave vectors on ξ 9 , which is the coefficient of (n 3 · v) 2 in equation (16). Using our coordinate system this dot product is independent of the azimuthal angle φ, and therefore we have that the terms in the critical wave vectors that are proportional to ξ 9 are also independent of φ. The second panel shows the dependence of the critical wave vectors on ξ 2 which is the coefficient of (n 3 · v) 2 . Both the ξ 9 and ξ 2 terms have been included in previous calculations. In the third panel we show the dependence of the critical wave vectors on ξ 6 , which is the coefficient of (n 1 · v)(n 3 · v), and in the last panel we show the dependence on ξ 8 , the coefficient of (n 1 · v) 3 (n 3 · v). (left), p 2 3crit (center) and p 5crit (right). The axes show x = cos(θ) ∈ (0, 1) and φ ∈ (0, π/2).