Nematic Skyrmions in Odd-Parity Superconductors

TL;DR

This paper investigates topological excitations called nematic Skyrmions in two-component nematic superconductors, especially in Cu$_x$Bi$_2$Se$_3$, revealing their structure, interactions, and potential for forming bound states.

Contribution

It introduces the concept of nematic Skyrmions as coreless vortices with topological charge in nematic superconductors, highlighting their interactions and bound state formation.

Findings

Identifies coreless vortices as nematic Skyrmions with topological charge.

Shows inter-Skyrmion forces are dipolar and orientation-dependent.

Demonstrates formation of multi-Skyrmion bound states.

Abstract

We study topological excitations in two-component nematic superconductors, with a particular focus on Cu$_x$Bi$_2$Se$_3$ as a candidate material. We find that the lowest-energy topological excitations are coreless vortices: a bound state of two spatially separated half-quantum vortices. These objects are nematic Skyrmions, since they are characterized by an additional topological charge. The inter-Skyrmion forces are dipolar in this model, i.e. attractive for certain relative orientations of the Skyrmions, hence forming multi-Skyrmion bound states.

Equations

$$H(\mathbf{r})$$

$$\displaystyle\mathcal{H}_{\mathrm{BCS}}=-\sum_{\sigma,\sigma^{\prime}}\int d^{3}r\bigg{[}\Psi^{{\dagger}}_{\sigma 1}(\mathbf{r})\Psi^{{\dagger}}_{\sigma^{\prime}2}(\mathbf{r})\Delta_{\sigma^{\prime}\sigma}(\mathbf{r})+\mathrm{h.c.}\bigg{]}\,,\leavevmode\nobreak\ \leavevmode\nobreak\$$

$$H_{\textrm{BdG}}(\mathbf{r})=\bigg{[}\begin{matrix}H(\mathbf{r})-\mu&\Delta(\mathbf{r})\\ \Delta^{+}(\mathbf{r})&-\sigma_{y}H^{*}(\mathbf{r})\sigma_{y}+\mu\end{matrix}\bigg{]}$$

$$\displaystyle\mathcal{F}=\sum_{s=\pm}\bigg{\{}-|\Delta_{s}|^{2}+|D_{x}\Delta_{s}|^{2}+|D_{y}\Delta_{s}|^{2}+\beta_{z}|D_{z}\Delta_{s}|^{2}$$

$$\displaystyle+\beta_{\perp}(D_{-s}\Delta_{s})^{*}D_{s}\Delta_{-s}+\frac{|\Delta_{s}|^{4}}{2}+\frac{\gamma}{2}|\Delta_{s}|^{2}|\Delta_{-s}|^{2}\bigg{\}}\,,$$

$$(D_x^2+D_y^2+{\beta }_z{D}_z^2)$$

$$=-{\beta }_{\perp }{D}_s^2{\Delta }_{-s},$$

$$\mathbf{N}\cdot\mathbf{D}\Delta_{s}+(\beta_{z}-1)N_{z}D_{z}\Delta_{s}+N_{s}D_{s}\Delta_{-s}\leavevmode\nobreak\ \bigg{|}_{\mathrm{surf}}=0\,,$$

$$\displaystyle-2{\bm{\nabla}}\times\mathbf{b}=\sum_{s=\pm}\bigg{\{}[\Delta_{s}^{*}\mathbf{D}\Delta_{s}+(\beta_{z}-1)\Delta_{s}^{*}\hat{z}D_{z}\Delta_{s}+\mathrm{c.c.}]$$

$$\displaystyle+(\hat{x}+is\hat{y})\beta_{\perp}[\Delta_{s}^{*}D_{s}\Delta_{-s}+\Delta_{-s}(D_{-s}\Delta_{s})^{*}]\bigg{\}}.$$

$$\mathcal{Q}(\mathbf{n})=\frac{1}{4\pi }{\int }_{{\mathrm{R}}^2}\mathbf{n}\cdot {\partial }_x\mathbf{n}\times {\partial }_y\mathbf{n}dxdy.$$

$$\displaystyle H(\mathbf{r})=v\tau_{z}\left[{\bm{\sigma}}\times\bigg{(}-i{\bm{\nabla}}-\frac{e}{c}\mathbf{A}(\mathbf{r})\bigg{)}\right]\cdot\hat{z}+v_{z}\tau_{y}\bigg{(}-i\partial_{z}-\frac{e}{c}\mathbf{A}(\mathbf{r})\bigg{)}+m\tau_{x}+(\mathbf{h}+\mathbf{m}\tau_{x})\cdot{\bm{\sigma}},$$

$$\Psi (\mathbf{r})=({\Psi }_{\uparrow ,1}(\mathbf{r}),{\Psi }_{\downarrow ,1}(\mathbf{r}),{\Psi }_{\uparrow ,2}(\mathbf{r}),{\Psi }_{\downarrow ,2}(\mathbf{r}){)}^{\mathrm{T}}.$$

$$E_{\pm }(\mathbf{p})=\pm \{v^2{p}_{\perp }^2+v_z^2{p}_z^2+m^2{\}}^{1/2}.$$

$${\mathcal{H}}_{\mathrm{int}}=-\sum _{\tau }\int d^{3}r{\Psi }_{\downarrow \tau }^{†}(\mathbf{r}){\Psi }_{\uparrow \tau }^{†}(\mathbf{r})U_{\tau }{\Psi }_{\uparrow \tau }(\mathbf{r}){\Psi }_{\downarrow \tau }(\mathbf{r})-\sum _{\sigma ,{\sigma }^{\prime }}\int d^{3}r{\Psi }_{\sigma 1}^{†}(\mathbf{r}){\Psi }_{{\sigma }^{\prime }2}^{†}(\mathbf{r})\lambda {\Psi }_{{\sigma }^{\prime }2}(\mathbf{r}){\Psi }_{\sigma 1}(\mathbf{r})$$

$${\mathcal{H}}_{\mathrm{BCS}}$$

$$-$$

