Pairing in graphene-based moir\'e superlattices

TL;DR

This paper systematically classifies and analyzes possible superconducting pairing states in graphene moiré superlattices, highlighting the conditions for singlet-triplet mixing, effects of symmetries, and implications for exotic superconductivity.

Contribution

It provides a comprehensive theoretical framework for understanding pairing instabilities and symmetry effects in graphene-based moiré superconductors, including singlet-triplet mixing mechanisms.

Findings

Triplet superconductivity signs in twisted double-bilayer graphene.

Conditions for singlet and triplet pairing coexistence and mixing.

Emergence of nematic mixed singlet-triplet states near mean-field parameters.

Abstract

We present a systematic classification and analysis of possible pairing instabilities in graphene-based moir\'e superlattices. Motivated by recent experiments on twisted double-bilayer graphene showing signs of triplet superconductivity, we analyze both singlet and triplet pairing separately, and describe how these two channels behave close to the limit where the system is invariant under separate spin rotations in the two valleys, realizing an SU(2)$_+$ $\times$ SU(2)$_-$ symmetry. Further, we discuss the conditions under which singlet and triplet can mix via two nearly degenerate transitions, and how the different pairing states behave when an external magnetic field is applied. The consequences of the additional microscopic or emergent approximate symmetries relevant for superconductivity in twisted bilayer graphene and ABC trilayer graphene on hexagonal boron nitride are described in detail. We also analyze which of the pairing states can arise in mean-field theory and study the impact of corrections coming from ferromagnetic fluctuations. For instance, we show that, close to the parameters of mean-field theory, a nematic mixed singlet-triplet state emerges. Our study illustrates that graphene superlattices provide a rich platform for exotic superconducting states, and allow for the admixture of singlet and triplet pairing even in the absence of spin-orbit coupling.

Tables

Table 1. Table 1: Summary of the different intervalley pairing states transforming under the trivial representation of the point group C 3 subscript 𝐶 3 C_{3} in the absence of a magnetic field. For notational convenience, we neglect the extra label m v = 0 subscript 𝑚 𝑣 0 m_{v}=0 to indicate intervalley pairing. λ 𝒌 subscript 𝜆 𝒌 \lambda_{\boldsymbol{k}} is a real-valued and Brillouin-zone-periodic function that is invariant under C 3 subscript 𝐶 3 C_{3} . To lowest order, we can take λ 𝒌 subscript 𝜆 𝒌 \lambda_{\boldsymbol{k}} to be independent of 𝒌 𝒌 \boldsymbol{k} . We also indicate the minimal number of nodes, which state it transforms to when setting J = 0 𝐽 0 J=0 [“SO(4) parent”] and reversing the sign of J 𝐽 J (“Hund’s partner”), and whether the state can be found in a single-band mean-field (MF) computation neglecting residual interactions and/or when the ferromagnetic (FM) fluctuation corrections discussed in Sec.  V dominate. In the last line, η 𝜂 \eta describes the temperature-dependent strength of admixing of the unitary triplet state.

Pairing	$M_{𝒌+}$	Nodes	SO(4) parent	Hund’s partner	MF/FM
$A^{1_s}$	${\lambda }_{𝒌}{\sigma }_0$	none	$A(1;0,0,0)$	$A^{3_s}(1,0,0)$	✓/✓
$A^{3_s}(1,0,0)$	${\lambda }_{𝒌}{\sigma }_x$	none	$A(1;0,0,0)$	$A^{1_s}$	✓/✗
$A^{3_s}(1,i,0)$	${\lambda }_{𝒌}\left({\sigma }_x+i{\sigma }_y\right)$	$\downarrow$ gapless/none	$A(1;1,0,0)$	$A^{1_s}+A^{3_s}(1,0,0)$	✗/✓
$A^{1_s}+A^{3_s}(1,0,0)$	${\lambda }_{𝒌}\left({\sigma }_0+\eta {\sigma }_x\right)$	none	$A(1;1,0,0)$	$A^{3_s}(1,i,0)$	✗/✓

Table 2. Table 2: Summary of possible pairing states transforming under the complex representation E 𝐸 E of C 3 subscript 𝐶 3 C_{3} . The labeling of the pairing states and their symmetry properties can be found in the main text. The states are ordered by pure singlet, triplet, and admixtures of singlet and triplet. The latter are only expected generically when the SU(2) × − {}_{-}\times SU(2) + symmetry is weakly broken. We use X 𝒌 subscript 𝑋 𝒌 X_{\boldsymbol{k}} and Y 𝒌 subscript 𝑌 𝒌 Y_{\boldsymbol{k}} to denote real-valued continuous functions on the Brillouin zone that transform as k x subscript 𝑘 𝑥 k_{x} and k y subscript 𝑘 𝑦 k_{y} under C 3 subscript 𝐶 3 C_{3} [see, e.g., Eq. ( 30 )]. The temperature-dependent coefficient η 𝜂 \eta describes the admixture of a triplet/singlet pairing at a second transition to a purely singlet/triplet one. Furthermore, a , b ∈ ℝ 𝑎 𝑏 ℝ a,b\in\mathbb{R} vary continuously with system parameters. The minimal number of nodes on any Fermi surface enclosing the Γ Γ \Gamma , K 𝐾 K , or K ′ superscript 𝐾 ′ K^{\prime} point is indicated in the column “Nodes”. As before, two states are referred to as Hund’s partners if they transform into each other under reversing the sign of the Hund’s coupling, see, e.g., Fig.  5 . As singlet and triplet mix for both δ > 0 𝛿 0 \delta>0 and δ < 0 𝛿 0 \delta<0 , there are no Hund’s partners for E 3 s ​ ( g ) superscript 𝐸 subscript 3 𝑠 𝑔 E^{3_{s}}(g) and E 3 s ​ ( h ) superscript 𝐸 subscript 3 𝑠 ℎ E^{3_{s}}(h) ; the corresponding mixed phases, contained in the last two lines of the table, are their own Hund’s partners.

Pairing	$M_{𝒌+}$	Nodes	Hund’s partner	MF/FM
$E^{1_s}(1,0)$	$X_{𝒌}{\sigma }_0$	2 points	$E^{3_s}\left(a\right)$	✗/✓
$E^{1_s}(1,i)$	$\left(X_{𝒌}+i{Y}_{𝒌}\right){\sigma }_0$	0	$E^{3_s}\left(d\right)$	✓/✓
$E^{3_s}\left(a\right)$	$X_{𝒌}{\sigma }_x$	2 points	$E^{1_s}(1,0)$	✗/✗
$E^{3_s}\left(b\right)$	$X_{𝒌}{\sigma }_x+Y_{𝒌}{\sigma }_y$	0	$E^{1_s}(0,i)+E^{3_s}\left(a\right)$	✓/✗
$E^{3_s}\left(c\right)$	$X_{𝒌}\left({\sigma }_x+i{\sigma }_y\right)$	$\downarrow$ gapless/2 points	$E^{1_s}(1,0)+E^{3_s}\left(a\right)$	✗/✓
$E^{3_s}\left(d\right)$	$\left(X_{𝒌}+i{Y}_{𝒌}\right){\sigma }_x$	0	$E^{1_s}(1,i)$	✓/✓
$E^{3_s}\left(e\right)$	$\left(X_{𝒌}+i{Y}_{𝒌}\right)\left({\sigma }_x+i{\sigma }_y\right)$	$\downarrow$ gapless/0	$E^{1_s}(1,i)+E^{3_s}\left(d\right)$	✗/✓
$E^{3_s}\left(f\right)$	$\left(X_{𝒌}+i{Y}_{𝒌}\right){\sigma }_x+\left(X_{𝒌}-i{Y}_{𝒌}\right){\sigma }_y$	2 points	$E^{1_s}(1,-i)+E^{3_s}\left(d\right)$	✗/✗
$E^{3_s}\left(g\right)$	$a\left(X_{𝒌}+i{Y}_{𝒌}\right)\left({\sigma }_x+i{\sigma }_y\right)+b\left(X_{𝒌}-i{Y}_{𝒌}\right){\sigma }_z$	0	—	✗/✗
$E^{3_s}\left(h\right)$	$a\left[\left(X_{𝒌}+i{Y}_{𝒌}\right){\sigma }_x+\left(X_{𝒌}-i{Y}_{𝒌}\right){\sigma }_y\right]+b{Y}_{𝒌}{\sigma }_z$	0	—	✗/✗
$E^{1_s}(0,i)+E^{3_s}\left(a\right)$	$i{Y}_{𝒌}{\sigma }_0+\eta X_{𝒌}{\sigma }_x$	0	$E^{3_s}\left(b\right)$	✓/✗
$E^{1_s}(1,0)+E^{3_s}\left(a\right)$	$X_{𝒌}\left({\sigma }_0+\eta {\sigma }_x\right)$	2 points	$E^{3_s}\left(c\right)$	✗/✓
$E^{1_s}(1,i)+E^{3_s}\left(d\right)$	$\left(X_{𝒌}+i{Y}_{𝒌}\right)\left({\sigma }_0+\eta {\sigma }_x\right)$	0	$E^{3_s}\left(e\right)$	✗/✓
$E^{1_s}(1,-i)+E^{3_s}\left(d\right)$	$\left(X_{𝒌}+i{Y}_{𝒌}\right){\sigma }_0+\eta \left(X_{𝒌}-i{Y}_{𝒌}\right){\sigma }_x$	2 points	$E^{3_s}\left(f\right)$	✗/✗
$E^{3_s}\left(g\right)+E^{1_s}(1,-i)$	$a\left(X_{𝒌}+i{Y}_{𝒌}\right)\left({\sigma }_x+i{\sigma }_y\right)+\left(X_{𝒌}-i{Y}_{𝒌}\right)\left(b{\sigma }_z+\eta {\sigma }_0\right)$	0	$E^{1_s}(1,-i)+E^{3_s}\left(g\right)$	✗/✗
$E^{3_s}\left(h\right)+E^{1_s}(1,0)$	$a\left[\left(X_{𝒌}+i{Y}_{𝒌}\right){\sigma }_x+\left(X_{𝒌}-i{Y}_{𝒌}\right){\sigma }_y\right]+b{Y}_{𝒌}{\sigma }_z+\eta X_{𝒌}{\sigma }_0$	0	$E^{1_s}(1,0)+E^{3_s}\left(h\right)$	✗/✗

Table 3. Table 3: Character table of the point group D 3 subscript 𝐷 3 D_{3} together with the corresponding basis functions and IRs of D 6 subscript 𝐷 6 D_{6} for singlet/triplet pairing.

	$E$	$2C_3$	$3C_{2y}$	Basis functions	IRs of $D_6$
$A_1$	$1$	$1$	$1$	$x^2+y^2$/$y\left(3x^2-y^2\right)$	$A_1$/$B_1$
$A_2$	$1$	$1$	$-1$	$z$/$x\left(x^2-3y^2\right)$	$A_2$/$B_2$
$E$	$2$	$-1$	$0$	$(2xy,x^2-y^2)$/$(x,y)$	$E_2$/$E_1$

Table 4. Table 4: Summary of the different intervalley pairing states classified by the IRs of the point group D 3 subscript 𝐷 3 D_{3} . The notation closely parallels that of Table  2 . Here, we use λ 𝒌 1 subscript superscript 𝜆 1 𝒌 \lambda^{1}_{\boldsymbol{k}} and λ 𝒌 2 subscript superscript 𝜆 2 𝒌 \lambda^{2}_{\boldsymbol{k}} to denote continuous functions on the Brillouin zone that are even and odd under ( k x , k y ) → ( k x , − k y ) → subscript 𝑘 𝑥 subscript 𝑘 𝑦 subscript 𝑘 𝑥 subscript 𝑘 𝑦 (k_{x},k_{y})\rightarrow(k_{x},-k_{y}) , respectively, and are both invariant under C 3 subscript 𝐶 3 C_{3} , λ 𝒌 j = λ C 3 ​ 𝒌 j subscript superscript 𝜆 𝑗 𝒌 subscript superscript 𝜆 𝑗 subscript 𝐶 3 𝒌 \lambda^{j}_{\boldsymbol{k}}=\lambda^{j}_{C_{3}\boldsymbol{k}} . Furthermore, X 𝒌 φ subscript superscript 𝑋 𝜑 𝒌 X^{\varphi}_{\boldsymbol{k}} and Y 𝒌 φ subscript superscript 𝑌 𝜑 𝒌 Y^{\varphi}_{\boldsymbol{k}} are rotated basis functions defined in Eq. ( 56 ); for instance, a possible choice for twisted bilayer graphene with Brillouin zone in Fig.  1 (b) is given by ( X 𝒌 , Y 𝒌 ) T = R ( π + φ ) / 2 ​ ( X 𝒌 ( 1 ) , Y 𝒌 ( 1 ) ) T superscript subscript 𝑋 𝒌 subscript 𝑌 𝒌 𝑇 subscript 𝑅 𝜋 𝜑 2 superscript subscript superscript 𝑋 1 𝒌 subscript superscript 𝑌 1 𝒌 𝑇 (X_{\boldsymbol{k}},Y_{\boldsymbol{k}})^{T}=R_{(\pi+\varphi)/2}(X^{(1)}_{\boldsymbol{k}},Y^{(1)}_{\boldsymbol{k}})^{T} with X 𝒌 ( 1 ) subscript superscript 𝑋 1 𝒌 X^{(1)}_{\boldsymbol{k}} , Y 𝒌 ( 1 ) subscript superscript 𝑌 1 𝒌 Y^{(1)}_{\boldsymbol{k}} in Eq. ( 30a ). To keep the notation short, each line with reference to φ 1 subscript 𝜑 1 \varphi_{1} or φ 2 subscript 𝜑 2 \varphi_{2} corresponds to two distinct states with φ 1 = 0 , π / 3 subscript 𝜑 1 0 𝜋 3 \varphi_{1}=0,\pi/3 and φ 2 = 0 , π / 2 subscript 𝜑 2 0 𝜋 2 \varphi_{2}=0,\pi/2 . The indicated number of nodal points refers to a Fermi surface enclosing the Γ Γ \Gamma point.

Pairing	$M_{𝒌+}$	Nodes around $\Gamma$	Hund’s partner	MF/FM
$A_1^{1_s}$	${\lambda }_{𝒌}^1{\sigma }_0$	none	$A_2^{3_s}(1,0,0)$	✓/✓
$A_2^{1_s}$	${\lambda }_{𝒌}^2{\sigma }_0$	6 points	$A_1^{3_s}(1,0,0)$	✓/✓
$A_1^{3_s}(1,0,0)$	${\lambda }_{𝒌}^2{\sigma }_x$	6 points	$A_2^{1_s}$	✓/✗
$A_2^{3_s}(1,0,0)$	${\lambda }_{𝒌}^1{\sigma }_x$	none	$A_1^{1_s}$	✓/✗
$A_1^{3_s}(1,i,0)$	${\lambda }_{𝒌}^2\left({\sigma }_x+i{\sigma }_y\right)$	$\downarrow$ gapless/6 points	$A_2^{1_s}+A_1^{3_s}(1,0,0)$	✗/✓
$A_2^{3_s}(1,i,0)$	${\lambda }_{𝒌}^1\left({\sigma }_x+i{\sigma }_y\right)$	$\downarrow$ gapless/none	$A_1^{1_s}+A_2^{3_s}(1,0,0)$	✗/✓
$A_1^{1_s}+A_2^{3_s}(1,0,0)$	${\lambda }_{𝒌}^1\left({\sigma }_0+\eta \left({\sigma }_x+i{\sigma }_y\right)\right)$	none	$A_2^{3_s}(1,i,0)$	✗/✓
$A_2^{1_s}+A_1^{3_s}(1,0,0)$	${\lambda }_{𝒌}^2\left({\sigma }_0+\eta \left({\sigma }_x+i{\sigma }_y\right)\right)$	6 points	$A_1^{3_s}(1,i,0)$	✗/✓
$E^{1_s}{(1,0)}_{{\phi }_1}$	$X_{𝒌}^{{\phi }_1}{\sigma }_0$	2 points	$E^{3_s}{\left(a\right)}_{{\phi }_1}$	✗/✓
$E^{1_s}(1,i)$	$\left(X_{𝒌}^0+i{Y}_{𝒌}^0\right){\sigma }_0$	0	$E^{3_s}\left(d\right)$	✓/✓
$E^{3_s}{\left(a\right)}_{{\phi }_1}$	$X_{𝒌}^{{\phi }_1}{\sigma }_x$	2 points	$E^{1_s}{(1,0)}_{{\phi }_1}$	✗/✗
$E^{3_s}\left(b\right)$	$X_{𝒌}^0{\sigma }_x+Y_{𝒌}^0{\sigma }_y$	0	${\left(E^{1_s}(0,i)+E^{3_s}\left(a\right)\right)}_{{\phi }_1}$	✓/✗
$E^{3_s}{\left(c\right)}_{{\phi }_1}$	$X_{𝒌}^{{\phi }_1}\left({\sigma }_x+i{\sigma }_y\right)$	$\downarrow$ gapless/2 points	${\left(E^{1_s}(1,0)+E^{3_s}\left(a\right)\right)}_{{\phi }_1}$	✗/✓
$E^{3_s}\left(d\right)$	$\left(X_{𝒌}^0+i{Y}_{𝒌}^0\right){\sigma }_x$	0	$E^{1_s}(1,i)$	✓/✗
$E^{3_s}\left(e\right)$	$\left(X_{𝒌}^0+i{Y}_{𝒌}^0\right)\left({\sigma }_x+i{\sigma }_y\right)$	$\downarrow$ gapless/0	$E^{1_s}(1,i)+E^{3_s}\left(d\right)$	✗/✓
$E^{3_s}{\left(f\right)}_{{\phi }_2}$	$\left(X_{𝒌}^{{\phi }_2}+i{Y}_{𝒌}^{{\phi }_2}\right){\sigma }_x+\left(X_{𝒌}^{{\phi }_2}-i{Y}_{𝒌}^{{\phi }_2}\right){\sigma }_y$	2 points	${\left(E^{1_s}(1,-i)+E^{3_s}\left(d\right)\right)}_{{\phi }_2}$	✗/✗
$E^{3_s}\left(g\right)$	$a\left(X_{𝒌}^0+i{Y}_{𝒌}^0\right)\left({\sigma }_x+i{\sigma }_y\right)+b\left(X_{𝒌}^0-i{Y}_{𝒌}^0\right){\sigma }_z$	0	—	✗/✗
$E^{3_s}{\left(h\right)}_{{\phi }_1}$	$a\left[\left(X_{𝒌}^{{\phi }_1}+i{Y}_{𝒌}^{{\phi }_1}\right){\sigma }_x+\left(X_{𝒌}^{{\phi }_1}-i{Y}_{𝒌}^{{\phi }_1}\right){\sigma }_y\right]+b{Y}_{𝒌}^{{\phi }_1}{\sigma }_z$	0	—	✗/✗
${\left(E^{1_s}(0,i)+E^{3_s}\left(a\right)\right)}_{{\phi }_1}$	$i{Y}_{𝒌}^{{\phi }_1}{\sigma }_0+\eta X_{𝒌}^{{\phi }_1}{\sigma }_x$	0	$E^{3_s}\left(b\right)$	✓/✗
${\left(E^{1_s}(1,0)+E^{3_s}\left(a\right)\right)}_{{\phi }_1}$	$X_{𝒌}^{{\phi }_1}\left({\sigma }_0+\eta {\sigma }_x\right)$	2 points	$E^{3_s}{\left(c\right)}_{{\phi }_1}$	✗/✓
$E^{1_s}(1,i)+E^{3_s}\left(d\right)$	$\left(X_{𝒌}^0+i{Y}_{𝒌}^0\right)\left({\sigma }_0+\eta {\sigma }_x\right)$	0	$E^{3_s}\left(e\right)$	✗/✓
${\left(E^{1_s}(1,-i)+E^{3_s}\left(d\right)\right)}_{{\phi }_2}$	$\left(X_{𝒌}^{{\phi }_2}+i{Y}_{𝒌}^{{\phi }_2}\right){\sigma }_0+\eta \left(X_{𝒌}^{{\phi }_2}-i{Y}_{𝒌}^{{\phi }_2}\right){\sigma }_x$	2 points	$E^{3_s}{\left(f\right)}_{{\phi }_2}$	✗/✗
$E^{3_s}\left(g\right)+E^{1_s}(1,-i)$	$a\left(X_{𝒌}^0+i{Y}_{𝒌}^0\right)\left({\sigma }_x+i{\sigma }_y\right)+\left(X_{𝒌}^0-i{Y}_{𝒌}^0\right)\left(b{\sigma }_z+\eta {\sigma }_0\right)$	0	$E^{1_s}(1,-i)+E^{3_s}\left(g\right)$	✗/✗
${\left(E^{1_s}(1,0)+E^{3_s}\left(h\right)\right)}_{{\phi }_1}$	$a\left[\left(X_{𝒌}^{{\phi }_1}+i{Y}_{𝒌}^{{\phi }_1}\right){\sigma }_x+\left(X_{𝒌}^{{\phi }_1}-i{Y}_{𝒌}^{{\phi }_1}\right){\sigma }_y\right]+b{Y}_{𝒌}^{{\phi }_1}{\sigma }_z+\eta X_{𝒌}^{{\phi }_1}{\sigma }_0$	0	${\left(E^{3_s}\left(h\right)+E^{1_s}(1,0)\right)}_{{\phi }_1}$	✗/✗

