Signatures of Nematic Superconductivity in Doped Bi$_2$Se$_3$ under Applied Stress

TL;DR

This paper investigates how applied stress affects nematic superconductivity in doped Bi$_2$Se$_3$, proposing experimental methods to verify the nematic nature of the order parameter and its impact on Josephson tunneling.

Contribution

It develops a Ginzburg-Landau theoretical framework for nematic superconductors under stress and suggests experimental probes to confirm nematic order.

Findings

Stress can be used to probe the nematic superconducting phase.

Applied stress influences Josephson tunneling, acting as an on/off switch.

Temperature dependence of critical current varies with stress, differing from conventional s-wave behavior.

Abstract

The $M_x \text{Bi}_2 \text{Se}_3$ family are candidates for topological superconductors, where $M$ could be Cu, Sr, or Nb. Two-fold anisotropy has been observed in various experiments, prompting the interpretation that the superconducting state is nematic. However, it has since been recognized in the literature that a two-fold anisotropy in the upper critical field $H_{c2}$ is incompatible with the na\"{i}ve nematic hypothesis. In this paper we study the Ginzburg-Landau theory of a nematic order parameter coupled with an applied stress, and classify possible phase diagrams. Assuming that the $H_{c2}$ puzzle is explained by a pre-existing "pinning field", we indicate how a stress can be applied to probe an extended region of the phase diagram, and verify if the superconducting order parameter is indeed nematic. We also explore the Josephson tunneling between the proposed nematic superconducting state and an s-wave superconductor. The externally applied stress is predicted to serve as an on/off switch to the tunneling current, and in certain regime the temperature dependence of the critical current can be markedly different from that between two conventional s-wave superconductors.

Equations

$$\vec{\eta }=\left(\begin{array}{c}{\eta }_x\\ {\eta }_y\end{array}\right)=\sqrt{D}\left(\begin{array}{c}\cos \theta \\ e^{i\varphi }\sin \theta \end{array}\right).$$

$$\begin{array}{rl}{\mathcal{F}}_0=& \alpha |\eta {|}^2+\left[{\beta }_1^{\prime }|\eta {|}^4+{\beta }_2^{\prime }\left(\vec{\eta }\cdot \vec{\eta }\right)\left({\vec{\eta }}^{\ast }\cdot {\vec{\eta }}^{\ast }\right)\right]\\ & +\left\{{\gamma }_1^{\prime }|\eta {|}^6+{\gamma }_2^{\prime }\left[({\eta }_1+i{\eta }_y{)}^3({\eta }_x^{\ast }+i{\eta }_y^{\ast }{)}^3+({\eta }_x-i{\eta }_y{)}^3({\eta }_x^{\ast }-i{\eta }_y^{\ast }{)}^3\right]+{\gamma }_3^{\prime }|\eta {|}^2\left(\vec{\eta }\cdot \vec{\eta }\right)\left({\vec{\eta }}^{\ast }\cdot {\vec{\eta }}^{\ast }\right)\right\}.\end{array}$$

$$\alpha =\kappa \left(\frac{T-T_0}{T_0}\right),$$

$${\mathcal{F}}_0=\alpha D+\beta D^2+\left[{\gamma }_1+{\gamma }_2\cos \left(6\theta \right)\right]D^3,$$

$$\begin{array}{rl}\vec{\epsilon }& =\left(\begin{array}{c}{\epsilon }_{xx}-{\epsilon }_{yy}\\ -2{\epsilon }_{xy}\end{array}\right)\equiv \left(\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\end{array}\right)\\ & =\epsilon \left(\begin{array}{c}\cos \left[2(-\Phi +{\Phi }_0)\right]\\ \sin \left[2(-\Phi +{\Phi }_0)\right]\end{array}\right),\end{array}$$

$$\vec{S}=$$

$$\vec{T}=$$

$$\vec{U}=$$

$${\mathcal{F}}_{\epsilon }=-g_0^{\prime }\vec{\epsilon }\cdot \vec{S}-g_1^{\prime }\vec{\epsilon }\cdot \vec{T}+g_2^{\prime }\vec{\epsilon }\cdot \vec{U}.$$

$$-g_0^{\prime }\vec{\epsilon }\cdot \vec{S}=-g_0\epsilon D\cos \left(2\theta -2\Phi \right)$$

$$-g_1^{\prime }\vec{\epsilon }\cdot \vec{T}=-g_1\epsilon D^2\cos \left(2\theta -2\Phi \right).$$

$$g_2^{\prime }\vec{\epsilon }\cdot \vec{U}=-g_2\epsilon D^2\cos \left(4\theta +2\Phi \right),$$

$$\begin{array}{rl}{\mathcal{F}}_{\epsilon }=& -g_0\cos \left(2\theta -2\Phi \right)\epsilon D\\ & -\left[g_1\cos \left(2\theta -2\Phi \right)+g_2\cos \left(4\theta +2\Phi \right)\right]\epsilon D^2.\end{array}$$

$$\mathcal{F}={\mathcal{F}}_0+{\mathcal{F}}_{\epsilon }=aD+b{D}^2+c{D}^3,$$

$$b=\beta -g_1\epsilon \cos (2\theta -2\Phi )-g_2\epsilon \cos (4\theta +2\Phi ).$$

$$\epsilon \ll \frac{\beta }{|g_1|},\frac{\beta }{|g_2|}.$$

$$a=\alpha -g_0\epsilon \cos (2\theta -2\Phi ).$$

$$\alpha ={\alpha }_1(\epsilon )\equiv g_0\epsilon ,$$

$$g_1<g_2.$$

$$T_1(\epsilon )=\kappa {\alpha }_1(\epsilon )+T_0$$

$$\Delta c_{v,1}=\frac{{\kappa }^2}{2T_1(\epsilon )}\frac{1}{\beta -\epsilon [g_1+g_2\cos (6\Phi )]}.$$

$${\left(\frac{\partial \mathcal{F}}{\partial D}\right)}_{D=\tilde{D}}=0.$$

$$\begin{array}{rl}\tilde{\mathcal{F}}& \equiv a\tilde{D}+b{\tilde{D}}^2+c{\tilde{D}}^3\\ & =k_0+k_2{\theta }^2+k_4{\theta }^4+\dots \end{array}$$

$$T_2=T_0\left[1-\frac{2}{3}\frac{\beta }{\kappa }\sqrt{\frac{g_0}{{\gamma }_2}}\sqrt{\epsilon }+\frac{g_0}{\kappa }\left(\frac{2{\gamma }_2-{\gamma }_1}{3{\gamma }_2}\right)\epsilon \right],$$

$$\Delta c_{v,2}=\frac{{\kappa }^2{T}_2(\epsilon )}{T_0^2}\frac{9{\gamma }_2{g}_0\epsilon }{2\left[g_0\epsilon (4{\gamma }_1-5{\gamma }_2+4\beta \sqrt{{\gamma }_2{g}_0\epsilon })\right]\left[\beta +({\gamma }_1+{\gamma }_2)\sqrt{g_0\epsilon /{\gamma }_2}\right]}.$$

$$\underset{\epsilon \to 0}{\lim }\left(\Delta c_{v,1}+\Delta c_{v,2}\right)=\Delta c_v,$$

