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

To investigate D-wave pairing conveniently, we simply consider the point contact pairing interaction of the form

$$\begin{equation}V(\mathbf{k},{\mathbf{k}}^{\,\prime })=\sum\limits _{m=-2}^{m=2}4\pi {g}_{\vert m\vert }{Y}_{2}^{m}(\hat{\mathbf{k}}){Y}_{2}^{m\ast }({\hat{\mathbf{k}}}^{\,\prime }),\end{equation} \tag{ 5 }$$

with $\hat{\mathbf{k}}=\mathbf{k}/k$. We take the z axis in the direction of the quantization axis of the spherical harmonic function ${Y}_{2}^{m}$. Different from the S-wave pairing, the pairing gap, which is also anisotropic, can be expressed as

$$\begin{equation}{\Delta}(\hat{\mathbf{k}})=\sum\limits _{m=-2}^{m=2}\sqrt{4\pi }{\delta }_{\text{m}}{Y}_{2}^{m}(\hat{\mathbf{k}}).\end{equation} \tag{ 6 }$$

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}_{2}^{0}$) 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}_{2}^{0}=\sqrt{5/16\pi }(3\,{\cos }^{2}\,\theta -1)$. Which can ensure the axial symmetry of ξ_Δ and ${\Delta}(\hat{\mathbf{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 ,

$$\begin{equation}{\delta }_{0}={g}_{0}\sum\limits _{{\mathbf{k}}^{\,\prime }}\frac{{\delta }_{0}4\pi \vert {Y}_{2}^{0}({\hat{\mathbf{k}}}^{\,\prime }){\vert }^{2}}{2\sqrt{{\xi }_{{\mathbf{k}}^{\,\prime }}^{2}+4\pi {\delta }_{0}^{2}\vert {Y}_{2}^{0}({\hat{\mathbf{k}}}^{\,\prime }){\vert }^{2}}}\left[1-f({E}_{{\mathbf{k}}^{\,\prime }}^{+})-f({E}_{{\mathbf{k}}^{\,\prime }}^{-})\right].\end{equation} \tag{ 7 }$$

Noting that the pairing gap δ₀ is a real number, we replace it by real parameter Δ for convenience hereafter. The pairing gap turn into ${\Delta}(\mathbf{k})=\sqrt{5/4}{\Delta}(3\,{\cos }^{2}\,\theta -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

$$\begin{align}\hfill {\xi }_{\mathbf{k}}& =\frac{{\mathbf{k}}^{2}}{2m}+\frac{{\mathbf{q}}^{2}}{2m}-\mu ,\hfill \\ \hfill \delta \varepsilon & =h-\frac{\mathbf{k}\cdot \mathbf{q}}{m}=h-\frac{kq}{m}\cos (\hat{\mathbf{k}\mathbf{q}})\hfill \\ \hfill & =h-\frac{kq}{m}\left[\mathrm{sin}\,{\theta }_{0}\,\mathrm{sin}\,\theta \cos (\phi -{\phi }_{0})+\cos \,{\theta }_{0}\,\cos \,\theta \right].\hfill \end{align} \tag{ 8 }$$

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 (θ₀, ϕ₀) 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 $\frac{\mathbf{k}\cdot \mathbf{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 ϕ₀, i.e., ϕ₀ can be eliminated by choosing a special spherical coordinate in which the orientation of 2q is represented as (θ₀, ϕ₀ = 0). Only θ₀ 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 θ₀ 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.

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

$$\begin{equation}\mathcal{F}={\Omega}+\mu \rho +h\delta \rho ,\end{equation} \tag{ 9 }$$

where δρ = ρ_a − ρ_b . The thermodynamic potential Ω takes the form of

$$\begin{equation}{\Omega}=U-T\mathit{\text{S}}.\end{equation} \tag{ 10 }$$

Where the entropy can be expressed as

$$\begin{equation}\mathit{\text{S}}=-\sum\limits _{\mathbf{k}}\left\{\,f({E}_{\mathbf{k}}^{+})\mathrm{ln}\,f({E}_{\mathbf{k}}^{+})+[1-f({E}_{\mathbf{k}}^{+})]\mathrm{ln}[1-f({E}_{\mathbf{k}}^{+})]+f({E}_{\mathbf{k}}^{-})\mathrm{ln}\,f({E}_{\mathbf{k}}^{-})+[1-f({E}_{\mathbf{k}}^{-})]\mathrm{ln}[1-f({E}_{\mathbf{k}}^{-})]\right\}\end{equation} \tag{ 11 }$$

in the mean-filed approximation and the internal energy of the grand canonical ensemble reads

$$\begin{equation}U=\sum\limits _{\mathbf{k}}\left[{\varepsilon }_{a}{n}_{a}(\mathbf{k})+{\varepsilon }_{b}{n}_{b}(\mathbf{k})\right]-\sum\limits _{\mathbf{k},{\mathbf{k}}^{\,\prime }}V(\mathbf{k},{\mathbf{k}}^{\,\prime }){\nu }^{{\dagger}}(\mathbf{k})\nu ({\mathbf{k}}^{\,\prime }).\end{equation} \tag{ 12 }$$

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(\omega )\mathrm{ln}\,f(\omega )+[1-f(\omega )]\mathrm{ln}[1-f(\omega )]=-\frac{\omega }{T}-\mathrm{ln}(1+{\text{e}}^{\frac{-\omega }{T}})$)

$$\begin{equation}{\Omega}=\frac{{{\Delta}}^{2}}{{g}_{0}}-\sum\limits _{\mathbf{k}}\left[{\xi }_{{\Delta}}-{\xi }_{\mathbf{k}}+T\,\mathrm{ln}(1+{\text{e}}^{-{E}_{\mathbf{k}}^{+}/T})+T\,\mathrm{ln}(1+{\text{e}}^{-{E}_{\mathbf{k}}^{-}/T})\right].\end{equation} \tag{ 13 }$$

The gap equation (7) and density (3) can be equivalently expressed as

$$\begin{equation}\frac{\partial {\Omega}}{\partial {\Delta}}=0,\quad \rho =-\frac{\partial {\Omega}}{\partial \mu },\quad \delta \rho =-\frac{\partial {\Omega}}{\partial h}.\end{equation} \tag{ 14 }$$

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.

For a fixed particle number system, the solution is equivalent to the minimum of $\mathcal{F}(\mu ,h,{\Delta},q,{\theta }_{0})$ in the hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0})$ defined by the two number densities equation in the parameter space spanned by μ, h, Δ, q, θ₀. The parameters μ and h can be treated as implicit functions of Δ, q, and θ₀ as defined by the density constraints (14). In hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0})$, the free parameters Δ, q, and θ₀ are determined by minimizing the free energy. In order to make the expression more clear, we define

$$\begin{equation}\frac{\mathcal{D}\mathcal{F}}{\mathcal{D}{x}_{i}}\equiv \sum\limits _{c}\frac{\partial \mathcal{F}}{\partial {\beta }_{c}}\frac{\partial {\beta }_{c}}{\partial {x}_{i}},\end{equation} \tag{ 15 }$$

where x_i = Δ, q, θ₀, and β_c = μ, h, Δ, q, θ₀. The derivatives $\frac{\partial \mu }{\partial {x}_{i}}$ and $\frac{\partial h}{\partial {x}_{i}}$ can be obtained from the implicit relations, i.e., density constraints (3). Whereas the $\frac{\partial {x}_{j}}{\partial {x}_{i}}={\delta }_{j}^{i}$. From the relation (9) between the free energy $\mathcal{F}$ and the thermodynamic potential Ω, it is easy to check that $\mathcal{D}\mathcal{F}/\mathcal{D}{x}_{i}=\partial {\Omega}/\partial {x}_{i}$, which develops into $\mathrm{d}\mathcal{F}/\mathrm{d}{\Delta}=\partial {\Omega}/\partial {\Delta}$ when q and θ₀ are absent. Accordingly, the solutions of $\frac{\mathcal{D}\mathcal{F}}{\mathcal{D}{x}_{i}}=0$, corresponding to the minimum of free energy, turn into

$$\begin{equation}\frac{\partial {\Omega}}{\partial {\Delta}}=0,\quad \frac{\partial {\Omega}}{\partial q}=0,\quad \frac{\partial {\Omega}}{\partial {\theta }_{0}}=0,\end{equation} \tag{ 16 }$$

with the density constraints

$$\begin{equation}\rho =-\frac{\partial {\Omega}}{\partial \mu },\quad \delta \rho =-\frac{\partial {\Omega}}{\partial h}.\end{equation} \tag{ 17 }$$

The stability of these solutions corresponds to positive definite Hessian matrix $\mathcal{H}$ of the free energy. The elements of $\mathcal{H}$ is

$$\begin{equation}{\mathcal{H}}_{ij}=\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{x}_{i}\mathcal{D}{x}_{i}}=\sum\limits _{bc}\frac{{\partial }^{2}\mathcal{F}}{\partial {\beta }_{b}\partial {\beta }_{c}}\frac{\partial {\beta }_{b}}{\partial {x}_{i}}\frac{\partial {\beta }_{c}}{\partial {x}_{j}},\end{equation} \tag{ 18 }$$

