Critical magnetic fields in a superconductor coupled to a superfluid

TL;DR

This paper investigates how coupling a superconductor to a superfluid affects its magnetic properties, revealing a shift in phase transition types and the potential for flux tube clustering, relevant for neutron star cores.

Contribution

It introduces a temperature-dependent Ginzburg-Landau model for coupled superfluid-superconductor systems and analyzes the impact on critical magnetic fields and phase transitions.

Findings

Transition regions between type-I and type-II become first order due to superfluid coupling.

Presence of superfluid alters the nature of phase transitions at critical magnetic fields.

Flux tube clustering may occur in neutron star cores with coupled superfluid and superconductor.

Abstract

We study a superconductor that is coupled to a superfluid via density and derivative couplings. Starting from a Lagrangian for two complex scalar fields, we derive a temperature-dependent Ginzburg-Landau potential, which is then used to compute the phase diagram at nonzero temperature and external magnetic field. This includes the calculation of the critical magnetic fields for the transition to an array of magnetic flux tubes, based on an approximation for the interaction between the flux tubes. We find that the transition region between type-I and type-II superconductivity changes qualitatively due to the presence of the superfluid: the phase transitions at the upper and lower critical fields in the type-II regime become first order, opening the possibility of clustered flux tube phases. These flux tube clusters may be realized in the core of neutron stars, where superconducting protons are expected to be coupled to superfluid neutrons.

Figures

Equations

$$\mathcal{L}={\mathcal{L}}_1+{\mathcal{L}}_2+{\mathcal{L}}_{\mathrm{int}}+{\mathcal{L}}_{\mathrm{YM}},$$

$${\mathcal{L}}_i$$

$${\mathcal{L}}_{\mathrm{int}}$$

$${\mathcal{L}}_{\mathrm{YM}}$$

$$⟨{\phi }_i⟩=\frac{{\rho }_i}{\sqrt{2}}{e}^{-i{\psi }_i}.$$

$$q\equiv q_1,q_2=0.$$

$$U(\vec{r})$$

$$G\equiv \frac{g_1+g_2}{2},g\equiv \frac{g_1-g_2}{2}.$$

$$\frac{\partial U}{\partial {\rho }_1}=0,\frac{\partial U}{\partial {\rho }_2}=0,$$

$${\rho }_1={\rho }_2=0,U_{\mathrm{NOR}}=0.$$

$${\rho }_1^2$$

$${\rho }_2^2$$

$${\rho }_1^2$$

$$U_{\mathrm{COE}}$$

$$\Omega(\mu_{1},\mu_{2},T)=U+T\sum_{i=1}^{6}\int\frac{d^{3}k}{(2\pi)^{3}}\ln\big{(}1-e^{-\epsilon_{ki}/T}\big{)}\,,$$

$$U(\vec{r})$$

$$m_{1,T}^2$$

$$m_{2,T}^2$$

$$h_T$$

$${\rho }_{0i}^2(T)={\rho }_{0i}^2(T=0)\left(1-\frac{T^2}{T_{ci}^2}\right),$$

$$T_{c1}^2$$

$$T_{c2}^2$$

$$F=\int d^{3}rU(\vec{r}).$$

$$\mathcal{G}=F-\frac{\vec{H}}{4\pi }\cdot \int d^{3}r\vec{B}.$$

$${\mathcal{G}}_{\mathrm{COE}}=V{U}_{\mathrm{COE}},$$

$${\mathcal{G}}_{\mathrm{SF}}=V\left(U_{\mathrm{SF}}-\frac{H^2}{8\pi }\right),$$

$$H_c=\sqrt{8\pi (U_{\mathrm{SF}}-U_{\mathrm{COE}})}=2\pi q\sqrt{2}\kappa \sqrt{1-\frac{h_T^2}{{\lambda }_1{\lambda }_2}}{\rho }_{01}^2.$$

$$\kappa =\frac{\ell }{\xi }=\sqrt{\frac{{\lambda }_1}{4\pi q^2}},$$

$$\ell =\frac{1}{\sqrt{4\pi q^2}{\rho }_{01}},\xi =\frac{1}{\sqrt{{\lambda }_1}{\rho }_{01}}.$$

$$H_{c2}=\frac{1}{q{\xi }^2}\left(1-\frac{h_T^2}{{\lambda }_1{\lambda }_2}\right)=\sqrt{2}\kappa \sqrt{1-\frac{h_T^2}{{\lambda }_1{\lambda }_2}}{H}_c,$$

$$\frac{{\mathcal{G}}_{\mathrm{COE}}}{V}=\frac{{\mathcal{G}}_{\mathrm{SF}}}{V}+{\lambda }_1⟨{\bar{\phi }}_1^4⟩\left[\frac{1}{2{\kappa }^2}-1+\frac{h^2}{{\lambda }_1{\lambda }_2}{\mathcal{I}}_1(p)\right],$$

$$p^2=\frac{2{\lambda }_2{\rho }_{\mathrm{SF}}^2}{q{H}_{c2}},$$

$${\mathcal{I}}_1(p)\equiv \frac{p{e}^{p^2/4}}{2\sqrt{2}}{\int }_{-\infty }^{\infty }dt{e}^{-t^2}\left\{e^{pt}\left[1-\mathrm{erf}\left(\frac{p}{2}+t\right)\right]+e^{-pt}\left[1-\mathrm{erf}\left(\frac{p}{2}-t\right)\right]\right\},$$

Full text

I Introduction

I.1 Goal

An array of magnetic flux tubes is created in certain superconductors for intermediate strengths of an external magnetic field. Superconductors with this property are said to be of type II. This is in contrast to type-I superconductors, where the magnetic field is either completely expelled or completely destroys the superconducting state, but never penetrates partially through quantized flux tubes. The Ginzburg-Landau parameter $\kappa$ – the ratio between the magnetic penetration depth and the coherence length of the superconducting condensate – predicts whether a superconductor is of type I or of type II.

The goal of this paper is to study the critical magnetic fields for the flux tube lattice in a two-component system, where the superconductor is coupled to a superfluid. We consider a system of two complex scalar fields and an abelian gauge field, with the two scalar fields coupled to each other and one of them coupled to the gauge field – the neutral scalar field is then indirectly coupled to the gauge field through the charged scalar field. Various aspects of this system will be discussed, such as the effect of different forms of the coupling between the scalar fields (density coupling vs. derivative coupling), effects of nonzero temperature, and the interaction between magnetic flux tubes. Special emphasis will be put on the transition region between type-I and type-II behavior, because this region is changed qualitatively by the presence of the superfluid, and one of the main results will be the topology of the phase diagram in this region.

I.2 Methods

Our calculations are based on a Ginzburg-Landau free energy for two condensates. We start, however, from a field-theoretical Lagrangian from which we compute the thermal fluctuations of the system. This is necessary in order to generalize the standard temperature-dependent coefficients of the Ginzburg-Landau potential to the situation of two coupled fields. We shall work in a relativistic formalism, but the main results hold for non-relativistic systems as well because we only consider the static limit. The coupled equations of motion for the two condensates and the gauge field – which yield the profile and energy of a single flux tube – are computed numerically. Nevertheless, where possible, we derive simple analytical results. For instance, when we compute the free energy of a flux tube array, we employ an approximation valid for sparse arrays, based on the numerical solution for a single flux tube, which is sufficient to derive certain aspects of the phase structure. For a complete study of the phase diagram a fully numerical calculation would be necessary. We believe that our results provide guidance and physical insights that can support and complement such a numerical calculation in future studies.

I.3 Astrophysical context

A superconductor that is coupled to a superfluid is expected to exist in the core of neutron stars in the form of superconducting protons which coexist with superfluid neutrons Bogoliubov (1958); Migdal (1959); Page et al. (2014); Sedrakian and Clark (2006); Haber and Schmitt (2017). Although we keep all our results as generic as possible, this is the application we have in mind when we make certain choices for the parameters of our model. It is also the main motivation for including a derivative coupling between the superconductor and the superfluid; for a calculation of the strength of this coupling in dense nuclear matter see for instance Ref. Chamel and Haensel (2006). Microscopic calculations – which have to be taken with care at these extreme baryon number densities – suggest that the proton superconductor turns from type II to type I as the density increases, i.e., as we move further towards the center of the star. In other words, a neutron star has a spatially varying $\kappa$, and the transition from type-II to type-I superconductivity might be realized as a function of the radius of the star Glampedakis et al. (2011). The resulting interface between the two superconducting phases might affect the evolution of the magnetic field in the star and is thus of potential relevance to observations. Even if this interface is not realized, be it because the central density is not sufficiently large or because quark matter is preferred before the necessary density is reached, it is important to understand the magnetic properties of the flux tube phase in the presence of the neutron superfluid.

The energy gaps from nucleon Cooper pairing depend strongly on density, varying non-monotonically along the profile of the star, with a maximum of the order of $1\,{\rm MeV}$ at intermediate densities and being much smaller at higher densities deep in the core Wambach et al. (1993). Therefore, the critical temperatures, which can be as high as $T_{c}\sim 10^{10}\,{\rm K}$, are very small in certain regions of the star. And, the critical magnetic fields for proton superconductivity, at their maximum about $H_{c}\sim 10^{16}\,{\rm G}$ – larger than the largest measured surface fields – become very small as well. (A very feeble superconducting pairing gap is neither robust against temperature nor against a magnetic field.) This motivates us to study the behavior of the superconductor at magnetic fields close to the critical fields, and it motivates us to include temperature. For predictions in the astrophysical context, the coefficients of our effective model should be made density-dependent, using results from more microscopic calculations (which, however, are prone to large uncertainties). In the present work we mainly focus on deriving general results and only mimic the situation of dense nuclear matter by varying our parameters in a way that is reminiscent of the situation in a neutron star.

There are other possible two- or multi-fluid phases in the core of a neutron star, where at least one of the components is charged. For instance, hyperon condensation may yield further condensate species Gusakov et al. (2009), a charged hyperon condensate in coexistence with a proton superconductor possibly forming a two-superconductor system. Two-component systems are also possible in dense quark matter. In the color-flavor locked (CFL) phase Alford et al. (1999), the pairing of all quarks is usually described by a single gap function. This is different in the presence of a magnetic field, and the study of color-magnetic flux tubes Iida (2005) or domain walls Giannakis and Ren (2003) in a Ginzburg-Landau approach shows striking similarities with our two-component system. The color-magnetic flux tubes in CFL are not protected by topology Eto et al. (2014), but if there is a mechanism to stabilize them, for instance an external magnetic field, they may have interesting implications for neutron star physics Glampedakis et al. (2012), like their analogues in 2SC quark matter Alford and Sedrakian (2010). In coexistence with a kaon condensate Bedaque and Schäfer (2002); Alford et al. (2008a), the CFL phase couples a color superconductor with a superfluid and represents another interesting system to which our results can be potentially applied.

I.4 Broader context

