The direction of the $d$-vector in a nematic triplet superconductor

We investigate the states of triplet pairing in a candidate nematic superconductor versus typical material parameters, using the mean field theory for two- and three-dimensional tight-binding models with local triplet pairing in the $E_u$ representation. In the two-dimensional model, the system favors the fully gapped chiral state for weaker warping or lower filling level, while a nodal and nematic $\Delta_{4x}$ state is favorable for stronger warping or higher filling, with the $d$-vector aligned along the principle axis. In the presence of lattice distortion, relative elongation along one of the principle axes, ${\bf a}$, tends to rotate the nematic $d$-vector orthogonal to ${\bf a}$, resulting in the nematic $\Delta_{4y}$ state at sufficient elongation. Three-dimensionality is seen to suppress the chiral state in favor of the nematic ones. Our results may explain the variety in the probed direction of the $d$-vector in existing experiments.


I. INTRODUCTION
Time-reversal-invariant (TRI) topological superconductors (TSC) attract sustained interests due to the interesting Majorana zero modes in vortices and itinerant Majorana fermions on the boundary. [1][2][3] It has been established that the key requirement for TSC in inversion symmetric systems is odd-parity pairing. 4,5 Although considerable efforts are made both theoretically and experimentally, [6][7][8][9][10] , definite evidence of a TRI-TSC is yet to be found. A best candidate to date is Cu x Bi 2 Se 3 , 11 which is made by intercalation of Cu into Bi 2 Se 3 . The parent compound is a topological insulator (TI). While the topology of the normal state is not necessary for TSC, 4 the strong spin-orbital coupling (SOC) in such a material makes TSC more likely. The maximum transition temperature T c observed is 3.8K. Early specific heat measurement 12 seems to indicate a full pairing gap. The upper critical field exceeds the Pauli limit, 13 suggesting triplet pairing. Assuming odd parity triplet pairing, a candidate pairing function is φ(k) = τ 2 σ 3 iσ 2 , 4 dubbed ∆ 2 pairing. Henceforth, τ i (σ i ) is the Pauli matrix in the orbital (spin) basis. The two orbitals are derived from the p z orbitals of Bi-atoms in a quintable layer of Bi 2 Se 3 . The ∆ 2 pairing is a nondegenerate representation of D 3d . The d-vector is along z and the SC state is fully gapped. The zero-bias conductance peak observed in the point-contact spectroscopy indicates the existence of unusual in-gap surface states. 14 However, the scanning-tunneling-microscopy (STM) measurement reveals a full gap with no sign of in-gap surface states. 15 With the spectroscopic uncertainty in mind, recent nuclear magnetic resonance (NMR) experiment 16 makes a breakthrough in this field. The observed Knight-shift K develops a two-fold oscillation as a function of the angle of the in-plane applied field H, with strongest (or no) suppression of K below T c for H along (or orthogonal to) one of the Se-Se bonds, the principle axes henceforth. This could be understood if the d-vector of the triplet aligns along a principle axis. Fu realized that the inplane NMR nematicity suggests the pairing function must belong to a doublet representation. 17 For local pairing, the desired pairing function is obvious, φ(k) = τ 2 σ 1,2 iσ 2 , dubbed (∆ 4x , ∆ 4y ) pairing, with the understanding that the d-vector is along and orthogonal to the principle axis, respectively. Note that ∆ 4x leads to nodal SC gap, protected by the D 3d group, while ∆ 4y could become fully gapped in the presence of warping effect in the normal state band structure. 17 While the nematicity is observed by various probes in Cu x Bi 2 Se 3 , the identified direction of the d-vector varies. 16,[18][19][20][21][22][23] Theoretically, 24 direct visualization of the d-vector is possible by quasi-particle interference (QPI) 25 and STM: the leading peak momentum in QPI at sub-gap energies should be along the d-vector, and the STM profile of the vortex at low energies should be elongated also along the d-vector. The agreement between the results in the momentum space (from QPI) and real space (from vortex profile) is a stringent constraint for the nematic triplet.
In real samples, there may be lattice distortions, 26 which could pin the direction of the d-vector. On the other hand, the extent of warping, the filling level, the thickness of the sample, as well as the strength of inter-layer hybridization, may vary from sample to sample, and their roles for the d-vector are to be unravelled. Here we study how the nematic pairing, and the direction of the d-vector in particular, depends on the above typical material parameters, in order to understand the variety in the probed d-vector direction in existing experiments. We use a mean field theory (MFT) based on a tight-binding model, assuming local triplet pairing in the E u representation.
Our main results are as follows. In the two-dimensional (2D) model, the system favors the fully gapped chiral state for weaker warping or lower filling level, while a nodal and nematic ∆ 4x state is favorable for stronger warping or higher filling. In the presence of lattice distortion, relative elongation along one of the principle axes, tends to rotate the nematic d-vector in favor of the nematic ∆ 4y state. In the 3D model, increasing inter-layer hybridization suppresses the chiral state in arXiv:1902.10915v1 [cond-mat.supr-con] 28 Feb 2019 favor of the nematic ones.
The rest of the paper is organized as follows. The model is described in Sec.II, the effect of warping in the 2D model is described in Sec.III, the effect of lattice distortion in Sec.IV, and the effect of inter-layer hybridization in Sec.V. Finally, Sec.VI is a summary of this work.