$$\mathcal{H}_{\textrm{BCS}}=\frac{1}{2}\int d^{3}r\Phi^{{\dagger}}(\mathbf{r})H_{\textrm{BCS}}(\mathbf{r})\Phi(\mathbf{r})+\int d^{3}r\bigg{[}\sum_{\tau}\frac{|\Delta_{\tau}(\mathbf{r})|^{2}}{U_{\tau}}+\sum_{\sigma,\sigma^{\prime}}\frac{|\Delta_{\sigma\sigma^{\prime}}(\mathbf{r})|^{2}}{\lambda}\bigg{]}$$

$$H_{\textrm{BCS}}(\mathbf{r})=\bigg{[}\begin{matrix}H(\mathbf{r})&\Delta(\mathbf{r})\\ \Delta^{{\dagger}}(\mathbf{r})&-\sigma_{y}H^{*}(\mathbf{r})\sigma_{y}\end{matrix}\bigg{]},$$

$${\sigma }_y{H}^{\ast }(\mathbf{r}){\sigma }_y=H(\mathbf{r}){|}_{\mathbf{A}(\mathbf{r})\to -\mathbf{A}(\mathbf{r})},$$

$$\Delta(\mathbf{r})=\frac{\Delta_{1}(\mathbf{r})+\Delta_{2}(\mathbf{r})}{2}+\frac{\Delta_{1}(\mathbf{r})-\Delta_{2}(\mathbf{r})}{2}\tau_{z}+\bigg{[}\begin{matrix}\Delta_{\uparrow\downarrow}(\mathbf{r})&\Delta_{\downarrow\downarrow}(\mathbf{r})\\ -\Delta_{\uparrow\uparrow}(\mathbf{r})&-\Delta_{\downarrow\uparrow}(\mathbf{r})\end{matrix}\bigg{]}\frac{\tau_{x}+i\tau_{y}}{2}-\bigg{[}\begin{matrix}\Delta_{\downarrow\uparrow}(\mathbf{r})&\Delta_{\downarrow\downarrow}(\mathbf{r})\\ -\Delta_{\uparrow\uparrow}(\mathbf{r})&-\Delta_{\uparrow\downarrow}(\mathbf{r})\end{matrix}\bigg{]}\frac{\tau_{x}-i\tau_{y}}{2}.$$

$$\Delta (\mathbf{r})=\frac{{\Delta }_1(\mathbf{r})+{\Delta }_2(\mathbf{r})}{2}+\frac{{\Delta }_1(\mathbf{r})-{\Delta }_2(\mathbf{r})}{2}{\tau }_z+\mathbf{\sigma }\cdot \mathbf{\Delta }(\mathbf{r}){\tau }_y,$$

$${\Delta }_x(\mathbf{r})=-\frac{i}{2}({\Delta }_{\uparrow \uparrow }(\mathbf{r})-{\Delta }_{\downarrow \downarrow }(\mathbf{r})),\text{leavevmode}\text{leavevmode}{\Delta }_y(\mathbf{r})=-\frac{1}{2}({\Delta }_{\uparrow \uparrow }(\mathbf{r})+{\Delta }_{\downarrow \downarrow }(\mathbf{r}))\text{leavevmode}\text{leavevmode}{\Delta }_z(\mathbf{r})=\frac{1}{2}({\Delta }_{\uparrow \downarrow }(\mathbf{r})+{\Delta }_{\downarrow \uparrow }(\mathbf{r}))$$

$$H_{\textrm{BCS}}(\mathbf{k})=\bigg{[}\begin{matrix}H(\mathbf{k})&\Delta\\ \Delta^{{\dagger}}&-H(\mathbf{k})\end{matrix}\bigg{]},$$

$$H(\mathbf{k})$$

$$\Delta$$

$${\mathcal{E}}^2(\mathbf{p})=(E_{+}(\mathbf{p})-\mu {)}^2+|\Delta {|}^2\left[1-\frac{(1-\cos \varphi )(m^2+v_z^2{p}_z^2)}{2{\mu }^2}\right],$$

$${\mathcal{E}}^2(\mathbf{p})=(E_{+}(\mathbf{p})-\mu {)}^2+|{\Delta }_z{|}^2(1-m^2/{\mu }^2).$$

$$\displaystyle\mathcal{E}^{2}(\mathbf{p})=(E_{+}(\mathbf{p})-\mu)^{2}+(\Delta_{x}^{2}+\Delta_{y}^{2})\bigg{(}1-\frac{m^{2}}{\mu^{2}}\bigg{)}-\frac{(vp_{x}\Delta_{y}-vp_{y}\Delta_{x})^{2}}{\mu^{2}}$$

$${\sin }^2(\varphi -\theta )=\frac{{\mu }^2-m^2}{v^2{p}_{\perp ,F}^2},$$

$$\mathbf{S}(\mathbf{r})=\underset{{\mathbf{r}}^{\prime }\to \mathbf{r}}{\lim }\text{Tr}\text{leavevmode}T\sum _n\frac{\mathbf{\sigma }}{2}G(\mathbf{r},{\mathbf{r}}^{\prime },{\omega }_n),$$

$${\mathbf{P}}_i(\mathbf{r})=\underset{{\mathbf{r}}^{\prime }\to \mathbf{r}}{\lim }\text{Tr}\text{leavevmode}T\sum _n\frac{\mathbf{\sigma }}{2}{\tau }_{i}G(\mathbf{r},{\mathbf{r}}^{\prime },{\omega }_n),$$

Full text

