Phase diagrams for two-component D-wave fermion superfluids with density imbalance

Finite-temperature phase structures for asymmetric two-component D-wave fermion superfluids in the weak coupling BCS regime are investigated. The Fulde–Ferrell–Larkin-Ovchinnikov (FFLO) state becomes non-degenerate with the orientation of Cooper pair momentum due to the anisotropic pairing gap for D-wave pairing. The new parameter θ 0, corresponding to the orientation of Cooper pair momentum, greatly enriches the phase structures. Especially, the phase separation (PS) region in the phase diagram is significantly reduced comparing to that for S-wave pairing. Theoretically, the stability condition of solutions is equivalent to the positive definite Hessian matrix of free energy in the hyperplane defined by the number density constraints in the parameter space spanned by the average chemical potential μ, relative chemical potential h, paring gap Δ, pair momentum q, and θ 0. In addition, the stability condition of the superfluid density is essentially included in the positive definite Hessian matrix. Moreover, the positiveness of imbalance number susceptibility Dδρ/Dh is always equivalent to the positiveness of the smaller eigenvalue λ − of the number susceptibility matrix. While the criteria ∂2Ω/∂Δ2 and ∂2Ω∂Δ2∂2Ω∂q2−∂2Ω∂Δ∂q2 against PS in previous studies are equivalent to λ − only under the prerequisite that the solution is a local minimum in the corresponding hyperplane.

In the exotic phenomena, stemming from the mismatched Fermi surface, the non-zero Cooper pair momentum corresponding to the anisotropic Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states proposed by Fulde and Ferrell [15] and Larkin and Ovchinnikov [16] in an S-wave superconductor with Zeeman splitting about sixty years ago is especially concerned. While due to the sensitivity to impurities [35,36] and orbital pair breaking [37], convincing evidence for experimental realization of FFLO state at low temperatures and high fields is still rather scare. However, the appearance has been suggested in the unconventional heavy fermion superconductors [35] and organic superconductors [38,39] as well as

Formalism
In the present paper, the physical system we interested in is an infinite system composed of two components of number density ρ a and ρ b and of chemical potential μ a and μ b , respectively. They interact via an anisotropic D-wave attractive interaction in three-dimensional free space. The full fermion Green's function is given by − μ a , and − μ b represents the single-particle energy of the two species of Fermions, respectively. ω m is odd Matsubara frequencies [49]. 2q corresponds to the nonzero center-of-mass momentum, i.e., the Cooper pair momentum in FFLO state. The off diagonal element, corresponding to the anomalous propagator, takes the form of F † = −Δ(k) (iω m +δε−ξ Δ )(iω m +δε+ξ Δ ) . The summation over the Matsubara frequencies provides the density matrix of particles in the condensate, i.e., the pairing probability, Accordingly, the gap equation is expressed as where E ± k = ξ Δ ± δε is the quasi-particle energy in the condensate and V(k, k ) represents the anisotropic point contact interaction. f(x) = 1 + exp( x T ) −1 is the well-known Fermi-Dirac distribution function. The fermion Green's function gives the number densities of the components with the density distributions For fixed number densities, it is essential that the coupled equations (2) and (3) should be solved self-consistently.

