Half-quantum vortices on c-axis domain walls in chiral p-wave superconductors

Chiral superconductors are two-fold degenerate and domains of opposite chirality can form, separated by domain walls. There are indications of such domain formation in the quasi two-dimensional putative chiral $p$-wave superconductor Sr$_2$RuO$_4$, yet no experiment has explicitly resolved individual domains in this material. In this work, $c$-axis domain walls lying parallel to the layers in chiral $p$-wave superconductors are explored from a theoretical point of view. First, using both a phenomenological Ginzburg-Landau and a quasiclassical Bogoliubov-deGennes approach, a consistent qualitative description of the domain wall structure is obtained. While these domains are decoupled in the isotropic limit, there is a finite coupling in anisotropic systems and the domain wall can be treated as an effective Josephson junction. In the second part, the formation and structure of half-quantum vortices (HQV) on such $c$-axis domain walls are discussed.


Introduction
The unconventional superconductivity in the quasi two-dimensional layered perovskite Sr 2 RuO 4 , discovered more than two decades ago [1], continues to attract considerable interest. Following early predictions of spin-triplet superconductivity in connection with the strong Hund's coupling between the Ru d-orbitals [2,3], a large number of experiments have pointed toward an odd-parity [4,5], spin-triplet [6][7][8] and time reversal symmetry-breaking [9][10][11] superconducting state. This makes chiral p-wave pairing the most probable candidate [12][13][14], although challenging discrepancies remain [15]. In particular, recent NMR-Knight shift data show discrepancies to earlier results being more consistent with spin singlet pairing [16][17][18]. It is unclear so far whether spin-orbit coupling and the multi-orbital nature of the electronic band structure would make the results compatible with spin triplet symmetry [19,20].
The chiral p-wave state is two-fold degenerate with phases of positive and negative chirality ± connected through time reversal. If full rotation symmetry around the z-axis is present in the electronic structure, the corresponding Cooper pair states can be attributed a definite angular momentum ± . The formation of domains of opposite chirality is possible upon the nucleation of the superconducting phase and depends on the cooling process. However, despite extensive experimental investigation no direct observation of the domains has been reported such that their structure and size is unknown to date.
Two typical domain geometries should be considered for a layered superconductor with a quasi-2D electronic structure as illustrated in figure 1. The first is the ab-plane (inplane) domain wall which separates regions of opposite chiralities within the same layer. These domain walls support chiral quasiparticle modes indicative of the non-trivial topological nature of the chiral state [21,22]. The details of the domain wall structure, such as the stable configurations and the behavior of the two chiral components, have been studied previously within the Ginzburg-Landau (GL) theory [23][24][25], as well as within a quasiclassical approximation [26][27][28]. The existence of these domains was proposed to explain the spontaneous internal Figure 1. Illustration of the two typical domain wall geometries for a layered superconductor with a quasi-2D electronic structure. For the ab-plane (inplane) domain wall, the relevant perspective is a single layer, with translational invariance along the out-of-plane direction. The c-axis domain wall discussed in this paper spans the ab-plane and thus the relevant perspective is a vertical cut, with translational invariance along the other inplane direction. magnetic field observed in μSR experiments of Sr 2 RuO 4 [9,26]. Their effect on the magnetic properties of extended Josephson junctions has been discussed theoretically [29] and the domain wall motion was connected to the unusual interference patterns reported in junctions between Sr 2 RuO 4 and a conventional superconductor [30,31], in different setups based on the Sr 2 RuO 4 -Ru eutectic system [32,33], and in ring structures fabricated from single crystals [34]. However, no direct observation of domain formation has been reported in real space scanning probes [35][36][37]. Furthermore, the size of the domains inferred from the existing measurements seems inconsistent [38].
The other geometry is the c-axis domain wall which spans the ab-plane. This type of domain wall is less studied. However, it is energetically more favorable in comparison to the inplane domain walls owing to the weak interlayer coupling. Experimentally, it has been proposed that c-axis domains may explain the observed absence of the spontaneous flux at the surface [36].
In the present work we discuss the structure and magnetic properties of c-axis domain walls. In the isotropic limit (full rotation symmetry around the c-axis), such domain walls decouple the chiral domains completely from each other, because the angular momentum of the Cooper pairs is a good quantum number and restrictive selection rules apply in Cooper pair tunneling. Thus, a supercurrent along the c-axis through the domain wall would not be possible. By introducing anisotropy into the electronic band structure, however, a finite coupling appears, and supercurrents can flow between the chiral domains. Here we present a comprehensive discussion of the c-axis domain wall from both a phenomenological GL and a quasiclassical Bogoliubov-de Gennes (BdG) view point. Both descriptions confirm the finite coupling away from the isotropic limit, and find that the phase shift across the junction depends on the sign of the anisotropy and exhibits a non-trivial periodicity of π. Like Josephson junctions such domain walls can host vortices which due to the π-periodicity of the phase are half-quantum vortices (HQV). We show that the HQVs on the c-axis domain wall can indeed be stable, and that both the maximal current across the domain wall and the characteristic length scales of the HQV depend on the magnitude of the anisotropy. By tuning the system toward the isotropic limit, where the two domains are decoupled, the HQV dissolves along the domain wall. For stronger anisotropies, when the HQV becomes smaller than the relevant screening length, non-local effects have to be taken into account.
We note that the HQVs considered here belong to the class of fractional vortices within the framework of multicomponent order parameters [39] and are fundamentally different from the HQVs based on the spin rotation of the Cooper pairs in spin-triplet superconductors. These have attracted a lot of interest for their non-trivial topological properties [40][41][42]. Different directions of the associated defect [43] and the stability conditions [44] have been discussed theoretically, and they have been connected to the experimental observations of half-height magnetization steps [45,46] and Little-Parks oscillations [47] in specially shaped samples of Sr 2 RuO 4 .
In the following, first, the phenomenological GL approach is presented in section 2.1, providing both an approximative analytical and a self-consistent numerical solution. This is complemented by a quasiclassical BdG approach in section 2.2. Next, we explore the HQV on the domain wall. Its general structure is described in section 3.1, while the full junction phenomenology is analyzed in section 3.2 within a sine-Gordon framework for both the isotropic and the non-local limit.