here x_i,j = Δ, q, θ₀, and β_b,c = μ, h, Δ, q, θ₀. The positive definite $3\times 3\ \mathcal{H}$ requires all the determinants of the sequential principal minor of $\mathcal{H}$ are positive, which imply that $\mathrm{det}\,{\mathcal{M}}_{t} > 0$, $\mathrm{det}\,{\mathcal{M}}_{q} > 0$, $\mathrm{det}\,{\mathcal{M}}_{{\Delta}} > 0$, and $\mathrm{det}\,{\mathcal{M}}_{h} > 0$ (see the appendix A for detail). $\mathrm{det}\,{\mathcal{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 ${\mathcal{H}}_{11}=\frac{{d}^{2}\mathcal{F}}{d{{\Delta}}^{2}}=\frac{{\partial }^{2}{\Omega}}{\partial {{\Delta}}^{2}}+\frac{{\partial }^{2}{\Omega}}{\partial {\Delta}\partial \mu }\frac{\mathrm{d}\mu }{\mathrm{d}{\Delta}}+\frac{{\partial }^{2}{\Omega}}{\partial {\Delta}\partial h}\frac{dh}{d{\Delta}}$, 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.

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 $\mathcal{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

$$\begin{equation}\mathcal{F}({\mathbf{v}}_{\text{s}})=\mathcal{F}(0)+{\mathbf{j}}_{\text{s}}\cdot {\mathbf{v}}_{\text{s}}+\frac{1}{2}{\rho }_{\text{s}}{\mathbf{v}}_{\text{s}}^{2}+\cdots \end{equation} \tag{ 19 }$$

Where

$$\begin{equation}\mathcal{F}(0)=\mathcal{F}{\vert }_{{\mathbf{q}}_{\text{s}}},\quad {\left.{\mathbf{j}}_{\text{s}}=m\frac{\mathcal{D}\mathcal{F}}{\mathcal{D}\mathbf{q}}\right\vert }_{{\mathbf{q}}_{\text{s}}},\quad {\left.{\rho }_{\text{s}}={m}^{2}\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\mathbf{q}}^{2}}\right\vert }_{{\mathbf{q}}_{\text{s}}}.\end{equation} \tag{ 20 }$$

The equation (16) can ensure the zero of the supercurrent j_s. The positive superfluid density corresponds to the positive definite matrix $\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\mathbf{q}}^{2}}$, which is also called as the Cooper pair momentum susceptibility [12]. One should note that in the previous works [10–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, $\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\mathbf{q}}^{2}}$ can be expressed as

