Next Article in Journal
Ultraviolet Finiteness or Asymptotic Safety in Higher Derivative Gravitational Theories
Next Article in Special Issue
Oscillating Magnetized Color Superconducting Quark Stars
Previous Article in Journal
Muon to Positron Conversion
Previous Article in Special Issue
1S0 Pairing Gaps, Chemical Potentials and Entrainment Matrix in Superfluid Neutron-Star Cores for the Brussels–Montreal Functionals
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Superconducting Phases in Neutron Star Cores

1
School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne NE1 7RU, UK
2
Institute of Space Sciences (ICE-CSIC), Campus UAB, 08193 Barcelona, Spain
3
Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain
*
Author to whom correspondence should be addressed.
Submission received: 14 February 2022 / Revised: 5 April 2022 / Accepted: 6 April 2022 / Published: 8 April 2022
(This article belongs to the Special Issue Superfluidity and Superconductivity in Neutron Stars)

Abstract

:
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 field strength. Combining analytical and numerical techniques, we find 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 flux is distributed inhomogeneously, with bundles of magnetic fluxtubes separated by flux-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 fluxtubes can coexist with non-superconducting regions, but only in a small part of the parameter space.

1. 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 condensates. We adopt the term “type-1.5” superconductivity to maintain consistency with earlier literature on the topic. We note that there are other mechanisms that might produce bundles of fluxtubes inside neutron stars, such as a macroscopic instability of the fluxtube lattice [14] or interactions with quantized vortices in the neutron superfluid [15]. However, these scenarios do not correspond to type-1.5 superconductivity in the sense used here.
The nature of superconductivity in the star’s core may have observational consequences. For example, there is observational evidence for long-period precession in some neutron stars [16,17,18,19], which places constraints on dissipation within the core. If the superfluid neutron vortices in the outer core are tightly pinned to a regular array of fluxtubes then precession becomes essentially impossible [20,21]. This contradiction could be resolved if the entire core were in a type-I state [21,22,23], or if instabilities resulting from precession are able to unpin the vortices [24]. Pinning between vortices and fluxtubes could also significantly affect the rise time, amplitude, and relaxation rate of rotational glitches in pulsars [25,26,27], as well as the long-term evolution of their rotation and magnetic field [28,29].
The goal of this paper is to extend previous studies to include all relevant couplings between the neutron and proton condensates, and to determine the precise conditions under which different types of superconductivity arise. We use a Ginzburg–Landau model for the condensates, which is the simplest phenomenological description that captures the quantization of vorticity and magnetic flux. Specifically, we work with the most general Ginzburg–Landau functional that permits a consistent treatment of entrainment, while also correctly satisfying Galilean invariance on small scales [30,31]. Using analytical and numerical techniques, we then find the ground state for the condensates by minimizing the free energy of the two-component system in the presence of an imposed magnetic field, and thus construct superconducting phase diagrams. Although our aim here is only to identify the ground state, and hence to determine the type of superconductivity, we note that this local ground state is likely to be a good approximation to the situation in real neutron stars on small scales. This is because young neutron stars cool very efficiently via neutrino emission [32,33], and after ∼ 10 4 yr their internal temperatures lie far below the critical temperatures for superconductivity and superfluidity [34,35,36].
The paper is organized as follows. In Section 2 we present our Ginzburg–Landau model, including a review of entrainment and Galilean invariance. In Section 3 we solve the Ginzburg–Landau equations numerically, and hence construct superconducting phase diagrams. Section 4 presents analytical results concerning the various phase transitions found numerically. Finally, we discuss the implications for neutron stars in Section 5. Further details of some of the calculations are provided in Appendix A, Appendix B and Appendix C.

2. The Ginzburg–Landau Formalism

Our goal is to formulate a simple, phenomenological model of the neutron and proton condensates that includes (a) their mutual entrainment, and (b) the coupling to the magnetic field. The simplest such model is the Ginzburg–Landau model, in which the free energy of the condensates is expressed in terms of complex scalar order parameters for the proton and neutron condensates, ψ p and ψ n , and the magnetic vector potential, A . We acknowledge that such a model can only be rigorously justified close to the condensation temperatures for the protons and neutrons, whereas mature neutron star cores are generally well below both of those temperatures. This limits the rigorous applicability of the Ginzburg–Landau model to the part of the core in which the condensation temperatures are similar [37], and to the narrow time window in which the condensation occurs. Furthermore, the Ginzburg–Landau model neglects spin–orbit interactions, and therefore does not include all of the microphysics of a neutron star. In particular, the neutrons are believed to form an anisotropic 3 P 2 superfluid throughout much of the core, which cannot be described by a single scalar order parameter, e.g., [38]. Nevertheless, the relative simplicity of the Ginzburg–Landau formalism makes it an appealing phenomenological model, and hence it has been widely used to model the superfluid components of neutron stars, e.g., [7,8,39,40,41,42,43]. Moreover, phenomenological Ginzburg–Landau models have been shown to reproduce many properties of laboratory 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 = 2 m 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.

2.1. 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.,
arg ψ p = 2 m p V p , arg ψ n = 2 m n V n ,
where V p and V n are the superfluid velocities. In the presence of entrainment, the velocity-dependent terms in the free-energy density, F vel , of the condensates must take the form
F vel = 1 2 ρ p | V p | 2 + 1 2 ρ n | V n | 2 1 2 ρ pn | V p V n | 2 ,
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
F ent = 1 2 ( h 1 + h 2 ) m n m p 1 / 2 ψ n ψ p + m p m n 1 / 2 ψ p ψ n 2 + 1 2 ( h 1 h 2 ) m n m p 1 / 2 ψ n ψ p m p m n 1 / 2 ψ p ψ n 2 + 1 4 h 3 ( ψ p ψ p ) 2 + 1 4 h 4 ( ψ n ψ n ) 2 ,
which includes four real, independent parameters h 1 , …, h 4 . In terms of the superfluid densities and velocities, we can write this as
F ent = h 1 ρ n 4 m p 2 ρ p 1 / 2 2 + ρ p 4 m n 2 ρ n 1 / 2 2 + ρ p ρ n 2 | V p V n | 2 + h 2 8 m p m n ρ p · ρ n + h 3 16 m p 2 | ρ p | 2 + h 4 16 m n 2 | ρ n | 2 .
So by comparison with Equation (2) the entrainment coefficient is
ρ pn = 2 2 h 1 ρ p ρ n .
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.

2.2. 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
F vel = 1 2 ρ pp | V p | 2 + 1 2 ρ nn | V n | 2 + ρ pn V p · V n ,
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 ρ x y 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 non-relativistic 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.

3. Superconducting Ground States

3.1. 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
F [ ψ p , ψ n , A ] = g pp 2 | ψ p | 2 n p 2 2 + g nn 2 | ψ n | 2 n n 2 2 + g pn | ψ p | 2 n p 2 | ψ n | 2 n n 2 + 1 8 π | × A | 2 + 2 4 m u 2 i e c A ψ p 2 + 2 4 m u ψ n 2 + h 1 2 i e c A ( ψ n ψ p ) 2 + h 2 h 1 2 ( | ψ p | 2 ) · ( | ψ n | 2 ) + h 3 4 ( | ψ p | 2 ) 2 + ( | ψ n | 2 ) 2 ,
where we have assumed that the proton Cooper pairs have charge 2 e . 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 position-independent densities | ψ p | 2 = n p / 2 and | ψ n | 2 = n n / 2 . Therefore the parameters n p and nn 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 pn 2 > 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 pn 2 < 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 2 m u g pp n p and ξ n 2 m u g nn n n ,
and the “bare” London length is defined as
λ m u c 2 4 π e 2 n p .
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 / ( 2 e ξ p ) , and the coupling coefficients h i in units of g pp ξ p 2 . 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 p 2 / 4 ,
F [ ψ p , ψ n , A ] = i A ψ p 2 + 1 ϵ ψ n 2 + κ 2 | × A | 2 + 1 2 ( 1 | ψ p | 2 ) 2 + R 2 2 ϵ ( 1 | ψ n | 2 ) 2 + α ϵ ( 1 | ψ p | 2 ) ( 1 | ψ n | 2 ) + h 1 ϵ i A ( ψ n ψ p ) 2 + ( h 2 h 1 ) 2 ϵ ( | ψ p | 2 ) · ( | ψ n | 2 ) + h 3 4 ( | ψ p | 2 ) 2 + 1 ϵ 2 ( | ψ n | 2 ) 2 ,
where we have defined the following parameters:
κ λ ξ p , R ξ p ξ n , ϵ n p n n , α g pn g pp .
Note that κ is equivalent to our dimensionless “bare” London length.

