A Landau theory of trigonal two-component superconductors

# Identification of nematic superconductivity from the upper critical field

## Abstract

Recent nuclear magnetic resonance and specific heat measurements have provided concurring evidence of spontaneously broken rotational symmetry in the superconducting state of the doped topological insulator CuBiSe. This suggests that the pairing symmetry corresponds to a two-dimensional representation of the crystal point group, and that CuBiSe is a nematic superconductor. In this work, we present a comprehensive study of the upper critical field of nematic superconductors within Ginzburg-Landau (GL) theory. Contrary to typical GL theories which have an emergent U(1) rotational symmetry obscuring the discrete symmetry of the crystal, the theory of two-component superconductors in trigonal crystals reflects the true crystal rotation symmetry. This has direct implications for the upper critical field. First, of trigonal superconductors with symmetry exhibits a sixfold anisotropy in the basal plane. Second, when the degeneracy of the two components is lifted by, e.g., uniaxial strain, exhibits a twofold anisotropy with characteristic angle and temperature dependence. Our thorough study shows that measurement of the upper critical field is a direct method of detecting nematic superconductivity, which is directly applicable to recently-discovered trigonal superconductors CuBiSe, SrBiSe, NbBiSe, and TlBiTe.

## I Introduction

Unconventional superconductors can be defined by superconducting order parameters that transform nontrivially under crystal symmetries. For a given superconductor, possible unconventional order parameters are classified by non-identity representations of the crystal point group. Such representations are either one-dimensional or multi-dimensional, and this distinction defines two classes of unconventional superconductivity (1); (2). The first class is exemplified by -wave superconductors in cuprates (3); (4), while the second class is exemplified by the -wave superconductivity in SrRuO (5), with two degenerate components at the superconducting transition temperature. Superconducting states in the second class spontaneously break lattice or time-reversal symmetry (6), in addition to the U(1) gauge symmetry, leading to novel thermodynamic and transport properties not seen in single-component superconductors. The search for new superconductors with multi-component order parameters is therefore of great interest.

The doped topological insulator CuBiSe, a superconductor with K (7); (8), has recently attracted a lot of attention as a promising candidate for unconventional superconductivity (9); (10); (11); (12); (13); (14); (15); (16); (17); (18); (19). Fu and Berg proposed that it may have an odd-parity pairing symmetry resulting from inter-orbital pairing in a strongly spin-orbit-coupled normal state (9). While previous surface-sensitive experiments (20); (21) drew disparate conclusions regarding the nature of superconductivity in this material, direct tests of the pairing symmetry in the bulk of CuBiSe are carried out only very recently. A nuclear magnetic resonance (NMR) measurement (22) found that despite the three-fold rotational symmetry of the crystal, the Knight shift displays a twofold anisotropy below as the field is rotated in the basal plane. The twofold anisotropy is also found in the specific heat of the superconducting state under magnetic fields down to T corresponding to (23). Both experiments found that the twofold anisotropy vanishes in the normal state, establishing that the superconducting state of CuBiSe spontaneously breaks the three-fold rotational symmetry. This is only possible when the order parameter belongs to the two-dimensional or representation of the point group. The pairing has been ruled out by comparing the theoretically expected gap structure (24) with specific heat data (8); (23). These results taken together strongly suggest that the pairing symmetry of CuBiSe is , an odd-parity pairing with two-component order parameters (9).

Spontaneous rotational symmetry breaking due to superconductivity is a rare and remarkable phenomenon. Superconductors exhibiting rotational symmetry breaking from multi-component order parameters can be called nematic superconductors (24), in analogy with the nematic liquid crystals and nematic electronic states in non-superconducting metals (25); (26). Nematic and chiral superconductivity, the latter breaking time-reversal symmetry, are the two distinct and competing states of multi-component superconductors, corresponding to real and complex order parameters respectively (6); (1). Broken rotational symmetry has previously been reported in the heavy-fermion superconductor UPt (29) under a magnetic field (28). In addition, the A phase in a narrow temperature range at zero field is likely rotational symmetry breaking, which however may be due to antiferromagnetic order already present in the normal state (31); (30). Thus the recent discovery of broken rotational symmetry in CuBiSe, without broken time-reversal symmetry, may potentially open a fruitful research direction.

