Strange quark mass turns magnetic domain walls into multi-winding flux tubes

In three-flavor quark matter with sufficiently small mismatch in the Fermi momenta of the different fermion species, Cooper pairing predominantly occurs in the spin-zero channel and in the antisymmetric anti-triplet channels in color and flavor space, ${\left[\bar{3}\right]}_{\mathrm{c}}$ and ${\left[\bar{3}\right]}_{\mathrm{f}}$. As pairing is assumed to occur between fermions of the same chirality, the flavor channel stands for either left-handed or right-handed fermions. Therefore, the order parameter for Cooper pair condensation can be written as

$$\begin{equation}{\Psi}={{\Phi}}_{ij}{J}_{i}\otimes {I}_{j}\in {\left[\bar{3}\right]}_{\mathrm{c}}\otimes {\left[\bar{3}\right]}_{\mathrm{f}},\end{equation} \tag{ 1 }$$

where the anti-symmetric 3 × 3 matrices ${\left({J}_{i}\right)}_{jk}=-i{{\epsilon}}_{ijk}$ and ${\left({I}_{j}\right)}_{k\ell }=-i{{\epsilon}}_{jk\ell }$ form bases of the three-dimensional spaces ${\left[\bar{3}\right]}_{\mathrm{c}}$ and ${\left[\bar{3}\right]}_{\mathrm{f}}$, respectively. As a consequence, we can characterize a color-superconducting phase by the 3 × 3 matrix Φ, which has one (anti-)color and one (anti-)flavor index. We shall put the colors in the order (r, g, b) (for red, green, blue) and the flavors in the order (u, d, s) (for up, down, strange). The color charges are only labels and thus their order is not crucial, but the order of the flavors matters since electric charge and quark masses break the flavor symmetry ² . In this convention, for instance, Φ₁₁ carries anti-indices $\bar{r}$ and $\bar{u}$ and thus describes pairing of gd with bs quarks and of gs with bd quarks.

We consider the following Ginzburg–Landau potential up to fourth order in Φ,

$$\begin{align}\hfill U& =-12\left\{\mathrm{T}\mathrm{r}\left[{\left({D}_{0}{\Phi}\right)}^{{\dagger}}\left({D}_{0}{\Phi}\right)\right]+{u}^{2}\enspace \mathrm{T}\mathrm{r}\left[{\left({D}_{i}{\Phi}\right)}^{{\dagger}}\left({D}^{i}{\Phi}\right)\right]\right\}\hfill \\ \hfill & \quad -4{\mu }^{2}\enspace \mathrm{T}\mathrm{r}\left[{{\Phi}}^{{\dagger}}{\Phi}\right]+4\sigma \enspace \mathrm{T}\mathrm{r}\left[{{\Phi}}^{{\dagger}}\left(\frac{{M}^{2}}{2{\mu }_{\mathrm{q}}}+{\mu }_{\mathrm{e}}Q\right){\Phi}\right]\hfill \\ \hfill & \quad +16\left(\lambda +h\right)\mathrm{T}\mathrm{r}\left[{\left({{\Phi}}^{{\dagger}}{\Phi}\right)}^{2}\right]-16h{\left(\mathrm{T}\mathrm{r}\left[{{\Phi}}^{{\dagger}}{\Phi}\right]\right)}^{2}+\frac{1}{4}{F}_{\mu \nu }^{a}{F}_{a}^{\mu \nu }+\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }.\hfill \end{align} \tag{ 2 }$$

Apart from the mass correction, proportional to the parameter σ, this is exactly the same potential, and the same notation, as in reference [17], where the starting point was the potential for Ψ, based on previous works [14, 25, 26, 28]. Due to the broken Lorentz invariance in the medium, the temporal and spatial parts of the kinetic term have different prefactors, $u=1/\sqrt{3}$. The covariant derivative is given by

$$\begin{equation}{D}_{\mu }{\Phi}={\partial }_{\mu }{\Phi}-ig{A}_{\mu }^{a}{T}_{a}^{T}{\Phi}+ie{A}_{\mu }{\Phi}Q,\end{equation} \tag{ 3 }$$

where g and e are the strong and electromagnetic coupling constants, respectively. The color gauge fields are denoted by ${A}_{\mu }^{a}$, a = 1, ..., 8, and A_μ is the electromagnetic gauge field. Furthermore, T_a = λ_a /2, with the Gell–Mann matrices λ_a , are the generators of the color gauge group SU(3), and the electric charge matrix for the Cooper pairs Q = diag(q_d + q_s, q_u + q_s, q_u + q_d) = diag(−2/3, 1/3, 1/3) is the generator of the electromagnetic gauge group U(1). Here, q_u, q_d, q_s denote the individual quark charges in units of the elementary charge e. The field strength tensors are ${F}_{\mu \nu }^{a}={\partial }_{\mu }{A}_{\nu }^{a}-{\partial }_{\nu }{A}_{\mu }^{a}+g{f}^{abc}{A}_{\mu }^{b}{A}_{\nu }^{c}$ for the color sector, where f^abc are the SU(3) structure constants, and F_μν = ∂_μ A_ν − ∂_ν A_μ for the electromagnetic sector. The constants in front of the quadratic and quartic terms in Φ are written conveniently as combinations of μ, λ and h, whose physical meaning will become obvious after performing the traces (see equation (4)).

We have incorporated a mass correction to lowest order, with the (fermionic) quark chemical potential μ_q and the mass matrix for the Cooper pairs M = diag(m_d + m_s, m_u + m_s, m_u + m_d) ≃ diag(m_s, m_s, 0), i.e. we shall neglect the masses of the light quarks, m_u ≃ m_d ≃ 0, and keep the strange quark mass m_s as a free parameter. We have also included the contribution of the electric charge chemical potential μ_e, which is of the same order for small quark masses, ${\mu }_{\mathrm{e}}\propto {m}_{\mathrm{s}}^{2}/{\mu }_{\mathrm{q}}$. The correction term is identical to the one used in references [30, 31], which is easily seen by an appropriate rescaling of Φ.

In the following we restrict ourselves to diagonal order parameters, ${\Phi}=\frac{1}{2}\enspace \mathrm{diag}\left({\phi }_{1},{\phi }_{2},{\phi }_{3}\right)$ with ${\phi }_{i}\in \mathbb{C}$, where ϕ₁ corresponds to ds pairing, ϕ₂ to us pairing, and ϕ₃ to ud pairing. If flavor symmetry was intact, off-diagonal order parameters could always be brought into a diagonal form by an appropriate rotation in color and flavor space. This is no longer true when flavor symmetry is explicitly broken, and thus our restriction to diagonal order parameters is a simplification of the most general situation [48]. Even in this simplified case, 2³ = 8 qualitatively different homogeneous phases have to be considered in principle, accounting for each of the three condensates to be either zero or nonzero. The restriction to diagonal matrices Φ allows us to consistently set all gauge fields with off-diagonal components to zero, ${A}_{\mu }^{1}={A}_{\mu }^{2}={A}_{\mu }^{4}={A}_{\mu }^{5}={A}_{\mu }^{6}={A}_{\mu }^{7}=0$, such that the only relevant gauge fields are the two color gauge fields ${A}_{\mu }^{3}$, ${A}_{\mu }^{8}$, and the electromagnetic gauge field A_μ . Moreover, we are only interested in static solutions and drop all electric fields, i.e. we only keep the spatial components of the gauge fields, giving rise to the magnetic fields B₃ = ∇ × A₃, B₈ = ∇ × A₈, and B = ∇ × A.

Within this ansatz, performing the traces in equation (2) yields

$$\begin{align}\hfill U& ={\left\vert \left(\nabla +\frac{ig{\mathbf{A}}_{3}}{2}+\frac{ig{\mathbf{A}}_{8}}{2\sqrt{3}}+\frac{2ie\mathbf{A}}{3}\right){\phi }_{1}\right\vert }^{2}+{\left\vert \left(\nabla -\frac{ig{\mathbf{A}}_{3}}{2}+\frac{ig{\mathbf{A}}_{8}}{2\sqrt{3}}-\frac{ie\mathbf{A}}{3}\right){\phi }_{2}\right\vert }^{2}\hfill \\ \hfill & \quad +{\left\vert \left(\nabla -\frac{ig{\mathbf{A}}_{8}}{\sqrt{3}}-\frac{ie\mathbf{A}}{3}\right){\phi }_{3}\right\vert }^{2}-\left({\mu }^{2}-{m}_{1}^{2}\right)\vert {\phi }_{1}{\vert }^{2}-\left({\mu }^{2}-{m}_{2}^{2}\right)\vert {\phi }_{2}{\vert }^{2}\hfill \\ \hfill & \quad -\left({\mu }^{2}-{m}_{3}^{2}\right)\vert {\phi }_{3}{\vert }^{2}+\lambda \left(\vert {\phi }_{1}{\vert }^{4}+\vert {\phi }_{2}{\vert }^{4}+\vert {\phi }_{3}{\vert }^{4}\right)\hfill \\ \hfill & \quad -2h\left(\vert {\phi }_{1}{\vert }^{2}\vert {\phi }_{2}{\vert }^{2}+\vert {\phi }_{1}{\vert }^{2}\vert {\phi }_{3}{\vert }^{2}+\vert {\phi }_{2}{\vert }^{2}\vert {\phi }_{3}{\vert }^{2}\right)+\frac{{\mathbf{B}}_{3}^{2}}{2}+\frac{{\mathbf{B}}_{8}^{2}}{2}+\frac{{\mathbf{B}}^{2}}{2}.\hfill \end{align} \tag{ 4 }$$

This potential can be viewed as a generalized version of a textbook superconductor which has a single condensate coupled to a single gauge field, see for instance reference [49]. Here we have three condensates with identical (bosonic) chemical potential μ, with self-coupling λ, cross-coupling h, and three different effective masses (squared)

$$\begin{equation}{m}_{1}^{2}=\sigma \left(\frac{{m}_{\mathrm{s}}^{2}}{2{\mu }_{\mathrm{q}}}-\frac{2{\mu }_{\mathrm{e}}}{3}\right),\quad {m}_{2}^{2}=\sigma \left(\frac{{m}_{\mathrm{s}}^{2}}{2{\mu }_{\mathrm{q}}}+\frac{{\mu }_{\mathrm{e}}}{3}\right),\quad {m}_{3}^{2}=\sigma \frac{{\mu }_{\mathrm{e}}}{3}.\end{equation} \tag{ 5 }$$

All three gauge fields couple to the condensates. We can simplify the potential by a suitable rotation of the gauge fields.

We apply the double rotation from reference [17], which, denoting the rotated gauge fields by ${\tilde {A}}_{\mu }^{3}$, ${\tilde {A}}_{\mu }^{8}$, and ${\tilde {A}}_{\mu }$, reads,

$$\begin{equation}\left(\begin{matrix}\hfill {\tilde {A}}_{\mu }^{3}\hfill \\ \hfill {\tilde {A}}_{\mu }^{8}\hfill \\ \hfill {\tilde {A}}_{\mu }\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill \mathrm{cos}\enspace {\vartheta }_{2}\hfill & \hfill 0\hfill & \hfill \mathrm{sin}\enspace {\vartheta }_{2}\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill -\mathrm{sin}\enspace {\vartheta }_{2}\hfill & \hfill 0\hfill & \hfill \mathrm{cos}\enspace {\vartheta }_{2}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathrm{cos}\enspace {\vartheta }_{1}\hfill & \hfill \mathrm{sin}\enspace {\vartheta }_{1}\hfill \\ \hfill 0\hfill & \hfill -\mathrm{sin}\enspace {\vartheta }_{1}\hfill & \hfill \mathrm{cos}\enspace {\vartheta }_{1}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {A}_{\mu }^{3}\hfill \\ \hfill {A}_{\mu }^{8}\hfill \\ \hfill {A}_{\mu }\hfill \end{matrix}\right),\end{equation} \tag{ 6 }$$

