# Topological defects in mixtures of superconducting condensates with different charges

## Abstract

We investigate the topological defects in phenomenological models describing mixtures of charged condensates with commensurate electric charges. Such situations are expected to appear for example in liquid metallic deuterium. This is modeled by a multicomponent Ginzburg-Landau theory where the condensates are coupled to the same gauge field by different coupling constants whose ratio is a rational number. We also briefly discuss the case where electric charges are incommensurate. Flux quantization and finiteness of the energy per unit length dictate that the different condensates have different winding and thus different number of (fractional) vortices. Competing attractive and repulsive interactions lead to molecule-like bound state between fractional vortices. Such bound states have finite energy and carry integer flux quanta. These can be characterized by topological invariant that motivates their denomination as skyrmions.

###### pacs:

67.85.Jk,74.25.Ha, 67.85.Fg## Introduction

Although recently there has been substantial interest in multicomponent superconductors, typically the research is restricted to fields which have same value of electric charge modulus. See, e.g., Ref. Herland et al., 2010 for recent work with a field overview. In the typical condensed matter systems the charge is set by a Cooper pair charge of for electronic systems or for protonic superconductors Babaev et al. (2004); Babaev and Ashcroft (2007). By contrast, multicomponent systems with different values of electric charge attracted much less attention, yet some were discussed in the literature. One example is liquid metallic deuterium where deuteron is a charge- boson which can Bose condense and coexist with the Cooper pairs of electrons and/or protons Oliva and Ashcroft (1984a, b); Babaev et al. (2004); Bedaque et al. (2011). Cooper pairs carry twice the charge of their constituent fermion, while a Bose-Einstein condensate of deuterons carries only once the charge of its (boson) constituent. This system is currently a subject of vigorous experimental pursuit Eremets and Troyan (2011).

Mixtures of condensates carrying different electric charges may also apply to ultra-cold atomic gases with synthetic gauge field. Recent progress both in theoretical understanding and experimental techniques to control such systems makes it promising that artificial dynamical gauge fields may be realized there Banerjee et al. (2012); Zohar et al. (2012). In that case, effectively a system could be described by multicomponent charged condensate models. Our present discussion may, in the future, find applications to these systems.

In superconductors with more than two components and broken symmetry there can be pairing transitions to paired states driven by proliferation of composite vortices Babaev (2004); Kuklov and Svistunov (2003); Smørgrav et al. (2005). Such pairing mechanisms can lead to charge- electronic superconducting systems as hypothesized recently in various contexts Agterberg and Tsunetsugu (2008); Berg et al. (2009); Herland et al. (2010), along with other recently discussed mechanisms for charge- superconductivity Moon (2012).

The above examples from superconductivity and superfluidity, along with the multicomponent gauge theories which appear as effective field theories in other condensed matter systems Sachdev (2008); Chen et al. (2013), calls for investigation of mixtures of charged condensates with arbitrary ratio of condensate charges. To this end we study a phenomenological Ginzburg-Landau model that accounts for such mixtures, regardless of their underlying microscopic origin. We discuss below that if the mixtures of charged condensates carrying different electric charges are realized, either in natural or artificial systems, their response to external applied magnetic field or rotation would be very different from those of mixtures with equal charges. This is because of the substantial difference in the topological excitations as compared to systems where condensates carry the same electric charge.

In Section I we introduce a mean-field model that accounts for mixtures of charged condensates, when the condensates carry different electric charges. Section II is devoted to elementary topological excitations, that is fractional vortices, and flux quantization in our Ginzburg-Landau model. In Section III, we investigate the physics of flux carrying topological excitations within the London limit where condensates are assumed to have constant densities.

## I Mixtures of charged condensates

A charged condensate can, under certain conditions, be described by mean field Ginzburg-Landau free energy that couples it to the vector potential of the magnetic field through the kinetic term

(1) |

Here is the electric charge of the condensate and its rest mass. and are respectively the reduced Planck constant and the speed of light in the vacuum. is the interacting potential. For example, for an ordinary superconductor, a Cooper pair has charge twice that of an electron. Thu, there . On the other hand, a Bose-Einstein condensate of singly charged bosons will have electric charge . Thus mixtures of charged condensates, should generically have different couplings to the vector potential. A mixture of condensates with different masses and charges is thus described by

(2) |

Note that since the theory is invariant under complex conjugation, it is sufficient to consider positive charge only, without losing generality, condensates with negative charge being obtained by complex conjugation of the one with the positive charge. Now, to get rid of superfluous parameters, we express the energy in units of and rescale the fields as