Bulk superconductivity in topological insulator materials has recently been observed in electron-doped $\mathrm{Cu_{x}Bi_{2}Se_{3}}$ in Refs. Hor et al. (2010); Wray et al. (2010); Kriener et al. (2011) and unusual superconducting states in this system were theoretically considered in Ref. Fu and Berg (2010). There, it was argued that the fully-gapped single order parameter superconductor, which has a spin-triplet pairing with odd parity and that possesses topologically protected gapless surface states, is favored over other superconducting states. It was later put forward Fu (2014) that the interplay of the crystal lattice anisotropy and the nematic superconductivity might be consistent with Knight-shift anisotropy measurements in $\mathrm{Cu_{x}Bi_{2}Se_{3}}$, that show spontaneous breaking of the spin-rotation symmetry below the superconducting transition temperature Matano et al. (2016). In this model, the nematic superconducting state has an odd-parity spin-triplet pairing and is described by a two-component order parameter that spontaneously breaks the rotational symmetry in the basal plane of the lattice. This scenario of nematic superconductivity is supported by the recent observation of two fold rotational symmetry of the magnetic field in specific heat and upper critical field measurements of the superconducting state Yonezawa et al. (2017). Bulk superconductivity was also reported in $\mathrm{Sr_{x}Bi_{2}Se_{3}}$ Shruti et al. (2015); Liu et al. (2015); Pan et al. (2016) and in magnetically doped $\mathrm{Nb_{x}Bi_{2}Se_{3}}$ Qiu et al. ; Asaba et al. (2017), where upper critical field Pan et al. (2016) and magnetic torque measurements Asaba et al. (2017) reveal signatures of rotational symmetry breaking in the amplitude of the superconducting gap.

Nematic superconducting states have interesting properties that were recently theoretically addressed. For example, such states show a specific anisotropy of the upper critical magnetic field Venderbos et al. (2016a), undergo phase transitions to superconducting states that break the time-reversal symmetry as a result of the interplay of ferromagnetism and superconductivity Yuan et al. (2017); Chirolli et al. (2017), and host Majorana fermions at the surface Sasaki et al. (2011); Venderbos et al. (2016b). This raises the question of the properties of topological excitations in this kind of materials. Recent work Wu and Martin (2017) presented an ansatz-based investigation of Kramers pairs of Majorana fermions bound inside a specific type of composite vortices that do not carry magnetic flux. In the model that we consider below, such an ansatz describes unstable vortex solutions.

In this Letter we address the question of the nature of the lowest-energy vortex excitations in nematic superconductors. To this end, we investigate vortex solutions in a two-component Ginzburg-Landau (GL) model consistently derived from the microscopic theory. We find that the lowest-energy topological excitations are coreless and consist of two spatially separated half-quantum vortices (HQVs) 111Half-quantum vortices are of substantial current interest in various contexts due to proposals that they host Majorana fermions, see for example Refs. Read and Green (2000); Ivanov (2001). , such that the total superconducting density has no zeros. These excitations can be characterized by a Skyrmion topological index, which is zero for singular vortices. Heuristically, the Skyrmion terminology follows from the fact that the coreless vortex can be seen as a texture of a unit vector that fully covers the target two-sphere. The unit vector that maps to the target two-sphere is defined as a projection of the superconducting degrees of freedom onto the vector composed of the Pauli matrix set.

Such a coreless vortex, which is a bound state of two HQVs, shows as a dipolelike configuration of the relative phase between the components of the order parameter, and thus can mediate a long-range dipole interaction between the Skyrmions, thus binding them together into a multi-Skyrmion bound state.

We consider a model of a three-dimensional topological insulator in the presence of a magnetic field, having in mind $\mathrm{Bi_{2}Se_{3}}$ as a particular material candidate, which is a narrow gap semiconductor with a layered crystal structure. The system is described by the Hamiltonian $\mathcal{H}=\int\Psi^{{\dagger}}(\mathbf{r})H(\mathbf{r})\Psi(\mathbf{r})d^{3}r$, with

[TABLE]

where $\mathbf{A}(\mathbf{r})$ is the vector potential, $e<0$ is the electric charge, and $m$ describes the coupling between the orbitals of $\mathrm{Bi_{2}Se_{3}}$. Here, $v$ and $v_{z}$ are the Fermi velocities that characterize the anisotropic dispersion of the massive Dirac fermion in the absence of the magnetic field: $E_{\pm}(\mathbf{p})=\pm\{v^{2}(p_{x}^{2}+p_{y}^{2})+v_{z}^{2}p_{z}^{2}+m^{2}\}^{1/2},$ where $\mathbf{p}=(p_{x},p_{y},p_{z})$ is the momentum of a particle. The Pauli matrices $\sigma_{a}$ and $\tau_{a}$ (with $a=x,y,z$), respectively, describe the real spin ($\uparrow,\downarrow$) and the orbital pseudospin ($1,2$) degrees of freedom. The electron operator is given by $\Psi(\mathbf{r})=(\Psi_{\uparrow,1}(\mathbf{r}),\Psi_{\downarrow,1}(\mathbf{r}),\Psi_{\uparrow,2}(\mathbf{r}),\Psi_{\downarrow,2}(\mathbf{r}))^{\mathrm{T}}$, and $\hbar=1$ units are used here, and spin and pseudospin indices are omitted for clarity of notation throughout the Letter. The Zeeman contribution of the magnetic field to the Hamiltonian (1) is neglected compared to that the orbital effect. We also note that, although there is strong spin-orbit interaction in each orbital, the inversion symmetry of the system is preserved.

As demonstrated in Ref. Fu and Berg (2010) the electron-phonon interaction might lead to several distinct $s$-wave superconducting instabilities in this system: intraorbital spin-singlet, interorbital spin-singlet, and interorbital spin-triplet. Motivated by the experimental signatures for the nematic superconductivity, we focus here on the interorbital spin-triplet pairing, which is described by the interaction Hamiltonian within the Bardeen–Cooper–Schrieffer (BCS) approximation:

[TABLE]

where $\Delta_{\sigma^{\prime}\sigma}(\mathbf{r})=\lambda\langle\Psi_{\sigma^{\prime}2}(\mathbf{r})\Psi_{\sigma 1}(\mathbf{r})\rangle$ with the interaction constant $\lambda>0$. Note that we consider the zero harmonic of the electron-phonon interaction potential, which shall give a higher temperature of the superconductor-metal phase transition than higher harmonics. To proceed we introduce the Bogolyubov-de Gennes (BdG) Hamiltonian $\mathcal{H}_{\textrm{BdG}}=\frac{1}{2}\int\Phi^{{\dagger}}(\mathbf{r})H_{\textrm{BdG}}(\mathbf{r})\Phi(\mathbf{r})d^{3}r$, where

[TABLE]

