Chiral p-wave superconductors have complex coherence and magnetic field penetration lengths

TL;DR

This paper demonstrates that chiral p-wave superconductors with broken time reversal symmetry exhibit coupled magnetic and order parameter modes, leading to complex, hybridized decay lengths and oscillatory magnetic field behavior.

Contribution

It reveals that in anisotropic, time-reversal-breaking superconductors, magnetic and order parameter excitations are mixed, resulting in unique decay modes not characterized by traditional lengths.

Findings

Magnetic and order parameter modes are coupled in p+ip superconductors.

Decay lengths are hybridized and can be complex, causing oscillatory decay.

Normal modes have their own distinct decay lengths, differing from conventional superconductors.

Abstract

We show that in superconductors that break time reversal symmetry and have anisotropy, such as p+ip materials, all order parameters and magnetic modes are mixed. Excitation of the gap fields produces an excitation of the magnetic field and vice versa. Correspondingly the long-range decay of the magnetic field and order parameter are in general given by the same exponent. Thus one cannot characterize p+ip superconductors by the usual coherence and magnetic field penetration lengths. Instead the system has normal modes that are associated with linear combinations of magnetic fields, moduli of and phases of the order parameter components. Each such normal mode has its own decay length that plays the role of a hybridized coherence/magnetic field penetration length. On a large part of the parameter space these exponents are complex. Therefore the system in general has damped oscillatory decay of the magnetic field accompanied by damped oscillatory variation of the order parameter fields.

Equations

$$\mathcal{F}=\frac{1}{2}{Q}_{ij}^{\alpha \beta }{D}_i{\psi }_{\alpha }\overline{D_j{\psi }_{\beta }}+\frac{1}{2}{B}^2+F_p$$

$${\psi }_{\alpha }={\rho }_{\alpha }{e}^{i{\theta }_{\alpha }}\alpha =1,2$$

$$F_p=V({\rho }_1,{\rho }_2)+\frac{\eta }{8}{\rho }_1^2{\rho }_2^2\cos 2{\theta }_{12}$$

$$Q_{ij}^{\alpha \beta }{D}_i{D}_j{\psi }_{\beta }$$

$${\partial }_j({\partial }_j{A}_i-{\partial }_i{A}_j)$$

$$J_i:=\mathrm{Im}(Q_{ij}^{\alpha \beta }{\overline{\psi }}_{\alpha }{D}_j{\psi }_{\beta }).$$

$${\theta }_{\Sigma }:=\frac{1}{2}({\theta }_1+{\theta }_2).$$

$$p_i:=A_i-{\partial }_i{\theta }_{\Sigma }$$

$${\epsilon }_{\alpha }:={\rho }_{\alpha }-u_{\alpha },{\theta }_{\Delta }:=\frac{1}{2}({\theta }_{12}-{\theta }_0).$$

$$D_i{\psi }_1$$

$$D_i{\psi }_2$$

$${\partial }_j({\partial }_j{p}_i-{\partial }_i{p}_j)$$

$$=-Q_{ij}^{11}{u}_1^2(p_j-{\partial }_j{\theta }_{\Delta })-Q_{ij}^{22}{u}_2^2(p_j+{\partial }_j{\theta }_{\Delta })$$

$$-u_1{u}_2\cos {\theta }_0\{Q_{ij}^{12}(p_j+{\partial }_j{\theta }_{\Delta })-Q_{ij}^{21}(p_j-{\partial }_j{\theta }_{\Delta })\}$$

$$-\sin {\theta }_0\{Q_{ij}^{12}{u}_1{\partial }_j{\epsilon }_2-Q_{ij}^{21}{u}_2{\partial }_j{\epsilon }_1\},$$

$$-\left({\partial }_1^2+{\partial }_2^2\right)p_i+{\partial }_i{\partial }_j{p}_j-L_{ij}{\partial }_j{\theta }_{\Delta }+K_{ij}{p}_j$$

$$\pm Q_{ij}^{12}\left(u_1{\partial }_j{\epsilon }_2-u_2{\partial }_j{\epsilon }_1\right)$$

$$\pm Q_{ij}^{12}\left(u_2{\partial }_i{\partial }_j{\epsilon }_1+u_1{\partial }_i{\partial }_j{\epsilon }_2\right)-K_{ij}{\partial }_i{\partial }_j{\theta }_{\Delta }$$

$$+L_{ij}{\partial }_i{p}_j+2\eta u_1^2{u}_2^2{\theta }_{\Delta }$$

$$-Q_{ij}^{11}{\partial }_i{\partial }_j{\epsilon }_1\pm Q_{ij}^{12}{u}_2{\partial }_i\left(p_j+{\partial }_j{\theta }_{\Delta }\right)+{\mathcal{H}}_{1\beta }{\epsilon }_{\beta }$$

