Superconducting Phases in Neutron Star Cores

: Using a phenomenological Ginzburg–Landau model that includes entrainment, we identify the possible ground states for the neutron and proton condensates in the core of a neutron star, as a function of magnetic ﬁeld strength. Combining analytical and numerical techniques, we ﬁnd that much of the outer core is likely to be a “type-1.5” superconductor (instead of a type-II superconductor as often assumed), in which magnetic ﬂux is distributed inhomogeneously, with bundles of magnetic ﬂuxtubes separated by ﬂux-free Meissner regions. We provide an approximate criterion to determine the transition between this type-1.5 phase and the type-I region in the inner core. We also show that bundles of ﬂuxtubes can coexist with non-superconducting regions, but only in a small part of the parameter space.


Introduction
In the core of a neutron star, neutrons and protons are both expected to form condensates via Cooper pairing, leading to a neutron superfluid and a proton superconductor, respectively. Within the proton superconductor magnetic flux is quantized into microscopic filaments called fluxtubes, and the macroscopic dynamics of the star's core magnetic field depends on the microscopic dynamics of these fluxtubes. In the conventional picture [1][2][3], the inner core is a type-I superconductor, meaning that fluxtubes are mutually attractive and therefore coalesce into macroscopic patches of magnetic flux, whereas the outer core is a type-II superconductor, meaning that fluxtubes are mutually repulsive and therefore form a regular array. The transition between these two superconducting regimes occurs where the ratio of the London penetration length, λ, and proton coherence length, ξ p , the so-called Ginzburg-Landau parameter κ ≡ λ/ξ p , takes the value κ = 1/ √ 2. The inner and outer core are therefore defined to be the regions in which κ < 1/ √ 2 and κ > 1/ √ 2, respectively. However, the microscopic behavior of fluxtubes in the core is affected by the coupling between the proton and neutron condensates, and in particular by their mutual entrainment, which could significantly change the location of this transition, and might even result in entirely different types of superconductivity. For example, Buckley et al. [4,5] showed that strong density coupling between the condensates could produce type-I superconductivity throughout the whole core, although only if the coupling is much stronger than is generally expected [6]. A more detailed study by Alford and Good [7], incorporating density and density-gradient coupling, found that in some cases the transition from type-I to type-II is mediated by domains of "type-II(n)" superconductivity, wherein each fluxtube carries n magnetic flux quanta. Subsequently, Haber and Schmitt [8] argued that these type-II(n) fluxtubes are generally unstable, and instead there is a regime of "type-1.5" superconductivity, in which the fluxtubes form bundles with a preferred separation. Such unconventional behavior has also been discussed in the context of terrestrial multi-band superconductors [9][10][11][12][13]. We will show later that inhomogeneous flux distributions arise in our model under a wide range of parameter conditions, as a result of coupling between the two superconductors, even well below the condensation temperature [44]. We will therefore take a similar approach to model the neutron star core.
In what follows, we normalize the order parameters such that |ψ p | 2 , for example, is the number density of proton Cooper pairs, and so ρ p = 2m p |ψ p | 2 is the mass density of the proton condensate, where m p is the mass of a proton. The order parameters are therefore two-particle mean-field wave functions for the condensates. As noted earlier, the core temperature in mature neutron stars lies far below the critical temperature for superconductivity, so in the absence of magnetic flux practically all of the proton matter would reside in the condensed state. In the presence of magnetic flux, however, there will be normal proton matter present in the cores of fluxtubes and in any non-superconducting regions, where the proton condensate is absent. We are concerned here only with the ground state for the condensates, wherein interactions with any "normal" components of the core, which include electrons, thermal excitations and normal protons and neutrons, are suppressed. Hence, we can safely disregard the normal matter in what follows.

Entrainment and Local Phase Invariance
The residual strong interaction between neutrons and protons in the core leads to a non-dissipative drag between their condensates known as entrainment [39,45]. To illustrate the dynamical effect of entrainment, we temporarily neglect any coupling to the magnetic field by setting A = 0; the magnetic field will be reintroduced later by invoking gauge invariance. The hydrodynamical momenta of the condensates are then proportional to the gradient of the phases of the order parameters, i.e., where V p and V n are the superfluid velocities. In the presence of entrainment, the velocitydependent terms in the free-energy density, F vel , of the condensates must take the form where ρ p and ρ n are the true mass densities of the condensates and the coefficient ρ pn , which determines the strength of entrainment, is generally negative [46]. This form of the free energy-with an interaction term that depends only on the relative velocity-is necessary to ensure Galilean invariance [31] (as well as the more general constraint of "local phase invariance" [30]). Equation (2) demonstrates that the primary effect of entrainment is to disfavor any relative flow between the two condensates by imposing an energetic penalty wherever V p = V n . A more subtle but equally important consequence is that the condensates' hydrodynamical momenta, given in Equation (1), are no longer proportional to their mass fluxes, defined as ∂F vel /∂V x for x ∈ {p, n}. This has significant consequences for the structure of fluxtubes and vortices, and for their mutual interactions [39], but in the present work we are concerned only with fluxtubes, i.e., topological defects in the proton condensate. For this reason, we will neglect the star's rotation, meaning that no neutron vortices are present in the ground state. While neglecting rotation is certainly a simplification, it is justified when deriving a microscale model of the neutron star interior, in which the proton fluxtube density is many orders of magnitude larger than that of the neutron vortices, e.g., [47].
Within the Ginzburg-Landau mean-field framework, entrainment first enters the free-energy density at fourth order in the order parameters, and at second order in their derivatives [39]. The most general such term that satisfies global U(1) symmetry in each condensate is a linear combination of the quantities |ψ x | 2 |∇ψ y | 2 , ψ x ψ y ∇ψ x · ∇ψ y , ψ x ψ y ∇ψ x · ∇ψ y , ψ x ψ y ∇ψ x · ∇ψ y , where x, y ∈ {p, n} and indicates the complex conjugate. With the additional constraint of Galilean invariance (2), the most general form of the entrainment term is found to be which includes four real, independent parameters h 1 , . . . , h 4 . In terms of the superfluid densities and velocities, we can write this as So by comparison with Equation (2) the entrainment coefficient is The parameters h 2 , h 3 , h 4 only provide density-gradient coupling, and the simplest model of entrainment would therefore set these parameters to zero. On scales much larger than the fluxtube cores the superfluid densities are approximately constant, and so these terms will have negligible effect. However, we will show later that the density-gradient terms play a significant role in the transition between type-I and type-II superconductivity, and therefore must be included in the construction of phase diagrams.