Motivated by the recent experimental progress, in this work we study the upper critical field of trigonal nematic superconductors within the framework of Ginzburg-Landau (GL) theory. Such GL theory admits a new trigonal gradient term which is not allowed in hexagonal crystals (27). We relate the gradient terms to Fermi surface and gap function anisotropies by a microscopic calculation of the GL coefficients. Building on and generalizing the previous work (27), we show that the upper critical field generically displays a sixfold anisotropy within the basal plane of trigonal crystals. We further show that a uniaxial strain acts as a symmetry-breaking field in nematic superconductors, which directly couples to the bilinear of the two-component superconducting order parameter. As a result, in the basal plane exhibits a twofold anisotropy with a distinctive angle and temperature dependence, similar to theoretically expected results for UPt in the presence of anti-ferromagnetic order (46); (47). Our findings suggest that measurement of the upper critical field is a direct method of detecting nematic superconductivity. In particular, this method may shed light on the pairing symmetries of other superconducting doped topological insulators SrBiSe (32); (33), NbBiSe (34) and TlBiTe (35), which have yet to be determined.

## Ii Ginzburg-Landau theory

We start by constructing the GL theory of odd-parity two-component superconductivity in crystals with point group and strong spin-orbit coupling. The pairing potential , which is a -dependent matrix in spin space, takes the following form

 ^Δ(→k)=η1^Δ1(→k)+η2^Δ2(→k). (1)

The pairing potential is a linear superposition of two degenerate components , the basis functions of the two-dimensional pairing channel (specific gap functions are given in the Supplementary Material, Sec. III). For odd-parity superconductors the pairing components satisfy . As basis functions of , the two partners transform differently under the mirror symmetry , i.e., is even whereas is odd. A key property of (doped) BiSe materials is strong spin-orbit coupling that locks the electron spin to the lattice. The two complex fields define the superconducting order parameters . In contrast, in case of triplet superconductors in spin-rotation invariant materials the order parameter components are vectors in spin space.

The GL theory of two-component superconductivity is formulated in terms of the order parameters and the GL free energy is the sum of a homogeneous term and a gradient term given by , where and are the corresponding free energy densities. In addition, the free energy contains a Maxwell term , which for our purposes can be taken as a constant. The free energy densities and are polynomial expansions in the order parameter fields and their gradients, and consist of all terms invariant under the symmetry group of the crystal. For two-component trigonal superconductors the homogeneous contribution is the same as the corresponding expression for hexagonal symmetry (6); (1),

 fhom=Aη†η+B1(η†η)2+B2|η∗1η2−η∗2η1|2, (2)

to fourth order in , and we have defined . The coefficients and are phenomenological constants of the GL theory. The sign of GL coefficient determines the nature of the superconducting state, selecting either chiral or nematic order (24); (36).

Spatial variation of the superconducting order parameter is captured by the gauge-invariant gradient , with the electromagnetic vector potential and . In case of multicomponent order parameters, there generally exist multiple independent gradient terms which are allowed by crystal symmetry. It is insightful to present all gradient terms in order of “emergent symmetry”. For crystals with a principal rotation axis along the direction, such as the three- and sixfold rotations of trigonal and hexagonal crystals, four gradient terms with full continuous in-plane rotational symmetry are present and given by (1); (38); (39)

 fD=J1(Diηa)∗Diηa+J2ϵijϵab(Diηa)∗Djηb+J3(Dzηa)∗Dzηa+J4[|Dxη1|2+|Dyη2|2−|Dxη2|2−|Dyη1|2+(Dxη1)∗Dyη2+(Dyη1)∗Dxη2+(Dxη2)∗Dyη1+(Dyη2)∗Dxη1] (3)

(summation understood, , ), and are the phenomenological GL coefficients. The first three terms are invariant under independent U(1) rotation of coordinates and order parameters, and thus have an emergent U(1)U(1) symmetry, whereas the gradient term with coefficient is invariant under arbitrary joint rotations of coordinates and order parameters, i.e., an emergent U(1) symmetry. Therefore, does not reflect the discrete rotational symmetry of the crystal. However, a new gradient term , which we call trigonal gradient term, is uniquely present in crystals with trigonal symmetry, but not allowed in hexagonal crystals (27). It is given by the expression

 fD,trig=J5[(Dzη1)∗Dxη2+(Dzη2)∗Dxη1+(Dzη1)∗Dyη1−(Dzη2)∗Dyη2+c.c.]. (4)