Angle-dependent gap and the orientation of Cooper pair momentum
To investigate D-wave pairing conveniently, we simply consider the point contact pairing interaction of the form withk = k/k. We take the z axis in the direction of the quantization axis of the spherical harmonic function Y m 2 . Different from the S-wave pairing, the pairing gap, which is also anisotropic, can be expressed as Generally, all components δ m should be taken into account in solving the gap equations. However, the thermodynamical stability investigation of D-wave pairing base on the Ginzburg-Landau energy functional indicates [50] that the stable D-wave state corresponds to cyclic D-wave phase or single component Y 2 2 pairing phase in the weak coupling limit when the coupling constant g |m| are degenerate with m. For asymmetric case, the other single component (such as Y 0 2 ) may as well turn into stable at certain asymmetries in weak coupling region. In addition, because of the crystal lattice, the pairing interaction of the heavyfermion superconductors [51] could be different with various m. Under such circumstances, the most energetically favored component may exclude the others. Therefore, it make sense to consider the pairing of single component. For simplicity, we suppose g m=0 g m=other , i.e., we only take m = 0 component pairing into account in the present work corresponding to Y 0 2 = 5/16π(3 cos 2 θ − 1). Which can ensure the axial symmetry of ξ Δ and Δ(k) along the z axis (the symmetry axis of ADG). Substituting equations (5) and (6) in to gap equation (2), one can obtain the equation for δ m , Noting that the pairing gap δ 0 is a real number, we replace it by real parameter Δ for convenience hereafter. The pairing gap turn into Δ(k) = 5/4Δ(3 cos 2 θ − 1), which remains the rotational symmetry along the axis (θ = 0, φ = 0). Once the pairing gap is anisotropic, the quasi-particle energy becomes non-degenerate with the orientation of Cooper pair momentum. To account this effect, the symmetric and asymmetric part of the spectrum can be expressed Where the average and relative chemical potentials are defined as μ = (μ a + μ b )/2 and h = (μ a − μ b )/2, respectively. The orientations of 2q and k are related to (θ 0 , φ 0 ) and (θ, φ) in spherical coordinate, respectively. One should note that the properties of spherical harmonics have been employed to expand the cosine of the angle between 2q and k. It is obvious that the term k·q m due to the Cooper pair momentum can reduce the suppression from the mismatched Fermi surface h in certain direction. We stress here that the gap equation (7) and the densities (3) are independent of the constant φ 0 , i.e., φ 0 can be eliminated by choosing a special spherical coordinate in which the orientation of 2q is represented as (θ 0 , φ 0 = 0). Only θ 0 as a parameter describes the direction of Cooper pair momentum. Different to the S-wave pairing, the FFLO states are non-degenerate for the orientation of Cooper pair momentum. And both the value q and the direction of Cooper pair momentum θ 0 should be determined by minimizing the energy of the system. These two values are also embodied in the gap equation (7) via the symmetric and asymmetric part of the spectrum.