$$-Q_{ij}^{22}{\partial }_i{\partial }_j{\epsilon }_2\mp Q_{ij}^{12}{u}_1{\partial }_i\left(p_j-{\partial }_j{\theta }_{\Delta }\right)+{\mathcal{H}}_{2\beta }{\epsilon }_{\beta }$$

$$K_{ij}=Q_{ij}^{11}{u}_1^2+Q_{ij}^{22}{u}_2^2,$$

$$L_{ij}=Q_{ij}^{11}{u}_1^2-Q_{ij}^{22}{u}_2^2,$$

$${\mathcal{H}}_{\alpha \beta }={\frac{{\partial }^2{F}_p}{\partial {\rho }_{\alpha }\partial {\rho }_{\beta }}|}_{\left(u_1,u_2,\pm \frac{\pi }{2}\right)}$$

$$p_i$$

$${\theta }_{\Delta }$$

$$b=\frac{-\mbox{\boldmath{$n$}}}{\mbox{\boldmath{$n$}}\cdot K\mbox{\boldmath{$n$}}}\cdot\left(K\mbox{\boldmath{$n$}}^{\perp}\,a-L\mbox{\boldmath{$n$}}\,\theta_{\Delta}^{\prime}\pm Q^{12}\mbox{\boldmath{$n$}}(u_{1}\varepsilon_{2}^{\prime}-u_{2}\varepsilon_{1}^{\prime})\right)$$

$$\mathcal{A}{\vec{w}}^{\prime \prime }+\mathcal{B}{\vec{w}}^{\prime }+\mathcal{C}\vec{w}=0,$$

$$\vec{w}(X)=\sum _{i=1}^8{c}_i{\vec{v}}_i{e}^{-{\mu }_{i}X}$$

$$\det \left({\mu }^2\mathcal{A}-\mu \mathcal{B}+\mathcal{C}\right)=0.$$

$${\lambda }_i=\frac{1}{{\mu }_i}.$$

$$({\epsilon }_1,{\epsilon }_2,{\theta }_{\Delta },a)\sim c{\vec{v}}_{r,\ast }{e}^{-\mathrm{Re}({\mu }_{\ast })X}\cos \mathrm{Im}({\mu }_{\ast })X,(X\to \infty )$$

$$\displaystyle Q^{11}=\left(\begin{array}[]{cc}3+\nu&0\\ 0&1-\nu\end{array}\right),\quad Q^{22}=\left(\begin{array}[]{cc}1-\nu&0\\ 0&3+\nu\end{array}\right),$$

$$\displaystyle\quad Q^{12}=\left(\begin{array}[]{cc}0&1-\nu\\ 1-\nu&0\end{array}\right)$$

$$F_p$$

$$G={\int }_{\Omega }\mathcal{F}-H{\int }_{\Omega }B+{\int }_{\partial \Omega }{\mathcal{F}}_{\text{mbox}surf}$$

$${\mathcal{F}}_{\text{mbox}surf}$$

Full text

I Introduction

Superconducting states that spontaneously break time reversal symmetry (BTRS) are a subject of intense experimental pursuit. Two types of BTRS state that attract particular interest are chiral $p$-wave superconductors where the most intense discussions were focused on $Sr_{2}RuO_{4}$ Mackenzie (2017); Mackenzie and Maeno (2003), and $s+is$ or $s+id$ superconducting states, evidence for which was recently found in iron-based superconductors Grinenko et al. (2017, 2018). BTRS states are described by an order parameter that has at least two components, because they break at least $U(1)\times Z_{2}$ symmetry. Also, rather generically, there is anisotropy in such superconducting states. In this work we investigate the most basic property of that state: the magnetic field penetration and coherence lengths.