is written in the Nambu notations: $\Phi^{{\dagger}}(\mathbf{r})=(\Psi^{{\dagger}}(\mathbf{r}),\Psi^{\textrm{T}}(\mathbf{r})(-i\sigma_{y}))$. We assume the Fermi level to be in the conduction band and hence set $\mu>|m|$. In what follows we consider the interorbital spin-triplet pairing of the form $\Delta(\mathbf{r})={\bm{\sigma}}\cdot\bm{\Delta}(\mathbf{r})\tau_{y}$, where $\bm{\Delta}(\mathbf{r})=[\Delta_{x}(\mathbf{r}),\Delta_{y}(\mathbf{r}),0]$, such that $\Delta_{x}(\mathbf{r})=-i[\Delta_{\uparrow\uparrow}(\mathbf{r})-\Delta_{\downarrow\downarrow}(\mathbf{r})]/2$, and $\Delta_{y}(\mathbf{r})=-[\Delta_{\uparrow\uparrow}(\mathbf{r})+\Delta_{\downarrow\downarrow}(\mathbf{r})]/2$. The strong spin-orbit interaction locks electron spin to the momentum, thus fixing the orientation of vector $\Delta$ to $[\Delta_{x}(\mathbf{r}),\Delta_{y}(\mathbf{r}),0]$.

In order to investigate the structure of the topological excitations, we derive microscopically the GL free energy functional for the two-component order parameter $\Delta_{\pm}=(\Delta_{x}\pm i\Delta_{y})/\sqrt{2}$ (for details of the derivation, see Supplemental Material [ThedetailsofthederivationaregivenintheSupplementalMaterial(SM)][]Supplementary ). The scaled Ginzburg-Landau free energy functional reads as $4\pi F=h^{2}_{0}l^{3}_{m}\int\left(\mathcal{F}(\mathbf{R})+[{\bm{\nabla}}\times\mathbf{a}(\mathbf{R})]^{2}\right)d^{3}R$, where

[TABLE]

where $\Delta_{\pm}=|\Delta_{\pm}|\text{e}^{i\varphi_{\pm}}$. Here we have defined the order parameter $\Delta_{s}(\mathbf{R})$ in the scaled form. The explicit scaling transformation for coordinate $\mathbf{R}$, vector potential $\mathbf{a}(\mathbf{R})$, and operator $D_{\pm}=D_{x}\pm iD_{y}$ in which $\mathbf{D}=-(i/\kappa)\bm{\nabla}_{R}+\mathbf{a}(\mathbf{R})$ is given in Supplemental Material Sup . In the dimensionless units, the coupling constant $\kappa$ parametrizes the magnetic field penetration length relative to a characteristic length scale associated with density variation. Its explicit expression through the parameters of the microscopic model is given in the Supplemental Material Sup [(F.73)], as well as the expression for other coefficients $h_{0}$, $l_{m}$, $\beta_{\perp,z}$, and $\gamma$. For $\kappa\ll 1$ the vortex core energy is relatively large and the thermodynamically stable Skyrmions do not form.

The anisotropy of the electronic spectrum results in the anisotropic gradient terms on the first line of the GL functional (4). The interplay of the two-component order parameter and strong spin-orbit interaction gives rise to the mixed gradient term with the coefficient $0<\beta_{\perp}<1$ on the second line in Eq. (4). The GL free energy functional density is invariant under the joint rotation of coordinates and the components of the order parameter 222For example, for a rotation of coordinates $(x,y)\rightarrow(y,-x)$ the components of the order parameter rotate $(\Delta_{x},\Delta_{y})\rightarrow-(\Delta_{y},-\Delta_{x})$.. Finally, the value of $\gamma$ determines whether the superconductor is in the nematic $0<\gamma<1$ or in the chiral $\gamma>1$ phase. Indeed, in the spatially homogeneous case, when $0<\gamma<1$ the free energy is minimal if the order parameter has the form $\Delta=\Delta_{0}(\cos\theta,\sin\theta,0)$ (for some real constant $\theta$ and $|\Delta_{0}|=1/\sqrt{1+\gamma}$), while it reads as $\Delta=\tilde{\Delta}_{0}(1,\pm i,0)$ (where $|\tilde{\Delta}_{0}|=1$) when $\gamma>1$.

Time-reversal symmetry is preserved in the nematic state. The chiral state, on the other hand is characterized by the nonzero electron spin polarization $\propto|\bm{\Delta}(\mathbf{r})\times\bm{\Delta}^{*}(\mathbf{r})|\neq 0$. Interestingly, nontrivial pseudospin polarization, antiferromagnetic spin orientation in two orbitals, shows up in the first gradients of the order parameter in both chiral and nematic cases, see Supplemental Material Sup .

The GL equation for the component $\Delta_{s}$ of the order parameter is obtained by the functional variation of Eq. (4) with respect to $\Delta_{s}^{*}$:

[TABLE]

and is supplemented by the boundary condition

[TABLE]

where $N_{s}=N_{x}+isN_{y}$ and $\mathbf{N}$ is the unit vector directed normal to the surface. Note that here we do not consider the effects of the localized surface states. Finally, the vector potential $\mathbf{a}(\mathbf{R})$ of the magnetic field $\mathbf{b}(\mathbf{R})={\bm{\nabla}}_{R}\times\mathbf{a}(\mathbf{R})$ satisfies

[TABLE]

We now turn to the investigation of the nature of topological excitations in the nematic superconductor. In two-component models, due to the coupling of the components to the vector potential $\mathbf{a}(\mathbf{R})$, the only solutions with finite energy per unit length have the same phase winding in both components of the order parameter, that is, a bound state of vortices in the different components, each carrying a fraction of magnetic flux that adds up to a single flux quantum. In the current model of the nematic superconductor vortices in each component of the order parameter carry half of a magnetic flux quantum; hence, they are half-quantum vortices (HQVs) 333This follows from that both components have the same ground state density. In general this ’symmetry’ can be broken, leading to a more generic fractional quantization of the flux, see detailed discussion in Babaev (2002).. Typically, the magnetic interaction between HQVs favors cocentricity of the vortex cores in the different components (see, e.g., a detailed discussion in Ref. Babaev (2002)). On the other hand, the model considered here also features mixed gradients and biquadratic density-density terms that result in the repulsion between the cores of the half-quantum vortices. Provided the latter dominate, the competition between those forces may result in a bound state of nonoverlapping half-quantum vortices, thus breaking the axial symmetry of the solution.