Domain wall structure
We start with the analysis of the basic structure of c-axis domain walls in spin-triplet chiral p-wave superconductors for a systems with an anisotropic electronic structure. For a system with full rotation symmetry around the c-axis this state possesses a definite angular momentum L z = ±1 with a gap function d = Δ 0ẑ (k x ± ik y ). Much of the phenomenology of this pairing state can be transferred to other chiral superconducting phases with the same angular momentum property, such as the spin-singlet chiral d-wave state with the gap function ψ(k) = Δ 0 k z (k x ± ik y ). While we present here our discussion for the chiral p-wave state, all qualitative results also apply to related chiral states.

Phenomenological Ginzburg-Landau approach
The GL theory allows for a very efficient symmetry based approach to inhomogeneous structures of a superconducting order parameter. It will provide us with the essential ingredients for the study of the HQV in the second part of this paper, section 3.
The GL free energy of chiral p-wave superconductors is constructed from a two-component order parameter belonging to the two-dimensional irreducible representation E u of the full tetragonal point group D 4h , as used to describe the odd-parity spin-triplet pairing state given by d(r, k) =ẑ(η x (r)f x (k) + η y (r)f y (k)) in the d-vector notation, where by d ẑ corresponds to spin configuration of inplane equal-spin pairing [39]. Here { f x (k), f y (k)} are basis functions of E u , odd in k, and η = (η x , η y ) denotes the two-component order parameter which in the bulk takes the form η = η b (1, ±i) with the two chiralities ±.
The GL free energy functional is a scalar under all symmetry operations and, thus, given by [39] where V p is the superconducting region, T c the critical temperature, {a, b i , K i } the GL coefficients, A the vector potential, and D = ∇ − iγA the gauge-invariant derivative, where γ = 2π/Φ 0 with Φ 0 the flux quantum. We do not resolve the individual RuO 2 -layers along the c-axis as in a Lawrence-Doniach type of model of weak interlayer coupling [48][49][50], because in Sr 2 RuO 4 the coherence length ξ c along the c-direction is considerably longer than the interlayer distance s [12].