The basic fundamental length scales of superconductors were first discussed by Fritz and Heinz London London and London (1935) and Ginzburg and Landau Landau and Ginzburg (1950) in an ordinary superconductor. This was done in the model for the simplest superconductor that breaks $U(1)$ symmetry, described by a single complex field $|\Psi|$ neglecting crystal anisotropies. The London magnetic field penetration length $\lambda$ is the power in the exponential law of decay of the magnetic field: ${\bf B}={\bf B}_{0}e^{-r/\lambda}$. The coherence length $\xi$ is the scale associated with the exponential law describing how the modulus $|\Psi|$ of the complex field, describing the order parameter, restores its ground state value $\bar{|\Psi|}$ away from a perturbation: $|\Psi|{(r)}\approx\bar{|\Psi|}-\rm{const}\ e^{-r/\xi}$. The microscopic BCS theory of superconductivity related the modulus of the order parameter field $|\Psi|$ to a superconducting gap $\Delta$ in the single-electron spectrum. The definition of the coherence length in the context of superconductivity has an extra factor $\sqrt{2}$ which we absorb for brevity in the definition of $\xi$. Often the coherence length is assessed only approximately, and it is important to remember the limitations of these approximate definitions. For example while in the simplest Ginzburg-Landau model coherence length is often estimated via vortex core size or slope of the order parameter near the centre of the vortex core, such estimates are known to fail even in the simplest models at low temperates Gygi and Schlüter (1991). Another indirect way to assess coherence length assumes its inverse proportionality to the gap function $\Delta$ in the BCS expression $\xi_{0}\propto 1/\Delta$. Likewise this expression has very limited validity. It cannot serve as an estimate at strong coupling or in the multi-component case. For example in the case of several gaps that would give unphysical divergence of coherence length where a gap is closing (i.e. at the crossover from $s_{++}$ to $s_{\pm}$ states where all coherence lengths should be finite because there is no symmetry breaking and no accidental degeneracies). Similarly that estimate would miss the divergence of coherence length when a superconductor transitions from ordinary to BTRS state i.e. $s$ to $s+is$ or $s$ to $s+id$ state, whose existence is dictated by symmetry. These examples shows that accurate coherence length calculations are required while simple estimates can be highly misleading. Calculations of coherence and magnetic field penetration lengths have been made for isotropic multicomponent models for general interactions both in phenomenological and microscopic models Babaev et al. (2010); Carlström et al. (2011a); Silaev and Babaev (2011); Carlström et al. (2011b); Garaud et al. (2018). The multi-component nature of these systems strongly affects only the coherence lengths, while the magnetic field penetration length is merely renormalized by intercomponent couplings. The situation was found to be very different in $U(1)$ multiband superconductors if different bands have different anisotropies. While usually the magnetic (London) modes decouple from other normal modes of the system, such as density and phase difference (Leggett) modes, having different anisotropies in different bands results in a hybridization of the London mode with the phase difference mode Silaev et al. (2018); Winyard et al. (2019a, b). For a system with $N$ bands that means that magnetic field decay is described by several modes with different exponents and there could be up $N+1$ such modes in the systems considered in Silaev et al. (2018); Winyard et al. (2019a, b). Furthermore the powers in the corresponding exponents under certain conditions are complex leading to a damped oscillatory decay of the magnetic field.

That raises the question: what is the behaviour of the magnetic field and what are the coherence lengths in $p+ip$ superconductors, since such systems are inherently both multicomponent and anisotropic? The important difference with the systems considered in Silaev et al. (2018); Winyard et al. (2019a, b), as discussed below, is the fact that such a superconducting state has spontaneously broken time reversal symmetry.

The standard Ginzburg-Landau model for a $p+ip$ superconductor can we written in dimensionless units as

[TABLE]

where the greek indices enumerate components of the order parameter and latin indices stand for space directions. Summation over repeated indices is implied, $D_{i}=\partial_{i}-iA_{i}$ is the covariant derivative with the gauge field $A_{i}$, and the complex fields

[TABLE]

represent the different superconducting components. We consider here a quasi-two-dimensional system or a configuration of a three dimensional system that is translation invariant in the $z$ direction. The magnetic field $\bm{B}=(0,0,B)=(0,0,\partial_{1}A_{2}-\partial_{2}A_{1})$ is directed so that the spatial indices take only the values $1,2$. $F_{p}$ represents the potential terms which, by gauge invariance, may depend only on $\rho_{\alpha}$ and $\theta_{12}:=\theta_{1}-\theta_{2}$. For the standard $p+ip$ superconductor, the $\theta_{12}$ dependence enters only via a term of the form $(\psi_{1}\psi_{2}^{*})^{2}+c.c.$, Heeb and Agterberg (1999); Agterberg (1998); Vadimov and Silaev (2013) that is,

[TABLE]

with $\eta>0$. Then the ground states (minima) of $F_{p}$ are degenerate occurring with $\theta_{12}=\pm\pi/2$. The ground state of this system is not gauge equivalent to its complex conjugate. Hence the system exhibits broken time reversal symmetry. Note that, although our focus below will be on the example of $p+ip$ superconductors, the model is very general, also describing other BTRS states such as $s+is$ and $s+id$ superconductors Garaud et al. (2017); Vadimov and Silaev (2018). Our results obtained below apply also to such states when anisotropy is present.

The anisotropy of the system enters through the parameters $Q^{\alpha\beta}_{ij}$, which must satisfy $Q^{\beta\alpha}_{ji}=(Q^{\alpha\beta}_{ij})^{*}$ to ensure $F$ is real. Henceforth we assume, as is standard, that all $Q^{\alpha\beta}_{ij}$ are real.

II Calculation of length scales

The spatial dependence of the fields at equilibrium is governed by the Ginzburg-Landau (Euler-Lagrange) equations for the functional $F=\int_{{\mathbb{R}}^{2}}{\cal F}$,