(3) |

Dropping the symbol further on, a mixture of charged condensates is thus described by the free energy density

(4) |

Here are complex fields that stand for the charged condensates (with different indices ). The condensates are coupled together through the electromagnetic interactions mediated by the vector potential in the kinetic terms . Since are condensates that are essentially different, they should be independently conserved. It results that the potential has global invariance ensuring independent conservation of both particle numbers

(5) |

In the condensed phase, are negative parameters while . The model exhibits gauge invariance under local transformations. That is, for arbitrary , the energy (4) is unchanged under the transformations

(6) |

As will be explained later on, a necessary condition for finite energy flux carrying configurations is that the charge of the condensates should be commensurate. That is, the ratio of the coupling constants is a rational number: . To capture this constraint, it is convenient to parametrize the gauge couplings as where are integer numbers (). Moreover, the two integers and should be relatively prime (their greatest common divisor must be ). Within our parametrization, is an arbitrary number that uniquely parametrizes the London penetration length defined below. If we apply this model to liquid metallic deuterium, the couplings are and denotes the deuteron condensate while (carrying twice the electric charge of ) denotes electronic Cooper pairs.

Functional variation of the free energy (4) determines the Euler-Lagrange equations of motion. That is, variation with respect to complex fields gives the Ginzburg-Landau equation for the charged condensates, while variation with respect to the vector potential defines Ampère’s law

(7) |

with the supercurrent

(8) |

In the ground state, the condensates have constant densities and is a pure gauge. The length scales at which the condensates recover their ground state value after infinitesimal perturbations, the coherence lengths, are . The penetration depth of the magnetic field is consistently derived in section III.

When considering vortex matter, we restrict ourselves to field configurations varying in the plane only and with normal magnetic field, that is, field configurations describing both two-dimensional systems and three-dimensional system invariant under translations along the normal direction. To investigate the physical properties of topological excitations in our model for mixtures of charged condensates, we numerically minimize the free energy (4) within a finite element framework provided by the Freefem++ library Hecht (2012). For technical details, see the discussion in Appendix B.

## Ii Topological defects

Because we consider several condensates, the elementary topological excitations are fractional vortices, that is, field configurations with phase winding of a single condensate (e.g. has winding while ). A fractional vortex carries a fraction of the flux quantum. This can be seen by deriving the quantization condition for the magnetic flux. The supercurrent (8), defined from Ampère’s equation , reads as

(9) |

Here we defined the weighted density . Since the supercurrent is screened, it decays exponentially and the magnetic flux thus reads as

(10) |

Since the condensates are complex fields, their phases wind integer number of time. The couple denotes the field configurations with winding of the condensate . The elementary excitations are fractional vortices and with unit winding in each component. A given fractional vortex in the condensate thus carries a fraction of the magnetic flux . Here is the flux quantum. For the magnetic flux to be quantized, as long as , it is necessary that different condensates have different winding number . This follows from

(11) |

When both condensates carry the same electric charge, , the quantization condition (11) reduces to the quantization conditions for multiband/multicomponent superconductors Babaev (2002).

Thus, each condensate has to wind times, so that the resulting composite vortex carries one flux quantum . Fractional vortices have logarithmically divergent energy per unit length and thus cannot form in bulk systems Babaev (2002). This can be seen by rewriting the free energy into charged and neutral modes. For this, and using (9), the kinetic term can be rewritten

(12) |

with again the weighted density . Defining the weighted phase difference , the free energy (4) reads as

(13a) | ||||

(13b) | ||||

(13c) |

Since it decouples from the gauge field, the term (13c) is called neutral mode. This is the kinetic energy of the relative motion of the two condensates. That is, the co-directed (counter-directed) motion of particles with opposite (alike) charges . When has a winding, this neutral mode has logarithmically divergent energy. Indeed, asymptotically each phase is well approximated by , where are the (integer) vorticities and the polar angle. The condition for the logarithmic divergence to be absent and thus for the energy to be finite thus reads as

(14) |

For a configuration carrying a single flux quantum (and since
and are integers), the absence of winding in the weighted phase
difference dictates that . Thus the configurations which have
no logarithmic divergence winds times in and
times in . Note that this condition also implies that the
total flux is integer (11).
On the other hand, fractional vortices have logarithmically divergent
energy. This is because they have winding in the weighted phase
difference. That is, when , since the screening is
incomplete, the energy of the vortex grows with the system size.
Such vortices are typically thermodynamically unstable in bulk systems.
When both condensates carry the same electric charge , the
condition (II) is automatically satisfied, provided both
condensates have the same winding .
That is, only co-centered composite vortices with same winding in both
condensates have finite energy. Since fractional vortices have logarithmically
divergent energy per unit length, they cannot form in bulk systems
^{1}

