Vestigial nematic order and superconductivity in the doped topological insulator Cu$_{x}$Bi$_{2}$Se$_{3}$

If the topological insulator Bi$_{2}$Se$_{3}$ is doped with electrons, superconductivity with $T_{{\rm c}}\approx3-4\:{\rm K}$ emerges for a low density of carriers ($n\approx10^{20}{\rm cm}^{-3}$) and with a small ratio of the superconducting coherence length and Fermi wave length: $\xi/\lambda_{F}\approx2\cdots4$. These values make fluctuations of the superconducting order parameter increasingly important, to the extend that the $T_{c}$-value is surprisingly large. Strong spin-orbit interaction led to the proposal of an odd-parity pairing state. This begs the question of the nature of the transition in an unconventional superconductor with strong pairing fluctuations. We show that for a multi-component order parameter, these fluctuations give rise to a nematic phase at $T_{{\rm nem}}>T_{c}$. Below $T_{c}$ several experiments demonstrated a rotational symmetry breaking where the Cooper pair wave function is locked to the lattice. Our theory shows that this rotational symmetry breaking, as vestige of the superconducting state, already occurs above $T_{c}$. The nematic phase is characterized by vanishing off-diagonal long range order, yet with anisotropic superconducting fluctuations. It can be identified through direction-dependent para-conductivity, lattice softening, and an enhanced Raman response in the $E_{g}$ symmetry channel. In addition, nematic order partially avoids the usual fluctuation suppression of $T_{c}$.

The electron doped topological insulator Bi 2 Se 3 has been reported to exhibit a low carrier density (n ≈ 10 20 cm −3 ), together with a small ratio ξ/λ F ≈ 2 · · · 4 [1 -3]. Recently, NMR Knight-shift measurements [4] and measurements of the angular-dependent specific heat in the magnetic field [5] revealed spontaneous symmetry breaking of the superconducting state in addition to the global U (1)-symmetry. The threefold symmetry of the underlying lattice is broken. Similar nematic superconductivity was observed in Sr x Bi 2 Se 3 [6,7] and in the closely related Nb-doped Bi 2 Se 3 [8,9]. Early on, Fu and Berg made the proposal that Cu x Bi 2 Se 3 may have an odd-parity two-component superconducting order parameter [10]. Rotational symmetry breaking below T c is then a possible consequence of this pairing state which is, by symmetry, in the E u representation of the point group D 3d , see Ref. [10][11][12][13]. This is the representation that transforms like the in-plane coordinates x = (x, y).
In this paper we show that superconducting fluctuations induce a phase transition to a nematic state. We find that these fluctuations either give rise to a nematic phase transition at T nem > T c or drive the superconducting transition weakly first order. Our quantitative analysis prefers the former scenario, where nematicity is a vestigial precursor phase of superconductivity. This is due to the pronounced two-dimensional electronic structure seen in ARPES measurements [14] that is induced by Cu-intercalation. In distinction to the usual expectation where fluctuations suppress T c , we find that nematic order largely off-sets this suppression, i.e. strengthens pairing compared to the case without nematic order. In the nematic state, the overall superconducting phase averages out to zero, yet the relative orientation of the two components of the Cooper pair field condenses in a long ranged ordered state with broken Z 3 or three-states Potts Figure 1: Schematic phase diagram of CuxBi2Se3 with data points (blue bullets) taken from [3]. We predict a purely nematic phase above the superconducting phase (indicated in red), where superconducting fluctuations create an ordered state that breaks the threefold rotational symmetry. Following Ref. [15], we expect for low temperatures tricritical points (red bullets) below which the transitions should be joint first order. model symmetry at T nem . Superconductivity sets in at a temperature slightly below T nem . The resulting phase diagram for doped Bi 2 Se 3 is sketched in Fig.1. Because of the locking of the Cooper pair wave function to the lattice, the elastic constant c Eg together with the sound velocity along certain high-symmetry directions are reduced at the upper temperature T nem . As the nematic transition of a clean system turns out to be weakly first order, the elastic constant will however not completely vanish. Weak disorder changes the transition to become second order giving rise to a vanishing elastic constant. The nematic state above T c can also be identified through anisotropic paraconductvity T c < T < T nem . We determine this anisotropy from the fluctuation spectrum of the Cooper pair field. The nematic order discussed here has several parallels to spin-induced Ising nematic order above a striped magnetic state of the iron-based superconductors [15][16][17][18] or to time-reversal symmetry breaking proposed for chiral superconductors in the context of SrRuO 4 [19], revealing the universality of the underlying principle of composite or intertwined order [20].
Before we discuss the details of our analysis we summarize the key idea of this paper. The low energy Hamiltonian that describes the superconducting state of doped Bi 2 Se 3 in the band basis is of the form with the pseudo spin index s = {1, 2} and ψ † k = ψ † k,1 , ψ † k,2 denoting fermionic creation operators in the conduction band (cf. Methods section for details). The pairing is given by the d k -vector as where ∆ = ∆ (x) , ∆ (y) T form a two-component order parameter in the E u representation. The broken rotation symmetry below T c naturally implies that this is the appropriate pairing state. Within the D 3d point group, the only alternative would be pairing in the even-parity state E g . This corresponds to d x 2 −y 2 + d yz , d xy + d xz superconductivity. Most of our analysis would proceed without changes if this were the case. Nematic order can generally be characterized in terms of a symmetric trace-less second rank tensor. The nematic tensor in our problem isq where we use the Pauli matricesτ = (τ z ,τ x ). Throughout this paper, we use hat symbols for 2 × 2 matrices and bold symbols for vectors. The expectation value ∆ (µ) * ∆ (ν) measures Cooper pair correlations. While such an expectation value does not break an additional symmetry in the case of a single-component pairing state, we show below that q αβ = 0 breaks another symmetry. For example, q xx = 0 implies ∆ (x) * ∆ (x) = ∆ (y) * ∆ (y) , while in the high temperature phase both expectation values are equal. We show that the additional symmetry is separately broken at a distinct temperature. The emerging nematic phase is not superconducting but induced by superconducting fluctuations. It is a vestige of the superconducting phase. This is only possible because ∆ is a two-component order parameter, i.e. the irreducible representation E u has dimensionality two. To be precise, q αβ transforms according to the representation E g of the point group. Since E u ⊗ E u = A 1g ⊕ A 2g ⊕ E g such a composite order, made up of a bilinear combination, is indeed allowed. Since E u ⊗ E u = E g ⊗ E g the analysis of nematic order does not change for even-parity multi-component superconductivity. The other non-trivial bilinear form is q y ≡ µν ∆ (µ) * τ y µν ∆ (ν) , which transforms under A 2g and breaks time reversal symmetry. In what follows we focus on the nematic order parameter q αβ . It takes the general form The amplitude q 0 sets in at the nematic transition temperature T nem . The unit vector n = (cos θ, sin θ) T is the director of the nematic state that determines the eventual orientation of the superconducting order parameter Thus, the superconducting Cooper pair field acts as nematogen that enables a rotational symmetry breaking, even without superconducting long-range order. Finally, the lattice symmetry of Cu x Bi 2 Se 3 allows for three distinct values of the angle θ = 0, π 3 , 2π 3 (cf. Fig.3). The statistical mechanics of the nematic state then corresponds to a three state Potts model.