A mixture of a superconductor with a superfluid is conceivable not only in neutron stars but also in the laboratory, for example in ultra-cold atomic systems, where Bose-Fermi mixtures have been produced Ferrier-Barbut et al. (2014); Delehaye et al. (2015). Atoms are, of course, neutral, and thus this is actually a mixture of two superfluids. However, at least for a single atomic species, the coupling to a “synthetic magnetic field” has been realized, including the observation of analogues of magnetic flux tubes Lin et al. (2009); Dalibard et al. (2011); Goldman et al. (2014). Therefore, future experiments may well allow for the creation of a laboratory version of a coupled superconductor/superfluid system.

If we relax the condition of exactly one of the two components being charged, we find more realizations. Systems of two superconducting components have been discussed in the literature Carlstrom et al. (2011); Brandt and Das (2011); Babaev and Silaev (2013); Wu et al. (2016) and can be realized in the form of two-band superconductors, or even in liquid metallic hydrogen Babaev et al. (2004). Two coexisting superfluids, besides atomic Bose-Fermi mixtures, are conceivable in 3He – 4He mixtures Tuoriniemi et al. (2002); Rysti et al. (2012), although in this case it is experimentally challenging to have both components in the superfluid state simultaneously.

I.5 Relation to previous work

Our study makes use of and extends various results of the literature. The model we are using is a gauged version of the one of Ref. Haber et al. (2016), where two-stream instabilities in a system of coupled superfluids were discussed. Magnetic flux tubes from proton superconductivity in neutron stars have been studied extensively in the literature, usually with an emphasis on phenomenological consequences. More microscopic approaches often do not include a consistent treatment of both components and rather put together separate results from the proton superconductor and the neutron superfluid (which may be a good approximation for certain quantities because of the small proton fraction in neutral, $\beta$-equilibrated nuclear matter). Studies relevant to our work that do include both components within a single model can be found in Refs. Alpar et al. (1984); Alford and Good (2008); Kobyakov and Pethick (2017); Sinha and Sedrakian (2015). In Ref. Alford and Good (2008), flux tube profiles and energies are computed, results that we reproduce and utilize in the present paper. Our calculation of the interaction between flux tubes is performed within an approximation valid for large flux tube separations, based on old literature for a single-component superconductor Kramer (1971); for a different method leading to the same result see Ref. Speight (1997). Extensions to a system of a superconductor coupled to a superfluid can be found in Refs. Buckley et al. (2004a, b), where the results were restricted to the symmetric situation of approximately equal self-coupling and cross-coupling strengths of the scalar fields (which is unrealistic for neutron star matter Alford et al. (2005)), and no derivative cross-coupling was taken into account. Interactions between flux tubes have also been computed, based on the same approximation, in the context of cosmic strings for one-component Bettencourt and Rivers (1995) and two-component MacKenzie et al. (2003) systems. Our study is also related to so-called type-1.5 superconductivity, predicted to occur in systems with two superconducting components Babaev and Speight (2005); Babaev et al. (2010); Carlstrom et al. (2011). Although in our study only one component is charged, we shall find very similar effects, for instance the possibility of flux tube clusters.

I.6 Structure of the paper

In Sec. II, we present the model, compute the free energy densities of the various phases at vanishing magnetic field, and introduce effects of nonzero temperature. In Sec. III, we derive the expressions for the critical magnetic fields $H_{c}$, $H_{c1}$, and $H_{c2}$ for our two-component system and use the flux tube - flux tube interaction to point out the possibilities of first-order phase transitions. Our numerical results, most of them in the form of phase diagrams, are presented in Sec. IV, together with a discussion of the type-I/type-II transition region. We give our conclusions and an outlook in Sec. V. Throughout the paper, we use natural units $\hbar=c=k_{B}=1$ and Gaussian units for the electromagnetic fields, such that the elementary charge is $e=\sqrt{\alpha}\simeq 0.085$ with the fine structure constant $\alpha$.

II Model

II.1 Lagrangian and basic phase structure

Our calculation will essentially be a mean-field Ginzburg-Landau study, and we could thus, as a starting point, simply state the Ginzburg-Landau potential. We choose a slightly more general field-theoretical language, mainly because it provides us with the framework of thermal field theory to introduce temperature. Starting from a Ginzburg-Landau potential directly, this would be less straightforward in our two-component system. In the following, we thus start with a Lagrangian for two complex scalar fields, and the zero-temperature Ginzburg-Landau potential simply is the tree-level potential of this Lagrangian. This is Eq. (5). Temperature is then introduced in an approximation based on the thermal excitations of the system, providing a simple temperature dependence for the Ginzburg-Landau coefficients, given in Eqs. (13) and (14).

The Lagrangian is

[TABLE]

where

[TABLE]

with the covariant derivative $D_{\mu}\varphi_{i}=(\partial_{\mu}+iq_{i}A_{\mu})\varphi_{i}$, where $A_{\mu}$ is the gauge field and $q_{1}$, $q_{2}$ the electric charges, with the complex scalar fields $\varphi_{1}$, $\varphi_{2}$, the mass parameters $m_{i}\geq 0$, the self-coupling constants $\lambda_{i}>0$, and the field strength tensor $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. We have included two types of cross-couplings between the fields: a density coupling with dimensionless coupling constant $h$, and a derivative coupling which allows for two different structures with coupling constants $g_{1}$ and $g_{2}$ of mass dimension $-2$. Due to this derivative coupling, the model is non-renormalizable and an ultra-violet cutoff is required in general. However, in our Ginzburg-Landau-like study we are only interested in an effective potential for which the only occurring momentum integral is made finite by nonzero temperature. Therefore, the non-renormalizability will not play any role in the following. The chemical potentials $\mu_{1}$ and $\mu_{2}$ are introduced in the usual way, they can be formally included in the Lagrangian as temporal components of the gauge fields in the covariant derivatives, $q_{i}A_{0}\to-\mu_{i}$, including the covariant derivatives in the coupling terms Haber et al. (2016). In isolation, each of the fields would form a Bose-Einstein condensate if $\mu_{i}>m_{i}$. We parametrize the condensates by their moduli $\rho_{i}$ and their phases $\psi_{i}$,

[TABLE]

Since we are interested in a superconductor coupled to a superfluid, we assume only one of the fields to be charged, say field 1, and the second to be neutral,

[TABLE]

Moreover, we are only interested in static solutions and thus drop all time derivatives. Then, the zero-temperature tree-level potential $U=-{\cal L}_{\varphi_{i}\to\langle\varphi_{i}\rangle}$ is

[TABLE]

where we have reduced the Yang-Mills contribution to a purely magnetic term, $\vec{B}=\nabla\times\vec{A}$, and where we have introduced the abbreviations

[TABLE]

Boundedness of the tree-level potential requires $h+g\mu_{1}\mu_{2}<\sqrt{\lambda_{1}\lambda_{2}}$. In the remainder of the paper, we shall set $g=0$, mainly for the sake of simplicity111In Ref. Alford and Good (2008) the terms proportional to $g$ were not included from the beginning. In Ref. Haber et al. (2016), which did not discuss vortex solutions, only the tree-level potential with $\nabla\rho_{i}=0$ was used, such that $G$ dropped out and $g$ was the only relevant derivative coupling. . Some of our results would become more complicated with a nonzero $g$, for example the large-temperature expansion in Sec. II.2. Also, by reducing the number of parameters, the parameter space of our model becomes a little less unwieldy. On the other hand, at least for zero temperature, $g$ does not play an important role for the magnetic flux tube profiles because our solutions will not include any circulation of the neutral condensate, $\nabla\psi_{2}=0$, in which case we see from Eq. (5) that $g$ appears merely as a modification of the density coupling $h$.

To define the possible phases of the system and establish the notation for their condensates, we start with the simplest case of spatially uniform condensates in the absence of a magnetic field, $\nabla\rho_{1}=\nabla\rho_{2}=\nabla\psi_{1}=\nabla\psi_{2}=\vec{A}=0$. As a consequence of these assumptions, the potential becomes independent of $G$. The local minima of the potential yield the possible phases, i.e., we need to solve the algebraic equations

[TABLE]

which allow for the following solutions.

• •

In the normal phase (“NOR”), neither the charged nor the neutral field condenses,

[TABLE]

• •

In the (pure) superconductor (“SC”), only the charged field forms a condensate, whereas the condensate of the other field is zero,

[TABLE]

• •

In the (pure) superfluid (“SF”), only the neutral field forms a condensate, while the charged fields remains uncondensed,

[TABLE]

• •

In the coexistence phase (“COE”), both condensates exist simultaneously. Without coupling, the coexistence phase is realized if and only if both chemical potentials are larger than the corresponding masses. The coupling favors ($h>0$) or disfavors ($h<0$) the COE phase. The condensates and the free energy density are

[TABLE]

The ground state is then found by determining the global minimum of $U$. The resulting phase diagram in the $\mu_{1}$-$\mu_{2}$ plane is shown in Fig. 1, for both signs of the coupling $h$. The figure also contains the phase transitions at nonzero temperature, which we discuss now.

II.2 Introducing temperature

We intend to include temperature $T$ into the potential (5) in an effective way. In Ginzburg-Landau models this is usually done by introducing $T$-dependent coefficients, with a $T$-dependence that is strictly valid only close to the critical temperature. In our system, the form of these coefficients is not obvious because we have two fields and hence (at least) two critical temperatures. We thus proceed by introducing temperature in our underlying field theory and derive an effective potential. This will be done in a high-temperature approximation, assuming the condensates to be uniform, and without background magnetic field. Once we have derived the $T$-dependent Ginzburg-Landau potential, we shall reinstate the magnetic field for our discussion of the phase diagram and allow for spatially varying condensates and gauge fields in a flux tube. Neglecting zero-temperature quantum corrections, the one-loop potential is

[TABLE]

where the sum is taken over all 6 quasiparticle excitations $\epsilon_{ki}$. Without condensation, each of the complex scalar fields yields 2 excitations (both massive if $m_{i}>0$), corresponding to particle and anti-particle excitations, while the gauge field has two massless excitations, corresponding to the two possible polarizations of massless photons. These are 6 modes in total. In the coexistence phase both scalar fields condense. As a consequence, there is one Goldstone mode from the neutral field and one would-be Goldstone boson from the charged field, which becomes a third mode of the now massive gauge field. Together with the two massive modes from the scalar fields and the two original modes of the gauge field – which are now massive as well – these are again 6 modes. The excitations $\epsilon_{ki}$ are computed from the tree-level propagator. Their expressions are very complicated, but for the high-$T$ approximation we only need their behavior at large momenta. All details of this calculation are deferred to appendix A. From a field-theoretical perspective our high-$T$ approximation is very crude, and for a quantitative evaluation of the model for all temperatures more sophisticated methods are needed, such as the two-particle irreducible formalism Alford et al. (2014) or functional renormalization group techniques Fejos and Hatsuda (2016). These methods are beyond the scope of the present work because, firstly, if applied to our present context of magnetic flux tubes and their interactions, they would render the calculation much more complicated and purely numerical methods would be required. Secondly, having in mind the application of our model to nuclear matter, the next step towards a more sophisticated description should probably be to employ a fermionic model, rather than improving the bosonic one (note for instance that our bosonic system has well-defined quasiparticle excitations for all energies, while a fermionic one has a continuous spectral density for energies larger than twice the energy gap from Cooper pairing).