3.2. 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,
F = F ,
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,
G = F 2 κ 2 H · × A = F 2 κ 2 H · B .
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 c 1 < | H | < H c 2 the magnetic flux organizes into a hexagonal lattice of discrete fluxtubes. The transitions at the lower critical field, H c 1 , and the upper critical field, H c 2 , are both second-order, because fluxtubes appear with infinite separation at | H | = H c 1 , and the superconductor density becomes vanishingly small at | H | = H c 2 . In our dimensionless units, H c 2 = 1 and H c 1 = 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 c 2 , 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 c 1 < H c 1 and | H | = H c 2 > H c 2 , 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.

3.3. 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:
κ 2 × ( × A ) = Im ψ p ( i A ) ψ p + h 1 ϵ ψ n ψ p ( i A ) ( ψ n ψ p ) , 2 ψ n = R 2 ( | ψ n | 2 1 ) ψ n + α ( | ψ p | 2 1 ) ψ n
h 1 ψ p ( + i A ) 2 ( ψ p ψ n ) ψ n 2 h 2 h 1 2 | ψ p | 2 + h 3 2 ϵ | ψ n | 2 , ( i A ) 2 ψ p = ( | ψ p | 2 1 ) ψ p + α ϵ ( | ψ n | 2 1 ) ψ p
h 1 ϵ ψ n ( i A ) 2 ( ψ n ψ p ) ψ p 2 h 2 h 1 2 ϵ | ψ n | 2 + h 3 2 | ψ p | 2 .
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 x y -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
A ( x + L ) = A ( x ) + N π L x L y e z × L ,
ψ p ( x + L ) = ψ p ( x ) exp i N π L x L y e z × L · x ,
ψ n ( x + L ) = ψ n ( x ) ,
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:
F 1 N x = 0 L x y = 0 L y F d x d y .
As an example, in Figure 2 we plot F as a function of the area per magnetic flux quantum,
a L x L y N = 2 π B ¯ ,
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).
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:
a ¯ = ϕ a 1 + ( 1 ϕ ) a 2 , and F ¯ = ϕ F 1 + ( 1 ϕ ) F 2 .
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
G ¯ = G / a = F g ( a ) a 4 π κ 2 a | H | ,
and the minimum of G ¯ over all a can be interpreted graphically from the plots in Figure 2. Since F g ( a ) 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. The two ranges of a for which the function F g ( a ) is linear give rise to two critical values of | H | (i.e.,  H c 1 and H c 2 ) at which the ground state changes discontinuously. At these values there are first-order transitions between a hexagonal fluxtube lattice with finite mean field and either the Meissner state (at H c 1 ) or the non-superconducting state (at H c 2 ).
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.

3.4. 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 c 1 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 c 1 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 c 2 , as the proton condensate density vanishes discontinuously; in the experiment with B ¯ , there is a critical value B ¯ = B c 2 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.

4. 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 c 1 and H c 2 , 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.

4.1. 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
F = a 2 1 α 2 ϵ R 2 + ( 2 π κ ) 2 a .
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
G ¯ = 1 2 1 α 2 ϵ R 2 κ 2 | H | 2 .
Thus, assuming that we have a type-I superconductor, there is a first-order transition between the Meissner and the non-superconducting state at
| H | = H c 1 2 κ 1 α 2 ϵ R 2 .
This defines the critical field, H c , in agreement with previous studies [8,37].

4.2. The Lower Critical Field, H c 1 vs.  H c 1

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 c 1 = 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 c 1 = F / ( 4 π κ 2 ) , as in the case of a single-component superconductor, where
F = lim a F ( a ) .
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 ) , A = A θ ( r ) e θ .
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
f | ψ p | , g | ψ n | , χ arg ψ n , V ( arg ψ p ) A .
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
0 = f ( f 2 1 ) + α ϵ f ( g 2 1 ) 2 f + f | V | 2
+ h 1 ϵ [ f | g | 2 + f g 2 | V χ | 2 · ( g 2 f ) ] h 2 2 ϵ f 2 ( g 2 ) h 3 2 f 2 ( f 2 ) ,
0 = R 2 ϵ g ( g 2 1 ) + α ϵ g ( f 2 1 ) 1 ϵ 2 g + 1 ϵ g | χ | 2
+ h 1 ϵ [ g | f | 2 + g f 2 | V χ | 2 · ( f 2 g ) ] h 2 2 ϵ g 2 ( f 2 ) h 3 2 ϵ 2 g 2 ( g 2 ) ,
0 = · [ g 2 χ + h 1 f 2 g 2 ( χ V ) ] ,
0 = f 2 V + h 1 ϵ f 2 g 2 ( V χ ) + κ 2 × ( × V ) .
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:
1 + h 1 ϵ + h 3 2 δ f + h 2 ϵ 2 δ g = 2 δ f + 2 α ϵ δ g ,
1 + h 1 + h 3 ϵ 2 δ g + h 2 2 δ f = 2 R 2 δ g + 2 α δ f ,
κ 2 × ( × δ V ) = h 1 ϵ δ χ 1 + h 1 ϵ δ V ,
1 + h 1 2 δ χ = h 1 · δ V ,
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
δ V = V 0 K 1 r λ e θ ,
for some constant coefficient V 0 , where K 1 is a modified Bessel function and λ is the effective London length,
λ = κ 1 + h 1 ϵ 1 / 2 .
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-coherence-length structure, i.e.,
δ f = f 1 K 0 2 r / ξ 1 + f 2 K 0 2 r / ξ 2 ,
δ g = g 1 K 0 2 r / ξ 1 + g 2 K 0 2 r / ξ 2 ,
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
1 + h 1 ϵ + h 3 f i ξ i 2 + h 2 ϵ g i ξ i 2 = f i + α ϵ g i ,
1 + h 1 + h 3 ϵ g i ξ i 2 + h 2 f i ξ i 2 = R 2 g i + α f i .
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
F F 2 π κ 2 V 0 2 i 0 K 0 | x i | / λ 2 π j = 1 , 2 f j 2 + 2 α ϵ f j g j + R 2 ϵ g j 2 ξ j 2 i 0 K 0 2 | x i | / ξ j ,
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 ,
F F κ 2 V 0 2 λ d e d / λ f + 2 + 2 α ϵ f + g + + R 2 ϵ g + 2 ξ + 2 ξ + 2 d e 2 d / ξ +
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 c 1 , 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 c 1 , 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 c 1 = H c 1 (and B c 1 = 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.

4.3. The Upper Critical Field, H c 2 vs.  H c 2

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 c 2 , 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 B ¯ e 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
1 + h 1 ϵ 1 + α R 2 1 2 i B ¯ e z × x 2 δ ψ p = 1 α 2 ϵ R 2 δ ψ p .
As known from single-component systems, e.g., [52], bounded solutions of this equation first appear when
B ¯ = H c 2 1 α 2 ϵ R 2 1 + h 1 ϵ 1 + α R 2 .
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 c 2 . However, if the function F ( a ) is not convex at this point then the upper transition will occur not at | H | = H c 2 , but at a higher value | H | = H c 2 , 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 c 2 , 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
2 H c 2 + h 3 1 α 2 ϵ R 2 2 1 κ 2 H c 2 3 π β ϵ R 2 = i exp ( 1 2 H c 2 | x i | 2 ) ( h 2 h 1 ) 1 2 H c 2 | x i | 2 + h 1 + α 1 H c 2 2 1 R 2 + α + h 3 ϵ R 2 1 2 H c 2 | x i | 2 + 1 H c 2 ,
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,
β | ψ p | 4 ¯ ( | ψ p | 2 ¯ ) 2 = i exp ( 1 2 H c 2 | x i | 2 ) .
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 c 2 and H c 2 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 c 2 , H c 2 and B c 2 all coincide. In the case shown in Figure 5, this yields the value κ 1.58 , in agreement with our numerical results.

5. 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 density-gradient 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 c 1 ), 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 c 2 < B ¯ < H c 2 ) 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 / f m 3 , n n = 0.258 1 / f m 3 , ξ p = 31.4 f m , ξ n = 84.5 f m , h 1 = 85 MeV f m 5 , h 2 = 323 MeV f m 5 , and h 3 = 219 MeV f m 5 . This means that the magnetic field in Figure 5 is measured in units of c / ( 2 e ξ p 2 ) 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 c 2 < B ¯ < H c 2 —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 c 2 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
κ = 1 2 + h 1 2 ϵ .
However, our results demonstrate that the transition to the type-I inner core generally occurs where the critical fields H c and H c 2 are equal. Although there is no analytical expression for the latter, it is generally quite close to H c 2 , and so we have the approximate condition H c H c 2 . In terms of our dimensionless variables, this condition becomes
κ 1 + h 1 ϵ 1 + α R 2 2 2 α 2 ϵ R 2 1 / 2 .
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 / c m 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.