The appearance of this new gradient term, which has symmetry, can be understood from angular momentum, since in trigonal symmetry is equivalent to . Indeed, in momentum space () the trigonal gradient term can be expressed as , where and similarly for . The relative phases between () and () are determined by mirror symmetry: () is even (odd) under . It follows from the structure of that the spatial variation of the order parameter in the basal plane is coupled to spatial variation in the -direction, which is in sharp contrast to hexagonal and tetragonal crystals. In the rest of this work we map out the consequences of trigonal crystal anisotropy in the GL theory for the upper critical field.

## Iii Upper critical field in the basal plane

The angular dependence of was first proposed as a method to establish the multicomponent nature of unconventional superconductors in the context of heavy-fermion superconductors (41); (42); (43). The key idea is as follows. For the class of single-component (e.g., -wave) superconductors with trigonal, tetragonal, and hexagonal symmetry, is always isotropic within the GL theory, due to the emergence of rotational symmetry to second order in the gradients. In case of multicomponent superconductors, effects of crystal anisotropy can appear in the GL theory, removing the emergent U(1) symmetry, but this crucially depends on crystal symmetry. For instance, hexagonal systems with multicomponent order parameters do not show in-plane -anisotropy due to the emergent rotational symmetry of Eq. (3), whereas tetragonal symmetry can give rise to an angular dependence of with fourfold symmetry (42). In trigonal crystals, can exhibit a sixfold anisotropy in the basal plane (27) as of Eq. (4). Here we map out the basal plane upper critical field of trigonal superconductors for general GL gradient coefficients.

Within GL theory, the upper critical field is calculated by solving the GL equations obtained from , keeping only terms linear in since the order parameter is small at . Therefore, the calculation also applies to chiral superconductors. The resulting system of GL equations, which is given by

 −Aηa=J1(D2x+D2y)ηa+J3D2zηa+J2ϵab[Dx,Dy]ηb+J4[(D2x−D2y)τzab+{Dx,Dy}τxab]ηb+J5[{Dz,Dx}τxab+{Dz,Dy}τzab]ηb, (5)

can be solved as a two-component harmonic oscillator problem, leading to a Landau-level spectrum from which is determined as the lowest Landau-level solution. The coupling of the two harmonic oscillators is determined by the structure of the GL equations, and is in general complicated by the presence of multiple gradient terms. In hexagonal and tetragonal systems, straightforward or even exact analytical expressions for can be found (42). In contrast, the trigonal gradient term of Eq. (4) couples basal plane gradients to gradients in the orthogonal direction, giving rise to a different set of harmonic oscillator equations to which previous methods do not apply. A special limiting case was considered in Ref. (27). We generalize this result by solving the GL equations in the presence an in-plane magnetic field for general gradient coefficients. In deriving the general solution we adopt an operator based approach and exploit that harmonic oscillator mode operators corresponding different cyclotron frequencies can be related by squeezing operators. Here we present and discuss the main results, and give a detailed account of the lengthy calculations in the Supplemental Material (SM). For convenience, below we will refer to the appropriate section of the SM.

To demonstrate the key features of in trigonal crystals, we will focus the discussion on the most physical case, where trigonal anisotropy effects may be considered weak and can be treated as perturbation. We take the magnetic field in the basal plane to be given by , which corresponds to a vector potential . It is convenient to rotate the basal plane GL gradients according to the transformation

 (D∥D⊥)=(cosθsinθsinθ−cosθ)(DxDy), (6)

such that is along the field and is perpendicular to the field. These operators satisfy , and and define the magnetic algebra . Writing Eq. (5) in terms of and , and setting (i.e., no modulation along the field), one obtains

 −Aηa=(J1D2⊥+J3D2z)ηa−J4D2⊥(cos2θτzab+sin2θτxab)ηb+J5{Dz,D⊥}(−cosθτzab+sinθτxab)ηb. (7)