$$\begin{array}{rl}\tilde{\mathcal{F}}\approx & |\alpha {|}^{3/2}\left[3{\gamma }_1+3{\gamma }_2\cos (6\theta )\right]\\ & +\frac{|\alpha |\left\{\beta -\epsilon \left[g_1\cos (2\theta -2\Phi )+g_2\cos (4\theta +2\Phi )\right]\right\}}{3\left[{\gamma }_1+{\gamma }_2\cos (6\theta )\right]}\\ & +O\left(|\alpha {|}^{1/2}\right).\end{array}$$

$${\epsilon }_c=\frac{36g_0{\gamma }_3}{(g_1-4g_2{)}^2}.$$

$${\left(\frac{{\partial }^2}{\partial {\sigma }^2}\tilde{\mathcal{F}}\left(\alpha ={\alpha }_2,\epsilon ,\sigma ,\Phi =0\right)\right)}_{\sigma =0}=0.$$

$$\begin{array}{rl}\tilde{\mathcal{F}}({\alpha }_2({\epsilon }_{\ast })& ,{\epsilon }_{\ast },\sigma =0,\Phi =0)\\ & =\tilde{\mathcal{F}}({\alpha }_2({\epsilon }_{\ast }),{\epsilon }_{\ast },\sigma =\pi /2,\Phi =0)\end{array}$$

$$\vec{p}=p\left(\begin{array}{c}\cos (-2{\Phi }_p)\\ \sin (-2{\Phi }_p)\end{array}\right)\equiv \left(\begin{array}{c}{p}_1\\ p_2\end{array}\right)$$

Full text

I Introduction

$\text{Bi}_{2}\text{Se}_{3}$ is a topological insulator Xia et al. (2009); Zhang et al. (2009); Liu et al. (2010). When intercalated with electron-donating atoms such as CuHor et al. (2010); Kriener et al. (2011), SrShruti et al. (2015) or NbQiu et al. (2015); Asaba et al. (2017), it becomes superconducting at low temperature. Following their identification as candidates of topological superconductors Fu and Berg (2010); Fu (2014); Venderbos et al. (2016a), the family of materials has recently garnered a lot of interestYonezawa (2018). Experimentally, two-fold anisotropy in the superconducting phase has been observed in NMR Knight’s shiftMatano et al. (2016), specific heatYonezawa et al. (2017); Willa et al. (2018); Sun et al. (2019), magnetic torqueAsaba et al. (2017), magneto-transport and upper critical fieldPan et al. (2016); Yonezawa et al. (2017); Du et al. (2017); Smylie et al. (2018). This is incompatible with the crystal lattice symmetry, and has led to the hypothesis that the superconducting state is nematic in nature, i.e. it spontaneously breaks both lattice rotational and gauge symmetries. While the superconductivity is suppressed by an applied hydrostatic pressureBay et al. (2012); Nikitin et al. (2016), the two-fold anisotropy is observed as far as superconductivity holdsNikitin et al. (2016).

The parent material $\text{Bi}_{2}\text{Se}_{3}$ forms a rhombohedral crystal, with a lattice point group $D_{3d}$; see FIG 1. According to the nematic hypothesis, the complex superconducting order parameter $\vec{\eta}$ is believed to transform under the irreducible representation $E_{u}$Fu (2014), which is parity-odd and two-dimensional, and, in the absence of other symmetry breaking effect, spontaneously breaks the crystal rotation symmetry together with the $U(1)$ symmetry. It can be thought of as a nematic director in the basal plane, since $\vec{\eta}$ and $-\vec{\eta}$ can be identified up to a global $U(1)$ phase shift of $\pi$. Experiments that detect the anisotropy in the basal plane all reported a two-fold symmetry, consistent with the presence of a nematic superconducting order parameter. This naïve picture, however, is contradicted by several pieces of experimental facts.

First of all, the upper critical field $H_{c2}$ is experimentally observed to be two-fold anisotropicPan et al. (2016); Yonezawa et al. (2017); Du et al. (2017); Smylie et al. (2018); Sun et al. (2019) with the applied field parallel to the basal plane, seemingly inline with other measurements. However, this result is actually incompatible with the naïve nematic hypothesis, which indicates a six-fold anisotropy for $H_{c2}$Agterberg and Walker (1995a); Krotkov and Mineev (2002); Venderbos et al. (2016b).

To understand this, one considers the onset of superconductivity in high field. Instead of a finite $\vec{\eta}$ providing a preferred orientation, the magnetic field and the lattice anisotropy together determine the orientation of the infinitesimal $\vec{\eta}$ when it first emerges. This consideration leads to the prediction of a six-fold anisotropy associated with the six-fold improper rotation symmetry of $D_{3d}$.

The robustness of the two-fold nematic direction observed in experiments is another puzzle unexplained by the naïve nematic scenario. It was noted that the two-fold anisotropy for a given sample is always pinned at the same direction. In Pan *et. al.*Pan et al. (2016) part of the experiment was performed with the sample repeatedly cycled up to $5$K (about $1.6T_{c}$) and cooled back down to superconducting phase again before every field sweep111A. de Visser, private communication. In a separate account222M. P. Smylie, private communication, the same two-fold direction for any given sample persists as it was cycled in and out of the superconducting phase repeatedly, sometimes all the way to room temperature and back, with or without a background magnetic field. To account for this behavior, the rotational symmetry must have been already broken in some way at the room temperature.

Indeed, for the case of Sr, evidence of normal-state nematic response was reported by Kuntsevich *et. al.*Kuntsevich et al. (2018) and Sun *et. al.*Sun et al. (2019). On the other hand, the Argonne groupWilla et al. (2018); Smylie et al. (2018) reported isotropic normal-state response for Sr. This is in contrast with the case of Cu, where no nematicity was ever observed above $T_{c}$Matano et al. (2016); Yonezawa et al. (2017) to our knowledge.

To reconcile the $E_{u}$ pairing scenario and the observed robust pinning effect, some nematic symmetry-breaking field (SBF) must already be present in the sampleWilla et al. (2018). The precise nature of this SBF is yet unknown, though mechanical stress or strain were named as possible candidates Fu (2014). This SBF pins all the two-fold nematic responses (including $H_{c2}$) in the same direction for a given sample.

We emphasize that any scenario that solely relies on spontaneous symmetry breaking at low temperature to supply the nematic direction is incompatible with experiments. The naïve nematic scenario does not work, nor does any scenario where the rotational symmetry is broken independent of the $U(1)$ symmetry, for exampleHecker and Schmalian (2018). A pre-existing SBF (“pining field”) is necessary.

A conceptual difficulty then arises: in the presence of a pinning field, an otherwise isotropic s-wave superconductor can in principle also appear two-fold anisotropic. Can we tell this so-called anisotropic s-wave scenario apart from the true nematic scenario? This distinction is more than merely academic, because the nematic superconductor is topologicalFu (2014); Venderbos et al. (2016a), while the s-wave case is not.

In light of these recent developments, a better understanding of how the proposed nematic superconductivity interacts with an explicit SBF is clearly desirable. In this work we consider mechanical stress or strain as the prototypical SBF, although an SBF with a different microscopic nature will have exactly the same phenomenology. We further consider the application of stress in addition to the existing, unknown SBF, as an experimental probe of the nature of the superconducting phase.