Thermodynamics
For asymmetric two-component fermion systems with fixed the total number densities ρ = ρ a + ρ b and the asymmetry α = (ρ a − ρ b )/ρ at temperature T, the thermodynamic quantity describing the system is the free energy defined as where δρ = ρ a − ρ b . The thermodynamic potential Ω takes the form of Where the entropy can be expressed as in the mean-filed approximation and the internal energy of the grand canonical ensemble reads The chemical potential of quasiparticles is also embodied in the first term of the internal energy. The second term of the internal energy represents the BCS mean-field interaction energy among the particles in condensate. Using the particle pairing probability (1), density distribution (4), and gap equation (7), the thermodynamic potential can be expressed in the following form (noting the property f(ω) ln f(ω) The gap equation (7) and density (3) can be equivalently expressed as If the Cooper pair momentum q = 0, one easily finds that dF/dΔ = ∂Ω/∂Δ. As will be discussed in section 2.3 the gap equation is equivalent to the statement that the free energy reaches a minimum.

Solutions and stability conditions of the pairing solution for fixed ρ and δρ
For a fixed particle number system, the solution is equivalent to the minimum of F (μ, h, Δ, q, θ 0 ) in the hyperplane P(Δ, q, θ 0 ) defined by the two number densities equation in the parameter space spanned by μ, h, Δ, q, θ 0 . The parameters μ and h can be treated as implicit functions of Δ, q, and θ 0 as defined by the density constraints (14). In hyperplane P(Δ, q, θ 0 ), the free parameters Δ, q, and θ 0 are determined by minimizing the free energy. In order to make the expression more clear, we define where x i = Δ, q, θ 0 , and β c = μ, h, Δ, q, θ 0 . The derivatives ∂μ ∂x i and ∂h ∂x i can be obtained from the implicit relations, i.e., density constraints (3). Whereas the From the relation (9) between the free energy F and the thermodynamic potential Ω, it is easy to check that DF/Dx i = ∂Ω/∂x i , which develops into dF /dΔ = ∂Ω/∂Δ when q and θ 0 are absent. Accordingly, the solutions of DF Dx i = 0, corresponding to the minimum of free energy, turn into with the density constraints The stability of these solutions corresponds to positive definite Hessian matrix H of the free energy. The elements of H is here x i, j = Δ, q, θ 0 , and β b,c = μ, h, Δ, q, θ 0 . The positive definite 3 × 3 H requires all the determinants of the sequential principal minor of H are positive, which imply that det M t > 0, det M q > 0, det M Δ > 0, and det M h > 0 (see the appendix A for detail). det M h is always positive due to the number density constraints. For the other three constraints, only one key constraint remains to play crucial role in determining the stability of the solution under the circumstances considered, which will be discussed in section 3. As an example, if the solution is confined to q = 0, the stability condition of solutions [10] corresponds to the element of Hessian matrix ∂Δ∂h dh dΔ , which is always positive when the solution of the gap equation with fixed particle number exists for S-wave pairing. However, it is not true for D-wave pairing as will shown in section 3.1.

Superfluid density
For a superfluid state to be stable, the superfluid density must be positive [11]. Using the linear response theory, the superfluid density can be calculated via the response of the system to an external vector potential. Since the thermodynamics of the system with fixed number densities is the free energy, the supercurrent j s and the superfluid density ρ s should be defined via the Taylor expansion of F near the solution q s = (q s , θ 0s , φ 0s ) (here supposing μ s , h s , Δ s , q s , θ 0s represent certain solution of equations (16) and (17)), responding to a small velocity v s Where The equation (16) can ensure the zero of the supercurrent j s . The positive superfluid density corresponds to the positive definite matrix D 2 F Dq 2 , which is also called as the Cooper pair momentum susceptibility [12]. One should note that in the previous works [10][11][12]34] for S-wave pairing, the Cooper pair momentum susceptibility mainly focus on the appearance of the Cooper pair momentum, i.e., the linear response of the system to the value of q deviating from 0. In the present work, the Cooper pair momentum susceptibility refers in particular to the same focus point. In the spherical coordinates, D 2 F Dq 2 can be expressed as Here all the derivative with respect to φ 0 vanishes due to the fact that the free energy is independent of azimuth angle φ 0 , indicating the superfluid density matrix (including theê φ 0ê φ 0 component) is independent of φ 0 as well. Under the circumstance q = 0 for uniform superfluid state, Due to the axial symmetry of the pairing gap we consider in the present work, the superfluid density tensor is no longer isotropic and can be decomposed into a transverse and a longitudinal part due to the O(2) symmetry. The superfluid density tensor can be expressed as in the spherical coordinates with where It should be stressed that the forms 1 Thus the superfluid density is also isotropic. Due to the angle dependence of the pairing gap for D-wave pairing, the discrepancy between ρ T and ρ L result in the anisotropic ρ s .
For the uniform superfluid state with q = 0, one essentially confines the gap equation ∂Ω ∂Δ = 0 and the density equation (17) to the hyperplane P(Δ) defined by the density constraint in the parameter space spanned by μ, h, Δ (the free parameter is Δ only). The stable solution corresponds to that the free energy reaches a minimum in hyperplane P(Δ). And the stability condition of the solution is D 2 F DΔ 2 = d 2 F dΔ 2 > 0. However, the stable solution in hyperplane P(Δ) may become unstable in hyperplane P(Δ, q, θ 0 ) (the free parameters are Δ, q, θ 0 ), namely, the stability condition d 2 F dΔ 2 > 0 is insufficient once the solution (Δ, q = 0, θ 0 ) is treated as a specific solution in hyperplane P(Δ, q, θ 0 ). Under this circumstances, the stability condition of the specific solution should be the positive definiteness of the 3 × 3 H, which develop into Where D 2 F DΔ 2 > 0 corresponds to the stability condition of solutions [10] in hyperplane P(Δ). The extra requirement D 2 F Dq 2 | q=0 = ∂ 2 Ω ∂q 2 > 0 determines the appearance of the pairing momentum. If ρ L > ρ T , it is obvious that the ∂ 2 Ω ∂q 2 = ρ T firstly reaches zero and changes sign along the orientation of θ 0 = π 2 resulting in the first appearance of the FFLO state with q perpendicular to the symmetry axis of ADG. And the positive definiteness of the superfluid density result in the identical conclusion as ∂ 2 Ω ∂q 2 > 0 when ρ L > ρ T . This is also true for the case ρ L < ρ T . From the perspective of the occurrence of the FFLO state, D 2 F Dq 2 | q=0 > 0 is equivalent to the positive definiteness of the superfluid density. As will discussed in section 3.2, D 2 F Dθ 2 0 corresponds to theê θ 0ê θ 0 component of the superfluid density and describes the stability of the orientation of Cooper pair momentum.