Connection with Previous Work
Our general expression (5) for the entrainment energy differs from corresponding expressions found in some previous works. To highlight the differences, we note that the velocity contributions to the free-energy density, given by Equation (2), can equivalently be expressed as where ρ pp ≡ ρ p − ρ pn and ρ nn ≡ ρ n − ρ pn represent "effective" proton and neutron mass densities. Therefore, on scales much larger than the vortex and fluxtube cores, for which the superfluid densities are approximately constant, entrainment can be described by including in the free-energy density a term proportional to V p · V n , and renormalizing the proton and neutron masses accordingly. However, in a microscale model that correctly includes density variations in the fluxtube cores, the dependence of the coefficients ρ xy on the condensate densities must be chosen carefully to preserve Galilean invariance [31,37] and additional density-gradient coupling terms have to be included. A number of previous studies, e.g., [39,43] do not treat entrainment on small scales consistently, because their entrainment interactions are incompatible with Equations (4) and (5). Haber and Schmitt [8] have introduced a relativistic model of density and derivative couplings that in the nonrelativistic limit is also incompatible with our Equation (5), cf. their Equation (5). The model of Alford and Good [7] is compatible with Equations (4) and (5), but it only includes the h 2 term. Hence their model actually has no entrainment at all, i.e., ρ pn = 0. This appears to be an oversight on their part, because they chose the value for h 2 based on prior estimates of ρ pn . Finally, the model of Kobyakov [37] has a similar but subtly different form to Equation (4), because it takes the entrainment term to be Im{|ψ n | 2 ψ p ∇ψ p − |ψ p | 2 ψ n ∇ψ n } 2 |ψ p | 2 |ψ n | 2 .
Although this quantity is Galilean invariant, it cannot be obtained from products of the order parameters, their conjugates and derivatives, and therefore cannot arise in our mean-field formalism. For the same reason, our model does not include a term proportional to ∇|ψ p | · ∇|ψ n |, unlike the model of Kobyakov [37]. In what follows, we highlight where our results reproduce those of earlier studies, in appropriate limits.

The Free Energy Density
The total free-energy density in our model is obtained by adding the entrainment terms (4) to the usual free energy of a two-component superfluid, and introducing the magnetic vector potential A by minimal coupling. To reduce the number of parameters in the model, in what follows we assume that the protons and neutrons have equal masses, m p = m n = m u , and that h 3 = h 4 , which means that the entrainment interaction (4) is symmetric in the condensates. Using Gaussian c.g.s. units, the free energy density in its most compact form can then be expressed as where we have assumed that the proton Cooper pairs have charge 2e. The coefficients g pp and g nn measure the self-repulsion of the condensates, and g pn measures their mutual repulsion or attraction. In the absence of magnetic field (i.e., if B ≡ ∇ × A = 0) we expect the ground state to be a uniform mixture of proton and neutron condensates, with positionindependent densities |ψ p | 2 = n p /2 and |ψ n | 2 = n n /2. Therefore the parameters n p and n n represent the expected number density of superfluid protons and neutrons, respectively, in the absence of magnetic field. For reasons explained earlier, in a mature neutron star these will have values very close to the total (i.e., normal plus condensate) particle number densities. We will hereafter refer to this state, with uniform condensate densities and vanishing magnetic field, as the Meissner state, and we note that it has F = 0, according to Equation (9). However, if the mutual attraction/repulsion between the condensates is too strong, such that g 2 pn > g pp g nn , then this Meissner state becomes unstable. Such behavior is not expected in the neutron star core, where the two condensates are believed to be only weakly attractive [6], and so in what follows we will always assume that g 2 pn < g pp g nn . Moreover, we expect a neutron condensate to be present even in non-superconducting regions, where ψ p = 0, and this implies the further restriction g pn > −g nn n n /n p . The consequences of violating these restrictions have been discussed in detail by Haber and Schmitt [8], for instance. In most studies of two-component condensates, the coefficient g pn represents the principle interaction between the two components, e.g., [48][49][50][51], and its effect on superconductivity in the neutron star core has been studied extensively [7,8,37]. In the present work, however, our main focus is on the effect of entrainment and other higher-order coupling terms. In the numerical results we present later, we therefore take g pn = 0, and instead study the effect of the h i parameters on the superconductor. For completeness, and to facilitate comparison with earlier studies, we retain g pn in all of our analytical results.
In the absence of coupling between the condensates (i.e., for g pn = 0 and h i = 0) the "bare" coherence lengths are defined as ξ p ≡ 2m u g pp n p and ξ n ≡ 2m u g nn n n , (10) and the "bare" London length is defined as We will describe later how the effective coherence lengths and London length are modified by the coupling between the condensates, including their mutual entrainment.
In order to simplify the mathematical model, we now nondimensionalize the free-energy density (9) by measuring ψ p and ψ n in units of n p /2 and √ n n /2, respectively, lengths in units of ξ p , A in units of c/(2eξ p ), and the coupling coefficients h i in units of g pp ξ 2 p . To improve the readability of our equations, we avoid introducing specific notation for dimensionless quantities and instead point out that, from here on, all parameters refer to dimensionless quantities. The dimensionless free-energy density is then, in units of g pp n 2 p /4, where we have defined the following parameters: Note that κ is equivalent to our dimensionless "bare" London length.