where the mixing angles ϑ₁ and ϑ₂ are given by

$$\begin{equation}\mathrm{sin}\enspace {\vartheta }_{1}=\frac{e}{\sqrt{3{g}^{2}+{e}^{2}}},\quad \mathrm{cos}\enspace {\vartheta }_{1}=\frac{\sqrt{3}g}{\sqrt{3{g}^{2}+{e}^{2}}},\end{equation} \tag{ 7a }$$

$$\begin{equation}\mathrm{sin}\enspace {\vartheta }_{2}=\frac{\sqrt{3}e}{\sqrt{3{g}^{2}+4{e}^{2}}},\quad \mathrm{cos}\enspace {\vartheta }_{2}=\frac{\sqrt{3{g}^{2}+{e}^{2}}}{\sqrt{3{g}^{2}+4{e}^{2}}}.\end{equation} \tag{ 7b }$$

This rotation is the most convenient choice for our purpose of calculating flux tube profiles in the 2SC phase: the first rotation, with the usual '2SC mixing angle' ϑ₁, ensures that in the homogeneous 2SC phase, where only the condensate ϕ₃ is nonzero, only the ${\tilde {\mathbf{B}}}_{8}$ field is expelled. The other two rotated fields penetrate the superconductor unperturbed (assuming zero magnetization from the unpaired quarks). If we were only interested in the homogeneous 2SC phase, this rotation would be sufficient. However, we will allow for ϕ₁ and ϕ₂ to be induced in the core of the flux tube. Therefore, we apply a second rotation with mixing angle ϑ₂. This rotation leaves ${\tilde {\mathbf{B}}}_{8}$ invariant and creates a field, namely $\tilde {\mathbf{B}}$, which is unaffected by the superconductor even if all three condensates are nonzero. Thus, $\tilde {\mathbf{B}}$ simply decouples from the condensates and can be ignored in the calculation of the flux tube profiles. For g ≫ e both mixing angles are small and thus in this case ${\tilde {A}}_{\mu }^{3}$ and ${\tilde {A}}_{\mu }^{8}$ are 'almost' gluons with a small admixture of the photon, while ${\tilde {A}}_{\mu }$ is 'almost' the photon with a small gluonic admixture. For consistency, we shall work with the rotated fields (6) throughout the paper, including the discussion of the homogeneous phases in section 3.

Applying the gauge field rotation and writing the complex fields in terms of their moduli and phases,

$$\begin{equation}{\phi }_{i}\left(\mathbf{r}\right)=\frac{{\rho }_{i}\left(\mathbf{r}\right)}{\sqrt{2}}\enspace {\text{e}}^{i{\psi }_{i}\left(\mathbf{r}\right)},\quad i=1,2,3,\end{equation} \tag{ 8 }$$

the Ginzburg–Landau potential (4) becomes

$$\begin{equation}U={U}_{0}+\frac{{\tilde {\mathbf{B}}}_{3}^{2}}{2}+\frac{{\tilde {\mathbf{B}}}_{8}^{2}}{2}+\frac{{\tilde {\mathbf{B}}}^{2}}{2},\end{equation} \tag{ 9 }$$

with

$$\begin{align}\hfill {U}_{0}& ={\left(\nabla {\psi }_{1}+{\tilde {q}}_{3}{\tilde {\mathbf{A}}}_{3}+{\tilde {q}}_{81}{\tilde {\mathbf{A}}}_{8}\right)}^{2}\frac{{\rho }_{1}^{2}}{2}+{\left(\nabla {\psi }_{2}-{\tilde {q}}_{3}{\tilde {\mathbf{A}}}_{3}+{\tilde {q}}_{82}{\tilde {\mathbf{A}}}_{8}\right)}^{2}\frac{{\rho }_{2}^{2}}{2}\hfill \\ \hfill & \quad +{\left(\nabla {\psi }_{3}-{\tilde {q}}_{83}{\tilde {\mathbf{A}}}_{8}\right)}^{2}\frac{{\rho }_{3}^{2}}{2}+\frac{{\left(\nabla {\rho }_{1}\right)}^{2}}{2}+\frac{{\left(\nabla {\rho }_{2}\right)}^{2}}{2}+\frac{{\left(\nabla {\rho }_{3}\right)}^{2}}{2}-\frac{{\mu }^{2}-{m}_{1}^{2}}{2}{\rho }_{1}^{2}\hfill \\ \hfill & \quad -\frac{{\mu }^{2}-{m}_{2}^{2}}{2}{\rho }_{2}^{2}-\frac{{\mu }^{2}-{m}_{3}^{2}}{2}{\rho }_{3}^{2}+\frac{\lambda }{4}\left({\rho }_{1}^{4}+{\rho }_{2}^{4}+{\rho }_{3}^{4}\right)\hfill \\ \hfill & \quad -\frac{h}{2}\left({\rho }_{1}^{2}{\rho }_{2}^{2}+{\rho }_{1}^{2}{\rho }_{3}^{2}+{\rho }_{2}^{2}{\rho }_{3}^{2}\right),\hfill \end{align} \tag{ 10 }$$

where we have introduced the rotated charges

$$\begin{equation}{\tilde {q}}_{81}\equiv \frac{g}{2\sqrt{3}}\enspace \mathrm{cos}\enspace {\vartheta }_{1}+\frac{2e}{3}\enspace \mathrm{sin}\enspace {\vartheta }_{1}=\frac{3{g}^{2}+4{e}^{2}}{6\sqrt{3{g}^{2}+{e}^{2}}},\end{equation} \tag{ 11a }$$

$$\begin{equation}{\tilde {q}}_{82}\equiv \frac{g}{2\sqrt{3}}\enspace \mathrm{cos}\enspace {\vartheta }_{1}-\frac{e}{3}\enspace \mathrm{sin}\enspace {\vartheta }_{1}=\frac{3{g}^{2}-2{e}^{2}}{6\sqrt{3{g}^{2}+{e}^{2}}},\end{equation} \tag{ 11b }$$

$$\begin{equation}{\tilde {q}}_{83}\equiv \frac{g}{\sqrt{3}}\enspace \mathrm{cos}\enspace {\vartheta }_{1}+\frac{e}{3}\enspace \mathrm{sin}\enspace {\vartheta }_{1}=\frac{\sqrt{3{g}^{2}+{e}^{2}}}{3},\end{equation} \tag{ 11c }$$

$$\begin{equation}{\tilde {q}}_{3}\equiv \frac{g}{2}\enspace \mathrm{cos}\enspace {\vartheta }_{2}+\frac{e}{2}\enspace \mathrm{cos}\enspace {\vartheta }_{1}\enspace \mathrm{sin}\enspace {\vartheta }_{2}=\frac{g}{2}\frac{\sqrt{3{g}^{2}+4{e}^{2}}}{\sqrt{3{g}^{2}+{e}^{2}}}.\end{equation} \tag{ 11d }$$

In equation (9) we have separated the quadratic contributions of the magnetic fields, which is notationally convenient for the following.

We shall be interested in the phase structure at fixed external magnetic field H, which we assume to be homogeneous and along the z-direction, H(r) = H e_z . Therefore, we need to consider the Gibbs free energy density

$$\begin{equation}G=\frac{1}{V}\int {\mathrm{d}}^{3}\mathbf{r}\enspace \left(U-\mathbf{H}\cdot \mathbf{B}\right),\end{equation} \tag{ 12 }$$

where V is the volume of our system. We can obviously assume that all induced magnetic fields have only z-components as well. Denoting the z-components of the rotated fields by ${\tilde {B}}_{3}$, ${\tilde {B}}_{8}$, and $\tilde {B}$, we have $\mathbf{H}\cdot \mathbf{B}=H\left[\mathrm{sin}\enspace {\vartheta }_{1}\enspace {\tilde {B}}_{8}+\mathrm{cos}\enspace {\vartheta }_{1}\left(\mathrm{sin}\enspace {\vartheta }_{2}\enspace {\tilde {B}}_{3}+\mathrm{cos}\enspace {\vartheta }_{2}\enspace \tilde {B}\right)\right]$. Since $\tilde {\mathbf{B}}$ does not couple to any of the condensates, it remains homogeneous even in the presence of flux tubes. Consequently, the equation of motion for $\tilde {\mathbf{A}}$ is trivially fulfilled, and we determine $\tilde {B}$ by minimizing the Gibbs free energy, which, using equation (9), yields

$$\begin{equation}\tilde {B}=H\enspace \mathrm{cos}\enspace {\vartheta }_{1}\enspace \mathrm{cos}\enspace {\vartheta }_{2}.\end{equation} \tag{ 13 }$$

Reinserting this result into G, we obtain

$$\begin{align}\hfill G& =-\frac{{H}^{2}\enspace {\mathrm{cos}}^{2}\enspace {\vartheta }_{1}\enspace {\mathrm{cos}}^{2}\enspace {\vartheta }_{2}}{2}\hfill \\ \hfill & \quad +\frac{1}{V}\int {\mathrm{d}}^{3}\mathbf{r}\enspace \left[{U}_{0}+\frac{{\tilde {B}}_{3}^{2}}{2}+\frac{{\tilde {B}}_{8}^{2}}{2}-H\left(\mathrm{sin}\enspace {\vartheta }_{1}\enspace {\tilde {B}}_{8}+\mathrm{cos}\enspace {\vartheta }_{1}\enspace \mathrm{sin}\enspace {\vartheta }_{2}\enspace {\tilde {B}}_{3}\right)\right].\hfill \end{align} \tag{ 14 }$$

The potential U₀ (10) depends on the parameters μ, λ, h, σ. The discussion of the homogeneous phases in section 3 turns out to be sufficiently simple to keep these parameters unspecified and to investigate the general phase structure. Our main results, however, require the numerical calculation of the flux tube profiles, and a completely general study would be extremely laborious. Therefore, for the results in section 5, we employ the weak-coupling values of these parameters [14, 15, 28, 30, 31, 38],

$$\begin{align}\hfill & {\mu }^{2}=\frac{48{\pi }^{2}}{7\zeta \left(3\right)}{T}_{\mathrm{c}}^{2}\left(1-\frac{T}{{T}_{\mathrm{c}}}\right),\quad \lambda =\frac{72{\pi }^{4}}{7\zeta \left(3\right)}\frac{{T}_{\mathrm{c}}^{2}}{{\mu }_{\mathrm{q}}^{2}},\quad \hfill \\ \hfill & h=-\frac{36{\pi }^{4}}{7\zeta \left(3\right)}\frac{{T}_{\mathrm{c}}^{2}}{{\mu }_{\mathrm{q}}^{2}},\quad \sigma =-\frac{24{\pi }^{2}}{7\zeta \left(3\right)}\frac{{T}_{\mathrm{c}}^{2}}{{\mu }_{\mathrm{q}}}\enspace \mathrm{ln}\enspace \frac{{T}_{\mathrm{c}}}{{\mu }_{\mathrm{q}}},\hfill \end{align} \tag{ 15 }$$