Consider now the simplest case where the condensate carries
single charge while carries double : . The
resulting composite vortex carrying one flux quantum is the bound state
of one (fractional) vortex in and two vortices in . To
reduce the cost of kinetic energy in the neutral sector, the interaction
through the neutral sector binds vortices together. It is minimal when
cores in different condensates coincide (see discussion below in
Section III). However the two vortices in quite naturally repel
each other, as would two Abrikosov vortices do in single band systems.
As a result we can expect that if the magnetic repulsion is strong enough,
despite the attractive channel through the neutral sector, fractional
vortices will not overlap.
This is indeed the case, as shown in Fig. 1. The regime
shown on the first line has , so the condensates have
single and double winding respectively. The resulting topological defect,
carrying a single flux quantum, is a bound state of fractional vortices
that are not co-centred and the configuration looks like an elongated rod.
We refer to such bound state as a vortex molecule
^{2}

Bound states of non-overlapping fractional vortices feature special topological properties that motivate their denomination as skyrmions. The terminology follows from the fact that two-component models can be mapped to easy-plane nonlinear -models that are associated with a topological invariant Babaev et al. (2002); Babaev (2009); Garaud et al. (2014). In systems with the same charge for both condensates, the two-component model is rewritten in term of the total current , the total density , and the pseudo-spin . The pseudo-spin unit vector is the projection of the superconducting condensates on spin- Pauli matrices :

(15) |

When all condensates have the same couplings , the spinor is defined as the two-vector of complex condensates . The finiteness of the energy dictates that is asymptotically a constant vector, while vanishing of neutral modes implies that has no winding.

To obtain a similar projection for incommensurate charges, the spinor should be chosen so that the pseudo-spin does not wind asymptotically (when neutral mode vanish). There are several possibilities to realize such a projection and we choose

(16) |

The projection (15) of (16) maps to the two-sphere target space. It determines the pseudo-spin that we use along the paper for the visualization of the pseudo-spin texture. Note that (16) is a non-holomorphic map, so it hard to justify the quantization of the associated invariant. For this we introduce another map which is holomorphic

(17) |

The associated projection is a map from the one-point compactification of the plane () to the two-sphere target space spanned by . That is , which is classified by the homotopy class . This defines the integer valued topological charge

(18) |

If everywhere, is an integer number. Roughly speaking, counts the number of times the texture (15) covers the target sphere. For practical purpose, and since it gives much better accuracy, we compute the degree of the maps (16) or (17), instead of computing the formula (18). Numerically calculated topological charge for the various configurations we constructed is indeed found to be integer (with a negligible error of order ).

The map (17) provides a rigorous justification of the topological invariant (18), but the associated texture turns out to be very difficult to visualize. So we still use (16) that gives more straightforward physical interpretation for the visualization of the pseudo-spin texture. Interestingly, when computed numerically both definitions give similar result and accuracy of the topological charge.

In our parametrization of the gauge couplings , the charges of the condensates are commensurate. That is their ratio is a rational number and the simplest case we discussed is . There, the unit flux quantum excitation is a molecule made of two fractional vortices in bound together by a single vortex in the condensate. For different ratio of the gauge couplings, the unit flux molecule-like bound states assume very rich structures as shown in Fig. 1 for and . There, depending on the vorticities, fractional vortices can either be completely split apart or partially overlapping, as for example in . The topological properties when fractional vortices overlap are essentially different than when they do not. Indeed the topological charge (18) is quantized only if there is no core overlap. To emphasize the richness in unit flux structures, we display in Fig. 2 configurations with other ratios of the commensurate charges. All of these have non-trivial, very different signature of the magnetic field. The structure of the molecule-like bound state depends not only on the ratio of the gauge couplings, but also on the ratio of the densities associated with each condensate. This can be seen from additional regimes we displayed in Appendix A.

As discussed later on in Section III, depending on the ratio of densities in both condensates, the structure of the single flux quantum topological defect can be quite different. Indeed for substantial disparity in densities the ‘symmetric molecule’ structure, dominated by long range quadrupolar mode of the relative phase , is no longer preferred. Instead, an ‘asymmetric molecule’ is formed and it is characterized by a longer range dipolar mode of the relative phases. Such regimes with disparity in the condensate densities can be seen in Fig. 3. More regimes are given as additional material in Appendix A.

