Stability of superfluid vortices in dense quark matter
Abstract
Superfluid vortices in the colorflavorlocked (CFL) phase of dense quark matter are known to be energetically disfavored relative to wellseparated triplets of “semisuperfluid” color flux tubes. However, the shortrange interaction (metastable versus unstable) has not been established. In this paper we perform numerical calculations using the effective theory of the condensate field, mapping the regions in the parameter space of coupling constants where the vortices are metastable versus unstable. For the case of zero gauge coupling we analytically identify a candidate for the unstable mode, and show that it agrees well with the results of the numerical calculations. We find that in the region of the parameter space that seems likely to correspond to realworld CFL quark matter the vortices are unstable, indicating that if such matter exists in neutron star cores it is very likely to contain semisuperfluid color flux tubes rather than superfluid vortices.
pacs:
12.38.t, 25.75.NqI Introduction
The densest phase of matter according to standard model physics is the colorflavorlocked (CFL) phase (Alford et al., 1998). The CFL condensate breaks the baryon number symmetry of the theory, hence the CFL phase is a superfluid, and CFL matter in the core of a spinning neutron star will carry angular momentum in the form of superfluid vortices. Unlike the fermions in terrestrial superfluids, quarks interact via a nonAbelian gauge group; the structure of the vacuum manifold (Balachandran et al., 2006) is and this permits the existence of a nonAbelian vortex configuration which is three times lower in energy than the usual superfluid vortex. This configuration consists of three widelyseparated semisuperfluid flux tubes, each carrying color magnetic flux. At separations much larger than the size of the core any two semisuperfluid flux tubes strongly repel each other (Nakano et al., 2009), and it has therefore been conjectured that CFL superfluid vortices will spontaneously decay into triplets of semisuperfluid flux tubes. However, the shortrange interaction between the flux tubes has not been calculated (Eto et al., 2014), leaving open the possibility that the vortices might be metastable, in which case the decay rate, occurring via barrier penetration, could be extremely low. Such metastability has already been established for vortices in an analogous system, a three component BoseEinstein condensate (Cipriani and Nitta, 2013).
In this paper we address these unresolved questions. To probe the stability of the CFL superfluid vortices we solve the classical field equations for the CFL condensate on a twodimensional lattice, analogously to previous calculations done for SU(2) YangMillsHiggs theory (Gleiser and Thorarinson, 2009),(Heinz et al., 1997). This approach gives us a full understanding of the decay process, far exceeding the insight that can be gained from asymptotic methods. The results presented in this paper are:

We map the regions in the parameter space of the couplings where the superfluid vortex is metastable (Section IV).

In the unstable regions we numerically extract the unstable mode (Section IV).

We analytically construct an unstable mode arising in the case of zero gauge coupling (Section III), and show that it is very similar to the numerically extracted unstable mode.