We also simplify the result by only keeping the leading order contribution from the derivative coupling $G$. As a result, all temperature corrections can be absorbed into thermal masses and a thermal density coupling, and we can work with the effective potential

[TABLE]

where

[TABLE]

For the following, we can thus simply take Eqs. (8) – (11) and replace the masses and the density coupling by their thermal generalizations. The effect of nonzero temperature on the phase structure is shown in Fig. 1. Before we use the potential to compute the critical magnetic fields, we briefly comment on the critical temperatures of the coexistence phase without external magnetic field. In the presence of a derivative coupling $G$ the resulting expressions are very lengthy and not very insightful. Therefore, we set $G=0$ for the moment, such that the only effect of temperature is a modification of the masses $m_{1}$ and $m_{2}$. Inserting the thermal masses into Eqs. (11), we compute the $T$-dependent condensates

[TABLE]

where the critical temperatures $T_{c1}$ and $T_{c2}$ indicate the phase transitions to the SF and SC phases,

[TABLE]

In the limit $h=0$, Eq. (16a) reduces to the well-known result for a single charged field, see for instance Eq. (4.24) in Ref. Kapusta (1981) (in this reference Heaviside-Lorentz units are used, i.e., our charge $q$ has to be divided by $\sqrt{4\pi}$ to match that result exactly). If we set $h=0$ in Eq. (16b) the result becomes independent of the charge $q$, as it should be because field 2 is neutral and couples to the gauge field only indirectly through field 1.

The critical temperatures (16) and their more complicated versions with nonzero $G$ are interesting in themselves. For instance, they can be used to analyze systematically in which regions of parameter space the COE phase is superseded by the SF phase at high temperature (i.e., the charged condensate melts first, $T_{c1}<T_{c2}$) or by the SC phase (i.e., the neutral condensate melts first, $T_{c2}<T_{c1}$). Or, they can be used to identify regions in the parameter space where one or both critical temperatures squared become negative, indicating that one or both condensates “refuse” to melt. This interesting observation – although it may be an artifact of our approximation – has been pointed out previously in the literature, see for instance appendix C in Ref. Alford et al. (2008b) and references therein. Here we shall not further analyze the critical temperatures and proceed with our main concern, phases at nonzero external magnetic field. None of the parameter sets we shall use in the following show this unusual behavior, i.e., we choose parameters such that $T_{c1}$ and $T_{c2}$ exist.

III Critical magnetic fields

The free energy can be computed from the potential (13),

[TABLE]

Since we are interested in the phase structure at fixed external (and homogeneous) magnetic field $\vec{H}=H\vec{e}_{z}$, we need to consider the Gibbs free energy

[TABLE]

To determine the complete phase diagram, we would have to compute the Gibbs free energy for all possible phases at each point in the phase space given by the thermodynamic variables $(\mu_{1},\mu_{2},T,H)$. The possible phases are the NOR, SF, SC, and COE phases listed above, and for the phases that are superconducting (SC and COE) we have to distinguish the Meissner phase, in which the magnetic field is completely expelled, $\vec{B}=0$, from the flux tube phase, where a lattice of magnetic flux tubes is formed, admitting part of the applied magnetic field in the superconductor. We shall simplify this problem by not computing the Gibbs free energy for the flux tube phase in full generality, which would require us to determine the spatial profile of the condensate and the magnetic field, including the preferred lattice structure, fully dynamically. Instead – following the usual textbook treatment Tinkham (2004) – we shall compute the critical magnetic fields $H_{c1}$, $H_{c2}$, and $H_{c}$, although they do not provide complete information of the phase diagram, not even for a single-component superconductor. To interpret their meaning for the phase diagram (in particular in our two-component system) it is important to precisely recall how they are computed, and thus we start each of the following three subsections with the definition of the corresponding critical magnetic field before we compute them for our system. In general, when we speak of the superconducting phase, this can be either the COE or the SC phase, while the normal-conducting phase can either be NOR or SF. The concrete calculations will always be done for the most interesting case, where both charged and neutral condensates exist in the superconducting phase (COE) and the normal conductor is the pure superfluid (SF). The critical magnetic fields for the transition between the COE and NOR and between the SC and NOR phases are not needed for our main results, but can be computed analogously. The latter appears to be the standard textbook scenario. However, in our two-component system it is conceivable that in the SC phase a neutral condensate is induced in the center of a flux tube Forgacs and Luk ács (2016); Forga cs and Luk ács (2016). Therefore, the pure superconductor SC might acquire a superfluid admixture, which can affect the critical magnetic fields for the transition to the completely uncondensed phase (NOR). In the present paper, we shall only consider flux tube solutions that approach the COE phase, not the SC phase, far away from the center of the flux tube.

III.1 Critical magnetic field $H_{c}$

Definition. The critical magnetic field $H_{c}$ is the magnetic field at which the Gibbs free energies of the superconducting phase in the Meissner state and the normal-conducting phase are identical, resulting in a first-order phase transition between them.

The Gibbs free energy of the COE phase with complete expulsion of the magnetic field is

[TABLE]

where $V$ is the total volume of the system and $U_{\rm COE}$ is the free energy density from Eq. (11b), with the masses $m_{1}$, $m_{2}$ and the coupling $h$ replaced by their thermal generalizations $m_{1,T}$, $m_{2,T}$, $h_{T}$. We neglect any magnetization in the normal-conducting phases, and thus $\vec{B}=\vec{H}$ in the SF phase, which yields the Gibbs free energy

[TABLE]

with $U_{\rm SF}$ from Eq. (10). Note that the $H^{2}$ term is a sum of the magnetic energy $\propto B^{2}$ and the term $\propto HB$ in the Legendre transformation from the free energy $F$ to the Gibbs free energy ${\cal G}$. Therefore, the critical magnetic field, defined by ${\cal G}_{\rm COE}={\cal G}_{\rm SF}$, becomes

[TABLE]

Here we have introduced the Ginzburg-Landau parameter

[TABLE]

with the magnetic penetration depth $\ell$ and the coherence length $\xi$,

[TABLE]

III.2 Critical magnetic field $H_{c2}$

Definition. Suppose there is a second-order phase transition between the superconductor in the flux tube phase and the normal-conducting phase, such that the equations of motion can be linearized in the charged condensate. Then, the critical magnetic field $H_{c2}$ is the maximal magnetic field allowed by the equations of motion. $H_{c2}$ is a lower bound for the actual transition from the flux tube phase to the normal-conducting phase because it does not exclude a first-order transition at some larger $H$. We call the critical field for such a first-order transition $H_{c2}^{\prime}$.

By definition, as we approach $H_{c2}$, the charged condensate approaches zero and the neutral condensate approaches the condensate of the SF phase. For magnetic fields $H$ close to and smaller than $H_{c2}$, we can write the condensates and the gauge field as their values at $H_{c2}$ plus small perturbations. Then, for the calculation of $H_{c2}$ itself the equations of motion linear in the charged condensate are sufficient. We are also interested in checking whether and in which parameter regime the flux tube phase is energetically preferred just below $H_{c2}$. This is done within the same calculation, but taking into account higher order terms in the equations of motion and the free energy. This calculation is somewhat lengthy and is explained in appendix B. Here we summarize the results. The critical magnetic field becomes

[TABLE]

where the second expression relates $H_{c2}$ to $H_{c}$ by using Eq. (21). At zero temperature, $H_{c2}$ does not depend on the gradient coupling $G$. However, the difference in Gibbs free energies between the superconducting and the normal-conducting phases does depend on $G$, see Eq. (94). For $G=0$ we have

[TABLE]

where $\langle\bar{\varphi}_{1}^{4}\rangle$ is the spatial average of the charged condensate (87), where

[TABLE]

and where

[TABLE]

with the error function erf.

In the limit of a single superconductor, $h=0$, we recover the standard result: in that case, Eq. (24) shows that the critical fields $H_{c}$ and $H_{c2}$ coincide at $\kappa^{2}=1/2$, and Eq. (25) shows that the flux tube phase is preferred, ${\cal G}_{\rm COE}<{\cal G}_{\rm SF}$, if and only if $\kappa^{2}>1/2$. In the coupled system the situation is more complicated. Now, from Eq. (24) we see that $H_{c}$ and $H_{c2}$ coincide at a larger value of $\kappa$ (since $h^{2}<\lambda_{1}\lambda_{2}$ to ensure the boundedness of the potential for $h>0$ and to ensure the existence of the COE phase for $h<0$, the square root is always real and smaller than 1). This appears to take away phase space from the flux tube phase. However, from Eq. (25) we see that the difference in Gibbs free energies between the COE and the SF phases changes sign at a different point, and this point is given not just by the coupling constant $h$, but also depends on $p$, i.e., on the magnitude of the neutral condensate $\rho_{\rm SF}$ compared to the square root of the critical magnetic field $H_{c2}$. Despite this dependence we can make a general statement: we find $0\leq{\cal I}_{1}(p)<1$, and thus the factor ${\cal I}_{1}(p)$ weakens the effect of the term $h^{2}/(\lambda_{1}\lambda_{2})$. At the value of $\kappa$ where $H_{c}$ and $H_{c2}$ are equal, the superconducting phase is preferred and – for all $p$ – remains preferred along $H_{c2}$, until the smaller $\kappa$ defined through Eq. (25) is reached. This observation is indicative of the complications at the transition between type-I and type-II superconductivity in the two-component system, and we shall find further discrepancies to the standard scenario when we compute the critical field $H_{c1}$.

Anticipating the numerical results in Sec. IV, let us comment on a possible first-order phase transition at $H^{\prime}_{c2}$, as mentioned in the definition at the beginning of this section. Suppose we are in a parameter region where the flux tube phase is favored just below $H_{c2}$, i.e., let $\kappa$ be larger than the critical $\kappa$ defined through Eq. (25). Then, any phase transition from the flux tube phase to the normal phase at a critical field smaller than $H_{c2}$ is excluded because we know that the system prefers to be in the flux tube phase just below $H_{c2}$ (here we ignore the very exotic possibility that the system quits the flux tube phase and then re-enters it below $H_{c2}$). A phase transition at a critical magnetic field larger than $H_{c2}$ – instead of the one at $H_{c2}$ – is however possible. This phase transition must be of first order because by definition $H_{c2}$ is the largest magnetic field at which a second-order transition may occur. Putting these arguments together leads to the conclusion that $H_{c2}$ is a lower bound for the transition from the flux tube phase to the normal phase, possibly replaced by a first order transition at $H^{\prime}_{c2}>H_{c2}$. Our numerical results will indeed suggest such a first-order phase transition. However, we shall find $H^{\prime}_{c2}<H_{c2}$, which, as we will explain, is an artifact of the approximation we apply for the interaction between flux tubes. Nevertheless, our result will allow us to speculate about the correct critical field $H^{\prime}_{c2}$, obtained in a more complete calculation that goes beyond our approximation.

