{'fraction_non_alphanumeric': 0.03917936240218984, 'fraction_numerical': 0.013078220477230625, 'mean_word_length': 4.03241964658411, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Flux tubes and the type-I/type-II transition in a superconductor coupled to a superfluid

We analyze magnetic flux tubes at zero temperature in a superconductor that is coupled to a superfluid via both density and gradient (``entrainment'') interactions. The example we have in mind is high-density nuclear matter, which is a proton superconductor and a neutron superfluid, but our treatment is general and simple, modeling the interactions as a Ginzburg-Landau effective theory with four-fermion couplings, including only s-wave pairing. We numerically solve the field equations for flux tubes with an arbitrary number of flux quanta, and compare their energies. This allows us to map the type-I/type-II transition in the superconductor, which occurs at the conventional kappa = 1/sqrt(2) if the condensates are uncoupled. We find that a density coupling between the condensates raises the critical kappa and, for a sufficiently high neutron density, resolves the type-I/type-II transition line into an infinite number of bands corresponding to ``type-II(n)'' phases, in which n, the number of quanta in the favored flux tube, steps from 1 to infinity. For lower neutron density, the coupling creates spinodal regions around the type-I/type-II boundary, in which metastable flux configurations are possible. We find that a gradient coupling between the condensates lowers the critical kappa and creates spinodal regions. These exotic phenomena may not occur in nuclear matter, which is thought to be deep in the type-II region, but might be observed in condensed matter systems.

saturation density, where the critical temperatures for the superfluid and superconductor are of order MeV [8][9][10]. The temperature of a compact star drops below this value within minutes of its formation in a supernova, and is at or below the keV range after the first 1000 years [11]. When we discuss type-I versus type-II behavior we are referring to the response of the system to a magnetic field at the lower critical value, at T = 0.
As far as we know, there has been no previous work on how a flux tube in a superconductor is affected by a gradient coupling to a co-existing superfluid. However, there has been work on possible knot solitons [12], vortices in the SO(5) model of high-temperature superconductivity [13], and on the complementary situation, a superfluid vortex with gradient coupling to a co-existing superconductor. There the coupling leads to the "entrainment" or Andreev-Bashkin effect [14] whereby the proton condensate is dragged along with the neutron condensate, producing a non-zero proton current around the vortex, dressing it with some magnetic flux [15]. It is interesting to note that this flux is not a multiple of the flux quantum for proton flux tubes. This is possible because of the difference between the energetics of a neutron vortex and a proton flux tube. The flux tube has energy density localized to the vicinity of its core. Far from the core the energy density must vanish, which means the proton field must change in phase by a multiple of 2π, and the the vector potential must cancel the resultant gradient, leading to a quantized magnetic flux. A neutron vortex, by contrast, has gradient energy that is not localized to the vicinity of the vortex, and the total energy per unit length diverges in the infinite volume limit. The vector potential is therefore not constrained to cancel any gradient in the proton field, and takes on a value that minimizes the overall energy, with no quantization condition on the resulting magnetic flux.
Returning to the situation that we study, a proton flux tube in a neutron superfluid background, we do not expect a similar behavior. This is because the proton flux tube's energy density is localized around its core, giving it (unlike the neutron vortex) a finite energy per unit length. If the neutron condensate were entrained, and developed non-zero circulation around the flux tube, it would acquire a non-localized energy density, leading to an infinite energy per unit length for the flux tube, which is clearly energetically disfavored. We will see below that the effect of gradient couplings on the proton superconductor is more subtle: it leads to metastable regions near the type-I/type-II boundary.

II. STABILITY OF FLUX TUBES
Our aim is to explore the response of the proton superconductor to an applied critical magnetic field at zero temperature. We will therefore construct a phase diagram in the space of the coupling constants of the Ginzburg-Landau effective theory. We would like to be able to specify when it is of type II (at the lower critical magnetic field, flux tubes appear, and remain separate, i.e they repel) and when it is of type I (at the critical magnetic field, macroscopic normal regions appear, i.e. the flux tubes attract and coalesce). The simplest way to do this is to calculate the energy per unit length E n of a flux tube containing n flux quanta. The same approach has been used for vortices in the SO(5) model [16]. It is convenient to work in terms of the energy per flux quantum, When B n is negative the n-quantum flux tube is stable against fission into many single quantum flux tubes, and it is energetically favorable for n single quantum flux tubes to coalesce into one n-quantum flux tube. When B n is positive the n-quantum flux tube is unstable against fission, and coalescence is energetically disfavored. If one calculates B n for all n then the energetically favored value of n is the one that minimizes B n . In a traditional type I superconductor, small flux tubes attract each other and amalgamate into large ones and ultimately into macroscopic normal regions, so we would expect to find B n < 0 with its value dropping monotonically as n rises. In a type II superconductor we would expect B n > 0, with its value rising monotonically with n. Our calculations confirm these results for a single superconductor, but we will see that B n shows more complicated behavior when the superconductor feels interaction with a co-existing superfluid.
Calculations of B n are straightforward because they always occur in a cylindrically symmetric geometry, so the problem is one-dimensional. For a more detailed understanding of flux tube interactions, one would have to consider two single-quantum flux tubes a distance d apart. Their total energy is U (d), where U (0) = E 2 and U (∞) = 2E 1 , so B 2 = 1 2 (U (0) − U (∞)). As expected, B 2 < 0 means that the flux tubes have lower energy when they amalgamate, and B 2 > 0 means that the flux tubes have lower energy when they separate. If U (d) is monotonic, we can conclude that flux tubes either coalesce (B 2 < 0) or repel to infinite separation (B 2 > 0), corresponding to type-I or type-II behavior respectively. However, if there is a minimum in U (d) at some favored intermediate separation d = d * then irrespective of the sign of B n , one has a new variety of type II superconductor with some favored Abrikosov lattice spacing d * . Such behavior has been found to arise from a φ 6 term [17] and in the case of two charged condensates [7]. Calculating U (d) in the current context is an interesting but demanding problem which we leave for future work. In this paper we assume that U (d) is monotonic, so to analyze the attractiveness/repulsiveness of the flux tube interactions it is sufficient to calculate B n , or equivalently E n /n.

III. FLUX TUBES IN THE GINZBURG-LANDAU MODEL
A. Ginzburg-Landau model We start by writing down the zero-temperature Ginzburg-Landau effective theory of proton and neutron condensates in the presence of a magnetic field [5,18]. We denote the proton condensate field by φ p , the neutron condensate field by φ n , and the magnetic vector potential by A. The free energy density is where m c is twice the nucleon mass, q is twice the proton charge, U ent is the entrainment free energy density (see [18]) and σ is a parameter characterizing the strength of the gradient coupling, µ p and µ n are the chemical potentials of the proton and neutron condensate excitations, and a pp , a nn , and a pn are the GL quartic couplings. In zero magnetic field, the condensates would have position-independent bulk densities φ p 2 and φ n 2 obtained by minimizing the free energy. This allows us to eliminate the chemical potentials µ p , µ n by writing µ p = a pp φ p 2 + a pn φ n 2 µ n = a nn φ n 2 + a pn φ p 2 (5) so up to constants involving φ p and φ n , the potential V can be expressed in terms of the deviations of the condensate fields from their bulk values: In a neutron star, electrical neutrality keeps the proton fraction small, in the 5% to 10% range [19,20]; we will take φ p 2 / φ n 2 ≈ 0.05. As we now argue, a typical value for the entrainment coupling is σ ∼ 10 −1 . We first relate our formalism to the hydrodynamic limit of the free energy, following [18]. We focus on the phases of the fields, φ p = φ p exp(iχ p ) and φ n = φ n exp(iχ n ), and assume the fields have constant magnitude, and their phases have gradients The free energy density (2) then reduces to the hydrodynamic form where the symmetric matrix ρ of superfluid densities has elements Our entrainment parameter σ is therefore related to the parameter ǫ of Ref. [21][22][23] by σ = ǫ φ n / φ p . Since ǫ is of order 0.03, and φ n 2 / φ p 2 ∼ 20, we expect σ ∼ 10 −1 . This is consistent with the estimate ρ pn ≈ − 1 2 ρ pp used by [18]. In terms of the Andreev-Bashkin parametrization [14], ρ 12 = −ρ pn , ρ 1 = ρ pp + ρ pn , ρ 2 = ρ nn + ρ pn , so ρ 1 /ρ 12 ∼ 1, ρ 2 /ρ 12 ∼ 40. All the interactions in (8), including the entrainment, have their ultimate origin in the strong interaction between the nucleons, which is isospin symmetric, and hence does not distinguish protons from neutrons.

B. Flux tube solutions
To study a flux tube containing n flux quanta, we assume a cylindrically symmetric field configuration in which the proton condensate field winds (in a covariantly constant way) around the z-axis with a net phase 2πn, We have defined φ n as a real field, because, as noted above, any net phase change in the neutron condensate when it circles the flux tube would cost an infinite energy per unit length. Inserting the ansatz in (2) we obtain Generating the Euler-Lagrange equations using the standard procedure, we obtain a set of coupled differential equations for f , g and a: At this point we recall the definition of the Ginzburg-Landau parameter κ = λ/ξ, where the London penetration depth λ and superconducting coherence length ξ are (see [24]) To further simplify the equations, we then change variables to a dimensionless radial coordinater = r/ξ, obtaining The free energy per unit length of the flux tube, in terms of the variabler, is  (18). With no coupling between the condensates (β = σ = 0), the neutrons are undisturbed (δρn = 0). With a non-zero density coupling β, the neutron condensate (broken lines) is significantly perturbed by the flux tube. Note that the neutron δρn's are multiplied by 10 (not by 100 as in Fig. 2) to make them visible. The other parameters are κ = 3.0, σ = 0.0, and φp 2 / φn 2 =.05.
In addition to the system of equations, we require boundary conditions on the fields at the origin and at ∞. Far from the flux tube core, the fields will go to their uniform condensate value, so f (∞) = g(∞) = a(∞) = 1. Near the origin, f (r) ∝ r n , a(r) ∝ r 2 and g(r) is a constant. Therefore we have the conditions f (0) = 0, a(0) = 0 and g ′ (0) = 0. To obtain the energy of a flux tube we numerically solve the ODE system for the neutron and proton condensate and magnetic potential profile functions, then calculate the free energy of the system by inserting the results into (17) and integrating.
The system has five independent parameters: a pp , a nn /a pp , a pn /a pp , σ, and φ n / φ p . In neutral nuclear matter, the density of protons (neutrons) is proportional to φ p 2 ( φ n 2 ), and the proton density is approximately 5% of the total baryon number density [18], so we set φ p 2 / φ n 2 = .05 in most of our analysis. Following [5,6] we set a nn = a pp , and use (15) to exchange the parameter a pp for κ, which is the conventional parameter used in condensed matter studies of superconductivity. Our reduced set of parameters is therefore κ, the proton-neutron gradient coupling σ, and the proton-neutron density coupling β ≡ a pn /a pp . We also study some effects of varying φ p 2 / φ n 2 .

A. Flux tube solutions
For given values of φ p 2 / φ n 2 , κ, the proton-neutron gradient coupling σ, and the proton-neutron amplitude coupling β ≡ a pn /a pp we numerically solved the equations of motion (16) giving the field profiles for flux tubes with various numbers n of flux quanta. We obtained the solutions using a finite-element relaxation method, which is much less sensitive to initial conditions than the traditional "shooting" method, and better suited to repeatedly solving the equations for different sets of parameters. Next, we insert the solution for each profile into our expression for the free energy (17) and numerically integrate it to obtain a value for E n .
To estimate the numerical errors in our results, we varied the convergence criterion in the finite-element relaxation calculation, the spacing of the radial grid of points, and the radius out to which the grid extended. We found that the resultant variation in E n /n was of order 10 −6 , so numerical errors are invisible on the scale of the plots shown in Fig. 3.
Having obtained E n we then plot the series B n to determine whether the system is type I or type II for the chosen point in parameter space. In this way we find the points in parameter space where the system changes from a type I state to a type II state. Taking various slices through the parameter space, we can generate phase diagrams that  (18). With no coupling between the condensates (β = σ = 0), the neutrons are undisturbed (δρn = 0). With a non-zero gradient coupling σ, the neutron condensate (broken lines) is slightly perturbed by the flux tube. Note that the neutron δρn's are multiplied by 100 to make them visible. The other parameters are κ = 3.0, β = 0.0, and φp 2 / φn 2 =.05.
show the boundary curves between the various phases.
Figs. 1 and 2 each show a profile for a flux tube with a single flux quantum n = 1 on the left, and a profile for a flux tube with 100 flux quanta on the right. Fig. 1 shows the effect of non-zero density coupling β and Fig. 2 shows the effect of non-zero gradient coupling σ. We have plotted the normalized difference in density of the pair fields from their condensate values,

No coupling to neutrons
We do not show a plot of the flux tube profile for a simple superconductor, since this is well known: in a core region whose area rises as the number of flux quanta n, the proton condensate is suppressed; in a wall region the condensate returns to its vacuum value. At the Bogomolnyi point [25], κ = 1/ √ 2, the energy per flux quantum is independent of n [26], but on either side of this value there are area and perimeter contributions to the energy [27], so for κ close to 1/ √ 2 we expect the energy of a flux tube in a simple superconductor to have the following dependence on n, This is an expansion around n = ∞, but our numerical results will show that it works down to n = 1. We define δκ ≡ κ − 1/ √ 2. E Bog is the energy per unit flux at δκ = 0. By convention we take the parameter M , which has dimensions of energy, to be positive. The value of c 1 2 is then positive, ensuring that for δκ > 0, n = ∞ is disfavored (type-II), and for δκ < 0, n = ∞ is favored (type-I). We will see this behavior in our numerical results (Sec. IV B 1 and upper left plot of Fig. 3).

Density coupling to neutrons
For positive β, which corresponds to positive a np , equations (2) and (4) indicate that there is a repulsion between the neutron and proton condensates, so in the center of the flux tube, where the proton condensate is suppressed, the neutron condensate will be enhanced. That is exactly what we see in Fig. 1, where the dashed curve, showing the perturbation to the neutron density ρ n , rises inside the flux tube. For negative β there is attraction between the two condensates, and the neutron condensate is suppressed inside the flux tube (dash-dotted line). We therefore expect that the leading correction due to the interaction will be proportional to the core area, i.e proportional to n. The energy of an n-quantum flux tube is then where E (sc) n (κ) is the energy for an n-quantum flux tube in a pure superconductor, with no coupling to a superfluid (19). The leading correction is −M β n, which should be negative and quadratic in β for small β (see Sec. IV A 4), so the interaction energy parameter M β is positive and proportional to β 2 . The sub-leading term proportional to √ n arises from the energy cost of the gradient in ρ n at the edge of the flux tube, where it must return to its vacuum value, so we expect this term to be positive: b 1 2 > 0. We do not have an a priori expectation for the sign of the sub-sub-leading term b 1 .

Gradient coupling to neutrons
For positive σ, we expect from (2) and (3) that the positive gradient in ρ p at the wall of the flux tube will induce a positive gradient in ρ n in the same range of radii, which lowers the energy of the system. This is exactly what we see in Fig. 2, where the dashed curve showing the perturbation to ρ n has a positive slope in the range of radii where the solid curve (ρ p ) has the largest positive slope. On either side of that region it has a negative slope, as it returns to its unperturbed value. For negative σ the effect is reversed: the dash-dotted curve shows ρ n having a negative slope where ρ p has the largest positive slope.
We therefore expect that in the presence of a gradient coupling, the correction to the energy of a flux tube has a dominant core-perimeter term proportional to √ n, The energy correction is negative and quadratic in σ for small σ (see Sec. IV A 4), so the interaction energy parameter M σ is proportional to σ 2 ; choosing it to be positive by convention requires s 1 2 to be positive. We do not have an a priori expectation for the sign of s 1 .

Symmetry under change of sign of couplings
It is clear from Figs. 1 and 2 that for couplings β and σ of order 0.5 the modification of the field configuration due to the interaction between the condensates is extremely small, so it is reasonable to treat its effects perturbatively. (At the end of Sec. IV C we will discuss the limit of small neutron condensate, where the perturbative approach becomes questionable.) When we evaluate the perturbative correction to the energy of the flux tube, there is no linear term in β and σ. Such a term would arise from evaluating the β and σ terms from the Hamiltonian in the unperturbed field configuration. But in that configuration the neutron condensate sits at its vacuum value, so both terms evaluate to zero (g = 1, g ′ = 0 in (17)).
We therefore expect the change in the energy of the flux tube to be quadratic in the couplings β and σ. Firstly, this correction must be negative. This is a well-known result from perturbation theory: the second-order correction arises from the change in the configuration in response to the perturbation, which only occurs because it is driven by a resultant lowering of the energy. Secondly, the change in the energy will in general contain β 2 , σ 2 , and βσ terms. This means it will be even in β when σ = 0 and even in σ when β = 0, so we expect M β ∝ β 2 and M σ ∝ σ 2 in Eqs. (20) and (21).
However, if both β and σ are nonzero, then the βσ terms spoil the symmetry of the energy under negation of the couplings. This is clear from Figs. 1 and 2. For example, suppose that as well as non-zero β we have a very small non-zero σ. Now consider sending β → −β. From Fig. 1 we see that this changes the sign of the slope of ρ n in the wall region where ρ p has positive slope. If σ is nonzero then these two configurations will have different energies, since the gradient of ρ n is then coupled to the gradient of ρ p .

B. Energetic stability of flux tubes
In Fig. 3, the energy per flux unit (E n /n) is plotted against n for various values of the Ginzburg-Landau parameters, namely κ, the density coupling β, and the gradient coupling σ. We fixed φ p 2 / φ n 2 = 0.05 (Sec. III A).

No coupling to neutrons
The upper left plot of Fig. 3 shows E (sc) n (κ)/n, the energy per flux quantum when there are no interactions between the neutron and proton pairs. We see that the only possible phases are the standard type I and type II, with a transition at the Bogomolnyi point, κ = 1/ √ 2, where the favored value of n jumps from 1 to infinity.. The lower line (κ just below 1/ √ 2) corresponds to type-I, where the lowest energy/flux is at n = ∞, so flux tubes attract. The upper line (κ just above 1/ √ 2) corresponds to type-II, where the lowest energy/flux is at n = 1, so flux tubes always repel each other. The middle line corresponds to the transition point (κ = 1/ √ 2), where there is no interaction between flux tubes [25]. Our numerical results are consistent with the expected form (19): when δκ > 0 the asymptotic value of E n /n is increased, and E n /n rises monotonically towards that asymptotic value, and conversely when δκ < 0 the asymptotic value of E n /n is decreased, and E n /n falls monotonically towards that asymptotic value. It is clear that c 1 2 in (19) must be positive to obtain this behavior at large n. From fits to our numerical calculations we find that c 1 is always positive, so it "fights against" the leading c 1 2 / √ n term, but for all n 1 it is overwhelmed. In fact, we find that (19) gives an excellent fit to our results down to n = 1, without any higher order terms.
In the remaining panels of Fig. 3, we explore the effect of density and gradient couplings between the proton superconductor and the neutron superfluid.

Density coupling to neutrons
The upper right panel of Fig. 3 shows the effect of a density coupling between the condensates. From (19) and (20) we expect The first point to notice is that the density coupling shifts the critical κ to a larger value. The transition between type-I and type-II occurs when the asymptotic behavior at large n changes from rising to falling, i.e. when the coefficient of the 1/ √ n term changes sign. This occurs for some positive value of δκ which rises as β 2 because M , M β , b 1 2 , and c 1 2 are all positive, and M β ∝ β 2 when σ = 0 (Sec. IV A 4). Thus in the upper right panel of Fig. 3 we had to increase κ from around 0.707 to around 0.818 in order to find the transition.
The other important point is the presence of a minimum in E n /n when κ is just above the new type-I/type-II boundary, indicating that the favored value of n may be neither 1 (standard type-II) nor infinity (type-I) but some intermediate value. This is consistent with (22), as long as we assume that the coefficient b 1 from (20) is either positive, or negative and of sufficiently small magnitude, so that the 1/n term in (22) has a positive coefficient (recall that M β , M , and c 1 are all positive, and δκ is also positive in this region). The minimum will then arise from competition between the positive 1/n term, which dominates at smaller n, giving a negative slope, and the 1/ √ n term which has a negative coefficient (because δκ is just above the new critical value) and dominates at larger n giving a positive slope. However, as δκ is reduced the negative coefficient of 1/ √ n becomes smaller and smaller, and the minimum moves out to arbitrarily large n, so the energetically favored value of n does not jump suddenly from 1 to ∞ as in the standard case, but increases in steps from 1 to infinity as we lower κ through a range of values down to the new critical value. This creates an infinite number of "type-II(n)" phases, each with a different flux in the favored flux tube, and when that flux becomes infinite the superconductor becomes type-I. This behavior is seen in our numerical results (Fig. 4).

Gradient coupling to neutrons
The lower left panel of Fig. 3 shows the effect of a gradient interaction with the superfluid. From (19) and (20) we expect Here we see that the gradient coupling shifts the critical κ to a smaller value. The transition between type-I and type-II occurs when the coefficient of the 1/ √ n term changes sign, which in this case happens for small negative δκ, which is proportional to −σ 2 because M , M σ , s 1 2 , and c 1 2 are all positive, and M σ ∝ σ 2 when β = 0 (Sec. IV A 4). The other important feature of this plot is the presence of a maximum in E n /n when κ is close to the type-I/type-II boundary. This is consistent with (24), as long as we assume that the coefficient s 1 from (21) is either positive, or negative and of sufficiently small magnitude, so that the 1/n term in (24) has a negative coefficient. The maximum will then arise from competition between the negative 1/n term, which dominates at smaller n, giving a positive slope, and the 1/ √ n term, which dominates at larger n giving a negative slope. The presence of this maximum allows for the possibility of metastable flux configurations. If we scan down in κ, we start in a type-II region where E n /n has its minimum at n = 1 and rises monotonically with n. But at some point a metastable minimum at n = ∞ appears, which drops to become degenerate with the minimum at n = 1. At this point there is a first-order transition: at the critical field, n = 1 flux tubes would co-exist with macroscopic normal regions (i.e. flux tubes with n = ∞) but not with flux tubes of intermediate size. Reducing κ further, the n = 1 flux tube becomes energetically metastable, and finally unstable.

Density and gradient coupling to neutrons
The lower right panel of Fig. 3 shows the effect of a combination of gradient and density interactions. As κ is decreased, a metastable energy minimum emerges at finite n; it drops and becomes a new global minimum at n = n * , yielding a sharp transition from n = 1 type-II to n = n * type-II. As κ is reduced further the favored number of flux quanta in a flux tube rises in integer steps from n * to infinity, at which point the superconductor becomes type-I.  1, but no gradient coupling (σ = 0). The right panel is for gradient coupling σ = 0.1, but no density coupling (β = 0). In both cases, we see that the type-I/type-II transition converges to κ = 1/ √ 2 as the neutron condensate disappears. For the case of a density coupling, as the neutron condensate decreases, the type-I/type-II boundary changes at φn 2 / φp 2 ∼ 10 from a narrow region of type-II(n) phase bands (thick line) to wider metastable regions.

C. Phase diagrams
Figures 4-7 illustrate the additional structure in the phase diagram of the superconductor induced by the couplings to a superfluid. Each diagram is a two-dimensional slice through the parameter space. Figure 4 shows the consequences of a density coupling β between the superfluid and superconductor. We see that the density coupling, irrespective of its sign, favors type-I superconductivity, pushing the the critical κ for the type-I/type-II transition up to higher values, forming a parabolic phase boundary in the β-κ plane, as expected from (23). This can be thought of as arising from the fact that nonzero β lowers the energy per flux of the core of large flux tubes (see (22)), which favors type-I superconductivity.
In the right panel we zoom in on the transition line near β = 0.5 to show the substructure in the phase transition region that is invisibly small in the left panel. As one would expect from our discussion of Figure 3 (upper right panel), on the type-II side of the transition there is a series of bands distinguished by the number of flux quanta n in the energetically favored flux tube. "Type-II (n = 1)" is the standard type-II superconductor. With decreasing κ we find transitions to Type-II (n = 2), Type-II (n = 3), and on up to n = ∞ which is a type-I superconductor.
In Figure 5 we show the consequences of a gradient coupling σ between the superfluid and superconductor. We see that the gradient coupling, irrespective of its sign, favors type-II superconductivity, pushing the critical κ for the type-I/type-II transition down to lower values, forming an inverted parabolic phase boundary in the σ-κ plane, as expected from (25). It also makes the phase transition first order, with spinodal lines where the unfavored phase becomes metastable. Both these effects arise from the lowering of the energy of the wall of the vortex, as explained in Sec. IV B 3.
In Figure 6 we show phase diagrams for the combination of both density and gradient couplings, fixing σ = 0.5 and varying β. As discussed in Sec. IV A 4, we expect that when σ = 0 the β → −β symmetry is now broken. In the right panel we magnify the transition region near β = 0.5, illustrating that on the type-II side as κ decreases the number of flux quanta in the favored flux tube jumps from 1 to a finite value n = 5, and then there is a sequence of bands in which n rises, reaching infinity when the superconductor becomes type I. This is the expected behavior, based on our discussion in Sec. IV B 4.
Finally, in Figure 7, we anticipate one direction in which this work could be extended, by exploring the consequences of varying the ratio of the superfluid density to the superconductor density, which up to now was fixed to φ n 2 / φ p 2 = 20, an appropriate value for neutral beta-equilibrated nuclear matter, of the type we expect to find inside neutron stars. Figure 7 shows phase diagrams in the plane of κ and φ n 2 / φ p 2 for a system with a density coupling (left panel) and with a gradient coupling (right panel).
For the case of a density coupling we use a negative value of the coupling, because this corresponds to an attractive interaction, which gives smooth behavior in the limit where the neutron condensate disappears, φ n 2 / φ p 2 → 0. As is clear from the plot, the type-I/type-II transition then converges to the standard value for a single-component superconductor, κ = 1/ √ 2. For a repulsive interaction, the φ n 2 / φ p 2 → 0 limit is singular: we discuss this in more detail below. It is interesting to note that the effects of the density coupling change dramatically with the relative densities of the neutrons and protons. At φ n 2 / φ p 2 10 the density coupling produces a thin region of multi-fluxquantum "type-II(n)" phases, as was illustrated in Fig. 4. But for lower values, it has a similar effect to a gradient coupling, inducing metastable regions on either side of the type-I/type-II boundary. This should be understandable in terms of the dependence of the coefficients b 1 2 and b 1 (Eqn. (20)) on φ n 2 / φ p 2 . In Sec. IV B 2 we argued that if b 1 is large enough then the E n /n curve has a minimum at finite n, yielding a type-II(n) phase. We conjecture that as φ n 2 / φ p 2 gets smaller, b 1 becomes sufficiently negative that this is no longer the case, and instead there is a maximum, leading to metastability of the n = 0 and n = ∞ states in spinodal regions around the type-I/type-II boundary. This is a topic for future investigation.
For the case of a gradient coupling (right panel of Fig. 7) the effects of varying φ n 2 / φ p 2 are less dramatic. It is interesting that, as for a density coupling, the variation is non-monotonic. Again, we conjecture that this could be understood in terms of variation of the coefficients s 1 2 and s 1 (Eqn. (21)) with φ n 2 / φ p 2 . As the superfluid density drops to zero, its effects become negligible, and the critical value of κ converges towards 1/ √ 2 as one would expect. Finally, we discuss the singularity of the φ n 2 / φ p 2 → 0 limit for a positive β, i.e. a repulsive density coupling between the neutron and proton condensates. From (6) we see that the expectation value of the neutron condensate is φ n + 1 2 β( φ p − φ p ), so far from the flux tube, where φ p is φ p , it is φ n . But in the core of the condensate it is larger (there is less proton condensate to repel it). In fact, even if the parameter φ n 2 were zero or slightly negative, there would be a positive neutron condensate in the core of the flux tube. This shows that for positive β the neutrons do not decouple and become irrelevant in the limit φ n → 0. We note two consequences of this. Firstly, for small φ n the β → −β symmetry discussed in Sec. IV A 4 is no longer present, because the effect of the flux tube on the neutron condensate is no longer a small perturbation. Secondly, in a system where φ n 2 is small and negative (i.e. the neutrons just barely fail to condense in the presence of the proton condensate) flux tubes could have superfluid cores, which is another topic that we leave for future investigation.

V. CONCLUSION
We conclude that coupling a superconductor to a co-existing superfluid causes significant modification of the energetics of the flux tubes. On the basis of calculations restricted to the cylindrical geometry of n-quantum flux tubes, we conclude that a coupling between the densities of the condensates shifts the type-I/type-II boundary to larger κ, and, if the superfluid density is high enough, appears to create an infinite number of new "type-II(n)" phases whose most stable flux tubes contain multiples of the basic flux quantum. A gradient coupling between the condensates leads to metastable regions surrounding the transition between type-I and type-II superconductivity.
As discussed in Section II, our calculation corresponds to comparing the energy at zero and infinite separation of flux tubes with varying numbers of flux quanta. This leaves open the possibility that there might be additional minima at finite separation. It is therefore possible that in parts of the phase diagram there might be a different phase from the ones we identify, namely an alternative type of type-II superconductor in which the spacing between flux tubes is fixed by the microscopic physics rather than by the strength of the applied field. To resolve this question will require calculation of the free energy of a pair of flux tubes at arbitrary separation. Such calculations have been performed for large separation [28][29][30], and by perturbing about the Bogomolnyi point [17] and by numerical computation [31]. In particular, the numerical methods that have been used recently to follow the interaction and annihilation or vortex-antivortex pairs [32] would be readily applicable to the simpler time-independent calculation of the interaction potential of flux tubes. Another natural generalization of our calculation would be to allow for non-s-wave pairing, such as the 3 P 2 pairing that is believed to occur in the neutron superfluid in the core of a neutron star.
Our results add another example to the class of two-component Ginzburg-Landau models with non-standard superconducting behavior. Previous work in this area includes the SO(5) model of high-temperature superconductivity, which has flux tubes described by a two-component GL model, where each component carries a different U (1) charge, and only one of them condenses in the vacuum [13]. Another example is the case of a two-component GL model where both components have electric charge, very different mass, and nearly the same Fermi energy. This system was found to have non-monotonic E(n)/n and intermediate minima in the interaction potential [7].
The exotic phenomena that we predict are localized to the region around the type-I/type-II transition, so they may not turn out to be relevant for the inner core of a neutron star, which is believed to be well inside the type-II regime [1]. However, given the extremely impressive recent progress in creating exotic systems such as multi-component superfluids of trapped cold atoms, it seems quite conceivable that a material that is both a superconductor and a superfluid might be created in the laboratory, and could be studied under controlled conditions. Our results would be directly relevant to such a material.