$$\begin{equation}\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\mathbf{q}}^{2}}=\left(\begin{matrix}\hfill \ \ \ \ \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{q}^{2}}\hfill & \hfill \ \ \frac{\mathcal{D}}{\mathcal{D}q}\left(\frac{1}{q}\frac{\mathcal{D}\mathcal{F}}{\mathcal{D}{\theta }_{0}}\right)\hfill & \hfill \quad \quad 0\hfill \\ \hfill \frac{\mathcal{D}}{\mathcal{D}q}\left(\frac{1}{q}\frac{\mathcal{D}\mathcal{F}}{\mathcal{D}{\theta }_{0}}\right)\hfill & \hfill \frac{1}{q}\frac{\mathcal{D}\mathcal{F}}{\mathcal{D}q}+\frac{1}{{q}^{2}}\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\theta }_{0}^{2}}\hfill & \hfill \quad \quad 0\hfill \\ \hfill \quad 0\hfill & \hfill \quad \ \ 0\hfill & \hfill \frac{1}{q}\frac{\mathcal{D}\mathcal{F}}{\mathcal{D}q}+\frac{\cos \,{\theta }_{0}}{{q}^{2}\,\mathrm{sin}\,{\theta }_{0}}\frac{\mathcal{D}\mathcal{F}}{\mathcal{D}{\theta }_{0}}\hfill \end{matrix}\right).\end{equation} \tag{ 21 }$$

Here all the derivative with respect to ϕ₀ vanishes due to the fact that the free energy is independent of azimuth angle ϕ₀, indicating the superfluid density matrix (including the ${\hat{\mathbf{e}}}_{{\phi }_{0}}{\hat{\mathbf{e}}}_{{\phi }_{0}}$ component) is independent of ϕ₀ as well.

Under the circumstance q = 0 for uniform superfluid state, $\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\mathbf{q}}^{2}}{\vert }_{q=0}=\frac{{\partial }^{2}{\Omega}}{\partial {\mathbf{q}}^{2}}$, the superfluid density developed into

$$\begin{equation}{\rho }_{\text{s}}={m}^{2}\frac{{\partial }^{2}{\Omega}}{\partial {\mathbf{q}}^{2}}{\vert }_{{\mathbf{q}}_{\text{s}}=0}.\end{equation} \tag{ 22 }$$

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

$$\begin{equation}{\rho }_{\text{s}}=\left(\begin{matrix}\hfill {\rho }_{\text{T}}\,{\mathrm{sin}}^{2}\,{\theta }_{0}+{\rho }_{\text{L}}\,{\cos }^{2}\,{\theta }_{0}\hfill & \hfill ({\rho }_{\text{T}}-{\rho }_{\text{L}})\cos \,{\theta }_{0}\,\mathrm{sin}\,{\theta }_{0}\hfill & \hfill 0\hfill \\ \hfill ({\rho }_{\text{T}}-{\rho }_{\text{L}})\cos \,{\theta }_{0}\,\mathrm{sin}\,{\theta }_{0}\hfill & \hfill {\rho }_{\text{T}}\,{\cos }^{2}\,{\theta }_{0}+{\rho }_{\text{L}}\,{\mathrm{sin}}^{2}\,{\theta }_{0}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\rho }_{\text{T}}\hfill \end{matrix}\right)\end{equation} \tag{ 23 }$$

in the spherical coordinates with

$$\begin{align}\hfill {\rho }_{\text{T}}=m\rho +\int \frac{{k}^{4}dk}{8{\pi }^{2}}{\int }_{-1}^{1}\mathrm{d}\,\cos \,\theta \,{\mathrm{sin}}^{2}\,\theta \left[{f}^{\,\prime }({\tilde{E}}_{\mathbf{k}}^{+})+{f}^{\,\prime }({\tilde{E}}_{\mathbf{k}}^{-})\right]\\ \hfill {\rho }_{\text{L}}=m\rho +\int \frac{{k}^{4}\,\mathrm{d}k}{4{\pi }^{2}}{\int }_{-1}^{1}\,\mathrm{d}\,\cos \,\theta \,{\cos }^{2}\,\theta \left[{f}^{\,\prime }({\tilde{E}}_{\mathbf{k}}^{+})+{f}^{\,\prime }({\tilde{E}}_{\mathbf{k}}^{-})\right],\end{align} \tag{ 24 }$$