Equations

$${\mathcal{H}}_{\text{SC}}=\sum _{\mathbf{k}}{c}_{\mathbf{k}\sigma v}^{†}{\left({\Delta }_{\mathbf{k}}i{\sigma }_y{\tau }_x\right)}_{\sigma v,{\sigma }^{\prime }{v}^{\prime }}{c}_{-\mathbf{k}{\sigma }^{\prime }{v}^{\prime }}^{†}+\text{H.c.}.$$

$${\mathcal{G}}_1=C_3\times \text{SU}(2{)}_{+}\times \text{SU}(2{)}_{-}\times U(1{)}_v,$$

$${\mathcal{G}}_2=C_3\times \text{SU}(2{)}_s\times U(1{)}_v,$$

$$C_3:$$

$$\text{SU}(2{)}_s:$$

$$\text{SU}(2{)}_{\pm }:$$

$$\text{U}(1{)}_v:$$

$$\Theta c_{\mathbf{k}}{\Theta }^{†}=i{\sigma }_y{\tau }_x{c}_{-\mathbf{k}}.$$

$${\Delta }_{\mathbf{k}}=\sum _n\sum _{\mu =1}^{d_n}{\eta }_{\mu }^n{\chi }_{\mu }^n(\mathbf{k}),{\eta }_{\mu }^n\in \mathbb{C},$$

$${\Delta }_{\mathbf{k}}=\sum _{n_{C_3},n_v,n_s}\sum _{{\mu }_1=1}^{d_{n_{C_3}}}\sum _{{\mu }_2=1}^{d_{n_v}}\sum _{{\mu }_3=1}^{d_{n_s}}{\eta }_{{\mu }_1{\mu }_2{\mu }_3}^n{\chi }_{{\mu }_1}^{n_{C_3}}(\mathbf{k}){\chi }_{{\mu }_2}^{n_v}{\chi }_{{\mu }_3}^{n_s}.$$

$${\left({\Delta }_{\mathbf{k}}\right)}_{\sigma v,{\sigma }^{\prime }{v}^{\prime }}={\delta }_{v,v^{\prime }}{\chi }^A(\mathbf{k},v)\sum _{\mu =1}^{d_{n_s}}{\eta }_{\mu }^{n_s}{\left({\chi }_{\mu }^{n_s}(v)\right)}_{\sigma {\sigma }^{\prime }},$$

$${\mathcal{H}}_{\text{SC}}$$

$$M_{\mathbf{k}v}$$

$$M_{\mathbf{k}v}={\sigma }_y{M}_{-\mathbf{k}\bar{v}}^T{\sigma }_y.$$

$$\Theta :$$

$$\text{SU}(2{)}_{+}\times \text{SU}(2{)}_{-}:$$

$$\mathcal{F}$$

$$+\frac{b_2}{2}\text{tr}\left[{\Delta }_{+}^{†}{\Delta }_{+}^{}{\Delta }_{+}^{†}{\Delta }_{+}^{}\right]+\frac{b_3}{4}{\left|\text{tr}\left[{\sigma }_y{\Delta }_{+}^{}{\sigma }_y{\Delta }_{+}^T\right]\right|}^2.$$

$${\mathcal{H}}_{\text{SC}}=\sum _{\mathbf{k},v}{\lambda }_{\mathbf{k}v}^s{c}_{\mathbf{k}\sigma v}^{†}{\left(i{\sigma }_y\right)}_{\sigma ,{\sigma }^{\prime }}{c}_{-\mathbf{k}{\sigma }^{\prime }\bar{v}}^{†}+\text{H.c.},$$

$$\mathcal{F}\sim a(T){\mathbf{d}}^{†}\mathbf{d}+b_1^t{\left({\mathbf{d}}^{†}\mathbf{d}\right)}^2+b_2^t|{\mathbf{d}}^{\ast }\times \mathbf{d}{|}^2.$$

$${\mathcal{H}}_{\text{SC}}=\sum _{\mathbf{k},v}{\lambda }_{\mathbf{k}v}^t{c}_{\mathbf{k}\sigma v}^{†}{\left({\sigma }_{x}i{\sigma }_y\right)}_{\sigma ,{\sigma }^{\prime }}{c}_{-\mathbf{k}{\sigma }^{\prime }\bar{v}}^{†}+\text{H.c.},$$

$${\mathcal{H}}_{\text{SC}}=\sum _{\mathbf{k},v}{\lambda }_{\mathbf{k}v}^t{c}_{\mathbf{k}\sigma v}^{†}{\left(({\sigma }_x+i{\sigma }_y)i{\sigma }_y\right)}_{\sigma ,{\sigma }^{\prime }}{c}_{-\mathbf{k}{\sigma }^{\prime }\bar{v}}^{†}+\text{H.c.},$$

$$b_2^t=T\sum _{{\omega }_n}\int \frac{{\mathrm{d}}^2\mathbf{k}}{(2\pi {)}^2}\frac{|{\lambda }_{\mathbf{k}}^t{|}^4}{({\omega }_n^2+{\xi }_{\mathbf{k}+}^2{)}^2}.$$

$$\mathcal{F}\sim a(T)\left(|{\Delta }^s{|}^2+{\mathbf{d}}^{†}\mathbf{d}\right)+\delta a(T)\left(|{\Delta }^s{|}^2-{\mathbf{d}}^{†}\mathbf{d}\right),$$

$$\frac{\Delta T_c}{T_{c,0}}\sim \frac{|J|}{V^2\nu }.$$

$$\mathcal{F}$$

$$+(b_1+b_2)|{\Delta }^s{|}^4+(b_1+b_2){\left({\mathbf{d}}^{†}\mathbf{d}\right)}^2+b_2{\left|{\mathbf{d}}^{\ast }\times \mathbf{d}\right|}^2$$

$$+2(b_1+2b_2)|{\Delta }^s{|}^2{\mathbf{d}}^{†}\mathbf{d}+2b_2\text{Re}\left[{\left({\Delta }^s\right)}^2{\mathbf{d}}^{†}{\mathbf{d}}^{\ast }\right],$$

$$A_{m_v=0}(1;0,0,0)$$

$$A_{m_v=0}(1;1,0,0)$$

$$2c|{\Delta }^s{|}^2|\mathbf{d}(T){|}^2=(\delta a-a(T))|{\Delta }^s{|}^2,$$

$$T_c^{-}=T_{c,0}\left(1+\frac{c-|b_2|}{|b_2|}\delta \right),\delta \equiv \frac{\delta a}{\alpha T_{c,0}}\simeq \frac{J}{V^2\nu }.$$

$${\mathcal{F}}_M$$

$$+2\delta c_1{\mathbf{M}}_Z\cdot \text{Im}\left({\mathbf{d}}^{\ast }{\Delta }^s\right)+i{c}_2{\mathbf{M}}_Z\cdot {\mathbf{d}}^{\ast }\times \mathbf{d}$$

$$+(c_3{\mathbf{M}}_Z^2+c_4{\mathbf{M}}_O^2)\left(|{\Delta }^s{|}^2+{\mathbf{d}}^{†}\mathbf{d}\right)$$

$$+(\delta c_5{\mathbf{M}}_Z^2+\delta c_6{\mathbf{M}}_O^2)\left(|{\Delta }^s{|}^2-{\mathbf{d}}^{†}\mathbf{d}\right).$$

Full text

I Introduction

Experiments on twisted bilayers of graphene have recently revealed interaction-induced insulating phases and superconductivity when the relative angle between the layers is fine-tuned to yield almost flat moiré bands, which enhances the impact of electronic correlations Cao et al. (2018a, b); Yankowitz et al. (2019); Lu et al. (2019). Due to the strong-coupling nature of the problem, which is corroborated by tunneling spectroscopy measurements Kerelsky et al. (2019); Choi et al. (2019); Jiang et al. (2019); Xie et al. (2019); Scheurer (2019), the form and mechanism of the insulating and superconducting phases are still under debate, despite considerable theoretical effort Baskaran (2018); Wu et al. (2018a); Guo et al. (2018); Koshino et al. (2018); Yuan and Fu (2018); Po et al. (2018); Lian et al. (2019); Zou et al. (2018); Dodaro et al. (2018); Xu and Balents (2018); Thomson et al. (2018); Fidrysiak et al. (2018); Kang and Vafek (2018); Laksono et al. (2018); Kang and Vafek (2019); Seo et al. (2019); Isobe et al. (2018); Liu et al. (2018); Sherkunov and Betouras (2018); Su and Lin (2018); Venderbos and Fernandes (2018); Peltonen et al. (2018); Liu et al. (2019a); Kennes et al. (2018); Choi and Choi (2018); Wu et al. (2018b); Lin and Nandkishore (2019); Huang et al. (2019); González and Stauber (2019); Chen et al. (2019a, 2020a); Roy and Juričić (2019); Wu and Das Sarma (2019); Kozii et al. (2019); You and Vishwanath (2019); Ray et al. (2019); Alidoust et al. (2019). Another graphene-based moiré system that displays both superconducting and correlated insulating behavior is ABC-stacked trilayer graphene on hexagonal boron nitride Chen et al. (2019b, c). In this case, the moiré pattern results from the difference in lattice constants, and it can be controlled by application of a vertical electric field Chittari et al. (2019); Zhang and Senthil (2019).

The most recent member of the family of strongly correlated graphene superlattice systems is twisted double-bilayer graphene Shen et al. (2019); Liu et al. (2019b); Cao et al. (2019), where two individually aligned AB-stacked graphene bilayers are twisted with respect to one another. As theoretical calculations show Zhang et al. (2019); Chebrolu et al. (2019); Choi and Choi (2019); Lee et al. (2019); Koshino (2019); Liu et al. (2019); Haddadi et al. (2020), flat electronic bands can be realized by tuning the twist angle and a vertical electric field. Similar to the abovementioned graphene moiré systems, both correlated insulating Shen et al. (2019); Liu et al. (2019b); Cao et al. (2019) and superconducting Shen et al. (2019); Liu et al. (2019b) phases are observed in experiment. However, in stark contrast to twisted bilayer and trilayer graphene, the superconducting transition temperature is found to increase linearly with a weak in-plane magnetic field Liu et al. (2019b), which is a strong indication of triplet pairing Lee et al. (2019); Wu and Das Sarma (2019). Furthermore, the gap of the correlated insulating phase is seen to increase with an applied magnetic field, indicating ferromagnetic order Shen et al. (2019); Liu et al. (2019b); Cao et al. (2019). There are also clear experimental indications of ferromagnetism in twisted bilayer Sharpe et al. (2019); Lu et al. (2019); Zondiner et al. (2019) and ABC trilayer graphene Chen et al. (2020b).

In this paper, we study the possible pairing states in graphene moiré superlattices. Motivated by the recent experimental signatures of triplet pairing, we pay special attention to the triplet channel, and possible mixed singlet and triplet phases. While the weak spin-orbit coupling in graphene seems to disfavor the latter class of phases, projections of the Coulomb interaction on the relevant moiré bands evince that the interaction terms that couple the spin degrees of freedom of the two valleys, $v=\pm$, of the system are much weaker than other interaction terms that do not Koshino et al. (2018); Zhang and Senthil (2019); Lee et al. (2019). Together with the nearly valley-diagonal band structure, this indicates that the system is approximately invariant under independent spin rotations in the two valleys. As has been pointed out before Xu and Balents (2018); You and Vishwanath (2019), the associated SU(2)+ $\times$ SU(2)- symmetry renders the singlet and triplet pairing channels degenerate. This paper will address the questions: (i) under which conditions can singlet and triplet mix when the SU(2)+ $\times$ SU(2)- symmetry is only weakly broken, and (ii) which triplet state transforms into which singlet upon reversing the sign of the symmetry-breaking interactions? In this way, we map out all possible phase diagrams close to the SU(2)+ $\times$ SU(2)--invariant limit.

In light of the narrow bandwidth and strong correlations of the graphene-based moiré superlattices we are interested in, our analysis will begin with a comprehensive study of exact constraints resulting from symmetry. The symmetry-based classification will then be supplemented with energetics, by studying which of the pairing states can be realized in the weak-coupling limit and what changes in the presence of additional fluctuation corrections. As it has the smallest set of symmetries, we will begin our classification with twisted double-bilayer graphene: while the lattice is invariant under threefold rotation, $C_{3}$, perpendicular to the graphene sheets, and under a twofold in-plane rotation [see Fig. 1(a)], the latter is broken due to the vertical electric field that is applied to tune the band structure and to induce superconductivity. It seems currently unclear whether the superconducting state coexists with the likely ferromagnetic correlated insulator and whether, at least in part of the phase diagram, there is a thermal transition directly from the (paramagnetic) normal metal to superconductivity without any ferromagnetic order. For this reason, we will analyze two scenarios separately: (I) there is no ferromagnetic order around the critical temperature, $T_{c}$, of superconductivity, and (II) there is ferromagnetic order already at $T>T_{c}$ that coexists microscopically with superconductivity for $T<T_{c}$ (or at the minimum, the associated ferromagnetic moments couple significantly to the superconducting order parameter). We will begin with the analysis of the superconducting states transforming under the IRs of the point group $C_{3}$ assuming time-reversal symmetry in the high-temperature phase—this is relevant for case (I) above. In order to capture scenario (II), we will later add the coupling to the time-reversal-symmetry breaking magnetic moments and examine how it affects the superconducting transition. This allows us to determine which of the pairing states are compatible with the linear increase of the critical temperature with small magnetic field, $B$, and to describe the possible phase diagrams in the temperature-$B$ plane.

We also generalize our discussion to twisted bilayer graphene and ABC trilayer graphene on hexagonal boron nitride. Here, we have to take into account an additional twofold rotation symmetry, $C_{2}$, perpendicular to the plane of the system and an in-plane rotation symmetry $C_{2y}$ along the $y$-axis; these symmetries are either realized as exact microscopic symmetries of the lattice or as approximate emergent symmetries Po et al. (2018); Zou et al. (2018); Zhang and Senthil (2019) of those systems, see Fig. 1(b) and (c).

I.1 Brief summary of the main results

Due to the length of the paper, here, we provide a very concise summary of the key results of this work for the convenience of the reader:

1. (1)

We analyze the consequences of the enhanced SU(2)+ $\times$ SU(2)- spin symmetry, taking into account the possibility of several consecutive superconducting transitions with their difference in transition temperatures vanishing in the limit where SU(2)+ $\times$ SU(2)- becomes exact. The resulting complete sets of possible phases for the relevant symmetry groups $C_{3}$ and $D_{3}$ (or, equivalently, $D_{6}$, see Sec. VI.1) are summarized in Tables 1, 2, and 4. 2. (2)

As follows from these tables, all point groups and all of their irreducible representations (IRs) allow for singlet-triplet admixed phases in the absence of any spin-orbit coupling. As opposed to the conventional mechanism of singlet-triplet admixture, which is based on a reduced symmetry Gor’kov and Rashba (2001), here, it results from (the proximity to) an enhanced spin symmetry. 3. (3)

To supplement these purely symmetry-based considerations with energetics, we analyze which of those states can be realized in single-band mean-field theory, i.e., whether there exists a form of the effective electron-electron interaction that can stabilize the superconducting state when treated within the mean-field approximation; the result is indicated in the last column in Tables 1, 2, and 4. This identifies the most important pairing states from a weak-coupling perspective. The presence of any of the remaining pairing phases—as might eventually be established in future experiments—must result from the strong-coupling and/or interband nature of superconductivity. 4. (4)

We also study corrections to mean-field theory coming from ferromagnetic fluctuations, within a simplified phenomenological approach in Sec. V that is justified microscopically in Appendix B.1. We analyze two limits. First, we consider the case of weak fluctuations in order to lift the residual degeneracies within mean-field theory. We find that out of the two possible phase diagrams for the IR $E$ close to mean-field theory, shown in Fig. 5, the one in part (b) [part (a)] is favored for spin (orbital) ferromagnetic fluctuations. This reveals that a nematic mixed singlet-triplet phase is a natural candidate pairing phase in graphene moiré superlattices. Second, we analyze which pairing states are favored in the case where the fluctuation corrections dominate over the mean-field contributions (see last column in the tables mentioned above). 5. (5)

We study the coupling of the superconducting states to the magnetic field, $B$, and examine which states can give rise to a linear increase of the critical temperature for small $B$: if SU(2)+ $\times$ SU(2)- is broken significantly, triplet pairing has to dominate for $B=0$ and there are only three possible triplet states as leading instabilities for $B\neq 0$. In the case where SU(2)+ $\times$ SU(2)- is an approximate symmetry, even singlet pairing at $B=0$ can yield a linear increase. For instance, the possible phase diagrams in the presence of a magnetic field for pairing in the trivial IR $A$ of $C_{3}$ are shown in Fig. 3.

I.2 Relation to other works

Let us briefly comment on the relation of this article to other works in the literature. While our classification also contains the pure singlet states, which have been subject to intense scrutiny in twisted bilayer graphene, we are mainly interested in elucidating the consequences of the enhanced SU(2)+ $\times$ SU(2)- spin symmetry with respect to subsequent transitions and the associated nontrivial interplay of singlet and triplet pairing.

In the context of twisted double-bilayer graphene, where, recently, signs of triplet pairing have been discovered, Ref. Wu et al., 2019 mainly focuses on the correlated insulating phase in this system, whereas Ref. Lee et al., 2019 also discusses pairing. We extend the work of Ref. Lee et al., 2019 by allowing for momentum-dependent pairing states, contrasting weakly and significantly broken SU(2)+ $\times$ SU(2)- symmetry, investigating admixed singlet and triplet phases, analyzing fluctuation corrections to mean-field theory, and mapping out the phase diagram in the presence of a magnetic field. In a follow-up work Samajdar and Scheurer (2020), we will complement the analysis of this paper by a microscopic energetic study specifically for twisted double-bilayer graphene. In Ref. Samajdar and Scheurer, 2020, we discuss which IR is expected to be favored, the form of the associated basis functions, and the impact of disorder on superconductivity.

I.3 Structure of the paper