[TABLE]

where the total supercurrent is

[TABLE]

Consider the behaviour of the system a long distance from some defect (e.g. a vortex, domain wall or material boundary). Since the fields are close to their ground state values, they should be well approximated by solutions of the linearization of the Euler-Lagrange equations about the ground state. That is, since the characteristic exponents, such as coherence lengths, define the exponential decay of a small perturbation of a field from its ground state, in order to calculate them one expands fields in the Euler-Lagrange equations around their ground state values (see e.g. Landau and Lifshitz (1980); Plischke and Bergersen (1989)). For a conventional superconductor the coherence length is obtained by expanding in small deviations of the field modulus $|\psi|$ Tinkham (1995), but that cannot a priori be done for our system involving multiple fields. Instead we should expand in small deviations in all degrees of freedom and see if there is a coupling between the fields arising at the lowest order. Because we are dealing with a superconductor we have a coupling to the gauge field $A$ and some care must be taken in handling the gauge invariance of the system. Let us define the phase field

[TABLE]

Note that $\rho_{\alpha}=|\psi_{\alpha}|$ and $\theta_{12}$ are gauge invariant, while $\theta_{\Sigma}$ and $A_{i}$ are not. The combination

[TABLE]

is gauge invariant, and our strategy is to reexpress the Euler-Lagrange equations in terms of $\rho_{\alpha}$, $\theta_{12}$ and $p_{i}$. Let us denote the ground state values of $\rho_{\alpha}$ and $\theta_{12}$ by $u_{\alpha}$ and $\theta_{0}$ respectively; for the $p+ip$ model (3), $\theta_{0}=\pm\frac{\pi}{2}$, but it is instructive to leave it general, for the time being. Then saying that the fields are close to their ground state values means precisely that $p_{i}$, $\varepsilon_{\alpha}$ and $\theta_{\Delta}$ are small, where

[TABLE]

In particular, the small quantities $\varepsilon_{\alpha},p_{i},\theta_{\Delta}$ should obey the linearization of (5) about $(p_{i},\rho_{\alpha},\theta_{12})=(0,u_{\alpha},\theta_{0})$. The left hand side is exactly $\partial_{j}(\partial_{j}p_{i}-\partial_{i}p_{j})$ which is already of linear order, but we must compute the supercurrent $J_{i}$ to linear order. This is straightforward once we recognize that $D_{i}\psi_{\alpha}$ is to linear order,

[TABLE]

so the linearization of (5) is

[TABLE]

Note that the left hand side of this equation is precisely the usual curl of the magnetic field ($p_{i}$ differs from $A_{i}$ by a gradient, so their curls coincide). The key observation is that, unless $\theta_{0}=0$ or $\pi$, that is, unless the ground state is phase locked or antilocked (or $Q^{12}_{ij}\equiv 0$) this PDE couples all the degrees of freedom together (through its final term), so that they all decay to zero with the same dominant length scale. Any other value of $\theta_{0}$ (including $\pm\frac{\pi}{2}$) corresponds to a ground state $(\psi_{1},\psi_{2})=(u_{1}e^{i\frac{1}{2}\theta_{0}},u_{2}e^{-i\frac{1}{2}\theta_{0}})$ which is not gauge equivalent to its complex conjugate, and hence breaks time reversal symmetry. Hence, the effects described below are generic when one has BTRS and spatial anisotropy.

To compute the length scales, we must linearize (5) in $(p_{i},\varepsilon_{\alpha},\theta_{\Delta})$ also. Henceforth, we specialize to the $p+ip$ case with potential (3), so that $\theta_{0}=\pm\frac{\pi}{2}$. Substituting (8),(9) into the Euler-Lagrange equations and discarding all terms nonlinear in small quantities yields

[TABLE]

where we have defined the matrix coefficients,

[TABLE]

and

[TABLE]

is the Hessian of the potential $F_{p}$ about the ground state, with respect to $(\rho_{1},\rho_{2})$.

Note that in the case where there are no mixed gradient terms $Q^{12}=0$ the linearized equations decouple into a pair for $(p_{i},\theta_{\Delta})$ and a pair for $(\varepsilon_{1},\varepsilon_{2})$. That means that small fluctuations in the density fields do not cause a perturbation of the phase difference and do not create magnetic field, as is indeed the case in ordinary superconductors, or in the class of anisotropic models studied in Silaev et al. (2018); Winyard et al. (2019b, a). However we see that for anisotropic superconductors that break time reversal symmetry, such as $p+ip$ superconductors, no such simplification takes place: all the gauge invariant fields $\varepsilon_{\alpha},\theta_{\Delta},p_{i}$ are coupled to one another, and when one changes all the others should change too. The implication of this is that systems like chiral $p+ip$ superconductors cannot be characterized by coherence and magnetic field penetration lengths in the usual sense, but the decay length scale of a small perturbation of the order parameter field and magnetic field is in general the same. Furthermore it implies that one cannot reliably use the London limit to calculate the magnetic field penetration length because density modes do not asymptotically decouple from magnetic modes. Below we calculate these length scales.