To address whether vortices are singular (cocentered HQV) or coreless (i.e., noncocentered HQV), we numerically construct vortex solutions by minimizing the free energy (4), starting by an initial configuration, in which both components $\Delta_{\pm}$ have the same winding. The theory is discretized within a finite-element formulation Hecht (2012), and minimized using a nonlinear conjugate gradient algorithm [Fordetaileddiscussiononthenumericalmethods; seeforexamplerelateddiscussionin:][]Garaud.Babaev.ea:16. Minimization procedure leads, after the convergence of the algorithm, to a vortex configuration that carries a number of flux quanta that is specified by the initial phase winding. Figure 1 shows such a single-quantum vortex configuration in the model (4) for the nematic superconductor. Note that the picture shows a close-up view, displaying only a small part of the simulated numerical grid, which is chosen to be large enough so that vortices do not interact with the boundaries. Clearly, the vortex solution is not axially symmetric. Inspection of the core structure reveals that HQVs in different components are spatially separated and thus that this bound state of HQVs is coreless; i.e., there is no singularity of total density of superconducting components: $[|\Delta_{+}(\mathbf{R})|^{2}+|\Delta_{-}(\mathbf{R})|^{2}]^{1/2}$. We simulated vortex solutions for various initial guesses that always converge to configurations as in Fig. 1. All investigated values of the parameter $\kappa$ led to coreless vortices in the type-II regime. The distance between HQVs is determined by the competition between magnetic attraction and repulsion mediated by other terms such as density-density interaction. This cannot be addressed analytically, but quantitatively it can be seen that increasing the value of $\kappa$ decreases binding the HQVs thus increasing their separation. Since the superconductor is substantially away from the type-I regime, and because these excitations are energetically cheaper than singular vortices, a lattice of Skyrmions will form in external field.

As singularities in both components do not overlap, there is a dipolelike configuration of the relative phase $\varphi_{-}-\varphi_{+}$ between the components. Importantly, the phase-difference gradients are very strong. This indicates that intervortex forces include torque and a long-range dipole interaction, that can lead to a long-range attraction between single-quanta vortices. As displayed in Fig. 2, by initially setting a double phase winding in each component, we find that indeed two single-quantum vortices form a bound state due to the dipolar forces. Note that the presence of dipolar interactions usually has clear signatures in structure formation. In particular, studies of Skyrmion solutions in other systems with dipolar inter-Skyrmion forces Jäykkä and Speight (2010); Garaud et al. (2014); Agterberg et al. (2014) show, for example, that hexagonal symmetry is unfavorable for the Skyrmion lattices. Moreover, the long-range dipolar interaction can also result in long-range attractive interaction between Skyrmions and boundaries of a superconductor, thus suggesting possible abundance of topological defects near boundaries in a weak applied magnetic field.

Bound states of nonoverlapping HQVs are coreless defects that can be called Skyrmions, and the reason for that terminology is that they exhibit additional topological properties, as compared to singular vortices. These can be seen by introducing the unit vector $\bf n$ defined as the projection of the superconducting degrees of freedom $\eta^{\dagger}=(\Delta_{+}^{*},\Delta_{-}^{*})$ onto spin-$1/2$ Pauli matrices $\tilde{\bm{\sigma}}$, as ${\bf n}=\eta^{\dagger}\tilde{\bm{\sigma}}\eta/\eta^{\dagger}\eta$. That is, the $x$ and $y$ components of the vector ${\bf n}$ depend on the phase difference, while the $z$ component is determined by the ratio of the moduli of the complex fields. The associated projection is a map from the one-point compactification of the plane ($\mathrm{R}^{2}\cup\{\infty\}\simeq S^{2}$) onto the two-sphere target space spanned by $\bf n$. That is ${\bf n}:S^{2}\to S^{2}$, which is classified by the homotopy class $\pi_{2}(S^{2})\in\mathrm{Z}$. This defines the integer-valued $\mathrm{C}\mathbb{P}^{1}$ topological invariant, as

[TABLE]

If $\eta\neq 0$ everywhere (coreless vortex), ${\cal Q}$ is an integer number. In a way, ${\cal Q}$ counts the number of times the texture of $\bf n$ covers the target two-sphere.

Figure 3 shows the texture of the unit vector $\bf n$ that corresponds to the vortices in the nematic phase. The left panel corresponds to the single-quantum vortex displayed in Fig. 1. It illustrates that the unit vector ${\bf n}$ wraps the target two-sphere (once), thus implying this configuration has unit Skyrmionic charge ${\cal Q}=1$. The right panel shows the texture corresponding to the bound state of Fig. 2 that originates in long-range dipolelike forces. Note that this illustrates the dipole nature of the long-range interaction. Indeed the pair of Skyrmions alternates the north (red) and south (blue) poles of the target sphere.

It is worth noting that various coreless vortices were also considered to exist in a number of models of multicomponent superconductivity. There are, however, substantial differences in the structure and properties of these solutions, and they should have distinct experimental manifestations. For example, in the framework of various models of $p$-wave superconductors, it was advocated that multiquanta coreless vortices may be favored over single-quanta singular vortices Sauls and Eschrig (2009); Ichioka et al. (2012); Knigavko and Rosenstein (1999). The recent numerical studies show that two-quanta Skyrmions in the GL models of a chiral $p$-wave superconductor are energetically favored, and hence single-quanta Skyrmions do not form in the ground state in an external field Garaud and Babaev (2015); Garaud et al. (2016); Fernández Becerra et al. (2016); Zhang et al. (2016); Fernández Becerra and Milošević (2017). Different types of chiral Skyrmions were also discussed for $s+is$ superconducting states Garaud et al. (2011, 2013); Garaud and Babaev (2015). The structure of chiral Skyrmions is significantly dissimilar compared to the nematic Skyrmions since the former have fractional vortices and magnetic flux pinned on domain walls between different time-reversal symmetry broken states Garaud and Babaev (2014); thus, they have very different magnetic field configurations.