Results
Collective nematic fluctuations. Our starting point is the well established microscopic Hamiltonian (23) for Bi 2 Se 3 , yet, for the sake of clarity we refer to the Methods section for a discussion. It leads to the Ginzburg-Landau expansion valid in the vicinity of the superconducting phase transition. In terms of the two-component order parameter ∆ = ∆ (x) , ∆ (y) T the action reads: with r 0 = 1 g − ρ F log ω0 T , where g and ω 0 are the strength and characteristic energy of the pairing interaction, and ρ F the density of states at the Fermi level, respectively. The gradient term is in momentum space given as with m 0 (p) = d p 2 x − p 2 y +d p y p z , and m 2 (p) = 2d p x p y +d p x p z , and characterized by four parameters [13]. We use the shorthand notation´p · · · = T´d 3 p (2π) 3 · · · . From the microscopic Hamiltonian one can determine the coefficients of the Ginzburg-Landau expansion at weak coupling λ = gρ F 1. This analysis yields u > 0 and v > 0, in full agreement with earlier calculations [12]. Crucial for our analysis is however not the applicability of this expansion, but only the signs of u and v. The key implication of the positive sign of v is that the superconducting order parameter is time reversal symmetric and can be written in the form of Eqn. (5). If v < 0 we would find time-reversal symmetry breaking and an analyis analogous to ours leads to a vestigial order parameter q y . The regime v < 0 was predicted in Ref. [21] for thin layers of doped Bi 2 Se 3 . Thus, no matter what the sign of v, one always has an accompanying symmetry breaking. This is true for all crystalline symmetries that allow for multi-component Cooper pair fields. For a superconducting order parameter, transforming according to a higher dimensional irreducible representation Γ, the product representation Γ * ⊗ Γ contains non-trivial irreducible represenations. For all point groups relevant to periodic systems these give rise to vestigial order parameters. Thus, there is no direct second-order superconducting phase transition for a multi-component Cooper pair field. Either there is a vestigial phase above T c , or the transition is of first order due to the coupling to the vestigial order parameter.
For an analysis of fluctuation effects and the description of the nematic ordering, it is efficient to express the interaction in terms of the quadrupolar tensorq αβ of Eqn. (3): where we introduced u = u + v andr = ∆ † ∆ τ 0 . We decouple the two terms in S (4) via Hubbard-Stratonovich transformations, e.g.´DQ e − 1 8v tr(QQ)− 1 2 tr(Qq) ∝ e v 2 tr(qq) , and obtain with the pairing susceptibilitŷ Here, we have expanded the matricesR = Rτ 0 andQ = Q ·τ in terms of the Pauli matrices with Q = (Q 1 , Q 2 ) T . Next, the superconducting order parameter fluctuations are integrated out in both regimes, T < T c and T > T c . To this end, we include Gaussian fluctuations of the pairing field, treated formally within a large-N expansion where T * is the transition temperature without nematic order. The positive value of Tc demonstrates its enhancement due to nematicity as compared to the case with fluctuations but without nematic order. We find a first order nematic phase transition at T = Tnem, followed by a second order superconducting phase transition at T = Tc. We also depict the inverse uniform pairing susceptibility χ −1 pair , which experiences a sudden drop at Tnem before it vanishes at Tc. Finally, we show the inverse nematic susceptibility χ −1 nem that reflects nematic fluctuations above Tnem. Plotted are the dimensionless quantities ∆0 d V 2 0 /(TcVṽ), q0/(2d ṽ), χ −1 pair /(d ṽ) and δr/(d ṽ) with definitions and parameter values given in the Methods section.
of the vector field ∆. We also allow for superconducting symmetry breaking with the condensed pairing field Using the saddle point approximation we finally obtain the five coupled equations of state The saddle point value of the collective variableQ αβ equals the desired order parameter q αβ . Vestigial nematicity is fluctuation induced. Without such superconducting fluctuations, equation (13) does not allow for a finite nematic order, q 0 = 0, above T c . Transition temperatures. The result of the numerical solution of the coupled set of equations (12)- (14) are shown in Fig.2. Here, we plot the order parameters as function of δr = r 0 − r * 0 ∝ (T − T * ), where T * denotes the transition temperature without nematic order present, i.e. r * 0 = −2u . We find that the nematic order parameter sets in above the superconducting transition temperature T c . The superconducting transition is of second order and T c can also be obtained from the divergence of the uniform pairing susceptibility χ −1 pair = (r 2 −Q 2 )/r with r = r 0 +R, see Fig.2. Note that χ pair denotes the largest eigenvalue of the matrixχ p for p = 0. The nematic transition at T nem is weakly first order. The origin of this behavior is the trigonal symmetry which allows for a cubic invariant. Up to fourth order terms, the action for the real order parameter Q reads: This is the well-known Landau expansion of a three states Potts model [22]. Expanding the coupled set of equations (12)-(14) for small Q i yields exactly this term with w ∝d 2 d > 0. Thus, overall the first order transition is expected to be weak, where it holds for the jump of the order parameter δq 0 ∝ w. Given the uncertainty in several parameters we cannot reliably predict T nem − T c in Kelvin. Our numerical analysis suggests however that it can be up to 10% of T c . In Fig.2 we see also the positive effect of nematic order on superconductivity. Fluctuations without nematic order suppress the transition temperature to T * T 0 c ≈ ω 0 e −1/λ as the coupling constant is reduced λ −1 → λ −1 + R/ρ F . With nematic order this effect is significantly weakened as now λ −1 → λ −1 + (R − |Q|) /ρ F . If we assume a more isotropic electronic structure we obtain instead two joint first order transitions, a trend that also occurs in other problems with vestigial precursor order [15]. Photoemission experiments [14] for Cu x Bi 2 Se 3 strongly support a very anisotropic Fermi surface, i.e. split transitions.
Degeneracy of superconducting and nematic ground states. Let us analyze the allowed orientation of the nematic director n, i.e. the allowed values for the angle θ. As shown in Fig.3, the fluctuation induced term Eqn. (15) picks three distinct values of the angle θ = 0, π 3 , 2π 3 , i.e. we have The presence of a finite nematic order parameter predetermines which of the degenerate superconducting ground states of Eqn. (5) will be realized. At T c no additional rotational symmetry breaking takes place and only the global U (1) symmetry of the superconductor and parity are broken. The one-to-one correspondence between the superconducting and the nematic order parameters follows from Eqn. (14), or equivalently, by determination of the angle θ in (4) and (5). We find that Figure 3: Illustration of the correspondence between the nematic and the superconducting order parameters. The green triangles denote the threefold degenerate nematic ground state (16). Each nematic ground state corresponds to two superconducting ground states (blue hexagons) which differ by an overall phase of π from one another, indicated by the purple arrows. The three insets visualize the respective superconducting ground state and the entailed threefold rotational symmetry breaking. The dz(x, y)-component of the order parameter (see Eqn. (2)) is plotted color coded on the ground with the hexagonal atomic structure of the unit cell on top. The gap affects the electronic bonds differently, leading to the aforementioned symmetry breaking.
the first solution of (16) leads to the superconducting ground states with θ n = {0, π}, while the second and third solutions lead to θ n = 2π 3 , 5π 3 and θ n = π 3 , 4π 3 , respectively (see Fig.3). To visualize the origin of the inplane anisotropy in real space, the three insets of figure  3 show the components of the triplet vector d z (x, y)(see Eqn. (2)) of the respective superconducting ground states in real space. We also show the Bi and Se atoms in the respective layers of the crystalline unit cell to demonstrate how the bonds are affected by the anisotropic superconducting gap. For the remainder of the work, we choose, without restriction, the first of the three degenerate nematic solutions, where θ n = 0.
Experimental implications. In the following we study the experimental implication in the nematic phase above T c and in the high temperature phase above T nem . Above T nem the onset of nematicity can be probed via renormalizations of the elastic moduli of the system. The elastic energy relevant for the transition is where we focus on in-plane distortions. The symmetry allowed coupling between the Cooper pair field and the elastic strainˆ αβ is We can now add an external stress to the energy and determine the renormalized elastic constants. Alternatively, we can add a conjugate field to the nematic degrees of freedom and obtain the nematic susceptibilitŷ where the Q i are again the expansion parameters of the nematic tensor in the Pauli basisτ = (τ z ,τ x ) that we have been using. As long as the lattice is purely harmonic we obtain the following relation between the renormalized elastic modulus c * Eg and its bare value c Eg : whereχ nem =χ nem (p → 0). A similar result for spininduced nematicity was previously derived in Ref. [18].
As T → T nem from above the nematic susceptibility rises, leading to a suppression of elastic constants. Within the Gaussian fluctuation regime the nematic susceptibility can be obtained explicitly and is given by: whereχ (0) nem,ij =´p tr (τ iχpτjχp ). In Fig.2 we also show the temperature dependence ofχ −1 nem,ii which displays a Curie-Weiss dependence. Since the transition is first order,χ nem will not diverge, and for the offset holdŝ χ −1 nem (T = T nem ) ∝ w 2 . However, as the first order transition is weak, the nematic susceptibility is significantly enhanced. Note thatχ −1 nem,ii ∝ r Q with r Q occuring in (15) andχ −1 nem,ij = 0 for i = j and T > T nem . To determine κ and the actual lattice softening one would need to know the change in lattice parameters deep in the superconducting state.χ nem is however directly observable via electronic Raman scattering [23,24] in the E g -channel.
Next, we study observables in the nematic phase, i.e. for T c < T < T nem , where the threefold symmetry is broken. As the nematic state is fluctuation induced, the most natural quantity to reflect this anisotropy is the paraconductivity of the system. Our calculation of the fluctuation contribution to the resistivity is a natural generalization of the classical works by Aslamasov and Larkin [25,26]. We obtain the conductivity: with the lattice constants a j , the velocity matrixV j,p = ∂χ −1 p /∂p j and the matrix of the pairing fieldχ p . The calculated temperature dependence of the resistivityρ αα = (σ −1 ) αα is plotted in figure 4. As expected, we find an anisotropy between the two in-plane components. Moreover, the sudden drop at T nem once again indicates the first order nature of the transition and evidences that fluctuation effects are more pronounced inside the nematic phase. For the chosen ground state, i.e. θ n = 0, the resistivity in y-directionρ yy is larger thanρ xx , since the fluctuating pairing amplitude along the x-direction is much larger than in the orthogonal direction, see Fig.3.
Effect of disorder on the nematic phase. Apart from the usual pair-breaking effects, disorder has a profound impact on states with vestigial order (see also [27]). A disorder configuration that locally changes a certain crystalline orientation will naturally nucleate a specific value of the nematic order parameter in its vicinity. Thus, ordinary potential scatters act as random field disorder for the vestigial order parameter, which is according to (15) a three states Potts variable. The random-field threestates Potts model was analyzed in [28,29]. Using the results of these papers the implication for our problem is that disorder changes the first order transition to become second order. Thus, the lattice softening should become more pronounced. Most importantly, weak disorder is not expected to destroy the nematic state.

