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

Dense quark matter is expected to behave as a type-II superconductor at strong coupling. It was previously shown that if the strange quark mass $m_s$ is neglected, magnetic domain walls in the so-called 2SC phase are the energetically preferred magnetic defects in a certain parameter region. Computing the flux tube profiles and associated free energies within a Ginzburg-Landau approach, we find a cascade of multi-winding flux tubes as"remnants"of the domain wall when $m_s$ is increased. These flux tubes exhibit an unconventional ring-like structure of the magnetic field. We show that flux tubes with winding numbers larger than one survive for values of $m_s$ up to about 20% of the quark chemical potential. This makes them unlikely to play a significant role in compact stars, but they may appear in the QCD phase diagram in the presence of an external magnetic field.


I. INTRODUCTION
Cold and dense matter is a color superconductor, in which certain color-magnetic fields are screened just like ordinary magnetic fields are screened in an electronic superconductor [1]. What happens if a color superconductor is placed in an external ordinary magnetic field? This question is, firstly, of theoretical interest because it addresses the phase structure of Quantum Chromodynamics (QCD). Especially for vanishing baryon density, where lattice calculations can be employed, the behavior of QCD in an external magnetic field has caught a lot of attention [2,3], motivated by the creation of large magnetic fields in heavy-ion collisions [4,5]. Secondly, the question is of phenomenological interest because color-superconducting quark matter may be found in the interior of neutron stars or in quark stars. Some of these compact stars (then called magnetars) are known to have huge magnetic fields on their surface. Whether the magnetic field in the core of the star is large enough to affect the physics on the QCD scale is unknown, but conceivable [6][7][8]. Possibly, although perhaps less likely, color superconductors may also be created in future collider experiments that aim to reach large baryon densities [9]. How a color superconductor responds to an external ordinary magnetic field strongly depends on the particular color-superconducting phase.
In this paper we mostly focus on the so-called 2SC phase [10], where only up and down quarks participate in Cooper pairing, while the strange quarks and all quarks of one color remain unpaired. In this phase, just like in the color-flavor-locked (CFL) phase [11], there is a certain combination of the photon and the gluons whose magnetic field penetrates the color superconductor unperturbed [12]. In other words, the Cooper pairs are neutral with respect to a "rotated charge", which is a combination of electromagnetic and color charges. As a consequence, a certain fraction of an external ordinary magnetic field can penetrate in the form of this rotated field. It turns out that at strong coupling (strong coupling constant much greater than the electromagnetic one) this fraction is in fact very large, such that only a small fraction of the magnetic field is expelled. Nevertheless, it may be energetically favorable to admit additional magnetic flux through certain magnetic defects. As in electronic superconductors, this scenario is referred to as type-II superconductivity. We expect the 2SC and CFL phases to be in the type-II regime for sufficiently large pairing gaps, which are estimated to be reached for realistic coupling strengths under compact star conditions [13][14][15][16].
The main idea of this paper is as follows. If the strange quark mass m s is neglected, then the usual 2SC phase is indistinguishable from the phase where only up and strange quarks form Cooper pairs. To avoid confusion, we refer to these phases as 2SC ud and 2SC us , respectively (the third possibility, 2SC ds , is different even for m s 0 because of the electric charges of the quarks). This enhanced flavor symmetry renders a domain wall configuration possible, which smoothly interpolates from 2SC ud to 2SC us . Indeed, there is a parameter regime where, due to the admission of additional magnetic flux through the wall, the formation of such domain walls becomes energetically favorable [16] 1 . At sufficiently large densities, m s can safely be neglected compared to the chemical potential, but this is no longer true in the context of compact stars (where we may still neglect the masses of up and down quarks). In this case, the free energies of 2SC ud and 2SC us phases become non-degenerate because there is an energy cost involved in including the massive strange quarks in the Cooper pair condensate. Therefore, the two phases can no longer coexist on two sides of a planar defect. Since we can view the domain wall as a flux tube with infinite radius, it can be expected that as we increase m s the planar defect turns into a flux tube with decreasing radius. As we shall see, this is realized in a descending sequence of multi-winding flux tubes, until we are left with flux tubes with winding number one at sufficiently large values of m s . The flux tubes in this sequence are 2SC ud flux tubes with a 2SC us core. Flux tubes with a core different from normal-conducting matter are possible in multi-component systems such as quark matter, for instance CFL flux tubes with a 2SC core, which are preferred over CFL flux tubes with a normal core for realistic parameters [16]. Related examples are two-component mixtures of superconductors with superfluids, where the core of the flux tubes is normal-conducting, but the superfluid condensate remains nonzero [18][19][20]. Our flux tubes are different from these examples in that the induced condensate in the core is zero far away from the center of the flux tube. They also differ from these cases with regard to the multi-winding solutions: While it has been pointed out that multi-winding flux tubes can become favored in CFL matter [21] and in superconductor/superfluid mixtures [18], this concerns a small region close to the transition point between type-I and type-II behavior of the superconductor, whose precise phase structure is further complicated by a first-order transition due to an attractive long-range interaction between the flux tubes [14,16,20,22]. The multi-winding 2SC flux tubes discussed here exist in a much larger parameter region within the type-II regime.
We shall work within a Ginzburg-Landau framework, which has often been used to study certain aspects of color superconductivity [13,[23][24][25][26]. In particular, our starting point is a direct extension of Ref. [16], making use of mass insertions considered previously [27][28][29]. While the Ginzburg-Landau approach in principle allows us to make model-independent predictions, it also has several shortcomings. Perhaps most importantly, it is purely bosonic and thus does not account for the fermionic constituents of the Cooper pairs. This is relevant since the magnetic fields considered here are large enough to potentially modify the microscopic structure of the Cooper pair condensates, as demonstrated in fermionic frameworks [30][31][32][33][34]. Moreover, for our numerical evaluation we shall employ the weakcoupling form of the Ginzburg-Landau parameters and extrapolate the results to large values of the strong coupling constant. Also, we do not attempt to construct any flux tube arrays, but rather focus on single flux tubes. This enables us to compute the critical magnetic field at which it becomes favorable to place a single magnetic flux tube into the system, H c1 in the usual terminology, but we cannot determine the precise structure of the phase above this critical magnetic field. Finally, we only include the effect of the quark mass to lowest nontrivial order and thus for large quark masses our results have to be considered with some care.
We emphasize that the flux tubes considered here are "pure" magnetic flux tubes, i.e., they have (quantized) magnetic flux, but zero baryon circulation. In color-superconducting matter, such flux tubes are not protected by topology [15,35]. We shall compute the parameter regime where they become energetically stable, i.e., where they cannot decay despite their non-topological nature. Since the 2SC phase is not a superfluid, pure magnetic flux tubes are the only possible line defects. For comparison, the CFL phase allows for defects with nonzero magnetic flux and nonzero baryon circulation [36,37]. As for ordinary superfluid vortices, their energy is formally divergent, while the magnetic flux tubes considered here have finite energy, just like magnetic flux tubes in an ordinary electronic superconductor.
The main difference of our magnetic flux tubes compared to the textbook scenario is the appearance of multiple condensates and gauge fields. Therefore, our study can also be put into the wider context of magnetic defects in unconventional superconductors. By considering the diagonal subsector of the order parameter in color-flavor space, we allow for three nonzero condensates, and as a consequence we can restrict ourselves to three gauge fields of the color and electromagnetic gauge group SU (3) × U (1). This is somewhat similar to electroweak strings with gauge group SU (2) × U (1) [38][39][40]. In the so-called semilocal approximation, where the SU (2) remains ungauged, these strings are described within an abelian Higgs model [41]. In this context, the flux tube profiles are often calculated at the transition point between type-I and type-II superconductivity, also referred to as the Bogomolny limit [42]. Our calculation is not restricted to this point and in fact we shall see that the multi-winding flux tubes are only stable away from this transition point. As in Ref. [41] we will see that for the multi-winding flux tubes the profiles of the magnetic field are ring-like, i.e., the maximum of the magnetic field is not at the center of the flux tube. This ring-like structure has also been observed in a model with a non-standard magnetic permeability [43], and it is similar to the experimentally observed structure of two-component superfluid vortices [44].
Our paper is organized as follows. We start with setting up the Ginzburg-Landau formalism in Sec. II, including the discussion of the mass terms and the rotated electromagnetism. In Sec. III we discuss the homogeneous phases as a preparation for the subsequent sections. The critical magnetic field H c and the upper critical field H c2 are computed in Sec. IV, and we set up the calculation of the flux tube profiles in Sec. IV C needed for the numerical calculation of the lower critical field H c1 . The main results are presented and discussed in Sec. V, which is divided into the discussion of the flux tubes themselves, Sec. V A, and the resulting phase diagram, Sec. V B. We give a summary and an outlook in Sec. VI. Our convention for the metric tensor is g µν = diag(1, −1, −1, −1). We work in natural units = c = k B = 1 and use Heaviside-Lorentz units for the gauge fields, in which the elementary charge is e = √ 4πα 0.3, where α is the fine-structure constant.

A. Ginzburg-Landau potential
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, [3] c and [3] 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 where the anti-symmetric 3 × 3 matrices (J i ) jk = −i ijk and (I j ) k = −i jk form bases of the three-dimensional spaces [3] c and [3] 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 2 . In this convention, for instance, Φ 11 carries anti-indicesr andū 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 Φ, Apart from the mass correction, proportional to the parameter , this is exactly the same potential, and the same notation, as in Ref. [16], where the starting point was the potential for Ψ, based on previous works [13,[23][24][25]. Due to the broken Lorentz invariance in the medium, the temporal and spatial parts of the kinetic term have different prefactors, u = 1/ √ 3. The covariant derivative is given by where g and e are the strong and electromagnetic coupling constants, respectively. The color gauge fields are denoted by A 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 a µν = ∂ µ A a ν − ∂ ν A a µ + gf abc A b µ A 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 Eq. (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, µ e ∝ m 2 s /µ q . The correction term is identical to the one used in Refs. [27,28], which is easily seen by an appropriate rescaling of Φ.
In the following we restrict ourselves to diagonal order parameters, Φ = 1 2 diag(φ 1 , φ 2 , φ 3 ) with φ i ∈ C, where φ 1 corresponds to ds pairing, φ 2 to us pairing, and φ 3 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 [45]. Even in this simplified case, 2 3 = 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, such that the only relevant gauge fields are the two color gauge fields A 3 µ , A 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 Within this ansatz, performing the traces in Eq. (2) yields 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 Ref. [46]. Here we have three condensates with identical (bosonic) chemical potential µ, with self-coupling λ, cross-coupling h, and three different effective masses (squared) All three gauge fields couple to the condensates. We can simplify the potential by a suitable rotation of the gauge fields.

B. Rotated electromagnetism and Gibbs free energy
We apply the double rotation from Ref. [16], which, denoting the rotated gauge fields byÃ 3 µ ,Ã 8 µ , andÃ µ , reads, where the mixing angles ϑ 1 and ϑ 2 are given by 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" ϑ 1 ensures that in the homogeneous 2SC phase, where only the condensate φ 3 is nonzero, only theB 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 φ 1 and φ 2 to be induced in the core of the flux tube. Therefore, we apply a second rotation with mixing angle ϑ 2 . This rotation leavesB 8 invariant and creates a field, namelyB, which is unaffected by the superconductor even if all three condensates are nonzero. Thus,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Ã 3 µ andÃ 8 µ are "almost" gluons with a small admixture of the photon, while A µ 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 Sec. III. Applying the gauge field rotation and writing the complex fields in terms of their moduli and phases, the Ginzburg-Landau potential (4) becomes with where we have introduced the rotated charges In Eq. (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) = He z . Therefore, we need to consider the Gibbs free energy density where V is the volume of our system. We can obviously assume that all induced magnetic fields have only zcomponents as well. Denoting the z-components of the rotated fields byB 3 ,B 8 , andB, we have H·B = H[sin ϑ 1B8 + cos ϑ 1 (sin ϑ 2B3 + cos ϑ 2B )]. SinceB does not couple to any of the condensates, it remains homogeneous even in the presence of flux tubes. Consequently, the equation of motion forÃ is trivially fulfilled, and we determineB by minimizing the Gibbs free energy, which, using Eq. (9), yields Reinserting this result into G, we obtain C. Parameter choices The potential U 0 (10) depends on the parameters µ, λ, h, . The discussion of the homogeneous phases in Sec. III 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 Sec. V, we employ the weak-coupling values of these parameters [13,14,25,27,28,35], 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 [47], 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 /µ 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 Ref. [48]. With the quark chemical potential in the core of a compact star of about µ q ∼ (400 − 500) MeV 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 Refs. [27,28] 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 Refs. [49,50]) With this result, it is convenient to trade the dimensionful parameter for the dimensionless "mass parameter" such that the complete dependence of our potential on the strange quark mass is absorbed in α, 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/(µ 2 /λ 1/2 ), and the ratio At weak coupling η = −1/2, as one can see from Eq. (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/(µ 2 /λ 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 twoloop QCD beta function (which should not be taken too seriously at such low densities), µ q 400 MeV corresponds to α s 1 and thus g = √ 4πα s 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 Eq. (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 Eq. (17) we will obtain an upper limit in m s /µ q for the presence of multi-winding 2SC flux tubes. Any T > 0 in Eq. (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 [20]. And, secondly, away from the asymptotic weak-coupling regime the phase transition becomes first order due to gauge field fluctuations [51], and thus at strong coupling the behavior just below the phase transition would have to be modified in a more sophisticated approach.

III. HOMOGENEOUS PHASES
Here we construct all possible homogeneous candidate phases within our ansatz and compute their Gibbs free energy density. Since we have aligned the z-axis with the magnetic fields, we may write the gauge fields asÃ 3 = xB 3 e y , A 8 = xB 8 e y with constantB 3 andB 8 . Furthermore, in this section we assume the condensates ρ i and their phases ψ i to be independent of r and thus all gradient terms in Eq. (10) vanish. Then, the equations of motion for the gauge fields and the condensates become and Since here we are only interested in homogeneous phases, the terms proportional to x 2 have to vanish separately in each equation (unless the equation is satisfied trivially by a vanishing condensate). The simplest case is the normal phase, where all three condensates vanish, ρ 1 = ρ 2 = ρ 3 = 0. Here we expect to have no induced color-magnetic fields and the external magnetic field should penetrate unperturbed, i.e., in the unrotated basis B 3 = B 8 = 0 and B = H. As a check, we can derive this result within our rotated basis. First we observe that the equations of motion (20) and (21) are trivially fulfilled for vanishing condensates. Then, from Eq. (10) we obtain U 0 = 0. Inserting this into the Gibbs free energy density (14) and minimizing the result with respect tõ B 3 andB 8 yields nonzero results for these rotated magnetic fields. By undoing the rotation one can check that these results together with the result forB (13) indeed give B 3 = B 8 = 0 and B = H, as expected. Then, substituting the magnetic fields at the minimum back into the Gibbs free energy density yields This result does not receive any mass corrections since, within our Ginzburg-Landau approach, unpaired fermions and any possible mass effects on them do not appear explicitly; the normal phase is the "vacuum" of our theory. We will now discuss the various superconducting candidate phases, which all do receive mass corrections. These mass corrections can be written in a general form using the masses m 1 , m 2 , m 3 . However, especially for the free energies this leads to some lengthy expressions, which are not particularly instructive. Thus we make use of Eq. (18) and express all free energies in this section to linear order in the mass parameter α.

A. CFL phase
In the CFL phase, all three condensates are nonzero. Here we obtainB 3 =B 8 = 0 and Eqs. (21) yield three coupled equations for the condensates, which have the solution where, in the second steps, Eqs. (18) and (19) have been used 3 . In the massless limit, α = 0, we recover the result of three identical condensates. The potential (10) becomes 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 where we have used the explicit form of the mixing angles (7).

B. 2SC phases
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 ρ 1 , ρ 2 , or ρ 3 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 ρ 3 . From Eq. (21c) we then obtain andB 8 = 0. By inserting this into the Gibbs free energy density G and minimizing with respect toB 3 we find and inserting this back into G yields Analogously, if only the condensate ρ 1 is nonzero, i.e., for the 2SC ds phase, we find as well asB 3 =B 8 = 0 from minimizing G, and we compute Finally, for the 2SC us phase, we find and the magnetic fieldsB This yields the Gibbs free energy density Comparing Eqs. (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.

C. fSC phases
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 ρ 2 = 0. In this case, there are gd/bs and gs/bd Cooper pairs from ρ 1 and ru/gd and rd/gu Cooper pairs from ρ 3 . 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 Ref. [28], we refer to this phase as dSC. Analogously, the phases with vanishing ρ 1 and ρ 3 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 ρ 2 = 0. With the help of Eqs. (20) we obtainB 3 =B 8 = 0. Then, Eqs. (21a) and (21c) yield two coupled equations for ρ 1 and ρ 3 , with the solution Thus we compute where again we have dropped the quadratic contribution in α. Analogously, if ρ 3 = 0, we find againB 3 =B 8 = 0, as well as and Finally, for ρ 1 = 0 we findB 3 =B 8 = 0, and Since the Gibbs free energy densities (35), (37), (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 [16], 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 Refs. [27,28].

A. Critical field Hc
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 η, Then, there are 4 transitions for which we can compute a critical magnetic field, The remaining transition, between the NOR and dSC phases formally has a critical field as well, 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 Ref. [16], where a different rotation was used. We see that if the magnetic field is given in units of µ 2 / √ λ, the parameter space has reduced to four dimensions, spanned by H, η, α, and g (for given e 0.30). In Fig. 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 This point also occurs in the right panel, where it has the coordinates The other triple point, where dSC, 2SC, and CFL phases meet, visible only in the left panel, is given by  Figure 1: Homogeneous phases at nonzero external magnetic field H in the plane spanned by H and the parameter ratio η = h/λ (describing the cross-coupling between the condensates) for fixed strong coupling constant g = 3.5 (left) and in the H-g plane for fixed η = −1/2 (right). In both panels, the solid curves are for a mass parameter α = 0.3, while the dashed curves correspond to the massless limit, α = 0.
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 c 14 µ 2 / √ λ, approximately corresponds to 4 H c 1.1 × 10 19 G. The main effects of the strange quark mass, manifest in Fig. 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 [28]. 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 [16], 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.

B. Upper critical field Hc2
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 4 Here we have used 1 G = 1 g 1/2 cm −1/2 s −1 and thus 1 G(c ) 3/2 = β eV 2 with the numerical factor β 0.06925. Together with √ 4π B HL = B G , where B HL and B G are the magnetic fields in the Heaviside-Lorentz and the Gaussian system of electromagnetism, we conclude that 1 eV 2 in natural Heaviside-Lorentz units corresponds to √ 4π/β G 51.189 G in the Gaussian system. 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 Ref. [16]. 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 Eq. (40) in Ref. [16], 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 κ = 1/ √ 2). Here, by equating H c2 with H c for the 2SC/NOR transition in Eq. (41) we find for this transition point where the weak-coupling expression of λ in Eq. (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 φ 3 is nonzero, to an inhomogeneous phase is slightly more complicated. Assuming a second-order transition, we linearize the Ginzburg-Landau equations in φ 1 and φ 2 . In the massless limit, the resulting two (decoupled) equations yield the same critical magnetic field [16]. In other words, as we approach the flux tube phase by decreasing H, both φ 1 and φ 2 become nonzero simultaneously (and continuously). This is different for nonzero m s , in which case the two relevant equations are where we have setÃ 8 = 0 sinceB 8 = 0 in the 2SC phase and where φ 3 = ρ 3 / √ 2 is the condensate in the homogeneous 2SC phase (26). With the usual arguments and using the 2SC relation betweenB 3 and H from Eq. (27), we obtain two different critical fields, The most relevant case for us is the one where both H (1) c2 and H (2) c2 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 (1) c2 from the homogeneous 2SC phase to a phase where both φ 1 and φ 3 are nonzero, which is an inhomogeneous version of the dSC phase (see Sec. III C). Then, as we reach the "would-be" H (2) c2 by further decreasing H, the approximation by which this critical field was computed is no longer valid, and thus the value for H (2) c2 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 (1) c2 with H c for the 2SC/CFL transition from Eq. (41). Dropping terms quadratic in α we find Since α depends on T c /µ q this is an implicit equation for T c /µ q . To lowest nontrivial order in m 2 s /µ 2 q the solution is 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 [16,20]. 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.

C. Flux tubes and lower critical field Hc1
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, 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 3 (∞) = 1 and f 1 (∞) = f 2 (∞) = 0. As in ordinary single-component flux tubes, we allow for a nonzero winding number n ∈ Z, such that the phases of the condensates are Here we have set the winding numbers for the "non-2SC" condensates f 1 and f 2 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 1 (∞) = f 2 (∞) = 0.) In such configurations, f 1 and/or f 2 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 [16] 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 Eq. (53). As a consequence, the boundary condition for the 2SC condensate in the core is f 3 (0) = 0, while f 1 (0) and f 2 (0) can be nonzero and must be determined dynamically.
As we have seen in Sec. II, after the rotation of the gauge fieldsÃ decouples from the condensates. Therefore, we are left with two nontrivial gauge fields, for which we introduce the dimensionless versionsã 3 andã 8 viã with the boundary conditionsã 3 (0) =ã 8 (0) = 0. This yields the magnetic fields where we have introduced the dimensionless variable and where prime denotes derivative with respect to R. Following Ref. [16], we have separated a term inÃ 3 for convenience, which gives rise to the nonzero fieldB 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,ã 3 (∞) =ã 8 (∞) = 0 (alternatively, this effect could have been implemented in the boundary condition forã 3 ). The behavior ofB 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: B, which fully penetrates the superconductor and thus is irrelevant for the calculation of the flux tube profiles;B 8 , which behaves analogously to the ordinary magnetic field in an ordinary flux tube; andB 3 , which is nonzero far away from the flux tube and is affected nontrivially by the flux tube profile. As a consequence, sinceB 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 Eqs. (13) and (55), we find where we have introduced the dimensionless magnetic field We can now express the potential U 0 (10) in terms of the dimensionless condensates and gauge fields, Here we have denoted the potential of the homogeneous 2SC phase by and we have abbreviated Inserting Eq. (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, and for the condensates, By taking the limit R → ∞ of Eq. (62b) we concludẽ while no condition forã 3 (∞) 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. The Gibbs free energy density can be written as where L is the size of the system in the z-direction, and with a closed integration contour encircling the flux tube at infinity, is the magneticB 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 F = πρ 2 2SC I with Written in this form, the free energy does not have any explicit dependence on the mass correction and is identical to the one in Ref. [16] (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 Eq. (65) is zero. This condition can be written as where Ξ c1 is the dimensionless version of H c1 via Eq. (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, Eq. (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 (1) c2 , which indicates the second-order transition from the homogeneous 2SC phase to an inhomogeneous phase where the condensate f 1 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 (1) c2 from Eq. (49a) (setting η = −1/2) in terms of the dimensionless field Ξ (58), which yields to linear order in α Solving this equation simultaneously with Eq. (68), i.e., setting Ξ (1) c2 = Ξ c1 , gives the intersection point. In the practical calculation, this is best done by inserting Eq. (69) into Eq. (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 Fig. 6.

V. NUMERICAL RESULTS AND DISCUSSION
We are now prepared to compute the flux tube profiles and the resulting critical fields H c1 , which we will put together with H c and H c2 from the previous section. As in the calculations of the critical fields H c and H c2 we use the bosonic masses (5), keep terms linear in α in Eqs. (63), and, as discussed in Sec. II C, we set η = −1/2 and g = 3.5 for all following results. Then we solve the coupled second-order differential equations (62) and (63) numerically without further approximations using the successive over-relaxation method, which has been used before in similar contexts [16,20,21,[52][53][54]. Each numerical run yields a flux tube profile for given m s /µ q , T c /µ q (from which α and λ are obtained), dimensionless external magnetic field Ξ, and winding number n. Since we are interested in the critical field Ξ c1 , we need to solve Eq. (68), for which we employ the bisection method. This requires solving the differential equations about 10 -20 times until a reasonable accuracy for Ξ c1 is reached. This whole procedure thus yields a critical field H c1 for given m s /µ q , T c /µ q , and n. As argued in the introduction and as we shall see below, solutions with high winding numbers are expected to play an important role. Therefore, in principle, the procedure has to be repeated for all n to find the preferred flux tube configuration for each point in the parameter space spanned by m s /µ q and T c /µ q . In practice, we have performed the calculation for the lowest few n, and for selected parameter sets for much larger n to check our conclusions. An additional complication arises because in certain parameter regions there is more than one solution to the set of differential equations. The single-component flux tube with f 1 ≡ f 2 ≡ 0 always exists. Configurations where the condensate f 2 is induced in the core of the flux tube, but not f 1 , turn out to be preferred over the single-component configuration whenever they exist and we shall discuss them in detail. Configurations where both f 1 and f 2 are induced in the core do exist as well, but only in a parameter region where H c2 indicates that the ground state is a more complicated flux tube array. We shall therefore not discuss these three-component configurations.

A. Flux tube properties
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 Ref. [16] (see Fig. 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 3 = 1, f 2 = 0, and on the other side f 2 = 1, f 3 = 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 Sec. III, and thus the domain wall configuration becomes unstable. In Fig. 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 critical fields as a function of the winding number for the parameters of Fig.  2 are shown in Fig. 3, which proves the successive decrease in the preferred winding from low to high strange quark mass.
Let us first discuss the profiles themselves in Fig. 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. Taking the point where the two condensates have the same value as a measure for this radius, we see that for the largest winding shown here, the radius is about R 6. With the help of Eq. (56) we can translate this number into physical units. Using Eqs. (26) and (15), we see that the translation depends on α, i.e., on the particular value of m s , and that, besides the ratio T c /µ q , it also depends on the quark chemical potential µ q itself (rewriting T c in λ (15) as µ q × T c /µ q ). One finds that the m s dependence for the values chosen here is very weak, such that, choosing µ q 400 MeV one obtains r 7.7 fm for R = 10, which corresponds to the radial domain shown in the figure. Consequently, the thickness of the high-winding flux tubes, say for the smallest mass chosen here, m s = 0.046µ q , is about r 4.6 fm. In other words, it only takes a relatively small strange quark mass to bring down the thickness of the flux tubes from infinity to a few fm.
The plots in the upper panels also show the profile of the function f 2 2 + f 2 3 , which is the radial coordinate of the space spanned by (f 2 , f 3 ), 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 Eq. R Figure 2: Upper panels: Flux tube profiles of the dimensionless condensates f2 and f3 (black) and f 2 2 + f 2 3 (red) as functions of the dimensionless radial coordinate R. The condensate f3 is the usual 2SC condensate of ud Cooper pairs and asymptotes to 1 for large R, while the us condensate f2 is induced in the core of the flux tube. For all plots g = 3.5 and Tc = 0.0856µq, while ms/µq assumes the values given in each panel, i.e., it increases from upper left to lower right. The masses are chosen such that the preferred configurations are flux tubes with winding numbers n = 10, . . . , 1 from upper left to lower right (and each plot shows the preferred configuration). For compact star conditions, R = 10 translates to about r 7.7 fm. Lower panels: Ratio of the induced magnetic field over the external magnetic field, B/H, in the X-Y plane perpendicular to the flux tube for the 10 configurations of the upper panels (X and Y in the same dimensionless units as R). The scale of the shading is adjusted for each plot separately, from black (maximal) to white (minimal). The magnetic field enters the superconductor in ring-like structures for small strange quark mass (large winding number) and turns into the conventional flux tube behavior for large mass (small winding number). This structure is reflected in the red curves of the upper panels.
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 Eq. (57) it is obvious that it is the small mixing angle ϑ 1 (sin ϑ 1 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 [41] and a one-component model with non-standard coupling between the condensate and the abelian gauge field [43] (see also Refs. [55,56] 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 Fig. 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 Fig. 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 Fig. 4 we show the value of the induced condensate in the core f 2 (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 2 (0) becomes more and more step-like for larger windings, i.e., f 2 (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.

B. Phase structure
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 Fig. 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 Eq. (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 Fig. 4. 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 Ref. [16], where it was computed explicitly by using a planar instead of a cylindrical geometry (see Figs. 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 Eq. (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 Fig. 5 with the critical fields H c from Sec. IV A and the critical fields H c2 from Sec. IV B. The result is shown in the left panel of Fig. 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 1 , following the notation of Ref. [16]. This configuration always has winding number 1. The region labeled by "CFL flux tubes" is bounded from above by H (1) c2 , which actually indicates a transition to dSC flux tubes (see discussion below Eqs. (49)). For the given parameters, the "would-be" critical magnetic field H (2) c2 is very close to H (1) c2 , 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.
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 Fig. 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 Fig. 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 1 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 c1 = H (see also the discussion around Eq. (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 Fig. 2 Figure 6: Left panel: Phase diagram in the plane of external field H and Tc/µq for ms = 0.15µq (black solid lines) compared to the massless result from Ref. [16] (red dashed lines). The numbers at the critical field for the emergence of 2SC flux tubes with us core indicate the winding number n, with n = ∞ indicating the domain wall in the massless case. The critical lines for the emergence of 2SC flux tubes without a us core (which always have winding n = 1) are marked with S1. Right panel: Phase diagram in the ms/µq-Tc/µq plane indicating the preferred flux tube configurations. In the shaded region dSC and CFL flux tubes appear, and the red arrows indicate the range where, at ms = 0, domain walls are preferred. This phase diagram shows one of the main results of the paper, namely that remnants of the domain wall in the form of multi-winding flux tubes survive up to ms 0.21µq. 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 Fig. 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 Eq. (45), we recall that the triple point in the left panel of Fig. 6, where NOR, 2SC, and 2SC flux tube phases meet, corresponds to about H = 1.1×10 19 G (note that the physical H also varies along the horizontal direction of this plot since µ 2 /λ 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 [15]. And, if color-magnetic flux tubes exist in compact stars, they help to sustain possible ellipticities of the star, resulting in detectable gravitational waves [57]. 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.

VI. SUMMARY AND OUTLOOK
Using a Ginzburg-Landau approach, we have investigated magnetic flux tubes in the 2SC phase of dense quark matter, i.e., flux tubes that asymptote far away from their center to a phase where u and d quarks form a Cooper pair condensate. Improving earlier studies, we have included the strange quark mass m s as a free parameter in the microscopic calculation of the flux tube profiles and their free energies. This is more realistic and makes a qualitative difference: In the presence of the strange quark mass, domain wall configurations that interpolate between ud and us condensates and that have previously been found to be energetically favored in a certain parameter region are no longer stable. We have shown that these domain walls -which can be viewed as flux tubes with infinite radiusturn into flux tubes with finite radius and high winding number as the strange quark mass is switched on. These unconventional flux tubes, which have a us condensate in their core, exhibit a ring-like structure of the magnetic field, i.e., the maximum of the magnetic field sits at a nonzero value of the radial distance from the center. Already for moderate values m s 0.05µ q , the flux tube radius is reduced to under 5 fm, and flux tubes with winding number larger than one survive up to about m s 0.2µ q . This roughly corresponds to the lower limit that is expected for the effective strange quark mass in a compact star environment. Therefore, at least within our approximations, we have concluded that the exotic multi-winding flux tubes -"remnants" of the domain walls from the massless limitprobably play no major role in compact stars, although flux tubes with a us core, but winding number one, survive up to much larger masses. Our study is also relevant for the QCD phase diagram at a nonzero external magnetic field and for a more general understanding of two-component magnetic flux tubes. Although the situation in dense quark matter is very specific due to the presence of color-magnetic fields and their mixing with the ordinary magnetic field, our 2SC flux tubes with ring-like magnetic fields are not unlike configurations previously found in different models [41,43].
Our results suggest a very rich structure of the 2SC flux tube phase for not too large values of m s , containing flux tubes of different winding numbers. It would be interesting to see if there are phase transitions between different flux tube arrays or if there are arrays composed of flux tubes with different winding. We have also pointed out that the upper critical field of the CFL flux tubes, i.e., the transition to a homogeneous 2SC phase, is affected qualitatively by the strange quark mass. Just below this upper critical field dSC flux tubes seem to be favored, such that also the CFL flux tubes phase appears to be more complicated than expected from the massless limit. Also the CFL flux tubes themselves may be revisited using our setup in order to systematically investigate the effect of the strange quark mass. In this case, there is no equivalent of the domain wall and thus the effect is probably less dramatic than in the 2SC phase.
We have used various approximations that can be improved in future studies. It would be interesting, but also very challenging, to start from a fermionic approach rather than from the bosonic Ginzburg-Landau model. This may modify our results since the magnetic fields considered here are sufficiently large to resolve the fermionic structure of the Cooper pairs. A more straightforward extension would be to explore the parameter space more systematically. We have restricted ourselves to the weak-coupling values for the parameters of the Ginzburg-Landau potential, and extrapolated the results to a large strong coupling constant realistic for compact stars. Also, we have neglected offdiagonal components in the order parameter matrix, which become potentially relevant due to the explicitly broken flavor symmetry. It would also be interesting to improve the treatment of temperature in the Ginzburg-Landau potential used here and in previous studies. Since there are several condensates which are not expected to melt at the same point, a more refined temperature dependence would include more than one critical temperature [20].