The Helmholtz and Gibbs Free Energies
We now seek the ground state for this system in the presence of an imposed magnetic field. There are two distinct thought-experiments that can be considered. In the first experiment, we control the magnetic flux density, B, by imposing a mean or net magnetic flux, and minimize the Helmholtz free energy, where the angled brackets represent some kind of integral over our physical domain, which could be finite or infinite. This experiment closely approximates the conditions in the core of a neutron star, which becomes superconducting as the star cools in the presence of a pre-existing magnetic flux. However, as we will discuss below, the ground state under these conditions can be inhomogeneous, i.e., macroscopic domains of distinct physical behavior can appear. For conceptual convenience, we can consider an alternative experiment in which the system is coupled to a thermodynamic external magnetic field, H, by minimizing the dimensionless Gibbs free energy, In an unbounded domain, the ground state in this experiment is guaranteed to be homogeneous, and hence the phase diagram is generally simpler. For later reference, we present in Figure 1 the phase diagrams for a single-component Ginzburg-Landau superconductor, i.e., we discuss its state as a function of the Ginzburg-Landau parameter, κ. (For more details, we refer the reader to standard textbooks on superconductivity, e.g., Tinkham [52]). For κ < 1/ √ 2, we have a type-I superconductor; when H is used as the control parameter, there is a first-order transition between the Meissner state (with B = 0) and the non-superconducting state (with B = H) at the critical value |H| = H c = 1/( √ 2κ) in our dimensionless units. When the mean magnetic flux, B, is used as the control parameter, this discontinuity resolves into an intermediate phase for 0 < B < H c , in which Meissner regions alternate with non-superconducting ones. For κ > 1/ √ 2, on the other hand, we have a type-II superconductor; for H c1 < |H| < H c2 the magnetic flux organizes into a hexagonal lattice of discrete fluxtubes. The transitions at the lower critical field, H c1 , and the upper critical field, H c2 , are both second-order, because fluxtubes appear with infinite separation at |H| = H c1 , and the superconductor density becomes vanishingly small at |H| = H c2 . In our dimensionless units, H c2 = 1 and H c1 = F ∞ /(4πκ 2 ), where F ∞ is the energy per unit length of a single fluxtube. When B is the control parameter, there is a similar second-order transition at B = H c2 , and a first-order transition between the intermediate and fluxtube states at κ = 1/ √ 2. For our two-component system, we anticipate that the phase diagram will be more complicated than that shown in Figure 1. In particular, Haber and Schmitt [8] have argued that the upper and lower transitions to and from the fluxtube state can become first-order in some cases, occurring at |H| = H c1 < H c1 and |H| = H c2 > H c2 , respectively. In that case, in the experiment with an imposed mean magnetic flux, B, the ground state can feature an irregular array of fluxtubes, even in an unbounded domain. Some aspects of this phase space can be determined analytically, as we describe in Section 4. However, in order to produce a complete phase diagram, it is necessary to solve the Euler-Lagrange equations arising from one of the functionals F or G numerically, as we describe in the next section.

The Numerical Model
Whether we choose to work with the Helmholtz free energy, F , or with the Gibbs free energy, G, we obtain the same system of Euler-Lagrange equations: However, the appropriate boundary conditions for these two experiments are different, and also depend on the particular size and shape chosen for the domain. Without loss of generality, we will assume from here on that the magnetic field is oriented in the z-direction, and that all variables are independent of z; so our domain will be some region within the xy-plane. We solve a discretized version of Equations (16)- (18), which are obtained by minimizing a discrete approximation to the free energy on a regular grid in x and y. The gauge field is included via a Peierls substitution, in order to maintain gauge invariance. The equations are solved using a simple relaxation method, and the grid resolution is repeatedly refined until a sufficient level of accuracy has been obtained. Additional details on the numerical algorithm can be found in Appendix A.
In the present study, we are not interested in the effect of physical boundaries on the phase diagram, and so we would ideally use an infinite domain, but for numerical calculations the domain must of course be finite. Moreover, we cannot use periodic boundary conditions, because in the presence of fluxtubes neither A nor ψ p is spatially periodic. Instead, we must use quasi-periodic boundary conditions [53], which involves specifying not only the size of the domain, L x × L y , say, but also the number, N, of magnetic flux quanta within the domain. Working in the symmetric gauge, the quasi-periodic boundary conditions for our dimensionless variables take the form where L represents either of the translational symmetries (L x , 0) or (0, L y ). The effect of these boundary conditions is to impose a total magnetic flux of 2πN within our domain (in our dimensionless units, the quantum of magnetic flux is 2π), and thus a mean magnetic flux of B = 2πN/(L x L y ). By solving Equations (16)- (18) subject to these boundary conditions, in domains of varying size, we can thus obtain the Helmholtz free energy, F , as a function of the mean magnetic flux. We note that the choice of domain aspect ratio affects the configuration of any fluxtube array that forms. In particular, we can impose either a square or a hexagonal lattice symmetry by using the following domain shapes: • for a square lattice, we take N = 1 and L x /L y = 1; • for a hexagonal lattice, we take N = 2 and L x /L y = √ 3.
In order to directly compare these two cases, we calculate the Helmholtz free energy per magnetic flux quantum per unit length: As an example, in Figure 2 we plot F as a function of the area per magnetic flux quantum, for one particular set of parameters, whose motivation is explained later, in Section 5. Note that in the case of a fluxtube lattice, a corresponds to the area of a single Wigner-Seitz cell. This plot was produced by finding the minimum value of F for both square and hexagonal lattices for domains of various sizes. We also plot the energy in the non-superconducting state, which has ψ p = 0 and a uniform magnetic field B = (0, 0, 2π/a), and is known analytically (see Section 4.1).
The energy in both the square (long-dashed, cyan) and hexagonal (solid, blue) lattice states matches smoothly onto the energy of the non-superconducting state (short-dashed, purple) at a 12.9 (region enlarged in the upper inset), and both converge to the same finite value as a → ∞. The dotted gray lines indicate the lower convex envelope, plotted separately in Figure 3, which is the true ground state in an unbounded domain. The two red dots indicate the values for the two simulations in Figure 4. We show an enlarged view of the left point in the lower inset.
However, as discussed in the previous section, in some cases the true ground state might be inhomogeneous, if this allows the average energy to be lower than that of any homogeneous state. Suppose, for example, that a fraction, φ, of the total magnetic flux is contained in regions with a = a 1 and F = F 1 , while the rest is in regions with a = a 2 and F = F 2 . In that case, the overall values of a and F are given by the lever rule: In this way, an energy F that is lower than F (a) can be achieved in any range of a for which the function F (a) is not convex. In fact, the true ground state in an unbounded domain is given by the lower convex envelope of all the homogeneous states, which is indicated by the dotted gray lines in Figure 2. We will use the notation F g (a) to refer to the true ground-state energy as a function of the area a. For this particular case, there are four distinct behaviors seen across the full range of a: • for 0 < a 12.7 the ground state is non-superconducting; • for 12.7 a 18.5 the ground state is a mixture of non-superconductor and a hexagonal fluxtube lattice; • for 18.5 a 26 the ground state is a hexagonal fluxtube lattice; • for a 26 the ground state is a mixture of a hexagonal fluxtube lattice and the Meissner state.
Once the function F g (a) is known, it is straightforward to also determine the minimum Gibbs energy as a function of H. In fact, the mean Gibbs energy density is and the minimum of G over all a can be interpreted graphically from the plots in Figure 2.
is a convex and monotonically decreasing function, each point on the curve F g (a) corresponds to a ground state with energy G = F g (a), and the corresponding value of |H| can be found by extrapolating the tangent line up to the F -axis, as shown in Figure 3.  We emphasize that the function F g represents the minimum free energy only for a hypothetical unbounded domain, free from any geometrical constraints. In any simulation with a finite domain size, the free energy in the ground state will generally exceed this value. Nevertheless, by using a large enough computational domain, and choosing values of a within the appropriate ranges, we can obtain examples of the inhomogeneous ground states described above. Figure 4 shows two such examples, with a = 14.5 and a = 52. These values of a were chosen so that in each case approximately half of the domain contains a hexagonal lattice. As shown in Figure 2, in each case the free energy is lower than that of a pure lattice, but still significantly higher than for the true ground state in an unbounded domain.  1. The second panel shows a case with N = 14 magnetic flux quanta, and a (dimensionless) area of aN = 52 × 14; this is a mixture of Meissner state and hexagonal fluxtube lattice. In both cases the aspect ratio is √ 3, which means that a pure hexagonal lattice is a possible state of the system, but is not the ground state.