The object that we find here should be interesting from the viewpoint of electronic states. Indeed, HQVs are known to possess Majorana modes. The spatial separation between two HQVs in the Skyrmion implies that individual HQVs may be rather easily stabilized in a mesoscopic sample.

In conclusion, we discussed the topological excitations in nematic superconductors. We showed that the topological excitations are nematic Skyrmions, each of which can be viewed as a bound state of two spatially separated half-quantum vortices. The nematic Skyrmions have orientation-dependent dipolar attractive forces and form multiquanta bound states, which could be expected to have clear experimental signatures in structure formation. Moreover, being coreless, the Skyrmions are expected to have unusual electronic core-state properties that could allow the identification of Skyrmionic states using scanning tunneling microscopy.

Supplemental Material to

”Nematic Skyrmions in odd-parity superconductors”

In the Supplemental Material, we discuss the spectrum of quasiparticles and derive the Ginzburg-Landau free energy functional of the nematic superconductor introduced in the main text. We also evaluate expression for the spin polarization in the nematic superconductor.

.1 Main definitions

We introduce the model of the superconductivity in the doped Dirac material. As a particular example we consider theoretical model of a doped $\mathrm{Bi_{2}Se_{3}}$. This is a spatial inversion and time reversal symmetric layered material with two orbitals within single quintuple layer. For concreteness, we assume that $z$-axis is normal to the plane of the layers. The low energy excitations in the presence of the magnetic field can be described by the Hamiltonian $\mathcal{H}=\int d^{3}r\Psi^{{\dagger}}(\mathbf{r})H(\mathbf{r})\Psi(\mathbf{r})$, where

[TABLE]

in which $v$ and $v_{z}$ are the in-plane and out of plane components of the Fermi velocity, $|m|$ is the mass of the Dirac fermion describing electron tunneling between the orbitals of the quintuple layer, $e<0$ is the charge of electron, ${\bm{\sigma}}$ and ${\bm{\tau}}$ are the Pauli matrices describing the spin and orbital degrees of freedom, such that

[TABLE]

The presence of two orbitals with opposite sign of the coefficient, $v$, supports inversion symmetry in the system. The Zeeman effect of the external magnetic field results in the spin polarization of electrons in the orbital $\mathbf{h}\cdot{\bm{\sigma}}$ and spin polarized inter-orbital coupling $\mathbf{m}\cdot{\bm{\sigma}}\tau_{x}$. Although we have written these contributions in the Hamiltonian for completeness we will ignore Zeeman effect further compared to the orbital effect of the magnetic field. The spectrum of particles in the absence of magnetic field is given by

[TABLE]

This spectrum is anisotropic and has a gap $2m$ at $p\equiv|\mathbf{p}|=0$.

We consider phonon mediated attractive interaction, which can lead to the instabilities in the inter-orbital and intra-orbital couplings. Part of the Hamiltonian, which describes attractive electron-electron interaction is given by

[TABLE]

We note that we consider the zero harmonic of the electron-phonon interaction potential in the intra-orbital, $U_{\tau}>0$, and inter-orbital channels, $\lambda>0$, which shall give higher temperature of the superconductor-metal phase transition compared to higher harmonics. It can be shown that the spectrum of the Bogolubov quasiparticles in the case of the inter-orbital spin-singlet pairing contains nodal line. We ignore this coupling mechanism here, since such a pairing is energetically unfavorable compared to the intra-orbital spin-singlet and inter-orbital spin-triplet pairings, which are either fully gapped or contain point nodes in the spectrum.

Thus, taking into account only orbital effects of the magnetic field, the BCS Hamiltonian is given by

[TABLE]

Here second line describes intra-orbital spin-singlet pairing: $\Delta_{\tau}(\mathbf{r})=U_{\tau}\langle\Psi_{\uparrow\tau}(\mathbf{r})\Psi_{\downarrow\tau}(\mathbf{r})\rangle$, $\Delta_{\tau}^{+}(\mathbf{r})=U_{\tau}\langle\Psi^{{\dagger}}_{\downarrow\tau}(\mathbf{r})\Psi^{{\dagger}}_{\uparrow\tau}(\mathbf{r})\rangle$, third line describes inter-orbital pairing: $\Delta_{\sigma^{\prime}\sigma}(\mathbf{r})=\lambda\langle\Psi_{\sigma^{\prime}2}(\mathbf{r})\Psi_{\sigma 1}(\mathbf{r})\rangle$, $\Delta_{\sigma\sigma^{\prime}}^{+}(\mathbf{r})=\lambda\langle\Psi^{{\dagger}}_{\sigma 1}(\mathbf{r})\Psi^{{\dagger}}_{\sigma^{\prime}2}(\mathbf{r})\rangle$, which might contain both spin-singlet and spin-triplet components; $\mu$ is the Fermi energy. The BCS Hamiltonian can be rewritten as

[TABLE]

where we have performed unitary transformation and introduced new operators $\Phi^{{\dagger}}(\mathbf{r})=(\Psi^{{\dagger}}(\mathbf{r}),\Psi^{\textrm{T}}(\mathbf{r})(-i\sigma_{y}))$, which are written in the Nambu space. The Hamiltonian is given by

[TABLE]

We note that here each block is a $4\times 4$ matrix and that Hamiltonian satisfies an equality

[TABLE]

Order parameters can be generally written as

[TABLE]

Neglecting spin-singlet inter-orbital coupling one can parametrize superconducting order parameters as

[TABLE]

where we introduce a vector $\bm{\Delta}(\mathbf{r})=(\Delta_{x}(\mathbf{r}),\Delta_{y}(\mathbf{r}),\Delta_{z}(\mathbf{r}))$, such that

[TABLE]

and assume that the inter-orbital spin-singlet pairing $\propto\Delta_{\uparrow\downarrow}(\mathbf{r})-\Delta_{\downarrow\uparrow}(\mathbf{r})=0$ vanishes.

.2 Eigenvalues of the BCS Hamiltonian

Here we consider the spectrum of quasiparticles in the Dirac superconductor in few limiting cases. We set external magnetic field to zero and consider coordinate independent order parameters. The BCS Hamiltonian in the momentum representation is given by