Anisotropy and GL coefficients
The GL expansion coefficients {a, b i , K i , K 5 } are material dependent parameters and can either be extracted from the corresponding microscopic Hamiltonian or experimental data. Their range is subject to the condition of stability of the free energy [39]. Relations between the sets a, {b i }, {K i } and K 5 can easily be determined from the linearized GL equations through the experimentally measured coherence lengths ξ ab and ξ c , the London penetration depths λ ab and λ c , and the GL parameter κ. This information is useful to write the free energy in dimensionless units such that the only experimental quantities entering are the superconducting anisotropy γ s = ξ ab /ξ c and the GL parameter κ ≡ κ ab = λ ab /ξ ab , which for Sr 2 RuO 4 are 20 and 2.6, respectively [13]. The ratios within the sets {b i } and {K i }, on the other hand, originate from further details of the electronic structure [51]. In a weak-coupling approach for the chiral p-wave state on a single band (e.g. γ-band of Sr 2 RuO 4 ) the coefficients b i are related through the band (Fermi surface) and gap structure as where · FS denotes the average over the Fermi surface. Introducing the anisotropy parameter of the electronic structure as [52] with ν ∈ ] − 1, 1 [, these further reduce to For the sake of definiteness we use f x (k) = v x (k) and f y (k) = v y (k), with v i the components of the Fermi velocity, such that the anisotropy of the gap function is identified with the anisotropy of the Fermi surface [28]. The constant b = 2b 1 + b 2 is chosen such that the amplitude of the bulk order parameter is given by The parameter ν is a measure for the anisotropy of superconducting phase as imposed by the electron band and gap structure [28]. In the isotropic limit, i.e. for a completely rotationally symmetric system with a cylindrical Fermi surface, the basis functions are f x (k) ∝ k x and f y (k) ∝ k y , resulting in ν = 0. Note that ν = ±1 corresponds to placing the free energy at the boundary of the stable region of the chiral p-wave state [53].
The analogous discussion for the coefficients of the gradient terms, {K i }, leads to K 1 = (3 + ν)K/4 and K 2 = K 3 = K 4 = (1 − ν)K/4 with the constant K = K 1 + K 2 . These inplane gradient coefficients expressed in terms of the anisotropy ν will only enter our discussion in the later part on the HQV, section 3. For the structure of the c-axis domain wall as discussed below, the order parameter is translationally invariant for inplane coordinates.

Phase shift across the domain wall
Domain walls involve a spatial change of the order parameter which can be decomposed into a variation of the amplitude of the order parameter components in the two domains and a shift of the overall phase across the domain wall. First we discuss which phase shift minimizes the free energy, depending on the anisotropy of the electronic structure, and then address the shape of the order parameter across the domain wall in the following section.
To analyze the domain wall structure it is convenient to use an order parameter representation η = (η + , η − ) = (η x − iη y )/2, (η x + iη y )/2 directly addressing the two degenerate chiral states. It is further useful to parameterize the complex order parameter components in terms of amplitude and phase enters the free energy, with z dw the position of the domain wall leading to an order parameter of the form η = (|η + |, |η − |e iϕ ), Here we neglect variations of the order parameter phase for |z − z dw | ξ , where one order parameter component vanishes, with ξ the width of the domain wall (we use z dw = 0). As the basic structure of the domain wall only depends on the out-of-plane z-direction, the inplane spatial dependence can be neglected completely.
It is immediately obvious that the free energy is minimized for ϕ = const., such that ∂ z ϕ = 0. Therefore we may model the c-axis domain wall using the boundary conditions with the bulk value η b given above. Straightforward symmetry considerations lead to the general structure η = (g(z), g(−z)e iϕ ), with g(z) being real, positive and asymmetric with respect to the domain wall The only remaining term depending on the phase shift ϕ is then Note that this term is concentrated on the domain wall region and vanishes for |z| ξ . The sign of the anisotropy of the electronic structure ν thus single-handedly determines the most energetically favorable choice of ϕ across the domain wall, The resulting free energy and, thus, the shape of the order parameter, on the other hand, only depend on the magnitude of the absolute value |ν|.
We would like to comment here on our approximation to reach the solution for ϕ in equation (8). First, we assumed that the coherence length along the c-axis covers many layers. If this is not the case and K 5 would be very small, the approximation, ∂ z ϕ = 0 is not justified. Then passing with ν through zero may not result in a discontinuous change from 0 to π/2. Rather ϕ could move continuously between these two minima for a certain range of ν around 0, because the ϕ-dependence of the free energy would be rather weak. Indeed our numerical treatment points toward such a behavior. However, for the sake of simplicity we do not analyze this rather special limit here, as it lies outside our scope.
In the isotropic limit (ν = 0), the free energy is fully degenerate for all ϕ, signaling the complete phase decoupling of the two domains, analogous to the Josephson coupling between two condensates. This limit describes a system with complete rotation symmetry around the c-axis, where the Cooper pair orbital angular momentum is a good quantum number and is conserved during the tunneling event. Pair tunneling is, thus, prohibited by this selection rule, consistent with our phenomenological result. Once the rotational symmetry is broken the selection rule is no longer valid and Cooper pairs can be transferred, in the present case on the level of a second order coupling, i.e. even numbers of Cooper pairs pass together. In this way the two domains are phase coupled and a supercurrent can flow across the domain wall. Given the non-trivial π-periodicity of the discrete allowed values for ϕ, realized through the cos(2ϕ) term in the free energy, we will explore the possibility of HQVs on the domain wall in section 3.

Shape of the order parameter across the domain wall
While the phase shift across the domain wall only depends on the sign of the anisotropy, the change in the amplitude of the order parameter components only depends on its magnitude |ν|. We will now tackle the domain wall problem first with an approximate variational solution, which we then compare with the exact numerical solution for the shape of the order parameter.
Following the approximations introduced in reference [25], the total amplitude of the order parameter is kept constant everywhere |η| 2 = |η b | 2 with the ansatz (η + , η − ) = η b sin(χ(z)), cos(χ(z))e iϕ . The boundary condition from equation (6) translates into χ(z → ∞) → 0 and χ(z → −∞) → π/2. In this way the free energy simplifies to where the phase shift ϕ given in equation (8) has already been implemented. Therefore this free energy only depends on the absolute value of the anisotropy |ν| and not on its sign. Variation with respect to χ leads to the sine-Gordon type differential equation Using ξ c = K 5 /(bη 2 b ), the standard solution is such that the width of the domain wall scales as ξ = 2ξ c / 1 − |ν|. The free energy density per unit inplane area is then found by inserting this solution, where F 0 is the bulk free energy without the domain wall. The cost of the domain wall is largest in the isotropic limit |ν| = 0, where the two domains are fully phase decoupled. On the other hand, at the stability boundary of the chiral p-phase, ν = ±1, the domain wall energy vanishes and the width ξ diverges. Compared to the inplane domain wall whose width is connected with ξ ab , the c-axis domain wall scales with ξ c ( ξ ab ) and is energetically much cheaper.
We minimize the free energy numerically setting T = 0 and without any restrictions on the amplitude of the order parameter. In the dimensionless formulation of the free energy the only remaining parameter here is the superconducting anisotropy γ s for which we use the value for Sr 2 RuO 4 mentioned above. Fixed values are used for ϕ, however, according to equation (8), ϕ = {0, π/2}. The numerical scheme follows a one-step relaxed Newton-Jacobi method for boundary value problems [54][55][56], details see reference [51].
The numerical results for the two components |η + | (dashed) and |η − | (dotted) are shown in figure 2 together with the total absolute value |η| (solid curve), all at ν = 0. In addition, the results for higher values  of ν are indicated by the thin lines which display the growing domain wall width for increasing |ν|. Unlike in our variational ansatz the value of |η| shows a dip at the domain wall which is weakened as |ν| grows and is entirely constant in the limit |ν| → 1. Finally, we note that the symmetry (η + , η − ) ∝ g(z), g(−z) is borne out in the numerical solution.
The numerical result for the free energy density is shown in figure 3 for ϕ = {0, π/2} (empty and filled dots), together with the approximative analytical solution from equation (12) (solid black line), as a function of the anisotropy ν. In addition, a linear fit to the numerical data for ν ≈ 0 is indicated (dashed line), which will be used in section 3.2. The crossing of the two energy branches at ν = 0 indicates a first order change for ϕ between the two sectors of different values as given in equation (8). The qualitative agreement between the numerical solution and the variational approximation is very good. Close to |ν| → 1 even a quantitative agreement is found, as in this limit the dip in |η| disappears in the numerical solution, in accordance to the simplifying assumption for the variational ansatz above. Finally we note the symmetry F [ν, ϕ] = F [−ν, ϕ + π/2] always holds.

Quasiclassical Bogoliubov-deGennes approach
In this section, we address the c-axis domain wall from a more microscopic viewpoint through a self-consistent BdG treatment. For simplicity, we restrict to a one-band spinless Fermion tight-binding model on a square lattice which is sufficient for a spin-triplet superconductor. The conclusions are, however, applicable to spinful and multi-band systems. In contrast to the GL analysis above, here the vector potential is neglected. Nevertheless, the solutions obtained show the same main behavior and give an insight on the role of the quasiparticle states at the domain wall. To simulate a c-axis domain wall in numerical BdG, we also use inplane translational invariance in a system of N z layers. The corresponding mean-field BdG formulation is given by where c † l,k (c l,k ) creates (annihilates) an electron in the lth layer with inplane momentum k and ξ k = −2t(cos k x + cos k y ) − 4t cos k x cos k y − μ denotes the inplane dispersion relation. The last two terms represent the superconducting pairing with chiral p-wave symmetry, η l,k = η x (l)f x (k) ± iη y (l)f y (k), where we will use two types of pairing states assuming nearest neighbor [( f x , f y ) = (sin k x , sin k y )] or next-nearest neighbor [( f x , f y ) = (sin k x cos k y , sin k y cos k x )] pairing interactions. Note the quasi-two-dimensional band structure of Sr 2 RuO 4 implies that the dispersion along the c-axis is very small, t z t, t . Without loss of generality, t = 1, t = 0.375t and t z = 0.03t are used throughout this section.
The c-axis domain wall is formed by allowing the two chiral components η x and η y to vary from layer to layer. Far from the domain wall the pairing states on the two sides correspond to the bulk states of opposite chirality, with an additional phase shift ϕ imposed between the two domains and Δ 0 the self-consistently determined bulk value of gap function amplitude. To perform the self-consistent calculation of the bulk and the domain wall structure for the two cases of pairing states we introduce either a purely nearest-neighbor or a purely next-nearest-neighbor pairing interaction. We start the iteration process for self-consistency with an initial configuration for a domain wall at the center of the system (l = 0), following reference [25], where χ(l) = π 4 (1 + tanh l λ ), with λ being a constant that defines the initial value for the extension of the domain wall along z. The structure of this initial configuration may or may not correspond to an actual energetically favorable solution. In the case of the latter, the system in general evolves into a stable state after sufficient iterative steps in our self-consistence procedure. For the calculations in figure 5, λ = 2 was used in the initial configuration.

Stable domain wall configurations
In lattice models, the anisotropy parameter ν is generically non-vanishing and shall depend on the details of the gap and band structure anisotropy. In figure 4 we plot the relevant parameters as a function of chemical potential for two different chiral p-wave gap functions on a square lattice. Consistent with the GL treatment, our variational BdG yields two sectors of stable domain walls, depending on the sign of ν.
Next we would like to understand the role of the phase shift ϕ. For this purpose we use a non-self-consistent approach of the BdG scheme where the gap function is taken with the fixed χ(l) as given in equation (15) while keeping ϕ as a free parameter. In this way we deduce the energy of the domain wall as a function of ϕ, relative to the bulk condensation energy which is approximated by per layer. Figure 5 displays the energies for the types of gap functions introduced above (NN-pairing in the upper and NNN-pairing in the lower panel) for varying chemical potential. The anisotropy ν corresponding to the different curves connect with the black dots given in figure 4. The minima are found at ϕ = 0, π for negative ν and at ϕ = π/2 (3π/2) for positive ν. The smaller |ν| the weaker the ϕ-dependence, as expected from our previous discussion.

Domain wall Andreev quasiparticle states
The emergence of subgap chiral quasiparticles at the ab-plane domain walls is a well-known feature whose origin lies in the topological nature of the superconducting phase [22,26]. The c-axis domain walls also host localized subgap quasiparticle states which, in contrast, do not directly reflect topological properties. Nevertheless, the chiral pairing does have its impact on the spectrum of these subgap Andreev states and their contribution to the coupling of the two domains. For our analysis we consider the situation from the viewpoint of a planar c-axis Josephson junction where the two connected superconductors are chiral, for simplicity described by a gap function on the Fermi surface of the form η k,± = η 0 e ±iθ k . We examine the electronic states based on a quasiclassical approach, focusing on particle trajectories between the two superconductors (see figure 6). Such a trajectory can be labeled by the momentum k = (k cos θ k , k sin θ k , k z ), near the Fermi surface. Following figure 6 the trajectory connects gap functions of the phase ϕ/2 − θ k in η − -domain with those of −ϕ/2 + θ k in the η + -domain. Consequently, for such a trajectory the corresponding phase difference between the two domains is ϕ − 2θ k . The set of quasiparticle states derived from each trajectory (including the reflected part) contributes to the coupling and the associated energy depends on the phase difference, E(ϕ − 2θ k , θ k ) and, in general, on the direction θ k . The total energy of the junction is the integral over all trajectories, here reduced to the angles θ k , as we restrict to momenta at the Fermi surface, Together with the standard periodicity E tot (ϕ) = E tot (ϕ + 2πn) this corresponds well to cos(2ϕ) in lowest order coupling. This simplified viewpoint allows us now to give a qualitative discussion of the role of anisotropy. For full rotation symmetry, E(φ, θ) = E(φ, θ + α) for an arbitrary angle α. Thus, E(ϕ − 2(θ + α), θ + α) = E(ϕ − 2θ, θ) with ϕ = ϕ − 2α leads in equation (17) immediately to E tot (ϕ) = E tot (ϕ ) such that the total junction energy is independent of the phase shift ϕ and the two chiral states are phase decoupled.
On the other hand, if we assume a four-fold rotation symmetry, then above relations for the energy are only true for α = π 2 (2n + 1). From this we find which is borne out to lowest order by a dependence like cos(2ϕ) as found previously in the GL formulation and is also consistent with the numerical result for the domain wall energy in figure 5. The absence of phase coupling for the isotropic case means that the critical current through a domain wall along the c-axis would vanish and would even be rather small for the anisotropic system.

Half-quantum vortex
We now turn to the question of vortices on a c-axis domain wall. As we will show, the π-periodicity of the phase shift ϕ suggests the existence of HQV. The structure and magnetic properties of a single HQV on the domain wall will be discussed in section 3.1. In section 3.2, the phenomenology of treating the c-axis domain wall as an effective Josephson junction is presented. The c-axis domain wall considered in this work supports such fractional vortices, but limited to carrying (integer multiples of) half of a standard flux quantum Φ 0 = hc/2e. This results from the non-trivial π-periodicity of the phase shift across the domain wall described in equation (8) and from the underlying cos(2ϕ) coupling term in equation (7), such that a π-kink is the smallest possible increase of the phase shift between two stable configurations of the domain wall, unlike the standard Josephson vortices corresponding to a 2π-kink of ϕ.

Stability, phase, and shape of the order parameter
We consider the case of ν < 0 where the stable domain walls possess the phase shifts ϕ = nπ. The c-axis domain wall shall be centered at z = 0 separating the phases η + for z < 0 and η − for z > 0. We introduce now a line defect on the domain wall at x = 0 by choosing the phase shift ϕ = 0 for x → −∞ and ϕ = π for x → +∞. As illustrated in figure 7, ϕ has to change by ±π/2 for z ≷ 0 along the x-axis to connect these two stable domain wall configurations. Away from the line defect, such a phase gradient does not cost any energy but can be absorbed by the proper gauge γA x (x) = ∂ x ϕ(x). The resulting states of the order parameter (η + , η − ) serving as the boundary conditions are all indicated in figure 7. In addition, the order parameter phase φ x is shown as a density plot in the background, for which we switch back to the basis, η = (η x , η y ) = (|η x |e iφ x , |η y |e iφ y ). We also note that the system is still translationally invariant along the y-direction.
For the subsequent analysis we neglect any surface effects, that is, we consider an infinite sample. Moreover, for simplicity we assume that to lowest order A y = 0, even though this prevents a fully self-consistent analysis because the spatial variation of the order parameter along the x-direction is associated with a small but finite A y component through the self-screening of the induced source current J y along the y-direction. The detailed structure of the magnetic flux pattern around the line defect is computed by minimizing the full GL free energy functional for the boundary conditions indicated in figure 7 numerically using a relaxed on-step Newton-Jacobi method, described in detail in reference [51]. To facilitate the computations, a less extreme value of γ s = 10 is used in this part, while the value of γ s = 20 as found in the literature for Sr 2 RuO 4 [13] was used in the first part.
Various quantities extracted from the computational result for the order parameter are shown in figure 8, with the position of the line defect at (x, z) = (0, 0). The relative phase between the x-and the y-component of the order parameter is shown in the top left panel. It is +π/2 for z > 0 and −π/2 for z < 0, in accordance with the setup and consistent with the presence of the domain wall. The absolute value |η| 2 of the order parameter is shown in the top right panel. As discussed above, it is suppressed at the domain wall. Now, it is additionally reduced at the line defect. The phase φ x = arg(η x ) is shown in the bottom left panel and behaves as anticipated in figure 7. Finally, the panel at the bottom right shows the π-kink in the phase shift ϕ(x) along the domain wall, defined here through ϕ = φ − − φ + . This quantity is the analogue to the Josephson phase when treating the domain wall as an effective Josephson junction.  The line defect is characterized by the winding of one of two order parameter components. In our case this can be extracted by representing the order parameter as The behavior of these two components is illustrated in figure 9 with a three-dimensional plot of their amplitude and a colored density plot of their phase beneath. While η 1 vanishes at the line defect, η 2 remains finite everywhere. Taking the overall phase structure into account, an unusual vortex carrying half a flux quantum (HQV) emerges from this line defect, as we explicitly show below.

Characteristic length scales and magnetic properties
Perpendicular to the domain wall, the magnetic field of the HQV is screened efficiently by inplane currents on the length scale λ ab . The extension of the magnetic flux distribution along the domain wall, however, depends on the coupling between the two domains and, thus, on the supercurrent which can flow across the domain wall. In order to understand the behavior of the domain wall it is helpful to view it as an effective Josephson junction. The critical current scales like J c ∝ |ν| for small |ν|, as we will point out in equation (24). The extension of the HQV, like the Josephson vortex, corresponds to the Josephson penetration depth which scales as λ J ∝ 1/ √ J c ∝ 1/ |ν|. Since the critical current vanishes in the isotropic limit, we also expect that the HQV will dissolve for ν → 0. On the other hand, for growing |ν| the flux distribution along the x-axis shrinks and the picture of the Josephson vortex is not entirely appropriate anymore as new effects come into play.
For layered superconductors there exist two different screening lengths λ ab and λ c = γ s λ ab λ ab due to screening currents parallel and perpendicular to the layers, respectively, in analogy to the anisotropy of the coherence length γ s = ξ ab /ξ c . While for conventional Josephson junctions usually λ J λ London , for the situation considered here, the two relevant length scales can be comparable in size already at very small values of ν ≈ 4%, see equation (30). Once λ J < λ c , non-local magnetic properties of the Josephson junction have to be considered, as reviewed in reference [57]. In this case the long-range screening behavior is more like that of an Abrikosov vortex on the length scale λ c , while the core remains Josephson-like, but has a new characteristic length l = λ 2 J /λ c < λ J [58]. This is derived in detail in section 3.2, while below we describe how to extract these characteristic length scales from the computational results for the HQV.
The structure of the magnetic flux line and the current pattern circulating around the HQV are shown in figure 10. A density plot of the magnetic field B y (x, z) is displayed at the top, with an inset zooming in on its center, where a vector plot of the current (J x , J z ) can be seen. A measure of the extension of the HQV along the x-axis can be derived by limiting the integral for the magnetic flux through the boundaries at x = ±w, whose result is shown in the bottom left panel. For w → ∞ we observe that the flux Φ(w) saturates at Φ 0 /2 as expected for an HQV. We now use this behavior to define the length w box through Φ(w box ) = 0.49Φ 0 (this somewhat arbitrary cutoff does not qualitatively influence the final results), as indicated in figure 10. Furthermore, the core size of the HQV can be estimated using the profile of the current across the domain wall, J z (x, 0), shown in the bottom right panel. The length scale of the core is defined as the position x max of the maximal current J max z , which in turn gives a measure for the critical current of the domain wall considering it as an effective junction between the two domains. In addition, −J z (−x) is shown (dashed), which indicates that the current is slightly asymmetric, while the total current ∞ −∞ dxJ z (x, 0) still integrates to zero within the numerical accuracy. From a detailed analysis of the full GL free energy functional, and also from symmetry considerations, it becomes apparent that for the spatial variation of the phase shift ϕ(x) = −ϕ(−x) + π. Specifically, the relaxation away from the HQV toward the stable 0 (or 2π) phase shift occurs on a slightly different length-scale than toward the π phase shift, because the inplane gradient terms are fundamentally 2π-periodic only. This small effect, however, does not affect our overall discussion.

Junction phenomenology
We now explore the behavior of the HQV for varying anisotropy ν by treating the domain wall as an effective Josephson junction. First, the current-phase relation and the critical current are discussed using our variational approach and comparing it to the computational result. Next, a sine-Gordon model is formulated both for the isotropic (small ν) limit and for the non-local (large ν) limit. Eventually, the characteristic length scales w box and x max for the core size and the full size of the HQV as introduced above are analyzed, based on the analytical estimates from the two limits, and compared to the values extracted from the computational results.

Current-phase relation
Treating the domain wall as an effective Josephson junction, the Josephson phase difference is associated with the phase shift ϕ(x) = φ − − φ + . Based on the π-periodicity of the phase shift, the current-phase relation behaves to lowest order as J z (ϕ) = J c sin(2ϕ), with the critical current a parameter to be determined, either from an analytical approach (see below), or extracted from our computational results (see figure 10). Deriving the current-phase relation from the full domain wall free energy self-consistently is beyond the scope of this paper. Instead we focus on the limit of small anisotropies and use the approximation presented in the first part, equation (12), assuming a constant total order parameter amplitude |η| throughout the system and neglecting the inplane gradient coupling terms. The phase-dependent domain wall free energy is then given by In analogy to the Josephson junction, the current density across the domain wall is therefore Approximating for small anisotropies |ν| 1, the lowest order expression proposed in equation (21) is confirmed, with the critical current given by recovering the behavior for the isotropic limit discussed above. When the domains are decoupled, no supercurrent can flow across the domain wall, and indeed J c (ν = 0) = 0.
The assumption of a constant total amplitude |η| 2 = |η b | 2 was found to be the least valid for small ν, while this is exactly the limit of interest here. We therefore propose a semi-analytical model, where instead of using equation (22) for the domain wall energy directly, the computational result near ν ≈ 0 is fitted linearly, indicated in figure 3 by the dashed line. The proportionality J c ∝ |ν| still holds, but with a different slope, J fit c ≈ 0.5J c . To put the critical current of the domain wall in relation to a situation without any domain wall, i.e. the upper limit, the standard procedure for the depairing current is followed [59]. The maximal current along the c-direction is then given by Figure 11. Maximum current across the domain wall at the HQV extracted from the computational results (black dots), compared to the critical current obtained from the approximative analytical solution (black line) and from a linear fit to the numerical solution for the free energy (orange). The ultimate limit is the depairing current (dashed). The inset zooms into the limit of small anisotropies.  Position of the maximum current across the domain wall at the HQV, extracted from the numerical results (black dots), compared toλ J (orange) for the isotropic limit |ν| ≈ 0 and to l (black) for the non-local limit. The fundamental screening lengthλ c (dashed) and the crossover anisotropy ν cross are also indicated.
approximated by an Abrikosov type of vortex of an anisotropic superconductor [48,62], whose field is given by ignoring the core region. Figure 12 shows the value of w box extracted from the numerical results (black dots) and the limiting w A computed from the above expression for the Abrikosov vortex (dashed line), where only the screening lengths λ ab and λ c enter. We find that for |ν| ≈ 0 the HQV expands as predicted, indicating that the critical current vanishes in the isotropic limit. On the other hand, for growing |ν| the vortex size is well described by the long-distance behavior of an Abrikosov-type vortex. Let us turn to the core size x max of the HQV, which we defined as the position of the maximum current across the domain wall (see figure 10). In figure 13 the computational results for x max (black dots) are shown, together with the positionλ J defined via the Josephson penetration depth (λ J = arcsinh(1)λ J ) (orange line), and the characteristic length of the non-local limit l (black line). For both, the critical current used is extracted from J max z , and inserted into equations (29) and (35), respectively. In addition, the relevant screening lengthλ c = arcsinh(1)λ c is displayed (dashed line), as well as the value ν cross of the anisotropy at whichλ J =λ c . The size of the vortex core in the isotropic limit |ν| ≈ 0 is indeed given by λ J , while it is determined by l for larger anisotropies, and changes in the crossover region around ν cross .
The variation of the extension of both the magnetic flux distribution and the core size are captured well within our sine-Gordon models. Restricting to lowest order current-phase relation gives consistent results valid even beyond the crossover region. Only at larger anisotropies bulk effects start to interfere and the HQV develops a normal core. While with increasing |ν| the order parameter amplitude |η| becomes less reduced at the domain wall (see figure 3), it actually shrinks more strongly at the center of the HQV. In this