where T_c is the critical temperature. The ratio T_c/μ_q can be understood as a measure of the pairing strength since T_c is closely related to the pairing gap. For instance, at weak coupling, which is applicable at asymptotically large densities, the zero-temperature pairing gap is exponentially suppressed compared to μ_q. It is related by a numerical factor of order one to T_c [50], and thus T_c/μ_q is also exponentially small. We shall extrapolate our results to strong coupling, having in mind applications to compact stars, where the densities are large, but not asymptotically large. In this case, model calculations as well as extrapolations of perturbative results suggest that T_c is of the order of tens of MeV, while we expect μ_q ∼ (400–500) MeV in the cores of compact stars, which yields the estimate T_c/μ_q ∼ 0.1. Besides the implicit dependence on T_c/μ_q our potential also depends on the ratio m_s/μ_q. Since m_s is medium dependent, its value at non-asymptotic densities is poorly known. It is expected to be somewhere between the current mass and the constituent mass within a baryon, m_s ∼ (100–500) MeV; for a concrete calculation within the Nambu–Jona-Lasinio model see for instance reference [51]. With the range of the quark chemical potential given above we thus expect m_s/μ_q ∼ (0.2–1). Finally, our potential depends on the electric charge chemical potential μ_e. In a fermionic approach, this chemical potential would be determined from the conditions of beta-equilibrium and charge neutrality. Since our Ginzburg–Landau expansion is formally based on small values of the order parameter, we follow references [30, 31] and use the value of μ_e in the completely unpaired phase. At weak coupling and to lowest order in the strange quark mass this value is (see for instance references [52, 53])

$$\begin{equation}{\mu }_{\mathrm{e}}=\frac{{m}_{\mathrm{s}}^{2}}{4{\mu }_{\mathrm{q}}}.\end{equation} \tag{ 16 }$$

With this result, it is convenient to trade the dimensionful parameter σ for the dimensionless 'mass parameter'

$$\begin{equation}\alpha \equiv \frac{\sigma {m}_{\mathrm{s}}^{2}}{{\mu }^{2}{\mu }_{\mathrm{q}}}=\frac{{m}_{\mathrm{s}}^{2}}{2{\mu }_{\mathrm{q}}^{2}}{\left(1-\frac{T}{{T}_{\mathrm{c}}}\right)}^{-1}\enspace \mathrm{ln}\enspace \frac{{\mu }_{\mathrm{q}}}{{T}_{\mathrm{c}}},\end{equation} \tag{ 17 }$$

such that the complete dependence of our potential on the strange quark mass is absorbed in α,

$$\begin{equation}{m}_{1}^{2}=\frac{{\mu }^{2}}{3}\alpha ,\quad {m}_{2}^{2}=\frac{7{\mu }^{2}}{12}\alpha ,\quad {m}_{3}^{2}=\frac{{\mu }^{2}}{12}\alpha .\end{equation} \tag{ 18 }$$

We shall see that if we are only interested in homogeneous phases, the phase structure is most conveniently calculated in the space spanned by α, g, the normalized dimensionless magnetic field H/(μ²/λ^1/2), and the ratio

$$\begin{equation}\eta \equiv \frac{h}{\lambda }.\end{equation} \tag{ 19 }$$

To ensure the boundedness of our potential, we must require η < 1/2. At weak coupling η = −1/2, as one can see from equation (15). Later, in our explicit calculation of the flux tube profiles and the resulting critical magnetic fields, we consider a fixed g and the parameter space spanned by H/(μ²/λ^1/2), T_c/μ_q, and m_s/μ_q. To choose a value of g, realistic for compact star conditions, we observe that according to the two-loop QCD beta function (which should not be taken too seriously at such low densities), μ_q ≃ 400 MeV corresponds to α_s ≃ 1 and thus $g=\sqrt{4\pi {\alpha }_{\mathrm{s}}}\simeq 3.5$. Of course, choosing such a large value for g in our main results is a bold extrapolation, given that we work with the weak-coupling parameters (15). Furthermore, we shall set T = 0 in equation (17). Strictly speaking this is inconsistent because the Ginzburg–Landau potential is an expansion in the condensates, and we use a value for μ_e (16) that is only valid very close to a second-order transition to the unpaired phase. Choosing a different, nonzero temperature, would not change our result qualitatively because it only enters the relation between m_s/μ_q and α. The definition of α (17) shows that the mass effect is smallest for zero temperature (i.e. in this case α is smallest for a given m_s/μ_q). Therefore, by our choice T = 0 in equation (17) we will obtain an upper limit in m_s/μ_q for the presence of multi-winding 2SC flux tubes. Any T > 0 in equation (17) will give a smaller m_s/μ_q up to which these exotic configurations exist. In any case, the temperature dependence in the present approach is somewhat simplistic to begin with because, firstly, in a multi-component superconductor there can be different critical temperatures for the different condensates, resulting in temperature factors different from the standard Ginzburg–Landau formalism [21]. And, secondly, away from the asymptotic weak-coupling regime the phase transition becomes first order due to gauge field fluctuations [54], and thus at strong coupling the behavior just below the phase transition would have to be modified in a more sophisticated approach.

In the CFL phase, all three condensates are nonzero. Here we obtain ${\tilde {B}}_{3}={\tilde {B}}_{8}=0$ and equation (21) yield three coupled equations for the condensates, which have the solution

$$\begin{equation}{\rho }_{1}^{2}=\frac{\lambda \left({\mu }^{2}-{m}_{1}^{2}\right)+h\left({\mu }^{2}+{m}_{1}^{2}-{m}_{2}^{2}-{m}_{3}^{2}\right)}{\left(\lambda -2h\right)\left(\lambda +h\right)}=\frac{{\mu }^{2}}{\lambda \left(1-2\eta \right)}\left(1-\frac{1}{3}\alpha \right),\end{equation} \tag{ 23a }$$

$$\begin{equation}{\rho }_{2}^{2}=\frac{\lambda \left({\mu }^{2}-{m}_{2}^{2}\right)+h\left({\mu }^{2}-{m}_{1}^{2}+{m}_{2}^{2}-{m}_{3}^{2}\right)}{\left(\lambda -2h\right)\left(\lambda +h\right)}=\frac{{\mu }^{2}}{\lambda \left(1-2\eta \right)}\left[1-\frac{\left(7-2\eta \right)\alpha }{12\left(1+\eta \right)}\right],\end{equation} \tag{ 23b }$$

$$\begin{equation}{\rho }_{3}^{2}=\frac{\lambda \left({\mu }^{2}-{m}_{3}^{2}\right)+h\left({\mu }^{2}-{m}_{1}^{2}-{m}_{2}^{2}+{m}_{3}^{2}\right)}{\left(\lambda -2h\right)\left(\lambda +h\right)}=\frac{{\mu }^{2}}{\lambda \left(1-2\eta \right)}\left[1-\frac{\left(1+10\eta \right)\alpha }{12\left(1+\eta \right)}\right],\end{equation} \tag{ 23c }$$

where, in the second steps, equations (18) and (19) have been used ³ . In the massless limit, α = 0, we recover the result of three identical condensates. The potential (10) becomes

$$\begin{equation}{U}_{0}\simeq -\frac{3{\mu }^{4}}{4\lambda \left(1-2\eta \right)}\left(1-\frac{2}{3}\alpha \right),\end{equation} \tag{ 24 }$$

where we have dropped the contribution quadratic in α, i.e. all terms quartic in m_s. This is consistent with our starting point, which only includes corrections to linear order in α. The Gibbs free energy density (14) thus becomes

$$\begin{equation}{G}_{\text{CFL}}\simeq -\frac{3{g}^{2}{H}^{2}}{2\left(3{g}^{2}+4{e}^{2}\right)}-\frac{3{\mu }^{4}}{4\lambda \left(1-2\eta \right)}\left(1-\frac{2}{3}\alpha \right),\end{equation} \tag{ 25 }$$

where we have used the explicit form of the mixing angles (7).

Next we discuss the three possible phases where exactly one condensate is non-vanishing. In each case, two flavors and two colors participate in pairing. According to the flavor structure of the condensates we term the phases 2SC_ds, 2SC_us, or 2SC_ud if ρ₁, ρ₂, or ρ₃ is non-vanishing, respectively. As we shall see, only the 2SC_ud phase is relevant for the phase structure since the other two phases turn out to be energetically disfavored. Therefore, we often use 2SC synonymously for 2SC_ud.

Starting with the most relevant phase, we first consider a nonzero ρ₃. From equation (21c) we then obtain

$$\begin{equation}{\rho }_{3}^{2}=\frac{{\mu }^{2}-{m}_{3}^{2}}{\lambda }=\frac{{\mu }^{2}}{\lambda }\left(1-\frac{\alpha }{12}\right),\end{equation} \tag{ 26 }$$

and ${\tilde {B}}_{8}=0$. By inserting this into the Gibbs free energy density G and minimizing with respect to ${\tilde {B}}_{3}$ we find

$$\begin{equation}{\tilde {B}}_{3}=\frac{3egH}{\sqrt{3{g}^{2}+{e}^{2}}\sqrt{3{g}^{2}+4{e}^{2}}},\end{equation} \tag{ 27 }$$

and inserting this back into G yields

$$\begin{equation}{G}_{{2\text{SC}}_{\text{ud}}}\simeq -\frac{3{g}^{2}{H}^{2}}{2\left(3{g}^{2}+{e}^{2}\right)}-\frac{{\mu }^{4}}{4\lambda }\left(1-\frac{\alpha }{6}\right).\end{equation} \tag{ 28 }$$

Analogously, if only the condensate ρ₁ is nonzero, i.e. for the 2SC_ds phase, we find

$$\begin{equation}{\rho }_{1}^{2}=\frac{{\mu }^{2}}{\lambda }\left(1-\frac{1}{3}\alpha \right),\end{equation} \tag{ 29 }$$

as well as ${\tilde {B}}_{3}={\tilde {B}}_{8}=0$ from minimizing G, and we compute

$$\begin{equation}{G}_{{2\text{SC}}_{\mathrm{d}\mathrm{s}}}\simeq -\frac{3{g}^{2}{H}^{2}}{2\left(3{g}^{2}+4{e}^{2}\right)}-\frac{{\mu }^{4}}{4\lambda }\left(1-\frac{2}{3}\alpha \right).\end{equation} \tag{ 30 }$$

Finally, for the 2SC_us phase, we find

$$\begin{equation}{\rho }_{2}^{2}=\frac{{\mu }^{2}}{\lambda }\left(1-\frac{7}{12}\alpha \right),\end{equation} \tag{ 31 }$$

and the magnetic fields

$$\begin{equation}{\tilde {B}}_{3}=\frac{3eg\left(3{g}^{2}-2{e}^{2}\right)H}{2{\left(3{g}^{2}+{e}^{2}\right)}^{3/2}\sqrt{3{g}^{2}+4{e}^{2}}},\quad {\tilde {B}}_{8}=\frac{9e{g}^{2}H}{2{\left(3{g}^{2}+{e}^{2}\right)}^{3/2}}.\end{equation} \tag{ 32 }$$

This yields the Gibbs free energy density

$$\begin{equation}{G}_{{2\text{SC}}_{\text{us}}}\simeq -\frac{3{g}^{2}{H}^{2}}{2\left(3{g}^{2}+{e}^{2}\right)}-\frac{{\mu }^{4}}{4\lambda }\left(1-\frac{7}{6}\alpha \right).\end{equation} \tag{ 33 }$$

Comparing equations (28), (30) and (33), we see that for all values of the magnetic field H and all nonzero values of α, the 2SC_ud phase has the lowest Gibbs free energy and thus we can ignore the other two 2SC phases. In the massless limit α = 0 the 2SC_ud and 2SC_us phases become degenerate, as discussed in the introduction.

There are three possible phases in which exactly two condensates are non-vanishing. As for the 2SC phases, we shall see that one of these phases is favored over the other two for all parameter values. This is the phase where ρ₂ = 0. In this case, there are gd/bs and gs/bd Cooper pairs from ρ₁ and ru/gd and rd/gu Cooper pairs from ρ₃. Therefore, all three colors of the d quark are involved in pairing (while only two colors of the u and s quarks are involved), and thus, following reference [31], we refer to this phase as dSC. Analogously, the phases with vanishing ρ₁ and ρ₃ will be termed uSC and sSC, respectively. (Hence the collective term fSC, where f stands for the flavor.)