Since the equations are anisotropic, to extract the length scales we must first select a direction (normal to the domain wall or material boundary, or radial from the vortex core, depending on context), denoted by a unit vector \mbox{\boldmath{n}}=(n_{1},n_{2}), and then reduce the equations to ordinary differential equations (ODEs) with $n$-dependent coefficients, by imposing translation invariance orthogonal to $n$. So, we demand that

[TABLE]

where $X=n_{i}x_{i}$ and \mbox{\boldmath{n}}^{\perp}=(-n_{2},n_{1}). Substituting (19) into (12)-(15), one obtains a coupled set of five ODEs. The two-vector valued ODE (12) implies a pair of scalar-valued ODEs, obtained by taking its scalar product with $n$ and \mbox{\boldmath{n}}^{\perp}. The $n$ component implies

[TABLE]

(where ${}^{\prime}\equiv d/dX$), which can be used to eliminate $b(X)$ from the other ODEs. We now have 4 coupled ODEs, forming a linear system, that describes the response of the system to a small perturbation about its ground state.

[TABLE]

where $\vec{w}=\left(\varepsilon_{1},\varepsilon_{2},\theta_{\Delta},a\right)^{T}$ and $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$ are certain constant $4\times 4$ real matrices. It is important to note that $\mathcal{A}$ and $\mathcal{C}$ are symmetric, while $\mathcal{B}$ is skew, and that all three depend on the choice of direction $n$. Their exact form is given in Appendix A.

Recall that (21) is the linearized system of field equations describing how a system recovers from a perturbation in the $n$-direction under the assumption of translation invariance orthogonal to $n$, for example, how the system behaves near the boundary of a superconductor subject to an external magnetic field. Its general solution is

[TABLE]

where $\mu_{1},\mu_{2},\ldots,\mu_{8}$, are the solutions of the degree 8 polynomial equation

[TABLE]

The constants $\mu_{i}$ should be interpreted as field masses which set the length scale $\lambda_{i}$ of spatial decay of the associated linear combination of fields via

[TABLE]

The quantities $\vec{v}_{1},\vec{v}_{2},\ldots,\vec{v}_{8}$, are the corresponding eigenvectors (by eigenvector we mean a unit length vector satisfying $(\mu_{i}^{2}\mathcal{A}-\mu_{i}\mathcal{B}+\mathcal{C})\vec{v}_{i}=\vec{0}$), and $c_{1},c_{2},\ldots,c_{8}$ are arbitrary constants, determined by boundary conditions and nonlinearities. Each exponential power is associated to a normal mode, determined by $\vec{v}_{i}$. In an ordinary superconductor the normal mode associated with the coherence length is the modulus of the order parameter, while the magnetic field penetration length is attributed to a massive vector field: the magnetic field. Instead, we see that in the chiral $p+ip$ superconductor the normal modes are associated with linear combinations of magnetic and matter degrees of freedom.

Indeed, the polynomial equation (23) has real coefficients and is quartic in $\mu^{2}$ (since $\mathcal{A},\mathcal{C}$ are symmetric, while $\mathcal{B}$ is skew); hence, if $\mu$ is a solution, so are $-\mu$, $\mu^{*}$ and $-\mu^{*}$. This demonstrates that complex length scales are caused by mixing, as this is the only way for multiple length scales to become linked and hence be complex conjugates of each other. Exactly half the eigenvalues, which we choose to label $\mu_{1},\ldots,\mu_{4}$ have positive real part, while the others have negative real part. We seek solutions that decay to [math] as $X\rightarrow\infty$; these are obtained by setting $c_{i}=0$ for $i\geq 5$ in equation (22).

The long-range behaviour of the fields, in direction $n$, is governed by the dominant eigenvector $\vec{v}_{i}=\vec{v}_{*}$, defined to be the eigenvector whose eigenvalue $\mu_{i}=\mu_{*}$ has smallest positive real part (hence the longest length scale $\lambda_{*}=1/\mu_{*}$ of spatial decay). Note that, in general, $\mu_{*}$ may be complex, in which case the fields at large $X$ are spatially oscillatory, behaving like

[TABLE]

where $c$ is some real constant and $\vec{v}_{r,*}$ is the real part of $\vec{v}_{*}$. The complex magnetic field penetration length implies oscilatory decay of the magnetic field as observed in anisotropic systems without BTRS Silaev et al. (2018); Winyard et al. (2019b, a). Here we find that in a $p+ip$ superconductor one cannot assume that a perturbation of the gap fields will decay with real exponents: i.e. there are no real coherence lengths in general. Note that for dirty isotropic multiband superconductors, the phase difference and density modes can be mixed even without breaking time reversal symmetry Garaud et al. (2018). Our findings of the complete mixing of the order parameters and magnetic modes would apply also for that case.

