Nematicity Arising from a Chiral Superconducting Ground State in Magic-Angle Twisted Bilayer Graphene under In-Plane Magnetic Fields

Recent measurements of the resistivity in magic-angle twisted bilayer graphene near the superconducting transition temperature show two-fold anisotropy, or nematicity, when changing the direction of an in-plane magnetic field [Cao \textit{et al.}, Science \textbf{372}, 264 (2021)]. This was interpreted as strong evidence for exotic nematic superconductivity instead of the widely proposed chiral superconductivity. Counter-intuitively, we demonstrate that in two-dimensional chiral superconductors the in-plane magnetic field can hybridize the two chiral superconducting order parameters to induce a phase that shows nematicity in the transport response. Its paraconductivity is modulated as $\cos(2\theta_{\bf B})$, with $\theta_{\bf B}$ being the direction of the in-plane magnetic field, consistent with experiment in twisted bilayer graphene. We therefore suggest that the nematic response reported by Cao \textit{et al.} does not rule out a chiral superconducting ground state.

Recent measurements of the resistivity in magic-angle twisted bilayer graphene near the superconducting transition temperature show two-fold anisotropy, or nematicity, when changing the direction of an in-plane magnetic field [Cao et al., Science 372, 264 (2021)]. This was interpreted as strong evidence for exotic nematic superconductivity instead of the widely proposed chiral superconductivity. Counter-intuitively, we demonstrate that in two-dimensional chiral superconductors the in-plane magnetic field can hybridize the two chiral superconducting order parameters to induce a phase that shows nematicity in the transport response. Its paraconductivity is modulated as cos(2θ B ), with θ B being the direction of the in-plane magnetic field, consistent with experiment in twisted bilayer graphene. We therefore suggest that the nematic response reported by Cao et al. does not rule out a chiral superconducting ground state.
Recent transport measurements in MATBG provided key features of the pairing symmetry of the superconducting state by revealing a two-fold anisotropy or nematicity in the resistivity around the superconducting transition temperature T c when changing the direction of a relatively-strong ( 0.5 T) in-plane magnetic field [40]. The transport response is still isotropic when the magnetic field is smaller. At first glance it appears that chiral superconductivity should be ruled out since it respects the three-fold rotation symmetry of the lattice. Nematic superconductivity-an exotic phase that breaks the lattice rotational symmetry but respects the transla-tional one-may be favored, which was phenomenologically proposed to be a complicated coexisting phase [41] or intrinsic correlated phase [42]. Nematic fluctuation in the correlated insulating phases was indeed observed in MATBG by scanning tunneling microscopy (STM) [43][44][45], also in twisted double bilayer graphene [46]. But it is not very clear whether the insulating correlated phase is directly related to, or purely competitive with, the superconducting one, in that superconductivity can survive even when the insulating state is completely suppressed [47]. Less attention was, however, paid to the possible role of the in-plane magnetic field for superconductivity. As one marked exception, it was proposed to provide a vector potential to induce the Z 2 symmetry-breaking phase transition in Sr 2 RuO 4 films of (p ± ip)-wave chiral superconductors [48]. However, this was not observed since Sr 2 RuO 4 may not be a p-wave superconductor, as suggested by recent investigations [49][50][51].
In this Letter, we formulate the phase transition of a chiral (d ± id)-wave superconductor driven by a critical in-plane magnetic field in a prototype honeycomb lattice of MATBG and demonstrate that the new phase is nematic with two-fold anisotropy in the transport response. We predict an emerging Hall effect [52] in the paraconductivity of the driven nematic phase, also with a two-fold anisotropy with respect to the direction of the in-plane magnetic field [ Fig. 1(a)]. In detail, the two degenerate chiral states, represented by ξ 1 and ξ 2 at the north and south poles of the Bloch sphere depicted in Fig. 1(b), are coupled via the vector potential of the magnetic field through angular-momentum conservation. When the magnetic field B is larger than a critical one arXiv:2101.01426v2 [cond-mat.supr-con] 21 Jun 2021 B c , the two chiral states are hybridized with equal contributions of the form [ξ 1 + exp(2iθ B )ξ 2 ]/ √ 2, as depicted by points at the equator of the Bloch sphere. The coefficient of this superposition is modulated by the direction of the magnetic field denoted by angle θ B . Near the superconducting transition temperature T c the critical field driving the transition becomes arbitrarily small B c → 0 as T → T c (e.g., B c ∼ 0.6 T when T ∼ 0.9T c ). Although the chiral states ξ 1 and ξ 2 are both isotropic, the driven state given by their superposition is nematic with an anisotropy axis modulated by cos(2θ B ), i.e., showing twofold nematic response with respect to the applied field. The consistence with experimental measurements [40] indicates that chiral superconductivity might be not ruled out in MATBG by these experimental findings. We propose that magnetoelectric transport measurements are useful tools also for other possible (p±ip)-wave chiral superconductors, such as UPt 3 [36][37][38] and UTe 2 [39] thin films, to engineer and identify the pairing symmetry.
Symmetry analysis for magnetic-field driven nematicity.-Although the microscopic mechanism of superconductivity in MATBG may be sensitive to the detailed structure of the flat bands and many-body interaction , the Ginzburg-Landau (GL) phenomenology is largely independent of these details by relying solely on the system's symmetry [1,4]. For MATBG, the (emergent) D 6 symmetry, which contains a six-fold rotation around the normal and two-fold rotations around in-plane axes, is successfully adopted to describe the band structure in the continuum model and tight-binding model with lattice relaxation [25,[53][54][55]. For our purpose, we shall focus on the d-wave superconductivity, which is allowed in D 6 group by the two-dimensional irreducible representation with basis ξ 1,2 in form of (d x 2 −y 2 ± id xy )-waves, with which one can represent the condensate Bose field via introducing the superconducting order parameters ψ 1,2 (r, t) under the basis ξ 1,2 [1,4,35]: The two basis functions ξ 1,2 span a two-dimensional space, depicted as a Bloch sphere in Fig. 1(b), with the d x 2 −y 2 ± id xy states being the north and south poles. To construct the GL Lagrangian, we build the symmetryallowed quadratic, quartic, and gradient terms of order parameters (refer to Supplemental Material for details [56]). We note that symmetry analysis by the D 3 , D 6 and D 6h groups gives the same Lagrangian [1,35,56,57].
The applied static and uniform in-plane magnetic field B induces an effective in-plane vector potential A x,y via an average over the thickness d of MATBG [48], i.e., where · · · = d 0 dz(· · · ) denotes the spatial average over the sample thickness and A = z × B with the surface normal of the film along theẑ-direction. Including this vector potential, the GL Lagrangian density follows as where {α, β 1,2 , γ, λ 1,2 } are real GL parameters, where Λ(r) is an arbitrary function. A constant vector potential is equivalent to a uniform supercurrent, to which case the conclusion drawn by magnetic field can be also applied. Without the magnetic field, the mass term α and stiffness β 1,2 are isotropic for either ψ 1 or ψ 2 , viz., they are isotropic phases without nematicity. The two chiral order parameters are coupled via the orbital effect of a magnetic field. Under the application of an in-plane magnetic field, the ground state can be changed, which is found by min-imizing the free energy Eq. (3), yielding . Without loss of generality, we assume the initial state to be ψ 1 before applying the magnetic field, i.e., a d + id superconductor. After applying the magnetic field, ψ 2 is also mixed into the ground state. However, this admixture is small when the field is weak. On the contrary, when the applied field is sufficiently strong the two order parameters ψ 1 and ψ 2 are driven to be the same in magnitude, and the ground state is no longer a chiral one. We call the complete loss of chirality a phase transition since the symmetry of the phase is completely changed. We note that B → 0 when α → 0 near T c . With the ansatz ψ 2 = ψ 1 e iϕ , we find which is suppressed by orbital effects when β > γ, as is the case in MATBG (see numerical results below). The order parameters with equal contribution of ψ 1,2 can be generally represented bỹ When decomposing the condensate Boson field Φ = ψ 1ξ1 +ψ 2ξ2 [Eq. (1)], we can define new basisξ 1,2 = (ξ 1 ± e iϕ ξ 2 )/ √ 2 forψ 1,2 , respectively, which lie at the equator of the Bloch sphere [ Fig. 1(b)]. Since ψ 2 = ψ 1 e iϕ , the new state driven by the in-plane magnetic field is exactlyψ 1 . We can thereby obtain its Lagrangian by substituting the transformation Eq. (7) into Eq. (8), yielding the linearized free energy ofψ 1 Here we define where the components of the stiffness arẽ c xx = β + γ cos(2θ B ), c xy =c yx = γ sin(2θ B ).
Intriguingly, these components are tunable by the direction of in-plane magnetic field, and are anisotropic along thex-andŷ-directions, viz., correspond to emerging nematicity that breaks the three-fold rotation symmetry in the magnetic-field-driven phase. Nematic paraconductivity.-The field-driven phase that is the ground state below T c shows a nematic transport response as addressed below. Slightly above the superconducting transition temperature, the conductivity, called paraconductivity, is mainly contributed by the superconductor order parameters since their fluctuation under thermal noise can carry a supercurrent [58,59]. The measurement of DC resistivity around T c , as performed in MATBG [40], can thus reflect the fluctuation of the superconducting order parameters and provide information about the pairing symmetry of the superconducting state. To calculate the response to an electric field, the vector potential in H(r, t) is increased by A E = −Et, which is the contribution of the applied electric field [59]. Then by the free energy (8), the time-dependent GL equation, augmented by thermal noise, is written as [58,59] Γ∂ tψ1 (r, t) = −H(r, t)ψ 1 (r, t) + f (r, t), where Γ is the damping rate for the superconducting order parameterψ 1 (r, t) and f (r, t) represents thermal noise. We assume that the thermal noise is white [58,59]. Via a Fourier transformation, the electric current reads where S is the sample area and Λ(q) ≡ ∂H(q, t)/∂A E | A E →0 . From this the paraconductivity is determined in linear response to be which is a tensor that exhibits a Hall response. This emerging Hall effect is unique to the magnetic-fielddriven phase since it is absent in the chiral superconducting phase without magnetic field, which follows an isotropic paraconductivity σ c ij = k B T e 2 Γ/(2π 2 α)δ ij . In particular, when the electric field is applied along thexdirection in the coordinate system defined in Fig. 1(a), the induced current along thex-direction is modulated as c xx , and there is a Hall response with the current along theŷ-direction being modulated asc yx . They are both nematic with a two-fold anisotropy [Eq. (9)].
Parameter estimation.-To be specific for the estimation of GL parameters and paraconductivity, here we consider the tight-binding model on the honeycomb lattice with two p-orbitals {p x , p y } on every site proposed by Yuan and Fu [55,60,61], as shown in Fig. 2(a). The superlattice has a point group of D 3 . Note, however, that our general conclusions rely on the above symmetry analysis only and are thus applicable beyond this specific microscopic model. The chiral basis is denoted by (â k±,σ ,b k±,σ ) T = (â kx,σ ±â ky,σ ,b kx,σ ±b ky,σ ) T / √ 2 on the "A" and "B" sites with σ = {↑, ↓} = {+, −} being the electron spin, the Hamiltonian in momentum space is divided into subspaces, given by Here, µ is the chemical potential; t 1 and {t 2 , t 3 } are the hopping parameters between nearest and fifth neighboring sites that are connected by c µ={1,2,3} and d µ={1,2,3} [ Fig. 2(a)], respectively [60]; f (k) = µ=1,2,3 e ik·cµ and g(k) = µ=1,2,3 e ik·dµ ; and gµ B B is the Zeeman splitting via the in-plane magnetic field. We disregard the Zeeman effect since its suppression of the superconductivity is weak when B → 0, while the orbital effect is still relevant when T → T c . We use the interaction Hamiltonian that allows to stabilize d-wave superconductivity [1,31,34,35,62], where the pairing potential Here N is the number of honeycomb lattice sites, and the coupling constant V < 0. For this specific model, the condensate boson field contains three components Φ = (φ 1 , φ 2 , φ 3 ) T when we consider the pairing between three nearest neighbors on the honeycomb lattice, with which the effective Lagrangian [63] (see Supplemental Material [56] for details), where M µµ determines the superconducting transition, and T µµ δγ controls the spatial fluctuations. In this model, . We first estimate the magnitude of the in-plane magnetic field to realize the phase transition to the nematic superconducting state in MATBG. With parameters |V | = t 1 = 4 meV, t 2 = 0.2t 1 , t 3 = 0.05t 1 and |c µ | = 14/ √ 3 nm for the Moiré honeycomb lattice [5,6], the critical temperature T c of chiral d-wave superconductivity is calculated by solving α(T c ) = 0. T c is shown in Fig. 2(b) and exhibits a dome with a peak at hole doping n h > n s /2 where n s = 4/Ω characterizes doping four holes in one Moiré unit cell of area Ω [5,40]. We note that this T c is well below the BKT transition within our meanfield framework [56]. This peak is not at the Van Hove points of the band that are characterized by two peaks in compressibility dn/dµ with our parameters. With a typical hole doping at µ = −1.5t 1 , we estimate T c ≈ 1.1 K, α = −10 −4 /|c µ | 2 meV −1 · m −2 , β = 2.6 meV −1 , γ = −0.5 meV −1 , and λ 1 = −3λ 2 = 0.02/|c µ | 2 meV −3 · m −2 at temperature T = 1 K [56,60]. With the thickness d ≈ 0.8 nm of TBG estimated by roughly twice the singlelayer one ∼ 0.34 nm [64], the two order parameters ψ 1,2 become close in magnitude when the applied magnetic field B 0.6 T, as shown in Fig. 2(c), which is close to the typical value B c ≈ 0.5 T found experimentally [40].
We then estimate the paraconductivity by choosing Γ = ηα /(k B T c ) with a factor η of order 1 [65,66], leading to a universal paraconductivity σ ij = ησ qcij / β 2 − γ 2 with σ q = e 2 /(2π ) = 3.87 × 10 −5 S being the conductance quantum. With the parameters at µ = −1.5t 1 and T = 1.2 K, we plot σ ij /σ q in Fig. 2(d) with B 0.6 T and η taken to be 1. The conductivities oscillate with respect to the magnetic field with period of π, thus exhibiting a two-fold anisotropy. The amplitude of the oscillation is determined by |γ| that may depend on the behind microscopic mechanism. The direction of the in-plane magnetic field tunes the sign of σ xy and hence the direction of the Hall current, which also could provide an intriguing functionality for future applications.
Discussion.-We have demonstrated nematic paraconductivity that emerges in two-dimensional chiral superconductors under an in-plane magnetic field. This effect is particularly instructive for the chiral d-wave superconducting state of the honeycomb lattice in that the driven phase shows two-fold anisotropy that breaks the three-fold one of the lattice. Furthermore, the magneticfield-driven nematic phase shows a Hall effect in a nonferromagnetic system. The underlying mechanism relies on the hybridization of chiral order parameters by an inplane vector potential with shifted the nodes of the gap in the driven phase, which could be directly tracked by STM [67]. Our work has direct implications for the pairing symmetry of superconductivity in MATBG. Since our purely symmetry-based mechanism applies in the general context of two-dimensional superconducting order parameters it might be relevant to experimental observations in other materials, such e.g., few-layer NbSe 2 reported recently [68,69]