### A closer look to liquid metallic deuterium-like system

We argued that our model captures topological aspects of mixtures of condensates with commensurate charges and that it qualitatively applies to various systems and for example to the superconducting state for liquid metallic deuterium (LMD), where one expects coexistence of electronic Cooper pairs and Bose condensate of deuterons. Since reliable microscopic parameters for this state are not available we will study a phenomenological Ginzburg-Landau model for such a mixture of charged condensates. The deuterium nucleus (deuteron) is a spin-1 particle that can condense in several states Oliva and Ashcroft (1984a, b). Here we ignore spin degrees of freedom and treat it as a scalar charged condensate carrying electric charge , while electronic Cooper pairs carry charge . The mass of the electronic Cooper pairs is while the mass of the deuteron is . Let and respectively denote the deuteronic and electronic condensates, so . Because of the electric neutrality at zero temperature we consider

(19) |

(note for this expression we restore the symbols from (3)). Thus in this regime the electronic condensate is responsible for of the screening of the flux. In Fig. 3, we show a single flux quantum topological defect for big disparity in ground state densities that are likely to occur for liquid metallic deuterium. There, the arrangement of fractional vortices is somewhat different from those displayed in Fig. 1. Indeed, unlike previously, some of the fractional vortices overlap. The relative phase corresponding to the regime Fig. 1 assumes quadrupolar structure. The relative phases corresponding to Fig. 3, instead shows a dipolar structure. This is shown in Fig. 4. Dipole modes are longer range than quadrupole modes. This change in the long range behaviour of the relative phases should result in important modification of the large scale vortex matter structures. Long range dipolar modes were shown to play important role, although in a different context, on large scale vortex structure formation Garaud et al. (2014). As discussed below in Section III, the modification of the long range modes is consistently reproduced in the London approximation.

Here two constituent fractional vortices overlap. As a result, the topological invariant (18) here is not quantized. This can be heuristically understood by the fact that the target sphere is not completely covered. With our crude estimates the ratio of condensate densities for liquid metallic deuterium, is about the same as for liquid metallic hydrogen. However because the commensurate charges are different, the topological excitations are completely different. The lowest energy topological excitation in liquid metallic hydrogen-like system is an axially-symmetric composite vortex, while for LMD it is a composite object of two co-centered vortices plus one satellite fractional vortex. The magnetic signature of the topological defect in LMD looks like a pair of vortices, while it is a single vortex for LMH.

### The case of incommensurate charges

The conditions for finite energy solutions (II) and the flux quantization (11) rely on the fact that charges are commensurate, that is, that their ratio is a rational number so that they can be parametrized as where are integer numbers. For generality, here we address the question of what changes if charged condensates have incommensurate electric charges. When stands for elementary charges of elementary particles, they are integer multiples of an elementary electric charge. With current progress in creation of artificial gauge fields it could not be ruled out that systems with incommensurate coupling to the gauge field may be artificially realized.

For this exercise, we have to relax the condition that are both integer numbers. Since have to be integer for the ’s to be single valued, the condition (II) ensuring both finite energy and flux quantization cannot be satisfied. As a result, the elementary vortices carry different flux that cannot be added together to add up to a flux quantum. More precisely, the total flux (II) is

(20) |

and it is impossible to (consistently) write this as an integer times a flux quantum . Correspondingly, when charges are incommensurate, it is not possible by any mean, to eliminate the winding in the neutral sector (13c). Thus there are no finite energy solutions (in infinite domain). Instead solutions have logarithmically divergent energy due to the (superfluid) mode associated with the neutral sector. That is, a phase gradient resulting from a phase winding always causes logarithmic divergence of vortex energy and cannot be fully compensated by the vector potential.

Let us now close this remark about mixtures of condensates with incommensurate charges, and focus on the case where the ratio of the coupling constants is a rational number. That is where are integers. We found that topological defects carrying integer flux are bound states of different fractional number of vortices in the different condensates. Because of the frustrated interactions that vortices in the same condensate repel, while they try to overlap with vortices in the other condensate, these bound states arrange into very complicated molecule-like structures. A large part of this exotic vortex structures can be captured by investigating the London limit, where fractional vortices can be mapped to Coulomb charges.

## Iii London limit