This paper is organized as follows: as described above, we start with twisted double-bilayer graphene. In Sec. II, we introduce the model and the action of the relevant symmetries. We first discuss pairing in the trivial IR of the point group of the system in Sec. III and then generalize to the complex IR $E$ in Sec. IV. Section V demonstrates how strong fluctuations can yield significant corrections to mean-field theory. We extend our analysis to twisted bilayer graphene and ABC trilayer graphene in Sec. VI, and explore the consequences of the additional microscopic and emergent symmetries relevant to those systems. A discussion of our results can be found in Sec. VII.

II Model and symmetries

We first focus on the (nearly flat) conduction band of twisted double-bilayer graphene which appears to host the superconducting phase observed experimentally Shen et al. (2019); Liu et al. (2019b), and later discuss the modifications for the related moiré systems, bilayer and trilayer graphene. Owing to the presence of a gap to other bands in the relevant parameter regime Choi and Choi (2019); Liu et al. (2019); Koshino (2019); Lee et al. (2019), it is reasonable to describe the superconducting instability in a single-band picture. We stress, however, that many of our conclusions are symmetry-based and thus, also apply when several bands are taken into account. Exceptions are provided by the energetic mean-field and fluctuation considerations, where we will specifically comment on the consequences of interband effects that might be present in these systems Xie et al. (2019).

Denoting the corresponding electronic creation and annihilation operators by $c_{\boldsymbol{k}\sigma v}$, where $\boldsymbol{k}$ is crystal momentum, $\sigma$ spin, and $v=\pm$ represents the valleys, the general pairing term can be written as

[TABLE]

Here and in the following, $\sigma_{j}$ and $\tau_{j}$ are Pauli matrices in spin and valley space, respectively, and the $4\times 4$ matrix $\Delta_{\boldsymbol{k}}$ is the superconducting order parameter. In Eq. (1), we have already made the assumptions that only Cooper pairs with zero net momentum form and that superconductivity preserves translational symmetry. Due to the proximity of superconductivity to ferromagnetic order Shen et al. (2019); Liu et al. (2019b), relaxing this assumption could be interesting, but we leave this for future work. Consequently, we need not consider IRs of the full space group but rather, can concentrate on the point group $\mathcal{G}$ of the system and time-reversal $\Theta$.

In this regard, we study two different point groups: an approximate point group,

[TABLE]

where $C_{3}$ is the crystalline point group, SU(2)± is spin rotation in valley $v=\pm$, and $U(1)_{v}$ corresponds to valley charge conservation. As argued in Ref. Lee et al., 2019, the intervalley “Hund’s” coupling $J$ is much smaller than the intravalley-density interaction $V$. In combination with the fact that the noninteracting band structure only has very small valley mixing, the system is invariant under Eq. (2) to a good approximation. In the presence of a finite Hund’s coupling, Eq. (2) is reduced to

[TABLE]

where SU(2)s is global spin rotation. To define these symmetries more precisely, we specify their representation on the electronic field operators:

[TABLE]

with $P_{\pm}=(\tau_{0}\pm\tau_{z})/2$ being the valley projection operators. Furthermore, time-reversal is represented by the antiunitary operator $\Theta$ with

[TABLE]

To classify superconductivity, we proceed as usual Sigrist and Ueda (1991) and express $\Delta_{\boldsymbol{k}}$ in Eq. (1) in terms of the IRs $n$ (with dimension $d_{n}$) of the point group as

[TABLE]

where $\chi_{\mu}^{n}(\boldsymbol{k})$, $\mu=1,\dots,d_{n}$, are partner functions transforming under the IR $n$. Within the minimal description of pairing in Eq. (1), which only involves one band per valley, $\chi_{\mu}^{n}(\boldsymbol{k})\in\mathbb{C}^{4\times 4}$ are matrices in spin and valley space.

In our case, the point group has the form $\mathcal{G}_{j}=C_{3}\times U(1)_{v}\times\mathcal{G}^{s}_{j}$ with $\mathcal{G}^{s}_{1}=\text{SU}(2)_{+}\times\text{SU}(2)_{-}\simeq\text{SO}(4)$ and $\mathcal{G}^{s}_{2}=\text{SU}(2)_{s}$. As a consequence, the IRs of $\mathcal{G}_{j}$ have the form $n=n_{C_{3}}\times n_{v}\times n_{s}$ where $n_{C_{3}}$, $n_{v}$, and $n_{s}$ are IRs of $C_{3}$, $U(1)_{v}$, and $\mathcal{G}_{j}^{s}$, respectively. We can thus rewrite Eq. (6) more explicitly as

[TABLE]

In order to classify superconducting states, we need to consider the different IRs of $C_{3}$, $U(1)_{v}$, and $\mathcal{G}_{j}^{s}$.

Let us begin our discussion of IRs with $U(1)_{v}$. While it has, in general, countably infinite IRs (one-dimensional and with character $e^{im_{v}\varphi}$, $m_{v}\in\mathbb{Z}$), only three are relevant here as all representations with $|m_{v}|>1$ cannot be realized with only two valleys. First, there is the trivial representation, $m_{v}=0$, with $\chi^{m_{v}=0}=a\tau_{0}+b\tau_{z}$ with a priori unknown $a,b$. Recalling the extra factor of $\tau_{x}$ in Eq. (1), this translates to purely intervalley pairing. Secondly, the pair of complex conjugate representations with $m_{v}=\pm 1$ has to be considered. Note that due to time-reversal symmetry, the complex representations cannot be discussed separately. Here, the basis functions read as $\chi^{m_{v}=\pm 1}=\tau_{x}\pm i\tau_{y}$; as such, this corresponds to purely intravalley pairing.

We thus see that $U(1)_{v}$ prohibits the mixing of inter- and intravalley pairing. As time-reversal (5) interchanges the valleys along with sending $\boldsymbol{k}\rightarrow-\boldsymbol{k}$ and we assume zero-momentum Cooper pairs, we will restrict our discussion to intervalley pairing, i.e., $m_{v}=0$ for the rest of the paper.

As is well known Dresselhaus et al. (2008), $C_{3}$ has the following IRs, both of which are one-dimensional: the trivial one, $A$, and the complex representation $E$ (and its complex conjugate partner). We analyze each of these IRs in Secs. III and IV, and in both cases, discuss the differences between $\mathcal{G}^{s}_{1}$ and $\mathcal{G}^{s}_{2}$; we will also see how the states “connect” once $\mathcal{G}^{s}_{1}$ is weakly broken to $\mathcal{G}^{s}_{2}$ due to a small but finite value of the Hund’s coupling.

III Trivial representation of the crystalline point group

For simplicity, we begin with the trivial representation $A$ of $C_{3}$, which is real and one-dimensional. In fact, the following discussion will not be modified as long as the IR is real and one-dimensional and there is no crystalline symmetry relating the two valleys. Interestingly, the last assumption is violated in twisted bilayer graphene and trilayer graphene on boron nitride; see Sec. VI for a detailed discussion of the associated modifications.

As already mentioned, we consider only intervalley pairing which corresponds to a real and one-dimensional IR as well. This means that the order parameter in Eq. (7) has the form

[TABLE]

where $\chi^{A}(\boldsymbol{k},v)$ is invariant under $\boldsymbol{k}\rightarrow g\boldsymbol{k}$ for all generators $g$ of the crystalline point group (here, we only have $g=C_{3}$).

III.1 Limit of exact SU(2)${}_{+}\times\,$SU(2)- symmetry

To proceed further, we have to inspect the scenarios for both $\mathcal{G}_{1}^{s}$ and $\mathcal{G}_{2}^{s}$. We start with the former, i.e., we assume that the Hund’s coupling is zero. Inserting Eq. (8) in the general pairing Hamiltonian (1), we obtain a pairing term of the form

[TABLE]

where $\bar{v}=\mp$ for $v=\pm$, and $M_{\boldsymbol{k}v}$ as well as $\Delta_{v}=\sum_{\mu}\eta_{\mu}^{n_{s}}\chi_{\mu}^{n_{s}}(v)$ are matrices in spin space. Fermi-Dirac statistics implies

[TABLE]

Rewriting pairing in terms of singlet and triplet as $M_{\boldsymbol{k}v}=\sigma_{0}\Delta^{s}_{\boldsymbol{k}v}+\boldsymbol{\sigma}\cdot\boldsymbol{d}_{\boldsymbol{k}v}$, Eq. (10) is equivalent to $\Delta^{s}_{\boldsymbol{k}v}=\Delta^{s}_{-\boldsymbol{k}\bar{v}}$ and $\boldsymbol{d}_{\boldsymbol{k}v}=-\boldsymbol{d}_{-\boldsymbol{k}\bar{v}}$, as expected.

We now study the stable superconducting phases in this channel by writing down the most general Ginzburg-Landau expansion constrained by the symmetries

[TABLE]

Due to the constraint (10) stemming from Fermi-Dirac statistics, we express the free energy in terms of one valley only (say $v=+$) as $\mathcal{F}=\mathcal{F}[M_{\boldsymbol{k}+}=\chi^{A}(\boldsymbol{k},+)\Delta_{+}]$, and the pairing in the other valley just follows from Eq. (10). The most general free energy to quartic order in $\Delta_{+}$, invariant under Eq. (11) and $\Delta_{v}\rightarrow e^{i\varphi}\Delta_{v}$, reads as

[TABLE]

Note that $|\text{tr}[\sigma_{y}\Delta_{+}^{\phantom{\dagger}}\sigma_{y}\Delta_{+}^{T}]|^{2}/2=\text{tr}[\Delta_{+}^{\phantom{\dagger}}\sigma_{y}\Delta_{+}^{T}\Delta_{+}^{*}\sigma_{y}\Delta_{+}^{\dagger}]$, so the latter is not an independent term to consider. It further holds that $|\text{tr}[\sigma_{y}\Delta_{+}^{\phantom{\dagger}}\sigma_{y}\Delta_{+}^{T}]|^{2}/2=(\text{tr}[\Delta_{+}^{\dagger}\Delta_{+}^{\phantom{\dagger}}])^{2}-\text{tr}[\Delta_{+}^{\dagger}\Delta_{+}^{\phantom{\dagger}}\Delta_{+}^{\dagger}\Delta_{+}^{\phantom{\dagger}}]$, which allows us to set $b_{3}=0$ in the following without loss of generality.

Using the singular-value decomposition of $\Delta_{+}$, it is straightforward to find all symmetry-inequivalent minima of Eq. (12). There are two different states depending on the sign of $b_{2}$ which we label by $A_{m_{v}=0}(\Delta^{s};\boldsymbol{d})$, where $\Delta^{s}$ and $\boldsymbol{d}$ refer to the singlet and the triplet vector, respectively, $A$ indicates the trivial IR of $C_{3}$, and $m_{v}=0$ signifies intervalley pairing (IR of U(1)v with $m_{v}=0$). If $b_{2}>0$, we get $\Delta_{+}\propto\sigma_{0}$, i.e., $M_{\boldsymbol{k},\pm}=\lambda_{\pm\boldsymbol{k}}\sigma_{0}$ with $\lambda_{C_{3}\boldsymbol{k}}=\lambda_{\boldsymbol{k}}$; according to the notation introduced above, this state will be labeled as $A_{m_{v}=0}(1;0,0,0)$. There are (infinitely) many different equivalent representations of this state since, for instance, the transformations in Eq. (11b) mix the singlet and triplet components—as described by the isomorphism $\text{SU}(2)_{+}\times\text{SU}(2)_{-}\simeq\text{SO}(4)$. However, for the sake of notational clarity, we will henceforth only show one convenient representative of each state. The $A_{m_{v}=0}(1;0,0,0)$ state preserves time-reversal symmetry and breaks SU(2)+ $\times$ SU(2)- down to SU(2)s [rotations of the total spin, i.e., $\boldsymbol{\varphi}_{+}=\boldsymbol{\varphi}_{-}$ in Eq. (11b)].

On the other hand, if $b_{2}<0$, we find $\Delta_{+}\propto\sigma_{0}+\sigma_{z}$, which corresponds to $A_{m_{v}=0}(1;1,0,0)$. For this phase, the order parameter in Eq. (9) assumes the form $M_{\boldsymbol{k},\pm}=\lambda_{\pm\boldsymbol{k}}\left(\sigma_{0}\pm\sigma_{z}\right)$ with $\lambda_{C_{3}\boldsymbol{k}}=\lambda_{\boldsymbol{k}}$. This state preserves time-reversal symmetry too, but it breaks SU(2)+ $\times$ SU(2)- down to O(2)s (with $\boldsymbol{\varphi}_{+}=\boldsymbol{\varphi}_{-}=\varphi\,\hat{\boldsymbol{e}}_{z}$), i.e., rotations of the total spin along a single axis.

III.2 Turning on the Hund’s coupling

In reality, there is, of course, a finite Hund’s coupling that reduces $\mathcal{G}^{s}_{1}=$ SU(2)+ $\times$ SU(2)- to only global spin rotations, $\mathcal{G}^{s}_{2}=$ SU(2)s, already in the high-temperature phase. In Ref. Lee et al., 2019, the Hund’s coupling $J$ has been estimated to be about 60 times smaller than the intravalley interaction $V$. Note, however, $J$ might be enhanced due to loop corrections. For this reason, we first classify the possible instabilities in the absence of an approximate SU(2)+ $\times$ SU(2)- symmetry and then, analyze how the different states “connect” for small values of $J$ and whether admixtures of singlet and triplet are possible.

To introduce our notation, we will begin with the classification for the reduced symmetry group $\mathcal{G}_{2}$ in Eq. (3); in that case, we have either singlet or triplet pairing:

Singlet:

This corresponds to the $d_{n_{s}}$ $=$ $1$ one-dimensional IR of $\mathcal{G}^{s}_{2}$ with $\chi^{n_{s}}$ $=$ $\sigma_{0}$ in Eq. (8). The pairing Hamiltonian simply has the form

[TABLE]

with $\lambda^{s}_{\boldsymbol{k}v}$ $=$  $\lambda^{s}_{-\boldsymbol{k}\bar{v}}$ and $\lambda^{s}_{C_{3}\boldsymbol{k}v}$ $=$ $\lambda^{s}_{\boldsymbol{k}v}$. All symmetries of the high-temperature phase are preserved. We refer to this state as $A^{1_{s}}_{m_{v}=0}$ with the $1_{s}$ referring to spin singlet, $m_{v}=0$ to intervalley pairing, and $A$ to the trivial representation of $C_{3}$.

Triplet:

This pairing channel is associated with the three-dimensional IR of $\mathcal{G}^{s}_{2}$. A possible choice of basis functions is $\chi^{n_{s}}_{\mu}(v)=\sigma_{\mu}$, $\mu=1,2,3$, in Eq. (8). As it is a multidimensional representation, the free energy has to be expanded beyond quadratic order. Writing $\boldsymbol{d}=(\eta_{1}^{n_{s}},\eta_{2}^{n_{s}},\eta_{3}^{n_{s}})$, we have up to quartic order

[TABLE]

Observe that $\lvert\boldsymbol{d}^{T}\boldsymbol{d}\rvert^{2}$ is not an independent quartic term since $\lvert\boldsymbol{d}^{T}\boldsymbol{d}\rvert^{2}=(\boldsymbol{d}^{\dagger}\boldsymbol{d})^{2}-\lvert\boldsymbol{d}^{*}\times\boldsymbol{d}\rvert^{2}$. The free energy in Eq. (14) has two stable minima. For $b^{t}_{2}>0$, we have $\boldsymbol{d}\propto(1,0,0)^{T}$ and the corresponding pairing term is

[TABLE]

with $\lambda^{t}_{\boldsymbol{k}v}=-\lambda^{t}_{-\boldsymbol{k}\bar{v}}$ and $\lambda^{t}_{C_{3}\boldsymbol{k}v}=\lambda^{t}_{\boldsymbol{k}v}$. As is easily seen, this term preserves time-reversal symmetry and breaks SU(2)s down to spin rotation along a single axis. This state will be referred to as unitary triplet and denoted by the symbol $A^{3_{s}}_{m_{v}=0}(1,0,0)$, where the three components just indicate the direction of the triplet vector. If, instead, $b^{t}_{2}<0$, we obtain $\boldsymbol{d}\propto(1,i,0)^{T}$, whence

[TABLE]

with $\lambda^{t}_{\boldsymbol{k}v}$ as above. This is a nonunitary triplet state. It breaks time-reversal symmetry and will be denoted by $A^{3_{s}}_{m_{v}=0}(1,i,0)$.

One might wonder what kind of interaction or band structure would favor $A^{3_{s}}_{m_{v}=0}(1,i,0)$ over $A^{3_{s}}_{m_{v}=0}(1,0,0)$ and vice versa. In mean-field theory, as detailed in Appendix A, it is straightforward to show by evaluation of a one-loop diagram that

[TABLE]

Here, $\omega_{n}$ are fermionic Matsubara frequencies and $\xi_{\boldsymbol{k}+}$ is the electronic band energy in valley $v=+$ of the nearly flat band hosting superconductivity. We observe that $b^{t}_{2}>0$ holds irrespective of microscopic details and hence, $A^{3_{s}}_{m_{v}=0}(1,0,0)$ is generically favored if we neglect corrections beyond mean-field theory (such as residual interactions or frequency dependence of pairing). Intriguingly, there have been experimental reports Quintanilla et al. (2010) of intrinsically nonunitary pairing in LaNiC2, i.e., nonunitary triplet pairing born out of a paramagnetic normal state. Thus, there is reason to believe that we cannot generically exclude this state, but we do not expect it to show up in any simple mean-field computation.

III.2.1 How do the states connect in the $J=0$ limit?

Next, we establish how the three possible states, $A^{1_{s}}_{m_{v}=0}$, $A^{3_{s}}_{m_{v}=0}(1,0,0)$, and $A^{3_{s}}_{m_{v}=0}(1,i,0)$, connect to the two derived in the previous subsection with enhanced SU(2)+ $\times$ SU(2)- symmetry, namely $A_{m_{v}=0}(1;0,0,0)$ and $A_{m_{v}=0}(1;1,0,0)$. To this end, we decompose the Ginzburg-Landau expansion (12) into singlet and triplet by writing $\Delta_{+}=\Delta^{s}+\boldsymbol{\sigma}\cdot\boldsymbol{d}$. Since $\text{tr}\left[\Delta_{+}^{\dagger}\Delta_{+}^{\phantom{\dagger}}\right]=|\Delta^{s}|^{2}+\boldsymbol{d}^{\dagger}\boldsymbol{d}$, singlet and triplet are degenerate at quadratic order in $\mathcal{F}$ as a consequence of the enhanced SU(2)+ $\times$ SU(2)- symmetry. For nonzero $J$, this degeneracy is lifted and we have

[TABLE]

where $\delta a$ can be made arbitrarily small as $J\rightarrow 0$. Neglecting, for now, the “back action” of the superconducting order parameter that condenses first on the second one (as described by higher-order terms in the Ginzburg-Landau expansion), we conclude that there are two superconducting transitions at $T_{c,0}^{\pm}$ $=$  $T_{c,0}\pm\Delta T_{c}$ with $\Delta T_{c}=|\delta a\,(T_{c,0})|/\alpha$, taking $a(T)\sim\alpha(T-T_{c,0})$ near $T_{c,0}$. The extra index [math] in $T_{c,0}^{\pm}$ highlights the fact that the aforementioned higher-order terms in the Ginzburg-Landau expansion can significantly affect the lower transition temperature, $T_{c}^{-}\neq T_{c,0}^{-}$; of course, this has no effect on the higher transition temperature, $T_{c}^{+}=T_{c,0}^{+}$.

Before analyzing these corrections, it is useful to estimate the temperature scale $\Delta T_{c}$. Using the expected result, $T_{c}^{\pm}\simeq\Lambda\exp(-1/[(V\pm J)\nu])$ of mean-field theory (from the linearized gap equations)—where $\Lambda$ is the cutoff and $\nu$ the density of states at the Fermi level—leads to