II. MODEL AND METHODS
According to the first-principles calculations, 27 the conduction and valence bands of Bi 2 Se 3 are superpositions of Se p z orbitals on the top and bottom layers of the unit cell, each of which mixed with it's neighboring Bi p z orbital. The angle-resolved photoemission spectroscopy (ARPES) experiment 28 shows that the band structure of Cu x Bi 2 Se 3 is quite similar to Bi 2 Se 3 . A two-orbital continuum model 17 has been proposed to describe the low energy physics of Cu x Bi 2 Se 3 . Here we use a similar minimal model defined on the layered triangular lattice, which is equivalent to the continuum model at low energy. The free part of the Hamiltonian can be written as, in the momentum space, where ψ k is a four-component spinor describing two orbitals and two spins, Here the summation is over the three inplane translation vectors d 1 = (1, 0), t w is the warping parameter allowed in a D 3d system, induced here by the form factor f 1,2,3 = (1, −1, 1); µ is the chemical potential; and m controls the topology of the normal state band structure. The 3D dispersion is controlled by the interlayer hybridization t z , and in the 2D limit we set t z = 0.
Finally, the coefficient α i = 1 unless specified explicitly otherwise (for the distorted lattice). Throughout this work we use arbitrary units for qualitative purposes, and we fix m = −0.5 for concreteness.
As an example, we show the band structure in the 2D limit in Fig.1(a) for µ = 3.4 and t w = 0.1, which is consistent with the ARPES measurement. 28 The corresponding density of states (DOS) is shown in Fig.1(b). The DOS has a Van Hove singularity near the Fermi level, and this makes the system sensitive to perturbations that we will be addressing.
The NMR nematicity in the SC state strongly implies triplet pairing in a doublet E u channel. While the underlying pairing mechanism is not yet clear, for our purpose it is sufficient and reasonable to start with an effective MF Hamiltonian with attractive pairing interaction in the E u channel, is the Nambu spinor, t means transpose, and V = (∆ 4x , ∆ 4y ) is the d-vector order parameter determined self-consistently by where V s < 0 is the pairing interaction and N is the number of lattice sites. The local triplet is made possible by the antisymmetric pairing between the two orbitals (or orbital singlet). The chiral pairing is defined by V×V * = 0, while the nematic pairing is signaled by V × V * = 0. The Bogoliubov-de Gennes (BdG) quasiparticles, with energy dispersion E k , can be obtained by diagonalizing the single-particle part of H MF . In the following we define the BdG gap at a polar angle as |E k | min for k along the corresponding radial direction in the inplane Brillouin zone. Usually this is just E k on the normal state Fermi surface (FS), but in the present model the unusual pairing may cause |E k | min to deviate from the FS.
The value of n in the converged MF solution depends on the initial condition. Both the chiral state and the six-fold degenerate ∆ 4x state can be captured by a phenomenological Landau free energy in the form With α < 0 and β > 0 as usual, the chiral phase is stable for β < 0 and γ = 0, 30 while the nematic phase is stable for β > 0 and γ < 0. Note the 6-th order γ-term is required to reproduce the discrete six-fold degeneracy in the nematic state. Interestingly, a similar six-fold degeneracy of helical triplet pairing was found theoretically in doped BiH, but a 12-th order term is needed instead. 29 By systematic MF calculations with V s = −1.5 and in the zero temperature limit, we obtain a phase diagram, Fig.2(b), for the pairing state in the (µ, t w ) parameter space. Note a higher value of µ corresponds to a higher electron filling. We see that the chiral state is favorable for lower warping or filling, and the nematic state is favored otherwise. We stress that the nematic state is of ∆ 4x -type. The ∆ 4y -type pairing is not obtained in the MF calculations here. The characteristic BdG gap veresus the polar angle is shown as insets in the respective phase, which is nodeless/nodal in the chiral/nematic phase.