Discussion
We showed that superconductivity with odd-parity pairing in the doped topological insulator Bi 2 Se 3 is dramatically affected by fluctuations. These fluctuations are important given the low carrier concentration. As a result, the U (1) gauge symmetry and the rotational symmetry are separately broken at the temperatures T c and T nem > T c , respectively. The intermediate nematic state is characterized by strong anisotropic superconducting fluctuations that give rise to an anisotropic paraconductivity. The symmetry breaking will certainly be inherited by the helical surface states of Bi 2 Se 3 . The three states Potts universality class of the vestigial order parameter implies that there should be three distinct domains of vestigial order that can be aligned by applying external stress in the E g symmetry, i.e. for finiteˆ xx −ˆ yy or xy . In addition, a lattice softening and enhanced Raman response are expected above T nem . This mechanism of composite order gives rise to an enhancement of the superconductivity if compared to the usual fluctuation suppression of T c . While the transition temperature should still be smaller than the mean-field temperature, it offers an explanation for the comparatively large transition temperature of Cu x Bi 2 Se 3 , given the low carrier concentration. This observation further suggests to search for similar states of composite order in other low-carrier superconductors with strong spin-orbit interaction.

Methods
Derivation of Ginzburg-Landau expansion. We start from the established electronic structure of Bi 2 Se 3 near the center of the Brillouin zone with the Hamiltonian [30] Here, c k = (c k,+↑ , c k,+↓ , c k,−↑ , c k,−↓ ) refers to the electron annihilation operators for momentum k, located in the two relevant p z -orbitals in the unit cell (±), and with spin (↑↓) [30]. The Pauli matricesτ i andσ j act in orbital and spin space, respectively. The last term h k = −λ k 3 + + k 3 − τ zσz with k ± = k x ± ik y takes into account the point symmetry of the hexagonal lattice. The origin of the Rashba-type spin-orbit interaction in Eqn. (23) is caused by the lack of local inversion symmetry of the Bi-Se layers. Globally, the system is inversion symmetric, hence the coupling toτ z , that is odd under parity, in orbital space.
The superconducting pairing states of doped Bi 2 Se 3 were classified in Ref. [10]. The state that is compatible with a rotational symmetry breaking has odd-parity and gives rise to the expectation value c † in the odd-pairity symmetry channel: The band structure of H 0 gives rise to four bands. In bulk, two pairs of Kramers degenerate bands are separated by a gap. In electron doped Bi 2 Se 3 , such as Cu x Bi 2 Se 3 , the Fermi energy is shifted to the upper two bands. Thus, we follow Ref. [12] and project into the conduction bands. Specifically, we use the manifestly covariant Bloch basis ψ ks , see Ref. [31], that respects the transformation behavior in coordinate and spin space. The index s = {1, 2} refers to the pseudo-spin that labels Kramers degeneracy. It follows H 0 ≈ ks ψ † k,s (ε k − µ 0 ) ψ k,s with ε k = m 2 + v 2 0 k 2 x + k 2 y + v 2 z k 2 z . The pair creation operator in this basis is given by where ψ † k = ψ † k,1 , ψ † k,2 and theσ l are Pauli matrices in pseudo-spin space. The full expressions for the form factors ϕ k,0 = 1 ε k (v z k z , 0, −v 0 k x ) as well as d (y) k,0 = 1 ε k (0, v z k z , −v 0 k y ), see also Ref. [12]. This fully defines the low-energy Hamiltonian of doped Bi 2 Se 3 . For the evaluation of the Ginzburg-Landau parameters in Eqns. (6) and (7) we use the parameters v 0 = 3.3 eVÅ, m = 0.28 eV as in [32], µ = 0.50 eV as in [33], a ≡ a x = a y = 4.1Å as in [34] and T c = 3.8 K. We chose g to reproduce the experimental transition temperature. As there is experimental evidence that the velocity v z depends on the amount of Cu substitution, we kept v z as a tuning parameter. Depending on the choice of v z we found either a joint first order transition, or the scenario depicted e.g. in figure 2 where we used (v z /a z )/(v 0 /a) = 1/20. The plotted dimensionless quantities from figure 2 read ∆ 0 d V 2 0 /(T c Vṽ), q 0 /(2d ṽ), χ −1 pair /(d ṽ) and δr/(d ṽ), where we definedṽ = (vT c )/(d 2 V 0 ) and V 0 = a x a y a z .