Starting again with the most relevant phase, we set ρ₂ = 0. With the help of equation (20) we obtain ${\tilde {B}}_{3}={\tilde {B}}_{8}=0$. Then, equations (21a) and (21c) yield two coupled equations for ρ₁ and ρ₃, with the solution

$$\begin{equation}{\rho }_{1}^{2}=\frac{\lambda \left({\mu }^{2}-{m}_{1}^{2}\right)+h\left({\mu }^{2}-{m}_{3}^{2}\right)}{{\lambda }^{2}-{h}^{2}}=\frac{{\mu }^{2}}{\lambda \left(1-\eta \right)}\left[1-\frac{\left(4+\eta \right)\alpha }{12\left(1+\eta \right)}\right],\end{equation} \tag{ 34a }$$

$$\begin{equation}{\rho }_{3}^{2}=\frac{\lambda \left({\mu }^{2}-{m}_{3}^{2}\right)+h\left({\mu }^{2}-{m}_{1}^{2}\right)}{{\lambda }^{2}-{h}^{2}}=\frac{{\mu }^{2}}{\lambda \left(1-\eta \right)}\left[1-\frac{\left(1+4\eta \right)\alpha }{12\left(1+\eta \right)}\right].\end{equation} \tag{ 34b }$$

Thus we compute

$$\begin{equation}{G}_{\text{dSC}}\simeq -\frac{3{g}^{2}{H}^{2}}{2\left(3{g}^{2}+4{e}^{2}\right)}-\frac{{\mu }^{4}}{2\lambda \left(1-\eta \right)}\left(1-\frac{5}{12}\alpha \right),\end{equation} \tag{ 35 }$$

where again we have dropped the quadratic contribution in α.

Analogously, if ρ₃ = 0, we find again ${\tilde {B}}_{3}={\tilde {B}}_{8}=0$, as well as

$$\begin{equation}{\rho }_{1}^{2}=\frac{{\mu }^{2}}{\lambda \left(1-\eta \right)}\left[1-\frac{\left(4+7\eta \right)\alpha }{12\left(1+\eta \right)}\right],\quad {\rho }_{2}^{2}=\frac{{\mu }^{2}}{\lambda \left(1-\eta \right)}\left[1-\frac{\left(7+4\eta \right)\alpha }{12\left(1+\eta \right)}\right],\end{equation} \tag{ 36 }$$

and

$$\begin{equation}{G}_{\text{sSC}}\simeq -\frac{3{g}^{2}{H}^{2}}{2\left(3{g}^{2}+4{e}^{2}\right)}-\frac{{\mu }^{4}}{2\lambda \left(1-\eta \right)}\left(1-\frac{11}{12}\alpha \right).\end{equation} \tag{ 37 }$$

Finally, for ρ₁ = 0 we find ${\tilde {B}}_{3}={\tilde {B}}_{8}=0$,

$$\begin{equation}{\rho }_{2}^{2}=\frac{{\mu }^{2}}{\lambda \left(1-\eta \right)}\left[1-\frac{\left(7+\eta \right)\alpha }{12\left(1+\eta \right)}\right],\quad {\rho }_{3}^{2}=\frac{{\mu }^{2}}{\lambda \left(1-\eta \right)}\left[1-\frac{\left(1+7\eta \right)\alpha }{12\left(1+\eta \right)}\right],\end{equation} \tag{ 38 }$$

and

$$\begin{equation}{G}_{\text{uSC}}\simeq -\frac{3{g}^{2}{H}^{2}}{2\left(3{g}^{2}+4{e}^{2}\right)}-\frac{{\mu }^{4}}{2\lambda \left(1-\eta \right)}\left(1-\frac{2}{3}\alpha \right).\end{equation} \tag{ 39 }$$

Since the Gibbs free energy densities (35), (37) and (39) of all three phases receive the same contribution from the magnetic field, the one with the lowest energy penalty from the strange quark mass term is preferred (and they all become degenerate for vanishing mass). We can therefore focus on the dSC phase and ignore the other two phases. While in the massless case none of these phases appears in the phase diagram [17], we find that a window for this phase opens up in the presence of a strange quark mass, which was already observed (without external magnetic field) in references [30, 31].

As we have just seen, the only relevant homogeneous phases within our ansatz are the CFL, 2SC_ud, dSC, and NOR phases. From their Gibbs free energy densities we can now easily obtain the critical magnetic fields H_c for the first-order transitions between them. In the usual terminology for ordinary superconductors, H_c is the critical magnetic field for the first-order transition between the superconducting phase and the normal-conducting phase for a type-I superconductor. For a color superconductor, the situation is more complicated because the various condensates can be broken sequentially until for sufficiently large magnetic field the normal phase is reached. As of now, we have not yet determined the transition from type-I to type-II behavior. To this end, the critical fields for the appearance of inhomogeneous phases have to be calculated, which we shall do in the subsequent sections.

By pairwise equating the Gibbs free energies (22), (25), (28), and (35) we obtain the conditions for 6 potential phase transitions. There is one transition which is independent of the magnetic field, namely the CFL/dSC transition, which we can express for instance in terms of a critical value for η,

$$\begin{equation}\eta =-\frac{6-7\alpha }{2\left(3+\alpha \right)}\quad \text{dSC/CFL}.\end{equation} \tag{ 40 }$$

Then, there are 4 transitions for which we can compute a critical magnetic field,

$$\begin{equation}\frac{{H}_{\mathrm{c}}^{2}}{{\mu }^{4}/\lambda }=\begin{cases}\left(1-\frac{1}{6}\alpha \right)\frac{3{g}^{2}+{e}^{2}}{2{e}^{2}}\quad \hfill & 2\text{SC/NOR}\hfill \\ \frac{3}{2\left(1-2\eta \right)}\left(1-\frac{2}{3}\alpha \right)\frac{3{g}^{2}+4{e}^{2}}{4{e}^{2}}\quad \hfill & \text{CFL/NOR}\hfill \\ \frac{1+\eta -\left(11+2\eta \right)\frac{\alpha }{12}}{1-2\eta }\frac{\left(3{g}^{2}+{e}^{2}\right)\left(3{g}^{2}+4{e}^{2}\right)}{9{e}^{2}{g}^{2}}\quad \hfill & 2\text{SC/CFL}\hfill \\ \frac{1+\eta -\left(4+\eta \right)\frac{\alpha }{6}}{2\left(1-\eta \right)}\frac{\left(3{g}^{2}+{e}^{2}\right)\left(3{g}^{2}+4{e}^{2}\right)}{9{e}^{2}{g}^{2}}\quad \hfill & 2\text{SC/dSC}\hfill \end{cases}.\end{equation} \tag{ 41 }$$

The remaining transition, between the NOR and dSC phases formally has a critical field as well,

$$\begin{equation}\frac{{H}_{\mathrm{c}}^{2}}{{\mu }^{4}/\lambda }=\frac{1-\frac{5}{12}\alpha }{1-\eta }\frac{3{g}^{2}+4{e}^{2}}{4{e}^{2}}\quad \text{dSC/NOR},\end{equation} \tag{ 42 }$$

but it turns out that this transition is never realized in the phase diagram.

The phase structure must of course be invariant under the rotation chosen for the gauge fields. As a check, in the massless limit α = 0 one recovers the results of reference [17], where a different rotation was used. We see that if the magnetic field is given in units of ${\mu }^{2}/\sqrt{\lambda }$, the parameter space has reduced to four dimensions, spanned by H, η, α, and g (for given e ≃ 0.30). In figure 1 we show two slices of this parameter space, namely in the η–H plane for the value of g used later, and in the g–H plane for the value of η used later. In each case we compare the massless result with the α = 0.3 result. For a typical value T_c/μ_q ∼ 0.1 this value of α corresponds to m_s/μ_q ∼ 0.5. We see that for nonzero α there are two triple points, i.e. points where three phases have the same free energy. Both triple points are realized in the left panel. The one where CFL, 2SC, and NOR phases meet, and which also exists in the massless case, is given by

$$\begin{equation}\left(\eta ,\frac{H}{{\mu }^{2}/\sqrt{\lambda }}\right)=\left(\frac{\frac{3}{8}\left(1+\frac{4}{3}\alpha \right){g}^{2}-\left(1-\frac{11}{12}\alpha \right){e}^{2}}{\left(1-\frac{\alpha }{6}\right)\left(3{g}^{2}+{e}^{2}\right)},\sqrt{\left(1-\frac{\alpha }{6}\right)\frac{3{g}^{2}+{e}^{2}}{2{e}^{2}}}\right).\end{equation} \tag{ 43 }$$

This point also occurs in the right panel, where it has the coordinates

$$\begin{equation}\left(g,\frac{H}{{\mu }^{2}/\sqrt{\lambda }}\right)=\left(\frac{2\sqrt{2}e}{\sqrt{3}}\sqrt{\frac{1-\frac{11}{12}\alpha +\eta \left(1-\frac{\alpha }{6}\right)}{1+\frac{4}{3}\alpha -8\eta \left(1-\frac{\alpha }{6}\right)}},\frac{3}{\sqrt{2}}\sqrt{\frac{1-\frac{5}{6}\alpha +\frac{{\alpha }^{2}}{9}}{1+\frac{4}{3}\alpha -8\eta \left(1-\frac{\alpha }{6}\right)}}\right).\end{equation} \tag{ 44 }$$

Zoom In Zoom Out Reset image size

The other triple point, where dSC, 2SC, and CFL phases meet, visible only in the left panel, is given by

$$\begin{equation}\left(\eta ,\frac{H}{{\mu }^{2}/\sqrt{\lambda }}\right)=\left(-\frac{6-7\alpha }{2\left(3+\alpha \right)},\frac{1}{6}\sqrt{\alpha \left(15+\frac{4{e}^{2}}{{g}^{2}}+\frac{9{g}^{2}}{{e}^{2}}\right)}\right).\end{equation} \tag{ 45 }$$

In the massless case, this point moves to the H = 0 axis and is no longer a triple point because the intermediate dSC phase disappears in this limit.

To get an idea of the strength of the magnetic field in physical units, let us assume μ_q ≃ 400 MeV, T_c ≃ 0.05μ_q, which is typical for the interior of compact stars. Then, for instance, the critical magnetic field between the NOR and 2SC phases, ${H}_{\mathrm{c}}\simeq 14\enspace {\mu }^{2}/\sqrt{\lambda }$, approximately corresponds to ⁴ H_c ≃ 1.1 × 10¹⁹ G.

The main effects of the strange quark mass, manifest in figure 1, are as follows. As expected and observed in many other calculations, the strange quark mass tends to disfavor CFL compared to 2SC, and it also tends to slightly disfavor the 2SC phase compared to the normal phase. The dSC phase only appears in the presence of a strange quark mass, and does so already for vanishing magnetic field [31]. Both panels show that for the values η = −1/2 and g = 3.5, which we shall use in our calculation of the flux tubes, the dSC phase plays no role. As already pointed out in the massless case [17], there is a regime of small g where the system transitions directly from the CFL phase to the NOR phase as we increase the magnetic field. Since our focus is on more realistic values of g applicable to compact stars, we deal with the scenario where the transition to the NOR phase occurs via the 2SC phase. This is still true with a nonzero strange quark mass.

Below a certain value of η, for instance η ≲ −0.85 if g = 3.5, the 2SC phase is the ground state even at vanishing magnetic field. This is worth mentioning since our main results concern 2SC flux tubes. For such a value of η, which goes beyond the weak-coupling results, the 2SC phase is bounded by the NOR phase for large magnetic fields, but not bounded by any other phase for low magnetic fields, and thus the region in which 2SC flux tubes can be expected is larger than for the weak-coupling value η = −1/2. In the weak-coupling scenario—extrapolated to g = 3.5—the 2SC phase is also bounded from below, namely by the CFL phase, and we shall see that this limits the region of 2SC flux tubes.