Historically, much of the phenomenology of a superconducting order parameter with non-trivial symmetry has been known as part of the lore of the unconventional superconductors (see Sigrist and Ueda (1991); Joynt and Taillefer (2002) and references therein.) Pertinent to the present discussion is the fact that, due to the competition of the SBF and the crystal field, the superconducting transition may be split into two if the order parameter is multi-componentHess et al. (1989); Machida and Ozaki (1989). In the heavy fermion compound $\text{UPt}_{3}$, the double transition is due to the antiferromagnetic order in the normal state, and has been experimentally observedFisher et al. (1989). The idea of a stress causing similar split has been put forward for $\text{Sr}_{2}\text{Ru}\text{O}_{4}$Sigrist et al. (1987); Walker and Contreras (2002), but to our knowledge has never been observedHicks et al. (2014); Watson et al. (2018).

In this paper, we formulate the Ginzburg-Landau (GL) theory of a nematic superconducting order parameter coupled to an external SBF, and explore the phenomenology as a guide to future experiments. In particularly, following the suggestion of Fu (2014), we consider stress as the prototypical SBF. The GL theory is written down in Section II, and then in Sections III, IV and V, we identify the splitting of the superconducting transition, as well as another possible transition at a lower temperature, and classify possible phase diagrams as function of stress and temperature. In section VI, we further explore the phenomenology when the system is subjected to two SBFS: a pre-existing pinning field and a stress applied as an experimental probe, and propose experimental signatures that may help distinguish a truly nematic superconducting order parameter from a single-component one.

In addition, in Section VII we consider the critical tunnel current in the Josephson junction between the proposed nematic superconducting state and another s-wave superconductor. The idea of using the Josephson effect as a probe for unconventional pairing symmetry has long been consideredTanaka and Kashiwaya (1997); Sumiyama et al. (1998); Tanaka and Kashiwaya (1999); Stefanakis (2002); Gouchi et al. (2012). The Ginzburg-Landau approach is also applicable to describe the Josephson tunnelingYip et al. (1990). If the nematic hypothesis does hold, the critical current shows a non-trivial anisotropic dependence on the SBF, as well as a temperature dependence that differs from the s-wave-to-s-wave scenario. This provides another experimental test for the nematicity.

II Ginzburg-Landau Free Energy

We will formulate the GL theory by treating the system as “almost $D_{3d}$-invariant”. That is, we consider the nematic superconducting order parameter living in the $E_{u}$ representation of $D_{3d}$, and write down a $D_{3d}$-invariant free energy. The pinning effect is introduced via coupling to a (possibly weak) explicit SBF.

This scenario is similar to the nematic pseudogap phase of YBCOAndo et al. (2002); Hinkov et al. (2008); Daou et al. (2010); Sato et al. (2017), where the order parameter would have broken the four-fold rotation symmetry of the CuO2 plane, except that this “symmetry” is already weakly broken by the orthorhombic lattice. Nonetheless the nematic response is greatly enhanced by the nematic pseudogap order parameter.

II.1 Free Energy without SBF

Let us first give a brief recap of the GL free energy for a nematic superconductor in the absence of any SBF.

We define the coordinate axes as follows: the $z$-axis is aligned with the principal $C_{3}$ axis of the lattice, the $x$-axis is along one of the $C_{2}$ axes in the basal plane, and the $y$-axis is chosen to form a right-handed set of coordinate system.

We assume that the superconducting order parameter is a complex two-component quantity $\vec{\eta}$. It will be parameterized as:

[TABLE]

Here $\eta_{x}$ is kept real by a suitable global $U(1)$ phase rotation. See Fig 1 for illustration.

The point group $D_{3d}$ can be built from three operations: $2\pi/3$-rotation about the $z$-axis (denoted $C_{3}$), $\pi$-rotation about the $x$-axis (denoted ${C_{2}}^{\prime}$), and mirror reflection about the $yz$-plane (denoted $M$). The odd-parity proposalFu (2014) would have $\vec{\eta}$ transforming in the $E_{u}$ representation, similar to $(x,y)$.

However, we may equally well consider the possibility where $\vec{\eta}$ transforms in the $E_{g}$ representation instead, similar to the pair $(yz,-xz)$. Gauge invariance requires $\vec{\eta}$ and $\vec{\eta}^{*}$ to always come in pairs in the GL expansion, and the sign of parity is irrelevant. The GL free energy and most of the subsequent analysis remain unchanged for the $E_{g}$ scenario, except that a minor modification is needed for the Josephson tunneling, which does not affect our qualitative conclusion. For either case, one can consistently define $\eta_{x}$ to be the component invariant under ${C_{2}}^{\prime}$.

The quantity $\vec{\eta}$ is to be considered as a headless vector, since a sign in $\pm\vec{\eta}$ can be regarded as a $U(1)$ phase factor. Therefore one may restrict $\theta$ to the interval $(-\pi/2,\pi/2]$, keeping in mind that $\theta$ and $\theta+\pi$ are always degenerate.

The GL expansion of free energy is obtained by writing down all the terms that are invariant under $D_{3d}$, global phase shift of $\vec{\eta}$, and time reversal. Up to $O(\eta^{6})$, the most generic form is:

[TABLE]

Here $\alpha$ is taken to be

[TABLE]

where $T_{0}$ is the superconducting transition temperature without external field. We ignore the temperature dependence of all other coefficients.

Let us first consider the free energy up to $O(\eta^{4})$. To this order $\mathcal{F}_{0}$ is in fact invariant under arbitrary rotations in the $xy$-plane. Two solutions for $\vec{\eta}$ are possible: the complex chiral state that is rotationally invariant but spontaneously breaks time-reversal symmetry, and the real nematic state which breaks the rotational symmetry but is time-reversal invariant. The sign of $\beta_{2}^{\prime}$ decides which phase is favoredVolovik (1988); Agterberg and Walker (1995b).

Since the nematic state is by assumption the true solution, we take the appropriate sign $\beta_{2}^{\prime}<0$. The minimum of $\mathcal{F}$ now lies along $\phi=0$, and in effect $\left(\vec{\eta}\cdot\vec{\eta}\right)\left(\vec{\eta}^{\;*}\cdot\vec{\eta}^{\;*}\right)=|\eta|^{4}$. Therefore we define $\beta\equiv\beta_{1}^{\prime}+\beta_{2}^{\prime}$, and require $\beta>0$ for stability.

Now let us move on to include terms of $O(\eta^{6})$. The $\gamma_{3}^{\prime}$ term essentially has the same $\theta$- and $\phi$- dependence as the $\beta_{2}^{\prime}$ term. Since a transition from nematic to chiral state at a lower temperature is never observed, we explicitly require $\gamma_{3}^{\prime}<0$. On the other hand, it can be shown that all extrema of the $\gamma_{2}^{\prime}$ term correspond to real values of $\vec{\eta}$. Consequently, one may assume $\vec{\eta}$ is real and simply set $\phi=0$. In terms of $\theta$ and $D$, the free energy reads:

[TABLE]

where $\gamma_{1}\equiv\gamma_{1}^{\prime}+\gamma_{3}^{\prime}$ and $\gamma_{2}=2\gamma_{2}^{\prime}$. For stability we require $\gamma_{1}>0$ and $\gamma_{1}>|\gamma_{2}|$.

The $\gamma_{2}$ term breaks the full rotational symmetry down to a discrete six-fold symmetry. This is the lowest order at which the crystal anisotropy enter the GL expansionAgterberg and Walker (1995b).

II.2 Stress as an SBF

Stress serves as our prototypical SBF. We will focus on stress in the $xy$-plane. Generally, planar stress is a rank-two symmetric tensor $\hat{\varepsilon}$ with three independent real parameters. Under $D_{3d}$, this tensor decomposes into two parts: the scalar $\operatorname{Tr}{\hat{\varepsilon}}$, and the traceless part which is organized to form the two-component quantity