Next, it is convenient to diagonalize the term proportional to . This is achieved by a the rotation of the order parameters given by

 (η1η2)=(cosθ−sinθsinθcosθ)(f1f2). (8)

In terms of the rotated order parameters the GL equations read

 −A(f1f2)=(J3D2z+(J1−J4)D2⊥00J3D2z+(J1+J4)D2⊥)(f1f2)+J5{Dz,D⊥}(−cos3θsin3θsin3θcos3θ)(f1f2). (9)

Note that only the term proportional to depends on the angle . We now describe solutions to Eq. (56) obtained by treating as a perturbation.

To start, let us consider taking both . Solving the GL equations then yields two degenerate series of Landau levels with cyclotron frequency , with the upper critical field given by (more details are provided in Sec. II B of the SM). Including the gradient contribution in Eq. (3) proportional to simply makes the cyclotron frequencies inequivalent, , and increases the upper critical field to . This defines the exactly solvable unperturbed system. Then, introducing trigonal perturbation parametrized by couples the two series of Landau levels with different frequencies in a nontrivial way: the coupling of in-plane and out-of-plane gradients implies a coupling of canonically conjugate operators of the form . To solve the system of GL equations we assume that crystal anisotropy effects are weak and use second order perturbation theory to obtain the correction to the cyclotron frequency . (The calculations are lengthy and described in detail in Sec. II B 3 of the SM.) The upper critical field then becomes with . We find to lowest order in as

 Hc2(θ)˜Hc2=1+J252ω2+⎡⎣cos23θ(1−ω−ω+)2+sin23θ1−ω−ω+F(ω−ω+)⎤⎦, (10)

where the frequencies are defined as . In the limit of small these frequencies become and . The function arises due to the coupling of two series of Landau levels with different cyclotron frequencies and oscillator eigenfunctions. It takes the form

 F(x) = (1−x2)52x2∞∑m≥0(2m)!(m!)24mx2m(2m−x21−x2)22m+x(2m+1) (11) = 1−xx[√1+x1−x2F1(12,a2;1+a2;x2)−1]

where and is a hypergeometric function. The function has the property , which implies that for (corresponding to ) no angular dependence of exists. The latter is a consequence of an emergent rotational symmetry of in Eq. (4): it is invariant under in-plane rotations of the order parameters and coordinates according to , . (Note that this is not a physical symmetry.)

In general, however, considering all regimes of gradient coefficients that satisfy the stability constraints of the free energy, exhibits a six-fold anisotropy in the basal plane of the crystal. For instance, the sixfold -anisotropy can be obtained starting from a solution of the GL equations derived from Eqs. (3) and (4) for and , and treating as a small perturbation. This case was considered in Ref. (27) and is described in Sec. II B 2 of the SM.

Figure 1 shows the angular dependence of the upper critical field for small to moderate and as obtained from Eq. (10). Note that in general, for materials with weak to moderate (crystal) anisotropy effects, one expects . To make the interplay between and explicit, we expand Eq. (10) for small and find

 Hc2(θ)Hc2(π4)=1+hcos6θ, (12)

where . This expression serves to highlight an important feature of the angular dependence of : , which is independent of system specific parameters. Here is defined by an axis orthogonal to a mirror plane.

Within weak coupling, the GL coefficients can be obtained in terms of Fermi surface and gap function properties using a microscopic mean-field Hamiltonian with pairing potential given by Eq. (1). The gradient coefficients , , , and are proportional to , where , , and are the Fermi energy, Fermi velocity, and correlation length respectively, and is the density of states. (The calculations are presented in detail in Sec. III of the SM.) We find that their relative strength depends on the crystal anisotropy of the Fermi surface and of the gap functions . In particular, is nonzero only when trigonal Fermi surface anisotropy is present, or when the gap function is composed of trigonal crystal spherical harmonics of the pairing channel (see Sec. III A of the SM), and is generally expected to be weak.

The general sixfold basal plane anisotropy of is a direct consequence of trigonal symmetry and a discriminating characteristic of two-component pairing symmetry. Indeed, single-component superconductivity corresponding to one-dimensional pairing channels of point group cannot exhibit sixfold anisotropy: the in-plane gradient term is given by and has emergent U(1) rotational symmetry. As a result, the sixfold anisotropy provides a clear experimental evidence for two-component pairing.