Importantly, as detailed below, in general one needs to retain contributions from the modes associated with shorter length scales.

Of course, our analysis should reproduce the usual picture of separate real length scales (the coherence length and magnetic penetration depth) in the case of a spatially isotropic system, where $Q^{\alpha\beta}_{ij}=\delta_{\alpha\beta}\delta_{ij}$, and should hold approximately for a small perturbation of this. In the near-isotropic regime, when $Q^{11},Q^{22}\approx I_{2}$ and $Q^{12}\approx 0$, the coupling between $\varepsilon_{\alpha},\theta_{\Delta},a$ is weak, the spectrum is real, and one of the eigenvectors, $\vec{v}_{4}$ say, is approximately $(0,0,0,1)$, while the others, $\vec{v}_{1},\vec{v}_{2},\vec{v}_{3}$, are approximately normal to $(0,0,0,1)$. We then recover the usual picture of separate length scales associated with the magnetic field, $\lambda_{mag}=\lambda_{4}$, and the condensates, $\lambda_{1,2,3}$. Consider the case where $\lambda_{*}=\lambda_{mag}$, that is, $\vec{v}_{4}$ is dominant. Although all fields do, strictly speaking, decay like (25) (with $\mathrm{Im}(\lambda_{*})=0$) at very large $X$, the coefficients in front of $\varepsilon_{1},\varepsilon_{2},\theta_{\Delta}$ are very small, while the coefficient in front of $a$ is of order of unity ($\vec{v}_{*}=\vec{v}_{4}\approx(0,0,0,1)$) so at intermediate range contributions from the subdominant eigenvectors are larger. This allows one to identify approximately $\lambda_{mag}$ as a penetration depth and $\min\{\lambda_{i}\}$ as a coherence length, and classify the system as type-2 (since $\lambda_{mag}$ is the largest length scale). Similar remarks apply if one of the condensate modes is dominant. This is consistent with the numerical solutions obtained earlier in such regimes Garaud et al. (2016). One therefore may approximately call the exponents associated to matter-field-dominated modes coherence lengths and those associated with magnetic-field dominated modes magnetic field penetration lengths. However this approximate picture disappears as one increases the anisotropy and magnetic and matter field couplings in (21) become significant.

In summary, the long range behaviour of spatially decaying solutions of our system is (25) where $\mu_{*}=1/\lambda_{*}$ is the solution of (23) with smallest positive real part. In general, $\mu_{*}$ depends on $n$, the direction along which we impose spatial decay, and may be complex, in which case the decay of both magnetic and gap fields is oscillatory. We have also shown that this coupling and the oscillations in all four fields is a direct result of BTRS and anisotropy.

In the next section we consider the implications of these findings for the Meissner state of a $p+ip$ superconductor.

III Meissner State in a $p+ip$ System

We consider a simple $p+ip$ model, such as the one discussed in the context of the debate of the nature of superconducting state in Sr2RuO4 in Bouhon and Sigrist (2010). This is of the form (1) with (after a trivial rescaling of fields which is shown in detail in appendix D)

[TABLE]

and potential,

[TABLE]

The model contains two unknown parameters: $-1<\nu<1$ which measures the anisotropy of the Fermi surface, and $V_{0}$, the overall strength of the potential (coinciding with $b/(\pi\gamma^{2}K^{2})$ in the notation of Ref. Bouhon and Sigrist (2010)). Its ground states are $(\rho_{\alpha},\theta_{12})=(1,\pm\pi/2)$, so $u_{1}=u_{2}=1$.

Consider a semi-infinite superconductor occupying the half-space $X\geq 0$ (where, as before, $X=n_{1}x_{1}+x_{2}n_{2}$), denoted $\Omega$, with the region $X<0$ occupied by an insulator. Denote by $\partial\Omega$ the boundary between these regions (where $X=0$). Note that \mbox{\boldmath{n}}=(\cos\varphi,\sin\varphi) is an inward pointing unit normal to this boundary. The system is subjected to a uniform external magnetic field $H$ in the $x_{3}$ direction. Provided $H$ is not too strong, the system will approach the ground state $\rho_{1}=\rho_{2}=1$, $\theta_{12}=\pi/2$ (say) in the bulk (as $X\rightarrow\infty$). To find the Meissner state, we minimize the Gibbs free energy

[TABLE]

over all fields in $\Omega$, assuming invariance under translations normal to $n$. Here, we use the standard boundary conditions, advocated in Sigrist and Ueda (1991), by including in the free energy the surface term

[TABLE]