III.3 Critical magnetic field $H_{c1}$

Definition. The critical magnetic field $H_{c1}$ is the magnetic field at which it becomes energetically favorable to put a single flux tube into the superconductor in the Meissner phase, resulting in a second-order phase transition from the Meissner phase into the flux tube phase. $H_{c1}$ is an upper bound for this transition because there can be a first-order transition at some smaller $H$, i.e., it can be favorable to directly form a flux tube lattice with a finite, not infinite, distance between the flux tubes. We call this first-order critical field $H_{c1}^{\prime}$.

According to the definition (18), the Gibbs free energy for the COE phase with a single magnetic flux tube is

[TABLE]

where $F_{\circlearrowleft}$ is the free energy of the flux tube, and where we have used

[TABLE]

with the winding number $n$ of the flux tube, the length of the flux tube $L$, and the fundamental flux quantum $\Phi_{0}=2\pi/q$. Placing a single flux tube into the system results in a loss in (negative) condensation energy, and thus the free energy increases. However, at fixed magnetic field $H$, there is an energy gain from allowing magnetic flux into the system. As a consequence, there is a competition between these two contributions of opposite sign in Eq. (28). At the critical point, the two contributions exactly cancel each other,

[TABLE]

The calculation of $H_{c1}$ thus amounts to the calculation of the free energy of a single flux tube $F_{\circlearrowleft}$, for which we can largely follow Ref. (Alford and Good, 2008). We work in cylindrical coordinates, $\vec{r}=(r,z,\theta)$, and make the following, radially symmetric, ansatz for the condensates,

[TABLE]

and the gauge field

[TABLE]

The profile functions $f_{i}$ and $a$ have to be computed numerically. Their boundary conditions are $f_{i}(\infty)=a(\infty)=1$, $f_{1}(0)=0$, and $\partial_{r}f_{2}(\infty)=\partial_{r}a(\infty)=0$, such that the condensates approach their homogeneous values $\rho_{0i}$ and the magnetic field vanishes far away from the center of the flux tube. The values of the neutral condensate and the gauge field at the center of the flux tube are determined dynamically. We have set the winding number of the neutral condensate to zero because the flux tube does not induce a superfluid vortex (Alford and Good, 2008).

We insert our ansatz into the potential (13) and separate the potential of the homogeneous COE phase,

[TABLE]

with

[TABLE]

To write the free energy of the flux tube in a convenient form, we introduce the dimensionless variable

[TABLE]

abbreviate the dimensionless gradient coupling by

[TABLE]

and the ratio of neutral over charged condensate by

[TABLE]

It is also useful to write $\mu_{1}^{2}-m_{1,T}^{2}=\lambda_{1}\rho_{\rm SC}^{2}=\lambda_{1}\rho_{01}^{2}-h_{T}\rho_{02}^{2}$ and $\mu_{2}^{2}-m_{2,T}^{2}=\lambda_{2}\rho_{\rm SF}^{2}=\lambda_{2}\rho_{02}^{2}-h_{T}\rho_{01}^{2}$, which follows from Eq. (11). Then, we obtain the free energy per unit length

[TABLE]

where prime denotes derivative with respect to $R$. This yields the equations of motion for $a$, $f_{1}$, $f_{2}$,

[TABLE]

We solve these equations numerically with a successive over-relaxation method. The profiles themselves have been discussed in detail in Ref. Alford and Good (2008)222Eqs. (39) are identical to Eqs. (16) in Ref. Alford and Good (2008) if we identify

$\frac{\Gamma}{2}\leftrightarrow\sigma\,,\qquad x\leftrightarrow\frac{\langle\phi_{n}\rangle}{\langle\phi_{p}\rangle}\,,\qquad\frac{h_{T}}{\lambda_{1}}\leftrightarrow-\frac{a_{pn}}{a_{pp}}\,,\qquad\frac{\lambda_{2}}{\lambda_{1}}\leftrightarrow\frac{a_{nn}}{a_{pp}}\,.$

, and we do not further comment on them. Instead we continue with the asymptotic solution, which will be needed later.

Far away from the center of the flux tube, all profile functions are close to one. Therefore, we write

[TABLE]

and linearize the profile equations (39) in $v$, $u_{1}$, and $u_{2}$,

[TABLE]

where

[TABLE]

We can decouple the equations for $u_{1}$ and $u_{2}$ by diagonalizing $M$,

[TABLE]

where $\nu_{\pm}$ are the eigenvalues of $M$ and $(\gamma_{\pm},1)$ its eigenvectors, given by

[TABLE]

where ${\cal Q}\equiv[(\lambda_{1}-\lambda_{2}x^{2})^{2}-2h_{T}\Gamma x(\lambda_{1}+\lambda_{2}x^{2})+x^{2}(4h_{T}^{2}+\Gamma^{2}\lambda_{1}\lambda_{2})]^{1/2}$. This yields two uncoupled equations for $\tilde{u}_{1}$ and $\tilde{u}_{2}$, where $\tilde{u}=U^{-1}u$, which we solve with the boundary condition $\tilde{u}_{1}(\infty)=\tilde{u}_{2}(\infty)=0$ (which leaves one integration constant from each equation undetermined). We undo the rotation with $u=U\tilde{u}$, and, together with the solution to Eq. (41a), insert the result into Eq. (40) to obtain the asymptotic solutions

[TABLE]

where $K_{0}$ and $K_{1}$ are the modified Bessel functions of the second kind, and the constants $C$, $D_{+}$, $D_{-}$ can only be determined numerically by solving the full equations of motion, including the boundary conditions at $R=0$. In deriving the linearized equations (41), we have not only used $u_{1},u_{2},v\ll 1$, but also $v^{2}\ll u_{1},u_{2}$, which implies $e^{-2R/\kappa}\ll e^{-\sqrt{\nu_{\pm}}R}$. This assumption is violated if $\kappa$ is sufficiently large compared to $1/\sqrt{\nu_{\pm}}$ (compared to $1/\sqrt{2}$ in a single superconductor), i.e., deep in the type-II regime. Later, when we use the solutions of the linearized equations for the interactions between flux tubes, we are only interested in the transition region between type-I and type-II behavior, where $1/\kappa\simeq\sqrt{\nu_{\pm}}$, i.e., for our purpose the linearization is a valid approximation.

III.4 Interaction between flux tubes and first-order phase transitions

If the phase transitions from the Meissner phase to the flux tube phase and from the flux tube phase to the normal-conducting phase were of second order we would be done. The critical magnetic fields of the previous sections would be sufficient to determine the phase structure. We shall see, however, that, due to the presence of the superfluid, first-order phase transitions become possible. To this end, we compute the Gibbs free energy of the entire flux tube lattice, rather than only of a single flux tube. We shall do so in an approximation of flux tube distances much larger than the width of a flux tube.

We generalize the Gibbs free energy (28) to a system with flux tube area density $\nu$ and add a term that takes into account the interaction between the flux tubes

[TABLE]

where we have eliminated $F_{\circlearrowleft}$ in favor of $H_{c1}$ with the help of Eq. (30), and where we have employed the nearest-neighbor approximation for the interaction term with the number of nearest neighbors $t$, and the dimensionless lattice constant $R_{0}$. For a hexagonal lattice, which we shall use in our explicit calculation, $t=6$ and $\nu=2/(\sqrt{3}R_{0}^{2})$. The interaction energy $F_{\rm int}^{\circlearrowleft}(R_{0})$ is defined by writing the total free energy of two flux tubes with distance $R_{0}$, say flux tubes $(a)$ and $(b)$, in terms of the free energy of the flux tubes in isolation plus the interaction energy,

[TABLE]

We calculate $F_{\rm int}^{\circlearrowleft}(R_{0})$ in appendix C in an approximation that is valid for large $R_{0}$. This calculation makes use of the method first employed in Ref. Kramer (1971), adapted to our two-component system with gradient coupling. All related references mentioned in Sec. I.5 are based on this method or an equivalent one, and our results reproduce the ones of those references in various limits. The result is

[TABLE]

As explained in the appendix in more detail, the integration can be reduced to an integral over the plane that separates the two Wigner-Seitz cells, which, in this simple setup, are two half-spaces. Since the integration along the direction of the flux tubes is trivial, we are left with a one-dimensional integral. As a consequence of the approximation, only the profile functions of a single flux tube appear in the integrand. In the derivation we have also assumed the asymptotic values of the condensates to be identical to the homogeneous values in the Meissner phase, $\rho_{01}$ and $\rho_{02}$. We shall later insert our numerical solutions $f_{1}$, $f_{2}$, and $a$ into Eq. (48) to compute the Gibbs free energy numerically. Before we do so we extract some simple analytical results with the help of the asymptotic solutions (45). Inserting them into Eq. (48) yields a lengthy expression which is not very instructive, especially due to the terms proportional to the gradient coupling. In appendix D we show that a simple expression can be extracted, even including the gradient coupling, if we restrict ourselves to the leading order contribution at large distances. Here we proceed with the simpler case of vanishing gradient coupling, $\Gamma=0$, to obtain straightforwardly

[TABLE]

where we have used $\gamma_{+}\gamma_{-}+x^{2}=0$ for $\Gamma=0$, which follows from Eqs. (44), the derivatives $K_{1}^{\prime}(x)=-K_{0}(x)-K_{1}(x)/x$, $K_{0}^{\prime}(x)=-K_{1}(x)$, and the integral

[TABLE]

The result (49) shows that there is a positive contribution, which makes the flux tubes repel each other due to their magnetic fields, and there is a negative contribution, which makes the flux tubes attract each other due to the lower loss of (negative) condensation energy if the flux tubes overlap. Let us first see how the case of a single superfluid is recovered by switching off the coupling $h$. (Since we have set $\Gamma=0$, there is no temperature dependence left in $h_{T}$ and we drop the subscript $T$ in this discussion.) As $h\to 0$, the quantities $\nu_{\pm}$ and $\gamma_{\pm}$ go to different limits, depending on the sign of $\lambda_{1}-\lambda_{2}x^{2}$. If $\lambda_{2}x^{2}>\lambda_{1}$, we have $\gamma_{+}\sim h$ and $\gamma_{-}\sim h^{-1}$. Numerically, we find that while $\gamma_{-}$ diverges, the product $D_{-}\gamma_{-}$ goes to a finite value. Moreover, $D_{+}$ goes to zero, such that the attractive terms reduce to $-D_{-}^{2}\gamma_{-}^{2}K_{0}(R_{0}\sqrt{2})$ since $\nu_{-}\to 2$ for $h\to 0$. In particular, all dependence on $x$, which contains the neutral condensate, has disappeared, as it should be. If, on the other hand, $\lambda_{2}x^{2}<\lambda_{1}$, we see from Eqs. (44) that now $\gamma_{+}\sim h^{-1}$ and $\gamma_{-}\sim h$, and it is the other term, $-D_{+}^{2}\gamma_{+}^{2}K_{0}(R_{0}\sqrt{2})$, which survives, again reproducing the correct result of a single superconductor. The result can be used to find the sign of the interaction at $R_{0}\to\infty$, i.e., to determine whether the flux tubes repel or attract each other at large distances. Since the Bessel functions fall off exponentially for large $R_{0}$, we simply compare the arguments of the Bessel functions of the negative and positive contributions. For the single superconductor, the long-distance flux tube interaction is thus attractive for $\kappa^{2}<1/2$ and repulsive for $\kappa^{2}>1/2$, i.e., the sign change appears exactly at the point where $H_{c}=H_{c2}$.

