Ginzburg-Landau free energy of crystalline color superconductors: A matrix formalism from solid-state physics
aa r X i v : . [ nu c l - t h ] D ec Ginzburg-Landau free energy of crystalline color superconductors:A matrix formalism from solid-state physics
Gaoqing Cao and Lianyi He ∗ The Ginzburg-Landau (GL) free energy of crystalline color superconductors is important for understandingthe nature of the phase transition to the normal quark matter and predicting the preferred crystal structure. Sofar the GL free energy at zero temperature has only been evaluated up to the sixth order in the condensate. Togive quantitative reliable predictions we need to evaluate the higher-order terms. In this work, we present a newderivation of the GL free energy by using the discrete Bloch representation of the fermion field. This derivationintroduces a simple matrix formalism without any momentum constraint, which may enable us to calculate theGL free energy to arbitrary order by using a computer.
PACS numbers: 12.38.-t, 21.65.Qr, 74.20.Fg, 03.75.Hh
I. INTRODUCTION
It is generally believed that the inhomogeneous Larkin-Ovchinnikov-Fulde-Ferrell (LOFF) phase appears in a super-conductor when the pairing between di ff erent fermion speciesis under the circumstances of mismatched Fermi surfaces [1–3]. The mismatched Fermi surfaces are normally inducedby the Zeeman energy splitting 2 δµ in a magnetic field [4].Early studies of the LOFF phase were restricted to 1D struc-tures, which include the Fulde-Ferrell (FF) state [2] with aplane-wave order parameter ∆ ( z ) = ∆ e iqz and the Larkin-Ovhinnikov (LO) state [1] with an antipodal-wave order pa-rameter ∆ ( z ) = ∆ cos(2 qz ). For s -wave pairing at weak cou-pling, it is known that the FF or LO state exists in a narrowwindow δµ < δµ < δµ , where the lower critical field δµ = . ∆ and the upper critical field δµ = . ∆ [1, 2] with ∆ being the pairing gap at vanishing mismatch. However,since the thermodynamic critical field is much lower than δµ due to strong orbit e ff ect, it is rather hard to observe the LOFFstate in ordinary superconductors [4]. In recent years, experi-mental evidences for the LOFF state in some superconductingmaterials have been reported [5–8].The studies of color superconductivity in dense quark mat-ter promoted new interests in the LOFF state [9–21]. Colorsuperconductivity in dense quark matter appears due to the at-tractive interactions in certain diquark channels [22–26]. Un-der the compact star constraints, di ff erent quark flavors ( u , d ,and s ) acquire mismatched Fermi surfaces because of the Betaequilibrium and the electric charge neutrality. On the otherhand, recent experiments on ultracold atomic Fermi gases pro-vided a controllable way to study the fermion superfluiditywith population imbalance [27, 28]. Quark color supercon-ductors under compact star constraints as well as atomic Fermigases with population imbalance are rather clean systems torealize the long-sought LOFF phase.In addition to the simple FF and LO states, there exist a ∗ E-mail address: [email protected]; Present address: Department of Physics,Tsinghua University, Beijing 100084, China variety of crystal structures. The general form of the crystalstructure of the order parameter can be expressed as [3] ∆ ( r ) = P X k = ∆ e iq ˆ n k · r . (1)A specific crystal structure corresponds to a multi-wave con-figuration determined by the P unit vectors n k ( k = , , ..., P ).Around the tricritical point in the temperature-mismatchphase diagram, the LOFF phase can be studied rigorously byusing the Ginzburg-Laudau (GL) analysis since both the gapparameter ∆ and the pair momentum q are vanishingly small[3]. It was found that the LO state is preferred near the tricriti-cal point [29–31]. However, the real ground state of the LOFFphase is still not quite clear due to the limited theoretical ap-proaches at zero temperature. Various theoretical approachessuggested that the LOFF state has a complicated crystal struc-ture near its phase transition to the normal state [10, 13, 32–35]. A recent self-consistent treatment of the 1D modulationshowed that a solitonic lattice structure is preferred near thephase transition to the BCS state [34].The GL free energy is important for us to understand thenature of the phase transition to the normal state and to searchfor the most preferred crystal structure near the phase transi-tion point. In a pioneer work, Bowers and Rajagopal investi-gated 23 crystal structures at weak coupling by using the GLapproach [10]. They evaluated the GL free energy up to thesixth order in ∆ , Ω GL ( ∆ ) N = P α ∆ + β ∆ + γ ∆ + O ( ∆ ) , (2)where N is the density of state at the Fermi surface. Thecoe ffi cient α is universal for all crystal structures and is givenby [10] α ( δµ, q ) = − + δµ q ln q + δµ q − δµ ! −
12 ln ∆ q − δµ ) . (3)Near the conventional second-order phase transition point δµ = δµ with the optimal pair momentum q = . δµ ,we have α ≃ ( δµ − δµ ) /δµ . The GL approach is meaningfulif the phase transition is of second order or weak first order.The GL coe ffi cients β and γ for 23 crystal structures wereevaluated in [10]. For most of the structures, we find β < γ >
0, the favored one seems to be the body-centered cubic (BCC) structure with P =
6. The order pa-rameter can be expressed as ∆ ( r ) = ∆ (cid:2) cos(2 qx ) + cos(2 qy ) + cos(2 qz ) (cid:3) . (4)Further, it was conjectured that the face-centered cubic (FCC)structure with P = β and γ are negative and their absolute values are thelargest [10]. The order parameter can be expressed as ∆ ( r ) = ∆ cos qx √ ! cos qy √ ! cos qz √ ! . (5)For the BCC structure, the GL free energy up to the sixthorder in ∆ predicts a strong first-order phase transition at δµ = δµ ∗ ≃ . ∆ with a large gap parameter ∆ ≃ . ∆ at δµ . δµ ∗ [10]. The prediction of a strong first-order phasetransition indicates that the GL approach is not valid or thehigher-order terms in ∆ are important. For the FCC structure,the GL free energy up to the sixth order in ∆ actually gives noprediction because both β and γ are negative.On the other hand, there exist some other approaches thatdo not use the GL approximation. Combescot and Mora em-ployed Eilenberger’s quasiclassical equation with a Fourierexpansion for the order parameter [32, 33]. This approachpredicted that the BCC-normal phase transition is of ratherweak first order : The upper critical field δµ ∗ is only 3 . δµ and ∆ ≃ . ∆ at δµ . δµ ∗ [33]. For the FCCstructure, it was found that its upper critical field is only 1 . δµ and hence it is less favored than BCC. An-other study employed a solid-state-like approach which cal-culated the free energy of the BCC structure by directly di-agonalizing the Hamiltonian matrix in the Bloch space [35].This approach also predicted that the BCC-normal phase tran-sition is of weak first order and the upper critical field of BCCis slightly higher than δµ .The contradiction between the predictions from the aboveapproaches and from the GL approach indicates that thehigher-order terms in the GL free energy are rather importantfor a quantitative prediction. Actually, for a first-order phasetransition, the higher-order expansions are crucial to deter-mine whether the first-order transition is weak or strong. Foran intuitive understanding, let us add the eighth-order term inthe GL free energy for the BCC structure. We have Ω GL ( ∆ ) N = P α ∆ + β ∆ + γ ∆ + η ∆ + O ( ∆ ) . (6)As a naive example, we assume η ≥
0. For η =
0, we obtaina strong first-order phase transition. However, the first-orderphase transition becomes weaker and weaker if we increasethe value of η . For η → + ∞ , the phase transition approachessecond order and δµ ∗ → δµ . On the other hand, if η is smallor even negative, the higher-order terms such as O ( ∆ ) maybecome important. Therefore, to give more precise predic-tions within the GL approach we need to evaluate the GL free energy to a su ffi ciently high order in ∆ . Actually, if we canevaluate the GL free energy to arbitrary order and determinethe convergence properties of the GL series, we may even treatthe strong first-order phase transition within the GL approach.The higher-order GL coe ffi cients can be evaluated by usingthe diagrammatic approach used in [10]. In this approach, toevaluate the 2 k -th order GL coe ffi cient one needs to sum allpossible configurations that satisfy the momentum constraint k X i = ( − i + q i = (7)for a set of 2 k wave vectors { q q ... q k } with q i = q ˆ n i . Forlarge k , the number of the configurations becomes very largefor the crystal structures with large number of waves P . More-over, one needs to introduce 2 k Feynman parameters to evalu-ate the integrals. Therefore, the calculation of the higher-orderterms is tedious and complicated in this formalism.In this work, we introduce a new derivation of the GL freeenergy based on a solid-state physics approach to the crystalstructures. This derivation provides a simple matrix formal-ism for the GL coe ffi cients without any momentum constraint.We find that the formalism used in [10] and our formalism canbe attributed to two di ff erent representations of the fermionfield: The diagrammatic approach [10] employs the usual con-tinuous momentum representation, while our derivation usesthe discrete Bloch representation. Since the matrix operationscan be easily realized by using a computer, this new formal-ism may enable us to calculate the GL free energy to arbitraryorder in ∆ .The paper is organized as follows. In Sec. II we brieflyreview the diagrammatic approach to the GL free energy. InSec. III we present our new derivation of the GL free energyby using the solid-state physics approach. The explicit matrixformalism of the GL coe ffi cients for some crystal structuresare presented in Sec. IV. We summarize in Sec. V. The naturalunits ~ = k B = II. CONTINUOUS MOMENTUM REPRESENTATION
We focus on the general two-flavor pairing at high density,low temperature, and weak coupling. In this case, the antipar-ticle degrees of freedom play no role. Therefore, we can startfrom a general e ff ective Lagrangian for two-flavor pairing athigh density. The Lagrangian density is given by [3] L e ff = ψ † (cid:2) i ∂ t − ε ( ˆ p ) + ˆ µ (cid:3) ψ + g ψ † σ ψ ∗ )( ψ T σ ψ ) , (8)where ψ = ( ψ u , ψ d ) T denotes the two-flavor fermion field, ε ( ˆ p ) is the fermion dispersion with ˆ p = − i ∇ , g is a contactcoupling which represents the attractive interaction, and σ isthe second Pauli matrix in the flavor space. We shall neglectall other internal degrees of freedom, such as color and spin.These degrees of freedom contribute a simple degenerate fac-tor and can be absorbed into the definition of the density ofstate N at the Fermi surface. The fermion chemical poten-tials are specified by the diagonal matrix ˆ µ = diag( µ u , µ d ) inthe flavor space, where µ u = µ + δµ, µ d = µ − δµ. (9)Here δµ plays the role of the mismatch between the Fermisurfaces. For ultra-relativistic fermion systems such as densequark matter we take ε ( p ) = | p | . The model also works fornonrelativistic fermion systems such as cold atomic Fermigases. In this case we take ε ( p ) = p (we set the fermionmass M = / ∆ at δµ = ε F ≃ µ . At weak coupling, the physical resultsshould be universal if we properly scale the physical quanti-ties by using the pairing gap ∆ and the density of states N atthe Fermi surface.To evaluate the GL free energy, we start from the partitionfunction Z of the system. In the imaginary time formalism, itis given by Z = Z [ d ψ ][ d ψ † ] e −S e ff (10)with the Euclidean action S e ff = − Z / T d τ Z V d r L e ff . (11)Here τ = it is the imaginary time and T is the temperature.Once the partition function Z is evaluated, the free energydensity is given by Ω = − TV ln Z . (12)Fermionic superconductivity is characterized by a nonzero ex-pectation value of the difermion fields ϕ ( τ, r ) = − g ψ u ( τ, r ) ψ d ( τ, r ) (13)The order parameter ∆ ( r ) = h ϕ ( τ, r ) i is static but inhomoge-neous in the LOFF state. At weak coupling and at low tem-perature, the order parameter fluctuation becomes negligible and we can employ the mean-field approach. With the help ofthe Nambu-Gor’kov (NG) spinor Ψ ( τ, r ) = ψ u ( τ, r ) ψ ∗ d ( τ, r ) ! , (14)the mean-field Lagrangian can be expressed as L MF = Ψ † ( τ, r )( − ∂ τ − H MF ) Ψ ( τ, r ) − | ∆ ( r ) | g , (15)where the Hamiltonian operator is given by H MF = ε ( ˆ p ) − µ − δµ ∆ ( r ) ∆ ∗ ( r ) − ε ( ˆ p ) + µ − δµ ! . (16)For convenience we first consider a finite system in a cu-bic box defined as x , y , z ∈ [ − L / , L /
2] and then set L → ∞ .Imposing the periodic boundary condition, the fermion mo-mentum becomes discrete and is given by p = π L (cid:16) l e x + m e y + n e z (cid:17) , l , m , n ∈ Z . (17)To convert to the momentum space, we use the Fourier trans-formation for the fermion fields Ψ ( τ, r ) = √ V X p X ω n ˜ Ψ ( i ω n , p ) e − i ω n τ + i p · r , Ψ † ( τ, r ) = √ V X p X ω n ˜ Ψ † ( i ω n , p ) e i ω n τ − i p · r , (18)where V = L is the system volume and ω n = (2 n + π T ( n ∈ Z ) is the fermion Matsubara frequency. For the generalcrystal structure of the order parameter given by (1), the mean-field action in the momentum space can be evaluated as S MF = VT P ∆ g − T X ω n ,ω n ′ X p , p ′ ˜ Ψ † ( i ω n , p ) (cid:16) i ω n δ ω n ,ω n ′ δ p , p ′ − δ ω n ,ω n ′ H p , p ′ (cid:17) ˜ Ψ ( i ω n ′ , p ′ ) , (19)where p , p ′ are the discrete momenta given by (17) and theHamiltonian matrix H p , p ′ reads H p , p ′ = ( ξ p − δµ ) δ p , p ′ ∆ ( F + ) p , p ′ ∆ ( F − ) p , p ′ ( − ξ p − δµ ) δ p , p ′ ! . (20)Here and in the following ξ p = ε ( p ) − µ . The matrices F ± that characterizes the crystal structure are given by( F ± ) p , p ′ = P X k = δ p − p ′ , ± q k , (21)where q k = q n k .The partition function in the mean-field approximation canbe evaluated in the momentum representation. We have Z MF = Z [ d ˜ Ψ ][ d ˜ Ψ † ] e −S MF . (22)Completing the Gaussian integral, we obtain the free energy Ω = Pg ∆ − TV X ω n Trln " ( S − ) p , p ′ T , (23)where the Trln acts in the momentum space and the NG space.The inverse fermion propagator S − defined in these spaces isgiven by ( S − ) p , p ′ ( i ω n ) = i ω n δ p , p ′ − H p , p ′ . (24)The free energy can be evaluated if we can diagonalize theHamiltonian matrix H p , p ′ . However, this is infeasible becausethe momenta p and p ′ become continuous in the large vol-ume limit. We therefore turn to the GL expansion of the freeenergy. The standard field theoretical approach is to use thederivative expansion. First, we separate the“BCS” self-energy Σ ∆ and write ( S − ) p , p ′ = ( S − ) p , p ′ − ( Σ ∆ ) p , p ′ , (25)where the inverse of the free fermion propagator reads( S − ) p , p ′ = ( S − + ) p , p ′
00 ( S − − ) p , p ′ ! , (26)with the matrix elements given by( S − ± ) p , p ′ = ( i ω n + δµ ∓ ξ p ) δ p , p ′ . (27)The self-energy term Σ ∆ is given by( Σ ∆ ) p , p ′ = ∆ F p , p ′ , (28) where the matrix F reads F p , p ′ = F + ) p , p ′ ( F − ) p , p ′ ! . (29)Using the derivative expansion, we obtain Ω = Ω N + Pg ∆ + TV X ω n ∞ X l = ∆ l l Tr h ( S F ) l i , (30)where Ω N is the free energy of the normal state, Ω N = − TV X ω n Trln ( S − ) p , p ′ T . (31)To obtain the GL free energy, we first complete the trace inthe NG space. It is easy to show the trace vanishes for odd l .After some manipulation, we obtain the GL free energy Ω GL ( ∆ ) = α ∆ + ∞ X k = α k k ∆ k . (32)The second-order GL coe ffi cient α is given by α = Pg + TV X ω n Tr [ S + F + S − F − ] . (33)The higher-order GL coe ffi cients with k ≥ α k = TV X ω n Tr h ( S + F + S − F − ) k i . (34)Note that the trace is now taken only in the momentum space.Next we complete the trace in the momentum space andobtain the integral expressions of the GL coe ffi cients. For thesecond order, we haveTr [ S + F + S − F − ] = P X a = P X b = X p X p , p , p δ p , p i ω n + δµ − ξ p δ p − p , q a δ p , p i ω n + δµ + ξ p δ p − p , − q b = P X a = P X b = X p X p ′ i ω n + δµ − ξ p δ p − p ′ , q a i ω n + δµ + ξ p ′ δ p ′ − p , − q b = P X a = X p i ω n + δµ − ξ p i ω n + δµ + ξ p − q a . (35)Therefore, at zero temperature the second-order GL coe ffi cient α can be expressed as α P = g + Z ∞−∞ dE π Z d p (2 π ) iE + δµ − ξ p iE + δµ + ξ p − q . (36)We notice that α / P is universal for all crystal structures. The above integral su ff ers from ultraviolet (UV) divergence. In theweak coupling limit, the pairing and hence the momentum integral is dominated near the Fermi surface. It is convenient to usethe regularization scheme that the momentum integral is restricted near the Fermi surface; i.e., − Λ < | p | − µ < Λ . Using the factthat δµ, q ≪ Λ ≪ µ , we can express α as α = PN α ( δµ, q ), where α ( δµ, q ) = gN + Z ∞−∞ dE π Z Λ − Λ d ξ Z d ˆ p π iE + δµ − ξ iE + δµ + ξ − p · q (37)with N = µ / (2 π ) being the density of state at the Fermi surface and ˆ p denoting the solid angle. The integrals can nowbe analytically worked out and the cuto ff dependence can be removed by using the pairing gap at vanishing mismatch, ∆ = Λ e − / ( gN ) . We finally obtain the analytical expression given by (3).For higher orders ( k ≥ k = h ( S + F + S − F − ) i = P X a = P X b = P X c = P X d = X p X p , p ,..., p δ p , p i ω n + δµ − ξ p δ p − p , q a δ p , p i ω n + δµ + ξ p δ p − p , − q b δ p , p i ω n + δµ − ξ p δ p − p , q c δ p , p i ω n + δµ + ξ p δ p − p , − q d ! = P X a = P X b = P X c = P X d = δ q a − q b + q c − q d , X p i ω n + δµ − ξ p i ω n + δµ + ξ p − q a i ω n + δµ − ξ p − q a + q b i ω n + δµ + ξ p − q a + q b − q c ! . (38)At weak coupling, the 2 k -th order GL coe ffi cient ( k ≥
2) can be generally expressed as α k = N X q , q ,..., q k J k ( q q . . . q k ) δ q s , , q s = k X i = ( − i + q i . (39)Here the summation over each q i means the summation over all P wave vectors q a ( a = , , ..., P ). The quantity J k ( q q . . . q k )for general k can be expressed in a compact form, J k ( q q . . . q k ) = Z ∞−∞ dE π Z ∞−∞ d ξ Z d ˆ p π k Y i = iE + δµ − ξ + p · k i iE + δµ + ξ − p · l i , (40)where the momenta k i and l i are given by (we define q = for convenience) k i = i − X n = ( − n + q n , l i = i − X n = ( − n + q n . (41)Note that we have set Λ → ∞ since the integral over ξ is freefrom UV divergence for k ≥ q with v F q where v F = √ µ is the Fermivelocity. The density of state N at the Fermi surface is re-placed by N = √ µ/ (4 π ). Therefore, at weak coupling, theGL free energy is universal for both the relativistic case andthe nonrelativistic case once we properly express the free en-ergy in terms of the density of state N at the Fermi surface,the pairing gap ∆ at vanishing mismatch, and the quantity v F q ( v F = ffi cients up to the sixth order ( k ≤
3) for 23crystal structures have been evaluated numerically by Bow-ers and Rajagopal [10]. They introduced Feynman param-eters to evaluate the integral J k ( q q . . . q k ). The advan- tage of the above formalism obtained by using the contin-uous momentum representation is that one can directly ap-proach the weak coupling limit by using the momentum cut-o ff scheme − Λ < | p | − µ < Λ . However, for higher-ordercoe ffi cients with large k , the calculation becomes complicatedand tedious. First, for general k , one needs to introduce 2 k Feynman parameters x , x , ..., x k and y , y , ..., y k . At large k , the integral over the Feynman parameters becomes com-plicated. Second, to obtain the GL coe ffi cients, one needs tosum over all possible configurations that satisfy the constraint P ki = ( − i + q i = . At large k , the number of these configura-tions becomes also large, which makes the calculation tedious.To the best of our knowledge, no results of the higher-orderGL coe ffi cients ( k ≥
4) have been reported so far.
III. SOLID-STATE PHYSICS APPROACH: DISCRETEREPRESENTATION
In this section, we turn to a discrete representation inspiredby solid-state physics. For a specific crystal structure givenby (1), it is periodic in coordinate space. The order parametercan be alternatively expressed as ∆ ( r ) = ∆ f ( r ) , (42)where f ( r ) = P X k = e iq ˆ n k · r (43)is a periodic function. Here we assume that the function f ( r )corresponds to a 3D lattice structure. The derivation of theGL free energy can be easily generalized to 1D and 2D latticestructures. For a 3D lattice structure, the unit cell is generatedby three linearly independent vectors a , a , and a . We canaccordingly define the reciprocal space generated by three lin-early independent vectors b , b , and b with b i obtained bythe relation a i · b j = πδ i j . The periodicity of the order param-eter means f ( r ) = f ( r + a i ). Therefore, the function f ( r ) canbe decomposed into a discrete set of Fourier components. Wehave f ( r ) = X G f G e i G · r = ∞ X l , m , n = −∞ f lmn e i G lmn · r , f ∗ ( r ) = X G f ∗ G e − i G · r = ∞ X l , m , n = −∞ f ∗ lmn e − i G lmn · r , (44)where the reciprocal lattice vector G is given by G = G lmn = l b + m b + n b , l , m , n ∈ Z . (45)The Fourier components f G and f ∗ G are given by f G = V c Z c d r f ( r ) e − i G · r , f ∗ G = V c Z c d r f ∗ ( r ) e i G · r , (46)where V c is the volume of the unit cell and the integration isrestricted in the cell volume. It is easy to show that X G | f G | = P . (47)Since the order parameter or pair potential ∆ ( r ) is periodic,the eigenvalue equation for the fermionic excitation spectrum, which is known as the Bogoliubov-de Gennes (BdG) equa-tion, is analogous to the Schr¨odinger equation of a quantumparticle moving in a periodic potential. The BdG equation forthe present system can be expressed as ε ( ˆ p ) − µ − δµ ∆ ( r ) ∆ ∗ ( r ) − ε ( ˆ p ) + µ − δµ ! φ λ ( r ) = E λ φ λ ( r ) . (48)According to the Bloch theorem, the eigenfunction φ λ ( r ) takesthe form of the Bloch function; i.e., φ λ ( r ) = e i k · r φ λ k ( r ) , (49)where the lattice momentum k is restricted in the Brillouinzone (BZ) and the function φ λ k ( r ) has the same periodicityas the order parameter ∆ ( r ). We therefore have the similarFourier expansion φ λ k ( r ) = X G φ G e i G · r = ∞ X l , m , n = −∞ φ lmn e i G lmn · r . (50)Substituting this expansion into the BdG equation, we finallyobtain a matrix equation in the G -space, X G ′ H G , G ′ ( k ) φ G ′ = E λ ( k ) φ G , (51)where the Hamiltonian matrix H G , G ′ ( k ) is given by H G , G ′ ( k ) = ( ξ k + G − δµ ) δ G , G ′ ∆ f G − G ′ ∆ f ∗ G ′ − G ( − ξ k + G − δµ ) δ G , G ′ ! . (52)For a given momentum k in the BZ, we can solve the eigen-values E λ ( k ) by diagonalizing the above matrix in the discrete G -space. This means that the fermionic excitation spectrumforms a band structure, in analogy to the energy spectrum of aquantum particle moving in a periodic potential.Now we turn to the field theory. We consider a finite sys-tem spanned by three vectors N a , N a , and N a and as-sume periodic boundary condition. Then the system con-tains N N N unit cells and the thermodynamic limit can bereached by setting N i → ∞ . In accordance with the Fourierexpansion (44) and the matrix equation (51), we expand thefermion field Ψ ( τ, r ) in terms of the Bloch function rather thanusing the usual momentum representation (18). We write Ψ ( τ, r ) = √ V X k ∈ BZ e i k · r X ω n X G ˜ Ψ ( i ω n , k , G ) e − i ω n τ + i G · r , Ψ † ( τ, r ) = √ V X k ∈ BZ e − i k · r X ω n X G ˜ Ψ † ( i ω n , k , G ) e i ω n τ − i G · r . (53)One can recover the usual momentum representation (18) by using the fact that k + G can generate all possible momentum p inthe momentum space. This expansion defines a new representation with two di ff erent quantum numbers k and G . We call thisBloch representation. In the Bloch representation, the mean-field action can be evaluated as S MF = VT ∆ g X G | f G | − T X ω n ,ω n ′ X k , k ′ ∈ BZ X G , G ′ ˜ Ψ † ( i ω n , k , G ) h i ω n δ ω n ,ω n ′ δ k , k ′ δ G , G ′ − δ ω n ,ω n ′ δ k , k ′ H G , G ′ ( k ) i ˜ Ψ ( i ω n ′ , k ′ , G ′ ) , (54)where the Hamiltonian matrix H G , G ′ ( k ) is given by (52). Inmathematics, the Bloch representation corresponds to a simi-larity transformation of the usual momentum representation,which makes the Hamiltonian matrix H p , p ′ block diagonalwith the blocks characterized by the lattice momentum k .The functional path integral can be worked out by performingintegrals over ˜ Ψ † ( i ω n , k , G ) and ˜ Ψ ( i ω n ′ , k ′ , G ′ ) for all possi-ble values of { i ω n , k , G } and { i ω n ′ , k ′ , G ′ } . Since the inversefermion propagator is diagonal in the frequency space and the k -space, the free energy can be expressed as Ω = Pg ∆ − TV X ω n X k ∈ BZ Trln " ( S − ) G , G ′ T , (55)where the Trln acts in the G -space and the NG space. Theinverse fermion propagator S − defined in the G -space andthe NG space is given by( S − ) G , G ′ ( i ω n , k ) = i ω n δ G , G ′ − H G , G ′ ( k ) . (56)The free energy can be evaluated if we can diagonalize theHamiltonian matrix H G , G ′ ( k ) to obtain all the eigenvalues orthe band spectrum { E λ ( k ) } . This is in principle feasible be-cause the G -space is discrete. Since the matrix has infinitedimensions, we should make a truncation − D ≤ l , m , n ≤ D ( D ∈ Z + ) to perform the diagonalization. The size of the ma-trix we need to diagonalize is 2(2 D + . To achieve conver-gence to the limit D → ∞ we normally need a large cuto ff D ,which leads to a large computing cost [35].Then we turn to the GL expansion based on the expression(55) of the free energy. The procedure is the same as we usedin Sec. II. First, for the inverse fermion propagator S − , weseparate the BCS self-energy and obtain( S − ) G , G ′ = ( S − ) G , G ′ − ( Σ ∆ ) G , G ′ . (57)Here the inverse of the free fermion propagator reads( S − ) G , G ′ = ( S − + ) G , G ′
00 ( S − − ) G , G ′ ! (58)with the matrix elements given by( S − ± ) G , G ′ = ( i ω n + δµ ∓ ξ k + G ) δ G , G ′ . (59)The self-energy term Σ ∆ can be expressed as( Σ ∆ ) G , G ′ = ∆ F G , G ′ , (60)where the matrix F is defined as F G , G ′ = F + ) G , G ′ ( F − ) G , G ′ ! . (61) Here the blocks F ± defined in the G -space are given by( F + ) G , G ′ = f G − G ′ , ( F − ) G , G ′ = f ∗ G ′ − G . (62)Using the derivative expansion, we obtain Ω = Ω N + Pg ∆ + TV X ω n X k ∈ BZ ∞ X l = ∆ l l Tr h ( S F ) l i , (63)where Ω N is the free energy of the normal state, Ω N = − TV X ω n X k ∈ BZ Trln ( S − ) G , G ′ T . (64)Using the fact that S − is diagonal and k + G generates allcontinuum momenta p , we recover the usual expression (31).Completing the trace in the NG space, we find that the tracevanishes for odd l . Then the GL free energy can be expressedas Ω GL = α ∆ + ∞ X k = α k k ∆ k . (65)The second-order GL coe ffi cient is given by α = Pg + TV X ω n X k ∈ BZ Tr( S + F + S − F − ) . (66)The higher-order GL coe ffi cients α k ( k ≥
2) read α k = TV X ω n X k ∈ BZ Tr h ( S + F + S − F − ) k i . (67)Note that the trace is now taken only in the G -space. At zerotemperature and in the large volume limit, we obtain α = Pg + Z ∞−∞ dE π Z BZ d k (2 π ) A ( iE , k ) . (68)and α k = Z ∞−∞ dE π Z BZ d k (2 π ) A k ( iE , k ) . (69)for k ≥
2. Here the quantity A k ( iE , k ) is defined as A k ( iE , k ) = Tr n [ S + ( iE , k ) F + S − ( iE , k ) F − ] k o (70)with ( S ± ) G , G ′ ( iE , k ) = δ G , G ′ iE + δµ ∓ ξ k + G . (71)The GL coe ffi cients should be independent of the represen-tation. Therefore, α / P is still universal for all crystal struc-ture. We can evaluate it from the FF state. We have α = P α FF2 . (72)For the FF state we have f ( r ) = f ( z ) = e iqz , which can beregarded as a 1D crystal structure with periodicity a = π/ q .We have the Fourier transformations f ( z ) = ∞ X n = −∞ f n e i π na z , f ∗ ( z ) = ∞ X n = −∞ f ∗ n e − i π na z , (73)where the Fourier components f n and f ∗ n are given by f n = f ∗ n = δ n , . (74)Therefore, the matrices F ± read( F + ) n , n ′ = δ n − n ′ , , ( F − ) n , n ′ = δ n ′ − n , . (75)The GL coe ffi cient α FF2 can be expressed as α FF2 = g + Z ∞−∞ dE π Z ∞ k ⊥ dk ⊥ π Z q − q dk z π A ( iE , k ⊥ , k z ) . (76)The trace in A can be worked out analytically. We have A ( iE , k ⊥ , k z ) = Tr [ S + F + S − F − ] = X n X n , n , n δ n , n iE + δµ − ξ n δ n − n , δ n , n iE + δµ + ξ n δ n − n , = X n iE + δµ − ξ n iE + δµ + ξ n − , (77)where ξ n is defined as ξ n ( k ⊥ , k z ) = h k ⊥ + ( k z + nq ) i ν − µ. (78)Here ν = / ν = k z + nq generates all continuousmomenta in the z direction, the final result for α can be ex-pressed as (36).For higher-order GL coe ffi cients α k ( k ≥ A k ( iE , k ) becomes tedious. However, be-cause the G -space is discrete, this procedure can be real-ized by using a computer with a proper code. Since the G -space has infinite dimensions, we first make a truncation − D ≤ l , m , n ≤ D to perform the matrix operations. Pre-cise results can be approached by using a su ffi ciently large D .Some of the GL coe ffi cients may su ff er from UV divergenceand a proper regularization scheme should be implemented.For ultra-relativistic dispersion ε ( p ) = | p | , α and α are di-vergent. For non-relativistic case with ε ( p ) = | p | , only α isdivergent. The GL free energy can be expressed as Ω GL ( ∆ ) N δµ = P ¯ α ∆ δµ ! + ∞ X k = ¯ α k k ∆ δµ ! k , (79)where the GL coe ffi cients become dimensionless. We have¯ α k = δµ k − N α k . (80) These dimensionless GL coe ffi cients are universal in the weakcoupling limit. Therefore, we can evaluate them by usingthe non-relativistic dispersion because all the coe ffi cients α k ( k ≥
2) are free from UV divergence. On the other hand, inthis matrix formalism, we need to treat the momenta k and G separately. Therefore, the usual momentum cuto ff scheme − Λ < ξ p < Λ is not proper to achieve the weak coupling limit.The dependence on the chemical potential µ becomes explicitin the matrix formalism. Weak coupling corresponds to thecase µ ≃ ε F ≫ δµ . In practice, we can vary the value of µ/δµ and approach the weak coupling limit. IV. FORMALISM FOR SOME CRYSTAL STRUCTURES
In the final part of this paper, we consider some crystalstructures of which the order parameters ∆ ( r ) are real. Wehave f ∗ G = f − G and hence( F + ) G , G ′ = ( F − ) G , G ′ ≡ F G , G ′ = f G − G ′ . (81)In the following, we list the explicit forms of the matrix F , thepropagators S ± , and the GL coe ffi cients α k ( k ≥ ε ( p ) = p . In practice, this leads to faster convergenceand hence is better for numerical calculations. A. LO
The LO state is a superposition of two antipodal planewaves ( P =
2) withˆ n = (0 , , , ˆ n = (0 , , − . (82)The function f ( r ) can be expressed as f ( r ) = f ( z ) = qz ) , (83)which forms a 1D crystal structure with periodicity a = π/ q .The Fourier decomposition is given by f ( z ) = ∞ X n = −∞ f n e niqz , (84)where the Fourier component f n reads f n = δ n , + δ n , − . (85)The matrix form of the BdG equation is given by X n ′ H n , n ′ ( k ) φ n ′ ( k ) = E λ ( k ) φ n ( k ) , (86)where the Hamiltonian matrix H n , n ′ ( k ) reads ( ξ n − δµ ) δ n , n ′ ∆ f n − n ′ ∆ f n − n ′ ( − ξ n − δµ ) δ n , n ′ ! (87)with ξ n ( k ⊥ , k z ) = k ⊥ + ( k z + nq ) − µ. (88)The BZ can be defined as k z ∈ [ − q , q ]. The GL coe ffi cients α k ( k ≥
2) can be expressed as α k = Z ∞−∞ dE π Z ∞ k ⊥ dk ⊥ π Z q − q dk z π Tr [ n ] h ( S + FS − F ) k i . (89)The matrix elements of S ± and F are given by( S ± ) n , n ′ = δ n , n ′ iE + δµ ∓ ξ n , F n , n ′ = δ n , n ′ + + δ n , n ′ − . (90) B. Square
The Square state is a superposition of four plane waves ( P =
4) with ˆ n = (1 , , , ˆ n = ( − , , , ˆ n = (0 , , , ˆ n = (0 , − , . (91)The function f ( r ) can be expressed as f ( r ) = f ( x , y ) = (cid:2) cos(2 qx ) + cos(2 qy ) (cid:3) , (92)which forms a 2D crystal structure. The unit cell is generatedby two linearly independent vectors a = a e x and a = a e y with the lattice spacing a = π/ q . The Fourier decompositionis given by f ( x , y ) = ∞ X l , m = −∞ f lm e i π a ( lx + my ) , (93)where the Fourier component f lm reads f lm = ( δ l , + δ l , − ) δ m , + δ l , ( δ m , + δ m , − ) . (94)The matrix form of the BdG equation is given by X l ′ m ′ H [ lm ] , [ l ′ m ′ ] ( k ) φ [ l ′ m ′ ] ( k ) = E λ ( k ) φ [ lm ] ( k ) , (95)where the Hamiltonian matrix H [ lm ] , [ l ′ m ′ ] ( k ) reads ( ξ lm − δµ ) δ l , l ′ δ m , m ′ ∆ f l − l ′ , m − m ′ ∆ f l − l ′ , m − m ′ ( − ξ lm − δµ ) δ l , l ′ δ m , m ′ ! (96)with ξ lm ( k x , k y , k z ) = ( k x + lq ) + ( k y + mq ) + k z − µ. (97)The BZ can be defined as k x , k y ∈ [ − q , q ]. The GL coe ffi cients α k ( k ≥
2) are given by α k = Z ∞−∞ dE π Z q − q dk x π Z q − q dk y π Z ∞−∞ dk z π Tr [ lm ] h ( S + FS − F ) k i . (98)The matrix elements of S ± and F read( S ± ) [ lm ] , [ l ′ m ′ ] = δ l , l ′ δ m , m ′ iE + δµ ∓ ξ lm , F [ lm ] , [ l ′ m ′ ] = ( δ l , l ′ + + δ l , l ′ − ) δ m , m ′ + δ l , l ′ ( δ m , m ′ + + δ m , m ′ − ) . (99) C. BCC
The BCC state is a superposition of six plane waves ( P = n = (1 , , , ˆ n = ( − , , , ˆ n = (0 , , , ˆ n = (0 , − , , ˆ n = (0 , , , ˆ n = (0 , , − , (100)The function f ( r ) can be expressed as f ( r ) = (cid:2) cos(2 qx ) + cos(2 qy ) + cos(2 qz ) (cid:3) , (101)which forms a 3D crystal structure. The unit cell is generatedby three linearly independent vectors a = a e x , a = a e y ,and a = a e z with the lattice spacing a = π/ q . The Fourierdecomposition is given by f ( r ) = ∞ X l , m , n = −∞ f lmn e iq ( lx + my + nz ) , (102)where the Fourier component f lmn reads f lmn = (cid:0) δ l , + δ l , − (cid:1) δ m , δ n , + δ l , (cid:0) δ m , + δ m , − (cid:1) δ n , + δ l , δ m , (cid:0) δ n , + δ n , − (cid:1) . (103)The matrix form of the BdG equation is given by X l ′ m ′ n ′ H [ lmn ] , [ l ′ m ′ n ′ ] ( k ) φ [ l ′ m ′ n ′ ] ( k ) = E λ ( k ) φ [ lmn ] ( k ) , (104)where the Hamiltonian matrix H [ lmn ] , [ l ′ m ′ n ′ ] ( k ) reads ( ξ lmn − δµ ) δ l , l ′ δ m , m ′ δ n , n ′ ∆ f l − l ′ , m − m ′ , n − n ′ ∆ f l − l ′ , m − m ′ , n − n ′ ( − ξ lm − δµ ) δ l , l ′ δ m , m ′ δ n , n ′ ! (105)with ξ lmn = ( k x + lq ) + ( k y + mq ) + ( k z + nq ) − µ. (106)The BZ can be defined as k x , k y , k z ∈ [ − q , q ]. The GL coe ffi -cients α k ( k ≥
2) are given by α k = Z ∞−∞ dE π Z q − q dk x π Z q − q dk y π Z q − q dk z π Tr [ lmn ] h ( S + FS − F ) k i . (107)The matrix elements of S ± and F read( S ± ) [ lmn ] , [ l ′ m ′ n ′ ] = δ l , l ′ δ m , m ′ δ n , n ′ iE + δµ ∓ ξ lmn , F [ lmn ] , [ l ′ m ′ n ′ ] = ( δ l , l ′ + + δ l , l ′ − ) δ m , m ′ δ n , n ′ + δ l , l ′ ( δ m , m ′ + + δ m , m ′ − ) δ n , n ′ + δ l , l ′ δ m , m ′ ( δ n , n ′ + + δ n , n ′ − ) . (108)0 D. FCC
The FCC state is a superposition of eight plane waves ( P =
8) with ˆ n = √ , , , ˆ n = √ − , − , − , ˆ n = √ , − , , ˆ n = √ − , , − , ˆ n = √ , − , − , ˆ n = √ − , , , ˆ n = √ , , − , ˆ n = √ − , − , . (109)The function f ( r ) can be expressed as f ( r ) = qx √ ! cos qy √ ! cos qz √ ! , (110)which forms a 3D crystal structure. The unit cell is generatedby three linearly independent vectors a = a e x , a = a e y , and a = a e z with the lattice spacing a = √ π/ q . The Fourierdecomposition is given by f ( r ) = ∞ X l , m , n = −∞ f lmn e i √ q ( lx + my + nz ) , (111)where the Fourier component f lmn reads f lmn = (cid:0) δ l , + δ l , − (cid:1) (cid:0) δ m , + δ m , − (cid:1) (cid:0) δ n , + δ n , − (cid:1) . (112)The matrix form of the BdG equation and the Hamiltonianmatrix H [ lmn ] , [ l ′ m ′ n ′ ] ( k ) take the same forms as (104) and (105)but with ξ lmn given by ξ lmn = k x + lq √ ! + k y + mq √ ! + k z + nq √ ! − µ. (113)The BZ can be defined as k x , k y , k z ∈ [ − q √ , q √ ]. The GLcoe ffi cients α k ( k ≥
2) are given by α k = Z ∞−∞ dE π Z q √ − q √ dk x π Z q √ − q √ dk y π Z q √ − q √ dk z π Tr [ lmn ] h ( S + FS − F ) k i . (114) The matrix elements of S ± and F can be expressed as( S ± ) [ lmn ] , [ l ′ m ′ n ′ ] = δ l , l ′ δ m , m ′ δ n , n ′ iE + δµ ∓ ξ lmn , F [ lmn ] , [ l ′ m ′ n ′ ] = ( δ l , l ′ + + δ l , l ′ − )( δ m , m ′ + + δ m , m ′ − ) × ( δ n , n ′ + + δ n , n ′ − ) . (115) V. SUMMARY
In this work, we presented a new derivation of the GL freeenergy of crystalline (color) superconductors and compared itwith the conventional diagrammatic derivation [10]. The di-agrammatic derivation and our derivation can be attributed tothe use of two di ff erent representations of the fermion field:The diagrammatic derivation employs the usual continuousmomentum representation, while our derivation uses the dis-crete Bloch representation. In either formalism, we need toevaluate the trace of the form Tr h ( S + F + S − F − ) k i . In the usualmomentum representation, taking this trace leads to a summa-tion of all possible configurations that satisfy the momentumconstraint (7). For large k , the number of these configurationsbecomes also large, which makes the calculation tedious. Inour formalism, this trace is taken in the discrete Bloch space.Therefore, the calculation of the GL coe ffi cients can be com-puterized: One can generate the matrices F ± and S ± and per-form the matrix operations by using a computer with a propercode. Once the computerization is realized, we may be able toevaluate the GL free energy to arbitrary order in ∆ . With theinformation of the higher-order terms in the GL free energy,we can give more reliable predictions for the phase transitions. Acknowledgments — The work of G. C. was supported bythe NSFC under Grant No. 11335005 and the MOST underGrant Nos. 2013CB922000 and 2014CB845400. The workof L. H. was supported by the US Department of Energy Top-ical Collaboration “Neutrinos and Nucleosynthesis in Hot andDense Matter”. L. H. also acknowledges the support fromFrankfurt Institute for Advanced Studies in the early stage ofthis work. [1] A. I. Larkin and Y. N. Ovchinnikov, Zh. Eksp. Teor. Fiz. ,1136 (1964).[2] P. Fulde and R. A. Ferrell, Phys. Rev. , A550 (1964).[3] R. Casalbuoni and G. Nardulli, Rev. Mod. Phys. , 263 (2004).[4] B. S. Chandrasekhar, Appl. Phys. Lett. , 7 (1962); A. M.Clogston, Phys. Rev. Lett. , 266 (1962).[5] K. Gloos, et. al. , Phys. Rev. Lett. , 501 (1993); R. Modler, et.al. , ibid , 1292 (1996); A. Bianchi, et. al. , ibid , 187004(2003). [6] J. L. O’Brien, et. al. , Phys. Rev. B61 , 1584 (2000).[7] L. Balicas, et. al. , Phys. Rev. Lett. , 067002 (2001); M. A.Tanatar, et. al. , Phys. Rev. B66 , 134503 (2002).[8] S. Kasahara, et. al. , Proc. Natl. Acad. Sci. , 16309 (2014).[9] M. Alford, J. Bowers, and K. Rajagopal, Phys. Rev.
D63 ,074016 (2001).[10] J. A. Bowers, K. Rajagopal, Phys. Rev.
D66 , 065002 (2002).[11] I. Shovkovy and M. Huang, Phys. Lett.
B564 , 205 (2003).[12] M. Alford, C. Kouvaris, and K. Rajagopal, Phys. Rev. Lett. , D70 , 054004 (2004).[14] M. Huang and I. Shovkovy, Phys. Rev.
D70 , 051501(R) (2004).[15] R. Casalbuoni, R. Gatto, M. Mannarelli, G. Nardulli, and M.Ruggieri, Phys. Lett.
B605 , 362 (2005).[16] K. Fukushima, Phys. Rev.
D72 , 074002 (2005).[17] I. Giannakis and H.-C. Ren, Phys. Lett.
B611 , 137 (2005).[18] E. V. Gorbar, M. Hashimoto, and V. A. Miransky, Phys. Rev.Lett. , 022005 (2006).[19] M. Mannarelli, K. Rajagopal, and R. Sharma, Phys. Rev. D73 ,114012 (2006); K. Rajagopal and R. Sharma, Phys. Rev.
D74 ,094019 (2006).[20] L. He, M. Jin, and P. Zhang, Phys. Rev.
D75 , 036003 (2007); A.Sedrakian and D. H. Rischke, Phys. Rev.
D80 , 074022 (2009);X.-G. Huang and A. Sedrakian, Phys. Rev.
D82 , 045029 (2010).[21] R. Anglani, R. Casalbuoni, M. Ciminale, R. Gatto, N. Ippolito,M. Mannarelli, and M. Ruggieri, Rev. Mod. Phys. , 509(2014).[22] M. Alford, K. Rajagopal, and F. Wilczek, Phys. Lett. B422 , 247(1998). [23] R. Rapp, T. Schaefer, E. V. Shuryak, and M. Velkovsky, Phys.Rev. Lett. , 53 (1998).[24] M. Alford, K. Rajagopal, and F. Wilczek, Nucl. Phys. B537 ,443 (1999).[25] D. T. Son, Phys. Rev.
D59 , 094019 (1999).[26] M. Alford, K. Rajagopal, T. Schaefer, and A. Schmitt, Rev.Mod. Phys. , 1455 (2008); D. H. Rischke, Prog. Part. Nucl.Phys. , 197 (2004); M. Buballa, Phys. Rep. , 205 (2005);I. A. Shovkovy, Found. Phys. , 1309 (2005).[27] M. W. Zwierlein, et. al. , Science , 492 (2006); G. B. Par-tridge, et. al. , ibid , 503 (2006).[28] D. E. Sheehy and L. Radzihovsky, Phys. Rev. Lett. , 060401(2006); Ann. Phys. (N. Y.) , 1790(2007).[29] A. I. Buzdin and H. Kachkachi, Phys. Lett. A225 , 341 (1997).[30] R. Combescot and C. Mora, Eur. Phys. J.
B28 , 397 (2002).[31] L. Jiang and J. Ye, Phys. Rev.
B76 , 184104 (2007).[32] R. Combescot and C. Mora, EPL , 79 (2004).[33] C. Mora and R. Combescot, Phys. Rev. B 71 , 214504 (2005).[34] D. Nickel and M. Buballa, Phys. Rev.
D79 , 054009 (2009).[35] C. Cao, L. He, and P. Zhuang, Phys. Rev.