[TABLE]

The large density of states, taken together with the estimated value of $V$—which is larger than even the bandwidth Lee et al. (2019) of the flat bands—and the relation $J\ll V$, implies that $\Delta T_{c}\ll T_{c,0}$ 111Taking $\nu\simeq 1/W$ with bandwidth $W\simeq 10\,\textrm{meV}$ and $V\simeq 35\,\textrm{meV}$, $J\simeq 0.6\,\textrm{meV}$ Lee et al. (2019), we estimate $\Delta T_{c}/T_{c,0}\simeq 0.5\%$.; the temperature/energy scale $\Delta T_{c}$ is most likely too small to be visible in experiments. While the estimate above is only based on mean-field theory, it indicates at least that it is important to study the behavior of superconductivity in the limit of small $\Delta T_{c}/T_{c,0}$ (and hence, weakly broken SU(2)+ $\times$ SU(2)- symmetry), accounting for the possibility of two transitions and mixing of singlet and triplet pairing (despite the absence of spin-orbit coupling). Moreover, we will see that nearly degenerate singlet and triplet pairing also has crucial consequences for the behavior of superconductivity in the presence of a magnetic field.

While we postpone the analysis of magnetic fields to Sec. III.3, here, we investigate the possibility of an admixture of singlet and triplet in the presence of time-reversal symmetry [relevant to scenario (I) defined in the introduction]. As anticipated above, this requires also considering the quartic terms of Eq. (12). We find

[TABLE]

neglecting corrections to the quartic terms coming from finite $J$.

Looking at the first transition with the higher transition temperature, we assess which of the two distinct triplet states, $A^{3_{s}}_{m_{v}=0}(1,0,0)$ and $A^{3_{s}}_{m_{v}=0}(1,i,0)$, and the singlet state can be stabilized by starting from $A_{m_{v}=0}(1;0,0,0)$ or $A_{m_{v}=0}(1;1,0,0)$ and turning on a finite Hund’s coupling $J$. For this purpose, we can neglect the coupling terms in the third line of Eq. (20). Clearly, if $\delta a<0$ (“anti-Hund’s coupling”), we get a singlet state for both $A_{m_{v}=0}(1;0,0,0)$ and $A_{m_{v}=0}(1;1,0,0)$. A straightforward way of establishing which of the triplet states is realized when $\delta a>0$ (“conventional” Hund’s coupling) proceeds by evaluating their respective free energy in Eq. (20). One finds that the state $A^{3_{s}}_{m_{v}=0}(1,0,0)$ is realized if $b_{2}>0$; otherwise, $A^{3_{s}}_{m_{v}=0}(1,i,0)$ is favored. This brings us to the conclusion that

[TABLE]

at the first transition (see the schematic phase diagram in Fig. 2). This result is just a consequence of the fact that the form $\Delta_{+}\propto\sigma_{0}$ for the $A_{m_{v}=0}(1;0,0,0)$ state we had chosen in the previous section can alternatively be written as $\Delta_{+}\propto\sigma_{x}$ due to the SU(2)+ $\times$ SU(2)- symmetry and thus, explicitly assumes the form of the unitary triplet state. Similarly, $\Delta_{+}\propto\sigma_{0}+\sigma_{z}$ used above for $A_{m_{v}=0}(1;1,0,0)$ can also be written as $\Delta_{+}\propto\sigma_{x}+i\sigma_{y}$. This is why it transitions into the nonunitary triplet state, upon turning on a nonzero Hund’s coupling.

In order to determine whether there is a second transition, we have to include the coupling terms between singlet and triplet in the third line of Eq. (20). To illustrate that these terms can be crucial, we consider the case $\delta a>0$ and $b_{2}>0$, i.e., the triplet state $A^{3_{s}}_{m_{v}=0}(1,0,0)$ condenses first. This leads to the coupling between singlet and triplet $2c|\Delta^{s}|^{2}|\boldsymbol{d}(T)|^{2}$, $c=b_{1}+b_{2}$, in the free energy, where we have made use of the fact that a relative phase of $\pi/2$ between singlet and triplet is energetically most favorable. As a result of $|\boldsymbol{d}(T)|^{2}=(\delta a-a(T))/(2c)$, which is valid as long as there is no additional singlet pairing, the growing triplet component induces the extra term

[TABLE]

which is always larger than the “bare” quadratic term of singlet pairing [in the first line of Eq. (20)]. Accordingly, there is no second transition (at least close to $T_{c,0}$ where our Ginzburg-Landau approach is valid) into a state that has a nonzero singlet component. We also checked that Eq. (20) does not allow for a first-order transition.

Similarly, all other cases can be scrutinized and one finds that if triplet dominates, there is no second transition. However, if singlet has a larger transition temperature ($\delta a<0$), there is a second transition into a phase with singlet and triplet pairing when $b_{2}<0$. This transition happens at the temperature

[TABLE]

The stability of the Ginzburg-Landau expansion only requires $c$ $>$ [math] and $c$ $>$ $-b_{2}$, so both $T_{c}^{-}$ $<$ $T_{c,0}^{-}$ and $T_{c}^{-}$ $>$ $T_{c,0}^{-}$ are possible. More importantly, unless $|b_{2}|/c$ is fine-tuned to be of order $\delta$, generically, $T_{c}^{-}\rightarrow T_{c,0}$ as $J\rightarrow 0$ and the two transitions, if present, are likely too close to be experimentally discernible. Due to the term $2b_{2}\,\text{Re}[\left(\Delta^{s}\right)^{2}\boldsymbol{d}^{\dagger}\boldsymbol{d}^{*}]$ in the free energy, we obtain the unitary triplet vector $\boldsymbol{d}=d_{0}(1,0,0)^{T}$ with $\Delta^{s}d_{0}^{*}\in\mathbb{R}$ (same phase). This is to be expected as $\Delta_{+}\propto\sigma_{0}+\sigma_{z}$ for the “parent” state $A_{m_{v}=0}(1;1,0,0)$.

A summary of these results is provided by the schematic phase diagrams in Fig. 2. We observe that the proximity to the enlarged symmetry in spin space, SU(2)+ $\times$ SU(2)-, favors the possibility of having a nonzero triplet component: for $b_{2}<0$, even a negative Hund’s coupling (anti-Hund’s) allows for $\boldsymbol{d}\neq 0$ and leads to the exotic possibility of significant ($d_{0}\simeq\Delta^{s}$ for $T_{c}^{-}-T>\Delta T_{c}$) singlet-triplet mixing in spite of the absence of spin-orbit coupling.

It is noteworthy that all the states are fully gapped (more precisely, they have no symmetry-enforced nodes) except for the nonunitary triplet $A^{3_{s}}_{m_{v}=0}(1,i,0)$, which is gapped for one spin species while the other is completely gapless. The admixture of singlet and unitary triplet has two unequal gaps for the two spin species both of which are finite as long as the magnitudes of singlet and triplet are not fine-tuned to be equal. All the states, along with their order parameters and properties, are summarized in Table 1.

We finally comment on the nature of the thermal phase transition for the different superconducting states once fluctuations of the order parameter are taken into account. Neglecting stray fields, the transition into the singlet phase $A^{1_{s}}$ is expected to be a BKT transition with quasi-long-range order of the complex-valued order parameter $\Delta^{s}$ below the transition temperature. For the triplet states, it is important to keep in mind that $\boldsymbol{d}$ cannot even have quasi-long-range order as it transforms as a three-component vector under spin-rotation. For the unitary triplet state [with order parameter manifold $(S_{2}\times S_{1})/\mathbb{Z}_{2}$] a BKT transition of the composite charge-$4e$ order parameter $\boldsymbol{d}^{T}\boldsymbol{d}$ is possible and is associated with the (un)binding of half vortices. This is different for the nonunitary state [with order parameter manifold $S_{3}/\mathbb{Z}_{2}\simeq\text{SO}(3)$] where $\boldsymbol{d}^{T}\boldsymbol{d}=0$ and no BKT transition into a quasi-long-range-ordered superconductor is expected. For the case of the two consecutive transitions in Fig. 2(b) with $\delta<0$, we first expect a BKT transition into a singlet phase followed by a crossover at which the triplet vector becomes nonzero.

However, we point out that, even in the simplest case of the singlet $A^{1_{s}}$, there are significant corrections to the BKT transition resulting from stray fields and mirror vortices Kogan (2007), which make the observation of a pristine BKT transition in a (charged) superconductor difficult. We believe that the current status of experiments does not allow one to exclude pairing phases that will not exhibit quasi-long-range order and a BKT transition in the limit of infinite system size.

III.2.2 Expectations within mean-field theory

Lastly, we evaluate what a naïve mean-field computation is expected to yield. In fact, from Eq. (17), we already know that the prefactor of the term $\left|\boldsymbol{d}^{*}\times\boldsymbol{d}\right|^{2}$ in Eq. (20) must be positive within mean-field theory and therefore, it holds that $b_{2}>0$. For completeness, we mention that in the mean-field approximation, $b_{1}=0$, as shown in Appendix A. Consequently, a single-band mean-field computation will generally favor Fig. 2(a) over (b); in other words, only half of the phases proposed in this section can be found in mean-field, which we also indicate in the last column of Table 1.

However, there is no fundamental mechanism prohibiting the mixing of singlet and triplet via two transitions (see, e.g., Ref. Mráz and Hlubina, 2005) and there are multiple reasons why we can effectively have $b_{2}<0$ (and $b_{1}>0$ to ensure stability): for instance, strong residual interactions and fluctuations have been shown to modify the values of the quartic terms in the free energy significantly Fernandes and Millis (2013); Kozii et al. (2019), thereby stabilizing phases that are otherwise not possible in the mean-field approximation. Given the small bandwidth and the underlying strong-coupling features of the problem Kerelsky et al. (2019); Choi et al. (2019); Jiang et al. (2019); Xie et al. (2019); Scheurer (2019), it is plausible that there are sizable corrections to mean-field theory. In addition, we recognize that there are other corrections arising from interband pairing, and that disorder can also dress the Ginzburg-Landau expansion. Moreover, it is unclear whether adding frequency dependence to the gap function could be of relevance.

In Sec. V, we will analyze the impact of ferromagnetic fluctuations, which are expected to be relevant for graphene moiré systems Shen et al. (2019); Liu et al. (2019b); Cao et al. (2019); Sharpe et al. (2019); Lu et al. (2019); Zondiner et al. (2019); Chen et al. (2020b), and find that these generically decrease the value of $b_{2}$; if sufficiently strong, these fluctuations will favor the phase diagram in Fig. 2(b).

III.3 In the presence of a magnetic field

We now generalize the Ginzburg-Landau expansion to also include the coupling to a Zeeman field $\boldsymbol{M}_{Z}=(M_{Z}^{x},M_{Z}^{y},M_{Z}^{z})$ and an (in-plane) orbital coupling $\boldsymbol{M}_{O}=(M_{O}^{x},M_{O}^{y})$. Both of these terms can either be due to an applied external magnetic field or due to the correlated insulating state. This enables us to discuss (i) the behavior of the superconducting critical temperature $T_{c}^{+}$ as a function of an external magnetic field in the absence of any ferromagnetic moments associated with the correlated insulating state [case (I) defined in the introduction]. At the same time, we can study (ii) how the transition temperature and the order parameter of superconductivity is affected by the potentially coexisting ferromagnetic order [case (II)].

III.3.1 Leading superconducting transition

We first turn our attention to the leading superconducting transition with the highest temperature $T_{c}^{+}$; potential subsequent superconducting transitions at lower temperatures are addressed later in Sec. III.3.2. For the goal of studying the first transition, we can restrict ourselves to quadratic order in the order parameter. Only keeping terms up to quadratic order in the magnetic field as well, we obtain

[TABLE]

While the prefactors $\delta a$, $\delta c_{1}$, $\delta c_{5}$, and $\delta c_{6}$ are necessarily zero in the limit $J\rightarrow 0$, where the SU(2)${}_{+}\times$ SU(2)- symmetry becomes exact, all remaining terms can be nonzero (and different in their values) at $J=0$. Notice that the third term has not been considered in Ref. Lee et al., 2019; this term arises only when both singlet and triplet are allowed for and leads to the admixture of a unitary triplet state with a singlet superconductor. The vanishing of $\delta a$ and $\delta c_{6}$ at $J=0$ is an obvious consequence of the enhanced SU(2)${}_{+}\times$ SU(2)- symmetry. To see that $\delta c_{1}$ also has to vanish as $J\rightarrow 0$, let us take $\boldsymbol{M}_{Z}$ along the $z$ direction; this breaks SU(2)${}_{+}\times$ SU(2)- down to $O(2)_{+}\times O(2)_{-}$, i.e., the system is only invariant under $c_{\boldsymbol{k}v}\rightarrow e^{i\varphi_{v}\sigma_{z}}c_{\boldsymbol{k}v}$. Performing this transformation with $\varphi_{+}=0$ and $\varphi_{-}=\pi/2$, we get $(\Delta^{s},d_{z})\rightarrow(id_{z},i\Delta^{s})$ and hence, $\delta c_{1}\rightarrow-\delta c_{1}$. With the same argument, it can be proven that $\delta c_{5}$ has to go to zero as $J\rightarrow 0$. In Appendix A, we show that $\delta c_{1}=0$ in mean-field theory within the single-band description, even when SU(2)${}_{+}\times$ SU(2)- is broken; this results from an emergent valley-exchange symmetry within the single-band mean-field approximation.

In discussing the highest critical temperature and the corresponding order parameter for $\boldsymbol{M}_{Z},\boldsymbol{M}_{O}\neq 0$, it is instructive to first look at the linear-in-field terms in Eq. (24). We find two different cases. If $|c_{2}M_{Z}|+\delta a>\sqrt{(\delta c_{1}M_{Z})^{2}+\delta a^{2}}$, one obtains a pure triplet state of the type $A^{3_{s}}_{m_{v}=0}(1,i,0)$. Choosing $\boldsymbol{M}_{Z}=M_{Z}\boldsymbol{e}_{x}$ with $M_{Z}>0$, the triplet vector is given by $\boldsymbol{d}=\left(0,1,\,\text{sign}(c_{2})i\right)^{T}$ and the critical temperature is

[TABLE]

Else, if $|c_{2}M_{Z}|+\delta a<\sqrt{(\delta c_{1}M_{Z})^{2}+\delta a^{2}}$, one finds an admixture of singlet and triplet with order parameter

[TABLE]

The transition temperature in this case is

[TABLE]

We see from Eq. (26) that there is an approximately equal mixing of singlet and triplet for $|\delta c_{1}|M_{Z}\gg|\delta a|$ while in the opposite limit, $|\delta c_{1}|M_{Z}\ll|\delta a|$, either singlet or triplet dominates depending on whether $\delta a<0$ or $\delta a>0$. The relative phase of $\pi/2$ between singlet and triplet makes the pairing state break time-reversal symmetry as is required in order to couple linearly to magnetic moments.

To understand how the approximate SU(2)+ $\times$ SU(2)- symmetry can naturally explain the linear-in-magnetic-field behavior, we first consider case (II), i.e., there is already microscopically coexisting ferromagnetic order (or there is at least a significant coupling between superconductivity and the ferromagnetic moments) at $T_{c}^{+}$. Then, $\boldsymbol{M}_{Z}$ and $\boldsymbol{M}_{O}$ should be thought of as the combination of the applied external magnetic field and the ferromagnetic order parameter. In this scenario, it is apt to assume $|\delta a|\ll\text{max}(\delta c_{1}M_{Z},c_{2}M_{Z})$ and we generically obtain a linear increase of the critical temperature with magnetic field [see Eqs. (25) and (27)]. If $c_{2}>\delta c_{1}$, we obtain the nonunitary triplet state with $\boldsymbol{d}\propto(1,i,0)^{T}$, which we expect close to the $J=0$ line, while $\delta c_{1}>c_{2}$ leads to the admixture of singlet and triplet with $\boldsymbol{d}\propto(1,0,0)^{T}$. As $\delta c_{1}$ vanishes for $J=0$ and in single-band mean-field theory (even when $J\neq 0$), we expect the former scenario to be more likely, which will favor the nonunitary triplet state as the leading instability.

In the case of scenario (I), we should view $\boldsymbol{M}_{Z}$ and $\boldsymbol{M}_{O}$ in Eq. (24) as resulting entirely from the Zeeman and orbital coupling of the external magnetic field alone. For large magnetic fields where $|\delta a|\ll\text{max}\,(\delta c_{1}M_{Z},c_{2}M_{Z})$, the same conclusions as above will apply and $T_{c}$ will generically vary linearly with the field. However, for sufficiently small magnetic fields, we have $|\delta a|\gg\text{max}(\delta c_{1}M_{Z},c_{2}M_{Z})$. In this limit, only $\delta a>0$ favoring the nonunitary triplet pairing $A^{3_{s}}_{m_{v}=0}(1,i,0)$ is consistent with the transition temperature changing linearly with magnetic field. Alternatively, the system could ultimately be in a singlet state at $M_{Z}=0$ (i.e., $\delta a<0$) but the magnitude of $\delta a$ is sufficiently small such that the “rounding off” of $T^{+}_{c}(M_{Z})$ at low $M_{Z}$ cannot be seen in experiment.

III.3.2 Quartic terms and sub-leading transitions

Having examined the first superconducting transition that takes place upon cooling the system down starting from the normal state, we now assess whether and what type of subsequent superconducting transitions can occur. In this context, we need to include terms quartic in the superconducting order parameter and extend Eq. (24) to

[TABLE]

where we kept only the terms linear in magnetic field, took $\boldsymbol{M}_{Z}$ along the $z$-axis, and re-expressed the triplet in the form $\boldsymbol{d}$ $=$ $d_{+}(1,i,0)/\sqrt{2}+d_{-}(1,-i,0)/\sqrt{2}+d_{0}(0,0,1)$. This parametrization is more convenient in the presence of a magnetic field than that used in Eq. (20). Additionally, we have neglected the impact of the magnetic field on the quartic terms.

Taking $\delta c_{1}$ $=$ [math] (as it has to vanish for $J$ $=$ [math]), the different possible phase diagrams are summarized in Fig. 3. The possibility illustrated in part (c) of Fig. 3 corresponds to the picture put forward by Ambegaokar and Mermin (1973) for He3 in the presence of a magnetic field, which might very well also apply to twisted double-bilayer graphene Liu et al. (2019b); Lee et al. (2019). The difference with Ref. Ambegaokar and Mermin, 1973 is that we do not get a third transition since we work with a one-dimensional IR of the spatial point group.

However, there are three other options, depicted in Fig. 3(a), (b), and (d), that we cannot easily exclude given the experimental data: owing to the strong-coupling properties of the problem at hand, a nonunitary triplet state might be dominant at $M_{Z}=0$, as seems to be the case in LaNiC2 Quintanilla et al. (2010) and is favored by our fluctuation approach of Sec. V; under this condition, only one transition is expected even when $M_{Z}\neq 0$ [see Fig. 3(d)]. It could also be that singlet dominates without a magnetic field instead. We can see in Fig. 3(a) and (b) that, in these two cases, triplet shows up and $T_{c}$ increases linearly when $M_{Z}>2|\delta a|/c_{2}$. The small value of $\Delta T_{c}/T_{c,0}$ estimated in Eq. (19) suggests that resolving this initial region, where $T_{c}$ is constant as a function of $M_{Z}$, is experimentally challenging.

III.3.3 Nonlinear couplings in a magnetic field