## Iv Nematic superconductivity and upper critical field

Within our GL theory, the rotational symmetry breaking superconducting state reported in Refs. (22); (23) corresponds to a real order parameter, i.e., . Up to fourth order [see Eq. (2)], the angle represents a continuous degeneracy. This degeneracy is lifted at sixth order by a crystal anisotropy term and leads to a discrete set of degenerate ground states (24); (36). In materials, such as CuBiSe, the remaining degeneracy may be further lifted by a symmetry-breaking pinning field, selecting a unique ground state. The origin of such pinning can be strain-induced distortions of the crystal (37), but in principle, any order with the same symmetry, electronic or structural, can pin the order parameter. In case of two-component superconductors, the symmetry-breaking (SB) pinning field couples linearly to order parameter in the following way

 fSB=g[(uxx−uyy)(|η1|2−|η2|2)+2uxy(η∗1η2+η∗2η1)], (13)

with coupling constant . The order parameter bilinears constitute a two-component subsidiary nematic order parameter (24) with the same symmetry as the symmetry-breaking field . For comparison, uniaxial strain in single-component superconductors couples to the gradient of the order parameter , taking the form different from Eq. (13). It is worth noting that the coupling considered here differs from the candidate theories proposed for the hexagonal superconductor UPt, in which case magnetic order couples quadratically, instead of linearly, to order parameter bilinears  (31); (44); (45); (38); (39); (29).

From a microscopic perspective, the origin of the order parameter pinning in Eq. (13) can be understood as a (strain-induced) Fermi surface distortion, leading to different Fermi velocities . A uniaxial distortion of this form couples to and has the effect of selecting either or by raising , resulting in a split transition. A quantitative calculation of the coupling constant , relating the order parameter bilinear to such Fermi surface distortion can be obtained within weak-coupling (see (40)). This effect of a Fermi surface distortion should be compared to uniaxial gradient anisotropies such as and , with the effect of the former being enhanced by a factor of  (40), where is the coherence length, , is a cutoff frequency, and is an effective interaction energy scale associated with the pairing. In addition, the effect of a uniaxial Fermi surface distortion on the shift of is enhanced by .

To address the effect of the SB field on in case of the trigonal nematic superconductors, we solve the linearized GL equations for small gradient coefficients in the presence of a uniaxial symmetry breaking term defined as , taking as a measure of the uniaxial anisotropy. Here we focus the discussion on the most salient features, for which we take , and relegate a more detailed account to the SM. A similar problem of upper critical field anisotropy was studied for split transitions in UPt (46); (47).

Setting in Eq. (7) and adding the contribution from the symmetry breaking field, the GL equations take the form

 −Aηa=(J1D2⊥+J3D2z)ηa+δτzabηb−J4D2⊥(cos2θτzab+sin2θτxab)ηb. (14)

The upper critical field is obtained by using the magnetic algebra of and , and projecting into the lowest Landau level. The upper critical field is then determined from the following implicit equation (see Sec. II D of the SM)

 −Aω=1l2b− ⎷J24J234ω4l4b+δ2ω2−J4J3δω3l2bcos2θ, (15)

(recall ) where the magnetic length is defined as . For we recover the result for in Eq. (10), to first order in (i.e., expanded to first order in ). For we simply find [see Eq. (12)], but with critical temperature with . This follows from comparing to , i.e., shifts the transition temperature and can be taken as a measure of . We define a dimensionless temperature by .

For general and nonzero we solve Eq. (15) for and show the representative results for and in Figs. 2(a) and 2(b). Two key characteristics of in the presence of a pinning field are evident in Fig. 2(a)–(b). First, the angular dependence of exhibits a distinct two-fold anisotropy, with a typical “peanut”-shape close to . This twofold anisotropy becomes more pronounced with increasing , as shown Fig. 2(b). Expanding the square root in Eq. (15) under the assumption of very small fields, i.e., , one finds (see Sec. II D of the SM). This “peanut”-shape of the profile should be contrasted with the profile of single-component superconductor where uniaxial gradient anisotropy leads to a weak elliptical angular dependence of , an effect which is parametrically smaller than the twofold anisotropy in the two-component case. Consequently, the twofold anisotropy of shown in Fig. 2, in particular the “peanut”-shape, is a discriminating property of two-component pairing.

