Properties of skyrmions and multi-quanta vortices in chiral p-wave superconductors

Chiral p-wave superconducting state supports a rich spectrum of topological excitations different from those in conventional superconducting states. Besides domain walls separating different chiral states, chiral p-wave state supports both singular and coreless vortices also interpreted as skyrmions. Here, we present a numerical study of the energetic properties of isolated singular and coreless vortex states as functions of anisotropy and magnetic field penetration length. In a given chiral state, single quantum vortices with opposite winding have different energies and thus only one kind is energetically favoured. We find that with the appropriate sign of the phase winding, two-quanta (coreless) vortices are always energetically preferred over two isolated single quanta (singular) vortices. We also report solutions carrying more flux quanta. However those are typically more energetically expensive/metastable as compared to those carrying two flux quanta.

the dominant component. It has a 4π winding of the relative phase between components, that follows from the Cooper pairs having nonzero internal orbital momentum 28 . Since the magnetic field lifts the degeneracy between chiralities, vortices and anti-vortices have different properties 29 .
It is experimentally seen that in an applied external field, Sr 2 RuO 4 exhibits vortex lattices with square symmetry at high fields [30][31][32] , and a transition to triangular lattice in lower fields 32,33 . Earlier theoretical calculations based on Ginzburg-Landau model for chiral p-wave superconductivity in Sr 2 RuO 4 [34][35][36] , are consistent with these observed transitions of the vortex lattice structure.
Besides single-quanta vortices, there also exists vortices carrying multiple quanta of the magnetic flux and that, as they are coreless, are essentially different from single-quanta vortices. For example as discussed in more details below, the component induced by a doubly quantized vortex in the dominant component has zero winding in subdominant one 29 . In this paper we demonstrate that the two-quanta (coreless) vortices, which can also be denoted as skyrmions, are energetically favoured as compared to (isolated) single-quanta vortices. Earlier works in the context of UPt 3 , even claim that lattices of similar two-quanta vortices may be energetically favoured as compared to those of single quanta 37,38 . The possible existence of lattices of different coreless vortices carrying single flux quantum in UPt 3 was also discussed recently 14 . It was also recently shown in the context of Sr 2 RuO 4 , based on solutions of microscopic Eilenberger equations, that lattices of two-quanta vortices are favoured for certain parameter sets 39 . Yet, such lattices of two-quanta vortices were never observed in Sr 2 RuO 4 . This motivates this work to further investigate the thermodynamic stability of skyrmions for broad parameter range.
In a previous work 40 , we reported isolated skyrmion solutions in a model for chiral p-wave superconductor. For the studied case of one of the chiralities, skyrmions can be energetically favoured as compared to vortices [Note that the Ginzburg-Landau model which was used in ref. 40 had slightly different coefficients in the potential terms compared to the standard GL model which follows from the weak-coupling mean-field theory. In this paper we use the same model as in ref. 36]. The skyrmions carrying two flux quanta are directly related to the two-quanta vortices mentioned above. However it was also demonstrated in ref. 40 that there are (meta-)stable skyrmions carrying larger number of flux quanta. In this model the energy and structure of vortices and skyrmions depends on the chirality. Equivalently, for a given chiral state, vortex/skyrmion solutions are not the same as anti-vortex/ anti-skyrmion. It thus calls for further investigation of vortex and skyrmion solutions (carrying two and more quanta), which we present below.