Going back to the full expression (49) for the two-component system, we compare $\nu_{-}$ with $1/\kappa^{2}$, because $\nu_{-}<\nu_{+}$, i.e., the term proportional to $K_{0}(R_{0}\sqrt{\nu_{-}})$ is less suppressed for $R_{0}\to\infty$. Therefore, the point at which the long-range interaction changes from repulsive to attractive is given by

[TABLE]

By comparing Eq. (51) with Eq. (24), we see that in the two-component system the long-distance interaction changes its sign at a point different from $H_{c}=H_{c2}$. This is made particularly obvious in the second line of Eq. (51), where we have expanded the result for large values of $x$, i.e., for large values of the neutral condensate compared to the charged one, $\rho_{02}/\rho_{01}\gg 1$. This limit is interesting for the interior of neutron stars, where protons are expected to contribute only about 10% to the total baryon number density333In Ref. Buckley et al. (2004b), the limit $x\gg 1$ was considered ($n_{1}/n_{2}\ll 1$ in the notation of that reference), and it was argued that the critical $\kappa$’s for $H_{c}=H_{c2}$ and the sign change of the long-range interaction are identical, in agreement with the leading-order contribution of our Eq. (51). Ref. Buckley et al. (2004b) only considered the near-symmetric situation $\lambda_{1}=\lambda_{2}\equiv\lambda$, $h=-\lambda+\delta\lambda$ with $0<\delta\lambda\ll\lambda$ (notice that $h<0$ here). In this case, our results show that $H_{c}=H_{c2}$ occurs at $\kappa^{2}\simeq\frac{\lambda}{4\delta\lambda}$ and the sign change in the long-range interaction energy at $\kappa^{2}\simeq\frac{\lambda}{4\delta\lambda}\frac{1+x^{2}}{x^{2}}$. Consequently, even in the near-symmetric situation the two critical $\kappa$’s are different and only become identical in the limit $x\gg 1$. . From Eq. (51) we recover $\kappa^{2}=1/2$ for $h=0$, but only if $\lambda_{2}x^{2}>\lambda_{1}$. The reason is that the limits $R_{0}\to\infty$ and $h\to 0$ do not commute in general: in deriving Eq. (51) we have fixed $h$ at a nonzero value and let $R_{0}\to\infty$, while in our above discussion of the single superconductor, we fixed $R_{0}$ while first letting $h\to 0$.

An attractive long-distance interaction between the flux tubes can have very interesting consequences. Recall that $H_{c1}$ is the magnetic field at which the phase with a single flux tube is preferred over the phase with complete field expulsion. In other words, at $H_{c1}$ the flux tube density is zero and increases continuously, while the flux tube distance decreases continuously from infinity at $H_{c1}$. If the interaction at infinite distances is attractive, the flux tubes do not “want” to form an array with arbitrarily small density. Assuming that the interaction always becomes repulsive at short range [which our numerical results confirm if we extrapolate Eq. (48) down to lower distances], there is a minimum in the flux tube - flux tube potential, which corresponds to a favored distance between the flux tubes. As a consequence, the transition from the Meissner phase to the flux tube phase occurs at a critical field lower than $H_{c1}$, which we call $H_{c1}^{\prime}$, at which the flux tube density jumps from zero to a nonzero value. An instructive analogy is the onset of nuclear matter as a function of the baryon chemical potential $\mu_{B}$. If the nucleon - nucleon potential was purely repulsive, there would be a second-order onset at the baryon mass, $\mu_{c}=m_{B}$. In reality, there is a binding energy $E_{b}$, and the baryon onset is a first-order transition at a lower chemical potential $\mu_{c}^{\prime}=m_{B}-E_{b}$. Here, the role of the chemical potential is played by the external field $H$, the role of the nucleons is played by the flux tubes with mass per unit length $H_{c1}=2qF_{\circlearrowleft}/(nL)$, and the binding energy is generated by the attractive interaction between the flux tubes.

In the single-component system, this first-order phase transition is not realized because it occurs in the type-I regime. More precisely, if we were to continue $H_{c1}$ into the type-I regime, then, at $H_{c1}$, it does not matter that the flux tube phase is made more favorable by an attractive interaction because the normal-conducting phase is the ground state (under the assumption that the gain in Gibbs free energy is not sufficient to overcome the difference to the normal phase). In the two-component system, however, the attractive interaction may exist in the regime where, at $H_{c1}$, the Meissner phase (and the phase with a single flux tube) is already preferred over the normal phase. Hence, any arbitrarily small binding energy will lead to a first-order phase transition at $H_{c1}^{\prime}<H_{c1}$. As we move along $H_{c1}$ towards smaller values of $\kappa$, i.e., towards the type-I regime, we hit the critical point given by Eq. (51), where the second-order transition turns into a first-order transition. Since our approximation is accurate for infinitesimally small flux tube densities, our prediction for this point is exact. If we then keep moving along $H_{c1}^{\prime}$, the flux tube density at the transition increases and our results have to be taken with care.

We can directly compute $H_{c1}^{\prime}$ by equating the Gibbs free energy of the flux tube phase (46) to the Gibbs free energy of the Meissner phase (19). In the flux tube phase we have to find the preferred flux tube distance $R_{0}$ (or, equivalently, the preferred flux tube density $\nu$), which is given by minimizing the Gibbs free energy. Hence, we compute $H_{c1}^{\prime}$ by solving the coupled equations

[TABLE]

for $H$ and $R_{0}$. We may use the same method to compute a potential first-order phase transition from the flux tube phase to the normal-conducting phase, i.e., in the free energy comparison we replace ${\cal G}_{\rm COE}$ with ${\cal G}_{\rm SF}$ from Eq. (20) and compute the resulting critical field $H_{c2}^{\prime}$.

IV Phase diagrams

IV.1 Taming the parameter space

The results in the previous sections have shown that the presence of the superfluid affects the transition from type-I to type-II superconductivity in a qualitative way, and we will make these results now more concrete by discussing the phase diagram of our model. To this end, we need to locate this transition in the parameter space. A priori, we have to deal with a large number of parameters, $m_{1}$, $m_{2}$, $\lambda_{1}$, $\lambda_{2}$, $q$, $h$, $G$, and the external thermodynamic parameters $T$, $H$, $\mu_{1}$, $\mu_{2}$. Having in mind a system of neutron and proton Cooper pairs we set $m_{1}=m_{2}\equiv m$ and $q=2e$, and express all dimensionful quantities in units of $m$. Many interesting results can already be obtained with a density coupling alone, and we shall therefore set the gradient coupling to zero, $G=0$, which implies $h_{T}=h$, for all numerical results. This leaves us with the 3 coupling constants $\lambda_{1}$, $\lambda_{2}$, $h$, plus 4 thermodynamic parameters. If we take the condition $H_{c2}=H_{c}$ as an indication for the location of the type-I/type-II transition, then Eq. (24) shows that the transition is, for $G=0$ and fixed $q$, given by a surface in the $\lambda_{1}$-$\lambda_{2}$-$h$-space. (This surface is independent of $\mu_{1}$, $\mu_{2}$, and $T$, but these parameters of course determine the favored phase, and thus, if embedded in the larger parameter space, not everywhere on that surface the COE phase is the preferred phase at $H=0$.) Therefore, the phase diagrams in Fig. 1, where $\lambda_{1}$, $\lambda_{2}$, and $h$ are fixed, are not very useful for our present purpose, and it is more suitable to start from the $\lambda_{1}$-$\lambda_{2}$ plane, where, for a given cross-coupling $h$, we obtain a nontrivial curve $H=H_{c2}$. Two phase diagrams in the $\lambda_{1}$-$\lambda_{2}$ plane at vanishing magnetic field are shown in the upper panels of Fig. 2, one for positive and one for negative cross-coupling $h$. We have chosen the chemical potentials to be larger than the common mass parameter, $\mu_{i}>m$, in which case it is always possible to find negative and positive values of $h$ such that at sufficiently low $T$ and $H$ there is a region in the phase diagram where the COE phase is preferred, cf. Fig. 1.

In the interior of a neutron star, as we move towards the center and thus increase the total baryon number, the system will take some complicated path in our multi-dimensional parameter space, under the assumption that the model describes dense nuclear matter reasonably well. Here we do not attempt to construct this path. But, we keep in mind that nuclear matter is expected to cross the critical surface $H=H_{c2}$ if we move to sufficiently large densities. Therefore, we now choose a path with this property. Starting from the diagrams in Fig. 2, the simplest way to do this is to choose a path in the $\lambda_{1}$-$\lambda_{2}$ plane with all other parameters held fixed. We parametrize the path by $\alpha\in[0,1]$, which is defined by

[TABLE]

with $\vec{\lambda}=(\lambda_{1},\lambda_{2})$. In Fig.  2 we show the paths for positive and negative $h$ that we shall use in the following. Both paths cross from a type-II region for small $\alpha$ into a type-I region for large $\alpha$. In a very crude way, $\alpha$ plays the role of the baryon density in a neutron star. Since our paths are chosen such that $\lambda_{1}$ decreases along them and the charge $q$ is fixed, the Ginzburg-Landau parameter $\kappa$ decreases as $\alpha$ increases.

IV.2 Phases at nonzero temperatures and magnetic fields

In the lower panels of Fig. 2 we show the zero-temperature critical magnetic fields $H_{c}$, $H_{c2}$, and $H_{c1}$, computed as explained in Secs. III.1 – III.3, and the critical temperatures at zero magnetic field, computed from Eqs. (16) for the transition between the COE phase and a single-condensate phase, and with the help of the condensates (9) and (10) together with the thermal masses (14) for the transitions from a single-condensate phase to the NOR phase. The horizontal axis is given by $\alpha$, i.e., we move through the $\lambda_{1}$-$\lambda_{2}$ plane along the paths shown in the upper panels of the figure. In principle, we can use the model straightforwardly to determine the phases in the entire $\alpha$-$H$-$T$-space. As a rough guide to this three-dimensional space notice that increasing the magnetic field at fixed $T$ will eventually destroy the charged condensate, i.e., if $H$ is sufficiently large only the SF and NOR phases survive, while increasing the temperature at fixed $H$ will eventually destroy all condensates, i.e., at sufficiently large $T$ only the NOR phase survives. Working out the details of the entire phase space might be interesting, but it is tedious and not necessary for the main purpose of this paper. Nevertheless, we emphasize that this possibility makes our model very useful for nuclear matter inside a neutron star. For instance, comparing our Fig. 2 with Fig. 1 in Ref. Glampedakis et al. (2011), we see that our results are – on the one hand – a toy version of more concrete calculations of dense nuclear matter, but – on the other hand – more sophisticated because they include all possible phases in a consistent way, not relying on any result within a single-fluid system.

