Vortex nucleation barriers and stable fractional vortices near boundaries in multicomponent superconductors

The magnetization process of a superconductor is determined by the potential barrier for vortex nucleation and escape. In multicomponent superconductors, fractional vortices with a winding in the phase of only one of the components can be stable topological solitons that carry a fraction of the flux quantum. While the formation of such objects in the bulk costs logarithmically or linearly divergent energy, these objects were shown to be stable near samples' boundaries in the two-component London model. Therefore, the conventional Bean-Livingston picture of magnetic flux entry does not apply to these superconductors, since the entry process can involve fractionalization of a vortex. In this paper, we address the nonlinear problem of determining the potential barrier for fluxoid penetration in a multicomponent superconductor, including the effects of various intercomponent couplings, by using the recently developed gauged string method. The method allows numerically exact (i.e., convergent) calculation of a sphaleron configuration in a gauge theory and thus the height of the nucleation barrier. We show how the fractionalized nucleation processes result in multiple sphalerons and intermediate states due to the complex shape of the energy landscape of multicomponent superconductors.


I. INTRODUCTION
Most superconductive systems of current interest have multiple components.That include multi-band superconductors with s-wave pairing, see e.g.[1][2][3][4][5][6], superconductors with broken time-reversal symmetry, which can be singlet s-wave superconductors [7][8][9][10], and superconductors materials with unconventional pairing [11].When a superconductor has multiple components, the vortex excitations have a composite character: they are bound states of fractional-flux vortices, which, individually, have logarithmically divergent energy due to intercomponent electromagnetic coupling, or linearly divergent density due to phase-difference-locking coupling, and therefore, form finite-energy bound states [12].This affects the magnetic properties of such materials since, to enter a superconductor, a vortex needs to overcome a potential barrier.
In the Bean-Livingston picture, the origin of the surface barrier for a vortex is described within the London model, as the effect of the competition between the attractive force of a mirrored image-vortex with opposite vorticity, and the repulsive force with the surface current induced by the external magnetic field [13].However, in general, the vortex position in the maximum energy configuration can be very close to the boundary, preventing the use of the superposition principle that holds only for linear theories.Therefore, a general description requires the solution of the full nonlinear problem, e.g. in a Ginzburg-Landau model.
These energy barriers correspond to infinitedimensional saddle points of the free energy landscape of the theory, and are called sphalerons, which are often associated with a transition between topologically distinct minima [14][15][16].These saddle points in the configuration space can be better understood as the local maximum of the minimum energy path of the nucleation process.That is defined as a path in the configuration space such that it crosses the minimum in the cotangent space of the path point by point.
In two-dimensional single-component systems, the number of processes featuring sphalerons is limited because the intervortex interaction is always repulsive, while the low dimensionality excludes processes involving complicated real-space topology changes like vortex reconnections.Indeed, sphalerons appear only in changes of the vortex lattice or in processes where vortices interact with the boundaries of the sample, like in the nucleation process.
When a second complex field comes into play, the energy landscape becomes more complicated and many new configurations, both stable states and sphalerons, are possible.In the nucleation process of a composite vortex, the vortices in the two components can experience attraction and repulsion to the boundary at different distances.The presence of different zones of attraction and repulsion can lead to the presence of barriers and metastable intermediate states.Vortex entry in a multicomponent superconductor is a much more complex process that, under certain conditions, involves the formation of a stable fractional vortex near the boundary.Stable fractional vortices in two-component systems were demonstrated both within the London model [17] and numerically in the Ginzburg-Landau model [18].
While many works, both analytical and numerical, focused on the calculation of the nucleation fields, that is the magnetic field needed to suppress the barrier [19][20][21][22][23], the calculation of the barrier height and the sphaleron configuration is very complicated and, until recently, there has been no analytical nor numerical fully-controllable approach to address this problem in gauge theories such as Ginzburg-Landau models.One of the many methods to compute the minimum energy path for other models is the string method [24,25], originally designed to compute saddle points in molecular dynamics problems.In Ref. 26 it has been proposed a generalization of the string method for classical gauge theories that allows calculating, in a numerically controllable way, vortex nucleation processes in a single-component Ginzburg-Landau model.
In this paper, we generalize to the multicomponent case the gauged string method presented in Ref. 26 that allows studying, extensively and free of approximations, the process of vortex nucleation and to calculate vortex nucleation barriers.The main objective is to study how nucleation is affected by the presence of a second component, looking for fractional vortex nucleation and edge pinning.Secondarily, this allowed the testing of the gauged string method for a system with more complicated gauge symmetry.

II. MODEL
Multicomponent Ginzburg-Landau theories have beeen microscopically derived for a variety of models [18,[27][28][29][30][31][32][33].Here, we consider a two-component model in which the order parameter Ψ is represented by two complex scalar fields (ψ 1 , ψ 2 ).The two component fields are coupled by some direct interaction V int , and to the electromagnetic field described by the gauge potential A. The system has a free energy where the free energy density functional reads, in natural units, ( Here the index k labels the axes of the coordinate system and the index α denotes the superconducting components, the vector field A is the gauge field while ψ α = |ψ α |e iθα are the two matter fields, and H is the external magnetic field.The parameters m αk are the masses of the two fields in the reference system.The operator D kα = −i∂ k − q α A k is the covariant derivative where q α is the coupling constant of the component to the gauge field that parameterize the magnetic field penetration depth.For the direct interactions term, we only consider the linear Josephson coupling, controlled by the parameters η, and the density coupling γ.Higher gradient terms, additional fields, and more general interactions, including the effects of time-reversal symmetry breaking and nematicity in the interaction potential can be straightforwardly included.
The vortex topological invariants of the model are represented by the two winding numbers defined as: If we neglect the direct interaction term V int (ψ 1 , ψ 2 ) the system described by the functional in Eq. 1 models two charged superfluids, which only interact through the gauge field.In the absence of phase-locking or phase separation, the system has U (1) × U (1) gauge symmetry.The magnetic flux carried by a vortex line, in the model that we consider, is where ψ α = −aα bα is the bulk equilibrium value of each order parameter, and N α is an integer.In conventional superconductors, q 1 = q 2 = q and therefore it simplifies to where Φ 0 = 2π/q is the magnetic flux quantum.
The elementary excitations in the model are fractional vortices, which have a nonzero winding only in one phase [12].Fractional vortices are thermodynamically unstable in the bulk [12] but they can be stabilized near the boundary of a sample by the interaction with surface currents [17].Only when the vortex carries an integer number of flux quanta, i.e.N 1 = N 2 , then the topological excitation is given the name of composite vortex and it is also stable in the bulk of the superconductor.
Another case, that we consider, is a system with condensates having different charges, that we parametrize as q i = g i q.With the simplifying assumption ψ 1 = ψ 2 , the flux carried by a vortex is Also in this case, the only stable soliton in the bulk is the one having an integer number of flux quanta, which are composite vortices with a (N 1 , N 2 ) = (g 1 , g 2 ) structure.We show later that, when the topological excitations do not satisfy the integer quantization of flux, they can exist only near the system's boundaries, because of a stabilization mechanism similar to the one for fractional vortices in conventional two-component superconductors.
In models with direct intercomponent interactions, i.e.V int = 0, it is not possible to associate separate coherence length to each component as the coherence lengths describe the decay of linear combinations of the matter fields [35,39,41].In this scenario, is more convenient to speak of a large-core component and a small-core component based on the typical dimension of the fractional vortex cores, which we determine numerically.
Here, other types of solitons can be present.For example, two well separated fractional vortices can stick together.These objects can be described as skyrmions of the pseudospin vector field S = Ψ † σΨ Ψ † Ψ where Ψ = (ψ 1 , ψ 2 ) T where and σ are the Pauli matrices vector [18,42].
An extremely useful tool to study the formation of solitons is the concept of minimum energy paths.A generic path in the configuration space of the system is a function q(s) = (A(s), ψ 1 (s), ψ 2 (s)), where the transition coordinate s ∈ [0, 1] has been introduced.A minimum energy path can be defined in a variational formulation as where F ⊥ is that the free energy functional F restricted to the cotangent space of the trajectory in the point q MEP (s) [43].Once the minimum energy path q MEP (s) is obtained, it is possible to compute the quantities of interest, like the free energy F or the winding numbers N i , and track their evolution along the path.The evolution along s can be considered a pseudodynamics, which mimics some features of the real time dynamics simulated by time-dependent models.We calculate minimum energy paths in the two-component Ginzburg-Landau model through the gauged string method introduced in Ref. 26.The details on the numerical implementation of this method can be found in the supplementary material of Ref. 26.For the reparametrization step of the method, we used the Euclidean distance calculated for the observable magnetic field B = ∇ × A, and currents 2m αk ψ α + c.c..In this work, we mainly study the minimum energy paths of the single vortex nucleation process.The energy of the principal sphaleron, the most important barrier in the process, is just the maximum of free energy along the minimum energy path F S = max s F (s).We define F M as the energy of the Meissner state, which is a state with no vortices in the bulk.When studying the nucleation of the first vortex in the sample, the nucleation barrier is therefore defined as ∆F n = F S − F M .
When looking at composite nucleation processes, the final point is a configuration when the vortex is in the center of the system.We denote the energy of this state as F C .If no surface bond states are present, the escape energy of the vortex is simply given by ∆F e = F S − F C .When considering fractional vortices, the most stable position is not in the center of the sample but usually near the surface, where the vortex forms a bond with surface currents.We label the free energy of this configuration by F B .
To characterize the surface-bonded fractional state, we are interested in evaluating the strength of the bond to the surface, by comparing the energy of the system with the same quantity when the vortex is in the bulk.However, the state with the vortex in the center of the bulk is not a stable state of the system but rather an unstable equilibrium point.For this reason, a minimum energy path cannot, in principle, be calculated.To circumvent this problem, we remove a single point from the superconductive domain to artificially introduce an extremely shallow local minimum, with negligible effects on the total energy of the field configuration.To minimize the perturbation of the free energy we choose point (0, 0) that is the farthest from the boundaries.We denote the energy of the configuration with the vortex pinned in the exact center of the domain by F C also in this case.
For the fractional nucleation curve, the escape barrier is then defined as the energy necessary to expel the vortex starting from the local minimum, ∆F e = F S − F B , while the bonding energy ∆F b = F C −F B is the energy required to move the vortex from the edge toward a region deep in the bulk.All these quantities are more clearly shown in Fig. 1.
For the simulations, we consider a system with computational domain Ω = [−L/2, +L/2] 2 with L = 12.We build the initial guesses for the component fields by moving the singularities from a point outside of the domain to the center.Notice that, the initial guess for the endpoint configurations at s = 0 and s = 1 are topologically equivalent to the ones in the converged minimum energy path as their motion follows a simple gradient descent.

III. RESULTS
The most general two-component Ginzburg-Landau model has a huge parameter space.For this reason, we To pin the vortex in the center of the domain the origin (x, y) = (0, 0) has been removed from the superconductive domain.This is necessary since the final configuration is not a minimum of the system.With this trick, the vortex is within a local minimum with negligible effects on the total energy.
show the results for some representative systems.For each series of simulations, we fix all parameters according to Tab.I except for one or a pair.First, we focus on the case where the intercomponent direct interaction is absent, i.e.V int = 0.The characteristic length scales for the systems considered are shown in Tab.II.We examine the nucleation in various systems: System E features two perfectly similar complex scalar fields in a type-2 regime, System U introduces different length scales but still belonging to the type-2 class.System H belongs to the type-1.5 class, while System D introduces condensates with different charges.In Systems U, H, and D, the two component fields have a different set of parameters and thus we identify ψ 1 as the large-core component while ψ 2 is the small-core one.

A. Surface bonded fractional vortices in type-2 systems
The first case we address is the vortex nucleation processes in isotropic two-components type-2 systems.
In the case of a superconductor with two matter fields, having different characteristic lengths, like System U, we observe a fractionalized nucleation.This means that the optimal nucleation process is, in general, split into two steps as the fractional vortices in each component enter the bulk in different stages.This is qualitatively in agreement with the earlier studies using the London model [17].This behavior can be observed by computing the winding number for each field along the path N i (s) and comparing this curve with the free energy F (s), as shown in Fig. 2. Indeed, the small-core vortex nucleates first, followed by a vortex in the large-core component.The two fractional vortices are bound together in a composite vortex which then, pushed by the repulsion of surface currents, moves towards the center of the domain.
If, instead, we consider the nucleation of the fractional vortex with nonzero winding in the large-core component only, we see that a surface bonded state can exist in a region of the parameters space.The minimum energy path of the fractional vortex nucleation of this kind is shown in Fig. 3.
By studying the changes in the minimum energy path for fractional vortex nucleations, as a function of the external magnetic field, we verify that the stabilization mechanism occurs due to the interaction with the surface Meissner currents.The fractional vortex is attracted by the surface of the domain and forms a bond.As shown in Fig. 4, an increase of the external magnetic field makes the minimum deeper.
It is possible to analyze the phase space of the system by extracting the free energy barriers ∆F n , ∆F e , ∆F b as a function of the external field.In this way, we can identify a region where fractional vortices are stabilized by surface currents.In the considered regime, the bonding energy ∆F b is approximately constant with respect to the magnetic field, meaning that the stability of the surface-bonded fractional state is weakly affected by the external magnetic field H, in the given field range.Fractional surface-bonded vortices can exist in the region between H = 0.7 and H = 1.4 while higher fields cause the nucleation of composite vortices.In this zone, the bonding energy is constant, while the escape barrier increases with the magnetic field intensity.We notice how, for fields higher than H = 1.2, the escape barrier exceeds the nucleation one, meaning that a fractional vortex is more stable than the Meissner state.This means that, for fields in the region starting from H = 1.2 to H = 1.4, a surface bonded fractional vortex phase is expected to be the equilibrium state of the system.This is consistent with the magnetization simulations of a slightly different model [18].In this zone, the bonding energy is constant while the escape barrier increase with the magnetic field intensity.Notice how, for fields higher than H = 1.2, the escape barrier exceeds the nucleation one, meaning that a fractional vortex is more stable than the Meissner state.The bonding energy is approximately constant with respect to the magnetic field.
It can be seen in Fig. 6, for the considered parameters, that, if the system has a vortex in the large-core component trapped near the surface, the activation energy for the nucleation of a vortex in the small-core component is usually low.The nucleation of a vortex in the other component frees the vortex from binding to the surface, as the composite vortex is repelled by the surface currents.This limits the lifetime of the surface-bonded state.
Comparison of nucleation minimum energy paths of a composite vortex (C), a fractional vortex (F) and the nucleation of a second component vortex starting from a fractional one (F → C), in System U at H = 1.1.The black line is the energy of the system with a fractional vortex with respect to the composite vortex free energy F C and the gray line is the energy of the Meissner state with respect F C .Notice how the nucleation of a fractional vortex (red line) has a barrier that is slightly smaller than the one of the composite vortex (blue line).Moreover, the nucleation process for a vortex in the second component starting from a fractional vortex bond to the surface has a very small barrier.This suggests that, in the given example, the fractional surface vortices are unstable to small fluctuations.Also in this case, the dotted line indicates the absence of winding in both components, the dasheddotted line signals the presence of a fractional vortex, while the continuous line means that there is a composite vortex in the system.

B. Anisotropy effects
The introduction of anisotropies in the model affects the nucleation process since it induces geometrical deformations in the vortices' cores.One way of introducing anisotropies in the model is to alter the ratio between the mass parameter of the fields.In the following, we consider the ratios m 2x /m 2y = 10, and m 2y /m 2x = 10.In all the simulations discussed in the paper, the vortex enters along the x direction.The two cases considered here can, indeed, be seen as rotations of the same system or equivalently nucleation processes at two different sides.
As it can be seen in Fig. 7, the nucleation barrier in the case m 2y > m 2x is lower and the state bonded to the surface is very stable and characterized by a fractional vortex with winding only in the isotropic field.On the contrary, if the heavy direction is parallel to the normal to the surface, the barrier is higher, and the vortex in the anisotropic component is the first to enter the sample.Hence, anisotropies provide a way to enhance the stability of vortices near the boundary.This can lead to the case where fractional vortices can be stable when laying near the edges of the system aligned with one of the principal axes while being unstable near edges oriented in the orthogonal direction.This situation is shown in Fig. 7, as only one of the minimum energy paths has a clearly defined minimum which is also lower than the initial state.In these simulations, we study the vortex entry from a surface aligned with the y direction.
In the case where the heavy direction is parallel to the surface, the bond state energy is enhanced by the vortex core deformation.To report the nucleation steps, we use the following scheme described in Fig. 2 C. Nucleation in type-1.5 systems In type-2 systems the nucleation process of consecutive vortices is not notably different from the nucleation of the first vortex, if the density of vortices is low enough.This is due to the repulsive interaction between composite vortices.
The situation is different in type-1.5 systems because of the long-range attractive interaction between vortices that tends to cluster them together.An additional nucleating vortex is, therefore, attracted by the composite vortices already present in the system.Hence, the nucleation is favored in the region of the boundary closer to the pre-existing vortex.
In Fig. 8, we compare the minimum energy path of the nucleation of the first composite vortex, with the nucleation of the second one.The nucleation barrier for the entry of a second composite vortex, in the considered example, is considerably lower than the one for the first vortex.Moreover, the presence of the attractive interaction between the vortices lowers the energy of the two-composite-vortex cluster.
Then, we consider the nucleation of a fractional vortex in the large-core component in the presence of a composite one in the bulk, i.e. (1, 1) → (2, 1).The nucleation barrier for such an event is comparable to the one for the nucleation of a second composite vortex, i.e.
(1, 1) → (2, 2) .We attribute this effect to the attractive interaction in the large-core component.The resulting soliton has a (2, 1) internal structure.Since it has a fractional flux, it is unstable in the bulk.However, it can be stabilized by the interaction with the boundary resulting in an unconventional surface bonded fractional soliton.Notice that this structure is more stable that the Meissner state, even if it is still metastable with respect to the two-vortex cluster.In the fractional (2, 1) soliton, the small core fractional vortex is shared between the two large-core vortices.The small energy barrier, combined with the limited increase in the system free energy for this configuration, increases the probability of seeing such solitons due to thermal fluctuations or by careful preparation of the state.This state can also be seen as containing a trapped skyrmion near the interface.

D. Condensates with different charges
All of the most common models for multicomponent superconductors assume that the charges of the condensate are equal.However, the multicomponent Ginzburg-Landau model can be extended to include systems where the charges of the two components are different but commensurate.These models describe a mixture of charged superfluids that can model more exotic systems like, e.g., the proposed coexisting electronic superconductivity with a Bose condensate of deuterons in liquid metallic deuterium [44] .In this case, unconventional vortex structures with a disparity in phase winding arise [45,46].
In these systems, the composite vortices have a complex structure, due to the different phase winding required to enclose a quantum of magnetic flux.We consider, as an example, a system in which one of two condensates has a double charge with respect to the other (System D).Here, the stable soliton is a cluster with a (2, 1) structure, composed of two fractional vortices in the first component, and one fractional vortex in the second component.Notice that, this structure, metastable in the conventional system, is now stable.While the (1, 1) composite vortex can now live only as a metastable solution near boundaries.
We study the nucleation process of such structure in Fig 9, where is it possible to observe how the nucleation process is strongly fractionalized.The entrance of each fractional vortex is, in fact, associated with an energy minimum.
Any deviation from the structure prescribed by Eq. ( 4) results in an unstable topological soliton in the bulk.However, unstable solitons can be stabilized by boundary effects, in a similar fashion to the surface bonded fractional vortices in the conventional two-component superconductors.We verified the stability of these structures by simulating the nucleation process of (1, 0) and (1, 1) vortices in Fig. 10.In both cases, a sizable barrier prevents the escape of the surface bonded solitons.

E. Effects of direct intercomponent interactions
In this section, we consider the effect of direct interactions, i.e.V int = 0, on the nucleation mechanism.As a start, we study the effect of bilinear Josephson on the nucleation process by simulating the nucleation of composite vortices in System E with the coupling term This term introduces phase-locking in the system such the spontaneously broken symmetry is reduced to U (1).For positive η, the minimum is obtained for θ 12 = ±π, while for negative values the minimum is for θ 12 = 0.In the following, we restrict the study to positive η as the sign can be absorbed in the phase of one field.
As show in Fig. 11 (a), the nucleation barrier gets higher with increasing coupling.The curves can be fitted reasonably well thought an exponential regression using ∆F ] as the ansatz, where C 1 , C 2 , and H n are the coefficients estimated for each value of η.In this way, we estimated the nucleation field shown in Fig. 11 (b).From the simulations, it appears that, in the considered regime, the nucleation field H n (η) increase linearly as a function of η.Note that the fractional vortex interaction potential is asymptotically linear in the presence of this coupling.
To check the effect of bilinear Josephson coupling on fractional surface vortices, we simulated the minimum energy path of fractional vortex nucleation for System U and H = 1.2.The bilinear Josephson coupling has a detrimental effect on surface-bond states.Already for η = 0.01, the fractional surface vortex becomes a metastable state with energy higher than the Meissner FIG. 10.Surface bond soliton in System D. Any soliton with a winding in the large core component that is smaller than the natural number (2 in this case), is unstable in the bulk and can be present only in a metastable state while bonded to the surface.The top panel reports the free energy for the first nucleation (green) and for the second nucleation (red) as a function of the parametrization parameter s.The contour plots report the system's order parameters ψ1 and ψ2, and magnetic field B at specified configurations in the minimum energy path marked by a blue dot in the top panel.
state.Increasing further η, we find that, at values around 0.04, the nucleation barrier for a vortex in the second component reaches zero such that the surface-bond fractional vortex state becomes unstable.
Let us now consider the density coupling this term does not break the U (1) × U (1) gauge symmetry.Depending on the sign of γ, this term promotes an attraction or repulsion between fractional vortices in different components.A repulsive (positive γ) interaction can lead to the formation of skyrmions in the system [18,47].In this case, the composite vortex formed by exactly coinciding vortex cores is a metastable soliton, while the stable skyrmion is formed by spatially separated bonded vortices.We observe that the presence of repulsive density coupling decreases the nucleation barrier and fields, as shown in Fig. 12.Since the fractionalized entry is favored, even in the case of System E, where the two fields are perfectly equal, the metastable field configuration featuring only one fractional vortex can be viewed as a partialskyrmion state.
We have tested the effect of density coupling for System U (with H = 1.2) through a simulation of the decay of the fractional surface vortex by nucleation of a second fractional vortex in the small-core component.We observe that the nucleation barrier for this process decreases as γ increases.Moreover, the metastable minimum state is found for a vortex further from the surface.For γ 0.5 the surface vortex states spontaneously decay and multiple vortices in the small-core component nucleate.An example of a skyrmion entry process is shown in Fig. 13, for system H with a magnetic field close to H c1 .Nonetheless, for some type-2 systems, the partialskyrmion can still be metastable even in the presence of a positive density coupling.Negative γ has also a detrimental effect on surface vortex states, but with a different mechanism.The direct interaction leads to a reduction of the escape barrier for the fractional vortex such that, already for γ = −0.2,we observe spontaneous decay to a Meissner state.

IV. CONCLUSIONS
In this paper, we studied the minimum energy path of the vortex nucleation process for the Ginzburg-Landau model of two-component superconductors.We obtained solutions that feature multiple sphalerons connected to the entry of a single fractional vortex.Even composite vortices enter the sample through fractionalized nucleation in many systems.In the type-2 and type-1.5 regime, in agreement with the previous studies based on the London-model [17], fractional vortices with winding in the large-core component can be stabilized by the interaction with boundary currents.We characterized the stability of these solutions and the conditions for their appearance in a fully nonlinear theory.
Even when surface fractional vortices do not appear in the ground state of the system, but only as metastable solutions, this does not preclude the observability of such states, as their lifetime is, in general, not zero.Moreover, these metastable states can be stabilized by additional effects.We found, in particular, that mass anisotropy provides the strongest stabilization mechanism.
The method that we developed allowed for the inclusion and analysis of the effect of direct interaction between the fields like the bilinear Josephson coupling and density coupling.We find that fractional vortices are metastable near boundaries even for nonzero Josephson coupling.Note that, in some of the considered cases, the barrier for the nucleation process of the vortex in the second component is very low.For this reason, we do not expect the surface vortex states to be generically observed if not actively searched in in well-controlled experiments.In type-1.5 superconductors, vortex clusters with a total flux that is not an integer multiple of the flux quantum may appear as metastable solitons near boundaries.
These surface-trapped fractional objects can be observed in scanning probes and can serve as smoking-gun evidence of multicomponent order parameters in multiband superconductors, superconductors that break multiple symmetries, including models proposed for superconductivity in magic-angle twisted bilayer graphene [29,48].Applied currents can move these solitons along boundaries of samples with a high-quality surface, making them distinguishable from ordinary vortices pinned to a defect.
The formation of solitons with fractional magnetic flux can also be affected by the surface effects discussed in [49].This opens up the possibility of the creation of a new platform for fluxonic information processing [50].
As a secondary result, we have shown how the gauged string method can be applied to a system with additional gauge symmetries, like the U (1)×U (1) model, confirming the flexibility of this method.This can be easily extended to systems with other gauge symmetries or used to study other processes like vortex clustering in type-1.5 systems.

FIG. 1 .
FIG.1.Example of a minimum energy path of the nucleation of a fractional vortex with the endpoint in the center of the superconductor domain.The picture shows the definition of nucleation barrier ∆Fn, escape barrier ∆Fe, and bonding energy ∆F b .In case of absence of surface bonded vortex state ∆F b = 0. To pin the vortex in the center of the domain the origin (x, y) = (0, 0) has been removed from the superconductive domain.This is necessary since the final configuration is not a minimum of the system.With this trick, the vortex is within a local minimum with negligible effects on the total energy.

FIG. 2 .FIG. 3 .
FIG.2.Minimum energy path of the nucleation process for a composite vortex in System U with H = 1.0.The minimum energy path shows an example of fractionalized nucleation: the nucleation process happens in two steps, with the smallcore vortex nucleating first, followed by the large-core one.To report the nucleation steps, we use the following scheme: the dotted line signals that the systems has zero global winding number in both the components (N1, N2) = (0, 0), dashdotted line means that only one component feature a nonzero winding number (N1, N2) = (0, 1), finally continuous line means both the components have a non-zero winding (N1, N2) = (1, 1).The insets show the matter fields configuration in three points along the minimum energy path.

FIG. 4 .
FIG. 4. Stabilization of the surface fractional bonded state in System U, as the external magnetic field is increased from H = 0.6 (blue line) to H = 1.45 (red line).The solid lines represent the minimum energy path of fractional vortex nucleation.The dashed lines, representing the free energy of the system with the fractional vortex pinned in the center of the system, are included as a visual aid.F M indicates the free energy of the Meissner state.

FIG. 5 .
FIG.5.Free energy barriers for a fractional vortex in System U as a function of the external field H.At H = 0.7, the surface fractional vortex states become metastable, as the escape barrier becomes nonzero.Fractional surface-bond vortices can exist in the region between H = 0.7 and H = 1.4,while higher fields cause the nucleation of composite vortices.In this zone, the bonding energy is constant while the escape barrier increase with the magnetic field intensity.Notice how, for fields higher than H = 1.2, the escape barrier exceeds the nucleation one, meaning that a fractional vortex is more stable than the Meissner state.The bonding energy is approximately constant with respect to the magnetic field.

FIG. 7 .
FIG.7.Minimum energy path of a composite vortex nucleation in a two-component superconductor with spatial anisotropy in one of the two components.The simulations show nucleation processes occurring in System E with various mass ratios at H = 1.0.In these simulations, we study the vortex entry from a surface aligned with the y direction.In the case where the heavy direction is parallel to the surface, the bond state energy is enhanced by the vortex core deformation.To report the nucleation steps, we use the following scheme described in Fig.2

1 FIG. 8 .
FIG. 8. Minimum energy path of the first vortex nucleation compared with the nucleation of a second composite and fractional vortex in System H at H = 1.0.The dashed black line indicates the energy of a composite vortex.The plots in the second row show the first sphaleron configuration corresponding nucleation of the (1, 0) fractional vortex.In the third row, we show the configuration of the surface bonded fractional (2, 1) soliton.

1 FIG. 9 .
FIG. 9. Minimum energy path of the nucleation of a (2, 1) soliton in System D at H = 1.2.The top panel reports the minimum energy path of the nucleation process.The dashed lines N1(s) and N2(s) are included to help identify the nucleation of the fractional vortices.The contour plots report the system's order parameters, ψ1 and ψ2, in specified configurations along the minimum energy path marked by the a,b,c-points in the top panel.The nucleation is fractionalized in three steps with (1, 0) and (1, 1) as intermediate metastable states.The global minimum is achieved when (N1, N2) = (2, 1).

FIG. 11 .
FIG. 11.Effect of bilinear Josephson coupling on the nucleation barrier ∆Fn for composite vortex nucleation in System E with various bilinear Josephson coupling strength η as a function of the external field H.The continuous lines are the estimated exponential regressions.(b) Effect of the bilinear Josephson coupling on the nucleation field Hn in System E for composite vortex nucleation as a function of the bilinear Josephson coupling η calculated by exponential fitting.The nucleation field increases approximately linearly with η.

FIG. 13 .
FIG.13.Nucleation process of a skyrmion in System H with γ = 0.3 at H = 0.6.The intercomponent direct interaction causes the repulsion between the two fractional vortices in the two components, favoring fractionalized entry.The blue dot on the top panel identifies the transition's sphaleron, whose order parameter configurations, as well as the pseudospin vector field, are displayed in the bottom panels.

TABLE I .
Parameters used in the simulations.Since the model has many independent parameters, for each series of simulations we fix all the parameters as in this table except for one or a pair which are mentioned.

TABLE II .
The table reports the characteristic parameters of the system described in the paper.Component coherence length ξi and magnetic penetration depth λ, equilibrium density ψi, thermodynamic critical field Hc.The first critical field Hc1 is defined as the external field that guarantee