Phase Diagrams
Using the procedure described above, we can determine the superconducting phase transitions that occur when B or |H| is used as the control parameter. The results shown in Figure 2 demonstrate that these transitions can be qualitatively different from those seen in a single-component Ginzburg-Landau superconductor. In order to directly compare the resulting superconducting phases with those shown previously in Figure 1, we produce phase diagrams in which κ is used as the independent parameter and all other parameters are fixed to the same values used in Figure 2. The result is shown in Figure 5. For κ 1.42 we have type-I superconductivity, and for κ 5.36 we have type-II superconductivity. However, within the range 1.42 κ 5.36 we have type-1.5 superconductivity: in the experiment with |H|, the transition from the Meissner to the fluxtube state at H c1 is first order, because the lattice first appears with a finite separation between fluxtubes; in the experiment with B, there is a critical value B = B c1 below which bundles of fluxtubes form alongside flux-free Meissner regions. This behavior results from the fact that fluxtubes are mutually attractive at large separation distances, and mutually repulsive at small separation distances, i.e., they have a preferred separation, as indicated by the minimum in the curve F (a) seen in Figure 2.
Over the narrow range 1.42 κ 1.58 an additional feature is present: in the experiment with |H|, the transition from the fluxtube state to the non-superconducting state is first order at H c2 , as the proton condensate density vanishes discontinuously; in the experiment with B, there is a critical value B = B c2 above which bundles of fluxtubes form alongside non-superconducting regions. This behavior can also be understood intuitively in terms of interactions between neighboring fluxtubes: for this range of parameters, the tendency towards a preferred separation distance is so strong, and the dip in the F (a) curve so pronounced, that it becomes energetically favorable to form a mixture of fluxtubes and non-superconducting regions, rather than a periodic lattice with the "wrong" separation distance.
The phase diagrams shown in Figure 5 resemble those hypothesized by Haber and Schmitt [8], and demonstrate that type-1.5 superconductivity can arise over a significant range of the Ginzburg-Landau parameter κ, when coupling between the proton and neutron condensates is present. However, to determine whether the behavior seen for this particular choice of coupling parameters is generic, we need to understand these phase transitions in more detail. In the following section we therefore seek analytical expressions for the locations of the various transitions in phase space, in order to determine which couplings between the two condensates give rise to inhomogeneous ground states.

Phase Transitions with H
In this section we will analyze in detail the phase transitions in the experiment involving the external magnetic field, H. As is clear from the previous section, the existence of first-order transitions at the lower and upper critical fields, H c1 and H c2 , results from the non-convexity of the free energy F (a) in the pure lattice state. To describe these phase transitions in general requires quite detailed knowledge of the function F g (a), which can only be determined with a 2D numerical model. However, some important features of the superconducting phase diagram can be determined either analytically, or from knowledge of the structure of a single fluxtube. This allows the phase diagram to be constructed more efficiently, because it reduces reliance on the 2D code, and it provides some physical insight into the origins of these first-order phase transitions. In the following sections, we describe those features of the function F (a) that can be determined analytically or semi-analytically.

The Critical Field, H c
The two simplest solutions of the Euler-Lagrange Equations (16)- (18) are the Meissner (flux-free) state, which has |ψ p | = |ψ n | = 1 and B = 0, and the non-superconducting state, which has |ψ p | = 0, |ψ n | = √ 1 + α/R 2 and a uniform magnetic flux density B with |B| = 2π/a. As discussed earlier, these two states can only be realized if the mutual repulsion/attraction g pn , and thus α, between the two condensates is sufficiently weak. Expressed in terms of the dimensionless parameter α, this implies the conditions α 2 < R 2 and α > −R 2 . We will assume from here on that both of these are satisfied.
To determine the thermodynamical critical field, H c , we equate the energy of the Meissner state with that of the non-superconducting state. By construction, the free-energy density (12) in the Meissner state is F = 0, and, hence, its mean Gibbs energy is G = 0. For a type-I system, the transition to the non-superconducting state therefore occurs when the corresponding energy density G becomes negative.
The purple, short-dashed curve in Figure 2 represents the free energy per unit length per quantum of flux in the non-superconducting state. From Equation (12), we obtain Substituting into Equation (25), we then find that the minimum mean Gibbs energy density, G, is achieved for a = 2π/|H|, as expected, and that this minimum is Thus, assuming that we have a type-I superconductor, there is a first-order transition between the Meissner and the non-superconducting state at This defines the critical field, H c , in agreement with previous studies [8,37].