Number susceptibility and stability of homogeneous two-component superfluid against phase separation
Generally, the stability condition for a two-component system against changes in the densities of its component can be described by the total free energy of the system [12,52] Supposing a tiny changes δρ i (r) deviate from uniform density distribution, due to the number conversation the first-order variation δF automatically vanishes. And the second-order variation δ 2 F is expressed as The stable homogeneous density distribution requires the 2 × 2 matrix ∂ 2 F /∂ρ σ ∂ρ τ = ∂μ τ /∂ρ σ to be positively definite. As shown in appendix B, the matrix ∂μ τ /∂ρ σ is the inverse matrix of the number With the relationships 2 is a unitary transformation matrix. Therefore the stable homogeneous density distribution corresponds to the positive definite matrix Which is equivalent to the statement that the eigenvalues λ ± of M 2×2 are both positive. Here with the trace Tr(M 2×2 ) = Dρ/Dμ + Dδρ/Dh. Fortunately, λ + is always positive which roughly corresponds to Dρ Dμ . Another eigenvalue λ − roughly corresponds to Dδρ Dh . The numerics shows λ − changes sign when Dδρ Dh changes sign for S-wave pairing [10]. Our calculations for D-wave pairing also confirm this. Actually, the authors [11,12] adopt χ = Dδρ Dh > 0 as the criterion for the stability of superfluid against the PS.
Since the appearance of det M [see equation (B3) for the definition of M] in the denominator of the expressions in equation (B7), each element of M 2×2 will change sign when det M approaches zero and changes sign. Under this circumstance, 1/ det M 2×2 approaches zero and changes sign as discussed in reference [10]. Note that 1/ det M 2×2 = 1/(λ + λ − ) and λ + is always positive, the stability condition that the eigenvalues of number susceptibility matrix are positive is equivalent to κ(Δ, q, In reference [10], the numerics confirm that the stability condition κ(Δ) = ∂ 2 Ω ∂Δ 2 > 0 is equivalent to the imbalance number susceptibility χ = Dδρ Dh > 0. However, as will be shown in section 3.4, it is true only when the solutions are the local minima in the considered parameter space.

Results and discussions
As is well known, the point contact interaction adopted in the current work will result in the ultraviolet divergence of the gap equation. To avoid this, we employ a hard cutoff k Λ for three dimensional momentum which is usually used in the study of color superconductivity at moderate baryon density and isospin asymmetric nuclear matter. In this scheme, the normalization of the coupling constantg 0 = g 0 4π/(k F m) and the cutoff Λ = (k Λ /k F ) 2 with k F = (3π 2 ρ) 1/3 are also used to avoid the specific choice of model parameters. The value of Δ depend both on the cutoff Λ and the coupling constantg 0 , however, the scaled quantities such as Δ/Δ 0 and h/Δ 0 (hereafter Δ 0 is the value of Δ in the symmetric case with h = 0) are regularization scheme independent in the weak coupling region [18,34,53,54]. In the present work, we take the weak coupling constantg 0 = 0.5 with the cutoff Λ = 4 as in the previous work [48,55].

ADG state
We first consider the case without the pair momentum, i.e., q = 0, and we call this ADG state. The stability condition of the nontrivial solution is det M Δ / det M h = d 2 F /dΔ 2 > 0, which is always true when the solution of the gap equation with fixed number density exists for S-wave pairing. However, for D-wave pairing, there exist two solutions at low temperature, corresponding to a saddle point and a local minimum . One should note that this unstable solution is quite analogous to the Sarma state [17]. Both this unstable state in the present paper and Sarma state correspond to the saddle point of the thermodynamic quantities. Different as Sarma state which appears as the saddle point of the thermodynamic potential, the unstable state here corresponds to the saddle point of the free energy. Besides, since there exist gapless excitations in momentum space for D-wave pairing, the criteria of Sarma state for S-wave are difficult to be applied in D-wave pairing.
The existence of the unstable solution at low temperature indicate the first order transition from the superconduction state to the normal state since Δ does not approach to zero when α → α c (α c represent the critical asymmetry where the δf become positive). This is also confirmed from the free energy curve. It is obvious that the free energy of ADG state cross the free energy of normal state as exhibited in the insert at T/T F = 0.01. However, the first transition only occur at low temperature, with increasing temperature, the domain of the unstable state is getting smaller, and disappear when T/T F > 0.0125. The insert at T/T F = 0.02 implies a second order transition from the superconducting state to normal state.
The stability condition against PS develop into κ(Δ) = ∂ 2 Ω/∂Δ 2 when q = 0. Our numerics for D-wave pairing indicate that ∂ 2 Ω/∂Δ 2 > 0 is equivalent to χ = Dδρ/Dh > 0 only if the solution is locally minimal in hyperplane P(Δ). Since the ground state corresponds to the locally minimal solution only, both ∂ 2 Ω/∂Δ 2 and χ can be adopted to distinguish the PS. The finite-temperature phase diagram for D-wave pairing without Cooper pair momentum are displayed in figure 2, which is significantly different from the phase diagram for S-wave pairing [10,12]. There are two narrow yellow regions corresponding to the PS, one is settled around α = 0.09 at very low temperature. Another is located near the boundary of α c . This difference from S-wave pairing may result from the existence of the nodes and zero-lines of the pairing gap which can provide phase-space for the excess particles. Under this circumstance, it might be unnecessary to  destroy the pairing state to settle the excess particle. Therefore, the superconducting state in certain region which is unstable against PS for S-wave pairing becomes stable against PS for D-wave pairing.
To investigate pair momentum susceptibility, the scaled (scaled by 3mρ/8) transverse and longitudinal part of superfluid density are exhibited in figure 2 as well, corresponding to the red short dashed and black solid thin lines, respectively. On the whole, the transverse part ρ T always reaches zero first, which indicates the FFLO state with pair momentum perpendicular to the symmetry axis of ADG appears earlier than that along the symmetry axis of ADG when the asymmetry α grows. Nevertheless, with increasing temperature, ρ L reaches zero earlier than ρ T when α get large, indicating a phase transition from the uniform superconducting state to the FFLO state with Cooper pair momentum along θ 0 = 0. Figure 3 shows ρ L reaches zero firstly at temperature T/T F = 0.0326. However, both ρ T and ρ L are positive for all the superconducting state when T/T F > 0.0332, indicating the vanishing of the FFLO state as shown in figure 7.

FFLO state with Cooper pair momentum perpendicular/parallel to the symmetry axis of ADG
Due to the anisotropic pairing gap, the FFLO state is non-degenerate with the orientation of Cooper pair momentum. Two specific direction [47], corresponding to θ 0 = 0 and θ 0 = π 2 , are firstly discussed for neutron-proton pairing in asymmetric nuclear matter. Actually, one can easily find that θ 0 = 0 ( π 2 ) satisfies the equation of θ 0 which corresponds to that the direction of q and the symmetry axis of ADG is parallel (orthogonal). And we call these two configurations, i.e., θ 0 = 0 and θ 0 = π 2 , the FFLO-ADG-P state and the FFLO-ADG-O state, respectively.
For the case θ 0 = 0 ( π 2 ) with nonzero Cooper pair momentum, ∂ 2 Ω ∂β c ∂θ 0 | β c =θ 0 = 0, the Hessian matrix turn into The positive definiteness corresponds to det M q > 0, det M Δ > 0, det M h > 0, and ∂ 2 Ω individually describes the stability of the orientation of the Cooper pair momentum. In addition to these constraints, the denominator in the number susceptibility (B7) ∂Δ∂q ) 2 which corresponds to the stability condition against PS. It is obvious the constraint from PS is coincident with the FFLO state for S-wave pairing in reference [10], i.e., In fact, on account of the vanishing of ∂ 2 Ω ∂β c ∂θ 0 β c =θ 0 , the hyperplane P(Δ, q, θ 0 = 0/ π 2 ) constitutes an individual subspace. If one confine the solution in hyperplane P(Δ, q, θ 0 = 0/ π 2 ), the results are coincident with that in the hyperplane P(Δ, q) defined by the number density constraints in the parameter space spanned by μ, h, Δ, q for S-wave pairing. Therefore, both the stability conditions of solutions and against PS are identical to that for S-wave pairing.
However, as discussed in section 2.4, the stability conditions are insufficient if the solutions are obtained in the subspace. And the superfluid density need to be supplemented to constrain the stability. Especially, the superfluid density in FFLO-ADG-P states reads where cos θ 0 q 2 sin θ 0 ∂Ω ∂θ 0 has been calculated by adopting the L'Hospital rule. Note that the positive definiteness of   hyperplane P(Δ, q, θ 0 ). For convenience, unless otherwise specified, the FFLO-ADG-P state refers to the FFLO-ADG-P-2 state hereafter.
We now discuss the phase transition from uniform ADG state to FFLO state. Since the transverse part of the superfluid density ρ T in ADG state reaches zero first in most case, a phase transition from the ADG state to FFLO-ADG-O state arises first. We call this phase transition the A-O phase transition for convenience in the present work. It is obvious that the FFLO-ADG-O-1 state is related to the emergence of FFLO state. Hence, the transverse part of the superfluid density ρ T (scaled by 3mρ/8) in ADG state, the scaled Cooper pair momentum q/k F in FFLO state, and the scaled difference of the free energy between superconducting and normal state δf = (F| s − F| n )/E 0 for the ADG (blue short dotted lines), FFLO-ADG-O-1 (black solid lines), and FFLO-ADG-O-2 (red short dash-dotted lines) states at different temperatures are exhibited in the upper, middle, and lower panel of figure 5, respectively. At low temperature T/T F = 0.001, ρ T reaches zero at 0.0643 and change sign when α > 0.0643 indicating that the solution with q = 0 is no longer locally minimal in hyperplane P(Δ, q, θ 0 ) when α > 0.0643. However, in a small region 0.0616 < α < 0.0643 there are two local minima. One corresponds to the ADG state while the other corresponds to the Except the two special directions of the Cooper pair momentum, the orientation between the two special direction may also survive as discussed in reference [48]. However, it might not be true for the    figure 7. But, however, the FFLO-ADG-P state reappears in a small temperature region 0.0318 < T/T F < 0.332 as displayed in the insert of figure 7. This can also be confirmed in right panel of figure 6, the orientation susceptibility for FFLO-ADG-P becomes positive when α > 0.153. The corresponding B-P phase transition at this small temperature region become the second order due to the vanishing of the FFLO-ADG-B-3 state at high temperature.
In addition, the stability against PS can be determined by the κ(Δ, q, θ 0 ) = det M > 0 or χ = Dδρ/Dh > 0. Together with the constraints from the ADG, FFLO-ADG-O, and FFLO-ADG-P state, the T − α phase diagram is finally exhibited in figure 7. For convenience, the name of FFLO-ADG-O, FFLO-ADG-B, and FFLO-ADG-P have been shortened to ADG-O, ADG-B, and ADG-P in figure 7. We stress here the A-P phase transition (phase transition from the ADG state to the FFLO-ADG-P state) emerges near the critical temperature where the FFLO states disappear due to the fact that ρ L < ρ T occurs at high temperature as displayed in figure 3. The A-P phase transition is of second order because the unstable FFLO-ADG-P-1 solution disappears at high temperature. In the T − α phase diagram, two narrow PS regions are located near the occurrence region of the FFLO-ADG-O state and the FFLO-ADG-B state, respectively. These PS regions are obviously reduced comparing to that for ADG phase diagram. Moreover, this T − α phase diagram significantly differentiates from that for S-wave pairing [10,12].