where ${\tilde{E}}_{\mathbf{k}}^{\pm }={E}_{\mathbf{k}}^{\pm }{\vert }_{q=0}$. It should be stressed that the forms $\frac{1}{{q}^{2}}\frac{\partial \mathcal{F}}{\partial {\theta }_{0}}$ and $\frac{1}{q}\frac{\partial \mathcal{F}}{\partial q}$ might be non-zero though the solutions require $\frac{\partial \mathcal{F}}{\partial q}=\frac{\partial \mathcal{F}}{\partial {\theta }_{0}}=0$ when q = 0. For S-wave pairing, Δ(k) is independent of θ, ρ_T = ρ_L. 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 $\frac{\partial {\Omega}}{\partial {\Delta}}=0$ and the density equation (17) to the hyperplane $\mathcal{P}({\Delta})$ 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 $\mathcal{P}({\Delta})$. And the stability condition of the solution is $\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{{\Delta}}^{2}}=\frac{{\mathrm{d}}^{2}\mathcal{F}}{\mathrm{d}{{\Delta}}^{2}} > 0$. However, the stable solution in hyperplane $\mathcal{P}({\Delta})$ may become unstable in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0})$ (the free parameters are Δ, q, θ₀), namely, the stability condition $\frac{{\mathrm{d}}^{2}\mathcal{F}}{\mathrm{d}{{\Delta}}^{2}} > 0$ is insufficient once the solution (Δ, q = 0, θ₀) is treated as a specific solution in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0})$. Under this circumstances, the stability condition of the specific solution should be the positive definiteness of the $3\times 3\ \mathcal{H}$, which develop into

$$\begin{equation}\mathcal{H}=\left(\begin{matrix}\hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{{\Delta}}^{2}}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{q}^{2}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\theta }_{0}^{2}}\hfill \end{matrix}\right).\end{equation} \tag{ 25 }$$

Where $\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{{\Delta}}^{2}} > 0$ corresponds to the stability condition of solutions [10] in hyperplane $\mathcal{P}({\Delta})$. The extra requirement $\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{q}^{2}}{\vert }_{q=0}=\frac{{\partial }^{2}{\Omega}}{\partial {q}^{2}} > 0$ determines the appearance of the pairing momentum. If ρ_L > ρ_T, it is obvious that the $\frac{{\partial }^{2}{\Omega}}{\partial {q}^{2}}={\rho }_{\text{T}}$ firstly reaches zero and changes sign along the orientation of ${\theta }_{0}=\frac{\pi }{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 $\frac{{\partial }^{2}{\Omega}}{\partial {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, $\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{q}^{2}}{\vert }_{q=0} > 0$ is equivalent to the positive definiteness of the superfluid density. As will discussed in section 3.2, $\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\theta }_{0}^{2}}$ corresponds to the ${\hat{\mathbf{e}}}_{{\theta }_{0}}{\hat{\mathbf{e}}}_{{\theta }_{0}}$ component of the superfluid density and describes the stability of the orientation of Cooper pair momentum.

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], $F=\int {\mathrm{d}}^{3}\mathbf{r}\,\mathcal{F}[{\rho }_{\tau }(\mathbf{r})]$ (τ = a, b). 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 δ² F is expressed as

$$\begin{equation}{\delta }^{2}F=\frac{1}{2}\int {\mathrm{d}}^{3}\mathbf{r}\sum\limits _{\sigma ,\tau =a,b}\frac{{\partial }^{2}\mathcal{F}}{\partial {\rho }_{\sigma }\partial {\rho }_{\tau }}\delta {\rho }_{\sigma }(\mathbf{r})\delta {\rho }_{\tau }(\mathbf{r}).\end{equation} \tag{ 26 }$$

The stable homogeneous density distribution requires the 2 × 2 matrix ${\partial }^{2}\mathcal{F}/\partial {\rho }_{\sigma }\partial {\rho }_{\tau }=\partial {\mu }_{\tau }/\partial {\rho }_{\sigma }$ to be positively definite. As shown in appendix B, the matrix ∂μ_τ /∂ρ_σ is the inverse matrix of the number susceptibility matrix $\mathcal{D}{\rho }_{\sigma }/\mathcal{D}{\mu }_{\tau }$, where we define

$$\begin{align}\hfill \frac{\mathcal{D}}{\mathcal{D}{\mu }_{a}}& \equiv \frac{\partial }{\partial {\mu }_{a}}+\sum\limits _{i}\frac{\partial {x}_{i}}{\partial {\mu }_{a}}\frac{\partial }{\partial {x}_{i}}\hfill \\ \hfill \frac{\mathcal{D}}{\mathcal{D}{\mu }_{b}}& \equiv \frac{\partial }{\partial {\mu }_{b}}+\sum\limits _{i}\frac{\partial {x}_{i}}{\partial {\mu }_{b}}\frac{\partial }{\partial {x}_{i}}.\hfill \end{align} \tag{ 27 }$$

With the relationships $\mathcal{D}/\mathcal{D}\mu =\mathcal{D}/\mathcal{D}{\mu }_{a}+\mathcal{D}/\mathcal{D}{\mu }_{b}$ and $\mathcal{D}/\mathcal{D}h=\mathcal{D}/\mathcal{D}{\mu }_{a}-\mathcal{D}/\mathcal{D}{\mu }_{b}$. One can find easily that

$$\begin{equation}\left(\begin{matrix}\hfill \frac{\mathcal{D}{\rho }_{a}}{\mathcal{D}{\mu }_{a}}\hfill & \hfill \frac{\mathcal{D}{\rho }_{a}}{\mathcal{D}{\mu }_{b}}\hfill \\ \hfill \frac{\mathcal{D}{\rho }_{b}}{\mathcal{D}{\mu }_{a}}\hfill & \hfill \frac{\mathcal{D}{\rho }_{b}}{\mathcal{D}{\mu }_{b}}\hfill \end{matrix}\right)=\frac{1}{2}A\left(\begin{matrix}\hfill \frac{\mathcal{D}\rho }{\mathcal{D}\mu }\hfill & \hfill \frac{\mathcal{D}\rho }{\mathcal{D}h}\hfill \\ \hfill \frac{\mathcal{D}\delta \rho }{\mathcal{D}\mu }\hfill & \hfill \frac{\mathcal{D}\delta \rho }{\mathcal{D}h}\hfill \end{matrix}\right)A,\end{equation} \tag{ 28 }$$