Author Contributions

Formal analysis, T.S.W. and V.G.; Investigation, T.S.W. and V.G.; Methodology, T.S.W. and V.G.; Visualization, T.S.W. and V.G.; Writing—original draft, T.S.W. and V.G.; Writing—review and editing, T.S.W. and V.G. All authors have read and agreed to the published version of the manuscript.

Funding

Part of TSW’s time was funded by EPSRC Grant EP/R024952/1 and by STFC grant ST/W001020/1. VG acknowledges partial support from a McGill Space Institute postdoctoral fellowship and the Trottier Chair in Astrophysics and Cosmology as well as the H2020 ERC Consolidator Grant "MAGNESIA" under grant agreement No. 817661 (PI: Rea) and Spanish National Grant PGC2018-095512-BI00. This work was also partially supported by the program Unidad de Excelencia María de Maeztu CEX2020-001058-M, and by the PHAROS COST Action (No. CA16214).

Data Availability Statement

The data used to produce the plots in this paper are available at https://0-doi-org.brum.beds.ac.uk/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 ψ p i , j and ψ n i , j , 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 x i + 1 / 2 , j and A y i , j + 1 / 2 . The gauge coupling between ψ p and A is implemented using a standard Peierls substitution, noting that for instance
x i A x ψ p = x exp ( i A x d x ) ψ p x i A x ψ p i + 1 / 2 , j 1 δ x exp ( i A x i + 1 / 2 , j δ x ) ψ p i + 1 , j ψ p i , j .
By including the gauge coupling in this way, we exactly preserve the gauge symmetry
ψ p i , j exp ( i ϕ i , j ) ψ p i , j , A x i + 1 / 2 , j A x i + 1 / 2 , j + ϕ i + 1 , j ϕ i , j δ x , A y i , j + 1 / 2 A y i , j + 1 / 2 + ϕ i , j + 1 ϕ i , j δ y .
This leads to a discrete approximation to the total free energy, F dis [ ψ p i , j , ψ n i , j , A x i + 1 / 2 , j , A y i , j + 1 / 2 ] . 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
ψ p i , j ψ p i , j N / 4 δ y / δ x + δ x / δ y F dis / ψ p , i , j 1 + h 1 ϵ | ψ n i , j | 2 + h 3 | ψ p i , j | 2 ,
ψ n i , j ψ n i , j N / 4 δ y / δ x + δ x / δ y F dis / ψ n , i , j 1 + h 1 | ψ p i , j | 2 + h 3 ϵ | ψ n i , j | 2 ,
A x i + 1 / 2 , j A x i + 1 / 2 , j N / 4 δ y / δ x + δ x / δ y F dis / A x i + 1 / 2 , j 2 κ 2 ,
A y i , j + 1 / 2 A x i , j + 1 / 2 N / 4 δ y / δ x + δ x / δ y F dis / A y i , j + 1 / 2 2 κ 2 .
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.

Appendix B. Long-Range Interaction between Fluxtubes

The dimensionless free-energy density (12), when written in terms of the real variables f, g, V and χ defined in Section 4.2, takes the form
F [ f , g , V , χ ] = 1 2 ( f 2 1 ) 2 + R 2 2 ϵ ( g 2 1 ) 2 + α ϵ ( f 2 1 ) ( g 2 1 ) + f 2 + f 2 V 2 + 1 ϵ g 2 + 1 ϵ g 2 χ 2 + κ 2 | × V | 2 + h 1 ϵ f 2 g 2 + g 2 f 2 + f 2 g 2 | V χ | 2 + h 2 2 ϵ f 2 · g 2 + h 3 4 1 ϵ 2 g 2 2 + f 2 2 .
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
δ δ Ψ F = 0 .
We assume that each fluxtube is located far inside its Wigner–Seitz cell, which allows us to make two approximations:
  • Within each cell, the solution can be approximated as a linear perturbation to the single fluxtube solution;
  • 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
δ 1 F = δ Ψ · δ δ Ψ F EL equations + · Q ,
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
δ 2 F = δ Ψ · δ δ Ψ F EL equations + 1 2 δ Ψ · δ 1 δ δ Ψ F linearized EL equations + · Q ( 2 ) ,
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 x y -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
F F = C 0 F ( all ) d x d y R 2 F ( 0 ) d x d y = C 0 F ( all ) i F ( i ) d x d y ,
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 ) :
F F C 0 δ 2 F [ δ Ψ ( i 0 ) ; Ψ ( 0 ) ] i 0 δ 2 F [ δ Ψ ( i ) ; Ψ ( none ) ] d x d y = C 0 Q ( 2 ) [ δ Ψ ( i 0 ) ; Ψ ( 0 ) ] i 0 Q ( 2 ) [ δ Ψ ( i ) ; Ψ ( none ) ] · d S C 0 Q ( 2 ) [ i 0 δ Ψ ( i ) ; Ψ ( 0 ) ] i 0 Q ( 2 ) [ δ Ψ ( i ) ; Ψ ( none ) ] · d S ,
where C 0 represents the boundary of cell C 0 , and d S 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
Q = 2 δ f f + 2 ϵ δ g g + 2 ϵ g 2 δ χ χ + 2 κ 2 δ V × ( × V ) + 2 h 1 ϵ g 2 δ f f + f 2 δ g g + f 2 g 2 δ χ ( χ V ) + 2 h 2 ϵ f g δ f g + δ g f + 2 h 3 f 2 δ f f + 1 ϵ 2 g 2 δ g g .
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
Q ( 2 ) = Q + g ( δ f ) 2 g + f δ f δ g g + f g δ f δ g .
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
C 0 i 0 2 δ f ( i ) δ g ( 0 ) + j 0 δ f ( i ) δ g ( j ) δ f ( i ) δ g ( i ) · d S = C 0 i , j i δ f ( i ) δ g ( j ) + i δ f ( i ) δ g ( 0 ) δ f ( 0 ) δ g ( i ) · d S = i C 0 · δ f ( i ) δ g ( 0 ) δ f ( 0 ) δ g ( i ) d x d y .
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):
i j i C 0 δ f ( i ) δ g ( j ) · d S = i j 0 C i δ f ( 0 ) δ g ( j ) · d S = j 0 R 2 δ f ( 0 ) δ g ( j ) · d S = 0 .
Here, R 2 is the boundary of the entire x y -plane, where the integrand is exponentially small.
By applying a similar procedure to the remaining terms in Equation (A12), we eventually find that
F F i C 0 δ Ψ ( i ) · L [ δ Ψ ( 0 ) ] δ Ψ ( 0 ) · L [ δ Ψ ( i ) ] d x d y ,
where δ Ψ = ( δ f , δ g , δ V , δ χ ) and
L [ δ Ψ ] = 1 + h 1 ϵ + h 3 2 δ f + h 2 ϵ 2 δ g 1 ϵ 1 + h 1 + h 3 ϵ 2 δ g + h 2 ϵ 2 δ f κ 2 × ( × δ V ) + h 1 ϵ δ χ 1 ϵ 2 δ χ + h 1 ϵ · ( δ χ δ V ) .
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
L [ δ Ψ ] 1 2 δ 1 δ δ Ψ F [ δ Ψ ; Ψ ( none ) ] ,
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 :
2 K 0 ( k | x x i , 0 | ) = k 2 K 0 ( k | x x i , 0 | ) 2 π δ ( 2 ) ( x x i , 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
F F δ f , δ g 2 π i 0 K 0 2 | x i | ξ j j = 1 , 2 1 + h 1 ϵ + h 3 f j 2 + 2 h 2 ϵ f j g j + 1 ϵ 1 + h 1 + h 3 ϵ g j 2 .
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 d S = ε d θ e r , and so
F F δ V κ 2 i 0 C 0 δ V ( i ) × ( × δ V ( 0 ) ) δ V ( 0 ) × ( × δ V ( i ) ) · d S = κ 2 V 0 i 0 0 2 π ε δ V ( i ) | x = 0 K 0 ε λ 1 λ + K 1 ε λ ( × δ V ( i ) ) z | x = 0 d θ .
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
F F δ V κ 2 V 0 i 0 0 2 π λ ( × δ V ( i ) ) z | x = 0 d θ = 2 π κ 2 V 0 2 i 0 K 0 | x i | λ ,
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:
  • 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.
  • 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 ) .
  • 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