Here we proceed with the discussion of the critical magnetic fields, and for the remainder of the paper we shall restrict ourselves to zero temperature.

IV.3 Type-I/type-II transition region

At first sight, the phase structure in Fig. 2 regarding the critical magnetic fields looks as expected from a single superconductor, only with a critical $\kappa$ that is shifted from the standard value. But, we already know from Sec. III.4 that the point at which the second-order onset of flux tubes turns into a first-order transition is different from the point where $H_{c}$ and $H_{c2}$ intersect. We have marked this point in both lower panels of Fig. 2. Moreover, in the presence of the superfluid, the three critical magnetic fields do not intersect in a single point. This is only visible on a smaller scale, and we discuss this transition region in detail now. With respect to that region, there is no qualitative difference between the two parameter sets chosen in Fig. 2, and therefore we will restrict ourselves to the set with $h<0$.

In Fig. 3, we present the critical magnetic fields in the region that covers their intersection point(s). In the left panel, we have, for comparison, set the coupling to the superfluid to zero, $h=0$, with all other parameters held fixed. As a result, we obtain the expected phase structure of an ordinary superconductor. All three critical magnetic fields intersect at one point – which can be viewed as a check for our numerical calculation of $H_{c1}$ – and this point corresponds to $\kappa^{2}=1/2$. For magnetic fields smaller than $H_{c}$ and $H_{c1}$ the superconductor expels the magnetic field completely, and magnetic fields larger than $H_{c}$ and $H_{c2}$ penetrate the system and superconductivity breaks down. In the open “wedge” between $H_{c1}$ and $H_{c2}$, an array of flux tubes (with varying flux tube density) is expected to exist, with second-order phase transitions at $H_{c1}$ and $H_{c2}$.

In the right panel we zoom in to the critical region of the lower left panel of Fig. 2. From our analytical results we know the following. $(i)$ The critical magnetic fields $H_{c}$, $H_{c2}$ intersect at a point given by Eq. (24), which corresponds for the chosen parameters to $\alpha\simeq 0.37182$. $(ii)$ Just below the curve $H_{c2}$ the flux tube phase is energetically favored over the normal-conducting phase (not necessarily over the Meissner phase) for all $\alpha<0.38265$, as we can compute from Eq. (25). This point is beyond the right end of the scale shown in Fig. 3. $(iii)$ The second-order phase transition from the Meissner phase to the flux tube phase turns into a first-order transition at the point given by Eq. (51), here $\alpha\simeq 0.29236$, which is beyond the left end of the scale of the plot. In the single superconductor, these three $\alpha$’s (or $\kappa$’s) coincide. Had we only computed $H_{c}$, $H_{c2}$, and $H_{c1}$, we would have obtained a puzzling collection of potential phase transition lines. However, together with the first-order phase transitions $H_{c1}^{\prime}$ and $H_{c2}^{\prime}$, computed from Eq. (52), a consistent picture of the phase structure emerges. Before we comment on this structure, we make the behavior at $H_{c1}^{\prime}$ more explicit by plotting the flux tube density $\nu$ and the Gibbs free energies in Fig. 4. The right panel of this figure includes the results for higher winding numbers. We see that they are energetically disfavored for the parameter set chosen here. In Ref. Alford and Good (2008) it was shown that higher winding numbers become important if the magnetic flux, instead of the external field $H$, is fixed. We did check that our numerical results indeed reproduce that observation, but we have not checked systematically whether and for which parameters flux tubes with higher winding numbers are favored in an externally given magnetic field $H$. This is an interesting question for future studies.

The most straightforward interpretation of the right panel of Fig. 3 is to simply ignore the second-order phase transition curves. Then, the topology of the critical region is the same as in the left panel, only with first-order instead of second-order transitions at the boundaries of the flux tube phase (with $H_{c1}^{\prime}$ turning into a second-order phase transition at $\alpha\simeq 0.29236$). However, this cannot be the complete picture. The reason is that after we have left the flux tube phase through $H_{c2}^{\prime}$ and keep increasing $H$ we reach $H_{c2}$, and we know that there should be flux tubes just below $H_{c2}$ for all $\alpha<0.38265$. In other words, our result contradicts the observation that $H_{c2}$ is a lower bound for the transition from the flux tube phase to the normal-conducting phase, as explained at the end of Sec. III.2. This contradiction is resolved when we remember the regime of validity of our approximation for the free energy of the flux tube lattice. Our approximation is accurate where $H_{c1}$ turns into $H_{c1}^{\prime}$ because the distance between the flux tubes is infinitely large at this critical point. As we move along $H_{c1}^{\prime}$ upon increasing $\alpha$, and then along $H_{c2}^{\prime}$ upon decreasing $\alpha$, our approximation becomes worse and worse. Within the present calculation we can thus not determine the phase structure unambiguously, but it is easy to guess a simple topology of the type-I/type-II transition region that is consistent with all our results and takes into account the shortcomings of our approximation. This conjectured phase structure is shown in Fig. 5.

The motivation for the conjecture is as follows. The existence of the first-order line $H^{\prime}_{c1}$ and its starting point is predicted rigorously in our approach. Let us move along that line assuming that we go beyond our approximation and know the complete result. As we move towards large $\alpha$, we will deviate from the line predicted by our approximation. At some value of $\alpha$, we will intersect the curve $H_{c}$. In order to resolve the contradiction of our phase structure, we expect this intersection to occur “on the other side” of the intersection between $H_{c2}$ and $H_{c}$. This implies that our approximation underestimates the binding energy of the flux tubes, i.e., we expect the flux tube phase to be more favored in the full result. We have not found a simple reason – other than the inconsistency of the phase structure – why our approximation distorts the full result in this, and not the other, direction. Now, at the new, correct, intersection of $H^{\prime}_{c1}$ and $H_{c}$, there must necessarily be a third line attached, namely $H_{c2}^{\prime}$ (just like in our approximation). The reason is that if we cross $H_{c1}^{\prime}$ we end up in the flux tube phase and if we cross $H_{c}$ we end up in the normal-conducting phase, and these two phases must be separated by a phase transition line. This critical field $H_{c2}^{\prime}$ might be larger than $H_{c2}$ for all $\alpha$ (below the $\alpha$ of the triple point where $H^{\prime}_{c1}$, $H^{\prime}_{c2}$ and $H_{c}$ intersect) or $H_{c2}^{\prime}$ might merge with $H_{c2}$, leading to an additional critical point. The latter is the scenario shown in the right panel of Fig. 5. One might ask whether $H_{c1}^{\prime}$ and $H_{c}$ intersect exactly at the point where $H_{c}$ and the second-order line $H_{c2}$ intersect. In this case, the entire upper critical line would be of second order and given by $H_{c2}$. However, this seems to require some fine-tuning of the interaction between the flux tubes since the second-order line $H_{c2}$ does not know anything about this interaction.

IV.4 Flux tube clusters

The first-order phase transitions with $H$ as an external variable translate into mixed phases if we fix the magnetic field $B$ (spatially averaged) instead. Again, this can be illustrated by the analogy to the onset of baryonic matter at small temperatures. As a function of $\mu_{B}$, this onset is a first-order transition with a discontinuity in baryon number density $n_{B}$. If we instead probe this onset with fixed $n_{B}$ (spatially averaged), we pass through a region of mixed phases, for example nuclei in a periodic lattice, until we reach the saturation density. These mixed phases are realized in the outer regions of a neutron star, and it would be an intriguing manifestation of this analogy if the mixed flux tube phases discussed here are realized in the core of the star. Each first-order transition in $H$ yields two critical magnetic fields $B$ which we compute as follows. At $H_{c1}^{\prime}$, the lower critical field is $B=0$, and the upper critical field is $\langle B\rangle=\Phi_{0}\nu$ [using Eq. (29)], where $\nu$ is the numerically computed flux tube area density as we approach the first-order transition from above; at $H_{c2}^{\prime}$, the lower critical field is $\langle B\rangle=\Phi_{0}\nu$, with $\nu$ now being the numerically computed density as we approach the first-order transition from below, while the upper critical field is $B=H_{c2}$; at $H_{c}$, the lower critical field is $B=0$, and the upper one is $B=H_{c}$. We perform this calculation with the parameters of Fig. 3. As discussed for the $H$-$\alpha$ phase diagrams above, also for the $B$-$\alpha$ phase structure we do not expect our approximation to yield quantitatively reliable results where the flux tube density is large. Therefore, our results reflect the topology of the $B$-$\alpha$ phase diagram correctly, but the precise location of the phase transition lines cannot be determined within our approach. The phase diagrams for the single superconductor and the two-component system are shown in Fig. 6. In a single superconductor, there is only one possible mixed phase: macroscopic regions in which the magnetic field penetrates, mixed with regions in which the magnetic field remains expelled Tinkham (2004). The geometric structure of these regions depends on the details of the system such as the surface tension, and it is beyond the scope of this paper to determine them. In the two-component system, two additional mixed phases are possible, both of which contain flux tube clusters. (Unrelated to the first-order phase transitions pointed out here, flux tube clusters have been suggested to exist in neutron stars in the vicinity of superfluid neutron vortices Sedrakian and Sedrakian (1995).) Firstly, at $H_{c1}^{\prime}$, flux tube clusters are immersed in a field-free superconducting region, as predicted for “type-1.5 superconductivity” Babaev and Speight (2005). Secondly, at $H_{c2}^{\prime}$, there is a mixed phase of flux tubes with the normal-conducting phase, i.e., superconducting regions that enclose flux tubes and that are themselves surrounded by completely normal-conducting regions.

V Conclusions

We have shown that the coupling to a superfluid can have profound effects on the magnetic properties of a superconductor. We have started from a microscopic model for two complex scalar fields, coupled to each other via density and gradient coupling terms, with one of the fields being electrically charged. By computing the thermal excitations of the system we have derived a Ginzburg-Landau-like effective potential for the charged and neutral condensates and the gauge field. This potential has then been evaluated at nonzero temperatures and external magnetic fields, computing the two condensates dynamically for all 4 possible phases: condensation of both fields (superconductor + superfluid), condensation of only one field (pure superconductor or pure superfluid), or no condensation. We have discussed the structure of the resulting phase diagram in the multi-dimensional parameter space, with the main focus on the transition region between type-I and type-II superconductivity. To this end, we have computed the critical magnetic fields $H_{c}$, $H_{c2}$ (analytically) and $H_{c1}$ (numerically, based on the profile functions of a magnetic flux tube). In contrast to the standard scenario of a single superconductor, these three magnetic fields do not intersect in a single point if the superconductor coexists with a superfluid. The phase structure around these intersection points is (at least partially) resolved by computing the first-order phase transitions $H_{c2}^{\prime}$ and $H_{c1}^{\prime}$. This has been done by employing a simple approximation for the free energy of a flux tube array that is valid for large flux tube distances and that effectively reduces the calculation to solving the equations of motion for a single flux tube. The new critical fields $H_{c}$, $H_{c2}^{\prime}$, $H_{c1}^{\prime}$ do intersect in a single point, restoring the topology of the transition region, with (segments of) the second-order transition lines replaced by first-order transitions. In particular, we have identified a new critical point – and derived an analytical expression for its location – where the second-order flux tube onset $H_{c1}$ turns into a first order transition $H_{c1}^{\prime}$. The presence of the first-order transitions allows for mixed phases with flux tube clusters, very similar to a type-1.5 superconductor, which consists of two charged fields coupled indirectly through the gauge field.