IV. LATTICE DISTORTION
In the previous section, we find the ∆ 4y pairing is not stabilized in the ideal 2D lattice model. However, this type of pairing state was reported in experiments. We then ask how it could be realized by perturbations away from the ideal limit. One possible perturbation is the lattice distortion, which is indeed observed in Sr x Bi 2 Se 3 by X-ray diffraction. 26 By symmetry, the effect of the lattice distortion on the hopping coefficients α 1,2,3 defined in Eq.1 can be decomposed into two symmetry channels, ∆α i ∝ x 2 i − y 2 i and ∆α i ∝ 2x i y i . Here (x i , y i ) = d i is the two components of the first-neighbor bond d i . We concentrate on the d x 2 −y 2 channel. (The other channel can be discussed similarly.) Although all of α 1,2,3 would be changed, we can perform rescaling so that α 2 = α 3 = 1, and α 1 = 1 − γ, to represent the qualitative effect of the strain. A positive γ corresponds to relative elongation of d 1 versus d 2,3 , and vice versa. (The rescaling may also change m and µ, but for a small distortion such changes could be ignored in the leading order approximation.) Fig.3(a) shows the MF solution versus γ in the parent nematic case with (µ, t w ) = (3.4, 0.4). The ∆ 4y component increases rapidly for positive γ, and the ∆ 4x component decreases to zero already for γ ≥ 0.008. We checked that ∆ 4y /∆ 4x remains real, so the pairing is still nematic. We see that the relative elongation along x can rotate the d-vector to the orthogonal direction, in favor of the ∆ 4y -pairing. For small positive γ, the d-vector interpolates between x-and y-directions. Such an intermediate state may have been observed in Ref. 23. Note for a negative γ, or relative compression of the lattice along x, the ∆ 4x state is the only stable one. This is reasonable since the system favors ∆ 4x already in the parent undistorted lattice for (µ, t w ) = (3.4, 0.4).
In the parent chiral case with (µ, t w ) = (3.4, 0.1), the effect of γ on the ∆ 4x,4y components of V is shown in Fig.3(b). The two components coexist within |γ| ≤ 0.05, and we checked that ∆ 4y /∆ 4x remains imaginary. So in this regime the pairing is still chiral. However, only ∆ 4y (∆ 4x ) is left for γ > 0.05 (γ < −0.05), hence the pairing becomes nematic. Although the critical value of γ for the transition between chiral and nematic states is relatively large, we again see that the relative elongation of the x-bonds in the 2D lattice tends to favor ∆ 4y pairing, and vice versa. By symmetry, the conclusion can be generalized to the relative elongation of any one of the principle axes of the triangle lattice.

V. DIMENSION CROSSOVER
It is known that at low doping, Cu x Bi 2 Se 3 has an ellipsoidal Fermi pocket centered at Γ, while at high doping the FS is open and cylinder-like along k z , as observed in the photoemission 31 and de Haas-van Alphen measurements. 32,33 In our tight-binding model, the shape of the 3D FS is controlled by the inter-layer hybridization t z . For small t z , the FS is open and cylinder-like. For t z > 2, the FS becomes closed around the Γ point. We should remark that in the case of an open cylinder-like FS, both types of nematic pairing are topologically trivial, irrespectively of whether they are nodal or nodeless, since the FS encloses two time-reversal invariant momenta. With this in mind, we continue to investigate the pairing state in the 3D limit. We model the 3D lattice by N z = 10 layers of triangle lattices and periodic boundary condition along all translation directions.
As an example of the 3D model, we consider µ = 3.4, t z = 1.0 and V s = −2.5. (In the 3D model, the DOS is lower. To get a sizable T c , a stronger pairing strength is used.) Fig.4(a) shows the FS in the inset, and the evolution of the MF solutions in the main panel, for the real (solid line) and imaginary (dotted line) parts of the ratio r = ∆ 4y /∆ 4x . We find in both cases of t w = 0 (blue lines) and t w = 0.1 (red lines), the solution converges to V = (∆ 4x , ∆ 4y ) ∝ (1, √ 3)/2, namely the ∆ 4x state up to a C 6 rotation.
By systematic calculations, we obtain the phase diagram in the (t z , t w ) parameter space, see Fig.4(b). The chiral state is limited to smaller t z or t w , while the nematic ∆ 4x phase prevails for larger t z or t w . In reality, the hopping along z is between quintuple layers, and should be weak. In this case, the result at large t w , e.g., t w = 0.4, is qualitatively similar to the 2D case discussed in Sec.III.
We note that the phenomenological Landau theory analysis in Ref. 34 reaches a similar conclusion that reduced 3D dispersion would favor the chiral state, however, it does not include the effect of warping, hence is unable to capture the nematic phase in the 2D limit. On the other hand, we have further checked the effect of lattice distortion in the 3D model. The qualitative effect is the same as in Fig.3(a), namely, a relative elongation along x would favor the ∆ 4y pairing.

VI. SUMMARY
To conclude, we studied the pairing states in a model for doped Bi 2 Se 3 superconductor with triplet pairing in the E u representation. In the 2D model, the fully gapped chiral state is favored if the warping parameter is small, while the nodal nematic triplet with the d-vector V = (∆ 4x , 0) is favored otherwise. Under lattice distortion, a relative elongation along x would favor a d-vector V = (0, ∆ 4y ). In the 3D model, the chiral state disappears for large interlayer hopping, in favor of the nematic ∆ 4x state, and the effect of lattice distortion is similar to the case of the 2D model. Taking the above material parameters into account, our results may provide a possible cause of the variety in the probed d-vector direction in existing experiments.