d F d ln a = C 0 [ 1 2 ( f 2 1 ) 2 + R 2 2 ϵ ( g 2 1 ) 2 + α ϵ ( f 2 1 ) ( g 2 1 ) κ 2 B z 2 ] d x d y ,
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:
f 1 + i ( f ( i ) 1 ) , g 1 + i ( g ( i ) 1 ) , B z i B z ( i ) .
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.

Appendix C. Weakly Nonlinear Lattice Solution

We consider a rectangular domain that contains an integer number of fluxtubes, N, with area a N . We have shown in Section 4.3 that the non-superconducting state is linearly unstable for a > 2 π / H c 2 , where H c 2 is given by Equation (49). Our goal is to compute weakly nonlinear solutions for a = 2 π / H c 2 + δ a , where δ a 1 . To do so we will roughly follow the same procedure as Abrikosov [56], except that by working with a finite domain and quasi-periodic boundary conditions we avoid having to manipulate products of infinite series.
It is convenient at this point to redefine our length scale, so that the area of the domain is normalized to unity, and remains fixed as the parameter a is varied. At the same time, we will also rescale A so that the boundary conditions have no dependence on a. Under the rescaling x ( a N ) 1 / 2 x and A ( a N ) 1 / 2 A , the free energy (12) becomes
F [ ψ p , ψ n , A ] = 1 2 ( 1 | ψ p | 2 ) 2 + R 2 2 ϵ ( 1 | ψ n | 2 ) 2 + α ϵ ( 1 | ψ p | 2 ) ( 1 | ψ n | 2 ) + 1 a N i A ψ p 2 + 1 ϵ a N ψ n 2 + κ 2 ( a N ) 2 | × A | 2 + h 1 ϵ a N i A ( ψ n ψ p ) 2 + ( h 2 h 1 ) 2 ϵ a N ( | ψ p | 2 ) · ( | ψ n | 2 ) + h 3 4 a N ( | ψ p | 2 ) 2 + 1 ϵ 2 ( | ψ n | 2 ) 2 .
Our domain now has the dimensions Γ × ( 1 / Γ ) , say, and we have the (quasi)periodic boundary conditions:
A ( x + L ) = A ( x ) + π N e z × L ,
ψ p ( x + L ) = ψ p ( x ) exp i π N e z × L · x ,
ψ n ( x + L ) = ψ n ( x ) ,
where L represents either of the translation symmetries ( Γ , 0 ) or ( 0 , 1 / Γ ) . The free energy per magnetic flux quantum (in the unscaled units) is
F ( a ) = a F ¯ ,
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 = π N e z × x . At the critical point a = 2 π / H c 2 , this solution becomes unstable to perturbations δ ψ p that lie in the kernel of the linear operator
L ( i π N e z × x ) 2 + 2 π N .
These perturbations have the form δ ψ p = e i π N y ( x + i y ) ϕ ( x + i y ) , 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:
ϕ ( z ) = j = 1 N exp ( 2 π i y j z ) ϑ 1 π Γ ( z z j ) | i Γ 2 ,
where the fluxtube locations, z j = x j + i y j , must satisfy j z j = ( m + 1 2 N ) Γ + ( n + 1 2 N ) i / Γ , for some m , n Z .
Just beyond the critical point, with a = 2 π / H c 2 + δ a , we anticipate that the solutions have regular asymptotic expansions of the form
A = A ( 1 ) + ( δ a ) A ( 2 ) + ,
ψ n = ψ n ( 1 ) + ( δ a ) ψ n ( 2 ) + ,
ψ p = ( δ a ) 1 / 2 ψ p ( 1 ) + ( δ a ) 3 / 2 ψ p ( 2 ) + .
Substituting this ansatz into the rescaled Euler–Lagrange equations, the leading-order contributions recover the non-superconducting state for A ( 1 ) and ψ n ( 1 ) , as well as the linear equation L ψ p ( 1 ) = 0 . Without loss of generality, we will take A ( 1 ) = π N e z × x and ψ n ( 1 ) = 1 + α / R 2 . For the moment we do not need to choose the particular form of ψ p ( 1 ) . However, in what follows we will make use of its quasi-periodicity, and of the (gauge-invariant) identities
i A ( 1 ) ψ p ( 1 ) = i e z × i A ( 1 ) ψ p ( 1 ) ,
1 2 π 2 ln | ψ p ( 1 ) | = N + j δ ( 2 ) ( x x j ) ,
where x j are the fluxtube locations, and δ ( 2 ) is the two-dimensional delta function.
Now proceeding to next order in the Euler–Lagrange equations, we eventually obtain the following:
1 α 2 ϵ R 2 1 H c 2 2 κ 2 π N × ( × A ( 2 ) ) = e z × | ψ p ( 1 ) | 2 , 1 | ψ n ( 1 ) | 2 + h 3 ϵ 2 4 π N R 2 H c 2 Re ψ n ( 1 ) ψ n ( 2 )
= h 2 h 1 2 2 2 π N h 1 + α H c 2 | ψ p ( 1 ) | 2 ,
2 Im ψ n ( 1 ) ψ n ( 2 ) = 0 ,
as well as a lengthy equation for ψ p ( 2 ) :
1 α 2 ϵ R 2 1 H c 2 L ψ p ( 2 ) = N 2 π H c 2 | ψ p ( 1 ) | 2 1 + α 2 ϵ R 2 ψ p ( 1 ) h 3 2 ψ p ( 1 ) 2 ( | ψ p ( 1 ) | 2 ) ψ p ( 1 ) h 2 h 1 ϵ 2 α ϵ 4 π N H c 2 Re ψ n ( 1 ) ψ n ( 2 ) + 1 α 2 ϵ R 2 1 H c 2 ( i A ( 1 ) ) · ( i A ( 2 ) ψ p ( 1 ) ) + i A ( 2 ) · ( i A ( 1 ) ) ψ p ( 1 ) h 1 ϵ ψ n ( 1 ) ψ n ( 2 ) ( i A ( 1 ) ) 2 ψ p ( 1 ) + ( i A ( 1 ) ) 2 ( ψ p ( 1 ) ψ n ( 1 ) ψ n ( 2 ) ) .
Equation (A35) can be integrated once to obtain B z ( 2 ) . We note that the boundary condition (A24) implies that A ( n ) is spatially periodic for all n > 1 , and so B z ( n ) ¯ = 0 . Hence we find that
B z ( 2 ) = 1 α 2 ϵ R 2 π N H c 2 2 κ 2 | ψ p ( 1 ) | 2 ¯ | ψ p ( 1 ) | 2 .
Now, in order for Equation (A38) to have regular solutions, its right-hand side must be orthogonal to ψ p ( 1 ) . To prove this, we note that L is self-adjoint with respect to the inner product
ψ , ϕ x = 0 Γ y = 0 1 / Γ ψ ϕ d y d x ,
provided that both arguments satisfy the quasi-periodic boundary condition (A25). All of the terms ψ p ( n ) satisfy these boundary conditions, and so
ψ p ( 1 ) , L ψ p ( 2 ) = L ψ p ( 1 ) , ψ p ( 2 ) = 0 .
Hence, taking the inner product of ψ p ( 1 ) with Equation (A38), we obtain the compatibility condition
0 = N 2 π H c 2 | ψ p ( 1 ) | 4 ¯ 1 α 2 ϵ R 2 | ψ p ( 1 ) | 2 ¯ + h 3 2 | ( | ψ p ( 1 ) | 2 ) | 2 ¯ + h 2 h 1 ϵ ( | ψ p ( 1 ) | 2 ) · Re { ψ n ( 1 ) ψ n ( 2 ) } ¯ + α ϵ 4 π N H c 2 | ψ p ( 1 ) | 2 Re { ψ n ( 1 ) ψ n ( 2 ) } ¯ + 1 α 2 ϵ R 2 1 H c 2 | ψ p ( 1 ) | 2 B z ( 2 ) ¯ + 4 π N h 1 ϵ | ψ p ( 1 ) | 2 Re { ψ n ( 1 ) ψ n ( 2 ) } ¯ .
Using Equation (A34), it can be shown that
| ( | ψ p ( 1 ) | 2 ) | 2 ¯ = 2 π N | ψ p ( 1 ) | 4 ¯ ,
and combining this with Equations (A36) and (A39), we can write the compatibility condition (A41) in a more symmetric form:
1 α 2 ϵ R 2 1 | ψ p ( 1 ) | 2 ¯ = 2 π H c 2 + π h 3 β + 1 α 2 ϵ R 2 2 π H c 2 3 κ 2 1 β R 2 ϵ 1 2 N 1 R 2 + α + h 3 ϵ R 2 γ 2 ¯ + 2 π H c 2 γ 2 ¯ ,
where we have defined
β | ψ p ( 1 ) | 4 ¯ | ψ p ( 1 ) | 2 ¯ 2 and γ ( x ) 2 Re ψ n ( 1 ) ψ n ( 2 ) | ψ p ( 1 ) | 2 ¯ .
If we now expand the free-energy density (A23) up to O ( δ a 2 ) , and use the identities derived above, we eventually find
F = a F ¯ = a 2 1 α 2 ϵ R 2 1 H c 2 ( δ a ) 2 2 π | ψ p ( 1 ) | 2 ¯ + ( 2 π κ ) 2 a + O ( δ a 3 ) .
The transition to the non-superconducting state is second order if (and only if) F ( a ) is convex in a neighborhood of the critical point a = 2 π / H c 2 , which requires that F < 0 and F > 0 . Using Equation (A45), and the compatibility condition (A43), these criteria become
( H c 2 κ ) 2 > 1 2 1 α 2 ϵ R 2 ,
2 H c 2 + h 3 1 α 2 ϵ R 2 2 1 H c 2 3 κ 2 π β ϵ R 2 > 1 2 N 1 R 2 + α + h 3 ϵ R 2 γ 2 ¯ + 2 π H c 2 γ 2 ¯ .
In a single-component Ginzburg–Landau superconductor these criteria both reduce to κ > 1 / 2 , so the upper transition in a type-II superconductor is always second-order [56]. Moreover, minimizing the free energy is equivalent to minimizing β , and hence the hexagonal lattice is energetically preferred [57]. In our more complicated two-component system, the second criterion (A47) cannot be evaluated analytically. However, for any particular fluxtube arrangement, we can in principle compute ψ p ( 1 ) , then solve Equation (A36) to obtain γ , and thus test this criterion numerically. In particular, for the case of a square or hexagonal lattice, Equation (A36) can be solved in Fourier space, eventually leading to the result (50). We note that the perturbation to the neutron condensate produced by the fluxtubes, which is represented by γ in Equation (A47), always acts to reduce the overall free energy, making a first-order transition more likely. The magnitude of γ depends, via Equation (A36), on the coupling parameters h 1 and h 2 , and so a first-order transition is guaranteed if these parameters are sufficiently large.
Interestingly, there are certain combinations of the parameters for which Equation (A36) can be solved analytically. In particular, if α = 0 and
h 2 h 1 1 = 1 R 2 1 + h 3 ϵ 1 + h 1 ϵ ,
then we find that γ | ψ p ( 1 ) | 2 . In that case, the O ( δ a 2 ) term in the free energy (A45) has the form
2 π 1 + h 1 ϵ + π h 3 π κ 2 1 + h 1 ϵ 3 1 ϵ R 2 h 1 2 1 + h 1 ϵ 1 + 1 2 R 2 1 + h 3 ϵ 1 + h 1 ϵ β .
If one of the parameters ϵ , R or κ is sufficiently small, then the quantity in square brackets will be negative. Taken at face value, this result seems to suggest that the free energy can be made arbitrarily small, because β can be made arbitrarily large in an unbounded domain. However, this singularity actually just reflects the breakdown of our weakly nonlinear analysis in the limit of an infinite domain. To illustrate how this breakdown occurs, we have calculated the free energy of various multiply-charged hexagonal and square lattice states for one particular set of parameters, chosen such that the quantity in Equation (A49) is slightly negative. In Figure A1, we plot the resulting free energy as a function of a, for several different lattice types. As the value of a is reduced from 14 towards the critical value 2 π / H c 2 10.3 , we find that the singly-charged hexagonal lattice is replaced by the doubly-charged hexagonal lattice, and then by the triply-charged square lattice, as the energetically preferred lattice state. As a is further reduced, we expect that even more highly-charged lattice states (with higher values of β ) will become energetically preferred. However, throughout this whole range of a, the true ground state for an unbounded domain is not a lattice at all, and is instead a mixture of the non-superconducting state and a singly-charged hexagonal lattice. Our weakly-nonlinear analysis does not apply to the true ground state, which develops as a nonlinear instability for a > 2 π / H c 2 9.9 . In fact, for all parameter ranges we have studied, the fluxtubes always adopt a singly-charged, hexagonal pattern whenever the domain is sufficiently large, as seen in Figure 4.
Figure A1. The Helmholtz free energy per flux quantum per unit length, F ( a ) , for various lattice states, with (dimensionless) parameters h 1 = 0.061 , h 2 = 0.422 , h 3 = 0.062 , R = 0.412 , ϵ = 0.096 , and κ = 1.17 . The true ground state (dotted, gray line) arises as a subcritical bifurcation from the non-superconducting state (short-dashed, purple line). The two insets show more detail of the cross-over regions between two lattice configurations of different charge.
Figure A1. The Helmholtz free energy per flux quantum per unit length, F ( a ) , for various lattice states, with (dimensionless) parameters h 1 = 0.061 , h 2 = 0.422 , h 3 = 0.062 , R = 0.412 , ϵ = 0.096 , and κ = 1.17 . The true ground state (dotted, gray line) arises as a subcritical bifurcation from the non-superconducting state (short-dashed, purple line). The two insets show more detail of the cross-over regions between two lattice configurations of different charge.
Universe 08 00228 g0a1