The Lower Critical Field, H c1 vs. H c1
As shown in Section 3.3, the lower critical field is a first-order phase transition if the Helmholtz free energy per unit length per flux quantum F (a) in the pure lattice state has a minimum at some finite value of a, since then a fluxtube lattice forms with finite separation. If this minimum is F min , say, then we have H c1 = F min /(4πκ 2 ). If, on the other hand, F is a monotonically decreasing function of a in the lattice state, then there is a second-order transition at H c1 = F ∞ /(4πκ 2 ), as in the case of a single-component superconductor, where In this limit, interactions between the fluxtubes vanish, and F ∞ is equivalent to the energy per unit length of a single fluxtube in an infinite domain. This energy can be computed efficiently by using polar coordinates r, θ centered on the fluxtube core, and seeking solutions of Equations (16)- (18) in the form [7] ψ p = f (r) e iθ , ψ n = g(r) , This ansatz assumes that the fluxtube carries a single quantum of magnetic flux, and that there is no corresponding phase defect in the neutron condensate. It further results in a system of ordinary differential equations that can be solved numerically, yielding the value of F ∞ to high accuracy. However, to determine whether the lower transition is second-order or first-order, we need to know whether the function F (a) tends to F ∞ from above or from below, which is equivalent to asking whether the long-range interaction between fluxtubes is repulsive or attractive. This can be derived rigorously using a method introduced by Kramer [54], which we describe in detail in Appendix B. In fact, the main result can be obtained heuristically by considering the perturbations produced by a fluxtube in the far-field, i.e., at a large distance from its core. It is convenient to work with the real variables The gauge invariance of the free energy density (12) guarantees that it can be rewritten in terms of these variables without loss of generality; for the full expression see Equation (A6). The corresponding Euler-Lagrange equations are then The far-field structure of the fluxtube can be determined by linearizing these about the uniform state with f = g = 1, and V = ∇χ = 0. This leads to the following system: where δ f , δg, δV, and δχ denote the linear perturbations. In the case of a single fluxtube, we are interested in the solution that is axisymmetric and decays at large distance from the origin. We deduce from Equations (38) and (39) that this solution has δχ = 0 and for some constant coefficient V 0 , where K 1 is a modified Bessel function and λ is the effective London length, Note that, compared to the London length λ in the absence of coupling, λ is made smaller by the parameter h 1 , because the effective mass of the protons is made smaller by the entrainment of neutrons [39]. From Equations (36) and (37) we find that, owing to the coupling between the condensates, the fluxtube in the far-field has a double-coherencelength structure, i.e., where K 0 is a zeroth-order modified Bessel function. The coefficients f i , g i and the effective coherence lengths ξ i satisfy the equations This is reminiscent of the fluxtube structure found in two-component superconductors, although in that case the two coherence lengths arise because a fluxtube is a phase defect in both condensates [55]. In our model, fluxtubes are phase defects in the proton condensate only, but the coupling coefficients α and h 2 ultimately achieve a similar effect. In some parameter regimes (as we have indeed seen in Section 3.4), we might therefore expect our model to exhibit type-1.5 superconductivity, which is common in two-component superconductors, and generally occurs when the London length lies between the two coherence lengths [13]. The physical reason is that the electromagnetic interaction between fluxtubes, which decays on the London length, is generally repulsive, whereas the density interactions are generally attractive. Thus type-1.5 superconductivity arises because fluxtubes are mutually attractive at large separation distances, but become mutually repulsive at shorter separations, producing bundles of fluxtubes with a preferred density, as indicated by the minimum in the F (a) curve.
As mentioned earlier, this heuristic argument can be put on a more rigorous basis by calculating the long-range interaction energy between fluxtubes in a lattice configuration, using the method first described by Kramer [54]. A similar method was used by Haber and Schmitt [8] to demonstrate the existence of type-1.5 superconductivity resulting from density and derivative couplings. As we evaluate the interaction energy for a more general case than [8], and thus obtain a different result, we present our analysis in full in Appendix B, although the steps closely follow those of Kramer [54] for a single-component superconductor. Our final result for the interaction energy is where x i is the location of the i-th lattice point, assuming that x 0 = 0. In the absence of any coupling between the two fluids (α = h 1 = h 2 = h 3 = 0) this reduces exactly to the result of Kramer [54]. In principle, we can now use this result to estimate the free energy of a particular lattice state using just the values of the coefficients f i , g i , V 0 , which can themselves be inferred from the nonlinear solution for a single fluxtube. However, this result is only strictly valid if the fluxtubes are very widely separated, and it becomes inaccurate once the fluxtubes are close enough to interact nonlinearly. Nevertheless, we can deduce that, in the asymptotic limit a → ∞, where d ∝ √ a is the lattice constant, and ξ + represents the larger of the two effective coherence lengths, which in practice is always larger than both of the bare coherence lengths. If λ > ξ + / √ 2 then, according to Equation (47), the long-range interaction energy is positive, implying that there is a second-order transition at the lower critical field H c1 , as for a type-II superconductor. Conversely, if λ < ξ + / √ 2 then the interaction energy is negative, implying that the F (a) curve has a minimum at a finite value of a. In that case, there is a first-order transition at the lower critical field H c1 , and we have a type-1.5 superconductor, which confirms the heuristic argument given above. In summary, we can identify the precise value of κ at which H c1 = H c1 (and B c1 = 0) by equating λ , which is given by Equation (41), with ξ + / √ 2, which is independent of κ. In the case shown in Figure 5, this yields the value κ 5.36, in agreement with our numerical results.

The Upper Critical Field, H c2 vs. H c2
As shown in Section 3.3, if the free energy F (a) is convex (and monotonically decreasing), then the transition between the fluxtube lattice state and the non-superconducting state is second-order. This means that the order parameter ψ p vanishes smoothly at the transition point, H c2 , which can thus be determined analytically by considering linear perturbations to the non-superconducting state. Working in the symmetric gauge, the non-superconducting state is given by |ψ p | = 0, |ψ n | = √ 1 + α/R 2 , and A = 1 2 Be z × x, where B = 2π/a is the (uniform) mean magnetic flux. From the linearized version of Equation (18), we find that the perturbation δψ p must satisfy the equation As known from single-component systems, e.g., [52], bounded solutions of this equation first appear when Therefore, the solid (blue) and long-dashed (cyan) lines in Figure 2, as highlighted in the upper inset, meet the short-dashed (purple) line at the point where a = 2π/H c2 . However, if the function F (a) is not convex at this point then the upper transition will occur not at |H| = H c2 , but at a higher value |H| = H c2 , and will be first order. To determine whether the function F (a) is convex at this point, we can seek weakly nonlinear solutions in the vicinity of a = 2π/H c2 , following the method of Abrikosov [56]. We present the details of this calculation in Appendix C. The main result is that the function F (a) becomes non-convex when where x i is the location of the i-th lattice point, as in Section 4.2, and β is the kurtosis of the proton order parameter, Equation (50) is valid for both hexagonal and square lattices of singly-charged fluxtubes. (Singly-charged here means that each fluxtube carries a single quantum of magnetic flux.) Although there are cases in which other fluxtube configurations have a lower free energy (see Appendix C), the point of transition between H c2 and H c2 is always determined by the singly-charged hexagonal lattice, for which β 1.1596 [57]. Thus from Equation (50) we can determine the value of κ at which H c2 , H c2 and B c2 all coincide. In the case shown in Figure 5, this yields the value κ 1.58, in agreement with our numerical results.