Before we turn to the flux tubes themselves, it is useful to compute their upper critical field H_c2. In the standard scenario of a single condensate, this is the maximum magnetic field which can sustain a nonzero condensate, under the assumption of a second-order transition to the normal phase. It is therefore the critical field below which an array of flux tubes is expected. In our multi-component system the situation is more complicated, and we have to calculate different critical fields H_c2 depending on which condensates melt. Having H_c and H_c2 at hand, we can then determine the parameter regime where the color superconductor is of type II, and in particular where we expect 2SC flux tubes.

The calculation of H_c2 is a generalization to nonzero strange quark mass of the analogous calculation done in reference [17]. That calculation, in turn, was a generalization of the standard single-component calculation which can be found in many textbooks. In a single-component superconductor, one linearizes the Ginzburg–Landau equations for a small condensate. The equation for the condensate then has the form of the Schrödinger equation for the harmonic oscillator, from which one reads off the maximal possible magnetic field H_c2, corresponding to the ground state energy. We know from the previous subsection that for strong coupling, as we decrease the magnetic field within the NOR phase, we encounter the 2SC phase. Consequently, for the corresponding critical field H_c2 we only have to take into account a single condensate, i.e. this case is analogous to the textbook scenario and leads to the simple generalization of equation (40) in reference [17],

$$\begin{equation}2\text{SC /NOR:}{H}_{\mathrm{c}2}=\frac{3\left({\mu }^{2}-{m}_{3}^{2}\right)}{e}=\frac{3{\mu }^{2}}{e}\left(1-\frac{\alpha }{12}\right).\end{equation} \tag{ 46 }$$

In this standard scenario all three critical magnetic fields H_c, H_c1, and H_c2 intersect at a single point (as a function of a model parameter, usually the Ginzburg–Landau parameter κ). Therefore, this intercept defines the transition between type-I and type-II behavior (usually at $\kappa =1/\sqrt{2}$). Here, by equating H_c2 with H_c for the 2SC/NOR transition in equation (41) we find for this transition point

$$\begin{equation}\frac{{T}_{\mathrm{c}}}{{\mu }_{\mathrm{q}}}=\frac{\sqrt{7\zeta \left(3\right)}}{12\sqrt{3}{\pi }^{2}}\sqrt{{g}^{2}+\frac{{e}^{2}}{3}}+\mathcal{O}\left(\frac{{m}_{\mathrm{s}}^{4}}{{\mu }_{\mathrm{q}}^{4}}\right),\end{equation} \tag{ 47 }$$

where the weak-coupling expression of λ in equation (15) has been used. It thus turns out that T_c/μ_q is a natural parameter to distinguish between type-I and type-II behavior—with large T_c/μ_q corresponding to type II. Interestingly, we see that there is no mass correction to the transition point within the order of our approximation. We shall see later that indeed H_c1 intersects H_c and H_c2 at the same T_c/μ_q. The reason is that in the vicinity of this point the system effectively behaves as a single-component system. Additional condensates can be induced in the cores of 2SC flux tubes—and our main results concern such unconventional flux tubes—but we will see that this is not the case close to the point (47).

The transition from the homogeneous 2SC phase, where ϕ₃ is nonzero, to an inhomogeneous phase is slightly more complicated. Assuming a second-order transition, we linearize the Ginzburg–Landau equations in ϕ₁ and ϕ₂. In the massless limit, the resulting two (decoupled) equations yield the same critical magnetic field [17]. In other words, as we approach the flux tube phase by decreasing H, both ϕ₁ and ϕ₂ become nonzero simultaneously (and continuously). This is different for nonzero m_s, in which case the two relevant equations are

$$\begin{equation}\left[{\left(\nabla +i{\tilde {q}}_{3}{\tilde {\mathbf{A}}}_{3}\right)}^{2}+\left({\mu }^{2}-{m}_{1}^{2}\right)+2h\vert {\phi }_{3}{\vert }^{2}\right]{\phi }_{1}\simeq 0,\end{equation} \tag{ 48a }$$

$$\begin{equation}\left[{\left(\nabla -i{\tilde {q}}_{3}{\tilde {\mathbf{A}}}_{3}\right)}^{2}+\left({\mu }^{2}-{m}_{2}^{2}\right)+2h\vert {\phi }_{3}{\vert }^{2}\right]{\phi }_{2}\simeq 0,\end{equation} \tag{ 48b }$$

where we have set ${\tilde {\mathbf{A}}}_{8}=0$ since ${\tilde {B}}_{8}=0$ in the 2SC phase and where ${\phi }_{3}={\rho }_{3}/\sqrt{2}$ is the condensate in the homogeneous 2SC phase (26). With the usual arguments and using the 2SC relation between ${\tilde {B}}_{3}$ and H from equation (27), we obtain two different critical fields,

$$\begin{equation}{H}_{\mathrm{c}2}^{\left(1\right)}=\frac{2\left(3{g}^{2}+{e}^{2}\right)}{3e{g}^{2}}\left[{\mu }^{2}\left(1+\eta \right)-\left({m}_{1}^{2}+\eta {m}_{3}^{2}\right)\right]=\frac{2{\mu }^{2}\left(3{g}^{2}+{e}^{2}\right)}{3e{g}^{2}}\left(1+\eta -\frac{4+\eta }{12}\alpha \right),\end{equation} \tag{ 49a }$$

$$\begin{equation}{H}_{\mathrm{c}2}^{\left(2\right)}=\frac{2\left(3{g}^{2}+{e}^{2}\right)}{3e{g}^{2}}\left[{\mu }^{2}\left(1+\eta \right)-\left({m}_{2}^{2}+\eta {m}_{3}^{2}\right)\right]=\frac{2{\mu }^{2}\left(3{g}^{2}+{e}^{2}\right)}{3e{g}^{2}}\left(1+\eta -\frac{7+\eta }{12}\alpha \right).\end{equation} \tag{ 49b }$$

The most relevant case for us is the one where both ${H}_{\mathrm{c}2}^{\left(1\right)}$ and ${H}_{\mathrm{c}2}^{\left(2\right)}$ are positive (a formally negative value indicates that the critical field does not exist, indicating that the homogeneous phase persists down to H = 0). This is the case for η = −1/2 and all reasonable, i.e. not too large, values of α. In this scenario, there is a transition at ${H}_{\mathrm{c}2}^{\left(1\right)}$ from the homogeneous 2SC phase to a phase where both ϕ₁ and ϕ₃ are nonzero, which is an inhomogeneous version of the dSC phase (see section 3.3). Then, as we reach the 'would-be' ${H}_{\mathrm{c}2}^{\left(2\right)}$ by further decreasing H, the approximation by which this critical field was computed is no longer valid, and thus the value for ${H}_{\mathrm{c}2}^{\left(2\right)}$ becomes irrelevant. Nevertheless, it can be expected that there will be some transition from an inhomogeneous dSC phase to an inhomogeneous CFL phase. The existence of an intermediate inhomogeneous dSC phase due to the nonzero strange quark mass is an interesting new observation, but it is beyond the scope of this paper to construct this phase explicitly.

We may again compute the transition point between type-I and type-II behavior. For η = −1/2 (where there is no homogeneous dSC phase), we equate ${H}_{\mathrm{c}2}^{\left(1\right)}$ with H_c for the 2SC/CFL transition from equation (41). Dropping terms quadratic in α we find

$$\begin{equation}\frac{{T}_{\mathrm{c}}}{{\mu }_{\mathrm{q}}}\simeq c\left(1-\frac{\alpha }{4}\right),\quad c\equiv \frac{\sqrt{7\zeta \left(3\right)}}{12\sqrt{2}{\pi }^{2}}\frac{g\sqrt{3{g}^{2}+4{e}^{2}}}{\sqrt{3{g}^{2}+{e}^{2}}}.\end{equation} \tag{ 50 }$$

Since α depends on T_c/μ_q this is an implicit equation for T_c/μ_q. To lowest nontrivial order in ${m}_{\mathrm{s}}^{2}/{\mu }_{\mathrm{q}}^{2}$ the solution is

$$\begin{equation}\frac{{T}_{\mathrm{c}}}{{\mu }_{\mathrm{q}}}=c\left(1+\frac{{m}_{\mathrm{s}}^{2}}{8{\mu }_{\mathrm{q}}^{2}}\enspace \mathrm{ln}\enspace c\right)+\mathcal{O}\left(\frac{{m}_{\mathrm{s}}^{4}}{{\mu }_{\mathrm{q}}^{4}}\right).\end{equation} \tag{ 51 }$$

Therefore, this transition point between type-I and type-II behavior does receive a correction quadratic in m_s (linear in α), in contrast to the transition point (47). The detailed phase structure around this point is expected to be complicated. This is due to the intermediate inhomogeneous dSC phase, as just discussed, but even without mass correction this transition point is affected in a nontrivial way by the multi-component nature of the system [17, 21]. Most importantly, if the lower boundary of the flux tube region H_c1 is computed in the usual way, i.e. assuming a second-order transition, it turns out that the three critical fields no longer intersect in a single point, and the situation becomes more complicated due to a first-order entrance into the flux tube phase. Here we do not have to deal with these complications, since the precise location of the transition point (51) and the phase transitions in its vicinity are not relevant for the 2SC flux tubes.

Having identified the parameter range where type-II behavior with respect to 2SC flux tubes is expected, we can now turn to the explicit construction of these flux tubes. We will restrict ourselves to the calculation of an isolated, straight flux tube, such that we can employ cylindrical symmetry and our calculation becomes effectively one-dimensional in the radial direction. This is sufficient to compute the critical field H_c1, which is defined as the field at which it becomes favorable to place a single flux tube in the system, indicating a second-order transition to a phase containing an array of flux tubes. Since the distance between the flux tubes goes to infinity as H_c1 is approached from above, the interaction between flux tubes plays no role. As explained in the introduction, our main goal is to determine the fate of the 2SC domain walls in the presence of a nonzero strange quark mass. Therefore, we focus exclusively on 2SC flux tubes, i.e. configurations which asymptote to the 2SC phase far away from the center of the flux tube.

In order to compute the profiles of the condensates and the gauge fields we need to derive their equations of motion and bring them into a form convenient for the numerical evaluation. We work in cylindrical coordinates (r, θ, z), where, as above, the z-axis is aligned with the external magnetic field H and thus with the flux tube. We introduce dimensionless condensates f_i (i = 1, 2, 3), which only depend on the radial distance to the center of the flux tube,

$$\begin{equation}{\rho }_{i}\left(\mathbf{r}\right)={f}_{i}\left(r\right){\rho }_{2\text{SC}},\end{equation} \tag{ 52 }$$

where we have denoted the condensate of the homogeneous 2SC phase (26) by ρ_2SC. Since we are interested in 2SC flux tubes, we impose the boundary conditions f₃(∞) = 1 and f₁(∞) = f₂(∞) = 0. As in ordinary single-component flux tubes, we allow for a nonzero winding number $n\in \mathbb{Z}$, such that the phases of the condensates are

$$\begin{equation}{\psi }_{1}\left(\mathbf{r}\right)={\psi }_{2}\left(\mathbf{r}\right)=0,\quad {\psi }_{3}\left(\mathbf{r}\right)=n\theta .\end{equation} \tag{ 53 }$$