We clarify the nature of interaction of semisuperfluid vortices at short distances in the unstable region.
Ii Vortices and flux tubes
The nonAbelian GinzburgLandau Lagrangian used to describe semisuperfluid vortices Balachandran et al. (2006) in the CFL phase of dense quark matter is given by
(1)  
where , . represents the gluonic gauge field. We choose the normalization for the generators. In the CFL phase represents the colorflavorlocked diquark condensate. It is a 33 complex matrix. An element of may be denoted as , where is a color index and is a flavor index. In the symmetrybreaking phase the field in the ground state has a nonzero vacuum expectation value,
(2) 
where
(3) 
and the mass spectrum contains the Goldstone boson and two massive Higgs fields associated with perturbations along the singlet and adjoint directions. Their masses are and respectively. Stability of the ground state requires and . There are also massive gluons of mass .
The superfluid vortex is
(4) 
where is the radial profile that satisfies
(5) 
with boundary conditions
(6) 
However, as noted above, the superfluid vortex is not the lowest energy configuration in the topological sector of configurations with net winding number 1 at radial infinity. The lowest energy configuration consists of three (red, green, and blue) semisuperfluid color flux tubes, each with global winding . A red semisuperfluid flux tube solution is
(7)  
(8)  
(9) 
Green and blue flux tubes are obtained by permuting the diagonal elements. The profile functions , and obey
(11)  
(12) 
with boundary conditions
(13)  
(14) 
To understand why a configuration of three wellseparated semisuperfluid flux tubes has lower energy than a single superfluid vortex, compare the energy densities far from the core:
(15) 
The energy density arises entirely from the scalar field gradient, which for each component is proportional to , where is the net winding of the field. In the superfluid vortex there is net winding of 1 in each of the three diagonal components, whereas in the semisuperfluid flux tube there is net winding of in each component. So the energy density of a single semisuperfluid flux tube is of the energy density of a superfluid vortex. This leads to a repulsive force between the flux tubes, since the further apart they are, the more of space is filled with the energetically cheaper semisuperfluid field configuration, as opposed to the energetically costlier superfluid vortex field configuration. To estimate the leading term in the resultant potential
(16) 
we note that the energy density of the three semisuperfluid vortices can be approximated as being the same as a superfluid vortex at , and being a superposition of the three individual energy densities at . If we assume that the cores of the flux tubes have radius , and neglect the core contributions (which are independent of ), we find
(17)  
For large separation , there is a strong repulsive force decaying as . This justifies our neglect of contributions that would be subleading in , such as the core energies.
Iii Stability analysis of the superfluid vortex
The GinzburgLandau energy functional for a twodimensional static field configuration is
(18) 
with the Hamiltonian density
(19)  
where has been defined in equation (3). We define the spatial behavior of the superfluid vortex (4) by
(20) 
Up to an additive constant, the energy density of the vortex is
(21) 
To investigate the stability of the vortex, consider a perturbation which only affects the quark condensate, leaving the gauge field unperturbed,
(22) 
(We use the normalization for the generators.) To second order, the change in energy density due to the perturbation is given by,
(23) 
Focusing on perturbations in the color direction, and decomposing , we can write
(24) 
It can be shown that, for (i.e. ) the vortex is unstable to a perturbation of the form
(25) 
which corresponds to a translation of the red and green components of the vortex a small distance in the direction, and translation of the blue component a distance in the opposite direction.
Using (24) we obtain the energy of the perturbation to order ,
(26) 
This is the main result of this section. We see that if (i.e. ) then the perturbation (25) lowers the energy of the vortex. At this point this is just a guess: there might be a lowerenergy perturbation that involves the gauge field or has a different spatial profile or color structure. However, we will see in Sec. IV that the numerically obtained unstable mode matches (23) very closely.
Iv Numerical results
To analyze the instability of the superfluid vortex and map the unstable/metastable boundary in the parameter space spanned by the three couplings , , and , we solve the classical field equations on a twodimensional lattice. For details see Appendix A. For all numerical calculations we chose the mass scale to be and use a lattice spacing . This provides an adequate resolution for the superfluid vortex solution on the lattice for all parameter values that we studied, since the size of the vortex depends only on , not on any of the couplings.
We use the ket as a convenient notation for the actual lattice configurations of matrices in 2dimensional position space. We define the following inner product and norm
(27) 
To characterize the degree to which field configurations resemble each other, we introduce the notion of an angle between two lattice configurations and ,
(28) 
so when the configurations are the same up to an overall multiplicative factor.
iv.1 Numerical analysis of the unstable mode
In parameter regions where the superfluid vortex is unstable, the vortex solution has an unstable mode, and spontaneously decays to three semisuperfluid flux tubes that repel each other. Snapshots of this process are shown in Figure 1 where we plot the total energy density of the system along the vertical axis. The instability arises from a direction in the highdimensional configuration space of perturbations to the superfluid vortex along which the potential possesses a negative curvature. This direction is the unstable mode , and its time evolution is
(29) 
In parameter regions where the vortex is metastable, there is no such unstable direction, and all perturbations oscillate but remain small.
To find out whether, for a given set of parameter values , , and , the superfluid vortex is unstable or metastable, we proceed as follows.
We first generate a lattice field configuration that is the exact superfluid vortex solution of the lattice field equations. We do this by transferring the continuum vortex solution (with radial profile obeying Eq. (5)) to the lattice and then allowing it to relax to the lowest energy state by evolving it using the Langevin equation with damping but no noise. We can do this even when the superfluid vortex is unstable because the initial approximate vortex configuration is proportional to the unit matrix in the colorflavor space of complex matrices, and the colorflavor symmetry of the Lagrangian guarantees that under time evolution it will remain proportional to the unit matrix, whereas decay would require the generation of nonsinglet color components.
To probe the stability of the equilibrated vortex we
add a small perturbation to the superfluid vortex, and
evolve it forward in time.
The exact form of the perturbation is not important, as
long as it has some overlap with the unstable mode (if any).
If the vortex is unstable then even a tiny initial
perturbation with some component along the unstable direction
will grow exponentially as we evolve forward in
time. It quickly dominates the other components of the initial perturbation.
We tried three different perturbations for the field (initial gauge links put to unity without perturbation):