When the electromagnetic repulsion is strong enough, integer vortices split to form a bound state of fractional vortices. The underlying physics describing the core splitting can be accurately captured within the London approximation where everywhere (except for a sharp cut-off at vortex core). There, the expression (13) further simplifies to

(21a) | ||||

(21b) |

The interaction energy of two non-overlapping fractional vortices is approximated in this London limit by considering charged (21a) and neutral modes (21b), separately. The energy of the charged sector (21a) reads as

(22) |

where the London penetration length is . The London equation for a (point-like) vortex placed at and carrying a flux is

(23) |

and its solution is

(24) |

where is the modified Bessel function of second kind. For two vortices located at and , respectively carrying fluxes and , the source term in the London equation reads as and the magnetic field is the superposition of two contributions . Thus

(25) |

and is the (self-)energy of the vortex . Finally, the interaction energy of two vortices in components reads as

(26) |

The interaction through the charged sector is thus screened interaction given by the modified Bessel function. When the couplings are parametrized such that they have the same sign (since the theory is invariant under complex conjugation, this is always possible), this interaction is always positive for any having the same sign of vorticity. It then gives, repulsive interaction between any kind of fractional vortices with co-directed winding. That is vortices repel while a vortex and an anti-vortex attract each other. On the other hand, the interaction through the neutral sector is attractive (resp. repulsive) for fractional vortices of the different (resp. same) condensate. The energy associated with the neutral mode (21b) reads as

(27) |

To evaluate the interaction between fractional vortices in different condensates and respectively located at and , the neutral sector is expanded

(28) |

At sufficiently large distance, a phase winding around some singularity located at a point , is well approximated by . Thus

(29) |

As a result, the interaction part in the neutral sector reads as

(30) |

Similarly, the interaction between two vortices in the same condensate is computed by requiring that the phase be the sum of the individual phases , while . Then the interaction reads as

(31) |

with here . To summarize, the interaction of vortices in different condensates is

(32) |

while interactions of vortices of similar condensates are

(33) |

with , and the sample size. Choosing the energy scale to be and defining the parameters

(34) |

The interaction between fractional vortices reads as

(35) |

Thus vortex matter in the London limit of a two-component superconductor with incommensurate charges is described by a 3-parameter family . This is illustrated in Fig. 5.

The interaction between vortices in the same condensates is repulsive. In multi-component superconductors where both condensates have the same number of vortices in each component, vortices in different condensates will attract each other to form a bound state of co-centered vortices that minimizes the energy cost of the neutral sector Babaev (2002); Smiseth et al. (2005). The situation here is more subtle. Because both condensates have different number of (fractional) vortices, the system has to compromise between vortices in similar condensates that repel each other and the fact that vortices in different condensates try to overlap. This explains why, beyond the London limit (see e.g. Fig. 1), integer vortices form a molecule-like bound state of split fractional vortices.

### Transition in the structure of skyrmion

In the context of the mapping to point charges, the finite energy condition (II) is equivalent to require charge neutrality. Thus we complete the mapping to point charges by defining the electric charge and . A neutral set of charged particle thus satisfy the charge neutrality

(36) |

We already know, from our solutions of the full non-linear model, that vortex solutions do exist and they have both finite energy and carry unit flux quanta. We also observed that, provided is small enough, that cores of fractional vortices are not superimposed. We now try to reproduce our results, using the equations (III) of the London limit. Here since a vortex configuration is a neutral set of discrete charge, it can be globally described by its dipole and quadrupole moments in two dimensions

(37) |

The minimum of the interaction energy (III) should describe the location of vortex cores. We apply this for a single vortex with . Thus it is described by a set of three point particles, two carrying a single positive charge and one twice negatively charged. According to the interaction energies, it forms a bound state of non-overlapping particles, see Fig. 6. Remarkably there is a transition of the ‘molecular’ structure at a given . For , the molecule is symmetric and it has no dipole moment. The long-range interactions are thus dominated by quadrupole modes. For sufficient disparity in densities, when , the least energetic arrangement is no longer symmetric and thus the vortex molecule develops a dipole moment that is long range and should dramatically alter the large scale structures. Indeed, quadrupole modes are smaller in amplitude and also decay faster than dipole modes. For examples of the effect of long-range dipolar interactions in large scale structures (in a different context), see Garaud et al. (2014).

In our estimates to describe liquid metallic deuterium, the London limit parameter . Thus according to the London limit picture Fig. 6, the single flux quantum vortex in the superconducting phase of liquid metallic deuterium should sit in the dipolar regime. As we illustrated in Fig. 3, this is indeed the case, and there is perfect agreement between the London limit picture and direct numerical simulations.

Note that since the model describes independently conserved condensates, it has invariance. Then both condensates have in principle different critical temperatures at which they condense. Then, there is in principle also always a regime where one of the condensates has much less density than the other. So there is always, at least, a regime where (resp. ) if (resp. ) condenses first. So there is always a phase dominated by the long range dipolar mode in the relative phases, that should dramatically influence large scale structuring of the vortex matter.

### The case of many particles

The physics of the topological defects carrying a single flux quantum is shown to be quite rich. Indeed even in the simplest case when , the single skyrmion has a transition in its internal structure (see Fig. 6). The emerging long-range modes should have very important influence on the many-skyrmion states. Investigating the many-skyrmion states, may give valuable information about the transport properties or about the lattice structures and their melting. Within the London limit, such properties can be investigated using molecular dynamics or Monte-Carlo simulations of the point particles interacting according to (III). Although the point-charge model does not completely capture all the underlying physics, it can reproduce several aspects of the structures obtained beyond the London limit. This is beyond the scope of the current paper, yet we can address few general comments about the case of many particles.

In order to investigate the many-body properties of our model, one approach is to model fractional vortices by point charges with elementary interactions (III), that is, to consider a set of particles corresponding to (fractional) vortices in the condensate . For a system containing integer flux, the number of particles should satisfy the relation . The interacting energy in the case of many particles reads as

(38) |

where denotes the position of the -th vortex of the condensate and interaction energies are given by (III). Note that this problem is related to the problem of unconventional plasma discussed in the context of quantum Hall states Bonderson et al. (2011); Herland et al. (2012, 2013).

The many-particle problem (III) can be investigated using different standard techniques such as molecular dynamics or Monte Carlo simulations. In the light of the complicated structure of the single skyrmions, one may expect very rich phases of the vortex matter there. However this deserves full investigation that is beyond the scope of the current paper.

## Iv Conclusions

We investigated physical properties of mixtures of charged (bosonic) condensates, carrying different electric charges. More precisely, we introduced a Ginzburg-Landau model that accounts phenomenologically for such mixtures. Disregarding the underlying microscopic theories that describe mixtures of charged condensates, this model is expected to qualitatively describe the topological excitations therein.

Elementary topological excitations are fractional vortices, that is, vortex configurations with winding in only one condensate. Because of the existence of a neutral mode, describing relative counter-directed motion of particles, fractional vortices have logarithmically divergent energy. The condition for having a finite energy solution in a mixture of condensates having commensurate electric charges and is that the phase should wind times in and times in . Because of the commensuration of electric charge, finite-energy configurations have different number of fractional vortices in different condensates.

Fractional vortices in the same condensate repel while fractional vortices in different condensates attract each other in order to reduce the energy cost associated with the counterflow of charge carriers. As a compromise, the topological excitation carrying an integer flux quantum can form a molecule-like bound state of fractional vortices, where there is no overlapping of vortices in contrast to systems with the same charges Babaev (2002); Garaud et al. (2014).

We also addressed the question of the underlying topology. There, we showed that two configurations carrying an integer flux quantum are differentiated from each other by a topological invariant. The topological excitations were explicitly constructed numerically and their structure, namely the spatial arrangement of constituent fractional vortices can be understood by investigating the London limit physics, where fractional vortices are mapped to point Coulomb charges.

The model we introduced and its topological excitations applies, at least qualitatively, to various physical systems where different condensates are formed and where they are commensurately coupled to the vector potential of the magnetic field. Namely it could effectively describe projected superconducting state of liquid metallic deuterium where deuterons form a charged Bose-Einstein condensate mixed with electronic Cooper pairs. This state of matter is currently a subject of experimental pursuit Eremets and Troyan (2011). Since these experiments are conducted in diamond anvil cell that can be equipped with a receiving coil, the finding which we report could help to confirm or rule out formation of this state. Mixtures of commensurately charged condensates might also be an interesting system to be realized in cold atoms with synthetic gauge fields.

We acknowledge fruitful discussions with Johan Carlström, Karl Sellin, Daniel Weston and especially with J.M. Speight. This work is supported by the Swedish Research Council, by the Knut and Alice Wallenberg Foundation through the Royal Swedish Academy of Sciences fellowship and by NSF CAREER Award No. DMR-0955902. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at National Supercomputer Center at Linköping, Sweden.

## Appendix A Additional material

## Appendix B Finite element energy minimization

We consider the two-dimensional problem (4) defined on a domain bounded by . In our simulations, we choose the domain to be a disk. The problem is supplemented by the boundary condition with the normal outgoing vector on . This condition physically implies that no current flows through the boundary. This is thus a superconductor/insulator or superconductor/vacuum boundary condition. Since this boundary condition is gauge invariant, additional constraint can be chosen on the boundary to fix the gauge. Our choice is to impose the radial gauge on the boundary (note that with our choice of domain, this is equivalent to ). This choice eliminates (most of) the gauge degrees and the boundary condition separates into two parts

(B.1) |

Note that these boundary conditions allow a topological defect to escape from the domain. To prevent this in simulations of individual skyrmions or skyrmion groups when no field is applied, the numerical grid is chosen to be large enough so that the attractive interaction with the boundaries is negligible. The size of the domain is then much larger than the typical interaction length scales. Thus, with this method one has to use large numerical grids, which is computationally demanding. This guarantees that the solutions are not boundary pressure artefacts. In particular this means that the observed core splitting cannot be attributed to finite size effects as in mesoscopic samples Chibotaru et al. (2007); Geurts et al. (2008); Chibotaru and Dao (2010).

The variational problem is defined for numerical computation using a finite element formulation provided by the Freefem++ library Hecht (2012). Discretization within finite element formulation is done via a (homogeneous) triangulation over , based on Delaunay-Voronoi algorithm. Functions are decomposed on a continuous piecewise quadratic basis on each triangle. The accuracy of such method is controlled through the number of triangles, (we typically used ), the order of expansion of the basis on each triangle (second order polynomial basis on each triangle), and also the order of the quadrature formula to compute the integral on the triangles.

### Initial guess

The initial field configuration carrying flux quanta is prepared by using an ansatz which imposes phase windings around spatially separated vortices in each condensate:

(B.2) |

where , is the ground-state density of a given condensate, and its ground-state phase. Because of invariance, both can be chosen to be zero. parametrizes the size of cores while the functions

(B.3) |

denotes the position of the singularity of the -th vortex of the condensate. The starting configuration of the vector potential is determined by solving numerically Ampère’s equation on the background of the superconducting condensates given by Equations (B)–(B).

For a given starting configuration, the free energy is then minimized with respect to all degrees of freedom, with the condition (B.1) that no current flows through the boundary. Here we used a non-linear conjugate gradient method. The algorithm was iterated until relative variation of the norm of the gradient of the functional with respect to all degrees of freedom was less than .

### Footnotes

- A reservation should however be made here. Mesoscopic systems, because they introduce a boundary cut-off on the divergent mode, can allow fractional vortices to be thermodynamically stable Chibotaru et al. (2007); Chibotaru and Dao (2010); Geurts et al. (2008). Also fractional vortices can also be thermodynamically stable near boundaries Silaev (2011). In our paper we however only consider bulk systems. In particular the simulation grid is chosen to be much larger than the size of vortices and thus boundary effects are irrelevant.
- Note that this physics of split vortices is entirely different from that of spatially separated fractional vortices that may also exist in different systems due to biquadratic density-density interaction Nitta et al. (2014); Agterberg et al. (2014); Kobayashi and Nitta (2014), dissipationless drag interaction Chung et al. (2007); Garaud et al. (2014) or coexistence of vortices and domain walls Garaud et al. (2013, 2011). In the case of mixtures of condensates with commensurate charges, the splitting occurs because of competing attractive interaction through the neutral sector and repulsive interaction mediated by the magnetic field. As discussed below, unlike for core splitting induced by bi-quadratic terms, this physics is well captured in the London limit.

### References

- E. V. Herland, E. Babaev, and A. Sudbø, Phys. Rev. B 82, 134511 (2010).
- E. Babaev, A. Sudbø, and N. W. Ashcroft, Nature 431, 666 (2004).
- E. Babaev and N. W. Ashcroft, Nat. Phys. 3, 530 (2007).
- J. Oliva and N. W. Ashcroft, Phys. Rev. B 30, 1326 (1984a).
- J. Oliva and N. W. Ashcroft, Phys. Rev. B 30, 5140 (1984b).
- P. Bedaque, M. Buchoff, and A. Cherman, Journal of High Energy Physics 2011, 1 (2011).
- M. I. Eremets and I. A. Troyan, Nat Mater 10, 927 (2011).
- D. Banerjee, M. Dalmonte, M. Müller, E. Rico, P. Stebler, U.-J. Wiese, and P. Zoller, Phys. Rev. Lett. 109, 175302 (2012).
- E. Zohar, J. I. Cirac, and B. Reznik, Phys. Rev. Lett. 109, 125302 (2012).
- E. Babaev, Nuclear Physics B 686, 397 (2004).
- A. B. Kuklov and B. V. Svistunov, Phys. Rev. Lett. 90, 100401 (2003).
- E. Smørgrav, E. Babaev, J. Smiseth, and A. Sudbø, Phys. Rev. Lett. 95, 135301 (2005).
- D. F. Agterberg and H. Tsunetsugu, Nature Physics 4, 639 (2008).
- E. Berg, E. Fradkin, and S. A. Kivelson, Nat Phys 5, 830 (2009).
- E.-G. Moon, Phys. Rev. B 85, 245123 (2012).
- S. Sachdev, Nat. Phys. 4, 173 (2008).
- K. Chen, Y. Huang, Y. Deng, A. B. Kuklov, N. V. Prokof’ev, and B. V. Svistunov, Phys. Rev. Lett. 110, 185701 (2013).
- F. Hecht, J. Numer. Math. 20, 251 (2012).
- E. Babaev, Phys. Rev. Lett. 89, 067001 (2002).
- A reservation should however be made here. Mesoscopic systems, because they introduce a boundary cut-off on the divergent mode, can allow fractional vortices to be thermodynamically stable Chibotaru et al. (2007); Chibotaru and Dao (2010); Geurts et al. (2008). Also fractional vortices can also be thermodynamically stable near boundaries Silaev (2011). In our paper we however only consider bulk systems. In particular the simulation grid is chosen to be much larger than the size of vortices and thus boundary effects are irrelevant.
- Note that this physics of split vortices is entirely different from that of spatially separated fractional vortices that may also exist in different systems due to biquadratic density-density interaction Nitta et al. (2014); Agterberg et al. (2014); Kobayashi and Nitta (2014), dissipationless drag interaction Chung et al. (2007); Garaud et al. (2014) or coexistence of vortices and domain walls Garaud et al. (2013, 2011). In the case of mixtures of condensates with commensurate charges, the splitting occurs because of competing attractive interaction through the neutral sector and repulsive interaction mediated by the magnetic field. As discussed below, unlike for core splitting induced by bi-quadratic terms, this physics is well captured in the London limit.
- E. Babaev, L. D. Faddeev, and A. J. Niemi, Phys. Rev. B65, 100512 (2002).
- E. Babaev, Phys. Rev. B79, 104506 (2009).
- J. Garaud, K. A. H. Sellin, J. Jäykkä, and E. Babaev, Phys. Rev. B 89, 104508 (2014).
- J. Smiseth, E. Smørgrav, E. Babaev, and A. Sudbø, Phys. Rev. B 71, 214509 (2005).
- P. Bonderson, V. Gurarie, and C. Nayak, Phys. Rev. B 83, 075303 (2011).
- E. V. Herland, E. Babaev, P. Bonderson, V. Gurarie, C. Nayak, and A. Sudbø, Phys. Rev. B 85, 024520 (2012).
- E. V. Herland, E. Babaev, P. Bonderson, V. Gurarie, C. Nayak, L. Radzihovsky, and A. Sudbø, Phys. Rev. B 87, 075117 (2013).
- L. F. Chibotaru, V. H. Dao, and A. Ceulemans, EPL (Europhysics Letters) 78, 47001 (2007).
- R. Geurts, M. V. Milošević, and F. M. Peeters, Phys. Rev. A 78, 053610 (2008).
- L. F. Chibotaru and V. H. Dao, Phys. Rev. B 81, 020502 (2010).
- M. A. Silaev, Phys. Rev. B 83, 144519 (2011).
- M. Nitta, M. Eto, and M. Cipriani, Journal of Low Temperature Physics 175, 177 (2014).
- D. F. Agterberg, E. Babaev, and J. Garaud, preprint (2014), arXiv:1403.6655 [cond-mat.supr-con] .
- M. Kobayashi and M. Nitta, Journal of Low Temperature Physics 175, 208 (2014).
- S. B. Chung, H. Bluhm, and E.-A. Kim, Phys. Rev. Lett. 99, 197002 (2007).
- J. Garaud, J. Carlström, E. Babaev, and M. Speight, Phys. Rev. B 87, 014507 (2013).
- J. Garaud, J. Carlström, and E. Babaev, Phys. Rev. Lett. 107, 197001 (2011).