We finally come back to the quadratic couplings to the magnetic field, associated with the terms with prefactors $c_{3,4}$ and $\delta c_{5,6}$ in Eq. (24). We first notice that this will lead to an additional quadratic suppression of the leading transition temperatures in Fig. 3; in particular, the transition temperature into the singlet state in part (a) and (b) will not be field-independent any more. More interestingly, the suppression of singlet and triplet is enforced to be nearly identical for small $J$ due to the SU(2)+ $\times$ SU(2)- symmetry, $\delta c_{5,6}\ll c_{3,4}$. Resultantly, if the effective $J$ relevant for superconductivity is indeed small, the nonlinear terms $\propto\boldsymbol{M}_{Z}^{2},\boldsymbol{M}_{O}^{2}$ are not expected to affect the competition between singlet and triplet significantly and the qualitative form of the phase diagrams in Fig. 3 is not modified.

IV Complex representation of $C_{3}$

In this section, we extend our previous analysis to the complex IR $E$ of the spatial point group $C_{3}$. Time-reversal symmetry necessitates treating the representation and its complex-conjugate partner on an equal footing. Alternatively, one can think of a two-dimensional (reducible) representation with partner functions transforming as $x$ and $y$ under $C_{3}$.

Akin to our discussion earlier, we first study the case of nonzero Hund’s coupling, $J\neq 0$, with point group $\mathcal{G}_{2}$ in Eq. (3), which enables us to distinguish between singlet and triplet pairing. After discussing all symmetry-allowed singlet and triplet states separately, we will derive the phase diagrams analogous to Fig. 2: we will examine how these states “connect” when adiabatically changing the Hund’s coupling from negative to positive values, and whether singlet and triplet can mix when $J$ is small and the SU(2)+ $\times$ SU(2)- symmetry is only weakly broken.

IV.1 Nonzero Hund’s coupling

To proceed with singlet pairing, we parametrize $M_{\boldsymbol{k}v}$ in Eq. (9) according to

[TABLE]

while $M_{\boldsymbol{k}-}$ is determined by the Fermi-Dirac constraint (10); $X_{\boldsymbol{k}}$ and $Y_{\boldsymbol{k}}$ are real-valued functions that are continuous on the Brillouin zone and transform as $k_{x}$ and $k_{y}$ under $C_{3}$. A one-parameter family of possible choices for the lowest-order functions (i.e., with minimal number of sign changes in the Brillouin zone) is given by

[TABLE]

Both $X_{\boldsymbol{k}}$ and $Y_{\boldsymbol{k}}$ have to vanish at $\Gamma$, $K$ and $K^{\prime}$ as these momenta are invariant under $C_{3}$. Further, both $X_{\boldsymbol{k}}$ and $Y_{\boldsymbol{k}}$ must have lines of zeros going through these high symmetry points; the orientation of these lines is, however, not fixed due to the absence of additional reflection or in-plane rotation symmetries—this is different from the situation for twisted bilayer and trilayer graphene in Sec. VI. For Eq. (30), the orientation of these zeros changes with $\phi$.

With the parametrization defined in Eq. (29), the relevant symmetries act as follows

[TABLE]

It readily follows from Eq. (31) that the most general free energy up to quartic order reads as

[TABLE]

The sign of $b^{s}_{2}$ therefore distinguishes between two different singlet phases: if $b^{s}_{2}>0$, we have $(\eta_{+},\eta_{-})=(1,0)$, which corresponds to

[TABLE]

Exactly as in Sec. III, we always show only one out of the many symmetry-equivalent representations of the order parameter—instead of using a general parametrization of a phase—to make the notation and the discussion of properties of the superconducting state more easily accessible. The state in Eq. (33) breaks time-reversal symmetry but preserves $C_{3}$ (and spin-rotation symmetry). We refer to this state as a chiral singlet superconductor and denote it by $E^{1s}(1,i)$ in the following. It is fully gapped (unless the Fermi surfaces go through the $\Gamma$, $K$, or $K^{\prime}$ point) and has been investigated extensively in the recent literature on pairing in twisted bilayer graphene You and Vishwanath (2019); Venderbos and Fernandes (2018); Wu and Das Sarma (2019); Kennes et al. (2018); Guo et al. (2018); Huang et al. (2019); Chen et al. (2020a); Lin and Nandkishore (2019); Liu et al. (2018); Fidrysiak et al. (2018).

Conversely, if $b^{s}_{2}<0$, we find that $|\eta_{+}|=|\eta_{-}|$ at the minimum of Eq. (32). As the relative phase $\varphi$ between $\eta_{+}$ and $\eta_{-}=\eta_{+}e^{i\varphi}$ is not fixed by Eq. (32), one might naively conclude that higher order terms have to be considered. In fact, in sixth order, there is indeed the contribution

[TABLE]

and the relative phase $\varphi$ will depend on $c_{1}/c_{2}$. However, upon reinserting $\eta_{-}=\eta_{+}e^{i\varphi}$ into Eq. (29), we notice that $\varphi\neq 0$ simply corresponds to rotating the basis functions $X_{\boldsymbol{k}}$ and $Y_{\boldsymbol{k}}$ into each other, which does not change their transformation behavior under $C_{3}$ [$\varphi$ is directly related to $\phi$ in Eq. (30a)]. Consequently, we can set $\varphi=0$ without loss of generality, which implies

[TABLE]

This state, which we call $E^{1s}(1,0)$, breaks $C_{3}$ but preserves time-reversal symmetry; this is the nematic singlet phase.

Within a single-band mean-field description (see Appendix A), we find $b^{s}_{1}=b^{s}_{2}/2>0$. As such, mean-field theory generically favors the chiral singlet superconductor over the nematic state $E^{1s}(1,0)$; this has been noted before in the context of twisted bilayer graphene You and Vishwanath (2019) and Ref. Kozii et al., 2019 discusses how strong fluctuations can stabilize the nematic phase.

Turning to triplet pairing, we now modify the parametrization (29) to

[TABLE]

where $X_{\boldsymbol{k}}$ and $Y_{\boldsymbol{k}}$ are defined exactly as before. For simplicity, we introduce the complex-vector notation, $\boldsymbol{d}_{\mu}=(\eta_{\mu,1},\eta_{\mu,2},\eta_{\mu,3})^{T}$, $\mu=\pm$. The representations of the symmetries now read as

[TABLE]

with $\mathcal{R}\in\text{SO}(3)$ and $\omega=e^{i\frac{2\pi}{3}}$. The most general free-energy expansion is given by

[TABLE]

up to quartic order, where $b^{t}_{j}$ $\in$ $\mathbb{R}$; the different symmetry-allowed phases follow from the stable minima of the free energy. When minimizing Eq. (38), we take into account that the relative phase between $\boldsymbol{d}_{+}$ and $\boldsymbol{d}_{-}$ can always be absorbed into a redefinition of the basis functions $X_{\boldsymbol{k}}$ and $Y_{\boldsymbol{k}}$, as for the singlet above. In total, we find eight distinct triplet states which we label by $E^{3_{s}}(a)$ through $E^{3_{s}}(h)$. Phase diagrams describing which of these phases is realized for a given configuration of the quartic couplings $b^{t}_{j}$ can be found in Appendix C; here, we list all the phases, describe their properties, and refer to Fig. 4 for an illustration of their respective spectra and densities of states:

1. (a)

This state, labeled as $E^{3_{s}}(a)$, can be represented by $\boldsymbol{d}_{+}$ $=$ $\boldsymbol{d}_{-}$ $=$ $\left(1,0,0\right)^{T}$ with the associated order parameter $M_{\boldsymbol{k}+}$ $=$ $X_{\boldsymbol{k}}\,\sigma_{x}$. More physically, it corresponds to a nematic unitary triplet phase. It preserves time-reversal symmetry, but breaks both SU(2)s spin-rotation symmetry [down to O(2)] and $C_{3}$ rotational symmetry. This state has two symmetry-enforced nodal points at each Fermi surface around the $K$, $K^{\prime}$, or $\Gamma$ point. Owing to the lack of any reflection symmetry (cf. the discussion of $D_{3}$ in Sec. VI below), the positions of these nodal points are not pinned to any specific direction. 2. (b)

One representative configuration of this phase is given by $\boldsymbol{d}_{+}$ $=$ $(1,-i,0)^{T}/2$ and $\boldsymbol{d}_{-}$ $=$ $(1,i,0)^{T}/2$; it can thus be seen as a helical triplet, consisting of two time-reversed copies of states with opposite chirality. The order parameter can be more explicitly written as $M_{\boldsymbol{k}+}$ $=$ $X_{\boldsymbol{k}}\sigma_{x}+Y_{\boldsymbol{k}}\sigma_{y}$, which can alternatively be thought of as a 2D analogue of the Balian-Werthamer state of the B-phase of superfluid 3He Vollhardt and Wolfle (2013). This state, denoted by $E^{3_{s}}(b)$ in the following, only has point nodes at $\Gamma$, $K$, and $K^{\prime}$, i.e., it is expected to exhibit a full gap for generic Fermi surfaces not going through these high-symmetry points. It preserves time-reversal symmetry. While this state breaks spin-rotation symmetry as well as $C_{3}$, the product of $C_{3}$ and a rotation in spin space along $\sigma_{z}$ with angle $2\pi/3$ is preserved; this can be viewed as the spontaneous formation of spin-orbit coupling. 3. (c)

Here, we can write $\boldsymbol{d}_{+}$ $=$ $\boldsymbol{d}_{-}$ $=$ $(1,i,0)^{T}/2$; hence, $M_{\boldsymbol{k}+}$ $=$ $X_{\boldsymbol{k}}(\sigma_{x}+i\sigma_{y})$. This is a nematic nonunitary triplet state which breaks time-reversal symmetry and $C_{3}$. One spin-species will be gapless while the other will have nodal lines (i.e., point nodes on the Fermi surface). 4. (d)

The triplet vectors in this phase can be written as $\boldsymbol{d}_{+}$ $=$ $(1,0,0)^{T}$, $\boldsymbol{d}_{-}$ $=$ [math] leading to $M_{\boldsymbol{k}+}$ $=$ $\left(X_{\boldsymbol{k}}+i\,Y_{\boldsymbol{k}}\right)\sigma_{x}$. As one of the two chiralities is preferred over the other ($|\boldsymbol{d}_{+}|\neq|\boldsymbol{d}_{-}|$), this state can be referred to as chiral unitary triplet. It is a 2D analogue of the A-phase of 3He Vollhardt and Wolfle (2013). It breaks SU(2)s spin-rotation symmetry [down to O(2)] and time-reversal, but preserves $C_{3}$. Except for $\Gamma$, $K$, and $K^{\prime}$, this state has no symmetry-imposed nodal points. In fact, its spectrum is identical to that of the helical triplet $E^{3_{s}}(b)$, which is why we group these two states together in Fig. 4. 5. (e)

For this state, we have $\boldsymbol{d}_{+}$ $=$ $(1,i,0)^{T}$, $\boldsymbol{d}_{-}$ $=$ [math], i.e., $M_{\boldsymbol{k}+}$ $=$ $\left(X_{\boldsymbol{k}}+i\,Y_{\boldsymbol{k}}\right)\left(\sigma_{x}+i\sigma_{y}\right)$. It consists of only one of the two time-reversed copies with opposite chirality of the $E^{3_{s}}(b)$ state discussed above and, thus, is a chiral nonunitary triplet state. This state can be seen as an analogue of the $A_{1}$-phase of 3He Vollhardt and Wolfle (2013). It preserves $C_{3}$, but breaks SU(2)s spin-rotation symmetry [down to O(2)] and time-reversal. Here, one of the spin components will be gapless while the other is fully gapped (as before, except for the high symmetry points $\Gamma$, $K$, and $K^{\prime}$ which are generically not on the Fermi surface). Note that although the spectrum of this state is not strictly identical to that of the nematic nonunitary triplet $E_{s}^{3_{s}}(c)$, we grouped them together in Fig. 4 as their respective plots are practically indistinguishable; this is related to the fact that, in both cases, the low-energy spectrum is dominated by the Fermi surface of one of the spin species. 6. (f)

In this phase, $\boldsymbol{d}_{+}$ $=$ $(1,0,0)^{T}$, $\boldsymbol{d}_{-}$ $=$ $(0,1,0)^{T}$, implying $M_{\boldsymbol{k}+}$ $=$ $\left(X_{\boldsymbol{k}}+i\,Y_{\boldsymbol{k}}\right)\sigma_{x}$ $+$ $\left(X_{\boldsymbol{k}}-i\,Y_{\boldsymbol{k}}\right)\sigma_{y}$. The state can, thus, be thought of as a superposition of two chiral unitary triplets with orthogonal spin polarizations or, when inserted into Eq. (9), as Cooper pairs of electrons with spin polarization $\downarrow\downarrow$ ($\uparrow\uparrow$) and orbital basis function $X_{\boldsymbol{k}}+Y_{\boldsymbol{k}}$ ($X_{\boldsymbol{k}}-Y_{\boldsymbol{k}}$). Time-reversal, $C_{3}$, and spin-rotation symmetry are all broken. The excitation spectrum is given by $E_{\pm}(\boldsymbol{k})=\sqrt{\xi_{\boldsymbol{k}+}^{2}+2(X_{\boldsymbol{k}}\pm Y_{\boldsymbol{k}})^{2}}$, so it is characterized by “two gaps”, given by $|X_{\boldsymbol{k}}\pm Y_{\boldsymbol{k}}|$, both of which are forced to vanish at two points for each Fermi surface enclosing $K$, $K^{\prime}$, and $\Gamma$. While the number of nodes of this state and of $E^{3_{s}}(a)$ are the same, the spin degrees of freedom on the Fermi surface have nodes at the same two momenta for $E^{3_{s}}(a)$. For $E^{3_{s}}(f)$, however, the two spin species have nodal points at different momenta. 7. (g)

Denoted by $E^{3_{s}}(g)$, this phase has $\boldsymbol{d}_{+}$ $=$ $\cos\,(\alpha)$ $(1,i,0)^{T}/\sqrt{2}$,  $\boldsymbol{d}_{-}=\sin(\alpha)\,(0,0,1)^{T}$, where the parameter $\alpha$ varies continuously with $b^{t}_{j}$ in the part of the phase diagram where this state is realized. The corresponding order parameter can be written as $M_{\boldsymbol{k}+}=\cos(\alpha)\,\left(X_{\boldsymbol{k}}+i\,Y_{\boldsymbol{k}}\right)(\sigma_{x}+i\sigma_{y})/\sqrt{2}+\sin(\alpha)\,\left(X_{\boldsymbol{k}}-i\,Y_{\boldsymbol{k}}\right)\sigma_{z},$ $0<\alpha<\pi$, and can be viewed as a superposition of a chiral nonunitary triplet state and a unitary state with opposite chirality. This state breaks time-reversal symmetry, spin-rotation invariance, and $C_{3}$ but preserves the product of $C_{3}$ and spin rotation by angle $2\pi/3$ along $\sigma_{z}$. So, similar to the state $E^{3_{s}}(b)$ above, this state spontaneously entangles rotations in spin and real space and its spectrum, see Fig. 4(g), is $C_{3}$ invariant. It is fully gapped (again, as long as the Fermi surfaces do not go through $\Gamma$, $K$, and $K^{\prime}$), with two different gaps $[(1\pm g_{\alpha})\left(X^{2}_{\boldsymbol{k}}+Y^{2}_{\boldsymbol{k}}\right)]^{1/2}$, where $g_{\alpha}=\cos\alpha\sqrt{1+\sin^{2}\alpha}$. 8. (h)

Finally, for the triplet phase $E^{3_{s}}(h)$, one has $\boldsymbol{d}_{+}$ $=$ $(\cos\alpha,0,i\sin\alpha)^{T}$, $\boldsymbol{d}_{-}=(0,\cos\alpha,-i\sin\alpha)^{T}$, which yields $M_{\boldsymbol{k}+}=\cos(\alpha)[(X_{\boldsymbol{k}}+i\,Y_{\boldsymbol{k}})\sigma_{x}+(X_{\boldsymbol{k}}-i\,Y_{\boldsymbol{k}})\sigma_{y}]-2\sin(\alpha)Y_{\boldsymbol{k}}\sigma_{z}$. It can be seen as a superposition of the states $E^{3_{s}}(a)$ and $E^{3_{s}}(f)$ to which it reduces for $\alpha=\pi/2$ and $\alpha=0$; it will have two nodal points for $\alpha$ close to these limiting cases, but can be fully gapped for other values of $\alpha$. For $\alpha\neq\pi/2$, this state breaks time-reversal, $C_{3}$, and spin-rotation symmetry.

In Appendix A, we show that $b^{t}_{1}=b^{t}_{3}/2=-b^{t}_{4}/2=-2b^{t}_{5}>0$ and $b^{t}_{2}=0$ within a single-band mean-field description. Minimizing Eq. (38) yields that the phases $E^{3_{s}}(b)$ and $E^{3_{s}}(d)$ have the lowest energy and are exactly degenerate for this configuration of quartic couplings. This degeneracy within mean-field theory, which was noted before in Ref. Xu and Balents, 2018, will be lifted by corrections resulting, e.g., from residual interactions. In Sec. V, we will find that $E^{3_{s}}(b)$ ($E^{3_{s}}(d)$) is favored in the presence of ferromagnetic spin (orbital) fluctuations. We will also see that significant fluctuations can stabilize phases other than the two, $E^{3_{s}}(b)$ and $E^{3_{s}}(d)$, favored in mean-field theory.

IV.2 Approximate SU(2)+ $\times$ SU(2)-

After having classified singlet and triplet separately, we now focus on small Hund’s coupling for which SU(2)+ $\times$ SU(2)- is an approximate symmetry, and singlet and triplet are nearly degenerate at the quadratic level of the free energy. This requires studying them on an equal footing and generalizing the parametrization in Eqs. (29) and (36) to include both singlet and triplet, i.e., extending the summation over $\nu$ in Eq. (36) to $\nu=0,1,2,3$. In analogy with Sec. III.1, we use $2\times 2$ matrices and write

[TABLE]

It is easy to see that the symmetries act according to

[TABLE]

where, recall, $\mathcal{G}^{s}_{1}\equiv\text{SU}(2)_{+}\times\text{SU}(2)_{-}$. Imposing $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ as an exact symmetry, the most general free energy up to quartic order reads as

[TABLE]

At first glance, one might think that there are additional terms with extra factors of $\sigma_{y}$, similar to the last term in Eq. (12). However, as before, all of them can be related to combinations of the terms already present in Eq. (41) as outlined in Appendix C.

Following the procedure applied in Sec. III to the one-dimensional IR $A$, we now add a small quadratic term, $\delta a\sum_{\mu}(|\Delta^{s}_{\mu}|^{2}-\boldsymbol{d}_{\mu}^{\dagger}\boldsymbol{d}^{\phantom{\dagger}}_{\mu})$, where $\Delta_{\mu}^{s}$ and $\boldsymbol{d}_{\mu}$ are the singlet and triplet component of $\Delta_{\mu}$ in Eq. (41), i.e., $\Delta_{\mu}=\sigma_{0}\Delta^{s}_{\mu}+\boldsymbol{\sigma}\cdot\boldsymbol{d}_{\mu}$. This term breaks $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ and hence, makes singlet and triplet inequivalent. It allows us to study which of the different singlet and triplet states defined above can mix, and to identify “Hund’s partners”, i.e., which states transform into each other when changing the sign of the Hund’s coupling $J$ and accordingly, of $\delta a$. This generalizes the phase diagrams in Fig. 2 and Table 1 to the complex representation.