[TABLE]

which transforms in the $E_{g}$ representation under $D_{3d}$.

Gauge invariance of the free energy requires that $\vec{\eta}$ and $\vec{\eta}^{\,*}$ to pair up, and all allowed combinations are parity-even. Therefore to couple to the nematic superconductor at all, an SBF must transform under the even-parity $E_{g}$ representation. It therefore shares the same phenomenology with stress (5). One can simply replace the stress with any pre-existing SBF of an arbitrary origin, and the resulting GL free energy is formally identical. In particular, we note that for externally applied strain, the form of GL free energy (14) remains exactly identical. The result in this paper equally applies to experiments which uses strain instead of stress.

We choose to parameterize $\vec{\varepsilon}$ as (5) because $\Phi$ would correspond to the physical rotation angle about the $z$-axis. When the stress is rotated by an arbitrary angle $\delta$, $\vec{\varepsilon}$ changes via $\Phi\rightarrow\Phi+\delta$. The angle $\Phi_{0}$ defines the direction relative to the $x$-axis that corresponds to $\Phi=0$, and we leave it unspecified for now, to be chosen for convenience later. The angle $\Phi$ defines a special direction in the basal plane, and subsequently we may refer to it as the “orientation” of the stress. This breaks the point group rotation symmetry.

For the sake of completeness, we state that one can form another $E_{g}$ pair using out-of-plane stress components: the pair $(\varepsilon_{yz},-\varepsilon_{xz})$ can be used in place of (5). Note that this is possible only for the trigonal $D_{3}d$ group. If B and C sites in Fig 1 were to be identified, the resultant hexagonal $D_{6h}$ point group wouldn’t have allowed the replacement.

For the remainder of this paper, we will always assume that $\operatorname{Tr}{\hat{\varepsilon}}$ is kept constant. The dependence of GL coefficients on $\operatorname{Tr}{\hat{\varepsilon}}$ can therefore be entirely disregarded.

We will only consider coupling to $\vec{\varepsilon}$ at linear order. As discussed above, gauge invariance requires $\vec{\eta}$ and $\vec{\eta}*$ to comes in pairs that transform under the $E_{g}$ representation. Up to order $O(\eta^{4})$, the following is the exhaustive list of possible combinations:

[TABLE]

The additional terms entering the free energy are

[TABLE]

The leading $\vec{\varepsilon}\cdot\vec{S}$ term was considered in the literatureVenderbos et al. (2016b); Willa et al. (2018). We go beyond quadratic order in $\vec{\eta}$ to study the interplay between the applied stress and crystal anisotropy.

All three terms in (9) can be shown to be extremized by real values of $\vec{\eta}$. Together with (2), every term in the free energy favors a real $\vec{\eta}$. One is again allowed to set $\phi=0$ from this point on.

Let us first look at $\vec{\varepsilon}\cdot\vec{S}$. Define $g_{0}=|g_{0}^{\prime}|>0$. By appropriately choosing $\Phi_{0}=0$ or $\pi/2$, this term can be written as

[TABLE]

The free parameter $\Phi_{0}$ has been used up to fix the sign of $g_{0}$; now the $\vec{T}$ term must be taken as-is. Define $g_{1}=\pm g_{1}^{\prime}$ with appropriate sign depending on the previous choice of $\Phi_{0}$, and one may write:

[TABLE]

Similarly, the $\vec{U}$ term becomes:

[TABLE]

with $g_{2}=\pm g_{2}^{\prime}$ appropriately defined to absorb the possible sign due to the choice of $\Phi_{0}$.

To sum up, the stress dependent part of GL free energy is

[TABLE]

The coefficient $g_{0}$ is made positive by appropriate choice of $\Phi_{0}$. The other two coefficients $g_{1}$ and $g_{2}$ can be of either signs.

It can be seen that (4), (10), (11) and (12) are all invariant under six-fold (improper) rotation about the $z$-axis. Consequently, we are free to impose the restriction $\Phi\in(-\pi/6,\pi/6]$, which amounts to redefining the $x$-direction relative to the external stress.

The full expression of GL free energy $\mathcal{F}$ employed in this paper is the sum of (4) and (13). Let us also introduce another notation:

[TABLE]

where $a$, $b$ and $c$ are functions of $\alpha$, $\theta$, $\Phi$ and $\varepsilon$.

II.3 Limit on the Magnitude of Stress

The explicit form of coefficient $b$ in $\mathcal{F}$ is

[TABLE]

If we take this expression at its face value, for $\varepsilon$ large enough, the minimum value of $b$ turns negative, and the transition from normal to superconducting state becomes first order.

While not implausible, this stress-induced first order transition has not been observed in any material to our knowledge. We therefore limit the range of $\varepsilon$ in our theoretical investigation to avoid this regime. The appropriate condition is:

[TABLE]

II.4 Sign of $\gamma_{2}$

In the absence of an SBF, up to order $D^{2}$ the free energy $\mathcal{F}_{0}$ is symmetric under arbitrary rotation around the z-axis. That is, $\eta_{x}$ and $\eta_{y}$ are completely degenerate. The six-fold symmetric $\gamma_{2}$ term breaks the degeneracy between $\eta_{x}$ and $\eta_{y}$. When $\gamma_{2}<0$, the ground state has non-zero $\eta_{x}$, while $\gamma_{2}>0$ results in non-zero $\eta_{y}$Agterberg and Walker (1995b).

Using a two-band lattice model, FuFu (2014) argued that $\eta_{y}$ is the correct superconducting order parameter for $\text{Cu}_{x}\text{Bi}_{2}\text{Se}_{3}$, which corresponds to $\gamma_{2}>0$. We will thus assume the positive sign for the remainder of the paper. However, we point out that all the result obtained subsequently can be easily mapped to the case where $\gamma_{2}<0$, should it turn out that way in experiment.

Let $\Delta\theta=\theta-\Phi$, and eliminate all occurrences of $\theta$ in favor of $\Delta\theta$ in (4) and (13). By shifting $\Phi\rightarrow\Phi+\pi/6$, one effectively reverses the signs of both $\gamma_{2}$ and $g_{2}$.

III Phase Diagrams

In this section we will describe the three possible phase diagrams, assuming the GL free energy (14). The reduced temperature $\alpha$, the magnitude of the stress $\varepsilon$, and the orientation of the stress $\Phi$ will span the three axes of the phase diagrams. The derivation and analysis of the features on these phase diagrams will be given in later sections.

We want to emphasize that any of the non-trivial structure of the phase diagram is a direct consequence of the nematic superconducting order parameter $\vec{\eta}$. An s-wave superconductor coupled to an explicit SBF may show anisotropic response, but will only exhibit a single superconducting phase, with $T_{c}$ depending on $\vec{\varepsilon}$ analytically.

First let us consider the stress-free case $\varepsilon=0$. The superconducting transition takes place at $\alpha=0$, and the order parameter is six-fold rotationally degenerate: $\theta=-\pi/6$, $\pi/6$, or $\pi/2$. We will refers to these directions as the “natural minima” of the free energy.

Let’s turn on the stress. By assumption, the size of stress $\varepsilon$ is such that $b$ remain positive for all values of $\theta$ and $\Phi$. The normal-to-superconducting transition is then solely controlled by the coefficient $a$ in (14):