For simplicity, we assume reflection from the boundary is specular, meaning that $\chi_{1}=\chi_{2}=\chi_{3}=\chi>0$. Having imposed translation invariance, the problem reduces to a one-dimensional variational problem on $[0,\infty)$, with natural boundary conditions at [math], which can be solved by a standard gradient-descent method. A more detailed discussion of the boundary conditions is given in appendix B. There is a caveat here. It has been demonstrated recently for $s$-wave superconductors, that boundary conditions can be different in superconductors from those based on the standard assumptions of Caroli-deGennes-Matricon type theory Samoilenka and Babaev (2019), which implies that the standard theory of boundary conditions for $p+ip$ should also be revised. However here we are interested not in the precise field values at the boundary but rather in the laws governing their decay away from the boundary. Therefore the precise form of the boundary conditions is not very important. In Appendix C we present results with the extra boundary terms omitted entirely, giving the same field decay behaviour.

The solutions depend on the unknown model parameters $\nu,V_{0},\chi$ as well as the applied field $H$ and the boundary orientation angle $\varphi$. We have run simulations for $\chi\in\{0,0.01,0.1,1,10\}$, finding no qualitative change in the physics we are focussed on. For that reason we fix $\chi=1$ for the remainder of this section and present a representative sample of the other parameters. We have included a plot, figure 4 for $\chi=0$ in Appendix C for comparison, to demonstrate that the oscillatory behaviour of the fields originates in complex coherence lengths, not from the boundary terms in eq. (28).

For $V_{0}=3$, $\nu=-0.95$, $H=0.3$, the Meissner states with boundary orientations $\varphi=0$ and $\varphi=\pi/3$ are presented in figure 1. Both exhibit oscillatory tails and field inversion of both $B$ and the condensates, consistent with exponential decay with a complex coherence and magnetic field penetration lengths. We shall return to this shortly. If the external field $H$ is increased further, the condensates separate more until, for some $\varphi$, one or other of the densities $\rho_{1}$ or $\rho_{2}$ hits zero. This produces a novel Meissner state depicted in figure 2 for $H=1$ (well below the lower critical field $H_{c1}=1.34$) for several angles $\varphi$. We see that neither matter field component vanishes for $\varphi=\pi/4$, whereas $\rho_{1}$ vanishes for $\varphi=0$, and $\rho_{2}$ vanishes for $\varphi=\pi/2$. As one condensate component goes to zero, the other achieves a maximum exceeding its ground state value, producing a stripe (orthogonal to $n$) of depletion of one condensate and surfeit of the other.

Returning to our main goal of testing the analysis of the previous section, it is straightforward to compute, for any given $\varphi$ (boundary orientation), $\nu$ (anisotropy parameter) and $V_{0}$ (potential energy scale in the GL energy) the dominant eigenvalue $\mu_{*}$, and hence map out the parameter set on which $\mu_{*}$ is complex. Figure 3 presents pictures of the $(\varphi,\nu)$ parameter plane, for a sequence of values of $V_{0}$, coloured to show the parameter domain where $\mu_{*}$ is complex. For $V_{0}$ small, the parameter domain of complex $\mu_{*}$ is small and confined to the edges where $|\nu|$ is close to $1$, but as $V_{0}$ increases, the domain swells, eventually covering the whole parameter space (when $V_{0}\approx 4$), predicting that the Meissner state should be spatially oscillatory for all anisotropies $\nu$ and all boundary orientations $\varphi$ if $V_{0}$ is around this value. Increasing $V_{0}$ still further, pockets of real $\mu_{*}$ return and gradually refill the whole parameter space for very high values of $V_{0}$. Turning to the parameter sets of figure 1, $V_{0}=3$, $\nu=-0.95$ and $\varphi=0,\pi/3$, we find in both cases that $\mu_{*}$ is complex, consistent with the nonlinear numerics ($\mu_{*}(0)=0.689+0.548i$, $\mu_{*}(\pi/3)=1.013+0.648i$). The oscillatory decay predicted by linear analysis predicts that the zeros of $B$, and of $\rho_{\alpha}-u_{\alpha}$ should be equally spaced with period $\pi/\mathrm{Im}(\mu_{*})$, as should successive extrema of these functions. These gap widths can easily be extracted from the nonlinear numerics, and are displayed, for these parameter sets, alongside the linear prediction, in figure 1. The agreement is remarkable. Finally, we have also chosen parameter sets for which $\mu_{*}$ is real, so that the linearization predicts non-oscillatory decay. While the solutions still generically exhibit a single peak in each field, after this initial overshoot in the nonlinear regime, the fields decay exponentially without oscillation as the linearization predicts.

IV Conclusions