We find that, out of the eight different triplet states $E^{3_{s}}(a)$ to $E^{3_{s}}(h)$, only two—$E^{3_{s}}(a)$ and $E^{3_{s}}(d)$—do not allow for a singlet-triplet admixture when reversing the sign of $J$ (or $\delta a$) so that singlet has the higher transition temperature. The reason for the absence of an admixture is the same as sketched by way of example in Sec. III.2: besides pure singlet and pure triplet terms, the quartic terms in Eq. (41) also contain couplings between singlet and triplet, as is readily seen by inserting the parametrization $\Delta_{\mu}=\sigma_{0}\Delta^{s}_{\mu}+\boldsymbol{\sigma}\cdot\boldsymbol{d}_{\mu}$, $\mu=\pm$ (the full expansion can be found in Appendix C). At the first transition, one of either singlet or triplet becomes nonzero and hence, “renormalizes” the quadratic term of the other channel. In some cases, this renormalization can prohibit the presence of a second transition. In the case of phases $E^{3_{s}}(a)$ and $E^{3_{s}}(d)$, we just obtain the pure singlets $E^{1_{s}}(1,0)$ and $E^{1_{s}}(1,i)$, respectively, without a second transition. The easiest way to interpret why we do not have an admixture in these cases is to look at the associated SO(4) parent states: the two triplets correspond to $(\eta_{+};\boldsymbol{d}_{+})=(\eta_{-};\boldsymbol{d}_{-})=(0;1,0,0)$ and $(\eta_{+};\boldsymbol{d}_{+})=(0;1,0,0)$, $(\eta_{-};\boldsymbol{d}_{-})=0$, respectively. Both of these configurations can be “rotated” into the pure singlets $(\eta_{+};\boldsymbol{d}_{+})=(\eta_{-};\boldsymbol{d}_{-})=(1;0,0,0)$ and $(\eta_{+};\boldsymbol{d}_{+})=(1;0,0,0)$, $(\eta_{-};\boldsymbol{d}_{-})=0$ via a $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ transformation.

For all other triplets, the Hund’s partner is an admixed phase. Specifically, as regards $E^{3_{s}}(b)$ and $E^{3_{s}}(c)$, the Hund’s partner is an admixture of a nematic singlet state and a nematic unitary triplet $E^{3_{s}}(a)$, with different relative phases and spatial orientations: for the former, the order parameter can be written as $i\,Y_{\boldsymbol{k}}\sigma_{0}+\eta X_{\boldsymbol{k}}\sigma_{x}$, where $\eta$ describes the temperature-dependent strength of mixing, while it is $X_{\boldsymbol{k}}(\sigma_{0}+\eta\sigma_{x})$ for the latter. On any Fermi surface around one of the high-symmetry points $\Gamma$, $K$, or $K^{\prime}$, these two states have zero and two nodal points, respectively. Again, the form of the admixed state can be understood from the representation of the triplet state in terms of $(\eta_{\mu};\boldsymbol{d}_{\mu})$. For instance, we have $(\eta_{+};\boldsymbol{d}_{+})=(0;1,-i,0)$, $(\eta_{-};\boldsymbol{d}_{-})=(0;1,i,0)$ for $E^{3_{s}}(b)$, which is equivalent to $(\eta_{+};\boldsymbol{d}_{+})=(1;1,0,0)$, $(\eta_{-};\boldsymbol{d}_{-})=(-1;1,0,0)$ after applying an appropriate $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ transformation.

Likewise, the Hund’s partners of $E^{3_{s}}(e)$ and $E^{3_{s}}(f)$ are admixtures of a chiral singlet and a unitary triplet state with the same and opposite chirality, respectively. The associated order parameters can be written as $(X_{\boldsymbol{k}}+i\,Y_{\boldsymbol{k}})(\sigma_{0}+\eta\,\sigma_{x})$ and $(X_{\boldsymbol{k}}+i\,Y_{\boldsymbol{k}})\sigma_{0}+\eta\,(X_{\boldsymbol{k}}-i\,Y_{\boldsymbol{k}})\sigma_{x}$. While the first of the two states has two fully established gaps, given by $(1\pm\eta)\sqrt{X^{2}_{\boldsymbol{k}}+Y^{2}_{\boldsymbol{k}}}$ (with $\pm$ referring to the spin species), the other has two gaps, $|X_{\boldsymbol{k}}|$ and $|Y_{\boldsymbol{k}}|$, with distinct momentum dependencies; it, thus, exhibits two point nodes per Fermi surface which occur at different positions for the two spin species, similar to the associated triplet phase $E^{3_{s}}(f)$.

In general, admixing a singlet component at a second transition to a triplet state is less likely to occur as a singlet state has less options to “adapt” (the order parameter comprises two complex numbers for $E$) than a triplet state (for which, the order parameter comprises six complex numbers). While this is not possible for the one-dimensional representation $A$ (see Fig. 2), the IR $E$ does allow for this scenario but only for the triplet states $E^{3_{s}}(g)$ and $E^{3_{s}}(h)$: for small $\delta a<0$, we find a second transition where an additional chiral (nematic) singlet component is admixed to $E^{3_{s}}(g)$ ($E^{3_{s}}(h)$). As both pure triplet states can be fully gapped, the same holds for the admixed phases. The admixture of the extra singlet component does not change the symmetries of $E^{3_{s}}(g)$ and $E^{3_{s}}(h)$ listed in Sec. IV.1 above. Reversing the sign of $\delta a$ to small positive values, we obtain the same admixed phase. The only difference is that the first transition is a singlet transition into a chiral (nematic) phase and the secondary triplet $E^{3_{s}}(g)$ [$E^{3_{s}}(h)$] becomes nonzero at a lower transition temperature.

The key results of this section, the pure triplet/singlet states and the possible admixed phases for small $J$ along with their order parameters and properties, are summarized in Table 2. As already discussed above, several states are degenerate within single-band mean-field theory. Depending on the form of the corrections to mean-field theory lifting this degeneracy, there are two possible phase diagrams, shown in Fig. 5. Interestingly, we observe that the chiral singlet, $E^{1_{s}}(1,i)$, is not the only possible phase close to mean-field theory for anti-Hund’s coupling: as can be seen in Fig. 5(b), a secondary phase transition into the nematic mixed singlet-triplet state $E^{1_{s}}(0,i)+E^{3_{s}}(a)$ is predicted. It is a fully gapped state with an anisotropic gap, $\sqrt{\eta^{2}X_{\boldsymbol{k}}^{2}+Y_{\boldsymbol{k}}^{2}}$, breaking rotational symmetry. Note that this route to a nematic superconducting state, indications of which are provided by recent experiments Cao et al. (2020), is distinct from that of other works Kozii et al. (2019); Chichinadze et al. (2019). Of course, sufficiently large corrections to mean-field theory can in principle yield any of the phases listed in Table 2; we will come back to these corrections in Sec. V below.

Let us finally discuss the impact of fluctuations of the order parameter on the thermal phase transitions. As readily follows from the respective order parameter manifolds, the singlet phases in Table 2 exhibit a conventional BKT transition, the triplets (a), (b), (d), (f), (g), and (h) will be charge-$4e$ superconductors where only spin-rotation invariant combinations of the triplet vector assume quasi-long-range order at finite temperature, and the triplets (c) and (e) will only display a crossover. However, as pointed out above, none of these three classes of transitions can currently be excluded based on the experimental data.

IV.3 Behavior in a magnetic field

Finally, we turn our attention to the behavior of the pairing states of the complex representation in the presence of a Zeeman field, $\boldsymbol{M}_{Z}$, and in-plane orbital coupling $\boldsymbol{M}_{O}$, along the same lines as Sec. III.3. From Eqs. (31) and (37), it follows that there are three possible coupling terms linear in the field and quadratic in the superconducting order parameter given by

[TABLE]

Notice that, exactly as for the IR $A$, there is no linear coupling to the in-plane orbital field, which is prohibited by time-reversal and $C_{3}$ rotation symmetry. While the first term in Eq. (42) is again forced to vanish for $J\rightarrow 0$ [for the same reason as $\delta c_{1}$ in Eq. (24)], the second singlet-triplet-mixing coupling, $c_{2}^{E}$, is not constrained to be zero for $J=0$. However, the emergent symmetry in the single-band mean-field description of Appendix A, leads to $c_{2}^{E}=0$, so it is natural to expect $c_{2}^{E}\ll c_{3}^{E}$ such that the last term in Eq. (42) describes the dominant linear coupling to the magnetic field—even when $J$ is small. As expounded in Appendix A, the expression for $c_{3}^{E}$ is identical in form to that for $c_{2}$ in Eq. (24). As such, the linear increase of the (first) superconducting transition temperature with small magnetic fields seen in experiment does not permit one to distinguish between the IRs $A$ and $E$.

There is one difference between the pairing states of the two IRs worth mentioning here: while the form of the leading triplet vector in a magnetic field is completely fixed to be $\boldsymbol{d}\propto(1,i,0)^{T}$ for the one-dimensional IR $A$, the complex IR allows for either the nematic nonunitary $E^{3_{s}}(c)$ or the chiral nonunitary $E^{3_{s}}(e)$ pairing for nonzero $\boldsymbol{M}_{Z}$. Which of the two is realized, depends on the value of the quartic terms in Eq. (38): if $b_{2}^{t}+b_{3}^{t}>0$, the state $E^{3_{s}}(e)$ will be preferred while the opposite sign corresponds to $E^{3_{s}}(c)$. Within single-band mean-field theory, we find $b_{2}^{t}=0$ and $b_{3}^{t}>0$, which leads to phase $E^{3_{s}}(e)$. In the next section, we will see that additional ferromagnetic fluctuations will further enlarge the positive value of $b_{2}^{t}+b_{3}^{t}$ and consequently, not affect the mean-field prediction that $E^{3_{s}}(e)$ is the leading triplet state with the highest transition temperature in the presence of a magnetic field.

V Fluctuation-induced superconductivity

Among the plethora of possible superconducting phases outlined in this paper, only a few can be realized in single-band mean-field theory (see Tables 1, 2, and 4). This originates from the fact that, within single-band mean-field theory, the ratio of the quartic terms is fixed and only one state or two degenerate states can occur for each IR. However, the presence of sizable correlations in the nearly flat bands of graphene moiré systems is expected to give rise to significant corrections to mean-field theory. This has recently been demonstrated for the case of charge-density-wave fluctuations in twisted bilayer graphene Kozii et al. (2019), and in the context of nematic fluctuations in the iron-based superconductors Fernandes and Millis (2013). In this section, we study how corrections associated with ferromagnetic fluctuations will split the mean-field degeneracies and, if sufficiently strong, realize phases distinct from mean-field theory.

To this end, we will first focus on spin fluctuations. This is prompted by experiments Shen et al. (2019); Liu et al. (2019b); Cao et al. (2019), which indicate a spin-polarized correlated insulating state in twisted double-bilayer graphene, and by the fact that the superconducting phase emerges when doping out of this polarized state. Likewise, we also expect ferromagnetic fluctuations to play an important role in twisted bilayer Sharpe et al. (2019); Lu et al. (2019); Zondiner et al. (2019) and trilayer graphene Chen et al. (2020b). In particular, in the latter two systems, however, these fluctuations will likely not only be of spin but also of orbital origin. This is why we will also discuss orbital fluctuations.

As it is known to capture the essential physics Kozii et al. (2019); Fernandes and Millis (2013), we focus in the main text on a phenomenological Ginzburg-Landau-like approach (that does not explicitly take into account fluctuations with nonzero momentum and frequencies), but provide a systematic microscopic derivation in Appendix B.1. Representing the ferromagnetic spin moment in valley $v=\pm$ by $\boldsymbol{m}_{v}$, we parametrize its contribution to the free energy as

[TABLE]

In this expression, $\hat{\chi}$ plays the role of the spin susceptibility (with $|\delta\chi|<\chi$ to ensure stability) and we expect $\delta\chi>0$ close to a phase where the spin moments in the two valleys are aligned. The ratio $\delta\chi/\chi$ controls how strongly the SU(2)${}_{+}\times$ SU(2)- symmetry is broken down to SU(2)s.

V.1 Trivial representation

Focusing first on the one-dimensional IR $A$ of $C_{3}$, the magnetic moments couple to the superconducting order parameter in Sec. III according to

[TABLE]

where we have retained only the couplings invariant under SU(2)${}_{+}\times$ SU(2)- and assumed that $\delta\chi\neq 0$ in Eq. (43) is the main symmetry-breaking perturbation. Upon making the association $\boldsymbol{M}_{Z}=\sum_{v}\boldsymbol{m}_{v}$, we notice that $c_{2}$ is the same prefactor as in Eq. (24). In the same vein as Ref. Kozii et al., 2019, we integrate out the massive fluctuations of $\boldsymbol{m}_{v}$. As a consequence of the coupling (44), this yields corrections to the terms quartic in the superconducting order parameters in Eq. (20), which can be conveniently split into two categories. First, there are corrections that preserve the SU(2)${}_{+}\times$ SU(2)- symmetry; these can be restated as renormalizations of the coefficients $b_{1}$ and $b_{2}$ in Eq. (20). Corrections of the second type break this symmetry, violating the form of the free-energy expansion (20). More explicitly, the renormalization of the free energy $\mathcal{F}$ in Eq. (20) due to the presence of ferromagnetic spin fluctuations can be compactly stated as

[TABLE]

where $\delta_{1}$ $=$ $-\delta_{2}$ $=$ $2c_{2}^{2}(\chi-\delta\chi)$ $>$ [math] and $\delta_{3}$ $=$ $2c_{2}^{2}\delta\chi$. As required by symmetry, the contribution $\delta_{3}$ of the second category breaking the SU(2)${}_{+}\times$ SU(2)- symmetry is proportional to $\delta\chi$.

We start with the limit $|\delta\chi|\ll\chi$, where the structure of Eq. (20) is asymptotically preserved and the form of the two possible phase diagrams in Fig. 2 is unchanged. Since $\delta_{2}<0$, strong ferromagnetic fluctuations will change the sign of $b_{2}$ from its positive mean-field value to negative and, as opposed to mean-field theory, favor the phase diagram in part (b) of Fig. 2 over part (a). We point out that naively taking Eq. (45) alone would render the quartic free-energy expansion unstable for large enough $\chi$. However, denoting the mean-field value of $b_{2}$ by $b_{2}^{0}$, there exists a regime, $b_{2}^{0}/2<c_{s}^{2}\chi<b^{0}_{2}$, for which $b_{2}<0$ due to fluctuation corrections and the free energy in Eq. (20) is stable. For larger values of $\chi$, we can imagine adding the sextic term $c(\text{tr}[\Delta_{+}^{\dagger}\Delta_{+}^{\phantom{\dagger}}])^{3}$ to the free energy to restore stability.

When $\delta\chi$ is of order $\chi$, the ferromagnetic fluctuations described by Eq. (43) induce considerable SU(2)${}_{+}\times$ SU(2)--symmetry-breaking interactions. The presumed sign $\delta\chi>0$ brings about a further enhancement of the term $-|\boldsymbol{d}^{*}\times\boldsymbol{d}|^{2}$ [as is obvious from Eq. (45)], which favors nonunitary triplet pairing relative to the SU(2)${}_{+}\times$ SU(2)--invariant form of the free energy in Eq. (20). Given that $\delta\chi<\chi$, strong ferromagnetic fluctuations are still expected to change the sign of $b_{2}$ relative to mean-field theory. The additional effect of $\delta\chi$ lies in effecting an additional first-order transition to a nonunitary triplet state in a third transition at lower temperatures for anti-Hund’s coupling in Fig. 2(b).

We have thus shown that significant ferromagnetic fluctuations can reverse the predictions of mean-field theory, and favor the nonunitary triplet state $A^{3_{s}}(1,i,0)$ and the admixed singlet-triplet phase $A^{1_{s}}+A^{3_{s}}(1,0,0)$ in Table 1.

V.2 Complex representation

The same analysis can be performed for the complex IR $E$ of Sec. IV. In this case, the most general SU(2)${}_{+}\times$ SU(2)--invariant coupling between the superconducting order parameter and the spin fluctuations allows for two independent coupling constants, $c_{\pm}$ $\in$ $\mathbb{R}$, and has the form

[TABLE]

Integrating out $\boldsymbol{m}_{v}$, we again obtain corrections to the free energy which are quartic in the superconducting order parameter. In the limit of SU(2)${}_{+}\times$ SU(2)- invariance, $\delta\chi=0$, these corrections can be represented by renormalizations of the couplings, $b_{j}\rightarrow b_{j}+\delta b_{j}$, in Eq. (41) with

[TABLE]

To study the ramifications of this result, we first consider the limit of weak fluctuations, for which $\delta b_{j}$ in Eq. (47) are much smaller in magnitude than the mean-field value of $b_{2}$. Albeit small, the corrections $\delta b_{j}$ are crucial here due to the exact degeneracy of the states $E^{3_{s}}(b)$ and $E^{3_{s}}(d)$ in mean-field theory observed earlier. From Eq. (41) with the replacement $b_{j}\rightarrow b_{j}+\delta b_{j}$, we find the free-energy difference of these two states to be

[TABLE]

thereby generically favoring $E^{3_{s}}(b)$ along with its Hund’s partner $E^{1_{s}}(0,i)+E^{3_{s}}(a)$, defined in Table 2; in other words, the phase diagram in Fig. 5(b) is favored over that in part (a). In the one-band description of Appendix B.2, it always holds that $c_{+}=c_{-}$, which is, in turn, a consequence of an emergent valley-exchange symmetry. However, multiband effects are expected to be present Xie et al. (2019) and to lead to nonzero $|c_{+}-c_{-}|\ll|c_{+}|$, which is enough to lift the degeneracy according to Eq. (48).

Next, we turn to the limit of strong ferromagnetic fluctuations, where the mean-field values of $b_{j}$ have to be treated as perturbations to the large $\delta b_{j}$ in Eq. (47). As $\chi\rightarrow\infty$, we find that, out of the triplet states in Table 2, $E^{3_{s}}(e)$ has the lowest energy unless $c_{+}=c_{-}$ or $c_{+}=-c_{-}$. We know that $c_{+}\simeq c_{-}$ and hence, can safely neglect the latter. For the former option, $E^{3_{s}}(e)$ is found to be degenerate with $E^{3_{s}}(c)$; however, for large but finite $\chi$, the additional contribution to $b_{j}$ from mean-field theory lifts this degeneracy, always selecting $E^{3_{s}}(e)$. Out of the multitude of possible pairing states in Table 2, strong ferromagnetic fluctuations thus favor the chiral nonunitary triplet state $E^{3_{s}}(e)$ and the mixed singlet-triplet phase $E^{1_{s}}(1,i)+E^{3_{s}}(d)$. Which of these two states is realized, depends on whether singlet or triplet has the higher transition temperature (the sign of $\delta a$).

V.3 Orbital fluctuations

Anticipating its relevance for twisted bilayer and trilayer graphene, here, we extend the previous analysis to the case of orbital ferromagnetic fluctuations. Due to the two-dimensional nature of the system, the in-plane orbital moments, $\boldsymbol{M}_{O}=(M_{O}^{x},M_{O}^{y})$, and the out-of-plane moment $M_{O}^{z}$ behave quite differently. Beginning with the complex representation, we already know from Sec. IV.3 that there is no linear coupling to $\boldsymbol{M}_{O}$; however, the superconductor can couple to $M_{O}^{z}$ as

[TABLE]

For concreteness, one might think of $M_{O}^{z}$ as valley fluctuations, associated with $\sum_{\boldsymbol{k}}c_{\boldsymbol{k}}^{\dagger}\tau_{z}c_{\boldsymbol{k}}^{\phantom{\dagger}}$, but our analysis is more general. Taking an energetic contribution quadratic in $M_{O}^{z}$ similar to Eq. (43) and integrating over $M_{O}^{z}$, we obtain a correction to the free-energy that can be conveniently expressed as

[TABLE]