[TABLE]

At finite stress, the transition occurs at

[TABLE]

and the order parameter is directed at $\theta=\Phi$. We will refer to this as the upper transition.

III.1 $\Phi$ = 0

We first examine the case when the direction of stress is fixed at $\Phi=0$. This is a plane in the full three-dimensional parameter space $(\Phi,\varepsilon,\alpha)$.

We first orient the stress along $\Phi=0$. Below $T_{c}$, the order parameter is locked to $\theta=0$ by symmetry. This is labeled as phase A. As discussed earlier, the $\gamma_{2}$ term favors $\theta$ along one of the natural minima. The competition between $g_{0}$ and $\gamma_{2}$ eventually leads to a second order phase transition from $\theta=0$ to $\theta\neq 0$ at a lower temperature $\alpha_{2}(\varepsilon)$. This will be referred to as the middle transition. The lower temperature phase will be called phase B. See Fig 2.

This can be explained in more physical terms. A hexagon has two types of high-symmetry directions. When the stress favors the $x$-type direction, but the naturally preferred orientation is of the $y$-type (or vice versa), the middle transition exists as a result of the competition.

This middle transition only exists for $\Phi=0$. Otherwise, $\theta$ will simply be pulled toward the nearer of $\pm\pi/6$ due to the reduced symmetry of the setup. At $\varepsilon=0$, the upper and middle transition merge into a single superconducting transition.

Below $\alpha_{2}(\varepsilon)$, the $\Phi=0$ plane is a first order coexistent surface, separating the $\theta>0$ and $\theta<0$ phases on either side. The middle transition is the critical end line of this surface.

Going still lower in temperature, another first order phase transition may occur at some $\alpha_{3}(\varepsilon)$, and forms the lower bound of the coexistence surface. It will be referred to as the lower transition. Below this lower transition, the orientation $\theta$ is fixed at $\pi/2$ (see FIG 2). The necessary and sufficient condition for the lower transition is

[TABLE]

The lower transition occurs at an temperature where $\vec{\eta}$ is no longer small. The sixth order $\gamma_{2}$ term dominates the free energy, and by the same argument $g_{0}$ is overwhelmed by $g_{1}$ and can be ignored. The $\gamma_{2}$ term has six degenerate minima (see equation (4),) and this degeneracy is lifted by $g_{1}$ and $g_{2}$. It will be shown subsequently that the difference $(g_{1}-g_{2})$ determines which of the six minima is most favorable, leading to the criterion (19).

All three lines of phase transition may continue indefinitely into higher stress, or alternatively the middle transition line may bend down and end when it merge with the lower transition at some value of stress $\varepsilon_{*}$. Overall, there are three different possible scenarios associated with this given $\mathcal{F}$ at $\Phi=0$, as shown in FIG 3.

III.2 $\Phi=\pi/6$

Immediately below the upper transition the order parameter has $\theta=\pi/6$, which is again locked by the symmetry. We shall refer to this as phase D. This orientation also minimizes the $\gamma_{2}$ term, however, and there is no middle transition.

On the other hand, the competition between $g_{0}$ and $g_{2}$ may be relevant at low temperature if (19) is satisfied, and then there is a lower transition. Crossing this transition from high to low temperature, the equilibrium state goes from $\theta=\pi/6$ to two-fold degenerate $(\theta-\pi/6)>0$ and $(\theta-\pi/6)<0$. Please see FIG 4 for illustration.

Purely from the symmetry standpoint, this transition can be of either first or second order. It will be shown that, for stress $\varepsilon$ smaller than some critical value $\varepsilon_{c}$, the lower transition is of first order; beyond this point it becomes second order. We will let $\alpha_{4}(\varepsilon)$ denote the line of this lower transition.

Similar to the middle transition at $\Phi=0$, lower transition here serves as the critical end line of a first order coexistent plane between the $\theta>\pi/6$ and $\theta<\pi/6$ phases on two sides of $\Phi=\pi/6$.

There are two possible scenarios at $\Phi=\pi/6$, one with and the other without a lower transition. The phase diagrams are as shown in FIG 5.

III.3 The Full Phase Diagram

A generic value of $\Phi$ breaks the six-fold rotational symmetry, and there cannot be a middle transition. If the (19) is satisfied, the resulting lower transition indeed span a first order coexistent surface, interpolating between the lines of lower transition we have identified at $\Phi=0$ and $\Phi=\pi/6$. Overall there are three possible scenarios, as shown in FIG 6.

As discussed earlier, the magnitude of stress in our discussion is limited by (16). Beyond the small-stress regime, the possible phase diagrams display very rich and complicated behaviors. But as noted above, this regime is unlikely to be physically relevant.

The middle transition at $\Phi=0$ marks the temperature below which $\gamma_{2}$ dominates over $g_{0}$. The lower transition is due to the competition between $g_{2}$ and $\gamma_{2}$. This already exhausts the list of possible competitions. Consequently we do not anticipate any other transition on the phase diagrams.

We have numerically verified the assertions made in this section. In the following sections, we will present analytical derivations the phase diagrams, and analyze their experimental implications.

IV The Split Superconducting Transition

IV.1 Upper Transition

The normal-to-superconducting transition occurs at $\alpha_{1}(\varepsilon)$, as given in (18). Equivalently, one may revert to using the physical temperature:

[TABLE]

is the corresponding critical temperature of the upper transition. This critical temperature is isotropic and independent of $\Phi$.

However, the fact that $\alpha_{1}$ depends linearly in $\varepsilon$ is itself an experimental signature for the nematic superconducting order. A single-component order parameter cannot couple linearly to $\vec{\varepsilon}$, as there is no way to form a combination that is invariant under three-fold rotation.

The specific heat jump will be anisotropic. Following standard analysis, one recover

[TABLE]

This six-fold anisotropy is another unique signature of the nematic superconducting state. It shares similar physical origin with the predicted six-fold anisotropy of $H_{c2}$ Agterberg and Walker (1995a); Krotkov and Mineev (2002); Venderbos et al. (2016b): the external field sees only the underlying lattice as the source of anisotropy. By mapping out the $\varepsilon$ and $\Phi$ dependence of $\Delta c_{v,1}$, one can in principle experimentally determine the Ginzberg-Landau coefficients involved. In particular, the the relative size of $g1$ and $g2$ can be used to predict whether a lower transition may exist.

IV.2 Middle Transition

The existence of a middle transition is another unique signature of the nematic superconducting state. It separates the $\theta=0$ phase from $\theta\neq 0$. Therefore $\theta$ itself is an appropriate order parameter to describe this transition. To this end, it is useful to first minizing $\mathcal{F}$ with respect to variation in $D$. Thus we define $\tilde{D}$ such that

[TABLE]

$\tilde{D}$ is an implicit function of $\theta$, $\alpha$, $\varepsilon$, and $\Phi$. We will assume that the stress orientation is tuned to $\Phi=0$ for the remainder of this part.

One can now derive a Ginzberg-Landau expansion for $\theta$ by replacing $D$ with $\tilde{D}$:

[TABLE]

Then the standard analysis of second order phase transition applies.

As noted in the previous section, the middle transition is due to the competition between $g_{0}$ and $\gamma_{2}$. It can be shown that, for asymptotically small $\varepsilon$, the leading expressions of coefficients $k_{2}$ and $k_{4}$ are independent of $g_{1}$ and $g_{2}$. Here we present these asymptotic formulas, though the exact analytic result can be obtained by carefully retaining all terms.