[TABLE]

where

[TABLE]

We will assume that the Fermi level is in the conduction band $\mu>|m|$ and for concreteness consider several representative cases of the order parameter.

1. Consider the intra-orbital pairing and assume $\Delta_{1}=\Delta e^{i\phi},\Delta_{2}=\Delta$, where $\phi$ is the phase difference of the order parameters in two orbitals. The spectrum of quasiparticles is described by

[TABLE]

which is gapped if the phase difference is $\phi=0$, while it has two point nodes at $\mathbf{p}_{\pm}=(0,0,\pm\sqrt{\mu^{2}-m^{2}}/v_{z})$ if the phase difference $\phi=\pi$.

2. Consider inter-orbital spin-triplet pairing in few limiting cases. In the case $\mathbf{\Delta}=(0,0,\Delta_{z})$ the spectrum is gapped

[TABLE]

3. In the case when $\bm{\Delta}=(\Delta_{x},\Delta_{y},0)$ and when both components $\Delta_{x,y}$ are real,the low energy excitations at the Fermi level are described by

[TABLE]

There exist point nodes in the plane $p_{z}=0$. For example, if $\bm{\Delta}=\Delta(\cos(\theta),\sin(\theta),0)$, then the positions of the nodes are determined by the solution of equation

[TABLE]

where we set $p_{x}=p_{\perp,F}\cos(\phi),p_{y}=p_{\perp,F}\sin(\phi)$ for the components of the in-plane momentum at the Fermi level.

.3 Spin polarization: magnetoelectric effect in nematic superconductor

In this section we will derive expression for the quasiparticle spin polarization in the superconductor in the absence of the magnetic field. The electron spin density $\mathbf{S}(\mathbf{r})$ can be written as follows:

[TABLE]

where $G(\mathbf{r},\mathbf{r}^{\prime},\omega_{n})$ is the Green function in Matsubara representation, $T$ is the temperature, and trace is taken over spin and layer degrees of freedom. One can also define pseudospin

[TABLE]

which describes relative spin polarization in different orbitals. One has two coupled equations for normal and anomalous Green functions

[TABLE]

supplemented by the self-consistency equation for the order parameter

[TABLE]

where $\lambda$ is the interaction constant. In order to find expression for the Green function $G(\mathbf{r},\mathbf{r}^{\prime},\omega_{n})$, we expand in powers of $\Delta(\mathbf{r})$ up to the second order. Thus, solving equations for the Green functions perturbatively in the order parameter

[TABLE]

we obtain

[TABLE]

where Green functions $G_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega_{n})$ and $\bar{G}_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega_{n})$ correspond to the case of the absence of the superconductivity and satisfy equations

[TABLE]

We can rewrite expression for the spin polarization in the form

[TABLE]

We are interested in the terms linear to the gradients of the order parameter.

[TABLE]

where we are using momentum representation of the Green function

[TABLE]

which is given by

[TABLE]

where $E_{\pm}(\mathbf{p})=\pm\{v^{2}p_{\perp}^{2}+v_{z}^{2}p_{z}^{2}+m^{2}\}^{1/2}$ is the spectrum of particles in the non-superconducting state and

[TABLE]

is the projector operator on electron or hole bands, which are defined by $s=\pm$. Operator $D_{s}(\mathbf{p})$ satisfies an equality $D_{s}^{2}(\mathbf{p})/4=D_{s}(\mathbf{p})/2.$ We also observe that

[TABLE]

and

[TABLE]

Now we consider the case when the Fermi energy is in the conduction band, $\mu>|m|$, in which it is enough to take into account contribution from the band with $s=+$. Noting that

[TABLE]

we find expression for the spin polarization

[TABLE]

We observe that only second term on the right hand side of Eq. C.35 contributes to the spin polarization through the first gradients of the order parameter. Consider simple current carrying expression for the order parameters $\Delta_{1,2}(\mathbf{r})=|\Delta|e^{i\phi_{1,2}(\mathbf{r})}$. We finally obtain expression for the spin polarization in the Dirac superconductor

[TABLE]

To summarize this section,

1. Finite term $\propto|\bm{\Delta}(\mathbf{r})\times\bm{\Delta}^{*}(\mathbf{r})|\neq 0$ describes spin polarization in the spontaneously time reversal symmetry broken state of the superconductor. Nontrivial pseudospin polarization $\mathbf{P}_{z}(\mathbf{r})$ shows up in the first gradients of the order parameter in both chiral and nematic cases. Interestingly, pseudospin polarization results in the antiferromagnetic spin orientation in two orbitals even in the nematic case.

2. Finite gradient of the phase difference of the order parameters in two orbitals ${\bm{\nabla}}[\phi_{1}(\mathbf{r})-\phi_{2}(\mathbf{r})]\neq 0$ breaks time reversal symmetry and thus leads to the spin polarization in the superconductor, $\propto[\hat{z}\times{\bm{\nabla}}\{\phi_{1}(\mathbf{r})-\phi_{2}(\mathbf{r})\}]$. This is the well known magneto-electric effect in superconductors. Contrary, if ${\bm{\nabla}}[\phi_{1}(\mathbf{r})-\phi_{2}(\mathbf{r})]=0$, while ${\bm{\nabla}}\phi_{1,2}(\mathbf{r})\neq 0$, then there exists finite pseudo-spin polarization, which corresponds to the antiferromagnetic in-plane spin orientation. Since there might be in-plane easy axis, anisotropy of the pseudospin polarization is expected.

.4 Derivation of Ginzburg-Landau equation

In this section we will derive Ginzburg-Landau equation and Ginzburg-Landau functional density for the case of inter-orbital spin-triplet pairing

[TABLE]

Due to strong anisotropy of the spectrum, it is convenient to separate the case $\bm{\Delta}(\mathbf{r})=(\Delta_{x}(\mathbf{r}),\Delta_{y}(\mathbf{r}),0)$ from $\bm{\Delta}(\mathbf{r})=(0,0,\Delta_{z}(\mathbf{r}))$. In the latter case the Ginzburg-Landau free energy functional has standard expression for the single order parameter. Thus, in what follows we set $\Delta_{z}=0$ and consider a vector $\bm{\Delta}(\mathbf{r})=(\Delta_{x}(\mathbf{r}),\Delta_{y}(\mathbf{r}),0)$, which allows to write