Discussion
Combining analytical and numerical techniques, we have determined the ground state for a mixture of proton and neutron superfluid condensates in the neutron star core, in the presence of magnetic flux. The condensates are coupled by density and densitygradient interactions, and by mutual entrainment. Our model extends the phenomenological Ginzburg-Landau framework used in previous works to consistently satisfy Galilean invariance on small scales. In addition to the three homogeneous phases (Meissner, fluxtube lattice, and non-superconducting) that are familiar from single-component superconductors, we have shown that inhomogeneous mixtures of any two of these phases can occur in the ground state, as a result of the coupling between the condensates. The novel mixed states occur because the fluxtubes perturb the neutron condensate density, which can reduce the overall free energy of the coupled system. As described in Section 4.2, coupling between the two condensates always increases the effective coherence length of fluxtubes, and hence under a range of conditions the fluxtubes have a preferred separation distance. Under these conditions, and for relatively weak mean field values (B < B c1 ), we have type-1.5 superconductivity, in which the fluxtubes form bundles with local hexagonal symmetry, with the rest of the domain in a flux-free Meissner state. In cases where the tendency towards fluxtube bundling is sufficiently strong, there may also be a range of mean field values (B c2 < B < H c2 ) for which there is a mixture of fluxtube and non-superconducting regions. This occurs if it is energetically favorable to confine the proton condensate within part of the domain, forming a lattice with the preferred separation, with the remaining magnetic flux concentrated into non-superconducting regions. In all cases, we find that fluxtubes preferentially adopt a singly-charged, hexagonal pattern, even in cases where a periodic lattice is not the ground state.
To relate our results to neutron stars, we must provide estimates for all of the dimensionless parameters in our model. The results presented in Section 3 assumed the following dimensional values: n p = 0.0251 1 /fm 3 , n n = 0.258 1 /fm 3 , ξ p = 31.4 fm, ξ n = 84.5 fm, h 1 = 85 MeV fm 5 , h 2 = 323 MeV fm 5 , and h 3 = 219 MeV fm 5 . This means that the magnetic field in Figure 5 is measured in units of c/(2eξ 2 p ) 3 × 10 15 G. These parameter values are reasonable for the core of a neutron star, reflecting realistic nuclear matter equations of state [58] and energy gaps [47,59]. In particular, we have determined h 1 following Chamel and Haensel [31] for the NRAPR equation of state [60], and we have chosen h 2 and h 3 to be of similar magnitude to h 1 . Our analytical results, presented in the previous section, indicate that the different phases seen in Figure 5 are generic, rather than an artefact of our particular parameter choice. For plausible choices of the coupling parameters α and h 2 , we always find type-1.5 superconductivity for a significant range of κ values. This suggests that much of the outer core is a type-1.5 superconductor, rather than type-II as is commonly assumed. However, the other notable feature of Figure 5-the existence of a mixed phase for B c2 < B < H c2 -is present only for sufficiently large values of the parameters h 1 or h 2 , which may explain why a phase of this kind was not seen in the results of Kobyakov [37]. In any event, such a phase could only be found in stars with a mean flux of order H c2 ∼10 15 G.
Since the pioneering work of Baym et al. [1], it has generally been accepted that the outer core of a neutron star is a type-II superconductor, in which fluxtubes are stable but mutually repulsive on all scales. In that case, fluxtubes are expected to form a lattice arrangement that expands over time, eventually leading to a flux-free state (albeit on a very long timescale). Our results suggest that fluxtubes have a preferred separation distance throughout a large part of the core, as a result of entrainment between the protons and neutrons, and will therefore form "hexagonal bundles" that can persist indefinitely.
Furthermore, previous works have often assumed, without providing justification, that the transition between the type-II outer core and type-I inner core occurs where the ratio of the effective London length, λ , to the proton coherence length, ξ p , takes the value 1/ √ 2, e.g., [39,47,61]. In terms of our dimensionless variables, this condition can be expressed as However, our results demonstrate that the transition to the type-I inner core generally occurs where the critical fields H c and H c2 are equal. Although there is no analytical expression for the latter, it is generally quite close to H c2 , and so we have the approximate condition H c H c2 . In terms of our dimensionless variables, this condition becomes An analogous result was obtained by Kobyakov [37]. In a neutron star we typically expect h 1 / to be of order unity and so, even if we neglect the density coupling parameter α, the type-I inner core will be significantly larger than predicted by Equation (52). We note that an alternative method to determine the point of transition to type-I superconductivity is to measure the "surface energy" of the interface between Meissner and non-superconducting regions [37,62,63]. The surface energy in a type-I superconductor is necessarily positive, because otherwise the interface is unstable to the formation of fluxtubes. However, this provides only a necessary condition for stability; the interface could still be unstable to nonlinear perturbations even if the surface energy is positive. Given that we have found the non-superconducting state to be metastable in some cases (i.e., unstable to nonlinear perturbations), we speculate that in those cases the boundary between the type-I and type-1.5 regimes does not occur exactly where the surface energy is zero. Further work will be required to confirm this however.
Although we have only considered the ground state for the condensates in this work, we expect the true microphysical state in the neutron star core to be qualitatively similar, provided that the temperature is well below the condensation temperature. However, the inhomogeneous ground states may be fragile, in the sense that the difference in energy density between the ground state and the fluxtube lattice state, (F (a) − F g (a))/a∼10 −3 , is small. (For the values of n p and ξ p quoted above, this corresponds to a difference in energy density of ∼2 × 10 26 erg /cm 3 .) To understand the time-dependent dynamics of the fluxtubes in detail, including effects of finite temperature, will require a more complete dynamical model than that presented here. In particular, in this work we have considered defects in the proton condensate only, i.e., we have neglected the presence of neutron vortices that will appear as a result of the superfluid's quantized rotation. This is justifiable when characterizing the system's ground state, because the fluxtubes outnumber the vortices by many orders of magnitude, but interactions between both types of defects are crucial for the dynamics of the rotation and magnetic field on larger scales. A fully dynamical description of the neutron star interior must also incorporate the (non-superfluid) electrons, whose scattering by fluxtubes is a dominant source of dissipation in the star's core [39,64,65]. Finding a consistent treatment of the electrons, whose mean free path far exceeds the typical distance between fluxtubes, is beyond the scope of the present work and is left for future study.