The critical temperature $T_{2}$ is

[TABLE]

and the specific heat jump is

[TABLE]

The formulas (24) and (25) are accurate in the limit where $(g_{1}\,\varepsilon/\beta)$, $(g_{2}\,\varepsilon/\beta)$, $g_{1}\sqrt{\varepsilon/g_{0}\,\gamma_{2}}$ and $g_{2}\sqrt{\varepsilon/g_{0}\,\gamma_{2}}$ are small.

We also note that $\Delta c_{v,2}$ vanishes in the limit $\varepsilon\rightarrow 0$. This is consistent with the sum rule

[TABLE]

where $\Delta c_{v}$ is the specific heat jump of the superconducting transition at zero stress.

V Lower Transition

The existence of a lower transition can be inferred by looking at the extremely low temperature $\alpha\rightarrow-\infty$ limit:

[TABLE]

The leading term has the six-fold rotational symmetry, favoring all natural minima equally. The sub-leading term lifts this degeneracy.

The sign of $(g_{1}-g_{2})$ controls the qualitative behavior of (27). First assuming a generic value of $\Phi$ that breaks the six-fold rotational symmetry. If $g_{1}>g_{2}$, out of the three minima, the one closest to $\Phi$ wins out. Conversely if $g_{1}<g_{2}$, the minimum furthest away from $\Phi$ becomes the lowest. On the other hand, the order parameter has $\theta=\Phi$ at the upper transition, and then initially drifts toward the nearest of the natural minima as the temperature is lowered. For $g_{1}<g_{2}$, there must be a first order transition separating the low-temperature asymptotic behavior from that just below the upper transition. This justifies (19) as the criterion for a lower transition.

The cases of $\Phi=0$ and $\pi/6$ require special attention. For the sake of clarity, we introduce $\sigma=(\theta-\Phi)$, and consider the range $-\pi/2<\sigma\leq\pi/2$. It can be shown that $\tilde{\mathcal{F}}$ exhibits the following properties at $\Phi=0$ or $\pi/6$:

$\tilde{\mathcal{F}}$ is symmetric under $\sigma\rightarrow-\sigma$ 2. 2.

$\tilde{\mathcal{F}}$ is always stationary at $\sigma=\pi/2$ and [math]. 3. 3.

$\tilde{\mathcal{F}}$ admits at most five stationary points within the range $-\pi/2<\sigma<\pi/2$ 333To see this, one notes that for a given $\sigma$, the relationship between $\tilde{D}$ and $\alpha$ is invertible. Instead of finding an upper bound on number of stationary points along the path of constant $\alpha$ and $\varepsilon$, one may solve for the much easier case of constant $\tilde{D}$ and $\varepsilon$. 4. 4.

Just below the upper transition, $\sigma=0$ is the global minimum of $\tilde{\mathcal{F}}$, and $\sigma=\pi/2$ is the global maximum

First let us look at the $\Phi=0$ case. Coming down from high temperature, $\sigma$ is initially fixed at [math], but drifts away from this high-symmetry direction below the middle transition. The previous argument for generic value of $\Phi$ therefore applies equally. However, below the lower transition, the location of global minimum in this case is pinned exactly at $\sigma=\pi/2$ by the enhanced symmetry.

Now we turn to the case of $\Phi=\pi/6$. As noted in the previous section, the lower transition here may be either first or second order.

If the lower transition is first order, $\tilde{\mathcal{F}}$ must develop a maximum-minimum pair on each side of $\sigma=0$ as the temperature is lowered. This already accounts for five stationary points within $-\pi/2<\sigma<\pi/2$, and $\sigma=0$ must always remain a local minimum. On the other hand, a second order lower transition implies the stationary point at $\sigma=0$ must revert its character as the temperature is lowered. Therefore the criterion separating first and second order transition is whether $(\partial^{2}\tilde{\mathcal{F}}/\partial\sigma^{2})_{\sigma=0}$ changes sign when the temperature is lowered. This gives the critical stress:

[TABLE]

If $\varepsilon>\varepsilon_{c}$, the lower transition becomes second order.

Finally, we will address the question of if and when the middle transition at $\Phi=0$ ends on the surface of lower transition. While a closed-form analytic solution is possible, the full expression is extremely long and unwieldy. We instead supply a recipe here.

The line of middle transition $\alpha_{2}(\varepsilon)$ is implicitly defined by

[TABLE]

Since the phase below the lower transition has exactly $\sigma=\pi/2$ at $\Phi=0$, the equation

[TABLE]

determines the stress $\varepsilon_{*}$ at which the middle and lower transitions meet if a solution exists. If there is no solution, then the line of middle transition extends indefinitely into large $\varepsilon$. Our numerical results indicate that, depending on the actual values of GL coefficients, $\varepsilon_{*}$ can be well within the limit (16), and the ending of the middle transition can be physically relevant.

VI Pre-existing Pinning Field and Signatures of Nematic Superconductivity

As noted in the introduction, there must be a pre-existing SBF in the $M_{x}$Bi2Se3 crystal to explain the pinning of the two-fold anisotropic direction. Thus far, we have assumed that mechanical stress is the only SBF present in our GL analysis. Should this be the case, one has direct control the SBF by application of another external stress, and the predictions in the previous sections can be directly tested.

However, the true nature of the pre-existing SBF is yet unknown, and it may not be mechanical stress or strain after all. In this section, we consider the crystal subjected to both a pre-existing pinning field $\vec{p}$ that is beyond our experimental control, and an artificially applied stress $\vec{\varepsilon}$ as an experimental probe.

The pinning field $\vec{p}$ must also transform under the $E_{g}$ representation. In the same notation developed for stress, $\vec{p}$ can be expressed as

[TABLE]

where the magnitude $p$ is positive. In the absence of stress or any other external SBFs, the orientation of $\vec{\eta}$ is pinned at $\theta=\Phi_{p}$ near the superconducting $T_{c}$.

As discussed in the introduction, the pinning effect is robust up to room temperature. Therefore we find it reasonable to model the $\vec{p}$ as essentially a given background, independent of temperature, near the superconducting $T_{c}$ ($\approx 3$K). While $\vec{p}$ may itself depends on the applied stress, one can formally expand

[TABLE]

and the higher order terms merely renormailzes the GL coefficients.

To our knowledge, the reported two-fold anisotropy always aligns with a lattice direction in all known casesPan et al. (2016); Yonezawa et al. (2017); Du et al. (2017); Smylie et al. (2018). We therefore conjecture that the pinning field, whatever it may be, breaks the lattice symmetry by favoring one of the three two-fold axes in the $\text{Bi}_{2}\text{Se}_{3}$ structure.

The two-fold axis of $H_{c2}$ in transverse field at zero applied stress can be identified experimentally. We can redefine this two-fold direction as the $x$-axis. Following the analysis of Venderbos et al. (2016b), one deduces that the order parameter $\vec{\eta}$ is pinned at either the $(1,0)$ or $(0,1)$ direction. The (unknown) sign of the GL coefficient to a certain gradient term determines the correct choice. The two cases respectively yield:

[TABLE]

Now we consider how the system behaves when both the pre-existing SBF $\vec{p}$ and the applied stress $\vec{\varepsilon}$ are present. We will focus on the split superconducting transition here, i.e. the pair of upper and middle transitions. To this end, one only needs to retain the coupling terms at order $O(\eta^{2})$. The modified GL coupling terms reads:

[TABLE]

where $S$ is given by (6).

We will define the so-called “total SBF” in the bulk of the superconductor as

[TABLE]

with the $\Phi_{b}$ the corresponding orientation angle. The GL free energy under consideration is:

[TABLE]

The ratio $g_{p}/g_{0}$ is unknown. However, one can always make $\vec{\varepsilon}\propto(1,0)$ so that it is parallel to $\vec{p}$. Physically, this amounts to keeping $\varepsilon_{2}=-2\varepsilon_{xy}=0$ while varying $\varepsilon_{1}=\varepsilon_{xx}-\varepsilon_{yy}$. The the change in magnitude $\Delta|P_{b}|$ is linearly dependent on the applied $\varepsilon_{1}$.

VI.1 Upper transition

Following (18) and (20), the critical temperature of the upper (normal-to-superconducting) transition will show a kink when one continuously varies the applied stress:

[TABLE]

This behavior is a unique signature of a nematic superconducting state. If the order parameter is single-component, it cannot couple linearly to either stress or the SBF. One would expect instead a quadratic dependence:

[TABLE]

In either case, the minimum of $T_{1}$ also marks the point where the two-fold direction is tilted by $\pi/2$.

If the kink in (37) can be identified, the predicted six-fold anisotropy purely due to the lattice should be restored at that point. This effect can be seen from the specific heat jump, following (21). Following Venderbos et al. (2016b), the upper critical field $H_{c2}$ should also exhibit the same six-fold anisotropy. This provides further verification for the nematic state.

VI.2 Middle Transition

The existence of the middle transition is another unique signature of the nematic superconducting state. Recall that it is due to the competition between the naturally preferred alignment of $\vec{\eta}$ and the explicit SBF.

We have been assuming that $\eta_{y}$ is the naturally preferred order parameter, or that the GL coefficient $\gamma_{2}$ is positive. The middle transition exists when the orientation of the total SBF obeys $\Phi_{b}=0$, or $\vec{P}_{b}\propto(1,0)$, which is satisfied on only one side of the $|P_{b}|=0$ kink. The second order phase transition may in principle be identified by the jump in specific heat (25).

We have made two assumptions here: the natural preferred direction, and the orientation of the pre-existing pinning field (33). First, it may well be that $\eta_{x}$ is the natural preference instead, but the conclusion that the middle transition exists on only one side of the $|P_{b}|=0$ kink stands unaltered.

It is trickier if the pre-existing $\vec{p}$ does not in fact align with a lattice direction. Should this turns out to be the case, one can instead apply both non-zero $\varepsilon_{1}$ and $\varepsilon_{2}$, keeping the ratio $\varepsilon_{1}/\varepsilon_{2}$ constant. The $\Phi_{b}=0$ condition is satisfied by exactly one value of the magnitude $\varepsilon$, and the middle transition exists only at that point. Nevertheless, just off that exact value, the sharp middle transition is smeared into a crossover, and one can still register a steep increase in specific heat as an experimental signature.

VII Josephson Junction with an S-wave superconductor

In this section, we change gear slightly and discuss the phenomenology of tunneling current between a nematic superconductor and an s-wave superconductor. As pointed out by Yip *et. al.*Yip et al. (1990), one may write down the effective Hamiltonian of a junction purely based on symmetry considerations, in the same spirit of the GL theory. The behavior of the critical current then reflects the symmetry of the nematic order parameter. Here the implicit assumption is that the interface itself has spin-orbit couplingMillis et al. (1988) so as to allow a tunneling supercurrent.

The theory for Josephson tunneling between two s-wave superconductors was worked out in the classic paper by Ambegaoker and BaratoffAmbegaokar and Baratoff (1963) (AB). For a junction between two different s-wave superconductors, right below the lower of the two critical temperatures, the critical current is proportional to $\sqrt{T_{c}-T}$, and shows no strong stress dependence. In this section we explore how the proposed nematic superconducting state gives rise to qualitative differences.

We consider the scenario where an s-wave superconductor (with a much higher $T_{c}$) is attached to the $M_{x}\text{Bi}_{2}\text{Se}_{3}$ sample being tested. The surface contact of the sample is cut perpendicular to the $z$-axis. We assume that the junction is in the tunneling limit.

The two-fold rotation ${C_{2}}^{\prime}$ (and space inversion) in the original $D_{3d}$ group is no longer a valid symmetry on the surface of contact. We will denote the reduced symmetry group as $G$. The junction coupling term in the effective Hamiltonian respects $G$, and of course the $U(1)$ gauge symmetry.

The s-wave superconductor is described by a complex scalar order parameter $\Psi=e^{i\chi}|\Psi|$. We make the gauge choice so that the order parameter $\vec{\eta}$ is real, and $\chi$ represents the phase difference between the two superconductors.

VII.1 The Intrinsic Contribution

Let us first discuss the leading order junction term that is independent of any SBF. One can form $G$-invariant combinations at order $O(\eta^{3})$, but these are not gauge-invariant in their own right, and must be coupled to $\Psi$. And then the time-reversal symmetry demands the coupling terms to be overall Hermitian. The leading coupling term, at $O(\Psi)$ and $O(\eta^{3})$, has the form:

[TABLE]

where $S$ is defined in (6). The coefficient $m_{0}$ must be real to satisfy time reversal symmetry.

Here, we assume that the junction is but a small perturbation to the bulk, and that the system is sufficiently close to the upper transition. One may therefore use the unperturbed solution for the nematic order parameter $\vec{\eta}$. Then the coupling (39) reduces to

[TABLE]

where $\theta$ is the orientation of the order parameter $\vec{\eta}$ in the bulk.

The supercurrent across the junction can be identifiedYip et al. (1990) as $2\,(\partial\!\mathcal{F}_{j0}/\partial\chi)$. The factor of two is due to that the order parameters in the GL theory describe Cooper pairs. Use also the fact that $D\propto|\alpha|$ just below the upper transition, one identify the critical current across the junction:

[TABLE]

In the hypothetical case where no SBF is present, (41) is the only contribution to the critical current. Our default scenario is that $\eta_{y}$ is naturally preferred, or that $\theta=\pi/2$ here. Then (41) shows a $3/2$ power-law temperature dependence that is markedly different from the AB theory. The other scenario is that $\theta=0$, and tunneling is completely disallowed, again a drastic departure from the AB theory. When the total SBF in the bulk $\vec{P}_{b}$ is non-vanishing, $\theta=\Phi_{b}$ in the immediate vicinity of the upper transition.

VII.2 The SBF-assisted Contribution

The quantity $S$ can be replaced by any SBF to yield an invariant coupling term. Here we will impose both a pre-existing $\vec{p}$ and an applied stress $\vec{\varepsilon}$. The coupling term is thus

[TABLE]

Again, $m_{p}$ and $m_{\varepsilon}$ are real by time reversal invariance.

Unfortunately, the ratio $m_{p}/m_{\varepsilon}$ is not necessarily the same as $g_{0}/g_{\varepsilon}$. One defines another “total” SBF for the junction:

[TABLE]

which is in general not aligned with $\vec{P}_{b}$ in the bulk.

Let $\Phi_{j}$ be the orientation angle corresponding to $\vec{P}_{j}$. Again using the bulk solution for $\vec{\eta}$, (42) becomes

[TABLE]

For a non-zero $|P_{j}|$, (44) is dominant in the immediate vicinity of the upper transition, where $D$ is still small. Taking this limit, the critical current is