Model
In the coordinate system where the crystal anisotropy axis is c z, the order parameter of the p x + ip y state is described by a complex two-dimensional vector η = (η x , η y ) 10,11,41 . Introducing the chiral order parameter components ( ) , the dimensionless Ginzburg-Landau free energy reads as (see e.g. [34][35][36]: There we use dimensionless units were the free energy is normalized to the condensation energy, and the lengths are given in units of . In these units, the gauge coupling g that appears in the covariant derivative D = ∇ + igA is related to the Ginzburg-Landau parameter g −1 = λ/ξ. The free energy (1) was derived from the weak coupling microscopic theory 34,35 . The anisotropy parameter v determines the anisotropy in the xy-plane ( ν | |<1 for the energy to be positively defined). It measures the tetragonal distortions of the Fermi surface, which has cylindrical geometry for v = 0, and is defined as (where 〈 ⋅ 〉 denote average over the Fermi surface). In the model Eq. (1), the dependence on the third coordinate is not considered (i.e. assuming two-dimensional system or translational invariance along z-axis).
The ground-state that minimizes the potential terms in (1) is degenerate and the solutions are (η + , η − ) = (1, 0) and (0, 1). The theory (1) is invariant under the (discrete) time-reversal operations  , as This invariance is spontaneously broken by the ground-state. The spontaneous breakdown of the discrete time-reversal symmetry dictates that the theory allows domain wall solutions that interpolate between regions with different ground-states. Such domain walls, rather generically created during phase transition where the discrete symmetry is broken 21 , carry a magnetic field perpendicular to the xy-plane 17,18 . The discrete degeneracy of the ground state is lifted by the magnetic field. Thus, depending on the direction of the external field, only one state is stable. Likewise, the vorticity of the superconducting condensates lifts the degeneracy between chiral (ground-)states.
Scientific RepoRts | 5:17540 | DOI: 10.1038/srep17540 As the components η + and η − behave differently for different sign of the winding, a complete study requires to consider both situations of counter-clockwise (positive) and clockwise (negative) vorticities. Note that this is equivalent to considering only positive vorticity but for both chiral states. For example, the configuration with a winding n + = + 1 on the ground-state (η + , η − ) = (1, 0) can be obtained by applying the time-reversal operation on the configuration whose ground state is (η + , η − ) = (0, 1) with the winding n − = − 1. In the following, we choose to fix the dominant component of the order parameter to be η − [i.e. the ground state is (η + , η − ) = (0, 1)] and thus need to investigate both positive and negative vorticities.
The asymptotic vorticity of the dominant component η − determines the sign of B z , as well as the vorticity of the subdominant component η + where θ is the polar angle. The relative phase between η + and η − (2), follows from the internal orbital momentum of Cooper pairs. In the Ginzburg-Landau model (1), this is the structure of mixed gradient that constraints the relative phase. Note that since the subdominant component η + vanishes asymptotically [i.e. it recovers its ground state value η + = 0], the winding n + can be located only in a close vicinity of the vortex core. Note also that the winding of the subdominant component does not affect the overall flux quantization, because the density of that component vanishes away from the vortex. From (2) it is rather straightforward to see that the vortex with the vorticity (n + , n − ) = (+ 3, + 1) and the (anti-)vortex with (n + , n − ) = (+ 1, − 1) have different core structures and thus different energy.

Results
In order to investigate the energetic properties of vortex matter, the fields η ± and A are discretized using a finite-element framework 42 and the free energy (1) is minimized using an nonlinear conjugate gradient algorithm (see the section Methods for details). In simulations of chiral p-wave superconductors on a finite domain, a special attention is required for boundary conditions in order to yield edge currents (see for example discussions in 19,41,43 ). Here, we are interested in the intrinsic energetic properties of isolated defects. Thus vortex configuration is created by an initial guess and placed them at the center of a large domain, with open boundary conditions, letting the fields freely reach the ground-state. By choosing a sufficiently large domain, this ensures that within numerical accuracy vortices will not interact with boundaries and thus we are able to probe their intrinsic structure and energy properties, without effects of boundary conditions. As it specifies the topological sector, a starting configuration with a given winding n − of the dominant component η − leads, after convergence of the algorithm, to a configuration that behaves as expected from (2). We systematically construct vortex solutions carrying one to four flux quanta for parameter space defined by wide range of values of the g and v. Figure 1 shows typical vortex solutions with different vorticities. Along this paper we also refer to vortices carrying multiple flux quanta, as skyrmions. The reason for that terminology is that they have additional topological properties, as compared to single quanta vortices. This is explained in more details by the end of the paper. The first and second blocks in Fig. 1 respectively show vortex solutions with B z > 0 and B z < 0. Vortices carrying from one to four flux quanta are displayed within each block. As expected from Eq. (2), single winding of the dominant component induces core structure of the subdominant component with different winding depending on that of η − (see the first column of each block). It is instructive to consider the last row in Fig. 1, that displays the relative phase between η − and η + . In agreement with (2), asymptotically the relative phase ϕ − − ϕ + = − 2θ, reflecting the orbital angular momentum difference between η − and η + . Moreover, the relative phase also indicates the position of the singularities of the components of the order parameter. Remarkably single quanta vortices are singular defects, since singularities in both components overlap (and thus η + = η − = 0). On the other hand, since both components never simultaneously vanish, two-quanta vortices are coreless defects. Interestingly the n − = − 2 configuration features a π jump of the relative phase when going outward from the vortex core. Inside the vortex core the time-reversed chiral state (η + , η − ) = (1, 0) is induced, while the (η + , η − ) = (0, 1) state is recovered asymptotically. The two quanta vortices thus feature a ringlike domain wall when going away from the center. The π jump of the relative phase for the n − = − 2 is located at this domain wall.
Like in conventional superconductors, the magnetic field for single quanta vortices is localized at the vortex core and screened at length scales determined by the penetration depth λ. Interestingly, the magnetic field for two-quanta vortices, and especially for n − = − 2, is not homogeneously distributed in the core. Rather it is localized at a given distance from the center and spread along the ring of the domain wall. Note that similar vortex configurations were also found to exist in the context of two-component model with  ( ) × ( ) × U 1 U 1 2 symmetry 44 . The ring-like distribution of the magnetic field for the two-quantum vortex can be understood as follows: B outside the vortex is screened by the (partial) currents in η − that run counter-clockwise, while inside the vortex currents in η + are responsible for the screening. Since η + vanishes away from the vortex core, it cannot contribute to the screening asymptotically. Conversely, η − vanishes at the vortex core and this is the induced subdominant component η + that screens B close to the center of the vortex core. The reason it can contribute to screening (inside the vortex) without having vorticity on its own, is only due to supercurrent produced by the vector potential (like the Meissner currents on the boundary of ordinary superconductors). Since those currents circulate clockwise, they compensate with the currents in η − so that at a certain distance (at the domain wall) there is no screening current. The magnetic field is thus localized at the domain wall. Although the core structure of single-quanta vortices are different depending on the sign of n − , their profile of the magnetic field looks quite similar. When considering vortices with > − n 1, both the core structure and the magnetic field profile are strikingly different and the skyrmions with negative n − do not resemble those with positive n − . Apart from the n − = − 2, − 3 skyrmions, the configurations that carry multiple flux quanta are far from being axially symmetric. Note that the n − = − 4 skyrmions resembles as some kind of bound state of two n − = − 2 skyrmions. As we will discuss below, this makes their decay into two n − = − 2 vortices rather easy.
Since the core structure is different, it is quite natural to expect that, unlike in conventional superconductors, vortices with opposite winding (and thus opposite directions of the magnetic field) are non-degenerate in energy. The (n + , n − ) = (+ 3, + 1) vortex has more total vorticity than the (n + , n − ) = (+ 1, − 1). Thus one could naively expect that the n − = − 1 vortices would be favoured as  Fig. 2(a), shows that the ratio of the energies of the single quanta vortices with n − = − 1 and n − = + 1, is always less than one. This implies that vortices n − = − 1 are always energetically favoured, as compared to those with n − = + 1. The first critical field of a vortex carrying a flux Φ is H c1 = E/2Φ , where E is its energy. As a result, Fig. 2(a) also implies that the n − = − 1 vortices also have lower first critical field H c1 in agreement with refs 36,45. Although both n − = ± 1 are perturbatively stable (i.e. they are minima of  ), only n − = − 1 is absolutely stable.
Note that the naive estimates based on counting the total vorticity provide the correct picture that (n + , n − ) = (+ 1, − 1) vortices are less energetic than (n + , n − ) = (+ 3, + 1) ones. It thus makes sense to apply the same arguments to configurations carrying more than one flux quantum. In the sector with negative n − , there are two possibilities to make a configuration that carries two flux quanta. Either to create two isolated (n + , n − ) = (+ 1, − 1) vortices carrying one flux quantum each or to create one (n + , n − ) = (0, − 2) two-quantum vortex. It turns out that a two-quantum vortex with smaller number of singularities is favoured as compared to two isolated single-quanta vortices. Figure 2(b) displays the ratio of the energies of two (isolated) n − = − 1 vortices and one n − = − 2 vortex. This ratio is always larger than one, thus implying that two-quanta vortices are energetically favoured as compared to two isolated single-quanta vortices. Note that the quantity displayed in Fig. 2(b), is actually also the ratio of first critical fields associated with single and double quanta vortices H c1 (n − = − 1)/H cl (n − = − 2). Note also that smaller H c1 for a higher-flux vortex does not necessarily imply that such vortices will nucleate first in low magnetic field. That is, due to higher winding they carry larger magnetic flux and thus can have a higher potential barrier to enter the sample (compared with the discussion of Bean-Livingston barrier in single component superconductors 46 ). The vortices (n − = − 1) and (n − = − 2) should interact differently with the Meissner currents and image charges, and thus even if the (n − = − 2) vortices have lower H cl , the interaction with the boundary may instead favour the entry of the vortices with (n − = − 1).
We also calculated the energy diagram similar to that in Fig. 2(b), but for vortices carrying three flux quanta n − = − 3 (data not shown). We found that unlike for n − = − 2, the n − = − 3 are not always stable. That is, in some regions of the parameter space the n − = − 3 is found, but in some other regions it decays into one single-quantum plus one double-quantum vortex. We find that for n − < 0, the n − = − 2 skyrmions are all energetically favoured. This behaviour can already be anticipated from the last column in Fig. 1 where the four quanta n − = − 4 skyrmion seems to be a bound state of two n − = − 2 vortices. As the skyrmions with n − < − 2 are more energetic than those with n − = − 2, one can easily see that the n − = − 4 configuration can decay into two n − = − 2 vortices and thus reduce its total energy. We find that in the regions where the n − = − 3 vortices exist, they are more energetic than three isolated single quanta vortices (or one single plus one two-quantum vortex).
Although being energetically unfavoured, it is still instructive to consider the properties of vortices carrying multiple flux quanta, with n − > 0. Diagrams in Fig. 3 show the ratio of the energies of multiple quanta vortices with n − > 0, compared to that of isolated vortices carrying smaller flux. The situation for n − > 0 is actually very different from that with n − < 0. Panels (a) and (b) in Fig. 3 display the ratio of the energies of isolated single-quanta vortices with that of vortices carrying two and three flux quanta. Depending on the anisotropy parameter v and on the gauge coupling g this ratio can either be smaller or larger than one and the solid lines on the diagram show when these are degenerate in energy. Below the solid line, isolated single-quanta vortices are energetically favoured as compared to multi-quanta vortices. Above this line, these are the vortices carrying two or three flux quanta which are favoured. Thus tetragonal distortions of the Fermi surface (i.e. larger ν | |) tend to favour n − = + 2 (and to a lesser  Figure 4 shows pseudo-spin texture for vortex solutions corresponding to those displayed in Fig. 1. Note that n is ill-defined for singular vortices, since there η =  0 at the core (i.e. singularities in both components overlap). Coreless vortices on the other hand have well defined pseudo-spin projection which is a map → n S : R 2 2 . Since at spatial infinity, n = (0, 0, − 1), the plane R 2 can be compactified to S 2 so that the pseudo-spin becomes a map → n S S : 2 2 . The homotopy invariants  π ( ) ∈ S 2 2 associated with such maps defines the integer-valued topological charge which can be used to classify various field configurations. Heuristically,  counts the number of times that the target sphere S 2 is wrapped while covering the xy-plane. Singular configurations for which the pseudo-spin is not everywhere well-defined, have  = 0. Non-singular solutions on the other hand, and in particular coreless vortices, have   π = Φ/ ∈ g 2 (where Φ is the flux). For example the two-quanta vortices, which are coreless, are characterized by  = 2. The fact that  = 0 for singular vortices can easily be seen from the plot of the pseudo-spin texture Fig. 4. There n never reaches the north pole and thus do not fully cover the unit sphere.

Discussion
Here, we reported a large-scale numerical investigation of the energy properties of isolated single and multiple quanta vortices/skyrmions in a Ginzburg-Landau model of chiral p-wave superconducting state. As pointed out previously, for a given ground-state chirality, vortices and anti-vortices are inequivalent. Thus we performed study for both orientations of the winding. The vortices with winding n − = − 1 in the dominant component are always preferred to those with winding n − = + 1. We also found that vortices carrying two flux quanta with n − = − 2 are always energetically favoured as compared to two isolated single-quanta vortices. Vortices with higher flux and negative n − , on the other hand, are either unstable or have higher energies per flux quantum. We also reported the structural and energetic properties of (meta-)stable skyrmions with various topological charge (i.e. for n − > + 1). The calculations show complicated sublinear scaling of the energy with the number of flux quanta that qualitatively agrees with previous works for a smaller parameter set in a related model 40 . Due to their very characteristic profiles of the magnetic field, their experimental observation, in e.g. scanning Hall and scanning SQUID experiments would provide a strong evidence of chiral p-wave superconductivity in the candidate materials described by the model (1). Note however that various aspects of microscopic physics may alter the form of the Ginzburg-Landau model (1), and in particular the balance between the different coefficients entering the free energy. This is currently a subject of ongoing studies, in connection with Sr 2 RuO 4 (see e.g. 19,20 ). Note added: after the completion of this work a study appeared reporting stable skyrmions as well as vortices, in this model, affected by mesoscopic effects in small samples 48 .

Methods
In this work we used the dimensionless two-component Ginzburg-Landau theory Eq. (1) that was previously microscopically derived in the weak coupling limit (see for example refs 34,35). In this work, we focus on the properties of vortex solutions in the xy-plane and neglect the dependence over the third coordinate z. This means that our solutions are either purely two dimensional, or describe bulk configuration, assuming translation invariance along z-axis (and thus neglecting possible surface effects). For the numerical investigation, the two-dimensional problem (1) is defined on a bounded domain Ω ⊂ R 2 . The boundary conditions for chiral p-wave superconductors can be very involved. Namely, in order to simulate chiral p-wave superconductors on a finite domain, a special attention has to be paid to boundary conditions to take into account edge currents properties. However, we are interested here in the intrinsic energetic properties of isolated defects. Thus we consider isolated vortices in large grids (such that there are no interactions with boundaries) and let the fields freely recover the ground-state. As a result, we probe the intrinsic structure and energy properties of vortices without any deformation originating from boundary behaviour. The simulation is run for a zero applied field (so that there are no Meissner currents), and the flux carrying solution is generated by a starting condition with a given winding of the dominant component. Because it enjoys topological protection, the (dominant) component cannot unwind by means of continuous transformations and thus topological properties (winding of the dominant component) are preserved by an energy minimization algorithm. Note that as simulations are run on a large but finite domain, there is still a possibility to change the topological sector, by moving the vortex across the boundary. This is possible, because without external fields there are no Meissner currents to prevent escape of a vortex. Note however that as we choose to work with large grids, the vortices in practice do not interact with boundaries, and thus they do not escape from the domain. The advantage of this choice is that it is guaranteed that obtained solutions are not affected by boundaries and that the calculated energies are those of isolated defects. The configurations displayed in the paper are close-up views of these defects.
For the actual numerical computation, the variational problem of minimizing the free energy is defined using a finite element formulation provided by the Freefem+ + library 42 . Discretization within finite element formulation is done via a (homogeneous) triangulation over Ω, based on Delaunay-Voronoi algorithm. Functions are decomposed over a continuous piecewise quadratic basis on each triangle. The accuracy of such method is controlled through the number of triangles, (we typically used 3 ~ 6 × 10 4 ), the order of expansion of the basis on each triangle (2nd order polynomial basis on each triangle), and also the order of the quadrature formula for the integral on the triangles. A nonlinear conjugate gradient algorithm is used to solve the variational nonlinear problem (i.e. to find the minima of  ). The algorithm is iterated until relative variation of the norm of the gradient of the functional  with respect to all degrees of freedom is less than 10 −8 (we verified that for this value, the configuration does not evolve and the energy remains constant).
For the minimization procedure to lead to a configuration that has the expected topological properties, the starting field configuration should exhibit itself those desired topological properties. Although strictly speaking there is no infinite energy barrier between different topological sectors in finite domains, the barrier is finite but large enough to prevent any unwinding. Thus typically gradient minimization converges to the configuration that has the topological properties of the starting guess. In order to have efficient numerics, it is also important that the starting field configuration is the closest as possible to the minimal energy configuration. The initial field configuration carrying N v flux quanta is prepared by using an ansatz which imposes phase winding of the dominant component (η − ) around a given point (x k , y k ): where  ( , ) = ( − ) + ( − ) x y x x y y k k k 2 2 and ξ a parametrizes the core size. The parametrization of η + , with nonzero density in the core enhances the convergence to form coreless defects. Finally, the starting configuration for the vector potential of the magnetic field A, is determined by solving Ampère's law equation ∇ × B + J = 0, for the supercurrent  δ δ = / J A specified by the superconducting condensates given by (4). Being an equation linear in A, this operation is rapidly solved. Once the starting configuration is constructed, all degrees of freedom are relaxed simultaneously.