There are several possible improvements and extensions of our work. Our approximation for the flux tube array can be improved for instance by determining dynamically the values of the condensates far away from the flux tubes instead of using the values of the homogeneous phase. To settle the precise location of the phase transition lines, it would be interesting to perform a brute force numerical calculation of the free energy of the flux tube phase, for which our results are a valuable guidance. There are several other interesting aspects of our model which we have mentioned but not worked out in detail. For instance, one could perform a more systematic study of the effect of the derivative coupling, which we have included in all our analytical results, but set to zero in the final numerical results of the phase diagrams. Or one could perform a more detailed study of flux tubes with higher winding numbers, which turned out to be energetically disfavored for the parameter regime we have studied, but which are known to potentially play a role in the two-component system. One can also study the phase structure at nonzero temperature in more detail and/or improve the large-temperature approximation on which our Ginzburg-Landau potential was based. Or one can include superfluid vortices, aiming at the phase structure at nonzero magnetic field and externally imposed rotation.

Our setup and our results are applicable to dense nuclear matter in the core of neutron stars. For instance, one can fit our model parameters, such as the density coupling and gradient coupling, to values predicted for nuclear matter and eventually compute the phase structure as a function of the baryon number density rather than of an abstract model parameter. One may also ask whether a potential phase of flux tube clusters would affect the transport properties of the core in a detectable way. Moreover, it would be interesting to employ our results in studies of the time evolution of the magnetic field in a neutron star. Here we have computed the ground state in equilibrium for given temperature, magnetic field and chemical potential, but for more phenomenological predictions one needs to know whether and on which time scale this ground state is reached.

Acknowledgements.

We would like to thank Mark Alford, Nils Andersson, Egor Babaev, Christian Ecker, Carlos Lobo, David Müller, and Andreas Windisch for valuable comments and discussions. We acknowledge support from the Austrian Science Fund (FWF) under project no. W1252, and from the NewCompStar network, COST Action MP1304. A.S. is supported by the Science & Technology Facilities Council (STFC) in the form of an Ernest Rutherford Fellowship.

Appendix A Derivation of the effective potential

In this appendix we compute an effective potential in a high-temperature approximation from the excitations of the system, taking into account the mixing of the two scalar fields with the photon. Elements of this derivation can be found in discussions of the standard abelian Higgs model, see for instance chapter 85 of Ref. Srednicki (2007).

It is convenient to split the complex scalar fields into their real and imaginary parts,

[TABLE]

Then, the Lagrangian (1), in the presence of chemical potentials $\mu_{1}$ and $\mu_{2}$, becomes

[TABLE]

where we have added a gauge fixing term,

[TABLE]

and where

[TABLE]

We allow for condensation of both fields by shifting $\phi_{1}\to\rho_{1}+\phi_{1}$, $\phi_{2}\to\rho_{2}+\phi_{2}$, i.e., we assume the condensates to be real, and from now on $\phi_{i}$ and $\chi_{i}$ are fluctuations about the condensates. The dispersion relations of the excitations are computed from the tree-level propagator in momentum space. To this end, we introduce the Fourier transformed fields via

[TABLE]

with the space-time four-vector $X=(-i\tau,\vec{r})$ and the four-momentum $K=(k_{0},\vec{k})$, where $k_{0}=-i\omega_{n}$ with the bosonic Matsubara frequencies $\omega_{n}=2\pi nT$, $n\in\mathbb{Z}$. In the imaginary time formalism, we have to replace $A_{0}\to iA_{0}$. The terms of second order in the fluctuations can then be written as

[TABLE]

with

[TABLE]

The inverse tree-level propagator is an $8\times 8$ matrix, which reads

[TABLE]

with the scalar field sector,

[TABLE]

where $\eta_{1/2}(\rho_{1},\rho_{2})\equiv-(\mu_{1/2}^{2}-m_{1/2}^{2})+\lambda_{1/2}\rho_{1/2}^{2}-h\rho_{2/1}^{2}$, the inverse gauge field propagator,

[TABLE]

where $\sigma\equiv 1-1/\xi$, and the off-diagonal blocks that couple the scalar fields to the gauge field,

[TABLE]

We are interested in an effective potential for the condensates $\rho_{1}$ and $\rho_{2}$, and thus we need to keep these condensates general. Nevertheless, it is instructive to first discuss the dispersions at the zero-temperature stationary point, i.e., we set $\rho_{1}=\rho_{01}$ and $\rho_{2}=\rho_{02}$ with the condensates in the coexistence phase $\rho_{01}$ and $\rho_{02}$ from Eq. (11). Let us first set the cross-coupling between the scalar field to zero, $h=G=0$. The dispersion relations $k_{0}=\epsilon_{k}$ are given by the zeros of ${\rm det}\,S^{-1}$. Since this is a polynomial of degree 8 in $k_{0}^{2}$, we obtain 8 dispersions, 6 of which are physical. The two unphysical ones are of the form $\epsilon_{k}=k$. These are the usual unphysical modes of the gauge field, whose contribution to the partition function is canceled by ghost fields. With the given gauge choice, ghosts do not couple to any of the fields and merely serve to cancel the unphysical modes. None of the modes depend on the gauge fixing parameter $\xi$, which only appears as a prefactor of the determinant ${\rm det}\,S^{-1}$ and thus does not have to be specified. The 6 physical dispersions are

[TABLE]

We have three gauge field modes with mass $\epsilon_{k=0}=\sqrt{4\pi}q\rho_{01}$ [the two modes of Eq. (75a) and the mode with the lower sign in Eq. (75b)], two more massive modes from the scalar fields, and the Goldstone mode [the mode with the lower sign in Eq. (75c)]. Let us now switch on the couplings $G$ and $h$ between the scalar fields. We find that the unphysical modes remain unaffected and all modes remain independent of the gauge fixing parameter $\xi$. The mass of the three gauge field modes and the entire dispersion (75a) for two of them is also unchanged. The expressions for the remaining dispersions become very complicated. They can easily be computed numerically, and we show the result in Fig. 7.

We now compute our effective potential by reinstating the general condensates, i.e., we need to compute the dispersions away from stationary point $\rho_{1}=\rho_{01}$, $\rho_{2}=\rho_{02}$. We restrict ourselves to the following large-momentum approximation, which is sufficient for the high-temperature approximation we are interested in,

[TABLE]

such that

[TABLE]

In general, the dispersions now do depend on the gauge fixing parameter $\xi$. However, in the limit (76) this dependence drops out, i.e., the coefficients $c_{1}$ and $c_{2}$ do not depend on $\xi$. Moreover, now the unphysical gauge modes no longer have the simple form $\epsilon_{k}=k$. Two of the physical gauge modes keep their simple form (75a), while for the other 4 physical modes the coefficients $c_{1}$ and $c_{2}$ are (at least some of them) very lengthy. However, adding up the result for all 6 physical modes yields a relatively compact result,

[TABLE]

where, in the second step, we have absorbed all terms that do not depend on $\rho_{1}$ or $\rho_{2}$ into “const.”, and dropped all higher-order terms in the derivative coupling (i.e., we assume $G\mu^{2}\ll 1$, where $\mu$ stands for all energy scales $\mu_{1}$, $\mu_{2}$, $m_{1}$, $m_{2}$, $\rho_{1}$ $\rho_{2}$). Dropping the constant contribution, we add the $T^{2}$ terms to the potential (5) and arrive at the potential (13) in the main text.

Appendix B Calculation of $H_{c2}$ and Gibbs free energy just below $H_{c2}$

Here we derive Eqs. (24) and (25). To this end, we need the equations of motion for the scalar fields and the gauge field. We go back to the Lagrangian (1), take the static limit and replace the parameters $m_{i}$ and $h$ by their $T$-dependent generalizations $m_{i,T}$ and $h_{T}$. This yields the potential

[TABLE]

and the equations of motion for $\varphi_{1}^{*}$, $\varphi_{2}^{*}$, and $\vec{A}$ become

[TABLE]

Since the transition from the flux tube phase to the normal-conducting phase is assumed to be of second order, the charged condensate becomes infinitesimally small just below $H_{c2}$, and we make the ansatz $\varphi_{1}=\bar{\varphi}_{1}+\delta\varphi_{1}$ with $\bar{\varphi}_{1}\propto(H_{c2}-H)^{1/2}$, and $\delta\varphi_{1}$ includes terms of order $(H_{c2}-H)^{3/2}$ and higher, i.e., is at least of order $\bar{\varphi}_{1}^{3}$. We also introduce perturbations for the neutral condensate and the gauge field, $\varphi_{2}=\bar{\varphi}_{2}+\delta\varphi_{2}$, $\vec{A}=(\bar{A}_{y}+\delta A_{y})\vec{e}_{y}$, where $\delta A_{y},\delta\varphi_{2}$ include terms of order $\propto H_{c2}-H$ and higher, i.e., they are at least of order $\bar{\varphi}_{1}^{2}$. As the magnetic field completely penetrates the superconductor at the phase transition, we can choose the unperturbed gauge field to be of the form $\bar{A}_{y}=xH_{c2}$, and we denote $\delta B=\partial_{x}\delta A_{y}$, such that $\vec{B}=(H_{c2}+\delta B)\vec{e}_{z}$. We assume all functions to be real and to depend on $x$ only, not on $y$ and $z$ (solutions with these properties are sufficient for our purpose, the derivation would also work without these restrictions but would be somewhat more tedious). We insert this ansatz into the equations of motion (80), and keep terms up to order $\bar{\varphi}_{1}^{3}$. Then, the linear contributions from Eqs. (80a) and (80b) yield two equations for $\bar{\varphi}_{1}$ and $\bar{\varphi}_{2}$,

[TABLE]

with

[TABLE]

while the subleading contributions from Eqs. (80a) and (80b) and the leading contribution from Eq. (80c) yield the following equations for the perturbations $\delta\varphi_{1}$, $\delta\varphi_{2}$, and $\delta A_{y}$,

[TABLE]

Inserting our ansatz into the potential (79), using partial integration and the equations of motion (81) and (83), and keeping terms up to order $\bar{\varphi}_{1}^{4}$, we find after some algebra the free energy

[TABLE]

We will first compute $H_{c2}$ from Eqs. (81) and afterwards compute the Gibbs free energy just below $H_{c2}$ from Eq. (84).

We assume the neutral condensate in the SF phase to be homogeneous, and thus Eq. (81b) yields $2\bar{\varphi}_{2}^{2}=\rho_{\rm SF}^{2}$, as expected. For the solution of Eq. (81a) we can simply follow the textbook arguments because it has the same structure as for a single-component superconductor. It reads