here $A={A}^{-1}=\left(\begin{matrix}\hfill 1\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill -1\hfill \end{matrix}\right)/\sqrt{2}$ is a unitary transformation matrix. Therefore the stable homogeneous density distribution corresponds to the positive definite matrix

$$\begin{equation}{\mathcal{M}}_{2\times 2}=\left(\begin{matrix}\hfill \frac{\mathcal{D}\rho }{\mathcal{D}\mu }\hfill & \hfill \frac{\mathcal{D}\rho }{\mathcal{D}h}\hfill \\ \hfill \frac{\mathcal{D}\delta \rho }{\mathcal{D}\mu }\hfill & \hfill \frac{\mathcal{D}\delta \rho }{\mathcal{D}h}\hfill \end{matrix}\right).\end{equation} \tag{ 29 }$$

Which is equivalent to the statement that the eigenvalues λ_± of ${\mathcal{M}}_{2\times 2}$ are both positive. Here

$$\begin{equation}{\lambda }_{\pm }=\frac{\mathrm{Tr}({\mathcal{M}}_{2\times 2})\pm \sqrt{\mathrm{Tr}{({\mathcal{M}}_{2\times 2})}^{2}-4\,\mathrm{det}({\mathcal{M}}_{2\times 2})}}{2}\end{equation} \tag{ 30 }$$

with the trace $\mathrm{Tr}({\mathcal{M}}_{2\times 2})=\mathcal{D}\rho /\mathcal{D}\mu +\mathcal{D}\delta \rho /\mathcal{D}h$. Fortunately, λ₊ is always positive which roughly corresponds to $\frac{\mathcal{D}\rho }{\mathcal{D}\mu }$. Another eigenvalue λ₋ roughly corresponds to $\frac{\mathcal{D}\delta \rho }{\mathcal{D}h}$. The numerics shows λ₋ changes sign when $\frac{\mathcal{D}\delta \rho }{\mathcal{D}h}$ changes sign for S-wave pairing [10]. Our calculations for D-wave pairing also confirm this. Actually, the authors [11, 12] adopt $\chi =\frac{\mathcal{D}\delta \rho }{\mathcal{D}h} > 0$ as the criterion for the stability of superfluid against the PS.

Since the appearance of $\mathrm{det}\,\mathcal{M}$ [see equation (B3) for the definition of $\mathcal{M}$] in the denominator of the expressions in equation (B7), each element of ${\mathcal{M}}_{2\times 2}$ will change sign when $\mathrm{det}\,\mathcal{M}$ approaches zero and changes sign. Under this circumstance, $1/\mathrm{det}\,{\mathcal{M}}_{2\times 2}$ approaches zero and changes sign as discussed in reference [10]. Note that $1/\mathrm{det}\,{\mathcal{M}}_{2\times 2}=1/({\lambda }_{+}{\lambda }_{-})$ and λ₊ is always positive, the stability condition that the eigenvalues of number susceptibility matrix are positive is equivalent to $\kappa ({\Delta},q,{\theta }_{0})=\mathrm{det}\,\mathcal{M} > 0$, which evolves into $\kappa ({\Delta})=\frac{{\partial }^{2}{\Omega}}{\partial {{\Delta}}^{2}} > 0$ when q = 0 and $\kappa ({\Delta},q)=\frac{{\partial }^{2}{\Omega}}{\partial {{\Delta}}^{2}}\frac{{\partial }^{2}{\Omega}}{\partial {q}^{2}}-{(\frac{{\partial }^{2}{\Omega}}{\partial {\Delta}\partial q})}^{2} > 0$ when q ≠ 0 for S-wave pairing. In reference [10], the numerics confirm that the stability condition $\kappa ({\Delta})=\frac{{\partial }^{2}{\Omega}}{\partial {{\Delta}}^{2}} > 0$ is equivalent to the imbalance number susceptibility $\chi =\frac{\mathcal{D}\delta \rho }{\mathcal{D}h} > 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.

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 $\mathrm{det}\,{\mathcal{M}}_{{\Delta}}/\mathrm{det}\,{\mathcal{M}}_{h}={\mathrm{d}}^{2}\mathcal{F}/\mathrm{d}{{\Delta}}^{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 of $\mathcal{F}$ in hyperplane $\mathcal{P}({\Delta})$, respectively. The locally minimal solution satisfy the stability condition whereas the saddle point solution correspond to the negative ${\mathrm{d}}^{2}\mathcal{F}/\mathrm{d}{{\Delta}}^{2}$ (as shown in figure 4(a)). The solution related to the local minimum (saddle point) is marked with red short dash doted (blue solid) line in figure 1. The normalized pairing gap Δ/Δ₀ and the scaled difference of the free energy between superconducting and normal states $\delta f=(\mathcal{F}{\vert }_{\text{s}}-\mathcal{F}{\vert }_{\text{n}})/{E}_{0}$ are shown in the upper and lower panel of figure 1, respectively. Where the constant ${E}_{0}={k}_{\text{F}}^{5}/(8{\pi }^{2}m)$ and the Fermi temperature ${T}_{\text{F}}={k}_{\text{F}}^{2}/(2m{k}_{\text{B}})$. 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.

Zoom In Zoom Out Reset image size

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 κ(Δ) = ∂²Ω/∂Δ² when q = 0. Our numerics for D-wave pairing indicate that ∂²Ω/∂Δ² > 0 is equivalent to $\chi =\mathcal{D}\delta \rho /\mathcal{D}h > 0$ only if the solution is locally minimal in hyperplane $\mathcal{P}({\Delta})$. Since the ground state corresponds to the locally minimal solution only, both ∂²Ω/∂Δ² 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.

Zoom In Zoom Out Reset image size

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. 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.

Zoom In Zoom Out Reset image size

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 and ${\theta }_{0}=\frac{\pi }{2}$, are firstly discussed for neutron-proton pairing in asymmetric nuclear matter. Actually, one can easily find that θ₀ = 0 $(\frac{\pi }{2})$ satisfies the equation of θ₀

$$\begin{equation}\frac{\partial {\Omega}}{\partial {\theta }_{0}}=\sum\limits _{\mathbf{k}}\frac{kq}{m}\left[f({E}_{\mathbf{k}}^{-})-f({E}_{\mathbf{k}}^{+})\right](\cos \,{\theta }_{0}\,\mathrm{sin}\,\theta \,\cos \,\phi -\,\mathrm{sin}\,{\theta }_{0}\,\cos \,\theta )=0,\end{equation} \tag{ 31 }$$

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 and ${\theta }_{0}=\frac{\pi }{2}$, the FFLO-ADG-P state and the FFLO-ADG-O state, respectively.

For the case θ₀ = 0 $(\frac{\pi }{2})$ with nonzero Cooper pair momentum, $\frac{{\partial }^{2}{\Omega}}{\partial {\beta }_{c}\partial {\theta }_{0}}{\vert }_{{\beta }_{c}\ne {\theta }_{0}}=0$, the Hessian matrix turn into

$$\begin{equation}\mathcal{H}=\left(\begin{matrix}\hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{{\Delta}}^{2}}\hfill & \hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\Delta}\mathcal{D}q}\hfill & \hfill 0\hfill \\ \hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}q\mathcal{D}{\Delta}}\hfill & \hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{q}^{2}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}\hfill \end{matrix}\right).\end{equation} \tag{ 32 }$$