Comparing the stability constraints κ and χ against phase separation
The stability of the polarized superfluid against PS requires the number susceptibility matrix must be positive definite. As discussed in section 2.5, both χ = Dδρ Dh > 0 and κ(Δ, q, θ 0 ) = det M can be adopted as the criteria of the stability of the polarized superfluid against PS, and they are equivalent to each other for S-wave pairing [10]. For D-wave paring, these two stability constraints might be different from that for S-wave pairing owing to the emergence of the unstable states, such as FFLO-ADG-O-1 and FFLO-ADG-O-3 states.
For the sake of verifying the validity of κ and χ, the scaled χ, κ, and the second eigenvalue λ − of number susceptibility matrix M 2×2 for FFLO state with θ 0 = π 2 are plotted in figure 8. The denominator in the number susceptibility (B7) evolves into κ(Δ, q) , κ is opposite to λ − and κ cannot provide available constraint against the PS. However, one may probably evaluate the stability condition against the PS after confirming the stability of the solution in the adopted parameter space. Under the prerequisite that the solution should at least be a local minimum of F in the adopted parameter space, the criterion κ is equivalent to the criteria χ and λ − . The solutions in the S-wave pairing with fixed number density is always stable, and the two criteria are equivalent to each other [10]. Nevertheless one should treat the D-wave pairing with carefulness.