in Eq. (41). It is easily seen that taking this as a small correction to mean-field theory will now favor the phase diagram in Fig. 5(a) over that in part (b). On the other hand, in the limit of strong orbital fluctuations, the chiral unitary, $E^{3_{s}}(d)$, and the chiral nonunitary triplet, $E^{3_{s}}(e)$, (along with their Hund’s partners) will be favored. This degeneracy will be lifted by the subleading ferromagnetic spin fluctuations, which favor the latter state, $E^{3_{s}}(e)$ (and its Hund’s partner), as readily follows from Eq. (47).

In the trivial representation, orbital fluctuations have no impact on which of the two possible phase diagrams in Fig. 2 is realized. This results from the fact that neither $\boldsymbol{M}_{O}$ (see Sec. III.3) nor $M_{O}^{z}$ can couple linearly to the superconducting states and their rotational invariant quadratic forms $(M_{O}^{z})^{2}$, $\boldsymbol{M}_{O}^{2}$ can only couple to $|\Delta^{s}|^{2}+\boldsymbol{d}^{\dagger}\boldsymbol{d}$. Consequently, the energetic correction obtained by integrating out the orbital fluctuations will also only depend via $|\Delta^{s}|^{2}+\boldsymbol{d}^{\dagger}\boldsymbol{d}$ on the superconducting states and, as such, not affect the value of $b_{2}$ in Eq. (20) and Fig. 2.

V.4 In a magnetic field

Finally, we come back to the impact of fluctuation corrections on the leading triplet phase in the presence of a magnetic field. As we have seen in Sec. IV.3, the superconducting state with the highest transition temperature in the presence of a sufficiently strong magnetic field will be a triplet phase due to the linear coupling in the second line of Eq. (42). At the mean-field level, $b_{2}^{t}+b_{3}^{t}>0$, which prefers $E^{3_{s}}(e)$ over $E^{3_{s}}(c)$ as the order parameter of this phase. Using the relations in Eq. (83), it is straightforward to rephrase the fluctuation corrections (47) and (50) of $b_{j}$ in terms of $b_{j}^{t}\rightarrow b_{j}^{t}+\delta b_{j}^{t}$ in Eq. (38). This yields $\delta b^{t}_{1}+\delta b^{t}_{2}=\chi(c_{+}-c_{-})^{2}>0$ and $\delta b^{t}_{1}+\delta b^{t}_{2}=4\delta b>0$ for spin and orbital fluctuations, respectively. We conclude that, as expected, ferromagnetic fluctuations do not change the mean-field prediction in this case and $E^{3_{s}}(e)$ is the dominant triplet order parameter in the presence of a magnetic field, for both strong and weak ferromagnetic fluctuations, and in their absence.

VI Adding further symmetries

In this section, we will analyze how the results presented above are modified once the additional symmetries, twofold rotation, $C_{2}$, perpendicular to the plane of the system, and in-plane rotation symmetry, $C_{2y}$, are added. As shown in Fig. 1(b) and (c), these symmetries are relevant as either exact microscopic or approximate emergent symmetries of twisted bilayer graphene and ABC trilayer graphene on hexagonal boron nitride, both of which exhibit superconductivity Cao et al. (2018b); Chen et al. (2019c).

VI.1 Consequences of a $C_{2}$ rotation symmetry

One crucial difference in twisted bilayer compared to twisted double-bilayer graphene is that the former has an approximate $C_{2}$ symmetry Po et al. (2018) that mixes the two valleys, i.e., the system is (approximately) invariant under

[TABLE]

To relate to our notation used above, we assume that it is sufficient to focus on a single band for describing superconductivity in twisted bilayer graphene as well. This is quite a natural assumption and, unless stated otherwise, we expect our conclusions to hold when additional bands are taken into consideration.

This (approximate) symmetry has attracted a lot of attention in the recent theory literature Kang and Vafek (2018); Yuan and Fu (2018); Koshino et al. (2018); Zou et al. (2018) of the system since it, combined with time-reversal and $C_{3}$, leads to a $C_{6}\Theta$ symmetry, which is responsible for not only the presence of (nearly gapless) Dirac cones at $K$ and $K^{\prime}$ but also the (approximate) vanishing of Berry curvature in twisted bilayer graphene. If the twist axis goes through the center of a hexagon, the system has $C_{6}$ rotation even as a microscopic symmetry. We note in passing that the (nearly) flat bands obtained in Refs. Lee et al., 2019; Koshino, 2019 for double-bilayer graphene do not feature any Dirac cones but have well-separated conduction and valence bands that are characterized by nonzero Chern numbers (at least in some parameter regime); this strongly indicates that $C_{2}$ is not an approximate symmetry in twisted double-bilayer graphene since $C_{2}\Theta$ would enforce zero Berry curvature.

In a similar fashion, Ref. Zhang and Senthil, 2019 has argued that the twofold symmetry (51) is also an approximate symmetry for ABC trilayer graphene on hexagonal boron nitride, although it is clearly not a microscopic symmetry of the system, as can be seen in Fig. 1(c).

All things considered, it is currently not known whether an approximate $C_{2}$ symmetry is relevant for superconductivity in twisted bilayer and ABC trilayer graphene. Therefore, we will now discuss what changes for the possible superconducting instabilities once we assume that the Hamiltonian is also invariant under the transformation in Eq. (51).

The $C_{2}$ transformation plays a special role in two dimensions as it is equivalent to $\boldsymbol{k}\rightarrow-\boldsymbol{k}$ and can, thus, significantly affect superconducting instabilities Scheurer et al. (2017). In graphene moiré superlattices, it also relates the two valleys and “interferes” with the Fermi-Dirac constraint (10): decomposing the pairing into singlet and triplet,

[TABLE]

Eq. (10) implies that $\lambda^{s}_{\boldsymbol{k}v}$ $=$ $\lambda^{s}_{-\boldsymbol{k}\bar{v}}$ and $\lambda^{t}_{\boldsymbol{k}v}$ $=$ $-\lambda^{t}_{-\boldsymbol{k}\bar{v}}$. Consequently, it holds (as long as the pairing matrix elements between different bands can be neglected) that

[TABLE]

i.e., all representations even (odd) in $C_{2}$ must be pure singlet (triplet) states and vice versa. This has a few implications worth mentioning. First, even if $C_{2}$ is not a good symmetry (say, it is significantly broken by interactions), SU(2)s spin-rotation invariance requires that the first transition must be into a pure singlet or triplet state and hence, the pairing must be either even or odd under $C_{2}$. In this sense, we can still distinguish between $p$-wave and $d$-wave pairing despite the presence of $C_{2}$-symmetry-breaking interactions. We emphasize that mixing will only be possible via multiple superconducting transitions (associated with admixtures of singlet and triplet) or interband pairing. The latter is expected to be quite weak given that the typical splitting between the bands at half-filling (at least a few meV Kerelsky et al. (2019)) is about or more than an order of magnitude larger than the superconducting critical temperature ($\simeq$ $0.15\,\textrm{meV}$ according to Ref. Cao et al., 2018b).

Secondly, if we do have an enhanced $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ symmetry (or are close to it), singlet and triplet are (nearly) degenerate. This forces the corresponding IRs of the spatial point group $D_{6}$ of the system, which behave identically under the subgroup $D_{3}$ but are even and odd under $C_{2}$, to be (nearly) degenerate at the quadratic level of the Ginzburg-Landau expansion. For instance, $A_{1}$ and $B_{1}$ of $D_{6}$ have to be degenerate, as summarized in Table 3. Without a Zeeman field, an extra $C_{2}$ symmetry with action in Eq. (53) also has no consequences for the higher-order terms in the free energy since spin-rotation invariance necessitates that all of these terms are even in the triplet vector. The only difference arises in the presence of a Zeeman field or magnetic fluctuations: with $C_{2}$ symmetry, it must hold that $\delta c_{1}=0$ in Eq. (24) and $\delta c_{1}^{E}=0$ in Eq. (42) even when the SU(2)${}_{+}\times$ SU(2)- symmetry is broken. Furthermore, a $C_{2}$ symmetry implies $c_{+}=c_{-}$ in Eq. (46). For this reason, weak ferromagnetic spin fluctuations do not lift the degeneracy of mean-field theory if we impose an exact $C_{2}$ symmetry and other types of fluctuations have to be considered. Recall, however, in both trilayer and twisted bilayer graphene, $C_{2}$ should only be considered as an approximate symmetry and terms breaking this symmetry will lead to $c_{+}\neq c_{-}$, thus lifting the degeneracy of mean-field theory by, e.g., favoring Fig. 5(b) over (a). Note that the $C_{2}$ symmetry also forces $c_{O}^{E}$ in Eq. (49) to vanish for $M_{O}^{z}$ corresponding to valley fluctuations. As such, the approximate $C_{2}$ symmetry does not specify whether spin or valley fluctuations are expected to be the dominant source of lifting the mean-field degeneracy. It only indicates that strong ferromagnetic fluctuations are most likely dominated by spin fluctuations.

In summary, when classifying superconducting states in twisted bilayer graphene or ABC trilayer graphene on hexagonal boron nitride in the absence of a Zeeman field, it is unimportant whether an approximate $C_{2}$ symmetry is relevant or not: singlet and triplet will always be even and odd under it. We can thus work with $D_{3}$ (instead of $D_{6}$) without loss of generality in the following. The only difference with twisted double-bilayer graphene (with finite displacement field) is an extra twofold rotation symmetry, $C_{2y}$, along the $y$-axes, see Fig. 1(b) and (c). Its action on the electronic operators reads as

[TABLE]

where $C_{2y}\boldsymbol{k}=(-k_{x},k_{y})$. The upshot of this additional symmetry for the possible superconducting instabilities is clarified in the next subsection.

VI.2 $D_{3}$ versus $C_{3}$

Due to the additional $C_{2y}$ symmetry, $D_{3}$ is a non-Abelian group and has three IRs—two one-dimensional and one two-dimensional representation (refer to the character table in Table 3). It is convenient to begin with the one-dimensional IRs $A_{1}$ and $A_{2}$ and take $J\neq 0$. Since $C_{2y}$ interchanges the valleys, its action on the intervalley pairing order parameter (52) can be written as

[TABLE]

So, we see that a singlet (triplet) state transforming under $A_{1}$ ($A_{2}$) has no nodes while a singlet (triplet) in the $A_{2}$ ($A_{1}$) channel has symmetry-imposed nodes on the line $k_{y}=0$ and along the directions rotated by $\pm\pi/3$. This creates six nodal points on any surface enclosing the $\Gamma$ point.

We can also readily understand from Eq. (55) how the one-dimensional representations “connect” at the SU(2)+ $\times$ SU(2)- point: at the high-symmetry point, $\lambda^{s}_{\boldsymbol{k}v}=\lambda^{t}_{\boldsymbol{k}v}$, ergo $A_{1}^{1_{s}}$ and $A_{2}^{3_{s}}$ or $A_{2}^{1_{s}}$ and $A_{1}^{3_{s}}$ must meet at the $J=0$ line in Fig. 2. We summarize these observations in Table 4.

In addition, for the case of the two-dimensional representation $E$ of $D_{3}$, the $C_{2y}$ symmetry has nontrivial consequences. Once again, we take $J\neq 0$ which permits us to study singlet and triplet independently. As singlet pairing has already been analyzed in detail for twisted bilayer graphene (see, e.g., Ref. Venderbos and Fernandes, 2018), we are chiefly concerned with the triplet states here. We parametrize the triplet pairing as in Sec. IV.1 with the sole distinction being that the basis functions $X_{\boldsymbol{k}}$ and $Y_{\boldsymbol{k}}$ are now constrained by the symmetries of $D_{3}$; we choose them to obey $X_{-C_{2y}\boldsymbol{k}}=-X_{\boldsymbol{k}}$ and $Y_{-C_{2y}\boldsymbol{k}}=Y_{\boldsymbol{k}}$, while transforming as $k_{x}$ and $k_{y}$ under $C_{3}$. A possible choice is given by Eq. (30) with $\phi=\pi/2$ for the Brillouin zone of twisted bilayer graphene in Fig. 1(b). With these conventions, the triplet vector transforms according to $(\boldsymbol{d}_{+},\boldsymbol{d}_{-})\rightarrow(\boldsymbol{d}_{-},\boldsymbol{d}_{+})$ under $C_{2y}$. This does not further constrain the quartic terms in the free energy (38), wherefore we can use the analysis of Sec. IV.1 for the point group $C_{3}$, bearing in mind the caveat that the relative phase, $\varphi$, between $\boldsymbol{d}_{+}$ and $\boldsymbol{d}_{-}$ cannot be absorbed in a redefinition of the basis functions $X_{\boldsymbol{k}}$ and $Y_{\boldsymbol{k}}$ any more due to the extra $C_{2y}$ symmetry. While $\varphi$ has no consequences for $E^{3_{s}}(d)$ or $E^{3_{s}}(e)$ and can be absorbed by performing a spin rotation for the phases $E^{3_{s}}(b)$ and $E^{3_{s}}(g)$, it describes different phases for all other stable minima of Eq. (38), and we have to go to higher order in the free-energy expansion to determine its value.

Consider $E^{3_{s}}(a)$ for instance. Writing $\boldsymbol{d}_{+}=(1,0,0)^{T}$ and $\boldsymbol{d}_{-}=e^{i\varphi}(1,0,0)^{T}$, it is easy to verify that the most general, $\varphi$-dependent sextic term to the free energy must have the form $c_{1}\cos(3\varphi)$ with $c_{1}\in\mathbb{R}$. This derives from Eq. (34) where the $C_{2y}$ symmetry forces $c_{2}$ to vanish. We thus find $\varphi=2\pi n/3$, $n\in\mathbb{Z}$, for $c_{1}<0$ and $\varphi=\pi/3+2\pi n/3$ when $c_{1}>0$. These two minima correspond to two different states, which can be compactly represented by defining the “rotated” basis functions

[TABLE]

with $X_{\boldsymbol{k}}$ and $Y_{\boldsymbol{k}}$ as introduced above. The order parameters are $M_{\boldsymbol{k}+}$ $=$ $X^{0}_{\boldsymbol{k}}\sigma_{x}$ and $M_{\boldsymbol{k}+}$ $=$ $X^{\frac{\pi}{3}}_{\boldsymbol{k}}\sigma_{x}$ $\equiv$ $(\sqrt{3}X_{\boldsymbol{k}}+Y_{\boldsymbol{k}})\sigma_{x}/2$ for $c_{1}<0$ and $c_{1}>0$, respectively. We denote these two states by $E^{3_{s}}(a)_{0}$ and $E^{3_{s}}(a)_{\frac{\pi}{3}}$, respectively. The first state, $E^{3_{s}}(a)_{0}$, preserves $C_{2y}$, but breaks $C_{3}$ rotation symmetry, and has a nodal line which is, as opposed to the states in Sec. IV.1, pinned to $k_{y}=0$. The other state, $E^{3_{s}}(a)_{\frac{\pi}{3}}$, however, breaks $C_{2y}$ and the nodal line is not pinned to the $k_{x}$ axis.

The remaining triplet states, $E^{3_{s}}(c)$, $E^{3_{s}}(f)$, and $E^{3_{s}}(h)$ of Sec. IV.1 can be analyzed in the same way. In all cases, we find two states corresponding to two different discrete values of the relative phase $\varphi$ between $\boldsymbol{d}_{+}$ and $\boldsymbol{d}_{-}$: for $E^{3_{s}}(c)$ and $E^{3_{s}}(h)$, we find $\varphi=0$ or $\varphi=\pi/3$ as before, whereas $E^{3_{s}}(f)$ requires even higher-order terms in the free energy expansion, yielding $\varphi=0$ or $\varphi=\pi/2$. In analogy to $E^{3_{s}}(a)_{\varphi}$, we label the states by $E^{3_{s}}(c)_{\varphi}$, $E^{3_{s}}(f)_{\varphi}$, $E^{3_{s}}(h)_{\varphi}$; their order parameters are the same as those of the corresponding states in Sec. IV.1 but with the rotated basis functions in Eq. (56) using the respective value of $\varphi$. Taken together, we obtain twelve triplet states for $D_{3}$, which are summarized in Table 4, instead of only eight for the point group $C_{3}$.

Finally, we can also ask how the different states behave for small $J$, i.e., whether singlet and triplet can mix and which phases are Hund’s partners. Exactly as illustrated above for the pure triplet phases, we have to consider higher-order terms that determine the relative phase between the chiral, $\mu=+$, and antichiral, $\mu=-$, basis functions. As this analysis closely parallels our previous discussions, we just present the result in Table 4. In total, there are ten symmetry-inequivalent mixed singlet and triplet phases. Seven of them are only possible if $\delta a<0$ (singlet dominates); the remaining three can be realized for either sign of $\delta a$.

VII Discussion and conclusion

In this work, we have presented a systematic classification and analysis of superconducting instabilities in graphene moiré systems. To this end, we have focused on zero-momentum Cooper pairs formed out of electrons in different valleys. Intervalley pairing is expected to be the dominant pairing channel as time-reversal relates the two valleys. We have first analyzed singlet and triplet pairing separately since spin-orbit coupling is expected to be very weak in graphene. However, theoretical estimates of the interaction terms of twisted bilayer Koshino et al. (2018), double-bilayer Lee et al. (2019), and trilayer Zhang and Senthil (2019) graphene indicate that these systems are approximately invariant under independent spin rotations in the two valleys, leading to an (approximate) SU(2)+ $\times$ SU(2)- symmetry and the (near) degeneracy of singlet and triplet pairing. For this reason, we have also classified the pairing instabilities close to this high-symmetry point, analyzing which triplet state transforms into which singlet phase upon changing the sign of the interactions breaking the SU(2)+ $\times$ SU(2)- symmetry. We have further derived the conditions under which singlet and triplet can mix despite the absence of spin-orbit coupling.

As it has the fewest symmetries, we first considered twisted double-bilayer graphene, for which there are also clear experimental indications of triplet pairing Shen et al. (2019); Liu et al. (2019b). Here, a displacement field, which is required to stabilize the superconducting state, reduces the point group to $C_{3}$. The pairing states and their properties associated with the real representation $A$ and the complex representation $E$ of $C_{3}$ are summarized in Tables 1 and 2, respectively.

Being one-dimensional and real, $A$ only allows for one singlet, a unitary and a nonunitary triplet phase, and one mixed phase. The latter is expected to be relevant only if SU(2)+ $\times$ SU(2)- is weakly broken and the two consecutive transitions in the schematic phase diagram of Fig. 2(b) are very close. Using the values of the coupling constants in Ref. Lee et al., 2019, we estimate the splitting to be about two orders of magnitude smaller than the critical temperature and hence, hard to see experimentally Note (1). Whether renormalization-group corrections could enhance the impact of these weak symmetry-breaking perturbations at energies of order of the transition temperature is an open question, which we leave for future work. The gap structure of the four phases transforming under $A$ is quite different: while the nonunitary triplet is gapless for one of the spin species, the singlet and unitary triplet have a single, fully established gap, and the mixed phase has two finite but distinct gaps for the two spin species. We have further shown that single-band mean-field theory will generically favor the phase diagram in Fig. 2(a) over Fig. 2(b). However, the small bandwidth and strong-coupling nature inherent in the problem makes the applicability of mean-field theory questionable and can lead to significant corrections which might eventually select other phases. We have illustrated these corrections for ferromagnetic fluctuations, expected to be relevant for twisted double-bilayer graphene Shen et al. (2019); Liu et al. (2019b); Cao et al. (2019), twisted bilayer Sharpe et al. (2019); Lu et al. (2019); Zondiner et al. (2019), and ABC trilayer graphene Chen et al. (2020b). We find that the resulting corrections will, as opposed to mean field, generally favor the phase diagram in part (b) of Fig. 2 over that in part (a).