Data Availability Statement:
The data used to produce the plots in this paper are available at https://doi.org/10.25405/data.ncl.c.5935723, accessed on 8 April 2022.

Acknowledgments:
The authors would like to thank the Institute for Nuclear Theory at the University of Washington for its kind hospitality and hosting INT Program INT-19-1a during which part of this work was carried out. We also thank Wynn Ho, John Miller, and Hayder Salman for helpful conversations related to this paper and Alexander Haber for providing feedback on our manuscript. Finally, we acknowledge the use of the following software: IPython [66], Matplotlib [67], NumPy [68][69][70], Pandas [71], and SciPy [72,73].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Minimization of the Free Energy
The dimensionless free-energy density, given by Equation (12), is approximated numerically on a regular 2D grid, with intervals δx and δy in the x and y directions. The order parameters ψ p and ψ n are defined on the gridpoints as ψ i,j p and ψ i,j n , where i and j denote the indices in x and y, respectively. The vector field A has only two components, (A x , A y ), which are defined on the corresponding links between the gridpoints, i.e., we have A i+1/2,j x and A i,j+1/2 y . The gauge coupling between ψ p and A is implemented using a standard Peierls substitution, noting that for instance By including the gauge coupling in this way, we exactly preserve the gauge symmetry This leads to a discrete approximation to the total free energy, . We obtain the ground state using a simple gradient-descent method, in which the step size is made as large as possible while maintaining numerical stability. Specifically, we use the iteration scheme We use an initial (dimensionless) resolution of δx, δy 0.5, and iterate until the change in energy drops below a threshold of 10 −7 . We then double the resolution in x and y and repeat the whole process until the value of F dis converges to at least five significant figures.
The reason for working with these real variables will be explained shortly. Recall that F = 0 in the absence of fluxtubes, vortices and magnetic flux, i.e., for f = g = 1 and V = ∇χ = 0. In what follows, we will refer to this as the Meissner solution.
Our goal is to determine the free energy per unit length per Wigner-Seitz cell, F (a), in the asymptotic limit of a widely-spaced fluxtube lattice, a → ∞. In this limit, we know that F converges to the value for a single fluxtube, i.e., F → F ∞ , and we wish to estimate the "interaction energy" F − F ∞ . If its value is positive, the lower transition is of second-order, whereas a negative value implies a first-order transition. To do so, we will first generalize the method introduced by Kramer [54] to the case of a free energy in the completely general form F = F[Ψ] , where Ψ represents the complete set of independent, real variables, and the angled brackets represent a domain integral. We will then apply the results to the particular case given by Equation (A6). We consider an infinite lattice of parallel fluxtubes that are widely separated, in the sense that the size of each Wigner-Seitz cell is large in comparison with both the coherence length and the London length. Any such lattice, characterized by Ψ, is a steady state, in that it is a solution of the Euler-Lagrange equations that are obtained from the functional derivative We assume that each fluxtube is located far inside its Wigner-Seitz cell, which allows us to make two approximations: 1.
Within each cell, the solution can be approximated as a linear perturbation to the single fluxtube solution; 2.
On the boundary of the cell, the solution can be approximated as the Meissner solution plus a superposition of linear, independent perturbations produced by the fluxtubes.
Note that we do not assume that the fluxtube is located exactly at the center of the cell, and we will verify later that the exact location of the fluxtube within the cell does not affect our result for the interaction energy. We will use the notation δ n F[δΨ; Ψ] to represent the n-th variation of F, i.e., the terms up to order (δΨ) n in the Taylor expansion of F[Ψ + δΨ] − F[Ψ]. If Ψ is a solution of the Euler-Lagrange Equation (A7), then δ 1 F[δΨ; Ψ] is an exact derivative. This follows directly from the definition of the functional derivative, which dictates that for some vector field Q[δΨ; Ψ]. Using the divergence theorem, the integral of δ 1 F over any domain can thus be expressed as an integral over the boundary of that domain.
Furthermore, if δΨ is a solution of the linearized Euler-Lagrange equations (i.e., linearized about Ψ), then it can be shown that δ 2 F[δΨ; Ψ] is also an exact derivative. In fact, we have where the vector field Q (2) [δΨ; Ψ] is given precisely by the terms up to order (δΨ) 2 in the Taylor expansion of Q[δΨ; Ψ + 1 2 δΨ]. Once we identify the functional form of δ 1 F, we deduce Q and thence Q (2) . We can then express the integral of δ 2 F over any domain as an integral over the boundary of that domain. The importance of this result in determining the interaction energy will become more obvious shortly. In what follows, we refer to the linear equations for δΨ as the Jacobi equations, by analogy with the equations defining Jacobi fields in Riemannian geometry [74].
Suppose the fluxtubes are aligned with the z-axis, at locations in the xy-plane indexed as x i . Without loss of generality, we will assume that x 0 = 0, and that x −i = −x i . We will also label the Wigner-Seitz cells as C i , such that x i ∈ C i . This allows us to express the interaction energy as where F (all) represents the free-energy density in the presence of the lattice, and F (i) represents the free-energy density in the presence of a single fluxtube at x = x i . To obtain the last line, we have used the fact that, due to the translational symmetry of the lattice, the energy density in cell C i resulting from a single fluxtube at x = 0 is equivalent to the energy density in cell C 0 resulting from a single fluxtube at x = x −i . We now apply assumptions 1 and 2 stated above. Within cell C 0 , we thus approximate Ψ (all) Ψ (0) + δΨ (i =0) , and Ψ (i) Ψ (none) + δΨ (i) for each i = 0, where δΨ (i =0) represents the perturbation produced by all fluxtubes external to C 0 and Ψ (none) is the Meissner solution. Furthermore, we approximate δΨ (i =0) ∑ i =0 δΨ (i) on the boundary of C 0 . We emphasize that δΨ (i) is the linear perturbation to the Meissner solution in the presence of a single fluxtube. In practice, this can be calculated from the Jacobi equations, as we demonstrated in Section 4.2. Under these approximations, we expand F (all) in Equation (A10) and cancel the zeroth-order term F (0) with the i = 0 term in the sum. Expanding the remaining i = 0 contributions and keeping in mind that the Meissner solution has F = 0, we are left with an expression for the interaction energy that only contains second variations of F, and can thus be written in terms of the vector field Q (2) : where ∂C 0 represents the boundary of cell C 0 , and dS is the boundary element on this boundary, with outward normal. The final step is to approximate Ψ (0) Ψ (none) + δΨ (0) on the boundary of C 0 , and retain only terms up to second order in δΨ; this calculation is straightforward once the functional form of Q (2) is known. In this way, we can express the interaction energy as an integral over the boundary of a single Wigner-Seitz cell, and so we do not need to consider the full domain volume to determine the nature of the phase transition at the lower critical field. Note that our assumptions 1 and 2 generally do not apply to the proton order parameter ψ p , because introducing an additional fluxtube causes a nonlinear change in the phase of the proton condensate throughout the domain. Thus, we cannot apply the above method directly to Equation (12). This is the reason for making the change of variables to the set Ψ ≡ ( f , g, V, χ), for which assumptions 1 and 2 do hold. Now taking first variations of the free energy F in the form of Equation (A6) and using integration by parts, we find that The interaction energy can now be calculated following the steps outlined above. For brevity let us just consider the representative term Q = 2 f g δ f ∇g. For this term, we find that Recalling that the Meissner solution has f = g = 1 and ∇χ = V = 0, we find that the contribution from this term to the interaction energy (A11) is In deriving the last equality, we have used the divergence theorem and the fact that the doubly-summed term vanishes, as can be shown using a similar argument to that leading to Equation (A10): Here, ∂R 2 is the boundary of the entire xy-plane, where the integrand is exponentially small.
By applying a similar procedure to the remaining terms in Equation (A12), we eventually find that where δΨ = (δ f , δg, δV, δχ) and We note the similarity between this linear operator L[δΨ] and the Jacobi Equations (36)- (39). In fact, for an arbitrary free-energy functional F[Ψ] it can be shown that the Formula (A14) for the interaction energy still holds, provided that we define i.e., L is defined by the Jacobi equations, implying that L[δΨ (i) ] = 0 at all points except x = x i (the center of the fluxtube), where δΨ (i) is not differentiable. Hence, the integrand in Equation (A14) vanishes at all points inside C 0 , except at x = 0, where it has the form of a delta function. It is for this reason that the exact location of the fluxtube within the Wigner-Seitz cell is immaterial in Equation (A14).
We can now evaluate the integral in Equation (A14) very much as for the simpler case of a single-component superconductor [54]. In the particular case given by Equation (A15), δΨ (i) is given by Equations (40) and (43), after substituting r → |x − x i |. The terms involving δ f and δg can be evaluated by using the following property of the Bessel function K 0 : where δ (2) is the two-dimensional delta function. Using Equations (44) and (45), we then find that all the terms cancel, apart from those involving delta functions, as expected in light of the comments below Equation (A16). The only remaining terms are Again using Equations (44) and (45) to simplify this result, we obtain the second contribution in Equation (46). To evaluate the remaining terms in Equation (A14), it is convenient to use the divergence theorem to rewrite the area integral over C 0 as a contour integral over ∂C 0 in order to avoid the singular behavior of the Bessel functions at r = 0. Remembering that the integrand of Equation (A14) vanishes everywhere but at x = 0, we can shrink the integration contour to a small circle of radius ε centered around the origin. We then have dS = ε dθ e r , and so Using the asymptotic behavior of the Bessel functions, i.e., K 0 (r) ∼ − ln(r) and K 1 (r) ∼ 1/r as r → 0, we observe that in the limit ε → 0 the first term vanishes, while the second one remains finite. More precisely, we find the first contribution in Equation (46). In principle, if we can numerically compute the nonlinear solution for a single fluxtube, from this we can determine the values of V 0 , f 1 and f 2 , and then use Equation (46) to calculate the interaction energy for any lattice of our choosing. In practice, however, it is difficult to obtain both f 1 and f 2 to sufficiently high accuracy to achieve quantitatively reliable results. Moreover, the assumptions made in obtaining this result are only valid in the asymptotic limit of a widely-spaced lattice, so caution is needed when applying this result to a lattice with finite separation between fluxtubes. For these reasons, we would like to have a more robust method for estimating the interaction energy. Haber and Schmitt [8] have suggested a possible approach: they followed essentially the same steps leading to Equation (A14), but chose to leave the result in the form of an integral over the boundary ∂C 0 . They then computed this integral numerically, approximating δΨ (i) using the solution obtained numerically for a single fluxtube. However, their approach has a number of shortcomings:

1.
Rather than computing the interaction energy for a lattice, they considered a pair of fluxtubes. However, this is generally not a steady state, i.e., it is not a solution of the Euler-Lagrange Equations (A7). This violates a basic assumption underlying the derivation.

2.
They chose to include some, but not all, of the higher-order terms in their calculation, leading to a result that lacks certain symmetries expected on physical grounds. By contrast, in deriving Equation (A14), we have consistently neglected all terms of higher order than (δΨ) 2 , and the result is antisymmetric between δΨ (0) and δΨ (i) .

3.
Their formula (C8) for the interaction energy depends on the location of the Wigner-Seitz cell boundary, relative to the fluxtube lattice. As we emphasized above, the location of the fluxtube within its cell is immaterial, and therefore should not change the result.
As an alternative approach, we suggest making use of the exact result dF d ln a = C 0 which we derive in Appendix C. Inside the integral, we can approximate the full solution by superposing the profiles of single fluxtubes with the Meissner solution: This formula is consistent with the rigorous result (46) in the asymptotic limit a → ∞, and because the integrand is spatially periodic by construction, it has none of the shortcomings described above. The formula will be accurate as long as the approximations (A22) hold, which in practice still requires that a 1. We have used this formula to independently verify some of the results from our 2D numerical model.
where the overbar represents the spatial average, which is equivalent to the area integral over the rescaled rectangular domain. Because the quantity a now appears only as a coefficient in the free energy, we can directly compute the derivative of F (a), which leads to the formula (A21).
In the rescaled units, the non-superconducting solution has the form |ψ p | = 0, |ψ n | = √ 1 + α/R 2 , A = πNe z × x. At the critical point a = 2π/H c2 , this solution becomes unstable to perturbations δψ p that lie in the kernel of the linear operator These perturbations have the form δψ p = e iπNy(x+iy) φ(x + iy), for some function φ, which must be chosen to match the quasi-periodic boundary condition (A25). The general solution can be expressed in terms of Jacobi theta functions: where the fluxtube locations, z j = x j + iy j , must satisfy ∑ j z j = (m + 1 2 N)Γ + (n + 1 2 N)i/Γ, for some m, n ∈ Z.