Summary and outlook
In summary, we have studied phase structure of a two-component D-wave fermion superfluid with density imbalance in the present work. The pairing gap becomes no longer isotropic owing to the anisotropic D-wave pairing interaction. As a result, the FFLO state develops into non-degenerate with the orientation of Cooper pair momentum. The extra freedom, i.e. the direction of Cooper pair momentum θ 0 , can significantly enrich the phase structures. Briefly, the FFLO state for S-wave pairing splits into three branches according to the possible orientations of Cooper pair momentum, i.e., the FFLO-ADG-O, FFLO-ADG-P, and FFLO-ADG-B states. The phase diagram significantly differentiate from that of S-wave pairing. To investigate the phase structure, the stability constraints are theoretically extended as well for D-wave pairing. The main conclusions are as follows.
(a) For a fixed density system, the stable solutions are essentially correspond to the local minima of the free energy F in the hyperplane defined by the density constraints in the parameter space. The solution is related to DF/Dx i with x i representing the free parameters such as Δ, q, θ 0 . The stability condition corresponds to that the Hessian matrix of the free energy H ij = D 2 F /Dx i Dx j is positive definite. (b) The superfluid density consists of two parts, i.e., the response to the value of Cooper pair momentum and the orientation of Cooper pair momentum, due to the anisotropic pairing gap for D-wave pairing.
These two parts correspond to the Cooper pair momentum susceptibility and orientation susceptibility. The constraints from superfluid density are effective only if the parameter space is sufficient. In fact, the stability condition of the super fluid density is essentially included in the positive definite Hessian matrix as long as the parameter space is complete. (c) The stability condition χ = Dδρ/Dh against PS is equivalent to the positiveness of eigenvalue λ − of the number susceptibility matrix. While the criteria κ(Δ) 2 , and κ(Δ, q, θ 0 ) = det M are equivalent to λ − only under the prerequisite that the Hessian matrix is positive definite in the adopted parameter space. (d) At low temperature, the transition from the uniform ADG state to normal state is of first order and becomes the second order with increasing temperature. For FFLO state, the A-O, O-B, B-P phase transitions are of the first order at low temperature, and become the second order transition with increasing temperature as well. The transition from the FFLO state to normal state is of second order. Near the critical temperature where FFLO state vanishes, there exists a second order A-P phase transition. These diversified phase transitions might be distinguishing feature of the considered D-wave pairing. The specific heat and thermal expansion may offer a circumstantial evidence of the FFLO state for D-wave pairing via analyzing the phase transition type. Additionally, the node structures of these ADG-FFLO states in momentum space might be distinctive. Therefore, the scanning tunneling microscopy quasiparticle interference [56][57][58][59] can perhaps be adopted to investigate the pairing structure of ADG-FFLO states. (e) Both the phase diagrams with and without Cooper pair momentum in the T − α plane show there are two narrow PS regions which are significantly different from those for S-wave pairing. This difference of the S-wave pairing may result from the existence of the nodes and zero-line of the quasi-particle spectrum near Fermi surface which can provide phase-space for the excess particles. Accordingly, it might be unnecessary to destroy the pairing state to settle the excess particles. And the superconducting state in certain region which is unstable against PS for S-wave pairing turns into stable against PS for D-wave pairing.
The current calculations are carried in the framework of mean-field level, the pseudogap effect [13] should also be take into account at finite temperature to obtain the realistic phase diagram in future. In addition, we adopt only the m = 0 component pairing interaction of D-wave pairing in the current calculation, it remains for future work that whether the other components result in the same phase diagram and whether the single component pairing is stable. the 4 × 4 matrix M q as and the 2 × 2 matrix M h as for convenience. Owing to the property , there are 15 independent elements of matrix M t , which read kq m (cos θ 0 sin θ cos φ − sin θ 0 cos θ) kq m (cos θ 0 sin θ cos φ − sin θ 0 cos θ)