Second, the angular dependence of is a function of temperature and has different shape in the vicinity of (i.e., small fields) as compared to far below (and high fields). This is in sharp contrast to the usual case, for instance Eq. (10), where only the overall magnitude of is temperature dependent. The unusual temperature dependence of can be more precisely captured by considering the upper critical field anisotropy ratio as function of temperature. In the vicinity of , the anisotropy ratio should exhibit temperature independent behavior given by (see Sec. II D of the SM). This is shown in Fig. 2(c), where the -anisotropy ratio is plotted for various values of . In contrast, using Eq. (15) we find that the -anisotropy ratio approaches unity for large temperature according to , which is independent of GL parameters. Within the model of Eq. (15), the temperature at which the transition between two behaviors occurs is given by . This “kink” feature was also found and discussed in the context of a hexagonal applicable to UPt (46); (47); (48). The distinctive temperature dependence of -anisotropy is uniquely associated with two-component pairing since single-component pairing with uniaxial gradient anisotropy leads to temperature independent -anisotropy.

## V Discussion and conclusion

To summarize, in this work we have addressed the magnetic properties of two-component superconductors in trigonal crystals with point group symmetry. Starting from a general GL theory of trigonal two-component superconductors, we find that the upper critical field exhibits a sixfold anisotropy in the basal plane, which is a discriminating property of two-component pairing. The sixfold anisotropy is a rare manifestation of discrete crystal symmetry, since effects of crystal anisotropy are typically obscured in GL theory by an emergent U(1) rotational symmetry. In addition, in this work we show that when a symmetry breaking field originating from, e.g., structural distortions selects a real order parameter, exhibits a twofold anisotropy with characteristic angular and temperature dependence.

The recent NMR and specific heat measurements on CuBiSe, which reported spontaneously broken rotational symmetry, indicate that this material belongs to the class of superconductors with two-component pairing symmetry. Prominent other examples of materials with trigonal symmetry, which have attracted increasing attention recently, are the doped BiSe superconductors SrBiSe, NbBiSe, and TlBiTe. Our theory of in-plane anisotropy of upper critical field stands to contribute to uncovering the pairing symmetry of these superconductors, which remains to be determined.

## Vi Acknowledgments

We thank Anne de Visser, Shingo Yonezawa, Yoshi Maeno for discussions. This work is supported by the David and Lucile Packard Foundation (L.F.), the DOE Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award No. DE- SC0010526 (V.K.), and the Netherlands Organization for Scientific Research (NWO) through a Rubicon grant (J.V.).

Supplemental material for

“Identification of nematic superconductivity from the upper critical field”

Jörn W. F. Venderbos, Vladyslav Kozii, and Liang Fu

Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

## Appendix A Landau theory of trigonal two-component superconductors

A Ginzburg-Landau (GL) theory of two-component superconductivity in a trigonal crystal with point group symmetry is formulated in terms of the superconducting order parameters, which in turn are obtained from pairing potential. In the presence of spin-orbit coupling, when spin-rotation symmetry is broken, and the symmetry group is the symmetry group of the crystal lattice, the superconducting pairing potential is decomposed into irreducible representations of . In most cases, one is interested in a single pairing channel, i.e., a single representation of the symmetry group, which may be one- or multi-dimensional. The pairing matrix corresponding to pairing channel takes the form (50); (49); (51)

 ^Δ(→k)=∑mηΓ,m^ΔΓ,m(→k), (16)

where labels the representation (i.e., pairing channel) and labels the components of the representation. The expansion coefficients are the superconducting order parameters and are complex scalars. In case of two-component superconductors the order parameter is the two-component complex number .

The GL free energy functional describing the two-component superconductor is then obtained as an expansion in the order parameters to given order,

 Missing or unrecognized delimiter for \left (17)

which consists of terms fully invariant under the symmetries of the crystal. Therefore, the GL free energy depends on the crystal symmetry group as well as the pairing channel. A trigonal crystal with point group symmetry admits two two-component representations, (even parity) and (odd parity), and in both cases the free energy density up to fourth order is given by

 f=F/V=A(|η1|2+|η2|2)+B1(|η1|2+|η2|2)2+B2|η∗1η2−η1η∗2|2, (18)