References

  1. Baym, G.; Pethick, C.; Pines, D. Superfluidity in Neutron Stars. Nature 1969, 224, 673–674. [Google Scholar] [CrossRef]
  2. Mendell, G. Superfluid Hydrodynamics in Rotating Neutron Stars. I. Nondissipative Equations. Astrophys. J. 1991, 380, 515. [Google Scholar] [CrossRef]
  3. Sedrakian, D.M.; Sedrakian, A.D.; Zharkov, G.F. Type I superconductivity of protons in neutron stars. Mon. Not. R. Astron. Soc. 1997, 290, 203–207. [Google Scholar] [CrossRef] [Green Version]
  4. Buckley, K.B.; Metlitski, M.A.; Zhitnitsky, A.R. Vortices and type-I superconductivity in neutron stars. Phys. Rev. C 2004, 69, 055803. [Google Scholar] [CrossRef] [Green Version]
  5. Buckley, K.B.; Metlitski, M.A.; Zhitnitsky, A.R. Neutron stars as type-I superconductors. Phys. Rev. Lett. 2004, 92, 151102. [Google Scholar] [CrossRef] [Green Version]
  6. Alford, M.; Good, G.; Reddy, S. Isospin asymmetry and type-I superconductivity in neutron star matter. Phys. Rev. C 2005, 72, 055801. [Google Scholar] [CrossRef] [Green Version]
  7. Alford, M.G.; Good, G. Flux tubes and the type-I/type-II transition in a superconductor coupled to a superfluid. Phys. Rev. B 2008, 78, 024510. [Google Scholar] [CrossRef] [Green Version]
  8. Haber, A.; Schmitt, A. Critical magnetic fields in a superconductor coupled to a superfluid. Phys. Rev. D 2017, 95, 116016. [Google Scholar] [CrossRef] [Green Version]
  9. Babaev, E.; Speight, M. Semi-Meissner state and neither type-I nor type-II superconductivity in multicomponent superconductors. Phys. Rev. B 2005, 72, 180502. [Google Scholar] [CrossRef] [Green Version]
  10. Moshchalkov, V.; Menghini, M.; Nishio, T.; Chen, Q.H.; Silhanek, A.V.; Dao, V.H.; Chibotaru, L.F.; Zhigadlo, N.D.; Karpinski, J. Type-1.5 Superconductivity. Phys. Rev. Lett. 2009, 102, 117001. [Google Scholar] [CrossRef] [Green Version]
  11. Babaev, E.; Carlström, J.; Speight, M. Type-1.5 superconducting state from an intrinsic proximity effect in two-band superconductors. Phys. Rev. Lett. 2010, 105, 067003. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Gutierrez, J.; Raes, B.; Silhanek, A.V.; Li, L.J.; Zhigadlo, N.D.; Karpinski, J.; Tempere, J.; Moshchalkov, V.V. Scanning Hall probe microscopy of unconventional vortex patterns in the two-gap MgB2 superconductor. Phys. Rev. B 2012, 85, 094511. [Google Scholar] [CrossRef] [Green Version]
  13. Babaev, E.; Carlström, J.; Silaev, M.; Speight, J.M. Type-1.5 superconductivity in multicomponent systems. Physica C 2017, 533, 20–35. [Google Scholar] [CrossRef] [Green Version]
  14. Muzikar, P.; Pethick, C.J. Flux bunching in type-II superconductors. Phys. Rev. B 1981, 24, 2533–2539. [Google Scholar] [CrossRef]
  15. Sedrakian, A.D.; Sedrakian, D.M. Superfluid core rotation in pulsars. I. Vortex cluster dynamics. Astrophys. J. 1995, 447, 305. [Google Scholar] [CrossRef]
  16. Stairs, I.H.; Lyne, A.G.; Shemar, S.L. Evidence for free precession in a pulsar. Nature 2000, 406, 484–486. [Google Scholar] [CrossRef]
  17. Shabanova, T.V.; Lyne, A.G.; Urama, J.O. Evidence for free precession in the pulsar B1642-03. Astrophys. J. 2001, 552, 321–325. [Google Scholar] [CrossRef] [Green Version]
  18. Haberl, F.; Turolla, R.; de Vries, C.P.; Zane, S.; Vink, J.; Méndez, M.; Verbunt, F. Evidence for precession of the isolated neutron star RX J0720.4-3125. Astron. Astrophys. 2006, 451, L17–L21. [Google Scholar] [CrossRef]
  19. Ashton, G.; Jones, D.I.; Prix, R. Comparing models of the periodic variations in spin-down and beamwidth for PSR B1828-11. Mon. Not. R. Astron. Soc. 2016, 458, 881–899. [Google Scholar] [CrossRef]
  20. Sedrakian, A.; Wasserman, I.; Cordes, J.M. Precession of isolated neutron stars. I. Effects of imperfect pinning. Astrophys. J. 1999, 524, 341–360. [Google Scholar] [CrossRef] [Green Version]
  21. Link, B. Constraining hadronic superfluidity with neutron star precession. Phys. Rev. Lett. 2003, 91, 101101. [Google Scholar] [CrossRef] [Green Version]
  22. Sedrakian, A. Type-I superconductivity and neutron star precession. Phys. Rev. D 2005, 71, 083003. [Google Scholar] [CrossRef]
  23. Charbonneau, J.; Zhitnitsky, A. Novel mechanism for type I superconductivity in neutron stars. Phys. Rev. C 2007, 76, 015801. [Google Scholar] [CrossRef] [Green Version]
  24. Glampedakis, K.; Andersson, N.; Jones, D.I. Stability of precessing superfluid neutron stars. Phys. Rev. Lett. 2008, 100, 081101. [Google Scholar] [CrossRef] [Green Version]
  25. Gügercinoğlu, E.; Alpar, M.A. Vortex creep against toroidal flux lines, crustal entrainment, and pulsar glitches. Astrophys. J. Lett. 2014, 788, L11. [Google Scholar] [CrossRef] [Green Version]
  26. Graber, V.; Cumming, A.; Andersson, N. Glitch rises as a test for rapid superfluid coupling in neutron stars. Astrophys. J. 2018, 865, 23. [Google Scholar] [CrossRef]
  27. Sourie, A.; Chamel, N. Vortex pinning in the superfluid core of neutron stars and the rise of pulsar glitches. Mon. Not. R. Astron. Soc. 2020, 493, L98–L102. [Google Scholar] [CrossRef]
  28. Alpar, M.A. Flux-vortex pinning and neutron star evolution. J. Astrophys. Astron. 2017, 38, 44. [Google Scholar] [CrossRef] [Green Version]
  29. Drummond, L.V.; Melatos, A. Stability of interlinked neutron vortex and proton flux-tube arrays in a neutron star—II. Far-from-equilibrium dynamics. Mon. Not. R. Astron. Soc. 2018, 475, 910–920. [Google Scholar] [CrossRef] [Green Version]
  30. Dobaczewski, J.; Dudek, J. Time-odd components in the mean field of rotating superdeformed nuclei. Phys. Rev. C 1995, 52, 1827–1839. [Google Scholar] [CrossRef] [Green Version]
  31. Chamel, N.; Haensel, P. Entrainment parameters in a cold superfluid neutron star core. Phys. Rev. C 2006, 73, 045802. [Google Scholar] [CrossRef] [Green Version]
  32. Yakovlev, D.G.; Kaminker, A.D.; Gnedin, O.Y.; Haensel, P. Neutrino emission from neutron stars. Phys. Rep. 2001, 354, 1–155. [Google Scholar] [CrossRef] [Green Version]
  33. Yakovlev, D.G.; Pethick, C.J. Neutron star cooling. Annu. Rev. Astron. Astrophys. 2004, 42, 169–210. [Google Scholar] [CrossRef]
  34. Yakovlev, D.G.; Levenfish, K.P.; Shibanov, Y.A. Cooling of neutron stars and superfluidity in their cores. Phys. Uspekhi 1999, 42, 737–778. [Google Scholar] [CrossRef] [Green Version]
  35. Kaminker, A.D.; Haensel, P.; Yakovlev, D.G. Nucleon superfluidity vs. observations of cooling neutron stars. Astron. Astrophys. 2001, 373, L17–L20. [Google Scholar] [CrossRef] [Green Version]
  36. Kaminker, A.D.; Yakovlev, D.G.; Gnedin, O.Y. Three types of cooling superfluid neutron stars: Theory and observations. Astron. Astrophys. 2002, 383, 1076–1087. [Google Scholar] [CrossRef] [Green Version]
  37. Kobyakov, D.N. Surface energy of magnetized superconducting matter in neutron star cores. Phys. Rev. C 2020, 102, 045803. [Google Scholar] [CrossRef]
  38. Bedaque, P.F.; Rupak, G.; Savage, M.J. Goldstone bosons in the 3P2 superfluid phase of neutron matter and neutrino emission. Phys. Rev. C 2003, 68, 065802. [Google Scholar] [CrossRef] [Green Version]
  39. Alpar, M.A.; Langer, S.A.; Sauls, J.A. Rapid postglitch spin-up of the superfluid core in pulsars. Astrophys. J. 1984, 282, 533–541. [Google Scholar] [CrossRef]
  40. Babaev, E. Unconventional rotational responses of hadronic superfluids in a neutron star caused by strong entrainment and a Σ hyperon gap. Phys. Rev. Lett. 2009, 103, 231101. [Google Scholar] [CrossRef] [Green Version]
  41. Sinha, M.; Sedrakian, A. Magnetar superconductivity versus magnetism: Neutrino cooling processes. Phys. Rev. C 2015, 91, 035805. [Google Scholar] [CrossRef]
  42. Kobyakov, D.N.; Pethick, C.J. Two-component superfluid hydrodynamics of neutron star cores. Astrophys. J. 2017, 836, 203. [Google Scholar] [CrossRef] [Green Version]
  43. Drummond, L.V.; Melatos, A. Stability of interlinked neutron vortex and proton flux tube arrays in a neutron star: Equilibrium configurations. Mon. Not. R. Astron. Soc. 2017, 472, 4851–4869. [Google Scholar] [CrossRef]
  44. Silaev, M.; Babaev, E. Microscopic derivation of two-component Ginzburg–Landau model and conditions of its applicability in two-band systems. Phys. Rev. B 2012, 85, 134514. [Google Scholar] [CrossRef] [Green Version]
  45. Andreev, A.F.; Bashkin, E.P. Three-velocity hydrodynamics of superfluid solutions. J. Exp. Theor. Phys. 1976, 42, 164–167. [Google Scholar]
  46. Allard, V.; Chamel, N. 1S0 pairing gaps, chemical potentials and entrainment matrix in superfluid neutron-star cores for the Brussels–Montreal functionals. Universe 2021, 7, 470. [Google Scholar] [CrossRef]
  47. Graber, V.; Andersson, N.; Hogg, M. Neutron stars in the laboratory. Int. J. Mod. Phys. D 2017, 26, 1730015. [Google Scholar] [CrossRef]
  48. Esry, B.D.; Greene, C.H.; Burke, J.P., Jr.; Bohn, J.L. Hartree–Fock theory for double condensates. Phys. Rev. Lett. 1997, 78, 3594–3597. [Google Scholar] [CrossRef]
  49. Law, C.K.; Pu, H.; Bigelow, N.P.; Eberly, J.H. “Stability Signature” in two-species dilute Bose–Einstein condensates. Phys. Rev. Lett. 1997, 79, 3105–3108. [Google Scholar] [CrossRef]
  50. Bashkin, E.P.; Vagov, A.V. Instability and stratification of a two-component Bose–Einstein condensate in a trapped ultracold gas. Phys. Rev. B 1997, 56, 6207–6212. [Google Scholar] [CrossRef]
  51. Riboli, F.; Modugno, M. Topology of the ground state of two interacting Bose–Einstein condensates. Phys. Rev. A 2002, 65, 063614. [Google Scholar] [CrossRef] [Green Version]
  52. Tinkham, M. Introduction to Superconductivity, 2nd ed.; Dover Press: New York, NY, USA, 2004. [Google Scholar]
  53. Wood, T.S.; Mesgarnezhad, M.; Stagg, G.W.; Barenghi, C.F. Quasiperiodic boundary conditions for three-dimensional superfluids. Phys. Rev. B 2019, 100, 024505. [Google Scholar] [CrossRef] [Green Version]
  54. Kramer, L. Thermodynamic behavior of type-II superconductors with small κ near the lower critical field. Phys. Rev. B 1971, 3, 3821–3825. [Google Scholar] [CrossRef]
  55. Svistunov, B.V.; Babaev, E.S.; Prokof’ev, N.V. Superfluid States of Matter; CRC Press: Boca Raton, FL, USA, 2015. [Google Scholar]
  56. Abrikosov, A.A. On the magnetic properties of superconductors of the second group. J. Exp. Theor. Phys. 1957, 5, 1174–1182. [Google Scholar]
  57. Kleiner, W.H.; Roth, L.M.; Autler, S.H. Bulk solution of Ginzburg-Landau equations for type II superconductors: Upper critical field region. Phys. Rev. 1964, 133, 1226–1227. [Google Scholar] [CrossRef]
  58. Chamel, N. Two-fluid models of superfluid neutron star cores. Mon. Not. R. Astron. Soc. 2008, 388, 737–752. [Google Scholar] [CrossRef]
  59. Ho, W.C.G.; Elshamouty, K.G.; Heinke, C.O.; Potekhin, A.Y. Tests of the nuclear equation of state and superfluid and superconducting gaps using the Cassiopeia A neutron star. Phys. Rev. C 2015, 91, 015806. [Google Scholar] [CrossRef]
  60. Steiner, A.W.; Prakash, M.; Lattimer, J.M.; Ellis, P.J. Isospin asymmetry in nuclei and neutron stars [review article]. Phys. Rep. 2005, 411, 325–375. [Google Scholar] [CrossRef] [Green Version]
  61. Glampedakis, K.; Andersson, N.; Samuelsson, L. Magnetohydrodynamics of superfluid and superconducting neutron star cores. Mon. Not. R. Astron. Soc. 2011, 410, 805–829. [Google Scholar] [CrossRef] [Green Version]
  62. Pippard, A.B. The surface energies of superconductors. Proc. Camb. Phil. Soc. 1951, 47, 617. [Google Scholar] [CrossRef]
  63. Lewis, H.W. Surface Energies in Superconductors. Phys. Rev. 1956, 104, 942–947. [Google Scholar] [CrossRef]
  64. Jones, P.B. Neutron superfluid spin-down and magnetic field decay in pulsars. Mon. Not. R. Astron. Soc. 1991, 253, 279. [Google Scholar] [CrossRef]
  65. Gusakov, M.E. Force on proton vortices in superfluid neutron stars. Mon. Not. R. Astron. Soc. 2019, 485, 4936–4950. [Google Scholar] [CrossRef] [Green Version]
  66. Perez, F.; Granger, B.E. IPython: A system for interactive scientific computing. Comput. Sci. Eng. 2007, 9, 21–29. [Google Scholar] [CrossRef]
  67. Hunter, J.D. Matplotlib: A 2D Graphics Environment. Comput. Sci. Eng. 2007, 9, 90–95. [Google Scholar] [CrossRef]
  68. Oliphant, T.E. A Guide to NumPy; Trelgol Publishing: Spanish Fork, UT, USA, 2006. [Google Scholar]
  69. van der Walt, S.; Colbert, S.C.; Varoquaux, G. The NumPy array: A structure for efficient numerical computation. Comput. Sci. Eng. 2011, 13, 22–30. [Google Scholar] [CrossRef] [Green Version]
  70. Harris, C.R.; Jarrod Millman, K.; van der Walt, S.J.; Gommers, R.; Virtanen, P.; Cournapeau, D.; Wieser, E.; Taylor, J.; Berg, S.; Smith, N.J.; et al. Array programming with NumPy. Nature 2020, 585, 357–362. [Google Scholar] [CrossRef]
  71. McKinney, W. Data Structures for Statistical Computing in Python. In Proceedings of the 9th Python in Science Conference, Austin, TX, USA, 28 June–3 July 2010; van der Walt, S., Millman, J., Eds.; pp. 56–61. [Google Scholar]
  72. Jones, E.; Oliphant, T.E.; Peterson, P. SciPy: Open Source Scientific Tools for Python. 2001. Available online: https://www.scipy.org/ (accessed on 16 March 2020).
  73. Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; et al. SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nat. Methods 2020, 17, 261–272. [Google Scholar] [CrossRef] [Green Version]
  74. Taub, A.H. Stability of general relativistic gaseous masses and variational principles. Commun. Math. Phys. 1969, 15, 235–254. [Google Scholar] [CrossRef]