Here we have set the winding numbers for the 'non-2SC' condensates f₁ and f₂ to zero. In principle, we might include configurations where these windings are nonzero. (The baryon circulation around the flux tube vanishes for arbitrary choices of the winding numbers as long as f₁(∞) = f₂(∞) = 0.) In such configurations, f₁ and/or f₂ would have to vanish far away from the flux tube and in the center of the flux tube, i.e. at best they would be non-vanishing in an intermediate domain. These configurations do not play a role in the massless limit [17] and there is no obvious reason why they should become important if a strange quark mass is taken into account. Therefore, we shall work with equation (53). As a consequence, the boundary condition for the 2SC condensate in the core is f₃(0) = 0, while f₁(0) and f₂(0) can be nonzero and must be determined dynamically.

As we have seen in section 2, after the rotation of the gauge fields $\tilde {\mathbf{A}}$ decouples from the condensates. Therefore, we are left with two nontrivial gauge fields, for which we introduce the dimensionless versions ${\tilde {a}}_{3}$ and ${\tilde {a}}_{8}$ via

$$\begin{equation}{\tilde {\mathbf{A}}}_{3}\left(\mathbf{r}\right)=\left[\frac{H\enspace \mathrm{cos}\enspace {\vartheta }_{1}\enspace \mathrm{sin}\enspace {\vartheta }_{2}}{2}r+\frac{{\tilde {a}}_{3}\left(r\right)}{r}\right]{\mathbf{e}}_{\theta },\quad {\tilde {\mathbf{A}}}_{8}\left(\mathbf{r}\right)=\frac{{\tilde {a}}_{8}\left(r\right)}{r}{\mathbf{e}}_{\theta },\end{equation} \tag{ 54 }$$

with the boundary conditions ${\tilde {a}}_{3}\left(0\right)={\tilde {a}}_{8}\left(0\right)=0$. This yields the magnetic fields

$$\begin{equation}{\tilde {\mathbf{B}}}_{3}=\left(H\enspace \mathrm{cos}\enspace {\vartheta }_{1}\enspace \mathrm{sin}\enspace {\vartheta }_{2}+\lambda {\rho }_{2\text{SC}}^{2}\frac{{\tilde {a}}_{3}^{\prime }}{R}\right){\mathbf{e}}_{z},\quad {\tilde {\mathbf{B}}}_{8}=\lambda {\rho }_{2\text{SC}}^{2}\frac{{\tilde {a}}_{8}^{\prime }}{R}{\mathbf{e}}_{z},\end{equation} \tag{ 55 }$$

where prime denotes derivative with respect to R, which is the dimensionless radial coordinate

$$\begin{equation}R=\frac{r}{{\xi }_{3}},\quad {\xi }_{3}\equiv \frac{1}{\sqrt{\lambda }{\rho }_{2\text{SC}}}.\end{equation} \tag{ 56 }$$

Here, ξ₃ is the coherence length; we have added a subscript 3 to indicate that in a 2SC flux tube it is the ud condensate f₃ whose asymptotic behavior is characterized by the coherence length. Following reference [17], we have separated a term in ${\tilde {\mathbf{A}}}_{3}$ for convenience, which gives rise to the nonzero field ${\tilde {\mathbf{B}}}_{3}$ in the 2SC phase, i.e. far away from the flux tube. Therefore, the dimensionless gauge fields do not create any additional magnetic fields at infinity, ${\tilde {a}}_{3}^{\prime }\left(\infty \right)={\tilde {a}}_{8}^{\prime }\left(\infty \right)=0$ (alternatively, this effect could have been implemented in the boundary condition for ${\tilde {a}}_{3}$). The behavior of ${\tilde {B}}_{3}$ is another qualitative difference of the 2SC flux tube to a textbook flux tube (besides potentially induced additional condensates and the existence of two gauge fields). In a standard flux tube, the magnetic field is expelled in the superconducting phase far away from the flux tube and penetrates through the normal-conducting centre. Here, we have three magnetic fields: $\tilde {B}$, which fully penetrates the superconductor and thus is irrelevant for the calculation of the flux tube profiles; ${\tilde {B}}_{8}$, which behaves analogously to the ordinary magnetic field in an ordinary flux tube; and ${\tilde {B}}_{3}$, which is nonzero far away from the flux tube and is affected nontrivially by the flux tube profile. As a consequence, since ${\tilde {B}}_{3}$ depends on the external field H, the flux tube profiles and flux tube energies also depend on H, which poses a technical complication. We should keep in mind, however, that it is the ordinary magnetic field B that dictates the formation of magnetic defects. A flux tube configuration, in which the condensation energy is necessarily reduced, may become favored if this energy cost is overcompensated by admitting magnetic B-flux in the system (this is the meaning of the −B ⋅ H term in the Gibbs free energy). It is therefore useful for the interpretation of our results to compute the unrotated field B from the profiles. Undoing the rotation (6) and using equations (13) and (55), we find

$$\begin{equation}\frac{B}{H}={\mathrm{cos}}^{2}\enspace {\vartheta }_{1}\left[1+\frac{e\enspace \mathrm{sin}\enspace {\vartheta }_{1}}{4{\Xi}}\left(\frac{{\tilde {a}}_{8}^{\prime }}{R}+\frac{3g}{\sqrt{3{g}^{2}+4{e}^{2}}}\frac{{\tilde {a}}_{3}^{\prime }}{R}\right)\right],\end{equation} \tag{ 57 }$$

where we have introduced the dimensionless magnetic field

$$\begin{equation}{\Xi}\equiv \frac{{\tilde {q}}_{3}H\enspace \mathrm{cos}\enspace {\vartheta }_{1}\enspace \mathrm{sin}\enspace {\vartheta }_{2}}{2\lambda {\rho }_{2\text{SC}}^{2}}=\frac{3e{g}^{2}}{4\sqrt{\lambda }\left(3{g}^{2}+{e}^{2}\right)}{\left(1-\frac{{m}_{3}^{2}}{{\mu }^{2}}\right)}^{-1}\frac{H}{{\mu }^{2}/\sqrt{\lambda }}.\end{equation} \tag{ 58 }$$

We can now express the potential U₀ (10) in terms of the dimensionless condensates and gauge fields,

$$\begin{align}\hfill {U}_{0}& ={U}_{2\text{SC}}+\frac{\lambda {\rho }_{2\text{SC}}^{4}}{2}\left[{f}_{1}^{\prime 2}+{f}_{2}^{\prime 2}+{f}_{3}^{\prime 2}+{f}_{1}^{2}\left(\frac{{f}_{1}^{2}}{2}-\frac{{\mu }^{2}-{m}_{1}^{2}}{{\mu }^{2}-{m}_{3}^{2}}\right)+{f}_{2}^{2}\left(\frac{{f}_{2}^{2}}{2}-\frac{{\mu }^{2}-{m}_{2}^{2}}{{\mu }^{2}-{m}_{3}^{2}}\right)\right.\hfill \\ \hfill & \quad \left.+\frac{{\left(1-{f}_{3}^{2}\right)}^{2}}{2}+\frac{{\left({\mathcal{N}}_{1}+{\Xi}{R}^{2}\right)}^{2}{f}_{1}^{2}+{\left({\mathcal{N}}_{2}-{\Xi}{R}^{2}\right)}^{2}{f}_{2}^{2}+{\mathcal{N}}_{3}^{2}{f}_{3}^{2}}{{R}^{2}}\right.\hfill \\ \hfill & \quad \left.-\eta \left({f}_{1}^{2}{f}_{2}^{2}+{f}_{1}^{2}{f}_{3}^{2}+{f}_{2}^{2}{f}_{3}^{2}\right)\right].\hfill \end{align} \tag{ 59 }$$

Here we have denoted the potential of the homogeneous 2SC phase by

$$\begin{equation}{U}_{2\text{SC}}=-\frac{{\left({\mu }^{2}-{m}_{3}^{2}\right)}^{2}}{4\lambda },\end{equation} \tag{ 60 }$$

and we have abbreviated

$$\begin{equation}{\mathcal{N}}_{1}\equiv {\tilde {q}}_{3}{\tilde {a}}_{3}+{\tilde {q}}_{81}{\tilde {a}}_{8},\quad {\mathcal{N}}_{2}\equiv -{\tilde {q}}_{3}{\tilde {a}}_{3}+{\tilde {q}}_{82}{\tilde {a}}_{8},\quad {\mathcal{N}}_{3}\equiv n-{\tilde {q}}_{83}{\tilde {a}}_{8}.\end{equation} \tag{ 61 }$$

Inserting equation (59) into the Gibbs free energy density (14) and using the expressions for the magnetic fields (55), we derive the equations of motion for the gauge fields,

$$\begin{equation}{\tilde {a}}_{3}^{{\prime\prime}}-\frac{{\tilde {a}}_{3}^{\prime }}{R}=\frac{{\tilde {q}}_{3}}{\lambda }\left[\left({\mathcal{N}}_{1}+{\Xi}{R}^{2}\right){f}_{1}^{2}-\left({\mathcal{N}}_{2}-{\Xi}{R}^{2}\right){f}_{2}^{2}\right],\end{equation} \tag{ 62a }$$

$$\begin{equation}{\tilde {a}}_{8}^{{\prime\prime}}-\frac{{\tilde {a}}_{8}^{\prime }}{R}=\frac{1}{\lambda }\left[{\tilde {q}}_{81}\left({\mathcal{N}}_{1}+{\Xi}{R}^{2}\right){f}_{1}^{2}+{\tilde {q}}_{82}\left({\mathcal{N}}_{2}-{\Xi}{R}^{2}\right){f}_{2}^{2}-{\tilde {q}}_{83}{\mathcal{N}}_{3}{f}_{3}^{2}\right],\end{equation} \tag{ 62b }$$

and for the condensates,

$$\begin{equation}0={f}_{1}^{{\prime\prime}}+\frac{{f}_{1}^{\prime }}{R}+{f}_{1}\left[\frac{{\mu }^{2}-{m}_{1}^{2}}{{\mu }^{2}-{m}_{3}^{2}}-{f}_{1}^{2}-\frac{{\left({\mathcal{N}}_{1}+{\Xi}{R}^{2}\right)}^{2}}{{R}^{2}}+\eta \left({f}_{2}^{2}+{f}_{3}^{2}\right)\right],\end{equation} \tag{ 63a }$$

$$\begin{equation}0={f}_{2}^{{\prime\prime}}+\frac{{f}_{2}^{\prime }}{R}+{f}_{2}\left[\frac{{\mu }^{2}-{m}_{2}^{2}}{{\mu }^{2}-{m}_{3}^{2}}-{f}_{2}^{2}-\frac{{\left({\mathcal{N}}_{2}-{\Xi}{R}^{2}\right)}^{2}}{{R}^{2}}+\eta \left({f}_{1}^{2}+{f}_{3}^{2}\right)\right],\end{equation} \tag{ 63b }$$

$$\begin{equation}0={f}_{3}^{{\prime\prime}}+\frac{{f}_{3}^{\prime }}{R}+{f}_{3}\left[1-{f}_{3}^{2}-\frac{{\mathcal{N}}_{3}^{2}}{{R}^{2}}+\eta \left({f}_{1}^{2}+{f}_{2}^{2}\right)\right].\end{equation} \tag{ 63c }$$

By taking the limit R → ∞ of equation (62b) we conclude

$$\begin{equation}{\tilde {a}}_{8}\left(\infty \right)=\frac{n}{{\tilde {q}}_{83}},\end{equation} \tag{ 64 }$$

while no condition for ${\tilde {a}}_{3}\left(\infty \right)$ can be derived, hence this value has to be determined dynamically. Due to the boundary value (64), the baryon circulation around the flux tube vanishes, as in a standard magnetic flux tube. We discuss the asymptotic behavior far away from the center of the flux tube in appendix A.

The Gibbs free energy density can be written as

$$\begin{equation}G={U}_{2\text{SC}}-\frac{{H}^{2}\enspace {\mathrm{cos}}^{2}\enspace {\vartheta }_{1}}{2}+\frac{L}{V}\left(\mathcal{F}-{\tilde {{\Phi}}}_{8}H\enspace \mathrm{sin}\enspace {\vartheta }_{1}\right),\end{equation} \tag{ 65 }$$

where L is the size of the system in the z-direction, and

$$\begin{equation}{\tilde {{\Phi}}}_{8}=\oint \mathrm{d}\mathbf{s}\cdot {\tilde {\mathbf{A}}}_{8}=\frac{2\pi n}{{\tilde {q}}_{83}},\end{equation} \tag{ 66 }$$

with a closed integration contour encircling the flux tube at infinity, is the magnetic ${\tilde {\mathbf{B}}}_{8}$-flux through the flux tube. Employing partial integration and the equations of motion (63), the free energy of a single flux tube per unit length is $\mathcal{F}=\pi {\rho }_{2\text{SC}}^{2}\mathcal{I}$ with

$$\begin{equation}\mathcal{I}\equiv {\int }_{0}^{\infty }\mathrm{d}R\enspace R\left[\frac{\lambda \left({\tilde {a}}_{3}^{\prime 2}+{\tilde {a}}_{8}^{\prime 2}\right)}{{R}^{2}}-\frac{{f}_{1}^{4}}{2}-\frac{{f}_{2}^{4}}{2}+\frac{1-{f}_{3}^{4}}{2}+\eta \left({f}_{1}^{2}{f}_{2}^{2}+{f}_{1}^{2}{f}_{3}^{2}+{f}_{2}^{2}{f}_{3}^{2}\right)\right].\end{equation} \tag{ 67 }$$

Written in this form, the free energy does not have any explicit dependence on the mass correction and is identical to the one in reference [17] (of course, the dependence on the mass enters implicitly through the equations of motion). The critical field H_c1 is defined as the field above which G is lowered by the addition of a flux tube, i.e. by the point at which the term in parentheses in equation (65) is zero. This condition can be written as

$$\begin{equation}{{\Xi}}_{\mathrm{c}1}=\frac{{g}^{2}\mathcal{I}\left({{\Xi}}_{\mathrm{c}1},n\right)}{8\lambda n},\end{equation} \tag{ 68 }$$

where Ξ_c1 is the dimensionless version of H_c1 via equation (58). In the ordinary textbook scenario, the free energy of a flux tube does not depend on the external magnetic field, and thus (68) would be an explicit expression for the (dimensionless) critical magnetic field. The free energy of a 2SC flux tube, however, does depend on the external magnetic field. Therefore, equation (68) is an implicit equation for Ξ_c1, which has to be solved numerically.

As we shall see in the next section, H_c1 may intersect with ${H}_{\mathrm{c}2}^{\left(1\right)}$, which indicates the second-order transition from the homogeneous 2SC phase to an inhomogeneous phase where the condensate f₁ is switched on. Since H_c1 becomes meaningless beyond this intercept, it is useful to compute this point explicitly. To this end, we write ${H}_{\mathrm{c}2}^{\left(1\right)}$ from equation (49a) (setting η = −1/2) in terms of the dimensionless field Ξ (58), which yields to linear order in α

$$\begin{equation}{{\Xi}}_{\mathrm{c}2}^{\left(1\right)}\simeq \frac{1}{4}\left(1-\frac{\alpha }{2}\right).\end{equation} \tag{ 69 }$$

Solving this equation simultaneously with equation (68), i.e. setting ${{\Xi}}_{\mathrm{c}2}^{\left(1\right)}={{\Xi}}_{\mathrm{c}1}$, gives the intersection point. In the practical calculation, this is best done by inserting equation (69) into equation (68) and solving the resulting equation for T_c/μ_q (for given strange quark mass and winding number). This calculation is relevant for the phase diagrams in figure 6.

We start by discussing individual flux tube profiles and the associated magnetic fields. We do so by choosing a fixed T_c/μ_q such that for vanishing strange quark mass there is a magnetic field at which it is energetically favorable to put a domain wall in the system. The values of T_c/μ_q for which this is the case are known from reference [17] (see figure 5 in that reference). The domain wall interpolates between the 2SC_ud and 2SC_us phases, i.e. on one side, and far away from it, we have f₃ = 1, f₂ = 0, and on the other side f₂ = 1, f₃ = 0. At nonzero m_s the free energies of 2SC_ud and 2SC_us are obviously no longer equal, as we have seen explicitly in section 3, and thus the domain wall configuration becomes unstable. In figure 2 we have chosen 10 different nonzero values of m_s/μ_q, such that the configuration with lowest H_c1 (to which we will refer as the 'energetically preferred' or simply the 'preferred' configuration) has winding number n = 10, ..., 1 as m_s/μ_q is increased. The profiles and magnetic fields are plotted at the corresponding minimal H_c1. The critical fields as a function of the winding number for the parameters of figure 2 are shown in figure 3, which proves the successive decrease in the preferred winding from low to high strange quark mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let us first discuss the profiles themselves in figure 2. For the smallest masses shown here the profiles of the condensates are reminiscent of a domain wall profile: the second condensate, which is induced in the core, assumes essentially the value of the homogeneous 2SC condensate. Of course, in contrast to a domain wall, the flux tubes have a finite radius, which decreases as the winding number decreases (with increasing strange quark mass). For a quantitative discussion of the size of the flux tubes we use the various characteristic length scales obtained from the asymptotic behavior of the flux tube profile, see appendix A. The penetration depth ℓ usually is characteristic for the decay of the magnetic field away from the center of the flux tube. In our case this scale corresponds to the decay of the ${\tilde {B}}_{8}$ field (in figure 2 we have restricted ourselves to plotting the ordinary magnetic field B). The coherence length ξ₃ introduced in equation (56) indicates the scale on which the ud condensate f₃ recovers its homogeneous value and corresponds to the usual coherence length of a single-component flux tube. Finally, and in contrast to the standard single-component scenario, our calculation in appendix A shows that there is a third scale ξ₂, characteristic for the induced us condensate f₂. Using equations (A8)–(A10) we can write these three length scales in physical units as

$$\begin{equation}\ell \simeq 0.94\left[1+\frac{1}{48}{\left(\frac{{m}_{\mathrm{s}}}{{\mu }_{\mathrm{q}}}\right)}^{2}\enspace \mathrm{ln}\enspace \frac{{\mu }_{\mathrm{q}}}{{T}_{\mathrm{c}}}\right]{\left(\frac{{\mu }_{\mathrm{q}}}{400\enspace \text{MeV}}\right)}^{-1}\enspace \text{fm},\end{equation} \tag{ 70a }$$

$$\begin{equation}{\xi }_{3}\simeq \frac{0.066}{{T}_{\mathrm{c}}/{\mu }_{\mathrm{q}}}\left[1+\frac{1}{48}{\left(\frac{{m}_{\mathrm{s}}}{{\mu }_{\mathrm{q}}}\right)}^{2}\enspace \mathrm{ln}\enspace \frac{{\mu }_{\mathrm{q}}}{{T}_{\mathrm{c}}}\right]{\left(\frac{{\mu }_{\mathrm{q}}}{400\enspace \text{MeV}}\right)}^{-1}\enspace \text{fm},\end{equation} \tag{ 70b }$$

$$\begin{equation}{\xi }_{2}\simeq 1.3{\left(\frac{H}{{\mu }^{2}/\sqrt{\lambda }}\frac{{T}_{\mathrm{c}}}{{\mu }_{\mathrm{q}}}\right)}^{-1/2}{\left(\frac{{\mu }_{\mathrm{q}}}{400\enspace \text{MeV}}\right)}^{-1}\enspace \text{fm},\end{equation} \tag{ 70c }$$

where we have used equation (15) with T = 0, the 2SC condensate (26), and g = 3.5. The length scale ξ₂ on which the induced condensate varies depends on the external magnetic field H, for which the dimensionless value (in units of ${\mu }^{2}/\sqrt{\lambda }$) can be inserted into equation (70c). For the parameters used in figure 2 the mass terms in equations (70a) and (70b) are negligibly small and setting μ_q = 400  MeV we find ℓ ≃ 0.94 fm, ξ₃ ≃ 0.77 fm, and ξ₂ ≃ (1.4–1.5) fm, while the maximal value of the radial coordinate in this figure R = 10 corresponds to r ≃ 7.7 fm. These values, characterizing the radial extent of the core of the flux tube do not depend on the winding number n. In the derivation in appendix A we have worked with a fixed and finite n while sending the radial coordinate to infinity. However, from the full numerical profiles we see that the core of the flux tube clearly grows with increasing n. For instance, taking the point where the two condensates have the same value as a rough measure for the size of the core, we see that for the largest winding shown here, n = 10, corresponding to the smallest mass, m_s = 0.046μ_q, this size is r ≃ 4.6 fm. This is about six times larger than the coherence length ξ₃, but it also shows that it only takes a relatively small strange quark mass to reduce the size of the flux tubes from infinity (domain wall) to a few fm.

The plots in the upper panels also show the profile of the function $\sqrt{{f}_{2}^{2}+{f}_{3}^{2}}$, which is the radial coordinate of the space spanned by (f₂, f₃), i.e. a domain wall can be understood as a rotation in this space from the horizontal to the vertical axis. This function illustrates the depletion of the 'combined' condensate in a cylindrical layer at nonzero radius for large winding numbers and in the center of the tube for low winding numbers. This depletion reflects the structure of the magnetic field, which is shown in the lower panels. We recall that the 2SC phase (just like the CFL phase) admits a fraction of the external field H even in the homogeneous phase, and this fraction is close to one for g ≫ e. Therefore, the energetic benefit of a 2SC flux tube is to admit additional magnetic flux on top of the flux already present. As a consequence of our choice g = 3.5, the ratio B/H is larger than 99% already in the homogeneous 2SC phase, as one can see for instance from equation (57). The magnetic field profiles in the figure, showing B/H between its minimal homogeneous 2SC value (white) and the maximal value (black) turn out to be in the range 0.9975 ≲ B/H ≲ 0.9995. Therefore, the magnetic field variations are unlikely to have any physical impact in the strong coupling regime. One might also wonder about the numerical accuracy needed to resolve such tiny variations of B/H. However, from equation (57) it is obvious that it is the small mixing angle ϑ₁ (sin ϑ₁ ≃ 0.05 for g = 3.5) that maps the variations into a tiny interval and thus the required numerical accuracy is not as high as one might think.

The magnetic profiles show interesting features. We see that the flux tubes with higher winding (small m_s) have their excess magnetic field concentrated in ring-like structures. This effect has been observed in the literature before with the help of a two-component abelian Higgs model [44] and a one-component model with non-standard coupling between the condensate and the abelian gauge field [46] (see also references [58, 59] for similar structures of magnetic monopoles in a non-abelian model). In these studies, the so-called Bogomolny limit was considered for simplicity, which corresponds to the transition point between type-I and type-II behavior. Here we observe the effect in a numerical evaluation of the general Ginzburg–Landau equations, and in particular we can point out regions in the phase diagram where the exotic ring-like structures are the preferred configuration (for a systematic discussion of the phase diagram see next subsection). The ring-like structure is easy to understand: in the domain wall (m_s = 0) the excess magnetic field is concentrated in the wall for symmetry reasons. As the strange quark mass is increased, the wall 'bends' to form a flux tube with finite radius, and it is obvious for continuity reasons that for very small masses the maximum of the magnetic field is still sitting in the transition region. Then, as the sequence in figure 2 shows, the magnetic rings gradually turn into ordinary flux tubes as m_s is increased (and the winding number decreases). Furthermore, we observe a double-ring structure, clearly visible in the multi-winding flux tubes. We have found that the double ring does not always occur. For instance, for the smaller value T_c = 0.08μ_q the double ring is replaced by a single ring.