The positive definiteness corresponds to $\mathrm{det}\,{\mathcal{M}}_{q} > 0$, $\mathrm{det}\,{\mathcal{M}}_{{\Delta}} > 0$, $\mathrm{det}\,{\mathcal{M}}_{h} > 0$, and $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}} > 0$. Our numerics show that $\mathrm{det}\,{\mathcal{M}}_{{\Delta}}$ and $\mathrm{det}\,{\mathcal{M}}_{h}$ are always positive, $\mathrm{det}\,{\mathcal{M}}_{q} > 0$ and $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}} > 0$ become the key stability condition. Where $\mathrm{det}\,{\mathcal{M}}_{q}=\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{{\Delta}}^{2}}\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{q}^{2}}-{(\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\Delta}\mathcal{D}q})}^{2}$ determines the stability of solutions in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0}=0/\frac{\pi }{2})$ and $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{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) evolves into $\frac{{\partial }^{2}{\Omega}}{\partial {{\Delta}}^{2}}\frac{{\partial }^{2}{\Omega}}{\partial {q}^{2}}-{(\frac{{\partial }^{2}{\Omega}}{\partial {\Delta}\partial 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., $\kappa ({\Delta},q)=\frac{{\partial }^{2}{\Omega}}{\partial {{\Delta}}^{2}}\frac{{\partial }^{2}{\Omega}}{\partial {q}^{2}}-{(\frac{{\partial }^{2}{\Omega}}{\partial {\Delta}\partial q})}^{2}$. In fact, on account of the vanishing of ${\left.\frac{{\partial }^{2}{\Omega}}{\partial {\beta }_{c}\partial {\theta }_{0}}\right\vert }_{{\beta }_{c}\ne {\theta }_{0}}$, the hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0}=0/\frac{\pi }{2})$ constitutes an individual subspace. If one confine the solution in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0}=0/\frac{\pi }{2})$, the results are coincident with that in the hyperplane $\mathcal{P}({\Delta},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

$$\begin{equation}{\rho }_{\text{s}}={m}^{2}\left(\begin{matrix}\hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{q}^{2}}\hfill & \hfill \ \ \ 0\hfill & \hfill \ \ \ 0\hfill \\ \hfill \ \ 0\hfill & \hfill \frac{1}{{q}^{2}}\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}\hfill & \hfill \ \ \ 0\hfill \\ \hfill \ \ 0\hfill & \hfill \ \ \ 0\hfill & \hfill \frac{1}{{q}^{2}}\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}\hfill \end{matrix}\right),\end{equation} \tag{ 33 }$$

where $\frac{\cos \,{\theta }_{0}}{{q}^{2}\,\mathrm{sin}\,{\theta }_{0}}\frac{\partial {\Omega}}{\partial {\theta }_{0}}$ has been calculated by adopting the L'Hospital rule. Note that the positive definiteness of the 2 × 2 matrix $\left(\begin{matrix}\hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{{\Delta}}^{2}}\hfill & \hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{\Delta}\mathcal{D}q}\hfill \\ \hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}q\mathcal{D}{\Delta}}\hfill & \hfill \frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{q}^{2}}\hfill \end{matrix}\right)$, corresponding to the stable solution in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0}=0)$, require the positiveness of $\frac{{\mathcal{D}}^{2}\mathcal{F}}{\mathcal{D}{q}^{2}}$. Thus, only the ${\hat{\mathbf{e}}}_{{\theta }_{0}}{\hat{\mathbf{e}}}_{{\theta }_{0}}$ component can provide an extra constraint $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}$ on orientation of the Cooper pair momentum, and we can call this the orientation susceptibility of the Cooper pair momentum similar as the Cooper pair momentum susceptibility. However, this constraint is essentially included in 3 × 3 Hessian matrix of the hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0})$. For FFLO-ADG-O state, the ${\hat{\mathbf{e}}}_{{\phi }_{0}}{\hat{\mathbf{e}}}_{{\phi }_{0}}$ component of the superfluid density becomes to zero due to the fact $\frac{\partial {\Omega}}{\partial {\theta }_{0}}=0$, the superfluid density might be positive semi-definite matrix. Nevertheless, the stability condition of the superfluid density can still provide the same constraint $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}$ as FFLO-ADG-P state. Summarily, the stability conditions of superfluid density are essentially included in the positive definite Hessian matrix if the adopted parameter space is complete.

