Angular Momentum Mixing in a Non-spherical Color Superconductor
aa r X i v : . [ h e p - ph ] N ov RU06-9-B
Angular Momentum Mixing in a Non-spherical Color Superconductor
Bo Feng, ∗ De-fu Hou, † and Hai-cang Ren
2, 1, ‡ Institute of Particle Physics, Huazhong Normal University, Wuhan, 430079, China Physics Department, The Rockefeller University, 1230 York Avenue, New York, NY 10021-6399 (Dated: October 30, 2018)We study the angular momentum mixing effects in the color superconductor with non-sphericalpairing. We first clarify the concept of the angular momentum mixing with a toy model for non-relativistic and spinless fermions. Then we derive the gap equation for the polar phase of denseQCD by minimizing the CJT free energy. The solution of the gap equation consists of all angularmomentum partial waves of odd parity. The corresponding free energy is found to be lower thanthat reported in the literature with p -wave only. PACS numbers: 26.30.+k, 91.65.Dt, 98.80.Ft
I. INTRODUCTION
The properties of quark matter at extreme conditions have been an active research area both theoretically andexperimentally. At high temperature, the quark-gluon plasma(QGP) has long been searched by colliding two nucleiat sufficiently high energy. On the other side, we expect that quark matter becomes color superconducting througha phase transition at high baryon density but low temperature[1, 2, 3, 4, 5, 6], which is the typical condition insidecompact stars.In a typical metallic superconductor, the electrons pair with equal chemical potential near the Fermi surface.The situation with a quark matter, however, is much more complicated. While the quark-quark interaction favorspairing between quarks of different flavors, the mass difference among u , d and s together with the charge neutralityrequirement induces a substantial mismatch among their Fermi momenta at the baryon density inside a compact starand thereby reduces the available phase space for Cooper pairing. A number of exotic color superconductivity phasesin the presence of mismatch have been proposed in the literature[7, 8, 9, 10, 11], but a consensus point of view of thetrue ground state has not been reached. The single flavor pairing[12, 13, 14, 15, 16], which is free from the Fermimomentum mismatch, is an interesting alternative in this circumstance and will be considered here. Since the quark-quark interaction is attractive in the color anti-triplet channel, the color wave function of the pair is anti-symmetric.For the equal helicity pairing to be considered in this article, the parity of the orbital wave function has to be oddas required by the Pauli principle. Except for the color-spin-lock phase examined in [12], the energy gap will not bespherical. The odd parity prevents the diquark wave function from realizing the full pairing potential. The energyscale of the color superconductivity is therefore reduced.At ultra-high baryon densities, asymptotic freedom of QCD ensures the validity of the weak coupling expansion,which has been carried out for CSC by a number of authors[6, 17, 18, 19, 20, 21, 22, 24, 25, 26, 27]. The dominantpairing interaction is mediated by one-gluon exchange and can be decomposed into partial waves, as is shown inEq.(59) below. Quantitative results of the transition temperature in the equal helicity channel with an arbitraryangular momentum has been obtained from the first principle [21, 22], and read T ( J ) c = 512 π ( 2 N f ) / µg exp h − π √ g + γ −
18 ( π + 4) + 3 c J i (1)where N f is the number of flavor, N is the number of colors, µ is the chemical potential, g is the running couplingconstant of QCD and γ (= 0 . ... ) is the Euler constant. The J -dependent constant c J = , for J = 0 , − J P n =1 1 n for J > . (2) ∗ Electronic address: fengbo@iopp. ccnu. edu. cn † Electronic address: hdf@iopp. ccnu. edu. cn ‡ Electronic address: [email protected]
The equal helicity pairing of odd parity picks up T (1) c as the transition temperature. We have T (1) c = e − T (0) c ≃ . × − T (0) c . Natural analogy is drawn with the superfluidity of He . But important difference between thepairing potential in quark matter and that in He has to be considered before ascertaining the angular dependence ofthe energy gap. The forward singularity of the one-gluon exchange renders the pairing strength equal for all partialwaves (the same leading order term inside the bracket of (1)). The paring potential in He , however, is entirely in thechannel of J = 1.The transition temperature was determined from the pairing instability of the diquark scattering amplitude in thenormal phase. In a perturbative treatment, the Dyson-Schwinger equation for the scattering amplitude is linear andthe partial wave decomposition in Eq.(1) is legitimate. This is not the case with the gap equation below T c , whichis nonlinear. A non-spherical solution of the gap equation, in general, will be a mixture of different partial wavesunless the the pairing is dominated in one angular momentum channel. The adjective ”non-spherical” refers to themagnitude of the energy gap of the quasi-particle. So the s -wave and the CSL gaps are spherical and are thereforefree from the angular momentum mixing. The gap equations in these cases are linear with respect to the angulardependence. A nonspherical gap function, φ M , for the equal helicity pairing at T < T c contains all spherical harmonicsof odd parity with the same azimuthal quantum number M . We have φ M = φ f M (ˆ p ) (3)where φ is the gap function of 2SC in the absence of the mismatch, ˆ p the direction of the relative momentum ofthe two quarks in a Cooper pair and the angular dependent factor f M (ˆ p ) = X J =1 , , ,... b J Y JM (ˆ p ) (4)with J the total angular momentum of the Cooper pair. Carrying the formulation of He over to QCD amounts todrop all higher multipoles except that of J = 1, which will not satisfy the gap equation of QCD. It was argued inthe literature that b = O (1) but b J = O ( g ) for J = 1, This, as will be shown below, is not the case. Instead, wefind that the function f M (ˆ p ) satisfies a nontrivial integral equation and thus b J = O (1) for all odd J ’s. Therefore theangular momentum mixing does occurs in the subleading order of the gap function. The angular momentum mixingwill modify all non-spherical ”spin-1” CSC examined in the literature, we shall focus our attention in this paper to theequal helicity pairing with zero azimuthal quantum number, i.e. the analog of the polar phase of He . The subscript M of φ M and f M will be suppressed below. Even though this phase is unstable, it is the simplest one to illustratethe mixing mechanism.The current work is organized as follows. In the next section, we shall clarify the concept of the angular momentummixing with a toy model of non-relativistic and spinless fermions. In the Sect. III the gap equation for the single flavorCSC will be derived by minimizing the CJT free energy of QCD. This gap equation will be reduced to an nonlinearintegral equation for the function f (ˆ p ) in the subsequent section and the numerical solution will be presented in theSect, V. We conclude the paper in the Sect. VI. Some technical details are deferred to the Appendices. Our unitsare ¯ h = c = k B = 1 and 4-vectors are denoted by capital letters, K ≡ K µ = ( k , ~k ) with k the Matsubara energy,which becomes continuous at T c = 0. Throughout the article, we shall follow the definition of the leading order andthe subleading order in [6]. Upon taking the logarithm of the transition temperature or the magnitude of the gapfunction, the order O ( g ) will be referred to as the leading one and the O (1) term to the subleading one. II. A TOY MODEL WITH ANGULAR MOMENTUM MIXING
To clarify the concept of the angular momentum mixing, we consider a toy model of nonrelativistic and spinlessfermions. In terms of the creation and annihilation operators, the model Hamiltonian reads H = X ~p ǫ p a † ~p a ~p − λ X ~p,~p ′ ′ V (ˆ p · ˆ p ′ ) a † ~p a †− ~p a − ~p ′ a ~p ′ (5)where ǫ p = p m − µ with m the mass and µ the chemical potential, λ > P ′ ~p,~p ′ extends to states with | ǫ p | < ω D and | ǫ p ′ | < ω D with ω D a UV cutoff (Debyefrequency for electronic superconductors). The angular dependent form factor V (ˆ p · ˆ p ′ ) can be expanded in series ofthe Legendre polynomials, V (ˆ p · ˆ p ′ ) = ∞ X J =0 (2 J + 1) v J P J (ˆ p · ˆ p ′ ) . (6)Introducing the order parameter χ ( ~p ) = < | a − ~p a ~p | > (7)with | > the ground state and expanding the interaction term of (5) to the linear order of the fluctuation a − ~p a ~p − < | a − ~p a ~p | > , we end up with the mean-field Hamiltonian H MF = 12 X ~p, | ǫ p | <ω D χ ∗ ( ~p ) φ ( ~p ) + X ~p ǫ p a † ~p a ~p − X ~p, | ǫ p | <ω D h φ ∗ ( ~p ) a − ~p a ~p + φ ( ~p ) a † ~p a †− ~p i , (8)where we have introduced the gap function via φ ( ~p ) = λ X ~p ′ , | ǫ p ′ | <ω D V (ˆ p · ~p ′ ) χ ( ~p ′ ) . (9)We have χ ( − ~p ) = − χ ( ~p ) and φ ( − ~p ) = − φ ( ~p ) following from their definitions. Upon a Bogoliubov transformation, wefind that χ ( ~p ) = φ ( ~p )2 ε p (10)with ε p = q ǫ p + | φ ( ~p ) | and the ground state energy E = X ~p,ǫ p < ǫ p + Ω F (11)with F the condensation energy density given by F = 12 X ~p, | ǫ p | <ω D h φ ∗ ( ~p ) φ ( ~p )2 ε p + | ǫ p | − ε p i . (12)Substituting (10) into (9), we obtain the gap equation φ ( ~p ) = λ X ~p ′ , | ǫ p ′ | <ω D V (ˆ p · ˆ p ′ ) φ ( ~p ′ ) ε p ′ . (13)In the weak coupling limit, ω D << µ and λD F << D F = m µ √ π the density of states on the Fermi surface,but the magnitude of V (ˆ p · ˆ p ′ ) remains of order one. We have X ~p, | ǫ p | <ω D = Ω Z | ǫ p | <ω D d ~p (2 π ) ≃ Ω D F π Z d ˆ p Z ω D − ω D dǫ. (14)Also, the support of the gap function extends only to a narrow band around the Fermi surface. We may ignore thedependence of φ on the magnitude p = | ~p | and switch the argument of φ from ~p to ˆ p . Following (14), the integrationsover p in (12) and (13) can be carried out readily and we end up with F = − D F π Z d ˆ p | φ (ˆ p ) | . (15)and φ (ˆ p ) = λD F π Z d ˆ p ′ V (ˆ p · ˆ p ′ ) φ (ˆ p ′ ) ln 2 ω D | φ (ˆ p ′ ) | . (16)The gap equation (16) is nonlinear because of the logarithm.In what follows, we consider two extremes of V (ˆ p · ˆ p ′ ), each of which gives rise to an exact solution to the gapequation (16). We present only the solution that is invariant under time reversal, i.e. the one with zero azimuthalquantum number. Case I : V (ˆ p · ˆ p ′ ) = 3 P (ˆ p · ˆ p ′ ) = 3ˆ p · ˆ p ′ . (17)It corresponds to the partial wave expansion (6) with v J = 1 and v J = 0 for J = 1. The angular dependence of thepairing force in He is of this type. The gap equation reads φ (ˆ p ) = 3 λD F π Z d ˆ p ′ ˆ p · ˆ p ′ φ (ˆ p ′ ) ln 2 ω D | φ (ˆ p ′ ) | (18)and its solution of zero azimuthal quantum number is given by φ (ˆ p ) = φ cos θ = φ P (cos θ ) (19)with θ the angle with respect to a prefixed direction in space and φ = 2 ω D e − λDF + . (20)The condensation energy density F = − ω D D F e − λDF + ≃ − . ω D D F e − λDF . (21)This solution corresponds to the polar phase of He [23]. Since the gap function contains only the partial wave of J = 1, there is no angular momentum mixing. The additional term in the exponent of (20), , comes from thelogarithm of (18). Case II : V (ˆ p · ˆ p ′ ) = 4 πδ (ˆ p − ˆ p ′ ) = X J (2 J + 1) P J (ˆ p · ˆ p ′ ) . (22)This corresponds to a singularity of the two body scattering amplitudes in the forward direction. We have v J = 1 forall J in (6). The last step of (22) follows from the addition theorem and the completeness of the spherical harmonics.The gap equation (16) becomes 1 = λD F ω D | φ (ˆ p ) | , (23)which implies a constant | φ (ˆ p ) | and yields a solution of odd parity and zero azimuthal quantum number. φ (ˆ p ) = 2 ω D e − λDF sign(cos θ ) = 2 ω D e − λDF ∞ X n =0 ( − n (4 n + 3) (2 n − n +1 ( n + 1)! P n +1 (cos θ ) . (24)The condensation energy density in this case reads F = − ω D D F e − λDF . (25)We refer to this case as the case with the angular momentum mixing because the gap function (24) contains all partialwaves. Carrying the solution of the case I to the case II amounts to drop all partial waves other than that of J = 1and would lead to a lower magnitude of the condensation energy (21).The case with QCD is similar to the case II above since the forward singularity of the diquark scattering renders thepairing strength of all partial waves equal to the leading order. The running coupling constant g of QCD correspondsto λ here and the angular momentum mixing shows up in the O (1) term of ln | φ | . Therefore we expect angularmomentum mixing to the subleading order of the angular dependence of the gap function. Besides being an ultrarelativistic system, the CSC of QCD differs from the toy model considered above in two aspects. The forwardsingularity of QCD also brings about the energy dependence of the gap, so the gap equation (13) will be replaced bythe Eliashberg equation derived by minimizing the CJT effective action of QCD. Secondly, the pairing strength ofeach partial wave does fall off with an increasing J in the sub-leading order of the pairing potential. It is this fallingoff that makes the amount of the angular momentum mixing numerically small for the solution considered in thisarticle. FIG. 1: The sunset diagram of Eq.(29). Straight line and wavy line denotes quark and gluon propagators respectively.
III. DERIVATION OF THE GAP EQUATION FROM THE CJT FREE ENERGY
The QCD Lagrangian for one flavor of massless quark is given by L = ¯ ψ ( iγ µ D µ + µγ ) ψ − G µνa G aµν + renormalization counterterms (26)where, ψ is the quark spinor in Dirac and color space and ¯ ψ = ψ † γ . The covariant derivative acting on the fermionfield is D µ = ∂ µ + igT a A aµ , where g is the running coupling constant, A aµ is the gauge potential, T a = λ a / a = 1 , ..., a -th SU (3) c generator with λ a the a -th Gell-Mann matrix. G aµν = ∂ µ A aν − ∂ ν A aµ + gf abc A bµ A cν is the fieldstrength tensor. Introducing the Nambu-Gorkov spinorΨ = (cid:18) ψψ C (cid:19) , ¯Ψ = ( ¯ ψ, ¯ ψ C ) (27)where ψ C = C ¯ ψ T is the charge-conjugate spinor and C ≡ iγ γ , the CJT effective action reads[6, 30]Γ[ D, S ] = − { Trln D − + Tr( D − D − − Trln S − − Tr( S − S − − [ D, S ] } (28)where D and S are the full gluon and quark propagators, D − and S − are the inverse tree-level propagators forgluons and quarks, respectively. Γ is the sum of all two-particle irreducible(2PI)vacuum diagrams built with D , S and the tree-level quark-gluon vertex ˆΓ. We haveΓ = −
14 Tr( D ˆΓ S ˆΓ S ) + ..., (29)where the first term corresponds to the sunset diagram of Fig.1 and the contribution from ... is beyond the subleadingorder of the gap function[6] .The stationary points of the CJT effective action are determined by δ Γ δD = 0 , δ Γ δS = 0 (30)which gives rise to the Dyson-Schwinger equation for gluons and quarks,Π ≡ − δ Γ δD , Σ ≡ δ Γ δS (31)where Π and Σ are the gluon and quark self-energy defined via D − = D − + Π and S − = S − + Σ. Instead ofsolving the two equations of (31) simultaneously, we shall reduce the CJT effective free energy with the aid of thefirst equation, leaving the gap function arbitrary. The gap equation ( which is the Nambu-Gorkov off diagonal partof the second equation of (31) ) will be derived after the CJT free energy is fully simplified under the weak couplingapproximation.Substituting the first equation of (31) into (28), the second term in Eq. (28) cancels the last term. We haveΓ[ D, S ] = −
12 [Trln D − − Trln S − − Tr( S − S − S − ≡ (cid:18) [ G +0 ] −
00 [ G − ] − (cid:19) (33)where [ G ± ] − = ( p ± µ ) γ − ~γ · ~p. (34)On writing the quark self-energy Σ ≡ (cid:18) Σ + Φ − Φ + Σ − (cid:19) (35)the full quark propagator, S = (cid:18) G + Ξ − Ξ + G − (cid:19) (36)can be obtained explicitly by inverting the matrix S − + Σ.For the single flavor pairing, the simplest choice of the off-diagonal block of Eq. (35) readsΦ + ( P ) = iφγ λ (37)and Φ + = Φ − (see Theorem 2 in[26]), where λ is the 2nd Gell-Mann matrix and φ is a function of the energy andthe momentum, i.e. φ = φ ( p , ~p ). φ is even in p and odd in ~p . By using the energy projectors of massless fermionsΛ ± p = (1 ± γ ~γ · ˆ p ) / G ± = p + ( p ∓ µ ) p − ( p ∓ µ ) − φ λ Λ + p γ + p ( p ± µ ) p − ( p ± µ ) − φ λ Λ − p γ (38)Ξ ± = iφλ γ p − ( p ± µ ) − φ λ Λ + p + iφλ γ p − ( p ∓ µ ) − φ λ Λ − p . (39)Because of the λ of (37), the excitation in the third color direction is ungapped.Now, we proceed to simplify the CJT free energy under the weak coupling approximation. Denote by Γ n the freeenergy density of the normal phase, we have Γ = Γ n + Ω F (40)where the condensate energy density F = −
12Ω [Trln D − − Trln D − n − Trln S − + Trln S − − Tr( S − S − D − − Trln[ D − n ] ≃ Tr[ D n δ Π] (42)where δ Π = Π − Π n , (43)with Π n the hard-dense-loop (HDL) resummed gluon self-energy in normal phase and D n the corresponding HDLgluon propagator. In the Coulomb gauge, the HDL gluon propagator is D n, ( K ) = D l ( K ) , D n, i ( K ) = D n,i = 0 , D n,ij = ( δ ij − ˆ k i ˆ k j ) D t ( K ) (44)where D l,t are the longitudinal and transverse propagators respectively and are diagonal in adjoint color space,i.e. D abl,t = δ ab D l,t . Consequently, we only need the 00-component, Π ( K ), and the transverse projection of theij-components, ( δ ij − ˆ k i ˆ k j )Π ij ( K ) = Π ii ( K ) − ˆ k i ˆ k j Π ij ( K ) (45)The gluon self-energy in super phase readsΠ µνab ( K ) = 12 T X P,P ′ Tr[ˆΓ µa S ( P )ˆΓ µb S ( P ′ )] (46)where K = P − P ′ and ˆΓ aµ ≡ (cid:18) Γ aµ
00 ¯Γ aµ (cid:19) (47)with Γ µa = γ µ T a and ¯Γ µa = − γ µ T Ta . Substituting Eq. (36) into Eq.(46), we find that Nambu-Gorkov space,Π µνab ( K ) = 12 T X P,P ′ (cid:8) Tr[Γ µ a G + (P)Γ ν b G + (P ′ )] + Tr[¯Γ µ a G − (P)¯Γ ν b G − (P ′ )]+Tr[Γ µa Ξ − ( P )¯Γ νb Ξ + ( P ′ )] + Tr[¯Γ µa Ξ + ( P )Γ νb Ξ − ( P ′ )] (cid:9) (48)Since the HDL gluon propagators are diagonal in color space, we only need the diagonal terms of Eq.(48) to deal withEq.(42). The explicit form of each diagonal term of (48) readsTr[Γ µa G + ( P )Γ νa G + ( P ′ )] = g T X P,P ′ Tr (cid:2) γ µ Λ + p γ γ ν Λ + p ′ γ (cid:3) w + a ( P, P ′ ) , (49a)Tr[¯Γ µa G − ( P )¯Γ νa G − ( P ′ )] = g T X P,P ′ Tr (cid:2) γ µ Λ − p γ γ ν Λ − p ′ γ (cid:3) w − a ( P, P ′ ) , (49b)Tr[Γ µa Ξ − ( P )¯Γ νa Ξ + ( P ′ )] = − g T X P,P ′ Tr (cid:2) γ µ γ Λ + p γ ν γ Λ − p ′ (cid:3) w a ( P, P ′ ) , (49c)Tr[¯Γ µa Ξ + ( P )Γ νa Ξ − ( P ′ )] = − g T X P,P ′ Tr (cid:2) γ µ γ Λ − p γ ν γ Λ + p ′ (cid:3) w a ( P, P ′ ) . (49d)where the repeated color indexes on LHS are not to be summed. The quantities w ± and w on RHS of Eqs.(49a-49d)are given by w ± a = p ± ǫ p p − ε p p ′ ± ǫ p ′ p ′ − ε p ′ , a = 1 , , h p ± ǫ p p − ε p p ′ ± ǫ p ′ p ′ − ǫ p ′ + p ± ǫ p p − ǫ p p ′ ± ǫ p ′ p ′ − ε p ′ i , a = 4 , · · · , p ± ǫ p p − ε p p ′ ± ǫ p ′ p ′ − ε p ′ + 13 p ± ǫ p p − ǫ p p ′ ± ǫ p ′ p ′ − ǫ p ′ , a = 8 (50a)and w a = − φ ( P ) φ ( P ′ )( p − ε p )( p ′ − ε p ′ ) , a = 1 , , , a = 4 , · · · , φ ( P ) φ ( P ′ )( p − ε p )( p ′ − ε p ′ ) , a = 8 (50b)where ǫ p = p − µ and ε p = p ( p − µ ) + φ ( P ). Since the dominant contributions in the weak coupling arise from thequasiparticles, we have ignored the contributions from the quasi-antiparticles in the calculations above. The traceover Dirac space is straightforwardTr (cid:2) γ Λ ± p γ γ Λ ± p ′ γ (cid:3) = − Tr (cid:2) γ γ Λ ± p γ γ Λ ∓ p ′ (cid:3) = 1 + ˆ p · ˆ p ′ , (51a) X i Tr (cid:2) γ i Λ ± p γ γ i Λ ± p ′ γ (cid:3) = X i Tr (cid:2) γ i γ Λ ± p γ i γ Λ ∓ p ′ (cid:3) = 3 − ˆ p · ˆ p ′ , (51b)Tr (cid:2) ~γ · ~k Λ ± p γ ~γ · ~k Λ ± p ′ γ (cid:3) = Tr (cid:2) ~γ · ~kγ Λ ± p ~γ · ~kγ Λ ∓ p ′ (cid:3) = (1 + ˆ p · ˆ p ′ ) ( p − p ′ ) k (51c)It can be shown that the contribution from Eq.(50a) to F is suppressed by an order g relative to that from Eq.(50b)and will be ignored here. We neglect also the dependence of the gap function on the magnitude of the momentum,but keep the dependence on the energy and the momentum orientation. Then the integrals over p and p ′ can becarried out easily. Since we are only interested in the zero temperature, the Matsubara sum becomes an integral overthe Euclidean energy. We findTr[ D n δ Π] = − g µ π Z dν Z dν ′ Z d ˆ p Z d ˆ p ′ φ ( ν, ˆ p ) φ ( ν ′ , ˆ p ′ ) q ( ν + φ ( ν, ˆ p ))( ν ′ + φ ( ν ′ , ˆ p ′ )) × (cid:2) D l ( ν − ν ′ , θ ) + D t ( ν − ν ′ , θ ) (cid:3) (52)where cos θ = ˆ p · ˆ p ′ . Making use of the Nambu-Gorkov formalism in Eq.(33-39), the rest terms of the condensateenergy density Eq.(32) can be evaluated readily.1Ω Trln S − − Trln S − = − µ (2 π ) Z dν Z d ˆ p h | ν | − p ν + φ ( ν, ˆ p ) i (53)1Ω Tr( S − S −
1) = − µ (2 π ) Z dν Z d ˆ p φ ( ν, ˆ p ) p ν + φ ( ν, ˆ p ) (54)The final expression of the condensation energy density reads F = − g µ π Z dν Z dν ′ Z d ˆ p Z d ˆ p ′ V ( ν − ν ′ , θ ) φ ( ν, ˆ p ) φ ( ν ′ , ˆ p ′ ) q ( ν + φ ( ν, ˆ p ))( ν ′ + φ ( ν ′ , ˆ p ′ ))+ 2 µ (2 π ) Z dν Z d ˆ p ν p ν + φ ( ν, ˆ p ) (55)where V contains the contribution from both magnetic and electric gluons, i.e. V = D l ( ν − ν ′ , θ ) + D t ( ν − ν ′ , θ ) (56)The gap equation can be derived by minimizing F with respect to the gap function φ ( ν, ˆ p ), δ Γ δφ = 0 (57)and we end up with φ ( ν, ˆ p ) = g µ π Z dν ′ Z d ˆ p ′ V ( ν − ν ′ , θ ) φ ( ν ′ , ˆ p ′ ) p ν ′ + φ ( ν ′ , ˆ p ′ ) (58)A consistent derivation of the gap equation up to the subleading order requires both the contribution from w ± a andthat from the diagonal block of (35) to be kept. The net result is to replace the first term inside the square root onRHS of (58) by ν ′ /Z ( ν ′ ) with Z ( ν ) the wave function renormalization of the normal phase. But it will not interferewith the angular dependence of the gap function to the subleading order as will be shown in the next section. IV. THE INTEGRAL EQUATION FOR THE ANGULAR DEPENDENCE OF THE GAP
Although the pairing strength are equal to the leading order of the QCD running coupling constant, similar tothe case II of the toy model, the subleading terms fall off with an increasing J . This makes the solution to the gapequation (58) highly nontrivial. In what follows, we shall isolate the energy dependence and the angle dependence ofthe pairing potential V . A differential equation with respect to the Matsubara energy will be derived from (58) thatfixes the gap function up to an arbitrary function of the angle. This function will be determined then by (58) with φ a known function of the Matsubara energy.Proceeding with the partial wave analysis, we expand V ( ν − ν ′ , θ ) in series of Legendre polynomials[22]: V ( ν − ν ′ , θ ) = 16 µ ln ω c | ν − ν ′ | X J (2 J + 1) P J (cos θ ) + 12 µ X J> (2 J + 1) c J P J (cos θ ) (59)where ω c = √ π µN f g and c J is given by Eq.(2). Using the completeness relation X J (2 J + 1) P J (cos θ ) = 4 πδ (ˆ p − ˆ p ′ ) (60)and the identity (proved in the Appendix A) Z d ˆ p ′ ∞ X J =1 (2 J + 1) c J P J (ˆ p · ˆ p ′ ) f (ˆ p ′ ) = 2 Z d ˆ p ′ f (ˆ p ′ ) − f (ˆ p ) | − ˆ p · ˆ p ′ | (61)with f (ˆ p ) an arbitrary function of ˆ p , the gap equation (58) becomes φ ( ν, x ) = ¯ g ω Z dν ′ n (cid:0) ln ω c | ν − ν ′ | + ln ω c | ν + ν ′ | (cid:1) φ ( ν ′ , x ) p ν ′ + φ ( ν ′ , x )+3 Z − dx ′ | x − x ′ | h φ ( ν ′ , x ′ ) p ν ′ + φ ( ν ′ , x ′ ) − φ ( ν ′ , x ) p ν ′ + φ ( ν ′ , x ) io (62)where ¯ g = g / (18 π ), x = ˆ p · ˆ z with ˆ z a fixed spatial direction and a UV cutoff, ω ∼ gµ is introduced. In deriving(62), we have assumed that the gap depends on x only, so the integration over the azimuthal angle of ˆ p ′ can be carriedout explicitly. The gap equation (62) can be further simplified by using the approximation of Son [17]ln ω c | ν − ν ′ | ≃ ln ω c | ν > | (63)with ν > = max( ν, ν ′ ). It is convenient to introduce ξ = ln ω c ν , a = ln ω c ω . (64)On writing φ = φ ( ξ, x ) and Φ( ξ, x ) ≡ ¯ g ∞ Z ξ dξ ′ φ ( ξ ′ , x ) q φ ( ξ ′ ,x ) ω c e ξ ′ (65)the gap equation(62) becomes φ ( ξ, x ) = ξ Φ( ξ, x ) − ξ Z a dξ ′ ξ ′ d Φ dξ ′ + 3 Z − dx ′ Φ( a, x ′ ) − Φ( a, x ) | x − x ′ | (66)Taking the derivative of both sides with respect to ξ , we find dφdξ = Φ( ξ, x ) (67)0which implies the boundary condition dφdξ → ξ → ∞ for all x . Another derivative of (67) yield the ordinary differential equation d φdξ + ¯ g φ q φ ω c e ξ = 0 (69)which is universal for all x . It follows from Eq.(66) that the gap equation is equivalent to a Φ( a, x ) − φ ( a, x ) + 3 Z − Φ( a, x ′ ) − Φ( a, x ) | x − x ′ | = 0 (70)The solution to (69) subject to the condition (68) contains an arbitrary function of x to be determined by (70). Nofurther approximation has been made up to now.The solution to the differential equation (69) proceeds in the same way as that for a spherical gap. To the leadingorder, the equation can be approximated by a linear one, d φ (0) dξ + ¯ g θ ( b − ξ ) φ (0) = 0 (71)where b ( x ) is to be determined by the condition φ ( b,x ) ω c e b = 1. Its solution that satisfies the boundary condition(68)and the continuity up to the first order derivative reads φ (0) ( ξ, x ) = φ ( x ) cos ¯ g [ b ( x ) − ξ ] , for ξ < b ( x ) ,φ ( x ) , for ξ ≥ b ( x ) . (72)where b ( x ) = ln ω c | φ ( x ) | (73)It follows from Eq.(67) then thatΦ( ξ, x ) = ¯ gφ ( x ) sin ¯ g [ b ( x ) − ξ ] , for ξ < b ( x ) , , for ξ ≥ b ( x ) . (74)The angle dependent factor f (ˆ p ) introduced in Eq.(3) is defined by f ( x ) ≡ φ ( x )∆ = O (1) (75)where ∆ is the s -wave gap given by π − ¯ g ln 2 ω c ∆ = 0 . (76)where the contribution from the wave-function renormalization is ignored. Up to the subleading order, the differentialequation(69) reads d φ (1) dξ + ¯ g θ ( b − ξ ) φ (1) = ¯ g h θ ( b − ξ ) − q φ ω c e ξ i φ (0) . (77)We find that φ (1) ( ξ, x ) = φ (1) ( ξ, x ) + A ( ξ, x ) u ( ξ, x ) − B ( ξ, x ) v ( ξ, x ) , (78)1where u ( ξ, x ) and v ( ξ, x ) are the two linearly independent solutions to the Eq.(71), u ( ξ, x ) = cos ¯ g [ b ( x ) − ξ ] , for ξ < b ( x ) , , for ξ ≥ b ( x ) . (79) v ( ξ, x ) = − sin ¯ g [ b ( x ) − ξ ] , for ξ < b ( x ) , ¯ gξ, for ξ ≥ b ( x ) . (80)and A ( ξ, x ) = ¯ g ∞ Z ξ dξ ′ h θ ( b − ξ ′ ) − q φ ( ξ ′ ,x ) ω c e ξ ′ i v ( ξ ′ , x ) φ (0) ( ξ ′ , x ) , (81a) B ( ξ, x ) = − ¯ g ∞ Z ξ dξ ′ h θ ( b − ξ ′ ) − q φ ( ξ ′ ,x ) ω c e ξ ′ i u ( ξ ′ , x ) φ (0) ( ξ ′ , x ) . (81b)At the point ξ = a , we have A ( a, x ) = 1 + O (¯ g ) , B ( a, x ) ≃ ¯ g ln2 (82)Therefore φ (1) ( a, x ) ≃ φ ( x ) h cos ¯ g ( b − a ) − ¯ g ln2 sin ¯ g ( b − a ) i (83)to the subleading order. Since ¯ g ( b − a ) = π/ O ( g ) according to Eq.(76), we haveΦ( a, x ) = ¯ gφ ( x ) + O ( g ) , (84)and φ (1) ( a, x ) ≃ φ ( x ) h π − ¯ g ( b − a ) − ¯ g ln2 i + O ( g ) (85)Substituting Eqs.(84) and (85) into Eq.(70), we obtain the gap equation to the subleading order − h π − ¯ g ln ω c | φ ( x ) | − ¯ g ln2 i φ ( x ) + Z − dx ′ φ ( x ′ ) − φ ( x ) | x − x ′ | = 0 (86)Then the integral equation for f ( x ), f ( x )ln | f ( x ) | − Z − dx ′ f ( x ′ ) − f ( x ) | x − x ′ | = 0 (87)follows from (76).Few comments are in order: 1) The spherical gap, f ( x ) = 1 is a trivial solution to Eq.(87) and there is no angularmomentum mixing. 2) The ”spin-1” gap, carried over from the polar phase of He , f ( x ) ∝ x , fails to satisfy thisequation. 3) Eq.(87) conserves the parity. In another word, its solution can be either an even or an odd functionof x . 4) If the wavefunction renormalization is restored, there will be an additional subleading term on RHS of (77)and an additional subleading term on RHS of B of eq.(82). This term, when substitute into Eq.(70), will cancel thecorresponding contribution to ∆ ∆ = πe − γ T (0) c = 2048 π µN f g e − π √ g − π . (88)leaving the integral equation (87) intact.2 V. THE NUMERICAL RESULTS OF THE ANGULAR DEPENDENCE
The solution to the integral equation Eq.(87) can be obtained from a variational principle. Upon substitutionof Eq.(72) with φ ( x ) = ∆ f ( x ) into Eq.(55), the condensate energy density becomes a functional of f (details inAppendix B), i. e. F = µ ∆ π F [ f ] , (89)where F [ f ( x )] = Z − dxf ( x ) h ln | f ( x ) | − i + 32 Z − dx Z − dx ′ (cid:2) f ( x ) − f ( x ′ ) (cid:3) | x − x ′ | = 2 Z dxf ( x ) h ln | f ( x ) | − i + 3 Z dx Z dx ′ n (cid:2) f ( x ) − f ( x ′ ) (cid:3) | x − x ′ | + (cid:2) f ( x ) + f ( x ′ ) (cid:3) x + x ′ o (90)with the last equality following from the odd parity of f ( x ), i. e. f ( − x ) = − f ( x ). Readers may easily verify that thevariational minimum of Eq.(90) does solve Eq.(87).Before the numerical solution, we consider a trial function f ( x ) = cx (91)and substitute it into the target functional (90). The minimization yields c = e − = e − ≃ . × − , (92)at which F [ f ] ≃ − . × − . (93)The trial function (91) is what people carried over from the polar phase of He . The ”-6” of the exponent of (92)comes from the pairing strength of the p -wave and the ” ” stems from the logarithm of (87). The latter contributionwas reported in [13]. The trial function (91) with p -wave alone is not optimal. The free energy will be lowered furtherby including higher partial waves of odd J as we shall see.To find the variational minimum, we discretize the integral of Eq.(90) by dividing the domain x ∈ (0 ,
1) into N ( >>
1) intervals with x j = ( j + 12 )∆ x, j = 0 , , , · · · , N − x = 1 /N . We have then F = lim N →∞ F N with F N = 2∆ x X j f j (cid:16) ln f j − (cid:17) + 6∆ x X j f j x j + 3∆ x X i,j,i = j h ( f i − f j ) | x i − x j | + ( f i + f j ) x i + x j i (95)where we have dropped the limit x ′ → x of the first term inside the curly bracket of Eq. (90). F N is a function of N variables. The stationary condition ∂ F ∂f j = 0 (96)yields f j ln f j + 3∆ x h x j + X i = j (cid:16) | x i − x j | + 1 x i + x j (cid:17)i f j − x X i = j (cid:16) | x i − x j | + 1 x i + x j (cid:17) f i = 0 (97)which is a discrete version of Eq. (87). Regarding f i ’s as given, the equation for f j is of the form(ln f j + a ) f j − b = 0 (98)3 J = 1 J = 3 J = 5 b J . × − − . × − . × − TABLE I: The first three expansion coefficients of the gap function according to Legendre polynomials. f ( x ) x=cos FIG. 2: The angular dependence of the gap function with angular momentum mixing. The dashed line and the solid line arethe initial configuration and the final numerical results respectively. with a and b positive. It has one and only one solution for f j > f j = e − x j as an initial configuration and update each f j by solving Eq.(98).This way we lower the value of the target functional F in each step and approach the solution to (87) eventually. Theprocess converges rapidly and our numerical solution to (87) is shown as solid line in Fig.2, which depart from thetrial function (dashed line) slightly. We find the minimum value of the target functional F [ f ] ≃ − . × − . (99)which drops from (93) by 3 . f ( x ) = X J =odd b J P J ( x ) . (100)The coefficients of the first three partial waves, J = 1 , ,
5, calculated by substituting the numerical solution into theformula b J = 2 J + 12 Z − dxf ( x ) P J ( x ) (101)are displayed in Table I. While the gap function contains all partial waves of odd J , the component of J = 1 is thebiggest. This is anticipated because the pairing strength of the all partial waves are equal in leading order but fall offwith an increasing J in the subleading order as is shown in the partial wave expansion (59).4 VI. CONCLUDING REMARKS
In summary we have explored the angular dependence of the gap function for a non-spherical pairing of CSC.Because of the equal strength of the pairing potential mediated by one-gluon exchange for all partial waves to theleading order of QCD running coupling constant and the nonlinearity of the gap equation, a non-spherical gap functioncannot be restricted to one angular momentum channel only. Other multipoles are bound to show up, which rendersthe angular dependence of the gap nontrivial. On the other hand, the pairing strength to the subleading orderdecreases with increasing angular momentum J . The mixing effect will not be as big as that in the soluble toy modelwe introduced for the purpose of clarification.For the single flavor CSC, we worked out the angular momentum mixing effect explicitly for the gap functionwith zero azimuthal quantum number at zero temperature. An nonlinear integral equation for the nontrivial angulardependence was derived and its solution was obtained numerically. The gap function in this case reads φ = ∆ f (ˆ p · ˆ z ) cos ¯ g (cid:16) ln ν ∆ | f (ˆ p · ˆ z ) | (cid:17) , for ν > ∆ | f (ˆ p · ˆ z ) | , ∆ f (ˆ p · ˆ z ) , for ν ≤ ∆ | f (ˆ p · ˆ z ) | . (102)where ∆ is given by Eq.(88) and f (ˆ p · ˆ z ) is plotted in Fig.2.The drop of the free energy of the modified polar phase by the mixing, however, is numerically small. The magnitudeof its condensation energy is smaller than that of the CSL phase by a factor of 1.48 instead of the factor 1.54 reportedin [13]. The CSL phase remains stable. In this sense our results at the moment is of theoretical values only. Thereare many other candidate pairing states between quarks of the same flavor [12, 13]. Among them are the states witha nonzero azimuthal quantum number and the pairing between quarks of opposite helicities. The former is analogousto the A phase of He and may be present in a compact star with a strong magnetic field. The pairing force in theunequal-helicity channel is stronger [12, 13, 22]. The angular momentum mixing effect is generic in all nonsphericalpairing states and the integral equation (87) can be readily generalized to these cases. There may be phenomenologicalimplications of the angular momentum mixing. A systematic survey of the angular momentum mixing effect in all”spin-1” CSC states covered in [13] will be reported in another paper.Another place where the angular momentum mixing shows up is the CSC-LOFF state in the presence of Fermimomentum mismatch. It has been speculated [28] that the forward singularity will increase the upper limit of themismatch value that supports a LOFF pairing. The new threshold was found in [29], motivated by the nearly equalpairing strength of all partial wave channels. The same mechanism works for the gap equation of LOFF pairing. Itsfree energy will be lowered by the angular momentum mixing and the lower edge of the LOFF window is expected tobe shifted to a lower value of the mismatch parameter. Acknowledgments
We would like to extend our gratitude to D. Rischke, T. Sch¨ a fer and A. Schmitt for stimulating discussions andvaluable comments. We are also benefitted from conversations with J.R. Li and Q. Wang. The work of D. F. H. andH. C. R. is supported in part by NSFC under grant No. 10575043 and by US Department of Energy under grantsDE-FG02-91ER40651-TASKB. The work of D. F. H. is also supported in part by Educational Committee of Chinaunder grant NCET-05-0675 and project No. IRT0624 APPENDIX A: THE DERIVATION OF EQUATION (61)
The integral formula of c J is[22] c J = Z − dx P J ( x ) − − x (A1)It is convenient to introduce c ǫJ = Z − dx P J ( x ) − − x + ǫ = Z − dx P J ( x )1 − x + ǫ − ln 2 + ǫǫ (A2)5where ǫ ( >
0) is an infinitesimal quantity. We have lim ǫ → + c ǫJ = c J . For the first term on RHS, we expand11 − x + ǫ = X J a J P J ( x ) (A3)according to a J = 2 J + 12 Z − dx P J ( x )1 − x + ǫ (A4)Therefore Eq. (A2) reads c ǫJ = 22 J + 1 a J − ln 2 + ǫ ∞ X J =1 (2 l + 1) c ǫJ P J (ˆ p · ˆ p ′ ) = ∞ X J =1 a J P J (ˆ p · ˆ p ′ ) − ∞ X J =1 (2 J + 1) P J (ˆ p · ˆ p ′ )ln 2 + ǫǫ = 21 − ˆ p · ˆ p ′ + ǫ − πδ (ˆ p − ˆ p ′ )ln 2 + ǫ f (ˆ p ), Z d ˆ p ′ ∞ X J =1 (2 J + 1) c ǫJ P J (ˆ p · ˆ p ′ ) f (ˆ p ′ ) = 2 Z d ˆ p ′ f (ˆ p ′ )1 − ˆ p · ˆ p ′ + ǫ − πf (ˆ p )ln 2 + ǫǫ = 2 Z d ˆ p ′ f (ˆ p ′ ) − f (ˆ p )1 − ˆ p · ˆ p ′ + ǫ + 2 Z d ˆ p ′ f (ˆ p )1 − ˆ p · ˆ p ′ + ǫ − πf (ˆ p )ln 2 + ǫǫ = 2 Z d ˆ p ′ f (ˆ p ′ ) − f (ˆ p )1 − ˆ p · ˆ p ′ + ǫ (A7)The Eq. (61) is obtained by taking the limit ǫ → + . APPENDIX B: THE CONDENSATION ENERGY DENSITY WITH THE ANGULAR MOMENTUMMIXING
In this appendix, we shall derive the expression Eq.(55)) of the condensation energy density with the angularmomentum mixing. Substituting Eq.(59) into the first term of Eq. (55), we find F = − g µ π Z dν ′ Z dν Z d ˆ p ′ Z d ˆ p h µ ln ω c | ν − ν ′ | ∞ X l =0 (2 l + 1) P l (ˆ p · ˆ p ′ )+ 12 µ ∞ X l =1 (2 l + 1) c l P l (ˆ p · ˆ p ′ ) i φ ( ν, ˆ p ) φ ( ν ′ , ˆ p ′ ) p ( ν + φ ( ν, ˆ p ))( ν ′ + φ ( ν ′ , ˆ p ′ ))= − ¯ g µ π n Z dν ′ Z dν Z d ˆ p ln ω c | ν − ν ′ | φ ( ν, ˆ p ) φ ( ν ′ , ˆ p ) p ( ν + φ ( ν, ˆ p ))( ν ′ + φ ( ν ′ , ˆ p )) − π Z dν ′ Z dν Z d ˆ p ′ Z d ˆ p − ˆ p ′ · ˆ p h φ ( ν, ˆ p ) p ν + φ ( ν, ˆ p ) − φ ( ν, ˆ p ′ ) p ν + φ ( ν, ˆ p ′ ) ih ( ν ↔ ν ′ ) io (B1)Because of the eveness of φ ( ν, ˆ p ) in ν , we have F = − ¯ g µ π n Z d ˆ p ω Z dν ′ ω Z ln ω c ν > φ ( ν, ˆ p ) φ ( ν ′ , ˆ p ) p [ ν + φ ( ν, ˆ p )][ ν ′ + φ ( ν ′ , ˆ p )] − π Z d ˆ p ′ Z d ˆ p ω Z dν ′ ω Z dν − ˆ p · ˆ p ′ h φ ( ν, ˆ p ) p ν + φ ( ν, ˆ p ) − φ ( ν, ˆ p ′ ) p ν + φ ( ν, ˆ p ′ ) ih ( ν ↔ ν ′ ) io (B2)6where the approximation (63) has been applied to the forward logarithm. For the gap function of zero azimuthalquantum number, φ ( ν, ˆ p ) depends only on x ≡ ˆ p · ˆ z . We find that F = − µ π ¯ g n Z − dx ∞ Z a dξ ∞ Z a dξ ′ ξ < d Φ( ξ, x ) dξ d Φ( ξ ′ , x ) dξ ′ − Z − dx Z − [Φ( a, x ′ ) − Φ( a, x )] | x − x ′ | o (B3)where, Φ( ξ, x ) has been defined in Eq.(65) and ξ and a have been defined in (64). The integral over ξ ′ followed bythe integral by part over ξ leads to F = µ π ¯ g n − a Z − dx Φ ( a, x ) − Z − dx ∞ Z a dξ Φ ( ξ, x ) + 32 Z − dx Z − dx ′ [Φ ( a, x ′ ) − Φ ( a, x )] | x − x ′ | o (B4)Making use of Eq. (65) and (85), we have ∞ Z a dξ Φ ( ξ, x ) = − ¯ g φ ( x ) h π − ¯ g ( b − a ) − ¯ g ln2 i + ¯ g ∞ Z a dξ φ ( ξ, x ) q φ ( ξ,x ) ω c e ξ (B5)and thus − a Φ ( a, x ) − ∞ Z a dξ Φ ( ξ, x ) = ¯ gφ ( x ) (cid:16) π − ¯ gb − ¯ g ln2 (cid:17) − ¯ g ω Z dν φ ( ν, x ) p ν + φ ( ν, x ) (B6)Substituting φ ( x ) = ∆ f ( x ) into (B4), we obtain that F = µ ∆ π n Z − dxf ( x )ln | f ( x ) | + 32 Z − dx Z − dx ′ (cid:2) f ( x ) − f ( x ′ ) (cid:3) | x − x ′ | o − µ π ω Z dν φ ( ν, x ) p ν + φ ( ν, x ) (B7)Then the condensate energy density with the angular momentum mixing reads F = F + µ π Z dν Z d ˆ p h | ν | − φ ( ν, ˆ p ) p ν + φ ( ν, ˆ p ) i = µ ∆ π n Z − dxf ( x ) h ln | f ( x ) | − i + 32 Z − dx Z − dx ′ (cid:2) f ( x ) − f ( x ′ ) (cid:3) | x − x ′ | o (B8)The minimization of this free energy give rise to Eq.(87) of the text. [1] B. Barrois, Nucl. Phys. B129 , 390(1977); S. C. Frautschi, in
Hadronic matter at extreme energy density , edited by N.Cabibbo and L. Sertorio (Plenum Press, 1980)[2] D. Bailin, and A. Love, Phys. Rep. , 325 (1984).[3] M. Alford, K. Rajagopal and F. Wilczek, Phys. Lett.
B422 , 247 (1998).[4] R. Rapp, T. Sch¨ a fer, E. V. Shuryak, and M. Velkovsky, Phys. Rev. Lett. , 53 (1998).[5] M. Alford, K. Rajagopal, and F. Wilczek, Nucl. Phys. B537 , 443 (1999)[6] Dirk H. Rischke , Prog. Part. Nucl. Phys. 52:197-296, 2004[7] Mei Huang, Pengfei Zhuang and Weiqin Chao, Phys. Rev.
D67 , 065015 (2003)[8] I. shovkovy and M. Huang, Phys. Lett.
B564 , 205 (2003)[9] M. Alford, C. Kouvaris, and K. Rajagopal, Phys. Rev. Lett. , 222001 (2004).[10] M. G. Alford, J. A. Bowers and K. Rajagopal, Phys. Rev. D63 , 074061 (2001)[11] P. F. Befaque, H. Caldas and G. Rupak, Phys. Rev. Lett. , 247002 (2003)[12] T. Sch¨ a fer, Phys. Rev. D62 , 094007 (2000)[13] A. Schmitt, Phys. Rev.
D71 , 054016 (2005); nucl-th/0405076 (Ph. D. thesis).[14] A. Schmitt, I. A. Shovkovy and Q. Wang, Phys. Rev. Lett. [15] A. Schmitt, Q. Wang and D. H. Rischke, Phys. Rev. Lett. D59 , 094019 (1999).[18] T. Sch¨ a fer and F. Wilczek, Phys. Rev. D60 , 114033 (1999).[19] R. D. Pisarski, and D. H. Rischke, Phys. Rev.
D61 , 074017 (2000).[20] D. K. Hong, V. A. Miransky, I. A. Shovkovy, and L. C. R. Wijewardhana, Phys. Rev.
D61 , 056001 (2000).[21] W. E. Brown, J. T. Liu, and H-C. Ren, Phys. Rev.
D61 , 114012 (2000).[22] W. E. Brown, J. T. Liu, and H-C. Ren, Phys. Rev.
D62 , 054016 (2000).[23] G. D. Mahan, Many Particle Physics, 1990, Plenum Press, New York.[24] W. E. Brown, J. T. Liu, and H-C. Ren, Phys. Rev.
D62 , 054013 (2000).[25] Qun Wang and Dirk H. Rischke, Phys. Rev.
D65 , 054005 (2002).[26] B. Feng, D-F. Hou, J-R. Li, and H-C. Ren, Nucl. Phys.
B754 , 351 (2006)[27] P. T. Reuter, Phys. Rev.
D74 , 105008 (2006)[28] A. K. Leibovich, K. Rajagopal and E. Shuster, Phys. Rev.
D64 , 094005 (2001).[29] I. Giannakis, J. T. Liu and H. C. Ren, Phys. Rev.
D66 , 031501 (2002)[30] J. M. Cornwall, R. Jackiw, and E. Tomboulis, Phys. Rev.