All profiles in figure 2 are 'unconventional' in the sense that they contain two condensates. As already mentioned above, these two-component solutions are not found everywhere in the parameter space. In the left panel of figure 4 we show the value of the induced condensate in the core f₂(0) as a function of the parameter T_c/μ_q for winding numbers n = 1, ..., 6 and for a nonzero strange quark mass, compared to the massless case. We see that as T_c/μ_q decreases, the induced condensate becomes smaller until it continuously goes to zero at a point that depends on the winding number. Interestingly, as the winding number is increased, this point seems to converge to the point (47) that distinguishes type-I from type-II superconductivity. In particular, the behavior of f₂(0) becomes more and more step-like for larger windings, i.e. f₂(0) is close to one until it sharply decreases near the type-I/type-II transition point. However, flux tubes with higher winding are energetically disfavored in the vicinity of the type-I/type-II transition point (as we shall see below), such that this interesting behavior does not seem to be physically relevant.

Zoom In Zoom Out Reset image size

The right panel of figure 4 is useful to understand why the induced condensate only appears in a certain parameter regime. In this plot we show the ratio of ξ₂ and ξ₃ given in equation (70). We see that the induced condensate vanishes if this ratio becomes too small. For a possible interpretation note that the coherence length is connected to the size of the Cooper pair. If we view ξ₃ as a typical size of the Cooper pair (also for us pairing), then the right panel of figure 4 suggests that the length scale for the variation of the induced condensate, which is determined by the external magnetic field, must be sufficiently large (at least by a factor of about 1.5) compared to the size of a Cooper pair. Otherwise the condensate 'does not fit' into the core of the flux tube.

We have seen that there are parameter choices for (g, T_c/μ_q, m_s/μ_q) where multi-winding 2SC flux tubes with a ring-like structure of the magnetic field are energetically preferred. In this section we investigate the parameter space more systematically. We should keep in mind that in QCD the parameters, g, T_c, and m_s are uniquely given by μ_q (at T = 0). Therefore, as we vary the quark chemical potential, the system will move along a unique, but unknown, curve in this three-dimensional parameter space. For instance, at asymptotically large μ_q, we start from small g ≪ 1, exponentially suppressed T_c and negligibly small m_s (compared to μ_q). Since we do not know the values of m_s and T_c at more moderate densities we keep our parameters as general as possible. In this sense, keeping g = 3.5 fixed for our results is a simplification, in an even more general calculation one might also vary g.

In figure 5 we compare the critical magnetic fields of flux tubes with different winding numbers for two different values of m_s as a function of T_c/μ_q. We plot the difference of the critical fields to the critical field of the n = 1 configuration because the different curves would be barely distinguishable had we plotted the critical fields themselves. It is instructive to start with the massless case (left panel). For T_c/μ_q ≃ 0.05 we observe the standard behavior at the type-I/type-II transition: in the type-I regime, the critical field H_c1 can be lowered by increasing the winding number, and as n → ∞ one expects H_c1(n) to converge to H_c, the critical field at which the NOR and 2SC phases coexist. Just above the critical value T_c/μ_q ≃ 0.05 the order of the different H_c1(n) is exactly reversed, and H_c1 is given by the flux tube with lowest winding, n = 1, the higher-winding flux tubes are disfavored. It is a good check for our numerics that all curves intersect at the same point, and this point is given by equation (47), whose derivation is completely independent of the numerical evaluation. This conventional behavior in the vicinity of the type-I/type-II transition point is expected because the second condensate plays no role here, at least for the lowest winding numbers, as we have seen in figure 4.

Zoom In Zoom Out Reset image size

Now, as T_c/μ_q is increased, the unconventional 2SC behavior becomes apparent. For any winding n > 1 there is a critical T_c/μ_q at which this higher winding becomes favored over n = 1. Similar to the type-I regime, there is a region in which H_c1(n) decreases monotonically with n until, for n → ∞, it is expected to intersect H_c1(1) at the critical point where domain walls set in (we have only plotted the curves for a few winding numbers, obviously the numerics become more challenging for large n). This point, T_c/μ_q ≃ 0.068, is taken from reference [17], where it was computed explicitly by using a planar instead of a cylindrical geometry (see figures 4 and 5 in that reference). Here we have marked the onset of domain walls by a vertical dashed line.

In the right panel we first note that the behavior in the vicinity of the type-I/type-II transition in the presence of a strange quark mass is qualitatively the same as in the massless case. The transition point itself, as we already know from equation (47), is even exactly the same, at least up to lowest nontrivial order in m_s. Now, however, the behavior for larger T_c/μ_q, where the second condensate does play a role, is qualitatively different. For the chosen value of m_s/μ_q, we find that for windings n ⩽ 6 there is a regime in which flux tubes with winding n are preferred, while flux tubes with n > 6 are never preferred. We see that the lowest H_c1 corresponds to flux tubes with winding numbers 1, 6, 5, 4, 3, 2, as T_c/μ_q is increased. These structures can be viewed as 'remnants' of the domain wall.

With the help of the results of this panel we can construct the phase diagram in the H–T_c/μ_q plane for this particular mass m_s = 0.15μ_q. To this end, we need to bring together the critical fields H_c1 from figure 5 with the critical fields H_c from section 4.1 and the critical fields H_c2 from section 4.2. The result is shown in the left panel of figure 6. This panel contains the first-order transition between NOR and 2SC phases (41), the second-order transitions from the NOR to the '2SC flux tube' phase (46) and from the 2SC to the 'CFL flux tube' phase (49a), and the second-order transition from the 2SC to the '2SC flux tube' phase just discussed, including the different segments corresponding to different winding numbers. We have included the phase transitions of the massless case for comparison, including the segment where domain walls are formed. The single-component flux tube configuration, i.e. the one without a us condensate in the core, is denoted by S₁, following the notation of reference [17]. This configuration always has winding number 1. The region labeled by 'CFL flux tubes' is bounded from above by ${H}_{\mathrm{c}2}^{\left(1\right)}$, which actually indicates a transition to dSC flux tubes (see discussion below equation (49)). For the given parameters, the 'would-be' critical magnetic field ${H}_{\mathrm{c}2}^{\left(2\right)}$ is very close to ${H}_{\mathrm{c}2}^{\left(1\right)}$, and thus we have, slightly inaccurately, termed the entire region 'CFL flux tubes', although we know that there is at least a thin slice where the inhomogeneous phase is not made of CFL flux tubes.

Zoom In Zoom Out Reset image size

Since our calculation only allows us to determine the critical fields for the entrance into the flux tube phases, we can only speculate about the structure of these inhomogeneous phases away from the second-order lines. It is obvious that the structure will be more complicated than in a standard single-component superconductor. This, firstly, concerns the transition between the phases labeled by 'CFL flux tubes' and '2SC flux tubes'. We have continued the H_c2 transition curve into the flux tube regime as a thin dotted line, but of course this curve should not be taken too seriously because it was calculated under the assumption of a homogeneous 2SC phase on one side of the transition. Secondly, and more closely related to our main results, the different winding numbers that occur upon entering the 2SC flux tube phase also suggest a complicated structure of the flux tube arrays. It is conceivable that there are transitions between 'pure' arrays of a single winding number or that there are arrays composed of flux tubes with different winding numbers.

We may now ask up to which values of the strange quark mass the multi-winding flux tubes survive. More precisely, at which value of m_s do only n = 1 flux tubes form at the entire phase transition curve H_c1 between the homogeneous 2SC phase and the '2SC flux tube' phase? This question can be answered by repeating the calculation needed for left panel of figure 6 for different values of m_s and determining the sequence of the preferred configurations as a function of T_c/μ_q in each case. The result is shown in in the right panel of figure 6. In this panel, the m_s/μ_q–T_c/μ_q plane is divided into different regions, each region labeled by the preferred flux tube configuration. One can check that a slice through this plot at m_s = 0.15μ_q reproduces exactly the sequence of phases shown in the left panel. In particular, we have indicated the transition between the standard flux tubes S₁ and the two-component flux tubes with winding number n = 1 by a dashed curve. We have restricted our calculation to winding numbers n ⩽ 10, but it is obvious how the pattern continues: as we move towards m_s = 0, say at fixed T_c/μ_q, the winding number of the preferred configuration increases successively and more rapidly. Eventually, the winding number and thickness of the flux tube go to infinity, which corresponds to the domain wall, whose range is indicated with (red) arrows. Therefore, the lines that bound the region of multi-winding configurations at small m_s from both sides—here calculated from n = 10—will slightly change if higher windings are taken into account and they are expected to converge to the arrows at m_s = 0. In the shaded region, dSC and CFL flux tubes start to become relevant, its boundary is given by the condition ${H}_{\mathrm{c}1}={H}_{\mathrm{c}2}^{\left(1\right)}$ (see also the discussion around equation (69)).

The main conclusion of this phase diagram is that the maximal strange quark mass for which remnants of the domain wall (i.e. 2SC flux tubes with winding number larger than one) exist, is m_s ≃ 0.21μ_q. This is on the lower end of the range expected for compact stars. Therefore, at least within our approximations, we conclude that it is conceivable, but unlikely that 2SC flux tubes with exotic structures as shown in figure 2 play an important role in the astrophysical setting. However, our phase diagram shows that two-component flux tubes (with n = 1) do survive for larger values of m_s. Since they have a us condensate in their core, they can also be considered as remnants of the domain wall, albeit with conventional magnetic field structure. We have checked that if the phase diagram in the right panel of figure 6 is continued to about m_s = 0.5μ_q, there is still a sizable range (of roughly the same size as for m_s = 0.25μ_q), starting at T_c/μ_q ≃ 0.083, where 2SC flux tubes with us core are preferred. Consequently, they may well have to be taken into account for the physics of compact stars.

In the discussion of the applicability of our results to compact stars, we should also return to the actual magnitude of the magnetic fields considered here. Using the conversion to physical units given below equation (45) and the expressions for μ and λ in equation (15), we find that the magnetic field in Gauss can be computed from

$$\begin{equation}H=1.6{\times}1{0}^{19}\frac{H}{{\mu }^{2}/\sqrt{\lambda }}\enspace \frac{{T}_{\mathrm{c}}}{{\mu }_{\mathrm{q}}}\left(1-\frac{T}{{T}_{\mathrm{c}}}\right){\left(\frac{{\mu }_{\mathrm{q}}}{400\enspace \text{MeV}}\right)}^{2}G.\end{equation} \tag{ 71 }$$

Consequently, setting T = 0, the triple point at T_c/μ_q ≃ 0.05 in the left panel of figure 6, where NOR, 2SC, and 2SC flux tube phases meet, corresponds to about H ≃ 1.1 × 10¹⁹ G, while the critical fields at T_c/μ_q = 0.08 are H_c1 ≃ 1.3 × 10¹⁹ G and H_c2 ≃ 2.9 × 10¹⁹ G (note that H in physical units also varies along the horizontal direction of this plot since μ²/λ^1/2 depends on T_c/μ_q). Magnetic fields of this magnitude are getting close to the limit for the stability of the star, and they are several orders of magnitude larger than observed at the surface. It is thus highly speculative, although possible, that these huge magnetic fields are reached in the interior of the star. However, exotic flux tubes as discussed here may be formed dynamically, for instance if the star cools through the superconducting phase transition at a roughly constant magnetic field [16]. And, if color-magnetic flux tubes exist in compact stars, they help to sustain possible ellipticities of the star, resulting in detectable gravitational waves [60]. Therefore, even though the equilibrium situation studied here may never be reached in an astrophysical setting, it is important to understand the various possible—non-standard—flux tube configurations.