In order to have an entire inspection of the solutions, the scaled gap Δ/Δ₀, the scaled (scaled by E₀) $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}$ and the dimensionless $\mathrm{det}\,{\mathcal{M}}_{q}$ (Ω has been normalized by E₀, μ, h, and Δ have been normalized by ${k}_{\text{F}}^{2}/(2m)$, and q has been normalized by k_F in the calculation) at T/T_F = 0.001 are displayed in figure 4. The FFLO-ADG-O and FFLO-ADG-P states are shown in figures 4(b) and (d), respectively. The red thick lines represent the scale pairing gap Δ/Δ₀, where as the black thin lines correspond to the key stability conditions in hyperplane $P({\Delta},q,{\theta }_{0}=0/\frac{\pi }{2})$ and the blue thin lines are related to the extra constraint $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}$. And the short dash-dotted and short dashed lines correspond to the stable solutions in hyperplane $P({\Delta},q,{\theta }_{0}=0/\frac{\pi }{2})$, while the solid and short dotted lines are related to the saddle point solutions. In addition, as discussed in section 3.1, there exists an unstable solution related to the saddle point for uniform ADG state in hyperplane $\mathcal{P}({\Delta})$, corresponding to the negative ${\mathrm{d}}^{2}\mathcal{F}/\mathrm{d}{{\Delta}}^{2}$ shown in figure 4(a).

Zoom In Zoom Out Reset image size