In conclusion we have shown that the normal modes in an anisotropic superconductor that breaks time reversal symmetry mix density and phase fields with the magnetic field. This precludes using the usual notion of coherence and magnetic field penetration lengths because long range decay of matter and magnetic fields is given by the same exponent. Additionally the fundamental length scales, associated with the normal modes, that mix the order parameters and magnetic field, are in general complex. We have also shown that this mode mixing requires BTRS along with anisotropy and that mixing is required for complex length scales. While systems exist with oscillations in magnetic field and phase difference due to anisotropy driven mixing between these two modes Silaev et al. (2018), all four fields having complex length scales can only happen in an anisotropic BTRS system. Calculating numerically the Meissner effect in a chiral $p+ip$ superconductor, we indeed find that application of an external magnetic field is screened in an oscillatory way and produces damped oscillatory decay of the order parameter fields. For strong anisotropy the effect should be detectable in muon spin relaxation experiments. Cutting sample boundaries under different angles relative to crystal axes and measuring magnetic field inversion can allow one to recover information about the order parameter. Finally we note that our analysis dictates that the effect is present for any inhomogeneous situation, including the domain wall excitations in $p+ip$ superconductors considered in Bouhon and Sigrist (2010); these should also exhibit oscillation and field inversion. That this was not observed in Bouhon and Sigrist (2010) might be an artifact of an overly restrictive ansatz. We plan to examine this further in a separate publication.

Appendix A A: The coupling matrices

Here we record the non-zero matrix elements of the $4\times 4$ matrices appearing in equation 21:

[TABLE]

Recall that $\mathcal{A}_{ij}\equiv\mathcal{A}_{ji}$, $\mathcal{B}_{ij}\equiv-\mathcal{B}_{ji}$ and $\mathcal{C}_{ij}=\mathcal{C}_{ji}$.

Appendix B B: Boundary Conditions

To compute the Meissner state in the region $\Omega$ numerically we must minimize the Gibbs free energy

[TABLE]

among all fields defined on $\Omega$. It is convenient to include a gauge-fixing term $\frac{1}{2}(\partial_{i}A_{i})^{2}$ in ${\cal F}$, and to denote the dynamical fields collectively as $\phi_{a}$, $a=1,\ldots,6$ (consisting of the real and imaginary parts of $\psi_{\alpha}$, and $A_{1}$, $A_{2}$). Then, under a variation $\delta\phi_{a}$, $G$ varies as

[TABLE]

where we have used the divergence theorem, and recalled that $n$ is an inward pointing normal to $\partial\Omega$. Demanding that $\delta G=0$ for all variations requires both these integrals vanish identically, and hence that $\phi_{a}$ satisfy the usual Euler-Lagrange equations in $\Omega$ together with the boundary conditions

[TABLE]

on $\partial\Omega$. For the model studied here, this reduces to

[TABLE]

Imposing the translationally invariant ansatz $\psi_{\alpha}=\psi_{\alpha}(X)$, $A_{i}=a(X)n_{i}^{\perp}+b(X)n_{i}$, where $X=n_{i}x_{i}$, this reduces further to

[TABLE]

These are the boundary conditions we impose at $X=0$. At $X=L$, large (our effective infinity), we demand that $b^{\prime}=a^{\prime}=0$, $\psi_{1}=u_{1}$ and $\psi_{2}=iu_{2}$ (the fields are in their ground state state).

Appendix C C: $\chi=0$ results

To confirm that the long-range decay behaviour holds for various boundary conditions, we include here a plot in figure 4 for the parameters used in our results section but with the boundary term removed, $\chi=0$.

Appendix D D: rescaling of fields

We have made use of the form of the potential argued for in Bouhon and Sigrist (2010), however we have made a few rescalings to rewrite the proposed model in a simpler fashion. The proposed model is,

[TABLE]

where $D_{i}=\partial_{i}-i\gamma A^{b}_{i}$ and some of the parameters are coupled such that,

[TABLE]

We will write our condensate fields as,

[TABLE]

and rescale our gauge field,

[TABLE]

Finally we can use a spatial rescaling,

[TABLE]

and then rescale the total energy to be,

[TABLE]

This finally gives the form of the energy given in equation 1 with potential,

[TABLE]

and anisotropy tensors,

[TABLE]

Where we have collected multiple parameters together,

[TABLE]

Note that without loss of generality we have reduced the number of parameters to two $(V_{0},\nu)$. This leads to the vacua and hence asymptotic values being $\theta_{12}=\pm\pi/2$ as required for BTRS and $\rho_{1}=\rho_{2}=1$ without loss of generality.

Acknowledgements

We thank Mihail Silaev for collaboration and discussions. The work of MS and TW is supported by the UK Engineering and Physical Sciences Research Council through grant EP/P024688/1. EB is supported by the Swedish Research Council Grants No. 642-2013-7837, 2016-06122, 2018-03659 and Göran Gustafsson Foundation for Research in Natural Sciences and Medicine. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.