[TABLE]

and thus is equivalent to the Schrödinger equation for the one-dimensional harmonic oscillator, $-\frac{\hbar^{2}}{2m}\psi^{\prime\prime}(x)+\frac{m}{2}\omega^{2}x^{2}\psi=E\psi$ with the identification $E/(\hbar\omega)=(\lambda_{1}\rho_{\rm SC}^{2}+h_{T}\rho_{\rm SF}^{2})/(2qH_{c2})$. Since the eigenvalues are $E_{n}=(n+\frac{1}{2})\hbar\omega$, the largest magnetic field for which the equation allows a physical solution is obtained by setting $n=0$,

[TABLE]

in agreement with Eq. (13) of Ref. Sinha and Sedrakian (2015). In the second expression we have rewritten the condensates $\rho_{\rm SC}$ and $\rho_{\rm SF}$ in terms of the charged condensate in the coexistence phase $\rho_{01}$, see Eq. (11), and used the definition of the coherence length $\xi$ from Eq. (23). Since the relevant eigenvalue of Eq. (85) is given by $n=0$, the corresponding eigenfunction is a Gaussian,

[TABLE]

where the exact value of the prefactor $C_{0}\propto(H_{c2}-H)^{1/2}$ is not relevant for the following. The result shows that, for $H$ just below $H_{c2}$, charged condensation with small magnitude of order $(H_{c2}-H)^{1/2}$ occurs in a slab confined in a direction perpendicular to the external magnetic field, here chosen to be the $x$-direction, with width $(qH_{c2})^{-1/2}$. Had we allowed for $y$ and $z$ dependencies of the condensate, we could have used this linearized approximation to discuss crystalline configurations and determine the preferred lattice structure. Here we continue by checking whether the solution (87) is energetically preferred over the normal-conducting phase for $H$ below and close to $H_{c2}$. To this end, we need to compute the Gibbs free energy, as defined in Eq. (18), from the free energy (84). We first solve Eq. (83c) with the boundary condition $\delta B(\pm\infty)=H-H_{c2}$ (since $B=H$ in the normal-conducting phase) to find

[TABLE]

Inserting this result into Eq. (84) and using Eq. (87) yields the Gibbs free energy

[TABLE]

with $\mathcal{G}_{\rm SF}$ from Eq. (20). It remains to compute $\delta\varphi_{2}$. We use Eq. (83b), which can be written as

[TABLE]

with the dimensionless variable $t=\sqrt{qH_{c2}}\,x$ and the dimensionless quantities

[TABLE]

where $p$ indicates the magnitude of the neutral condensate and $\gamma$ the magnitude of the gradient coupling $G$ relative to the density coupling $h_{T}$, both in units given by the critical magnetic field. With the boundary conditions $\delta\varphi_{2}(\pm\infty)=0$, this equation has the solution

[TABLE]

where we have abbreviated

[TABLE]

with the error function erf. Inserting Eq. (92) into Eq. (89) yields

[TABLE]

where $\langle\ldots\rangle$ denotes spatial average, and

[TABLE]

We discuss this result for the case without gradient coupling, $\gamma=0$, in the main text.

Appendix C Interaction between two flux tubes

In this appendix we derive the expression for the interaction energy Eq. (48). We start from the definition (47), i.e., we consider two parallel flux tubes $(a)$ and $(b)$ separated by the (dimensionless) distance $R_{0}$. We divide the total volume $V$ into two half-spaces $V^{(a)}$ and $V^{(b)}$, which are the simplest versions of two Wigner-Seitz cells: we connect the two flux tubes by a line with length $R_{0}$, and the plane in the center of and perpendicular to that line divides $V$ into $V^{(a)}$ and $V^{(b)}$. The interaction free energy is then computed from

[TABLE]

where, due to the symmetry of the configuration, we have restricted the integration to the half-space $V^{(a)}$, where $U^{(a)}_{\circlearrowleft}$, $U^{(b)}_{\circlearrowleft}$ are the free energy densities of the two flux tubes in the absence of the other flux tube, and where $U_{\circlearrowleft}^{(a)+(b)}$ is the total free energy of the flux tubes. (Recall that by definition $U_{\circlearrowleft}$ denotes the pure flux tube energy density, with the free energy density of the homogeneous configuration already subtracted.)

We assume $R_{0}$ to be much larger than the widths of the flux tubes, such that the contribution of flux tube $(b)$ to the free energy is small in $V^{(a)}$. Therefore, we will now compute the free energy density of a “large” contribution that solves the full equations of motion plus a “small” contribution that solves the linearized equations of motion. We shall do so in a general notation, not referring to the geometry of our two-flux tube setup. Only in Eq. (103), when we insert the results into the free energy (96), we shall come back to this setup and introduce a more explicit notation indicating the contributions of the two different flux tubes. Following Ref. Kramer (1971), we define

[TABLE]

and write

[TABLE]

The equations of motion for a single flux tube to leading order, $\delta\vec{Q}=\delta f_{1}=\delta f_{2}=0$, are (from now on, in this appendix, all gradients are taken with respect to the dimensionless coordinates)

[TABLE]

[equivalent to Eqs. (39) in the main text], and the equations of motion of first order in the corrections $\delta\vec{Q}$, $\delta f_{1}$, $\delta f_{2}$ become

[TABLE]

We denote the free energy density, up to second order and after using the equations of motions, by $U_{0}+\delta U$, where

[TABLE]

is the free energy density of a single flux tube from Eq. (38), and the first-order and second-order corrections can be written as a total derivative,

[TABLE]

Notice that any explicit dependence on the density coupling $h_{T}$ has disappeared, while the derivative coupling $\Gamma$ does appear explicitly.

We can now go back to the interaction free energy (96) and identify the full free energy $U_{\circlearrowleft}^{(a)+(b)}$ in the half-space $V^{(a)}$with $U_{0}+\delta U$. In $V^{(a)}$, $U_{\circlearrowleft}^{(a)}$ is given by setting $\delta\vec{Q}=\delta f_{1}=\delta f_{2}=0$ in $U_{0}+\delta U$ (which simply leaves $U_{0}$), and $U^{(b)}_{\circlearrowleft}$ is obtained by setting $\vec{Q}_{0}=0$, $f_{10}=f_{20}=1$ in $U_{0}+\delta U$ (which leaves various terms from $\delta U$). Consequently, we find

[TABLE]

where we have rewritten the volume integral as a surface integral and where we have made the contributions from the two flux tubes $(a)$ and $(b)$ explicit. Since the derivatives of all fields vanish at infinity, the integration surface is reduced to the plane that separates the two Wigner-Seitz cells. We now use the geometry of the setup to simplify this expression: we align the $z$-axis with flux tube $(a)$, such that this flux tube sits in the origin of the $x$-$y$ plane, with the $x$-axis connecting the two flux tubes. Therefore, $\vec{Q}^{(a)}$, $f_{10}^{(a)}$, $f_{20}^{(a)}$ are functions only of $R$, while $\delta\vec{Q}^{(b)}$, $\delta f_{1}^{(b)}$, $\delta f_{2}^{(b)}$ also depend on the azimuthal angle $\theta$. However, since we only need the functions and their gradients at the boundary between the two Wigner-Seitz cells and since this boundary is by assumption far away not only from flux tube $(b)$ but also from flux tube $(a)$, we can write ($i=1,2$)

[TABLE]

Note in particular that, at the relevant surface, $\delta f_{i}^{(b)}=\delta f_{i}^{(a)}$, but $d\vec{S}\cdot\nabla\delta f_{i}^{(a)}=-d\vec{S}\cdot\nabla\delta f_{i}^{(b)}$. Now, $\delta Q$ and $\delta f_{i}$ are functions only of $R$. Inserting Eqs. (104) into Eq. (103) yields

[TABLE]

We can employ this result by inserting the modified Bessel functions from Eq. (45),

[TABLE]

We may also extrapolate this result down to smaller distances by reinstating the full numerical functions through $\delta Q\to Q=-n(1-a)/R$ and $\delta f_{i}\to 1-f_{i}$, which yields the result (48) in the main text.

Appendix D Asymptotic approximation of flux tube interaction with gradient coupling

In the main text, we discuss the large-distance behavior of the flux tube interaction without gradient coupling. In the presence of a gradient coupling, the interaction is more complicated, but, as we show in this appendix, an equally compact expression can be derived if we are only interested in the leading order contributions, i.e., the exponential behavior.

We start by inserting the asymptotic solutions (45) into the expression for the interaction free energy (48). The result is an integral over a sum of many terms, each of which is a product of 2, 3, or 4 modified Bessel functions of the second kind. In each product, one factor is $K_{1}$ and the remaining ones are $K_{0}$. The integral over the terms with 2 Bessel functions that have the same argument can be expressed again as a Bessel function with the help of Eq. (50). For the integral over all other products we use the expansion,

[TABLE]

and only keep terms with the smallest exponential suppression. These terms are found as follows. With Eq.  (107) we approximate

[TABLE]

Then, we approximate each product of Bessel functions $K_{0}K_{1}$, $K_{0}K_{0}K_{1}$, $K_{0}K_{0}K_{0}K_{1}$ by the leading order term, and re-express the exponential as a product $K_{0}K_{1}$ with the help of Eq. (108). If we have started with a product $K_{0}K_{1}$ with different arguments, we arrive at an expression which we can integrate using Eq. (50). If we have started with a product of 3 or 4 Bessel function, we do not exactly reproduce the integrand of Eq. (50) because there is an additional factor $R^{-1/2}$ (for 3 Bessel functions) or $R^{-1}$ (for 4 Bessel functions). The resulting integral can be expressed in terms of the so-called Meijer G-function, which we expand again since we are anyway only interested in the asymptotic behavior. As a result, we obtain

[TABLE]

For each of the terms in the interaction energy we need to replace $\alpha_{i}$ by either $\sqrt{\nu_{+}}$ or $\sqrt{\nu_{-}}$. From Eq. (44) we see that $\sqrt{\nu_{+}}>\sqrt{\nu_{-}}$. Therefore, the largest contribution we obtain is $\exp(-\sqrt{\nu_{-}}R_{0})$, and this contribution is only created by the product of 2 Bessel functions with the same argument $\sqrt{\nu_{-}}$ because 2 Bessel functions with different arguments give rise to $\exp[-(\sqrt{\nu_{+}}+\sqrt{\nu_{-}})R_{0}/2]$, which is suppressed more strongly, 3 Bessel functions give rise to suppressions of at least $\exp[-3\sqrt{\nu_{-}}R_{0}/2]$ etc. The largest contributions are thus given by the terms where we can apply the integral (50), and we obtain

[TABLE]

Therefore, if $\gamma_{-}^{2}+x^{2}-\Gamma x\gamma_{-}>0$, one can use the same arguments as in the main text for the discussion of the attractiveness of the flux tube interaction at large distances, only with a more complicated eigenvalue $\nu_{-}$, which now depends on the gradient coupling $\Gamma$.