a random configuration
Here we used 8 complex random numbers (uniformly distributed between zero and ) as entries for the matrices at each lattice site. 
the analytically obtained unstable mode for
This mode is constructed using (25), which is controlled by two parameters, the direction and magnitude of the displacement, . Using the set as a complete basis, the radial profile serves as the position dependent coefficient of one or more of the basis elements , since is a stable direction. We usually used only, but one can equally well choose any other basis element or combinations thereof. 
the numerically obtained unstable mode (see below).
During an initial transient period (which lasts longest for the random initial perturbation) the unstable mode grows to become the dominant component. After the overlap of the growing mode with the original perturbation grows exponentially in time, following (29),
(30) 
The exponential growth ends after time when the amplitude is so large that nonlinearities become nonnegligible. Using (30) we can measure the growth rate without knowing the actual form of the unstable mode. If such an exponential growth is observed, the superfluid vortex is unstable. The growth rate is determined by the negative curvature of the energy in the unstable direction, and is independent of the initial “seed” perturbation.
For an unstable vortex we can numerically obtain the unstable mode up to an overall normalization factor, by computing
(31) 
where and similarly for . Once we know the unstable mode , we can go back and repeat the procedure described above using that mode as the initial seed perturbation, in which case the initial transient time is very short and exponential growth starts immediately.
During the epoch of exponential growth we use (28) to compute the angle between the growing mode and the initial perturbation. When using the numerically constructed unstable mode as a seed we find that is very small, typically in the range to 2 degrees. When using an unstable mode that has been constructed analytically by Eq. (25), is typically of the order of degrees. Thus, the analytic mode discussed in Section (III) seems to be identical to the unstable mode, at least as far as the case is concerned. Allowing for (larger) nonzero values of the gaugecoupling changes the picture, as indicated by Figure 4. The transition line between unstable and metastable configurations is not parallel to the axis, suggesting that the gauge field plays a more prominent role in the decay process as grows larger. The angle in the case of the random mode is slightly larger as compared to the case where we used the analytic mode to perturb the system. This discrepancy is most likely a numerical artifact, as we are operating with very small numbers to isolate and extract the unstable mode in this very highdimensional space.
iv.2 Parameter space scan
Using the procedure described above to determine whether the vortices for a gives set of couplings are unstable or metastable, we performed a scan of the parameter space for , , . The translation of these couplings to physical values is discussed in Section IV.3 below. Our main findings are summarized in Figures 2 and 3.
The critical surface separating the metastable from the unstable regime is depicted in Figure 2. We see that only a small subspace of the surveyed parameter space yields metastable superfluid vortices. The plot suggests that at higher values of the gauge coupling (), which is the relevant region for QCD, the superfluid vortices are unstable, spontaneously decaying in a short time. Larger values of the coupling also yield unstable vortices.
In Figure 3 we sliced the threedimensional parameter space along the physically most relevant coupling among the three, the gauge coupling . With increasing , the area of metastable solutions decreases rapidly. The dotted line in the plot corresponds to , which is the line along which the analytic analysis of Section III predicts the stable/unstable transition. For small values of , this seems to be approximately true. As the gauge coupling increases, however, the transition boundary deviates more and more from the predicted line, which is most likely due to the gauge field playing an increasingly important role in the decay process.
iv.3 Relating stable and unstable regions to ratio of masses
In Ref. Eto and Nitta (2009) it was suggested that the short range interaction between semisuperfluid flux tubes, and hence the unstable/metastable boundary for the superfluid vortex, is dictated by the hierarchy of the masses , and (Section II). In Figure 4 we investigate this proposal by choosing one slice of the parameter space along the plane and plotting the different mass hierarchies along with the unstable/metastable boundary. The vertical line corresponds to the line along which . We could not find a direct agreement between the unstable/metastable boundary and any of the boundaries derived from arranging the masses according to their mass hierarchy.
iv.4 Relation between the effective theory and QCD
At arbitrary density it is not possible to calculate the couplings in our effective theory in terms of the microscopic physics, namely QCD, because QCD is strongly coupled. However, in the ultrahigh density regime, where the coupling becomes weak, the parameters and have been calculated using the mean field approximation Iida and Baym (2001), Giannakis and Ren (2002) in terms of baryon chemical potential and the transition temperature of the CFL condensate,
(32) 
Using this expression, the range of values of and that we have explored in our study correspond to values of ranging from 400 MeV to 500 MeV and values of ranging from 10 MeV to 15 MeV. As we can see, in the weak coupling mean field calculation , i.e. . In our calculations the superfluid vortex is unstable for . We illustrate this in Figure 3 where the dashed straight line is where , and even at very small QCD coupling this line is in the unstable region. Increasing the value of the coupling constant shrinks the region of metastability away from the plane , and increasing the value of just takes us to larger and , so from Figure 3 it seems likely that the vortices will remain unstable at large and .
V Conclusion and discussion
We have studied the stability properties of the superfluid vortices in the CFL phase of dense quark matter. Using a GinzburgLandau effective theory of the condensate field discretized on a twodimensional spatial lattice, we evolved the vortex configuration in time and looked for an exponentially growing unstable mode. We scanned the parameter space of the couplings, mapped the regions where the vortices are unstable as opposed to metastable, and in unstable regions we identified the unstable mode. We found that the region where superfluid vortices are metastable is rather small. Vortices are metastable when the gauge coupling is sufficiently small. In that case, the transition line separating the metastable from the unstable direction is almost as predicted in equation (26), that is, , see upper left panel in Figure 3. At larger values of , the metastable region shrinks and vanishes around . We could not find any sign of metastability above values of . If we use meanfield weakcoupling calculations to relate the GinzburgLandau couplings to QCD parameters such as , , and then it seems likely that CFL vortices in neutron stars would be unstable rather than metastable, but these calculations are not valid in the density range of interest for neutron stars. If better calculations of the effective theory couplings become available, the physical region in our plots could be more precisely identified.
It is very interesting that a superfluid vortex can be rendered unstable by solely perturbing the quark condensate, in spite of the fact that in CFL matter gauge fields play an essential role in the later states of the decay process, in which the vortex separates into three color flux tubes that repel each other. In a theory without gauge fields there would be no such repulsive force, and our preliminary calculations suggest that in that theory the superfluid vortex decays into a moleculelike configuration of three separate vortices which remain bound at a fixed spacial separation.
This work opens up several avenues of future enquiry. We only investigated the early stages of the process of vortex decay. For regions of parameter space where the vortices are metastable it would be interesting to measure their lifetime and evaluate the height of the energy barrier. Our GinzburgLandau theory did not include entrainment (currentcurrent) interactions, and it would be interesting to study how that affect our results. The same methods that we used could be applied to vortices in the colorspinlocked phase of quark matter, and to study the stability of the proposed colormagnetic flux tubes in twoflavor color superconducting quark matter.
Acknowledgements.
This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award Number #DEFG0205ER41375. AW acknowledges support by the Schrödinger Fellowship J 3800N27 of the Austrian Science Fund (FWF). TV thanks Washington University for hospitality. TV acknowledges support by the U.S. Department of Energy, Office of Science, Office of High Energy Physics under Award number #DESC0013605.Appendix A Equations of motion on the lattice
We solved the equations of motion on a twodimensional spacial lattice whose lattice spacing we denote by . The matter field and its conjugate momentum are represented by complex matrices and with one color and one flavor index which live on the lattice sites. The gauge fields are unitary matrices with two color indices living on the links, and their time derivative is encoded in the electric field which is a Hermitian matrix. In what follows, and are discretized versions of spatial coordinates and , and represents time. represents a gauge link at time that emanates from the site along the direction. We perform our calculations in the temporal gauge, . One can construct a plaquette at time , denoted by , by starting at site and going in direction. After each step, the orientation is changed by a 90degree turn to the left, such that the plaquette closes after four steps. At any given time , starting at site , one can thus construct four different elementary plaquettes, depending on the direction of the initial step. The plaquette then corresponds to the product of the four link variables in the order of their appearance. For example, the plaquette becomes the matrix product .
The lattice energy functional is
(33) 
where
(34)  
(35)  
(37) 
The lattice equations of motion derived from the above energy functional are
(38) 
where
(39)  
(40) 
(41) 
(42) 
where ’s are the SU(3) generators. We choose the normalization .
(43) 
In order to obtain the initial superfluid vortex configuration on the lattice, we use the Langevin evolution method. The continuum profile of the superfluid vortex is first evolved using the Langevin approach, with temperature set to zero. We then take the final relaxed configuration and use it as the input configuration for our numerical studies on stability. In the Langevin implementation, equations (42) and (43) are modified as follows :
(44) 
(45) 
where is the coefficient of viscosity. By the fluctuationdissipation theorem the stochastic noise terms , , are independently drawn from the same Gaussian probability distribution
(46) 
where is a Gaussian random number of zero mean and unit variance. is the temperature of the thermal bath that is coupled to the system.
We tested two types of boundary conditions (BC): fixed and Neumann. Fixed BC consisted of locking the matter and gauge fields at the edge of the lattice to the values they would take when there is a superfluid vortex at the center of the lattice. With fixed BC the boundary affects the later stages of the decay of a superfluid vortex because the edges repel the semisuperfluid flux tubes. For Neumann BC we fixed the matter and gauge fields at the boundary to be the same as their neighbors one lattice spacing in. This sets the gradient of the field configuration to zero at the edge. With Neumann BC the later stages of the decay behave correctly, with the three semisuperfluid flux tubes leaving the lattice and disappearing across the boundary. However, with Neumann BC the early states of the decay were affected by a very slight attraction of the vortex to the boundary, so if the vortex lives for long enough it starts moving slowly towards the boundary. We therefore used fixed BC in our calculations.
In order to study finite size effects, we used different lattice sizes for our calculations. Calculations that led to the main results of this paper were performed on a lattice and were repeated on a lattice of size , where we found no significant discrepancy in the results. We furthermore varied the criteria in the code that are responsible for detecting the exponential growth of the unstable mode. We found our results to be robust under these changes as well.
References
 Alford et al. (1998) M. G. Alford, K. Rajagopal, and F. Wilczek, Phys. Lett. B422, 247 (1998), arXiv:hepph/9711395 [hepph] .
 Balachandran et al. (2006) A. P. Balachandran, S. Digal, and T. Matsuura, Phys. Rev. D73, 074009 (2006), arXiv:hepph/0509276 [hepph] .
 Nakano et al. (2009) E. Nakano, M. Nitta, and T. Matsuura, Phys. Lett. B672, 61 (2009), arXiv:0708.4092 [hepph] .
 Eto et al. (2014) M. Eto, Y. Hirono, M. Nitta, and S. Yasui, PTEP 2014, 012D01 (2014), arXiv:1308.1535 [hepph] .
 Cipriani and Nitta (2013) M. Cipriani and M. Nitta, Phys. Rev. A88, 013634 (2013), arXiv:1304.4375 [condmat.quantgas] .
 Gleiser and Thorarinson (2009) M. Gleiser and J. Thorarinson, Phys. Rev. D79, 025016 (2009), arXiv:0808.0514 [hepth] .
 Heinz et al. (1997) U. W. Heinz, C. R. Hu, S. Leupold, S. G. Matinyan, and B. Muller, Phys. Rev. D55, 2464 (1997), arXiv:hepth/9608181 [hepth] .
 Eto and Nitta (2009) M. Eto and M. Nitta, Phys. Rev. D80, 125007 (2009), arXiv:0907.1278 [hepph] .
 Iida and Baym (2001) K. Iida and G. Baym, Phys. Rev. D63, 074018 (2001), [Erratum: Phys. Rev.D66,059903(2002)], arXiv:hepph/0011229 [hepph] .
 Giannakis and Ren (2002) I. Giannakis and H.c. Ren, Phys. Rev. D65, 054017 (2002), arXiv:hepph/0108256 [hepph] .