in terms of the expansion coefficients and . When the GL coefficient is positive () the order parameter is real and breaks rotational symmetry. In that case, a continuous degeneracy remains at fourth order, which is lifted by crystal anisotropy effects at sixth order. The sixth order contribution to the free energy density reads as

 f(6)=C1[(η∗+η−)3+(η∗−η+)3]+C2(|η1|2+|η2|2)3+C3(|η1|2+|η2|2)|η∗1η2−η1η∗2|2, (19)

(where ) and it is the first term, proportional to , which is responsible for lifting the continuous degeneracy and selecting three degenerate ground states related by threefold rotation.

We extend the GL theory to include terms representing spatial inhomogeneity: the gradient terms. The superconductor is a charged superfluid with charge , and gradient terms are introduced by defining the covariant derivative as , where is the electromagnetic vector potential.

The gradient terms of the free energy expansion are obtained using the same representation theory recipe as before. The derivative components transform as , and transforms as of the point group. These representations are used to form products with , from which bilinears with the general structure

 DiXiaηa, (20)

are obtained. Here is a tensor transforming irreducibly under crystal symmetry (see Table 1). A gradient term is then simply given by , with an independent phenomenological stiffness constant for each distinct representation.

We now discuss all gradient terms of an odd-parity two-component superconductor in a trigonal crystal with symmetry. We start from the gradient terms which have a continuous U(1) rotational symmetry with respect to the coordinates and order parameters individually, thus giving rise to an emergent U(1)U(1) symmetry. They are given by

 fD,U(1)×U(1)=J1(Diηa)∗Diηa+J2ϵijϵab(Diηa)∗Djηb+J3(Dzηa)∗Dzηa, (21)

with gradient coefficients (and a sum over repeated and is understood, and ). As a result of the emergent U(1)U(1) symmetry, these gradient terms obscure the true discrete crystal rotation symmetry.

Whether or not the true crystal symmetry is reflected in the GL gradient expansion of multicomponent order parameters depends on the crystal system. For instance, in case of hexagonal symmetry (i.e., point group ) the gradient terms are obtained by considering all irreducible tensors of the hexagonal group. For the two-component representations the gradient terms are given by (50)

 fD=K1|Dxη1+Dyη2|2+K2|Dxη2−Dyη1|2+K3(|Dxη1−Dyη2|2+|Dxη2+Dyη1|2)+K4(|Dzη1|2+|Dzη2|2). (22)

A very similar result is obtained for the representations. This can be rewritten to obtain (see also, e.g., (52); (53))

 fD,U(1)×U(1)[J1,J2,J3]+fD,4[J4]=J1(Diηa)∗Diηa+J2ϵijϵab(Diηa)∗Djηb+J3(Dzηa)∗Dzηa+J4(τzijτzab+τxijτxab)(Diηa)∗Djηb, (23)

where summation over repeated indices is understood. The gradient constants and are related by

 J1=K1+K22+K3,J2=K1+K22−K3,J3=K4,J4=K1−K22. (24)

The term proportional to has an emergent U(1) rotational symmetry: it is invariant under joint rotations of coordinates and order parameters. (Note that the gradient term for symmetry is given by , which also possesses an emergent U(1) symmetry.) Consequently, effects of crystal symmetry do not appear in the GL theory of hexagonal two-component superconductors.

When crystal symmetry is lowered to trigonal symmetry, an additional gradient term arises which reflects the true symmetry of the crystal. Table 1 lists all irreducible tensors and their symmetries. With the help of this table it is straightforward to obtain all gradient terms. One observes that there are two distinct tensors with symmetry. As a result, a cross term of these two is an allowed gradient term. Specifically, the additional gradient term, defined as and referred to as trigonal anisotropy (gradient) term, takes the form (56)

 fD,trig[J5] = J5[(Dxη2)∗Dzη1+(Dxη1)∗Dzη2+(Dyη1)∗Dzη1−(Dyη2)∗Dzη2+c.c.], (25) = J5[τxab(Dxηa)∗Dzηb+τzab(Dyηa)∗Dzηb+c.