Figure 1. Phase diagrams for a one-component Ginzburg–Landau superconductor, for different values of the Ginzburg–Landau parameter, κ . The top panel shows the experiment with an imposed external field, | H | , in our nondimensional units. The first-order and second-order transitions at the different critical fields are indicated by solid and dashed black lines, respectively, and the resulting phases labeled accordingly. Shading of the respective regions is indicative of the magnetic flux distribution. The bottom panel shows the phase transitions in the experiment with an imposed mean flux, B ¯ . For more details see the text.
Figure 1. Phase diagrams for a one-component Ginzburg–Landau superconductor, for different values of the Ginzburg–Landau parameter, κ . The top panel shows the experiment with an imposed external field, | H | , in our nondimensional units. The first-order and second-order transitions at the different critical fields are indicated by solid and dashed black lines, respectively, and the resulting phases labeled accordingly. Shading of the respective regions is indicative of the magnetic flux distribution. The bottom panel shows the phase transitions in the experiment with an imposed mean flux, B ¯ . For more details see the text.
Universe 08 00228 g001
Figure 2. The Helmholtz free energy per flux quantum per unit length, F , as a function of the area per magnetic flux quantum, a, with the dimensionless parameters κ = 1.444 , R = 0.371 , ϵ = 0.097 , α = 0 , h 1 = 0.102 , h 2 = 0.387 , and h 3 = 0.263 . 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.
Figure 2. The Helmholtz free energy per flux quantum per unit length, F , as a function of the area per magnetic flux quantum, a, with the dimensionless parameters κ = 1.444 , R = 0.371 , ϵ = 0.097 , α = 0 , h 1 = 0.102 , h 2 = 0.387 , and h 3 = 0.263 . 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.
Universe 08 00228 g002
Figure 3. The ground-state Helmholtz free energy per unit length, F g ( a ) , can be used to infer the minimum Gibbs free-energy density, G ¯ , as a function of the external field, H . The tangent to each point on the curve F g ( a ) has slope G ¯ , and intersects the vertical axis at the point F = 4 π κ 2 | H | . The dotted lines indicate the transitions at H c 1 and H c 2 , both of which are first-order in this example. We have used the same dimensionless parameters as for Figure 2.
Figure 3. The ground-state Helmholtz free energy per unit length, F g ( a ) , can be used to infer the minimum Gibbs free-energy density, G ¯ , as a function of the external field, H . The tangent to each point on the curve F g ( a ) has slope G ¯ , and intersects the vertical axis at the point F = 4 π κ 2 | H | . The dotted lines indicate the transitions at H c 1 and H c 2 , both of which are first-order in this example. We have used the same dimensionless parameters as for Figure 2.
Universe 08 00228 g003
Figure 4. Inhomogeneous ground states for the parameters indicated by red dots in Figure 2. The brightness and hue indicate the density and phase of the proton order parameter, ψ p , respectively. As indicated by the color bar, the phase winds by 2 π around each fluxtube. The first panel shows a case with N = 24 magnetic flux quanta, and a (dimensionless) area of a N = 14.5 × 24 , corresponding to a mixture of non-superconducting protons and hexagonal fluxtube lattice. Approximately 2/3 of the magnetic flux is contained in the non-superconducting domain, visible as dark bands on both sides of the image, and hence only 8 fluxtubes are visible. The low brightness of the lattice domain indicates the low density of the proton condensate there, i.e., | ψ p | 2 1 . The second panel shows a case with N = 14 magnetic flux quanta, and a (dimensionless) area of a N = 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.
Figure 4. Inhomogeneous ground states for the parameters indicated by red dots in Figure 2. The brightness and hue indicate the density and phase of the proton order parameter, ψ p , respectively. As indicated by the color bar, the phase winds by 2 π around each fluxtube. The first panel shows a case with N = 24 magnetic flux quanta, and a (dimensionless) area of a N = 14.5 × 24 , corresponding to a mixture of non-superconducting protons and hexagonal fluxtube lattice. Approximately 2/3 of the magnetic flux is contained in the non-superconducting domain, visible as dark bands on both sides of the image, and hence only 8 fluxtubes are visible. The low brightness of the lattice domain indicates the low density of the proton condensate there, i.e., | ψ p | 2 1 . The second panel shows a case with N = 14 magnetic flux quanta, and a (dimensionless) area of a N = 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.
Universe 08 00228 g004
Figure 5. Phase diagrams for the two-component Ginzburg–Landau superconductor, for different values of the Ginzburg–Landau parameter, κ , and other parameters fixed to the same values as in Figure 2. Both figures are plotted in the same style as Figure 1. Shading of the respective regions is indicative of the magnetic flux distribution. In addition to the phases observed for the single-component case (we label the intermediate type-I phase as M/N), we also obtain inhomogeneous regimes where Meissner and fluxtube regions (M/F) as well as fluxtube and non-superconducting regions alternate. We refrain from labeling the latter to avoid overcrowding the plot. These are associated with the appearance of the critical fields H c 1 and H c 2 , respectively.
Figure 5. Phase diagrams for the two-component Ginzburg–Landau superconductor, for different values of the Ginzburg–Landau parameter, κ , and other parameters fixed to the same values as in Figure 2. Both figures are plotted in the same style as Figure 1. Shading of the respective regions is indicative of the magnetic flux distribution. In addition to the phases observed for the single-component case (we label the intermediate type-I phase as M/N), we also obtain inhomogeneous regimes where Meissner and fluxtube regions (M/F) as well as fluxtube and non-superconducting regions alternate. We refrain from labeling the latter to avoid overcrowding the plot. These are associated with the appearance of the critical fields H c 1 and H c 2 , respectively.
Universe 08 00228 g005aUniverse 08 00228 g005b
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wood, T.S.; Graber, V. Superconducting Phases in Neutron Star Cores. Universe 2022, 8, 228. https://0-doi-org.brum.beds.ac.uk/10.3390/universe8040228

AMA Style

Wood TS, Graber V. Superconducting Phases in Neutron Star Cores. Universe. 2022; 8(4):228. https://0-doi-org.brum.beds.ac.uk/10.3390/universe8040228

Chicago/Turabian Style

Wood, Toby S., and Vanessa Graber. 2022. "Superconducting Phases in Neutron Star Cores" Universe 8, no. 4: 228. https://0-doi-org.brum.beds.ac.uk/10.3390/universe8040228

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop