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

Andrea Maiani  Center for Quantum Devices, Niels Bohr Institute, University of
Copenhagen, DK-2100 Copenhagen, Denmark  Andrea Benfenati  Department of
Physics, The Royal Institute of Technology, Stockholm SE-10691, Sweden  Egor
Babaev  Department of Physics, The Royal Institute of Technology, Stockholm
SE-10691, Sweden

###### Abstract

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 inter-component 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 infinite-dimensional 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 [27, 18, 28, 29, 30, 31, 32, 33]. Here, we consider a two-
component model in which the order parameter $\Psi$ is represented by two
complex scalar fields $(\psi_{1},\psi_{2})$. The two component fields are
coupled by some direct interaction $V_{\mathrm{int}}$, and to the
electromagnetic field described by the gauge potential $\mathbf{A}$. The
system has a free energy
$F\quantity[\mathbf{A},\psi_{1},\psi_{2}]=\int\mathcal{F}\differential[3]{\mathbf{r}}$
where the free energy density functional reads, in natural units,

$\begin{split}\mathcal{F}[\mathbf{A},\psi_{1},\psi_{2}]&=\sum_{\alpha=1,2}\quantity[\sum_{k=x,y}\frac{\absolutevalue{D_{k\alpha}\psi_{\alpha}}^{2}}{2m_{\alpha
k}}+\frac{b_{\alpha}}{2}\quantity(\frac{a_{\alpha}}{b_{\alpha}}+\absolutevalue{\psi_{\alpha}}^{2})^{2}]+\frac{1}{2}\quantity(\curl\mathbf{A}-\mathbf{H})^{2}+V_{\mathrm{int}}\quantity(\psi_{1},\psi_{2})\,,\\\
V_{\mathrm{int}}\quantity(\psi_{1},\psi_{2})&=\frac{\eta}{2}\quantity(\psi_{1}\psi_{2}^{*}+\mathrm{c.c.})+\frac{\gamma}{2}\absolutevalue{\psi_{1}}^{2}\absolutevalue{\psi_{2}}^{2}\,.\end{split}$
(1)

Here the index $k$ labels the axes of the coordinate system and the index
$\alpha$ denotes the superconducting components, the vector field $\mathbf{A}$
is the gauge field while
$\psi_{\alpha}=\absolutevalue{\psi_{\alpha}}e^{i\theta_{\alpha}}$ are the two
matter fields, and $\mathbf{H}$ is the external magnetic field. The parameters
$m_{\alpha k}$ are the masses of the two fields in the reference system.  The
operator $D_{k\alpha}=-i\partial_{k}-q_{\alpha}A_{k}$ is the covariant
derivative where $q_{\alpha}$ 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 $\eta$, and the density coupling $\gamma$. 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:

$N_{i}=\int_{\partial\Omega}\gradient\theta_{i}\cdot\differential{l}.$ (2)