The complex representation allows for many more states: two pure singlets, eight triplets, and, if SU(2)+ $\times$ SU(2)- is only weakly broken, six distinct mixed phases. As compiled in Table 2, all of these three classes of states allow for nodal points and fully gapped phases. However, only the triplets can have nodal lines (residual ungapped Fermi surfaces of one spin species). Only one out of the two different triplet states of the IR $A$ allow for an admixture of singlet and triplet for weak anti-Hund’s coupling but, in contrast, six out of the eight triplets transforming under $E$ do so.

Out of the possible pairing states in Table 2, single-band mean-field theory favors the two triplet states $E^{3_{s}}(b)$ and $E^{3_{s}}(d)$ along with their respective Hund’s partners—the nematic mixed phase $E^{1_{s}}(0,i)+E^{3_{s}}(a)$ and the chiral singlet $E^{1_{s}}(1,i)$. We show the associated phase diagrams in the vicinity of mean-field theory in Fig. 5(a) and (b). We have discussed how additional weak ferromagnetic spin (orbital) fluctuations can lift the exact degeneracy of $E^{3_{s}}(b)$ and $E^{3_{s}}(d)$, generically favoring the former (latter) and, hence, the phase diagram in Fig. 5(b) [Fig. 5(a)]. In the limit of strong ferromagnetic fluctuations, we obtain the chiral nonunitary triplet $E^{3_{s}}(e)$ or, for weak anti-Hund’s coupling, the mixed singlet-triplet state $E^{1_{s}}(1,i)+E^{3_{s}}(d)$ as the dominant instability.

Motivated by the experimentally observed Liu et al. (2019b) linear increase of the transition temperature with an in-plane magnetic field in twisted double-bilayer graphene and signs of magnetism in bilayer and trilayer graphene, we have also mapped out the possible phase diagrams in the presence of a magnetic field. As expected, if the SU(2)+ $\times$ SU(2)- symmetry is significantly broken, the linear increase is only consistent with triplet pairing. For pairing in the $A$ channel, there are two possible phase diagrams, shown in Fig. 3(c) and (d), depending on which triplet state is realized in the absence of a magnetic field. The magnetic field fully determines the form of the leading triplet state to be $A^{3_{s}}(1,i,0)$ in the $A$ channel. For order parameters transforming under $E$, there are two possibilities for the leading triplet state, $E^{3_{s}}(c)$ or $E^{3_{s}}(e)$, in a magnetic field; which of the two is realized depends on the value of the quartic couplings in the free energy. Both mean-field theory and ferromagnetic fluctuations favor the $E^{3_{s}}(e)$ state. If, however, SU(2)+ $\times$ SU(2)- is only very weakly broken, singlet pairing as the dominant instability of the system is also consistent with the linear increase of the critical temperature; the two possible phase diagrams for the case of pairing in the IR $A$ are illustrated in Fig. 3(a) and (b).

We have also derived (within mean-field theory) the key couplings, $c_{2}$ in Eq. (24) and $c_{3}^{E}$ in Eq. (42), between the superconducting order parameter and the magnetic field $B$ that determine the slope of the increase, $\Delta T_{c}$, of the critical temperature with magnetic field. We found that they have the exact same mathematical form; as such, the behavior $\Delta T_{c}\simeq 2\mu_{B}B$, with Bohr magneton $\mu_{B}$, seen in experiment Liu et al. (2019b), is equally surprising for both pairing channels and does not favor one channel over the other. In both cases, this might either be accidental or due to quantum critical scaling Lee et al. (2019).

We have also studied, in Sec. VI, the changes in the classification when there is an extra in-plane rotation symmetry, $C_{2y}$, and a twofold rotation, $C_{2}$, perpendicular to the plane. These two symmetries are relevant (either as exact or emergent symmetries) to twisted bilayer graphene and ABC trilayer graphene. We find that while the $C_{2}$ symmetry has no consequences for the classification, $C_{2y}$ not only pins the nodes of certain pairing states along high-symmetry lines but also leads to more pairing states as summarized in Table 4.

This work further illustrates that graphene moiré systems provide a very rich playground for novel strongly correlated superconducting phases. We hope that our systematic analysis of pairing in the absence and presence of magnetic fields will help future theoretical and experimental studies to pinpoint the microscopic form of the superconducting state.

Acknowledgments

This research was supported by the National Science Foundation under Grant No. DMR-1664842. MS also acknowledges support from the German National Academy of Sciences Leopoldina through grant LPDS 2016-12. We thank Andrey Chubukov, Eslam Khalaf, Alex Kruchkov, Subir Sachdev, Harley Scammell, and Ashvin Vishwanath for helpful discussions.

Appendix A Microscopic Ginzburg-Landau expansion

In this appendix, we derive the prefactors of the various free-energy expansions in the main text within mean-field theory. Unless stated otherwise, we use a single-band description.

A.1 Without a magnetic field

We imagine performing a mean-field decomposition in the Cooper channel and keeping only the singlet and triplet pairing of the dominant IR. The ensuing mean-field Hamiltonian for the one-band model has the form

[TABLE]

where $\xi_{\boldsymbol{k}+}=\xi_{-\boldsymbol{k}-}$ due to time-reversal symmetry. In Eq. (57), we have omitted a constant term, which is quadratic in the superconducting order parameter and does not affect the quartic terms we derive below. Upon integrating out the fermions in Eq. (57) and expanding the resulting free energy in the superconducting order parameter, the Ginzburg-Landau expansion coefficients can be obtained order by order.

Starting with the one-dimensional real IR $A$ of $C_{3}$, we write $\Delta^{s}_{\boldsymbol{k}}=\lambda_{\boldsymbol{k}}^{s}\Delta^{s}$, $\boldsymbol{d}^{s}_{\boldsymbol{k}}=\lambda_{\boldsymbol{k}}^{t}\boldsymbol{d}$, where $\lambda_{\boldsymbol{k}}^{s}$ and $\lambda_{\boldsymbol{k}}^{t}$ are momentum-dependent basis functions that are invariant under $C_{3}$. Using the generalization of Eq. (20) to parametrize the free energy,

[TABLE]

which allows us to account for a nonzero $J$ making singlet and triplet nonequivalent, we find

[TABLE]

To keep the expressions compact, we have defined the functional

[TABLE]

When $J=0$, we have $\lambda_{\boldsymbol{k}}^{s}=\lambda_{\boldsymbol{k}}^{t}$ and hence, obtain

[TABLE]

which is compatible with the prefactors in Eq. (20) as required from the SU(2)+ $\times$ SU(2)- symmetry. On top, $\gamma_{1}=\gamma_{3}$ is an additional constraint arising from the mean-field approximation (and not related to an exact symmetry). In terms of the prefactors in Eq. (20), it sets $b_{1}=0$, as stated in the main text. The positive sign of the coefficients in Eq. (62) implies that mean-field theory always favors part (a) in the phase diagram in Fig. 2.

Similarly, we can study the complex representation of $C_{3}$ introduced in Sec. IV of the main text. Using the representation in Eq. (29) for the singlet pairing, $\Delta^{s}_{\boldsymbol{k}}=\sum_{\mu}\eta_{\mu}\left(X_{\boldsymbol{k}}+i\mu\,Y_{\boldsymbol{k}}\right)$, it is straightforward to show that

[TABLE]

for the coefficients $b^{s}_{1,2}$ in Eq. (32). Being positive, these coefficients favor the chiral superconductor $E^{1_{s}}(1,i)$ as was observed earlier as well You and Vishwanath (2019); Kozii et al. (2019).

Finally, repeating this procedure for the triplet state with parametrization (36), $\boldsymbol{d}_{\boldsymbol{k}}=\sum_{\mu}\boldsymbol{d}_{\mu}\left(X_{\boldsymbol{k}}+i\mu\,Y_{\boldsymbol{k}}\right)$, the coefficients in Eq. (38) evaluate to

[TABLE]

The triplet states $E^{3_{s}}(b)$ and $E^{3_{s}}(d)$ will have the lowest energy for this configuration of quartic coefficients as argued in the main text. The degeneracy between these two states is lifted by corrections beyond the mean-field approximation, such as the ferromagnetic fluctuations of Sec. V. In the presence of a magnetic field, Eq. (64) uniquely determines the chiral nonunitary triplet $E^{3_{s}}(e)$ as the leading instability (see Sec. IV.3).

A.2 Coupling to a magnetic field

In this subsection, we will analyze several important coupling terms between the superconductor and the magnetic field from a weak-coupling perspective. The microscopic form of the coupling to the Zeeman, $\boldsymbol{M}_{Z}$, and in-plane orbital field, $\boldsymbol{M}_{O}$, reads as

[TABLE]

where we have absorbed the $g$-factor of the Zeeman coupling into the definition of $\boldsymbol{M}_{Z}$. This is not possible for the orbital coupling, as its $g$-factor $\boldsymbol{g}_{v}(\boldsymbol{k})$ depends significantly on momentum. The form of $\boldsymbol{g}_{v}(\boldsymbol{k})$ is determined by microscopic details such as the Bloch states. All we need here is that $\boldsymbol{g}_{v}(\boldsymbol{k})=-\boldsymbol{g}_{\bar{v}}(-\boldsymbol{k})$, as follows from time-reversal symmetry (5), and we refer to Ref. Lee et al., 2019 for a microscopic derivation of its momentum dependence.

Let us first note that even when the actual interacting multiband system is not invariant under $C_{2}$, the single-band mean-field Hamiltonian, $\mathcal{H}_{\text{MF}}+\mathcal{H}_{B}$, is left invariant under the action of $C_{2}$ in Eq. (51) if we further set $\boldsymbol{d}_{\boldsymbol{k}}\rightarrow-\boldsymbol{d}_{\boldsymbol{k}}$ in Eq. (57) and $\boldsymbol{M}_{O}\rightarrow-\boldsymbol{M}_{O}$. This emergent symmetry is a consequence of the special role of $C_{2}$ in two dimensions as it acts on $\boldsymbol{k}$ in the same manner as time-reversal and, as such, can have crucial consequences for superconducting pairing Scheurer et al. (2017).

In the present case, this symmetry implies that the coupling terms $\delta c_{1}$ in Eq. (24) and $\delta c_{1}^{E}$, $c_{2}^{E}$ in Eq. (42) will vanish within single-band mean-field theory as is also readily confirmed by explicit calculation; we emphasize, however, that this is not an exact statement and we have checked that a multiband mean-field description allows for nonzero values. Nonetheless, we view the vanishing of these coupling in the weak-coupling single-band limit as an indication that they are likely small in the system.

Finally, the couplings of the Zeeman term to the triplet vector in Eqs. (24) and (42) are also not constrained by the emergent $C_{2}$ symmetry. We find these to be nonzero and given by

[TABLE]

respectively. Our main observation here is that the forms of $c_{2}$ and $c_{3}^{E}$ are identical: the nonuniversal part is a momentum integral which, in both cases, is weighted by a function that is invariant under $C_{3}$ and has no symmetry-imposed nodes on the Fermi surface. Accordingly, it is not possible to distinguish between the IRs $A$ and $E$ based on the slope of the increase of $T^{+}_{c}$ in small magnetic fields.

Appendix B Fluctuation corrections to mean-field

In this appendix, we provide further details on Sec. V.

B.1 Microscopic derivation

In this first part, we will derive, from a microscopic description of the system, that the prediction of the phenomenological approach of the main text provides the leading correction to the free energy of the superconductor in the limit where the mass of the fluctuations approaches zero.

To this end, we will use the field-theoretical formalism and describe the system by the action

[TABLE]

which consists of several contributions: first,

[TABLE]

is the free-electron contribution (with Grassmann fields $c$ and $c^{\dagger}$, in analogy to the operators in the main text), where $\int_{k}\dots\equiv T\sum_{\omega_{n}}\sum_{\boldsymbol{k}}\dots$ with fermionic Matsubara frequencies $\omega_{n}=\pi T(2n+1)$, and $k=(\boldsymbol{k},\omega_{n})$ comprising momentum and frequency. The second term,

[TABLE]

describes pairing, similar to the second line of Eq. (57), where we omitted the term proportional to the order parameter squared, since it is irrelevant for the free-energy contribution at quartic order in the superconducting state. The ferromagnetic fluctuations in valley $v$ with associated bosonic fields $\boldsymbol{\phi}_{qv}=(\phi^{x}_{qv},\phi^{y}_{qv},\phi^{z}_{qv})^{T}$ are described by the action

[TABLE]

where $q\equiv(\boldsymbol{q},\Omega_{n})$ is the bosonic analogue of $k$, i.e., $\Omega_{n}$ $=$ $2\pi Tn$ are bosonic Matsubara frequencies. Physically, $\hat{\chi}(\boldsymbol{q},i\Omega_{n})$ plays the role of the (analytic continuation of the) dynamical spin susceptibility [as compared to the static one in Eq. (43) of the main text]. We will focus here on the $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$-symmetric limit, where $[\hat{\chi}(q)]_{vv^{\prime}}=\delta_{v,v^{\prime}}\chi(q)$, and take the conventional Ornstein-Zernike form

[TABLE]

with velocity $v$ and mass $m$ (related to the correlation length, $\xi$, as $\xi=v/m$). Finally, the spin fluctuations couple to the electrons as described by the last contribution,

[TABLE]

with a coupling constant $g$ that, of course, can be absorbed into $\chi_{0}$ (or vice versa), but we will keep $g$ explicit here.

We follow Ref. Kozii et al., 2019 and focus on the leading-order correction of the coupling $g$ (second order, $\propto g^{2}$) to the free energy, $\mathcal{F}[\Delta_{\boldsymbol{k}},\boldsymbol{d}_{\boldsymbol{k}}]$, but emphasize that our expressions will differ from those of Ref. Kozii et al., 2019 as we consider a different type of fluctuations (centered around zero rather than finite momenta). These corrections are derived systematically by first integrating out the fermions, expanding the action to quadratic order in the bosonic fields $\boldsymbol{\phi}_{vq}$, which is sufficient to quadratic order in $g$, and integrating out the massive bosons.

We are interested in terms quartic in the superconducting order parameter, leaving us with the four distinct types of contributions in leading order in $g$ which are represented diagrammatically in Fig. 6. As indicated, the first diagram, in Fig. 6(a), only involves the zero-momentum and zero-frequency, $q=0$, fluctuations and it exactly captures the contributions of the simplified approach discussed in Sec. V of the main text.

The remaining diagrams—the self-energy in Fig. 6(b), the ladder in (c), and the vertex correction in (d) to the mean-field box-diagram—are fundamentally different: the loop-integrals involve integration over finite frequency and momentum of the bosonic fluctuations. As such, it is intuitively clear that they are less singular in the limit $m\rightarrow 0$ than the first diagram in Fig. 6(a), which is proportional to $\chi(\boldsymbol{q}=0,i\Omega_{n}=0)=m^{-2}$. In fact, these additional contributions can be shown to diverge with $\log(m)$ for small $m$. To illustrate this, let us consider the diagram in Fig. 6(b), which is proportional to

[TABLE]

where we introduced the function

[TABLE]

that only depends on the fermionic momenta and frequencies. Here, $\lambda_{\boldsymbol{k}}^{j}$ represent the basis functions of the involved superconducting order parameters. In the following, we will cut off the $\boldsymbol{q}$ integral by $\Lambda/v$ and expand $\xi_{\boldsymbol{k}+\boldsymbol{q}}\sim\xi_{\boldsymbol{k}}+\boldsymbol{v}_{\boldsymbol{k}}\cdot\boldsymbol{q}$, allowing us to write

[TABLE]

with the dimensionless velocity ratio $\hat{v}_{\boldsymbol{k}}:=|\boldsymbol{v}_{\boldsymbol{k}}|/v$. Since $|\omega_{n}|\geq\pi T$ and we work at finite temperature (given by the critical temperature of superconductivity), the term in the second line of Eq. (75) is finite in the limit $E\rightarrow 0$. The integral, thus, diverges as $\log(m)$ at small $E$ (infrared), as stated above. The other two diagrams in Fig. 6(c) and (d) can be analyzed in the same way and are also found to be subdominant as $m\rightarrow 0$ compared to the one in Fig. 6(a). This justifies the approach of Sec. V microscopically.

B.2 Enhanced symmetry in the one-band description

In the last part of this appendix, we discuss why $c_{+}\simeq c_{-}$ in Eq. (46) is expected. From the previous subsection of this appendix, we know that the results of the main text on fluctuation-induced superconductivity are captured by zero-momentum and zero-frequency fluctuations. Writing $\boldsymbol{m}_{v}:=\boldsymbol{\phi}_{q=0v}$ in Eq. (72) and generalizing to a momentum- and valley-dependent coupling constant, we here consider

[TABLE]

where $g^{m}_{v}(\boldsymbol{k})=g^{m}_{\bar{v}}(-\boldsymbol{k})$ as a consequence of time-reversal symmetry and we have, as before, assumed that we can focus on a single isolated electronic band. It is easy to see that $\mathcal{H}_{\text{MF}}+\mathcal{H}_{m}$, with $\mathcal{H}_{\text{MF}}$ in Eq. (57), is again invariant under the $C_{2}$ symmetry in Eq. (51) if we further replace

[TABLE]

While Eq. (44) is automatically invariant under Eq. (77), the coupling for the two-dimensional representation in Eq. (46) is invariant only if $c_{+}$ $=$ $c_{-}$. Consequently, multiband effects are required for nonzero $c_{+}-c_{-}$, wherefore we expect its value to be much smaller than $c_{+}+c_{-}$, as stated in the main text. We also checked by explicit calculation that $c_{+}\neq c_{-}$ is possible in a multiband description.

Appendix C Details for the complex representation

In this appendix, we present additional details of the different phases transforming under the complex representation $E$ of $C_{3}$.

As a starting point, it is helpful to chart out a phase diagram describing which of the triplet phases $E^{3_{s}}(a)$ to $E^{3_{s}}(h)$ is realized as a function of the quartic terms $b^{t}_{1,2,3,4,5}$ in Eq. (38). Upon recognizing that $b_{1}^{t}$ does not affect the form of the order parameter (but is assumed to be chosen so as to guarantee the stability of the expansion), we can conveniently display the phases as a function of $b_{j}^{t}/|b_{2}^{t}|$, $j=3,4,5$, discussing the two possible signs of $b_{2}^{t}$ separately. Such a phase diagram is drawn in Fig. 7.

As the main text contends, there are no independent terms involving $\sigma_{y}$ to add to the SU(2)${}_{+}\times$ SU(2)--invariant form of the free energy in Eq. (41). To see this, we note that it suffices to consider terms involving both $\Delta_{+}$ and $\Delta_{-}$ since terms with only $\Delta_{+}$ (or $\Delta_{-}$) have already been addressed in Sec. III.1. Among the terms that mix $\Delta_{+}$ and $\Delta_{-}$, the following are consistent with time-reversal and $C_{3}$ symmetry:

[TABLE]

However, all of these terms can be reformulated as

[TABLE]

so they do not constitute independent terms to add to Eq. (41).

In concluding this appendix, we present the explicit form of the free-energy (41) in terms of singlet and triplet components. Inserting $\Delta_{\mu}=\sigma_{0}\Delta^{s}_{\mu}+\boldsymbol{\sigma}\cdot\boldsymbol{d}_{\mu}$, $\mu=\pm$ in Eq. (41) and adding SU(2)+ $\times$ SU(2)- symmetry breaking only at the level of the quadratic terms, one arrives at

[TABLE]

where $\bar{\mu}=-$ for $\mu=+$ and vice versa. Due to the fewer number of independent parameters in Eq. (41), there are many relations between the different coefficients $\beta_{1},\dots,\beta_{12}$, namely:

[TABLE]

It is not difficult to observe that the five different purely triplet quartic terms, $\beta_{j=3,4,5,6,7},$ are all independent. Consequently, we can parametrize all twelve $\beta_{j}$ in terms of the five purely triplet terms and realize all of the triplet states of Sec. IV.1.