[TABLE]

We start from the self-consistency equation for the order parameter

[TABLE]

where $\lambda$ is the interaction constant. Anomalous Green function is defined by

[TABLE]

Thus in momentum representation (we will be include orbital effect of the magnetic field in the final expressions by introducing the gauge invariant derivatives) we obtain

[TABLE]

where we have neglected derivatives of the order parameter in the non-linear in powers of the order parameter terms. Using that the Green functions satisfy equality

[TABLE]

we obtain

[TABLE]

Here we also include fifth order term. Since we consider that the Fermi level is in the conduction band and the Fermi energy $\mu>|m|$ is the largest energy scale in the system, we can approximate the Green function as

[TABLE]

In the intermediate derivations one obtains terms which contain gradients of the Green function over momentum $\propto\partial_{\mathbf{p}}G_{0}(\mathbf{p},\omega_{n})$ under the integral over momentum. The natural limit of small temperature of the superconductor-metal phase transition $T_{c}/\mu\ll 1$, allows to approximate these terms as $\frac{D_{+}(\mathbf{p})}{2}\partial_{\mathbf{p}}\frac{1}{i\omega_{n}+\mu-E_{+}(\mathbf{p})}$.

We derive Ginzburg-Landau equation for the spin-triplet inter-orbital order parameter $\bm{\Delta}(\mathbf{r})=(\Delta_{x}(\mathbf{r}),\Delta_{y}(\mathbf{r}),0)$ in the form

[TABLE]

where

[TABLE]

Coefficients are given by:

[TABLE]

Here $\xi=E_{+}(\mathbf{p})-\mu,\leavevmode\nobreak\ \leavevmode\nobreak\ E_{+}(\mathbf{p})=\sqrt{m^{2}+v_{z}^{2}p_{z}^{2}+v^{2}p_{\perp}^{2}}$, and $(dp)=\frac{d^{3}p}{(2\pi)^{3}}=\frac{p_{\perp}dp_{\perp}}{2\pi}\frac{dp_{z}}{2\pi}$. The Ginzburg-Landau free energy functional is given by

[TABLE]

where $\bm{\Delta}\cdot\bm{\Delta}^{*}=|\Delta_{x}|^{2}+|\Delta_{y}|^{2}$ and $\bm{\Delta}\times\bm{\Delta}^{*}=\Delta_{x}\Delta_{y}^{*}-\Delta_{x}^{*}\Delta_{y}$. We have also introduced six order terms with coefficients

[TABLE]

Few remarks are in order.

1. In the superconducting state $\alpha<0$, i.e. at temperatures lower that the critical temperature, $T<T_{c}$.

2. Coefficient at the fourth order term $\gamma_{1}-\gamma_{2}>0$, while $\gamma_{2}$ can be either positive or negative. When $\gamma_{2}$ is positive the ground state of the superconductor is nematic, when $\gamma_{2}$ is negative the state is chiral. Note that $\gamma_{2}$ changes sign from positive at $v_{z}\sim v$ to negative at $v_{z}=0$. The latter case describes doped thin film of the three-dimensional topological insulator.

3. One observes following relation between coefficients $0<\beta_{2}/\beta_{1}\leq 1/2$.

4. We do not include in-plane crystal anisotropy into the model, six order terms do not lift the degeneracy of the nematic state and since they are small compared to the fourth order terms we will ignore them.

5. To have more insight, consider two limiting cases. First, spherically symmetric case. We set $v=v_{z}$ and obtain

[TABLE]

Here $\xi=E_{+}(\mathbf{p})-\mu,\leavevmode\nobreak\ \leavevmode\nobreak\ E_{+}(\mathbf{p})=\sqrt{m^{2}+v^{2}p^{2}},$ and $(dp)=\frac{p^{2}dp}{2\pi^{2}}$. Second, consider cylindrical symmetric case in which we set $v_{z}=0$. We find that

[TABLE]

Here $\xi=E_{+}(\mathbf{p})-\mu,\leavevmode\nobreak\ \leavevmode\nobreak\ E_{+}(\mathbf{p})=\sqrt{m^{2}+v^{2}p^{2}_{\perp}}$, and $(dp)=\frac{p_{\perp}dp_{\perp}}{2\pi}\frac{dp_{z}}{2\pi}$. Notice that $\gamma_{2}$ and $\delta_{2}$ change sign, which signals for the phase transition between nematic and chiral states of the superconductor.

.5 Symmetric form of the Ginzburg-Landau free energy functional

It is worth to rewrite Ginzburg-Landau free energy functional in the symmetric form. Let us introduce new variables

[TABLE]

Thus, one has

[TABLE]

Useful identities

[TABLE]

Thus, we rewrite Ginzburg-Landau free energy functional in the symmetric form

[TABLE]

.6 Scaled form of the Ginzburg-Landau free energy functional

Finally we rewrite Ginzburg-Landau free energy functional in the scaled form for temperatures lower than the temperature of the superconductor-metal phase transition $T<T_{c}$. We perform standard substitutions:

[TABLE]

We define lengths, which for a single order parameter would correspond to the coherence and magnetic lengths

[TABLE]

One measures lengths in units of $l_{m}$: $R=r/l_{m}$. We also introduce

[TABLE]

such that magnetic field $\mathbf{B}=\bm{\nabla}_{R}\times\mathbf{a}$ is measured in units of $h_{0}$, together with the following rewriting

[TABLE]

where

[TABLE]

We obtain

[TABLE]

where we introduce new coefficients:

[TABLE]

They are not independent. First is the symmetric case where we set $v=v_{z}$ and obtain

[TABLE]

treating $\kappa$ as the only free parameter. Second, is the cylindrical case where we set $v_{z}=0$ and obtain

[TABLE]

again treating $\kappa$ as the only free parameter. Finally, including Maxwell term $({\bm{\nabla}}_{R}\times\mathbf{a})^{2}$ into free energy one obtains

[TABLE]