If we neglect the direct interaction term
$V_{\mathrm{int}}\quantity(\psi_{1},\psi_{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)\crossproduct U(1)$ gauge symmetry.

The magnetic flux carried by a vortex line, in the model that we consider, is
$\Phi=\sum_{\alpha}N_{\alpha}\Phi_{\alpha}$ with $\Phi_{\alpha}=\frac{2\pi
q_{\alpha}\bar{\psi}_{\alpha}^{2}}{\sum_{\beta}q_{\beta}^{2}\bar{\psi}_{\beta}^{2}}$
where $\overline{\psi}_{\alpha}=\sqrt{\frac{-a_{\alpha}}{b_{\alpha}}}$ is the
bulk equilibrium value of each order parameter, and $N_{\alpha}$ is an
integer.

In conventional superconductors, $q_{1}=q_{2}=q$ and therefore it simplifies
to

$\Phi=\Phi_{0}\left[N_{1}\frac{\overline{\psi}^{2}_{1}}{\overline{\psi}^{2}_{1}+\overline{\psi}^{2}_{2}}+N_{2}\frac{\overline{\psi}^{2}_{2}}{\overline{\psi}^{2}_{1}+\overline{\psi}^{2}_{2}}\right]\,,$
(3)

where $\Phi_{0}=2\pi/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 $\overline{\psi}_{1}=\overline{\psi}_{2}$, the flux carried by a
vortex is

$\Phi=\Phi_{0}\left[\frac{g_{1}N_{1}+g_{2}N_{2}}{g_{1}^{2}+g_{2}^{2}}\right]\,.$
(4)

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.

The interaction of vortices between themselves and the boundary are determined
by the different length scales that are present in the model. In a isotropic
system without direct interaction, i.e. $m_{\alpha x}=m_{\alpha y}$ and
$V_{\mathrm{int}}=0$, we can identify the magnetic field penetration depth as
$\lambda=\quantity(\sum_{\alpha}\lambda_{\alpha}^{-2})^{-1/2}$ where
$\lambda_{\alpha}=\sqrt{\frac{m_{\alpha}}{q_{\alpha}^{2}\overline{\psi}_{\alpha}^{2}}}$,
and the components coherence lengths as
$\xi_{\alpha}=\frac{1}{2\sqrt{-m_{\alpha}a_{\alpha}}}$. If
$\xi_{1},\xi_{2}<\lambda$ the system behave as a type-2 superconductor, while
if $\xi_{1},\xi_{2}>\lambda$ it behaves as a type-1. Finally, the regime where
$\min(\xi_{1},\xi_{2})<\lambda<\max(\xi_{1},\xi_{2})$ is termed type-1.5
superconductivity [34, 35, 32, 36, 37, 6, 38, 39]. In this regime, the
interaction between vortices are long-range attractive and short-range
repulsive. When an anisotropic mass tensor is introduced the length scales
become direction-dependent increasing the number of possible regimes [40].

In models with direct intercomponent interactions, i.e. $V_{\mathrm{int}}\neq
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 [41, 39, 35]. 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
$\mathbf{S}=\frac{\Psi^{\dagger}{\bf\it\sigma}\Psi}{\Psi^{\dagger}\Psi}$ where
$\Psi=(\psi_{1},\psi_{2})^{T}$ where and ${\bf\it\sigma}$ are the Pauli
matrices vector [42, 18].

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
$\mathbf{q}(s)=\quantity(\mathbf{A}(s),\,\psi_{1}(s),\,\psi_{2}(s))$, where
the transition coordinate $s\in[0,1]$ has been introduced. A minimum energy
path can be defined in a variational formulation as

$\mathbf{q}_{\mathrm{MEP}}(s)=\arg\min F_{\bot}[\mathbf{q}(s)]\,,\quad\forall
s\,\in[0,1]\,,$ (5)

where $F_{\bot}$ is that the free energy functional $F$ restricted to the
cotangent space of the trajectory in the point $\mathbf{q}_{\mathrm{MEP}}(s)$
[43]. Once the minimum energy path $\mathbf{q}_{\mathrm{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
$\mathbf{B}=\curl\mathbf{A}$, and currents $j_{\alpha
k}=\psi_{\alpha}^{*}\frac{D_{k\alpha}}{2m_{\alpha
k}}\psi_{\alpha}+\mathrm{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_{\mathrm{S}}=\max_{s}F(s)$. We define $F_{\mathrm{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 $\Delta
F_{\mathrm{n}}=F_{\mathrm{S}}-F_{\mathrm{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_{\mathrm{C}}$. If no surface bond states are
present, the escape energy of the vortex is simply given by $\Delta
F_{\mathrm{e}}=F_{\mathrm{S}}-F_{\mathrm{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_{\mathrm{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_{\mathrm{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, $\Delta
F_{\mathrm{e}}=F_{\mathrm{S}}-F_{\mathrm{B}}$, while the bonding energy
$\Delta F_{\mathrm{b}}=F_{\mathrm{C}}-F_{\mathrm{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.

Figure 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 $\Delta F_{\mathrm{n}}$,
escape barrier $\Delta F_{\mathrm{e}}$, and bonding energy $\Delta
F_{\mathrm{b}}$. In case of absence of surface bonded vortex state $\Delta
F_{\mathrm{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.

For the simulations, we consider a system with computational domain
$\Omega=[-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 show the results for some representative systems.
For each series of simulations, we fix all parameters according to Tab. 1
except for one or a pair. First, we focus on the case where the intercomponent
direct interaction is absent, i.e. $V_{\mathrm{int}}=0$.

| $a_{1}$ | $b_{1}$ | $m_{1x}$ | $m_{1y}$ | $q_{1}$ | $a_{2}$ | $b_{2}$ | $m_{2x}$ | $m_{2y}$ | $q_{2}$
---|---|---|---|---|---|---|---|---|---|---
System E | -1 | 1 | 1 | 1 | -1 | -1 | 1 | 1 | 1 | -1
System U | -1 | 1 | 1 | 1 | -1 | -1 | 1 | 2.5 | 2.5 | -1
System H | -0.5 | 0.5 | 0.5 | 0.5 | -1 | -0.5 | 0.5 | 4 | 4 | -1
System D | -1 | 1 | 2 | 2 | -2 | -1 | 1 | 1 | 1 | -1
Table 1: 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.

The characteristic length scales for the systems considered are shown in Tab.
2. 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 $\psi_{1}$ as the large-
core component while $\psi_{2}$ is the small-core one.

| $\xi_{1}$ | $\xi_{2}$ | $\lambda$ | $\bar{\psi}_{1}$ | $\bar{\psi}_{2}$ | $H_{\mathrm{c}}$ | $H_{\mathrm{c1}}$
---|---|---|---|---|---|---|---
System E | 0.50 | 0.50 | 0.71 | 1.00 | 1.00 | 1.41 | 1.16
System U | 0.50 | 0.32 | 0.85 | 1.00 | 1.00 | 1.41 | 0.93
System H | 1.00 | 0.35 | 0.67 | 1.00 | 1.00 | 1.00 | -
System D | 0.35 | 0.50 | 0.58 | 1.00 | 1.00 | 1.41 | -
Table 2:  The table reports the characteristic parameters of the system
described in the paper. Component coherence length $\xi_{i}$ and magnetic
penetration depth $\lambda$, equilibrium density $\bar{\psi}_{i}$,
thermodynamic critical field $H_{\mathrm{c}}$. The first critical field
$H_{\mathrm{c1}}$ is defined as the external field that guarantee
$F_{M}(H_{\mathrm{c1}})=F_{C}(H_{\mathrm{c1}})$. Notice that the composite
vortex considered for System D has a $(2,1)$ structure.

###  III.1 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.

Figure 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 small-core 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
$(N_{1},N_{2})=(0,0)$, dash-dotted line means that only one component feature
a non-zero winding number $(N_{1},N_{2})=(0,1)$, finally continuous line means
both the components have a non-zero winding $(N_{1},N_{2})=(1,1)$. The insets
show the matter fields configuration in three points along the minimum energy
path.

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.

Figure 3: Free energy and winding numbers computed along the minimum energy
path of the nucleation process for a fractional vortex with winding in the
first component  in System U at $H=1.0$. If a vortex has non-zero winding only
in the large-core component, it remains bonded to the surface. The curves
flatten for $s\to 1$ because of the interaction of the vortex with edges in
the finite-size system. To report the nucleation steps, we use the following
scheme: the dotted line indicates that the systems’ winding numbers are
$(N_{1},N_{2})=(0,0)$, while the dash-dotted line means $(N_{1},N_{2})=(1,0)$.

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.

Figure 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_{\textrm{M}}$ indicates the free energy of
the Meissner state.

It is possible to analyze the phase space of the system by extracting the free
energy barriers $\Delta F_{\mathrm{n}}$, $\Delta F_{\mathrm{e}}$, $\Delta
F_{\mathrm{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 $\Delta F_{\mathrm{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].

Figure 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.

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.

Figure 6: Comparison of nucleation minimum energy paths of a composite vortex
($\mathrm{C}$), a fractional vortex ($\mathrm{F}$) and the nucleation of a
second component vortex starting from a fractional one
($\mathrm{F}\to\mathrm{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_{\textrm{C}}$ and the gray line is the energy of the
Meissner state with respect $F_{\textrm{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 dashed-dotted line signals the presence of a fractional
vortex, while the continuous line means that there is a composite vortex in
the system.

###  III.2 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.

Figure 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

###  III.3 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)\to(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)\to(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.

Figure 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.

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.

###  III.4 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.

Figure 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 $N_{1}(s)$ and $N_{2}(s)$ are included to help
identify the nucleation of the fractional vortices. The contour plots report
the system’s order parameters, $\psi_{1}$ and $\psi_{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
$(N_{1},N_{2})=(2,1)$.

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.

Figure 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 $\psi_{1}$ and $\psi_{2}$, and magnetic field $\mathbf{B}$ at
specified configurations in the minimum energy path marked by a blue dot in
the top panel.

###  III.5 Effects of direct intercomponent interactions

In this section, we consider the effect of direct interactions,  i.e.
$V_{\mathrm{int}}\neq 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

$V_{\mathrm{lin}}\quantity(\psi_{1},\psi_{2})=\frac{\eta}{2}\quantity(\psi_{1}\psi_{2}^{*}+\mathrm{c.c.})=\eta\absolutevalue{\psi_{i}}\absolutevalue{\psi_{2}}\cos(\theta_{12})\,.$
(6)

This term introduces phase-locking in the system such the spontaneously broken
symmetry is reduced to $U(1)$. For positive $\eta$, the minimum is obtained
for $\theta_{12}=\pm\pi$, while for negative values the minimum is for
$\theta_{12}=0$. In the following, we restrict the study to positive $\eta$ 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 $\Delta
F^{\textrm{(fit)}}_{\textrm{n}}(H)=C_{1}\quantity[1-\exp\quantity(-C_{2}(H-H_{\textrm{n}}))]$
as the ansatz, where $C_{1}$, $C_{2}$, and $H_{\textrm{n}}$ are the
coefficients estimated for each value of $\eta$. 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_{\mathrm{n}}(\eta)$ increase
linearly as a function of $\eta$. Note that the fractional vortex interaction
potential is asymptotically linear in the presence of this coupling.

Figure 11: (a) Effect of bilinear Josephson coupling on the nucleation barrier
$\Delta F_{\mathrm{n}}$ for composite vortex nucleation in System E with
various bilinear Josephson coupling strength $\eta$ 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 $H_{n}$ in System E for composite vortex nucleation as a function of the
bilinear Josephson coupling $\eta$ calculated by exponential fitting. The
nucleation field increases approximately linearly with $\eta$.

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 $\eta=0.01$, the fractional surface
vortex becomes a metastable state with energy higher than the Meissner state.
Increasing further $\eta$, 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

$V_{\mathrm{den}}\quantity(\psi_{1},\psi_{2})=\frac{\gamma}{2}\absolutevalue{\psi_{1}}^{2}\absolutevalue{\psi_{2}}^{2}\,,$
(7)

this term does not break the $U(1)\crossproduct U(1)$ gauge symmetry.
Depending on the sign of $\gamma$, this term promotes an attraction or
repulsion between fractional vortices in different components. A repulsive
(positive $\gamma$) 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 partial-skyrmion state.

Figure 12: (a) Effect of density coupling on the nucleation barrier $\Delta
F_{\mathrm{n}}$ of composite vortices in System E. The continuous lines are
exponential fitting of the nucleation barriers. Strong biquadratic Josephson
coupling changes the shape of the curve $H_{\mathrm{n}}(H)$ in a way that
reduces the quality of an exponential fitting. Panel (b) shows the composite
nucleation field $H_{n}$ as a function of the density coupling parameter
$\gamma$.

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 $\gamma$
increases. Moreover, the metastable minimum state is found for a vortex
further from the surface. For $\gamma\simeq 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_{\rm c1}$. Nonetheless, for some
type-2 systems, the partial-skyrmion can still be metastable even in the
presence of a positive density coupling.

Figure 13:  Nucleation process of a skyrmion in System H with $\gamma=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.

Negative $\gamma$ 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 $\gamma=-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 [48, 29]. 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)\times
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.

##  V Acknowledgments

We gratefully acknowledge the support of NVIDIA Corporation with the donation
of the Quadro P6000 GPU used for this research. A.M. was supported by the
Danish National Research Foundation, the Danish Council for Independent
Research |Natural Sciences. A.B. and E.B. were supported by the Swedish
Research Council Grants No. 2016-06122, 2018-03659, and Göran Gustafsson
Foundation for Research in Natural Sciences.

## References

  * Suhl  _et al._ [1959] H. Suhl, B. T. Matthias, and L. R. Walker, Phys. Rev. Lett. 3, 552 (1959).
  * Moskalenko [1959] V. A. Moskalenko, Fiz. Met. Met. 8, 2518 (1959).
  * Kamihara  _et al._ [2008] Y. Kamihara, T. Watanabe, M. Hirano, and H. Hosono, J. Am. Chem. Soc. 130, 3296 (2008).
  * Nagamatsu  _et al._ [2001] J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani, and J. Akimitsu, Nature 410, 63 (2001).
  * Kreisel  _et al._ [2020] A. Kreisel, P. J. Hirschfeld, and B. M. Andersen, Symmetry 12, 1402 (2020).
  * Svistunov  _et al._ [2015] B. V. Svistunov, E. S. Babaev, N. V. Prokof’ev, L. Bruno, and N. V. Prokof’ev,  _Superfluid States of Matter_ (CRC Press, 2015) pp. 1–546.
  * Grinenko  _et al._ [2017] V. Grinenko, P. Materne, R. Sarkar, H. Luetkens, K. Kihou, C. H. Lee, S. Akhmadaliev, D. V. Efremov, S.-L. Drechsler, and H.-H. Klauss, Phys. Rev. B 95, 214511 (2017).
  * Grinenko  _et al._ [2018] V. Grinenko, R. Sarkar, K. Kihou, C. H. Lee, I. Morozov, S. Aswartham, B. Büchner, P. Chekhonin, W. Skrotzki, K. Nenkov, and Others, arXiv Prepr. arXiv1809.03610 (2018).
  * Grinenko  _et al._ [2021] V. Grinenko, D. Weston, F. Caglieris, C. Wuttke, C. Hess, T. Gottschall, J. Wosnitza, A. Rydh, K. Kihou, C.-H. Lee, R. Sarkar, S. Dengre, I. Maccari, J. Garaud, A. Charnukha, R. Hühne, K. Nielsch, B. Büchner, H.-H. Klauss, and E. Babaev, Bosonic $Z_2$ metal: Spontaneous breaking of time-reversal symmetry due to Cooper pairing in the resistive state of Ba$_{1-x}$K$_x$Fe$_2$As$_2$ (2021), arXiv:2103.17190 [cond-mat.supr-con] .
  * Kivelson  _et al._ [2020] S. A. Kivelson, A. C. Yuan, B. Ramshaw, and R. Thomale, npj Quantum Mater. 5, 1 (2020).
  * Sigrist and Ueda [1991] M. Sigrist and K. Ueda, Rev. Mod. Phys. 63, 239 (1991).
  * Babaev [2002] E. Babaev, Phys. Rev. Lett. 89, 4 (2002), arXiv:0111192 .
  * Bean and Livingston [1964] C. P. Bean and J. D. Livingston, Phys. Rev. Lett. 12, 14 (1964).
  * Manton and Sutcliffe [2004] N. Manton and P. Sutcliffe,  _Topological solitons_ (Cambridge University Press, 2004) p. 403.
  * Rubakov [2009] V. Rubakov,  _Class. Theory Gauge Fields_ (Princeton University Press, 2009) pp. 1–444.
  * Năstase [2019] H. Năstase,  _Class. F. Theory_ (Cambridge University Press, 2019) pp. 1–462.
  * Silaev [2011] M. A. Silaev, Phys. Rev. B - Condens. Matter Mater. Phys. 83, 144519 (2011).
  * Agterberg  _et al._ [2014] D. F. Agterberg, E. Babaev, and J. Garaud, Phys. Rev. B - Condens. Matter Mater. Phys. 90, 64509 (2014), arXiv:1403.6655 .
  * de Gennes [1965] P. G. de Gennes, Solid State Commun. 3, 127 (1965).
  * Catelani and Sethna [2008] G. Catelani and J. P. Sethna, Phys. Rev. B 78, 224509 (2008).
  * Transtrum  _et al._ [2011] M. K. Transtrum, G. Catelani, and J. P. Sethna, Phys. Rev. B 83, 94505 (2011).
  * Liarte  _et al._ [2017] D. B. Liarte, S. Posen, M. K. Transtrum, G. Catelani, M. Liepe, and J. P. Sethna, Supercond. Sci. Technol. 30, 033002 (2017), arXiv:1608.00175 .
  * Pack  _et al._ [2020] A. R. Pack, J. Carlson, S. Wadsworth, and M. K. Transtrum, Phys. Rev. B 101, 144504 (2020), arXiv:1911.02132 .
  * E  _et al._ [2002] W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B - Condens. Matter Mater. Phys. 66, 523011 (2002), arXiv:0205527 [cond-mat] .
  * E and Vanden-Eijnden [2010] W. E and E. Vanden-Eijnden, Annu. Rev. Phys. Chem. 61, 391 (2010).
  * Benfenati  _et al._ [2020] A. Benfenati, A. Maiani, F. N. Rybakov, and E. Babaev, Phys. Rev. B 101, 220505 (2020), arXiv:1911.09513 .
  * Tilley [1964] D. Tilley, Proceedings of the Physical Society (1958-1967) 84, 573 (1964).
  * Wagner  _et al._ [2021] G. Wagner, H. S. Røising, F. Flicker, and S. H. Simon, Physical Review B 104, 134506 (2021).
  * Chichinadze  _et al._ [2020] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Phys. Rev. B 101, 224513 (2020).
  * Maiti and Chubukov [2013] S. Maiti and A. V. Chubukov, Physical Review B 87, 144511 (2013).
  * Garaud  _et al._ [2017] J. Garaud, M. Silaev, and E. Babaev, Physica C: Superconductivity and its Applications 533, 63 (2017).
  * Silaev and Babaev [2012] M. Silaev and E. Babaev, Phys. Rev. B 85, 134514 (2012).
  * Krohg and Sudbø [2018] F. N. Krohg and A. Sudbø, Phys. Rev. B 98, 014510 (2018).
  * Babaev and Speight [2005] E. Babaev and M. Speight, Phys. Rev. B - Condens. Matter Mater. Phys. 72, 180502 (2005), arXiv:0411681 .
  * Silaev and Babaev [2011] M. Silaev and E. Babaev, Phys. Rev. B 84, 094515 (2011), arXiv:1102.5734 .
  * Moshchalkov  _et al._ [2009] V. Moshchalkov, M. Menghini, T. Nishio, Q. H. Chen, A. V. Silhanek, V. H. Dao, L. F. Chibotaru, N. D. Zhigadlo, and J. Karpinski, Phys. Rev. Lett. 102, 1 (2009), arXiv:0902.0997 .
  * Babaev  _et al._ [2017] E. Babaev, J. Carlström, M. Silaev, and J. M. Speight, Phys. C Supercond. its Appl. 533, 20 (2017), arXiv:1608.02211 .
  * Babaev and Carlström [2010] E. Babaev and J. Carlström, Phys. C Supercond. its Appl. 470, 717 (2010).
  * Carlström  _et al._ [2011] J. Carlström, E. Babaev, and M. Speight, Phys. Rev. B 83, 174509 (2011), arXiv:1009.2196 .
  * Winyard  _et al._ [2019] T. Winyard, M. Silaev, and E. Babaev, Physical Review B 99, 064509 (2019).
  * Babaev  _et al._ [2010] E. Babaev, J. Carlström, and M. Speight, Physical review letters 105, 067003 (2010).
  * Babaev  _et al._ [2002] E. Babaev, L. D. Faddeev, and A. J. Niemi, Phys. Rev. B65, 100512 (2002).
  * Samanta and Weinan [2013] A. Samanta and E. Weinan, Commun. Comput. Phys. 14, 265 (2013).
  * Babaev  _et al._ [2004] E. Babaev, A. Sudbø, and N. Ashcroft, Nature 431, 666 (2004).
  * Garaud and Babaev [2014] J. Garaud and E. Babaev, Phys. Rev. B - Condens. Matter Mater. Phys. 89, 1 (2014).
  * Chatterjee  _et al._ [2020] C. Chatterjee, S. B. Gudnason, and M. Nitta, Journal of High Energy Physics 2020, 1 (2020).
  * Garaud  _et al._ [2014] J. Garaud, K. A. H. Sellin, J. Jäykkä, and E. Babaev, Phys. Rev. B 89, 104508 (2014).
  * Khalaf  _et al._ [2020] E. Khalaf, P. Ledwith, and A. Vishwanath, arXiv Prepr. arXiv2012.05915 (2020).
  * Samoilenka and Babaev [2020] A. Samoilenka and E. Babaev, Phys. Rev. B 101, 134512 (2020).
  * Golod  _et al._ [2015] T. Golod, A. Iovan, and V. M. Krasnov, Nature communications 6, 1 (2015).