In figure 4(b), four branches of solution [48] for the FFLO-ADG-O state are displayed. Our numerics show that $\mathrm{det}\,{\mathcal{M}}_{h}$ and $\mathrm{det}\,{\mathcal{M}}_{{\Delta}}$ are always positive for these four branches of solution. Thus only $\mathrm{det}\,{\mathcal{M}}_{q} > 0$ remains to determine the stability of the solutions in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0}=\frac{\pi }{2})$. Form this stability condition, one can find that the FFLO-ADG-O-1 and FFLO-ADG-O-3 states actually correspond to the saddle point due to the negative $\mathrm{det}\,{\mathcal{M}}_{q}$. This is also confirmed in the free energy surface $\mathcal{F}({\Delta},q)$ (here μ, h are obtained from the density constraint). Instead, the FFLO-ADG-O-2 and FFLO-ADG-O-4 are both related to the local minimum of $\mathcal{F}$ in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0}=\frac{\pi }{2})$. And the free energy of FFLO-ADG-O-2 (FFLO-ADG-O-4) state get smaller than that of FFLO-ADG-O-4 (FFLO-ADG-O-2) state in the region 0.615 < α < 0.161 (α > 0.161). But the orientation susceptibility of the Cooper pair momentum $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}$ becomes negative when α > 0.1562, indicating that the FFLO-ADG-O-4 state is always unstable in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0})$. Only the FFLO-ADG-O-2 state could be the candidate of the ground state, thus unless otherwise specified, the FFLO-ADG-O state refers to the FFLO-ADG-O-2 state hereafter. For the FFLO-ADG-P state shown in figure 4(d), two solutions corresponding to the saddle point and local minimum are labelled as the FFLO-ADG-P-1 and FFLO-ADG-P-2, respectively. Owing to the key stability condition $\mathrm{det}\,{\mathcal{M}}_{q}$ of solutions in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0}=0)$ and the orientation susceptibility $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}$, only the FFLO-ADG-P-2 state becomes the local minimum and the possible ground state when α > 0.155 in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{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 $\delta f=(\mathcal{F}{\vert }_{\text{s}}-\mathcal{F}{\vert }_{\text{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 $\mathcal{P}({\Delta},q,{\theta }_{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 FFLO-ADG-O-2 state. Between these two local minima, there exists a saddle point solution, i.e., the FFLO-ADG-O-1 state. The A-O phase transition must occur in this small region where the FFLO-ADG-O-1 exists, due to the fact that the ADG turn into unstable when α > 0.0643 and the FFLO-ADG-O state vanishes when α < 0.0616. The left vertical line in the middle panel at T/T_F = 0.001/0.003 corresponds to the A-O transition asymmetry. Thanks to the coexistence region of ADG and FFLO-ADG-O state, the A-O phase transition is of the first order as shown in figure 7. One should note that the stable ADG and FFLO-ADG-O states coexist in the region where the FFLO-ADG-O-1 state exists. With increasing temperature, the existence region of the FFLO-ADG-O-1 state is getting smaller and vanishes when T/T_F > 0.005 48. And the A-O phase transition turns into the second order as shown in figure 7.

Zoom In Zoom Out Reset image size

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 neutron–proton pairing in attractive ³ SD₁ channel [47]. To some extent, the possible orientation of the Cooper pair momentum depends on the configuration of the pairing gap. For the current pairing gap, the orientation between θ₀ = 0 and ${\theta }_{0}=\frac{\pi }{2}$ can survive in an appreciable region in the T − α phase diagram as displayed in figure 7. We call this state the FFLO-ADG-B state.

The stability condition of FFLO-ADG-B state require the Hessian matrix is positive definite, i.e., $\mathrm{det}\,{\mathcal{M}}_{t} > 0$, $\mathrm{det}\,{\mathcal{M}}_{q} > 0$, $\mathrm{det}\,{\mathcal{M}}_{{\Delta}} > 0$, and $\mathrm{det}\,{\mathcal{M}}_{h} > 0$. Fortunately, $\mathrm{det}\,{\mathcal{M}}_{q}$, $\mathrm{det}\,{\mathcal{M}}_{{\Delta}}$, and $\mathrm{det}\,{\mathcal{M}}_{h}$ are always positive for all the possible solutions. The key stability condition becomes $\mathrm{det}\,{\mathcal{M}}_{t} > 0$. At Low temperature, there are three branches of solutions labelled by FFLO-ADG-B-1, FFLO-ADG-B-2, and FFLO-ADG-B-3 in figure 4(c). It is obvious that only the FFLO-ADG-B-2 state is related to the local minimum of free energy in hyperplane $\mathcal{P}({\Delta},q,{\theta }_{0})$. While the other two branches correspond to the saddle point solution. Same as the FFLO-ADG-O state, unless otherwise specified, the FFLO-ADG-B state refers to the FFLO-ADG-B-2 state hereafter. One should note that the behavior of the FFLO-ADG-B-1 and FFLO-ADG-B-3 states are very similar to that of the FFLO-ADG-O-1 state, for instance, the FFLO-ADG-B-1 (FFLO-ADG-O-1) vanishes at the zero point of $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}{\vert }_{{\theta }_{0}=\pi /2}$ $(\frac{{\partial }^{2}{\Omega}}{\partial {q}^{2}}{\vert }_{q=0})$. This behavior indicates the first phase transition from the FFLO-ADG-O state to the FFLO-ADG-B state.

To understand the phase transition from the FFLO-ADG-B state to the FFLO-ADG-O (FFLO-ADG-P) state, the scaled pairing gap and the orientation of q for the FFLO-ADG-B state, and the scaled orientation susceptibility of the Cooper pair momentum $\frac{{\partial }^{2}{\Omega}}{\partial {\theta }_{0}^{2}}$ for the FFLO-ADG-O and FFLO-ADG-P states are exhibited in figure 6. At low temperature T/T_F = 0.012, the existence of the FFLO-ADG-B-1 solution implies a sudden change of θ₀ from the stable FFLO-ADG-O state to the stable FFLO-ADG-B state. This sudden change corresponds to the first phase transition from the FFLO-ADG-O state to the FFLO-ADG-B state. We call this transition the O-B phase transition. Similarly, we call the transition from the FFLO-ADG-B state to the FFLO-ADG-P state the B-P phase transition. The FFLO-ADG-B-3 state results in the first order B-P transition. With increasing temperature, the unstable FFLO-ADG-B-1 solution vanishes, indicating the second O-B phase transition. Interestingly, at temperature T/T_F = 0.022, the orientation susceptibility is negative in the whole existing region, indicating that the FFLO-ADG-P state vanishes as shown in 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition, the stability against PS can be determined by the $\kappa ({\Delta},q,{\theta }_{0})=\mathrm{det}\,\mathcal{M} > 0$ or $\chi =\mathcal{D}\delta \rho /\mathcal{D}h > 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].

The stability of the polarized superfluid against PS requires the number susceptibility matrix must be positive definite. As discussed in section 2.5, both $\chi =\frac{\mathcal{D}\delta \rho }{\mathcal{D}h} > 0$ and $\kappa ({\Delta},q,{\theta }_{0})=\mathrm{det}\,\mathcal{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 ${\mathcal{M}}_{2\times 2}$ for FFLO state with ${\theta }_{0}=\frac{\pi }{2}$ are plotted in figure 8. The denominator in the number susceptibility (B7) evolves into $\kappa ({\Delta},q)=\frac{{\partial }^{2}{\Omega}}{\partial {{\Delta}}^{2}}\frac{{\partial }^{2}{\Omega}}{\partial {q}^{2}}-{(\frac{{\partial }^{2}{\Omega}}{\partial {\Delta}\partial q})}^{2}$ in hyperplane $P({\Delta},q,{\theta }_{0}=\frac{\pi }{2})$. The solid, short dash-dotted, short dotted, and short dashed lines label the FFLO-ADG-O-1, FFLO-ADG-O-2, FFLO-ADG-O-3, and FFLO-ADG-O-4 solutions, respectively. Since λ₊ is always positive, λ₋ determines the positive definiteness of ${\mathcal{M}}_{2\times 2}$. One can realize from figure 8 that the imbalance number susceptibility χ and λ₋ always keep the same positiveness and change signs synchronously, signifying that the criterion $\frac{\mathcal{D}\delta \rho }{\mathcal{D}h} > 0$ is completely equivalent to the positive eigenvalue λ₋. In contrast, κ and λ₋ hold the same positiveness and alter signs only for the stable solutions in the adopted parameter space, i.e., the FFLO-ADG-O-2 and FFLO-ADG-O-4 solutions. In the FFLO-ADG-O-1 and FFLO-ADG-O-3 states corresponding to the saddle point in hyperplane $P({\Delta},q,{\theta }_{0}=\frac{\pi }{2})$, κ 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 $\mathcal{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.

Zoom In Zoom Out Reset image size