[TABLE]

If $\vec{\varepsilon}$ and $\vec{p}$ aren’t aligned, the angles $\Phi_{b}$ and $\Phi_{j}$ are different in general. This critical current therefore shows a two-fold anisotropy as one rotates the stress, and has the conventional square-root dependence on temperature. At first sight, this tunneling current may seem indistinguishable from that of the anisotropic s-wave scenario, but there are still some dramatic signatures that are direct consequences of the nematic superconductivity.

VII.3 Experimental Signatures

At zero applied stress and non-vanishing pre-existing pinning field, one has $\theta=\Phi_{b}=\Phi_{j}=\Phi_{p}$. We proceed to discuss the two cases in our conjecture (33) separately.

If $\Phi_{p}=0$, then at zero applied stress, both $I_{c0}$ and $I_{c1}$ vanishes, and there is no tunneling current allowed. An applied stress at a generic direction tilts both $\theta=\Phi_{b}$ and $\Phi_{j}$, and turns on the tunneling current.

Alternatively, one keeps $\varepsilon_{2}=0$ and applies a non-zero $\varepsilon_{1}$ as proposed in the previous section. $\Phi_{b}$ and $\Phi_{j}$ go from [math] to $\pi/2$ as $\vec{P}_{b}$ and $\vec{P}_{j}$ crosses zero, respectively. While this change in $\Phi_{j}$ does not affect the critical currents, the change in $\Phi_{p}$ switches on the tunneling current from zero.

On the other hand, if $\Phi_{p}=\pi/2$, tunneling current without applied stress is initially non-zero. One again applies $\varepsilon_{1}$ while keeping $\varepsilon_{2}=0$. When $\vec{P}_{b}$ crosses zero, the tunneling current is switched off.

Additionally, if one finds a region where $\vec{P}_{j}\approx 0$ while $\Phi_{b}=\pi/2$, then $\mathcal{F}_{j1}$ is anomalously suppressed by the fine-tuning of the applied stress, and the critical current will instead show the $3/2$ power-law temperature dependence of $I_{c0}$.

The above treatment amounts to the direct coupling of the two superconducting bulks. This is admittedly an over-simplification: for unconventional superconductors, surface depairing may occur depending on the exact detail of the gap functionYip et al. (1996). The extent of this surface effect is dictated by the coherence length. Therefore, near $T_{c}$ when the coherence length is large, surface depairing may substantially suppress the tunneling current; the overall temperature dependence will then have a higher powerTanaka and Kashiwaya (1997); Yip (1997) than our simple prediction of $3/2$. Nonetheless, the temperature dependence is clearly distinct from the s-wave-to-s-wave case.

VII.4 $\vec{\eta}$ in $E_{g}$ Representation

As advertised in the introduction, if $\vec{\eta}$ is instead even under space inversion, i.e. in the $E_{g}$ representation, the results in this section need minor modifications. This stems from the fact that $M$, the mirror reflection about the $yz$-plane, acts differently on $E_{u}$ and $E_{g}$. The coupling terms $\mathcal{F}_{j0}$ and $\mathcal{F}_{j1}$ becomes:

[TABLE]

The upshot is that sines are to be replaced by cosines in both $I_{c0}$ and $I_{c1}$:

[TABLE]

The experimental signatures proposed above are still available, but all angles involved are shifted by $\pi/2$. Without an applied stress $\vec{\varepsilon}$, now tunneling current is forbidden at $\Phi_{p}=\pi/2$, and a non-zero $\vec{\varepsilon}$ that tilts $\vec{P}_{j}$ away switch on the current. On the other hand, $\Phi_{p}=0$ initially allows a tunneling current that can be switched off by reversing $\vec{P}_{j}$. The temperature dependence of $I_{c0}$ can be observed by having $\Phi_{b}=0$ and $|P_{j}|\approx 0$.

VIII Conclusion

In this paper, we explore the coupling between the two-component superconducting order parameter $\vec{\eta}$ and the (traceless part of) stress $\vec{\varepsilon}$ in the basal plane, and map out possible phase diagrams allowed by symmetry constraints. Indeed the analysis is not restricted to mechanical stress: any SBF must couple in a similar manner. In later sections we consider the the case where there are a pre-exisiting pinning SBF in the sample, and an applied stress as an experimental probe.

In the presence of a SBF, we found that the superconducting transition splits into two (upper and middle transitions), and there may be another phase transition at an even lower temperature (lower transition), depending on the values of GL coefficients.

The critical temperature (20) of the upper transition does not show any anisotropy, but a six-fold anisotropy for the specific heat jump (21) is predicted as the SBF is rotated. The physics behind is similar to the theoretical angular dependence of upper critical field discussed in the literatureAgterberg and Walker (1995a); Krotkov and Mineev (2002); Venderbos et al. (2016b).

The middle transition only exist when the SBF is aligned along $\Phi=0$; the sharp transition becomes a crossover if the alignment is not exact. We give formulas for the transition temperature (24) and specific heat jumps (25). As the middle transition is connected to the superconducting transition without a symmetry breaking field, we expect that these results from our GL theory are still quantitatively accurate.

The same cannot be said for the lower transition, which is not connected to the superconducting transition without a SBF. The assumption that GL coefficients are temperature-independent may no longer be a quantitatively accurate approximation, and we focus our effort on the qualitative results. The existence of a lower transition hinges on the criterion (19). It is first order almost everywhere, except when $\Phi=\pi/6$ and $\varepsilon>\varepsilon_{c}$ it becomes second order. The critical value $\varepsilon_{c}$ is given in (28).

The observed two-fold anisotropy of $H_{c2}$, and the robust pinning of the two-fold direction over many cycles, suggest that a pre-existing pinning field explicitly breaks the rotational symmetry in the sample. We discuss how, at leading order in GL theory, the pinning field and an externally applied stress combine to form a “total” SBF felt by the sample. By varying the applied stress, one gains access to different regions of the three-dimensional phase diagram worked out above.

We point out two thermodynamics experimental signatures unique to the nematic state. First, the superconducting critical temperature is linearly proportional to the strength of the total SBF, as given in (37), as opposed to the quadratic relation if the order parameter is single-component. Second, the existence of the middle transition, and the associated finite crossover even if the total SBF isn’t exactly tuned to the right orientation, can be observed through calorimetry experiments.

When linked to another s-wave superconductor, the Josephson tunneling current also offers hints to the superconducting pairing symmetry. We discuss how an applied stress can switch the tunneling current on and off if the superconducting order parameter is indeed nematic. The unusual $(T_{1}-T)^{3/2}$ temperature dependence of critical current $I_{c0}$ (41) may also be seen in experiment if the required conditions are met. None of these peculiar behaviors can be seen if the superconductivity of $M_{x}$Bi2Se3 turns out to be s-wave.

We hope our findings here will guide future experimental effort in discerning the pairing symmetry of the $\text{Bi}_{2}\text{Se}_{3}$ family of superconductor, thereby helping to settle the debate on the topological nature of the superconductivity.

Acknowledgements.

The authors would like to thank Anne de Visser for sharing the details of experimental procedures carried out in Pan et. alPan et al. (2016), and Matthew Smylie for discussing the robustness of pinning effect seen in experiments. This work is supported by Ministry of Science and Technology, Taiwan under grant number MOST 107-2112-M-001-035-MY3, and PTH is supported under grant number MOST 107-2811-M-001-045.