$\eta$-pairing in correlated fermion models with spin-orbit coupling

TL;DR

This paper extends $ta$-pairing theory to Hubbard models with spin-orbit coupling, revealing exact pseudospin symmetry and conditions for $ta$-pairing states in various topological and semimetal systems.

Contribution

It generalizes $ta$-pairing to spin-orbit coupled Hubbard models, identifying exact symmetry conditions and analyzing phase stability in topological semimetals.

Findings

Exact $ta$-pairing conditions in SOC Hubbard models.

Identification of stable $ta$-pairing and charge-density-wave phases.

Existence of an exactly solvable multicritical line.

Abstract

We generalize the $\eta$-pairing theory in Hubbard models to the ones with spin-orbit coupling (SOC) and obtain the conditions under which the $\eta$-pairing operator is an eigenoperator of the Hamiltonian. The $\eta$ pairing thus reveals an exact $SU(2)$ pseudospin symmetry in our spin-orbit coupled Hubbard model, even though the $SU(2)$ spin symmetry is explicitly broken by the SOC. In particular, these exact results can be applied to a variety of Hubbard models with SOC on either bipartite or non-bipartite lattices, whose noninteracting limit can be a Dirac semimetal, a Weyl semimetal, a nodal-line semimetal, and a Chern insulator. The $\eta$ pairing conditions also impose constraints on the band topology of these systems. We then construct and focus on an interacting Dirac-semimetal model, which exhibits an exact pseudospin symmetry with fine-tuned parameters. The stability regions for the exact $\eta$-pairing ground states (with momentum $\bm{\pi}$ or $\bm{0}$) and the exact charge-density-wave ground states are established. Between these distinct symmetry-breaking phases, there exists an exactly solvable multicritical line. In the end, we discuss possible experimental realizations of our results.

Tables

Table 1. Table 1: This table summarizes the 16 local eigenstates and their respective eigenvalues of the bond Hamiltonians ( 8 ), for B = − t 𝐵 𝑡 B=-t . The 16 local bases are defined as follows: | 00 ⟩ ket 00 |00\rangle denotes the vacuum, | 22 ⟩ = c i ↑ † ​ c i ↓ † ​ c j ↑ † ​ c j ↓ † ​ | 00 ⟩ ket 22 superscript subscript 𝑐 ↑ 𝑖 absent † superscript subscript 𝑐 ↓ 𝑖 absent † superscript subscript 𝑐 ↑ 𝑗 absent † superscript subscript 𝑐 ↓ 𝑗 absent † ket 00 |22\rangle=c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}|00\rangle , | 20 ⟩ = c i ↑ † ​ c i ↓ † ​ | 00 ⟩ ket 20 superscript subscript 𝑐 ↑ 𝑖 absent † superscript subscript 𝑐 ↓ 𝑖 absent † ket 00 |20\rangle=c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}|00\rangle , | 02 ⟩ = c j ↑ † ​ c j ↓ † ​ | 00 ⟩ ket 02 superscript subscript 𝑐 ↑ 𝑗 absent † superscript subscript 𝑐 ↓ 𝑗 absent † ket 00 |02\rangle=c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}|00\rangle , | σ ​ σ ′ ⟩ = c i ​ σ † ​ c j ​ σ ′ † ​ | 00 ⟩ ket 𝜎 superscript 𝜎 ′ superscript subscript 𝑐 𝑖 𝜎 † superscript subscript 𝑐 𝑗 superscript 𝜎 ′ † ket 00 |\sigma\sigma^{\prime}\rangle=c_{i\sigma}^{\dagger}c_{j\sigma^{\prime}}^{\dagger}|00\rangle , | σ ​ 0 ⟩ = c i ​ σ † ​ | 00 ⟩ ket 𝜎 0 superscript subscript 𝑐 𝑖 𝜎 † ket 00 |\sigma 0\rangle=c_{i\sigma}^{\dagger}|00\rangle , | 0 ​ σ ⟩ = c j ​ σ † ​ | 00 ⟩ ket 0 𝜎 superscript subscript 𝑐 𝑗 𝜎 † ket 00 |0\sigma\rangle=c_{j\sigma}^{\dagger}|00\rangle , | σ ​ 2 ⟩ = c i ​ σ † ​ c j ↑ † ​ c j ↓ † ​ | 00 ⟩ ket 𝜎 2 superscript subscript 𝑐 𝑖 𝜎 † superscript subscript 𝑐 ↑ 𝑗 absent † superscript subscript 𝑐 ↓ 𝑗 absent † ket 00 |\sigma 2\rangle=c_{i\sigma}^{\dagger}c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}|00\rangle , and | 2 ​ σ ⟩ = c i ↑ † ​ c i ↓ † ​ c j ​ σ † ​ | 00 ⟩ ket 2 𝜎 superscript subscript 𝑐 ↑ 𝑖 absent † superscript subscript 𝑐 ↓ 𝑖 absent † superscript subscript 𝑐 𝑗 𝜎 † ket 00 |2\sigma\rangle=c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}c_{j\sigma}^{\dagger}|00\rangle .

	$h_{ij}^x$	$h_{ij}^y$
Eigenvalue	Eigenstate	Eigenstate
0	$|00⟩$	$|00⟩$
$\frac{U}{2}-\mu +4V$	$|22⟩$	$|22⟩$
$\frac{U}{4}-\frac{\mu }{2}\pm P$	$|20⟩\pm |02⟩$	$|20⟩\pm |02⟩$
$V-\frac{\mu }{2}$	$|\sigma {\sigma }^{\prime }⟩$	$|\sigma {\sigma }^{\prime }⟩$
$\pm t-\frac{\mu }{4}$	$|\sigma 0⟩\pm |0\overline{\sigma }⟩$	$|\sigma 0⟩\pm \sigma i|0\overline{\sigma }⟩$
$\pm t+\frac{U}{4}-\frac{3\mu }{4}+2V$	$|\sigma 2⟩\mp |2\overline{\sigma }⟩$	$|\sigma 2⟩\mp \sigma i|2\overline{\sigma }⟩$

Equations

$$H$$

$$H_0$$

$$H_U$$

$${\eta }_{\mathbf{Q}}^{†}=\sum _{\mathbf{k}}{c}_{\mathbf{k}\uparrow }^{†}{c}_{\mathbf{Q}-\mathbf{k}\downarrow }^{†}=\sum _j{e}^{i\mathbf{Q}\cdot j}{c}_{j\uparrow }^{†}{c}_{j\downarrow }^{†}.$$

$$[H_0,{\eta }_{\mathbf{Q}}^{†}]=\lambda {\eta }_{\mathbf{Q}}^{†},$$

$${\alpha }_{\mathbf{k}}={\alpha }_{\mathbf{Q}-\mathbf{k}},{\beta }_{\mathbf{k}}={\beta }_{\mathbf{Q}-\mathbf{k}},$$

$${\epsilon }_{\mathbf{k}}+{\epsilon }_{\mathbf{Q}-\mathbf{k}}+{\gamma }_{\mathbf{k}}-{\gamma }_{\mathbf{Q}-\mathbf{k}}=\lambda ,$$

$$[H_U,{\eta }_{\mathbf{Q}}^{†}]=U{\eta }_{\mathbf{Q}}^{†},$$

$$[H,{\eta }_{\mathbf{Q}}^{†}]=(\lambda +U){\eta }_{\mathbf{Q}}^{†}.$$

$$H^{\prime }=\sum _{x\text{-bonds}}{h}_{ij}^x+\sum _{y\text{-bonds}}{h}_{ij}^y,$$

$$h_{ij}^x$$

$$h_{ij}^y$$

$$V$$

$$U$$

$$V$$

$$U$$

$$P$$

$$U$$

$$V$$

$$B=-t,V=P=0,U=2\mu \le -8|t|.$$

$$H_0$$

$$[H_0,{\eta }_{\mathbf{Q}}^{†}]$$

$$S$$

$${\alpha }_{\mathbf{k}}={\alpha }_{\mathbf{Q}-\mathbf{k}},{\beta }_{\mathbf{k}}={\beta }_{\mathbf{Q}-\mathbf{k}},{\gamma }_{\mathbf{k}}-{\gamma }_{\mathbf{Q}-\mathbf{k}}=\lambda =0,$$

$$\mathcal{H}(\mathbf{k})=\mathcal{H}(\mathbf{Q}-\mathbf{k}),$$

$$H_{\eta }=\sum _{\mathbf{k}}{\Delta }_{\mathbf{k}}(c_{\mathbf{k}\uparrow }^{†}{c}_{-\mathbf{k}\downarrow }^{†}+c_{-\mathbf{k}\downarrow }{c}_{\mathbf{k}\uparrow }),$$

$$\mathcal{C}\hat{Q}{\mathcal{C}}^{-1}=-\hat{Q}.$$

$${\mathcal{C}}^2=(-1{)}^{\hat{Q}},$$

$${\mathcal{T}}_f=\mathcal{C}\mathcal{K},$$

$$\mathcal{C}=\left(\begin{array}[]{cc}i\tau_{y}&0\\ 0&\mathbf{I}\\ \end{array}\right)^{\otimes N},$$

$${\mathcal{C}}^2=(-1{)}^{{\hat{N}}_f}.$$

$$\mathcal{C}{c}_{j\sigma }{\mathcal{C}}^{-1}=-\sigma e^{i\mathbf{Q}\cdot j}{c}_{j,-\sigma }^{†},$$

$$[J_{+},{\Delta }_z]=-2{\Delta }_{+},[J_{-},{\Delta }_{+}]=-{\Delta }_z,[J_z,{\Delta }_{+}]={\Delta }_{+},$$

$$D_{+}(t,t^{\prime })=-\frac{i}{N}\Theta (t-t^{\prime })⟨[J_{+}(t),{\Delta }_z(t^{\prime })]⟩.$$

Full text

I introduction

Exact results or solutions of interacting many-electron models beyond one spatial dimension are rare and highly appreciated. One important advance in studying the Hubbard model is the so-called $\eta$-pairing theory, which is an exact result due to Yang Yang (1989). The $\eta$-pairing operator allows for the construction of exact eigenstates possessing off-diagonal long-range order (ODLRO) and reveals an $SU(2)$ pseudospin symmetry Zhang (1990) of the Hubbard model. More importantly, the $\eta$-pairing states can serve as exact superconducting ground states in a class of strongly correlated electronic models Essler et al. (1992, 1993); Shen and Qiu (1993); de Boer et al. (1995), which is beyond the framework of BCS theory. Studies inspired by the $\eta$-pairing theory have since been pursued, including the unification of various symmetry-breaking orders (and related multicritical behavior study) Zhang (1990, 1997); Fradkin et al. (2015); Roy and Foster (2018); Boyack et al. (2020); Szabo and Roy (2020), the construction of an $SU(2)$ gauge theory of the Hubbard model Kim (2006, 2007); Hermele (2007), the calculation of entanglement entropy of $\eta$-pairing states Fan and Lloyd (2005); Vafek et al. (2017), and the study of quantum thermalization physics and quantum many-body scars in Hubbard-like models Nandkishore and Huse (2015); Garrison et al. (2017); Veness et al. (2017); Yu et al. (2018); Moudgalya et al. (2020); Mark and Motrunich (2020).

It is now recognized that spin-orbit coupling (SOC) can produce a wealth of topological phases and physical systems with SOC are currently under intense investigation. In particular, a number of theoretical understandings have been obtained for strongly correlated electronic systems with SOC Pesin and Balents (2010); Witczak-Krempa et al. (2014a); Rau et al. (2016). Nevertheless, most of these theoretical investigations used various approximate methods such as mean field theory and more rigorous and clearer physical understanding of these systems deserves further studies. Therefore, exact results or solutions for spin-orbit coupled interacting electronic systems would be valuable.

In this work, we first show that Yang’s $\eta$-pairing theory can be generalized to a class of Hubbard models with SOC and derive the equations which make the $\eta$-pairing operator an eigenoperator of the Hamiltonian. The exact $SU(2)$ pseudospin symmetry can be then preserved in these systems, while the $SU(2)$ spin symmetry is explicitly broken by the SOC. We then exemplify these exact results with four concrete spin-orbit coupled systems, including a Dirac semimetal, a Weyl semimetal, a nodal-line semimetal, and a Chern insulator, in the noninteracting limit. Some consequences of the pseudospin symmetry and the relation between the present study and recent works are discussed. To convert the $\eta$-pairing eigenstates into ground states, additional interactions besides the Hubbard $U$ are needed. In particular, we construct and focus on an interacting Dirac-semimetal model on the square lattice. By making the Hamiltonian frustration-free in some parameter regions, we establish the exact $\eta$-pairing ground states and the exact charge-density-wave (CDW) ground states, where the $\eta$-pairing ground states with momentum $\bm{\pi}$ and $\bm{0}$ are identified with a pair-density-wave (PDW) phase and a uniform superconductor, respectively. Moreover, there exists an exactly solvable multicritical line between the three symmetry-breaking phases. In the end, we discuss possible experimental realizations of our results in optical lattices.

This paper is organized as follows: In Sec. II we introduce our model and generalize the $\eta$-pairing theory. In Sec. III we exemplify the generalized $\eta$-pairing theory with four concrete models. An interacting Dirac-semimetal model with exact $\eta$-pairing ground states is presented in Sec. IV and we identify an exactly solvable multicritical line in Sec. V. We discuss possible experimental realizations in Sec. VI and present a summary and discussion in Sec. VII. Certain details can be found in the appendices.

II Model and generalization of $\eta$ pairing

The generic Hamiltonian for our spin-orbit coupled Hubbard models in arbitrary spatial dimensions is given by

[TABLE]

The noninteracting part $H_{0}$ is written in momentum space, where $C_{\bm{k}}^{\dagger}=(c_{\bm{k}\uparrow}^{\dagger},c_{\bm{k}\downarrow}^{\dagger})$. In $H_{0}$, the three momentum-dependent real coefficients $\alpha_{\bm{k}},\beta_{\bm{k}},\gamma_{\bm{k}}$ with Pauli matrices $\sigma_{x,y,z}$ represent the SOC, and the term $\varepsilon_{\bm{k}}$ with the identity matrix $\sigma_{0}$ is the usual energy band without SOC. The interacting part $H_{U}$ written in real space is the on-site Hubbard interaction, where $n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$.

In the absence of SOC ($\alpha_{\bm{k}}=\beta_{\bm{k}}=\gamma_{\bm{k}}=0$), Yang finds that the so-called $\eta$-pairing operator $\eta_{\bm{\pi}}^{\dagger}=\sum_{j}e^{i\bm{\pi}\cdot j}c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}$ obeys the commutation relations $[H_{0},\eta_{\bm{\pi}}^{\dagger}]=-2\mu\eta_{\bm{\pi}}^{\dagger}$ and $[H_{U},\eta_{\bm{\pi}}^{\dagger}]=U\eta_{\bm{\pi}}^{\dagger}$, where $\mu$ is the chemical potential and $\bm{\pi}=(\pi,\pi,\pi)$ or $(\pi,\pi)$ for the cubic or square lattice. The total Hamiltonian and $\eta_{\bm{\pi}}^{\dagger}$ thus obey $[H,\eta_{\bm{\pi}}^{\dagger}]=(U-2\mu)\eta_{\bm{\pi}}^{\dagger}$, and one can construct many exact eigenstates $H(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle=m(U-2\mu)(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle$ ($m=0,1,...$) possessing ODLRO.

In Yang’s observation, the key point that makes the possible $\eta$ pairing is the fact that the energy band $\varepsilon_{\bm{k}}$ (for the cubic or square lattice without SOC) has the property $\varepsilon_{\bm{k}}+\varepsilon_{\bm{\pi}-\bm{k}}=-2\mu$ (independent of $\bm{k}$), which leads to the commutation relation $[H_{0},\eta_{\bm{\pi}}^{\dagger}]=-2\mu\eta_{\bm{\pi}}^{\dagger}$.

Therefore, to generalize the $\eta$-pairing theory to Hubbard models with SOC, e.g., Eq.(1), we need some new constraints on the coefficients $\alpha_{\bm{k}},\beta_{\bm{k}}$, and $\gamma_{\bm{k}}$. To achieve these new constraints, we first define the $\eta$-pairing operator with a generic momentum $\bm{Q}$ as

[TABLE]

It can then be shown that (see Appendix A)

[TABLE]

if $\alpha_{\bm{k}},\beta_{\bm{k}},\gamma_{\bm{k}},\varepsilon_{\bm{k}}$, and $\bm{Q}$ satisfy the following equations

[TABLE]

where $\lambda$ is a constant independent of $\bm{k}$. The commutator of the Hubbard interaction and the $\eta$-pairing operator reads

[TABLE]

which holds for arbitrary $\bm{Q}$. It follows from Eqs.(3) and (5) that the $\eta$-pairing operator is now an eigenoperator (note the different meaning of an eigenoperator introduced in the context of the rigged Hilbert space Roberts (1966a, b)) of the total Hamiltonian, namely

[TABLE]

Therefore, one can construct many exact eigenstates $H(\eta_{\bm{Q}}^{\dagger})^{m}|0\rangle=m(\lambda+U)(\eta_{\bm{Q}}^{\dagger})^{m}|0\rangle$ possessing ODLRO.

The $SU(2)$ spin symmetry is explicitly broken by the SOC terms in Eq.(1). However, from Eq.(6), we see that our spin-orbit coupled Hubbard model possesses an exact $SU(2)$ pseudospin symmetry when $\lambda+U=0$. The term pseudospin is attributed to the $SU(2)$ algebra formed by $\eta$-pairing operators, i.e., $[J_{\alpha},J_{\beta}]=i\epsilon_{\alpha\beta\gamma}J_{\gamma}\ (\alpha,\beta,\gamma=x,y,z)$, where $J_{x}=\frac{1}{2}(\eta_{\bm{Q}}^{\dagger}+\eta_{\bm{Q}}),J_{y}=\frac{1}{2i}(\eta_{\bm{Q}}^{\dagger}-\eta_{\bm{Q}})$, $J_{z}=\frac{1}{2}(\hat{N}_{e}-N)$, $\hat{N}_{e}=\sum_{i\sigma}n_{i\sigma}$, and $N$ is the number of lattice sites.

Physically, the pseudospin symmetry can be understood via a particle-hole transformation for the spin-down electrons: $f_{j\uparrow}^{\dagger}=c_{j\uparrow}^{\dagger},f_{j\downarrow}^{\dagger}=e^{-i\bm{Q}\cdot j}c_{j\downarrow}$. The pseudospin operators can now be rewritten in terms of the $f$ fermions as $J_{\alpha}=\frac{1}{2}\sum_{j}f_{j}^{\dagger}\tau_{\alpha}f_{j}$, where $\tau_{\alpha}$ are the Pauli matrices and $f_{j}^{\dagger}=(f_{j\uparrow}^{\dagger},f_{j\downarrow}^{\dagger})$. Thus, the pseudospin can be interpreted as the spin of $f$ fermions. This is reminiscent of the $SU(2)$ gauge structure in the fermionic-spinon representation of spin operators Affleck and Marston (1988); Dagotto et al. (1988); Wen (2002). It is important to emphasize that the $\eta$-pairing states do not have the pseudospin symmetry. In fact, they correspond to the ferromagnetic states of $f$ fermions with maximal total ”spin” $N/2$, i.e., $\bm{J}^{2}(\eta_{\bm{Q}}^{\dagger})^{m}|0\rangle=\frac{N}{2}(\frac{N}{2}+1)(\eta_{\bm{Q}}^{\dagger})^{m}|0\rangle$ and $J_{z}(\eta_{\bm{Q}}^{\dagger})^{m}|0\rangle=(m-\frac{N}{2})(\eta_{\bm{Q}}^{\dagger})^{m}|0\rangle$ (where $\bm{J}^{2}=\sum_{\alpha}J_{\alpha}^{2}$), which break the $f$ fermion’s spin symmetry and hence break the pseudospin symmetry.

III Application of the generalized $\eta$-pairing theory

Based on the above exact results for model (1), we find that the the $\eta$ pairing can exist in a variety of spin-orbit coupled systems. In particular, we apply these exact results to four concrete examples, which include three gapless and one gapped noninteracting topological phases due to SOC.

The first example is a Dirac-semimetal Hubbard model on a square lattice. Its Hamiltonian is given by Eq.(1) with $\alpha_{\bm{k}}=2t^{\prime}\sin k_{x},\beta_{\bm{k}}=2t^{\prime}\sin k_{y},\gamma_{\bm{k}}=0$, and $\varepsilon_{\bm{k}}=-2t(\cos k_{x}+\cos k_{y})-\mu$. The lattice spacing has been set to unity throughout our paper. There are four Dirac points $\bm{K}=(n_{1}\pi,n_{2}\pi)$ in momentum space at $U=0$, where $n_{1,2}=0$ or $1$. Specifically, near the doubly degenerate point (e.g., $\bm{K}=(0,0)$) of the energy bands, the Bloch Hamiltonian $\mathcal{H}(\bm{k})=\alpha_{\bm{k}}\sigma_{x}+\beta_{\bm{k}}\sigma_{y}+\gamma_{\bm{k}}\sigma_{z}+\varepsilon_{\bm{k}}\sigma_{0}$ to the leading order reads $\mathcal{H}(\bm{k})\simeq 2t^{\prime}k_{x}\sigma_{x}+2t^{\prime}k_{y}\sigma_{y}+(-4t-\mu)$, which is a 2-dimensional massless Dirac Hamiltonian Jackiw (1984); Haldane (1988); Melikyan and Tešanović (2007). Substituting the expressions of $\alpha_{\bm{k}},\beta_{\bm{k}},\gamma_{\bm{k}}$ and $\varepsilon_{\bm{k}}$ into Eq.(4), we find that the $\eta$-pairing momentum $\bm{Q}=\bm{\pi}$ and the constant $\lambda=-2\mu$.

Similarly, the next example we consider is a Weyl-semimetal Hubbard model on a cubic lattice Witczak-Krempa et al. (2014b). Its Hamiltonian is given by Eq.(1) with $\alpha_{\bm{k}}=2t^{\prime}\sin k_{x},\beta_{\bm{k}}=2t^{\prime}\sin k_{y},\gamma_{\bm{k}}=2t^{\prime}\sin k_{z}$, and $\varepsilon_{\bm{k}}=-2t(\cos k_{x}+\cos k_{y}+\cos k_{z})-\mu$. There are totally eight Weyl points $\bm{K}=(n_{1}\pi,n_{2}\pi,n_{3}\pi)$ at $U=0$, where $n_{1,2,3}=0$ or $1$. As can be seen from Eq.(4), the $\eta$-pairing momentum $\bm{Q}=\bm{\pi}$ and the constant $\lambda=-2\mu$. Notice that for our Dirac- or Weyl-semimetal Hubbard model, the original $\eta$-pairing theory Yang (1989) is recovered in the absence of SOC, e.g., $t^{\prime}=0$.

Our third example is a nodal-line-semimetal Hubbard model on a cubic lattice Roy (2017). Its Hamiltonian is given by Eq.(1) with $\alpha_{\bm{k}}=t(\cos k_{x}+\cos k_{y}-b),\beta_{\bm{k}}=t^{\prime}\sin k_{z},\gamma_{\bm{k}}=0$, and $\varepsilon_{\bm{k}}=-\mu$. At $U=0$, upon expanding $\alpha_{\bm{k}}$ and $\beta_{\bm{k}}$ near the $\bm{k}=\bm{0}$ point, we obtain a nodal ring with radius $\sqrt{2(2-b)}$ pinned at the $k_{z}=0$ plane. We notice that, from Eq.(4), the $\eta$-pairing momentum $\bm{Q}=(0,0,\pi)$ and the constant $\lambda=-2\mu$.

The last example we consider is a Chern-insulator Hubbard model on a triangular lattice Li et al. (2016). Its Hamiltonian is given by Eq.(1) with $\alpha_{\bm{k}}=-2t\cos(\frac{k_{x}}{2}-\frac{\sqrt{3}k_{y}}{2}),\beta_{\bm{k}}=-2t\cos(\frac{k_{x}}{2}+\frac{\sqrt{3}k_{y}}{2}),\gamma_{\bm{k}}=-2t\cos k_{x}$, and $\varepsilon_{\bm{k}}=-\mu$. This model breaks time-reversal symmetry and $H_{0}$ has a nonzero band Chern number $\pm 2$. In the large-$U$ limit, this system is described effectively by the triangular lattice Kitaev-Heisenberg spin model Li et al. (2016, 2015). We see that, from Eq.(4), the $\eta$-pairing momentum $\bm{Q}=\bm{0}$ and the constant $\lambda=-2\mu$. Note that the $\eta$-pairing operator with $\bm{Q}=\bm{0}$ is simply the $s$-wave Cooper-pairing operator in the usual BCS theory. In fact, the first generalization of $\eta$ pairing to a triangular lattice was obtained by adding a staggered $\pi/2$ flux through each triangle plaquette Zhai (2005), in the absence of SOC.

Now, a few remarks are in order. Firstly, as shown in Appendix B, the $\eta$ pairing conditions Eq.(4) impose constraints on the topological properties of spin-orbit coupled systems. For example, the number of Dirac points in a 2-dimensional Dirac-semimetal must be a multiple of four, and the band Chern number of a Chern insulator must be an even integer. These have been embodied in the above concrete models.

Secondly, we notice that the $SU(2)$ pseudospin symmetry group contains a finite subgroup, namely the charge-conjugation symmetry (see Appendix D), which leads to a vanishing $U(1)$ charge $\langle Q_{i}\rangle\equiv\langle n_{i}-1\rangle=0$ on each site, where $\langle A\rangle$ denotes the ensemble average of the operator $A$ at any temperature. For the above four examples of spin-orbit coupled Hubbard model (1), the exact $SU(2)$ pseudospin symmetry is respected when $\lambda+U=0$ or simply $\mu=U/2$, so the electron density is fixed at half filling (independent of the temperature) when $\mu=U/2$. Another consequence of the exact pseudospin symmetry is that the ODLRO and CDW long-range order must either be both absent or coexist Zhang (1990); Pu and Shen (1994). In particular, the PDW phase (see below) possesses these two orders, whose parent Hamiltonian Eq.(7) exhibits an exact pseudospin symmetry (see Appendix F).

Lastly, beyond these four concrete topological Hubbard models, the exact $SU(2)$ pseudospin symmetry has also been found in other Hubbard-like models, such as the extended Falicov-Kimball model Fujiuchi et al. (2019) and the Hubbard model with a Zeeman term Fava et al. (2020), which are in fact special examples of our general spin-orbit coupled Hubbard model (1). Therefore, the $\eta$ pairing structure and the associated $SU(2)$ pseudospin symmetry of the models in Refs.Fujiuchi et al. (2019); Fava et al. (2020) can be readily obtained by using our generalized $\eta$-pairing theory for Hubbard models with SOC.

IV Construction of the exact $\eta$-pairing ground states

In general, the exact $\eta$-pairing eigenstates of the Hubbard model (1) are not ground states Yang (1989). In fact, beyond one spatial dimension, the Hubbard model has no exact solution to this day. However, some extended Hubbard models containing additional interactions are exactly solvable by constructing the ground state wave functions explicitly Essler et al. (1992, 1993); de Boer et al. (1995); de Boer and Schadschneider (1995). Following the method used in Ref.de Boer and Schadschneider (1995), we would further show that the $\eta$-pairing states can be constructed as the exact ground states of extended Hubbard models with SOC.

For the sake of concreteness, let us consider a square lattice Dirac-semimetal model similar to the one discussed above. Besides the on-site Hubbard $U$, we now introduce additional bond interactions, and the resulting new Hamiltonian can be expressed as

[TABLE]

where the summation $\sum_{x\texttt{-bonds}}$ ($\sum_{y\texttt{-bonds}}$) runs over all the bonds along the $x$ ($y$) direction of the square lattice. The bond Hamiltonians are given by

[TABLE]

where $c_{i}^{\dagger}=(c_{i\uparrow}^{\dagger},c_{i\downarrow}^{\dagger})$, $n_{i}=n_{i\uparrow}+n_{i\downarrow}$, and $\bar{\sigma}=-\sigma$. The noninteracting part $t$ represents SOC which is responsible for the Dirac points in momentum space. The $B$ term in Eq.(8) is known as the bond-charge interaction Kivelson et al. (1987), which completely breaks the $SU(2)$ spin-rotation symmetry and can be thought of as due to SOC Montorsi and Rasetti (1991). Note that it differs from the usual $SU(2)$ spin symmetric bond-charge interaction $\sum_{\sigma}(c_{i\sigma}^{\dagger}c_{j\sigma}+c_{j\sigma}^{\dagger}c_{i\sigma})(n_{i\bar{\sigma}}+n_{j\bar{\sigma}})$ by spin flips on the bonds. In addition, the terms $V$ and $P$ denote the nearest-neighbor Coulomb interaction and the pair-hopping term, respectively.

We note that the $\eta$-pairing operator with momentum $\bm{\pi}$ can be an eigenoperator of our extended Hubbard model with SOC, and hence Eq.(7) exhibits an exact pseudospin symmetry with fine-tuned parameters, as shown in Appendix F. In the absence of the bond interactions ($B=V=P=0$), Eq.(7) reduces to the spin-orbit coupled Hubbard model (1), and one can show that, using the generalized $\eta$-pairing theory discussed above, the $\eta$-pairing operator with zero momentum is an eigenoperator.

To construct the exact ground states of Eq.(7), the basic idea is to identify some parameter regions where the Hamiltonian (7) is frustration-free: The ground states of Eq.(7) are simultaneous ground states of each and every local bond Hamiltonian $h_{ij}^{x,y}$ in Eq.(8).

To find the parameter region (i.e., restrictions on $U$, $B$, $V$, and $P$ in the form of equalities and inequalities) where the $\eta$-pairing states are frustration-free ground states, let us diagonalize the local bond Hamiltonians (8). Here, each bond Hamiltonian $h_{ij}^{x,y}$ has 16 local eigenstates and their respective energies, which are summarized in Table 1. Notice that we have set $B=-t$ in the calculation Essler et al. (1992, 1993); de Boer et al. (1995); de Boer and Schadschneider (1995). For $B\neq-t$, the local eigenstates become complicated and it is difficult to see whether there exists a global frustration-free ground state de Boer and Schadschneider (1995).

We note that the $\eta$-pairing states with momentum $\bm{\pi}=(\pi,\pi)$ , i.e., $(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle$, can be built completely from the local states $|00\rangle$, $|22\rangle$, and $|20\rangle-|02\rangle$. Thus, $(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle$ will be the common ground states of all the bond Hamiltonians $h_{ij}^{x}$ and $h_{ij}^{y}$, if $|00\rangle$, $|22\rangle$, and $|20\rangle-|02\rangle$ are local ground states. This requires that, from Table 1, $0=\frac{U}{2}-\mu+4V=\frac{U}{4}-\frac{\mu}{2}-P$ are the minimum eigenvalues, which yields the following constraints

[TABLE]

Therefore, inside this parameter region together with $B=-t$, the exact ground states of our interacting Dirac-semimetal model (7) are $(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle$ ($m=0,1,...,N$). The zero temperature ODLRO for this phase can be readily evaluated: $\langle c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}c_{j\downarrow}c_{j\uparrow}\rangle_{T=0}=\sum_{m}\frac{\langle c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}c_{j\downarrow}c_{j\uparrow}\rangle_{m}}{N+1}=\frac{1}{6}e^{i\bm{\pi}\cdot(j-i)}$, where $\langle c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}c_{j\downarrow}c_{j\uparrow}\rangle_{m}=\frac{m(N-m)}{N(N-1)}e^{i\bm{\pi}\cdot(j-i)}$ is the ODLRO for each ground state $(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle$ Yang (1989). It is this non-vanishing ODLRO $\langle c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}c_{j\downarrow}c_{j\uparrow}\rangle_{T=0}$ at large distance $|i-j|\rightarrow\infty$ that characterizes the superconducting phase Yang (1962) inside the region of Eq.(9). This phase breaks the $U(1)$ particle-number symmetry, say, the number of electrons $2m$ in each ground state $(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle$ is not fixed with changing $m$. In addition, the lattice translation symmetry is also broken due to the nonzero pairing momentum $\bm{\pi}$: Each ground state transforms as $(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle\rightarrow(-1)^{m}(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle$ under the translation by one lattice spacing along the $x$ or $y$ direction, but the factor $(-1)^{m}$ depends on the parity of $m$. We thus conclude that the above superconducting phase is in fact a PDW phase.

In the same way, one can establish the parameter region in which the $\eta$-pairing states with momentum $\bm{0}$ become the exact ground states. $(\eta_{\bm{0}}^{\dagger})^{m}|0\rangle$ is now built completely from the local states $|00\rangle$, $|22\rangle$, and $|20\rangle+|02\rangle$, which have to be made local ground states. Again from Table 1, we get the following constraints

[TABLE]

Similar to the $\eta$-pairing states with momentum $\bm{\pi}$, the phase inside this region also possesses a non-vanishing ODLRO $\langle c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}c_{j\downarrow}c_{j\uparrow}\rangle_{T=0}=\frac{1}{6}$ and is thus a superconducting phase. The global $U(1)$ symmetry is spontaneously broken, due to the fluctuating electron numbers in the degenerate ground states $(\eta_{\bm{0}}^{\dagger})^{m}|0\rangle$. However, this phase preserves the lattice translation symmetry since the pairing momentum is $\bm{0}$. So it is a spatially uniform superconductor. In fact, as shown in Appendix C, it is closely related to the conventional $s$-wave superconductor.

It is interesting to see that the CDW state $|\texttt{CDW}\rangle=\prod_{i\in\mathcal{A}}c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}|0\rangle$ can also be an exact ground state, where $\mathcal{A}$ denotes one sublattice of the square lattice. Notice that $|\texttt{CDW}\rangle$ can be constructed completely from the local states $|20\rangle$ and $|02\rangle$, we get the following constraints

[TABLE]

The CDW phase obviously breaks the lattice translation symmetry, accompanied by a charge density modulation with the same periodicity as that of the PDW phase. We notice that some important aspects of the physics of this simple CDW order in topological semimetals have been explored recently Sehayek et al. (2020).

Eqs.(9)-(11) are exact results that establish the stability regions of superconducting (i.e., the $\eta$-pairing states) and CDW ground states. We thus see that, for sufficiently small $U$, superconducting ground states are stabilized by attractive $V$ and finite $P$, whereas a repulsive $V$ (without $P$) favors a CDW ground state.

V Identification of an exact multicritical line

Given the three exactly solvable regions (9)-(11), it is desirable to determine the phase transitions between these three symmetry-breaking phases. We note that the three regions (9)-(11) have no overlaps with each other, however, overlaps exist in their boundaries. Specifically, as shown in Appendix G, any two of the three boundaries of regions (9)-(11) share a common phase boundary, say

[TABLE]

As a consequence, Eq.(12) is in fact an exact multicritical line, where the PDW, uniform superconductor, and CDW phases meet.

A complementary approach to confirm this is by checking the analyticity of ground-state energy near the line (12). As can be seen from Table 1, the ground-state energy (per bond) vanishes continuously from each of the three phases toward the line (12). However, its first derivatives (with respect to, e.g., $U$ and $\mu$) have a jump at line (12), namely a kink, implying a quantum phase transition at this line.

More interestingly, the multicritical line (12) is also an exactly solvable line: Within this regime, we see that, from Table 1, the local ground states are $|00\rangle,|20\rangle,|02\rangle$, and $|22\rangle$. In this case, the ground-state space of Eq.(7) is spanned by the bases $\prod_{i}(c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger})^{l_{i}}|0\rangle$ ($l_{i}=0,1$) and is highly degenerate, which contains both the $\eta$-pairing states (with arbitrary momentum) and the CDW states.

VI Experimental considerations

A few remarks are now in order concerning possible experimental realizations of our results in optical lattices. The essential ingredients of the present study are SOC and interaction. In general, the SOC can be implemented artificially in an optical lattice by introducing synthetic gauge fields Dalibard et al. (2011) and the interaction can be widely tuned through Feshbach resonances Chin et al. (2010); Giorgini et al. (2008). In particular, concrete schemes to realize spin-orbit coupled Hubbard models with controlled band structures and interactions have been proposed Duan et al. (2003); Liu et al. (2004). In fact, the Fermi-Hubbard model has been realized in an optical lattice Schneider et al. (2008); Jordens et al. (2008); Strohmaier et al. (2007), where the kinetic energy, the interaction strength, and the chemical potential can be varied independently Hart et al. (2015); Mazurenko et al. (2017); Parsons et al. (2016); Boll et al. (2016); Cheuk et al. (2016); Brown et al. (2017). Given these, it is plausible that the $SU(2)$ pseudospin symmetry can be realized via fine tuning of the model parameters, e.g., $\mu$ or $U$ so that $\mu=U/2$ as discussed previously, in cold-atom systems Bloch et al. (2008); Dalibard et al. (2011); Chin et al. (2010); Giorgini et al. (2008).

On the other hand, the exactly solvable interacting Dirac-semimetal model (7) contains off-site interactions, such as the $B$, $V$, and $P$ terms in Eq.(8), which can arise from the matrix elements of the low-energy effective potential between fermionic atoms in an optical lattice Werner et al. (2005); Rosch et al. (2008). Recent experiments in cold fermionic Rydberg atoms in an optical lattice have achieved a high degree of control over these interactions, such as adjusting the sign, strength, and range of the interactions Bloch et al. (2008); Li and Sarma (2015); Carroll et al. (2004); Viteau et al. (2011); Schauss et al. (2012); Labuhn et al. (2016); Chin et al. (2010); Giorgini et al. (2008). Therefore, the ability to control and tune interactions of extended Hubbard models in optical lattices has paved the way for the realization of Hamiltonian (7) over a wide range of parameters. Moreover, the $\eta$-pairing states with zero momentum are closely related to the $s$-wave superconductivity (see Appendix C), whose experimental evidence in an optical lattice has already been established and well studied Chin et al. (2006); Mitra et al. (2017). We are therefore optimistic that the exact $\eta$-pairing ground states of Eq.(7) could be realized in optical lattices in the near future.

VII Summary and discussion

In summary, we have generalized Yang’s $\eta$-pairing theory to a class of Hubbard models with SOC, where a hidden $SU(2)$ symmetry is uncovered. We apply these exact results to four concrete spin-orbit coupled Hubbard models, which include three gapless and one gapped topological phases in the noninteracting limit. We then construct and focus on an interacting Dirac-semimetal model and establish the stability regions for the exact $\eta$-pairing ground states (with momentum $\bm{\pi}$ or $\bm{0}$) and the exact CDW ground states. Between these distinct symmetry-breaking phases, there exists an exactly solvable multicritical line.

The present work may provide new guidance in exploring exotic quantum phases, for instance a PDW phase here, in a wider range of physical systems, such as spin-orbit coupled systems with $SU(2)$ pseudospin symmetry. The exactly solvable interacting Dirac-semimetal model we constructed deserves further studies: This model hosts various exact symmetry-breaking ground states and an exact multicritical line between them, which should serve as a clean platform to study intertwined orders Fradkin et al. (2015) and quantum phase transitions between the Dirac-semimetal and symmetry-breaking phases Classen et al. (2015); Roy et al. (2016); Sato et al. (2017); Janssen et al. (2018); Boyack et al. (2020); Torres et al. (2020).

Acknowledgements.

I am especially grateful to Jian-Xin Li and Zhenyu Zhang for helpful discussions and encouragement during the course of this work. I also thank Meng Cheng, Zhao-Long Gu, Zhao-Yang Dong, Feng Tang, Rui Wang, and Xiaoyu Zhu for enlightening conversations. This work was supported by the NSF of China under Grant No. 11704341 and the Open Research Program of the Laboratory of Solid State Microstructures (Nanjing University) under Grant No. M30020.

Appendix A Derivation of Eq.(4) in the main text

Here, we show how to obtain Eq.(4) in the main text. We first rewrite the noninteracting part $H_{0}$ in Eq.(1) as

[TABLE]

It is now convenient to use the momentum space expression of $\eta$-pairing operator $\eta_{\bm{Q}}^{\dagger}=\sum_{\bm{k}}c_{{\bm{k}}\uparrow}^{\dagger}c_{\bm{Q}-{\bm{k}}\downarrow}^{\dagger}$ to calculate its commutator with $H_{0}$. After a direct and lengthy calculation, we get the following result

[TABLE]

Now, let us consider the first summation term $S=\sum_{\bm{k}}\alpha_{\bm{k}}c_{\bm{k}\downarrow}^{\dagger}c_{\bm{Q}-{\bm{k}}\downarrow}^{\dagger}$ in Eq.(14), and we have the following identities

[TABLE]

where the first line follows from the anti-commutation relation of fermion operators, and the substitution $\bm{k}=\bm{Q}-\bm{k^{\prime}}$ is used in the second line. As can be seen from Eq.(15), if $\alpha_{\bm{k}}=\alpha_{\bm{Q}-\bm{k}}$, then we have $S=-S$ and hence $S=0$. The second summation term in Eq.(14) also vanishes in the same way when $\alpha_{\bm{k}}=\alpha_{\bm{Q}-\bm{k}}$. Similarly, the third and the fourth summation terms in Eq.(14) vanish separately when $\beta_{\bm{k}}=\beta_{\bm{Q}-\bm{k}}$. Finally, the last summation term in Eq.(14) is equal to $\lambda\sum_{\bm{k}}c_{{\bm{k}}\uparrow}^{\dagger}c_{\bm{Q}-{\bm{k}}\downarrow}^{\dagger}=\lambda\eta_{\bm{Q}}^{\dagger}$, if $\varepsilon_{\bm{k}}+\varepsilon_{\bm{Q}-\bm{k}}+\gamma_{\bm{k}}-\gamma_{\bm{Q}-\bm{k}}=\lambda$, where $\lambda$ is a constant independent of $\bm{k}$.

In conclusion, when $\alpha_{\bm{k}}=\alpha_{\bm{Q}-\bm{k}}$, $\beta_{\bm{k}}=\beta_{\bm{Q}-\bm{k}}$, and $\varepsilon_{\bm{k}}+\varepsilon_{\bm{Q}-\bm{k}}+\gamma_{\bm{k}}-\gamma_{\bm{Q}-\bm{k}}=\lambda$, i.e., Eq.(4), we have $[H_{0},\eta_{\bm{Q}}^{\dagger}]=\lambda\eta_{\bm{Q}}^{\dagger}$.

Appendix B Implications of $\eta$-pairing theory on the topological properties of spin-orbit coupled systems

Based on our generalized $\eta$-pairing theory, in the following we will show that Eq.(4) can have nontrivial restrictions on the properties of band topology of spin-orbit coupled systems.

Since we are interested in the topological properties, the usual energy dispersion term $\varepsilon_{\bm{k}}\sigma_{0}$ in $H_{0}$ can be ignored because it does not affect the band topology. Eq.(4) is thus reduced to

[TABLE]

where the last equation vanishes because $\lambda$ is independent of $\bm{k}$ and we can take $\bm{k}=\frac{\bm{Q}}{2}$ such that $\gamma_{\frac{\bm{Q}}{2}}-\gamma_{\bm{Q}-\frac{\bm{Q}}{2}}=\gamma_{\frac{\bm{Q}}{2}}-\gamma_{\frac{\bm{Q}}{2}}=0$. Eq.(16) implies that the corresponding Bloch Hamiltonian $\mathcal{H}(\bm{k})=\alpha_{\bm{k}}\sigma_{x}+\beta_{\bm{k}}\sigma_{y}+\gamma_{\bm{k}}\sigma_{z}$ has the property

[TABLE]

and the Bloch wave-functions also satisfy $\psi(\bm{k})=\psi(\bm{Q}-\bm{k})$ (up to an unphysical phase factor). From the point of view of topology, this implies that each point $\bm{k}$ is topologically identical to the point $\bm{Q}-\bm{k}$ in momentum space.

Let us now focus on 2-dimensional systems. If we translate the Brillouin zone (BZ) by vector $\frac{\bm{Q}}{2}$ and choose its center as the new origin of coordinate in momentum space, then the above arguments mean that each point $\bm{k}$ is topologically identical to the point $-\bm{k}$. We now bipartite the BZ into two halves, then each half BZ is topologically identical to a sphere which is a closed manifold (See Appendix A 3 in Ref.Li et al. (2016) for the proof).

Therefore, we arrive at the following conclusions:

(1). If the spin-orbit coupled system is a 2-dimensional Dirac-semimetal, then the number of Dirac points must be a multiple of four.

Proof: Firstly, the two half BZs host equal number of Dirac points due to the symmetry property of Bloch Hamiltonian, e.g., Eq.(17). Secondly, each half BZ is a closed manifold where the Dirac points must appear in pairs.

(2). If the system is a Chern insulator, then the band Chern number must be an even integer. [See Appendix A 3 in Ref.Li et al. (2016) for the proof.]

Appendix C The $\eta$-pairing states are topologically trivial

In this section, we argue that the $\eta$-pairing states, e.g., $(\eta_{\bm{0}}^{\dagger})^{m}|0\rangle$ ($m=0,1,...,N$), are closely related to the conventional BCS superconductivity and hence they are topologically trivial. Specifically, let us consider the $s$-wave BCS mean-field Hamiltonian $H_{BCS}=\sum_{\bm{k}\sigma}\xi_{\bm{k}}c_{\bm{k}\sigma}^{\dagger}c_{\bm{k}\sigma}+\sum_{\bm{k}}\Delta_{\bm{k}}(c_{{\bm{k}}\uparrow}^{\dagger}c_{-{\bm{k}}\downarrow}^{\dagger}+c_{{-\bm{k}}\downarrow}c_{{\bm{k}}\uparrow})$, where $\Delta_{\bm{k}}<0$. The intra-pair distance between the two electrons in the usual Cooper pair is finite, whereas it is zero (i.e., on-site singlet pairing) in the $\eta$-pairing states. Therefore, the ground state of $H_{BCS}$ is in general not a linear combination of $(\eta_{\bm{0}}^{\dagger})^{m}|0\rangle$.

However, $H_{BCS}$ is topologically equivalent to a simpler Hamiltonian

[TABLE]

since the energy spectrum remains fully gapped during the deformation $\xi_{\bm{k}}\rightarrow 0$ and there should be no topological phase transition. The ground state wave function (up to a normalization factor) of $H_{\eta}$ is now readily written as $|G\rangle=\prod_{\bm{k}}(1+c_{{\bm{k}}\uparrow}^{\dagger}c_{-{\bm{k}}\downarrow}^{\dagger})|0\rangle=\sum_{m=0}^{N}\frac{1}{m!}(\eta_{\bm{0}}^{\dagger})^{m}|0\rangle$, which is a superposition of $\eta$-pairing states $(\eta_{\bm{0}}^{\dagger})^{m}|0\rangle$ Shirakawa et al. (2020). Since $H_{BCS}$ is topologically trivial, $H_{\eta}$ and hence $(\eta_{\bm{0}}^{\dagger})^{m}|0\rangle$ are also topologically trivial.

An alternative way to see that the $\eta$-pairing states are topologically trivial may be as follows. As shown in the main text, the $\eta$-pairing states can be mapped to the ferromagnetic states through a local particle-hole transformation. A local transformation should not change the global topological property of a state. Therefore, the $\eta$-pairing states are topologically trivial since the ferromagnetic states are topologically trivial.

Appendix D Charge-conjugation symmetry as a subgroup of the $SU(2)$ pseudospin symmetry

As shown in the main text, the generators of the $SU(2)$ pseudospin symmetry group are the pseudospin operators $J_{x},J_{y}$ and $J_{z}$. The global pseudospin-rotation by $\pi$ about the $J_{y}$ axis, say $\mathcal{C}=e^{i\pi J_{y}}$, flips the sign of $J_{z}$: $\mathcal{C}J_{z}\mathcal{C}^{-1}=-J_{z}$. We recall that the pseudospin operator $J_{z}$ is proportional to the $U(1)$ charge operator $\hat{Q}\equiv\hat{N}_{e}-N$, say $J_{z}=\frac{1}{2}\hat{Q}$. Therefore, the pseudospin-rotation $\mathcal{C}=e^{i\pi J_{y}}$ is in fact a charge conjugation Hermele (2007); You et al. (2017), i.e.,

[TABLE]

In fact, $\mathcal{C}=e^{i\pi J_{y}}$ flips the sign of the local operator $J_{i}^{z}\equiv\frac{1}{2}Q_{i}\equiv\frac{1}{2}(n_{i}-1)$ on each site. So, if the model Hamiltonian has the $SU(2)$ pseudospin symmetry (and hence the charge-conjugation symmetry), then $\langle Q_{i}\rangle=\langle J_{i}^{z}\rangle=0$ on each site, where $\langle A\rangle$ denotes the ensemble average of the operator $A$ at any temperature.

It is interesting to see that the charge-conjugation operator $\mathcal{C}$ squares to the charge parity You et al. (2017)

[TABLE]

since $\mathcal{C}^{2}=e^{i2\pi J_{y}}=e^{i2\pi J_{z}}=(-1)^{2J_{z}}$ and $J_{z}=\frac{1}{2}\hat{Q}$. An alternative way to understand Eq.(20) is to make an analogy between the charge-conjugation $\mathcal{C}$ and the ”time reversal” $\mathcal{T}_{f}$. Specifically, let us focus on the $f$ fermions introduced by a particle-hole transformation for the spin-down electrons and we have $J_{\alpha}=\frac{1}{2}\sum_{j}f_{j}^{\dagger}\tau_{\alpha}f_{j}$, as shown in the main text. The charge-conjugation $\mathcal{C}=e^{i\pi J_{y}}$ rotates the spin of $f$ fermions $\mathcal{C}f_{j\uparrow}\mathcal{C}^{-1}=-f_{j\downarrow}$, $\mathcal{C}f_{j\downarrow}\mathcal{C}^{-1}=f_{j\uparrow}$. The ”time reversal” $\mathcal{T}_{f}$ acts on these $f$-fermion operators the same as $\mathcal{C}$, except that $\mathcal{T}_{f}$ is an antiunitary operator, e.g., $\mathcal{T}_{f}i\mathcal{T}_{f}^{-1}=-i$. So, we can define $\mathcal{T}_{f}$ to be

[TABLE]

where $\mathcal{K}$ denotes the complex conjugation operator. To be precise, consider the $f$-fermion occupation number basis $\prod_{j=1}^{N}(f_{j\uparrow}^{\dagger})^{l_{j\uparrow}}(f_{j\downarrow}^{\dagger})^{l_{j\downarrow}}|0\rangle$, where $|0\rangle$ is the vacuum state of $f$ fermions and $l_{j\sigma}\in\{0,1\}$. Now the operator $\mathcal{K}$ is defined in this basis, e.g., the base kets $\prod_{j=1}^{N}(f_{j\uparrow}^{\dagger})^{l_{j\uparrow}}(f_{j\downarrow}^{\dagger})^{l_{j\downarrow}}|0\rangle$ are invariant under the action of $\mathcal{K}$. In this representation, the charge-conjugation $\mathcal{C}=e^{i\pi J_{y}}$ acts as

[TABLE]

where $\tau_{y}$ is the second Pauli matrix and $\mathbf{I}$ is the $2\times 2$ identity matrix. The unitary operator $\mathcal{C}$ thus turns out to be a real orthogonal matrix in the above representation, say $\mathcal{C}^{\ast}=\mathcal{C}$. Using Eq.(21), we have $\mathcal{T}_{f}^{2}=\mathcal{C}\mathcal{C}^{\ast}=\mathcal{C}^{2}$. Now we recall that the ”time reversal” operator $\mathcal{T}_{f}$ squares to the fermion number parity $\mathcal{T}_{f}^{2}=(-1)^{\hat{N}_{f}}$, with $\hat{N}_{f}=\sum_{j,\sigma}f_{j\sigma}^{\dagger}f_{j\sigma}$. We finally conclude that

[TABLE]

Note that Eq.(23) is indeed equivalent to Eq.(20), since $\hat{Q}$ and $\hat{N}_{f}$ have the same parity: $\hat{N}_{f}=2J_{z}+2\sum_{j}f_{j\downarrow}^{\dagger}f_{j\downarrow}=\hat{Q}+2\sum_{j}f_{j\downarrow}^{\dagger}f_{j\downarrow}$, where the term $2\sum_{j}f_{j\downarrow}^{\dagger}f_{j\downarrow}$ always has an even parity.

For the sake of completeness, we outline the action of the charge-conjugation operator $\mathcal{C}$ on the original electron operators, say

[TABLE]

for which reason $\mathcal{C}$ is also called particle-hole transformation.

Appendix E Implications of $\eta$-pairing theory on the superconducting properties in Hubbard models with SOC

Following Zhang’s idea in Ref.Zhang (1990), we will show that the $SU(2)$ pseudospin symmetry of our spin-orbit coupled Hubbard models would lead to a pair of new collective modes together with the usual Goldstone mode, if the system is in the superconducting phase.

The crucial point in Zhang’s observation is the fact that the transverse components of the pseudospin operators look like a superconducting order, while the $z$ component looks like a CDW order, up to a local phase factor $e^{i\bm{Q}\cdot j}$. So the multiplet of the particle-particle and particle-hole operators should transform as a vector under the pseudospin rotation

[TABLE]

meaning that the pseudospin operators can rotate the superconducting order and the CDW order into each other, where the pseudospin operators $J_{+}=J_{x}+iJ_{y}=\eta_{\bm{Q}}^{\dagger}$, $J_{-}=J_{+}^{\dagger}$, $J_{z}=\frac{1}{2}(\hat{N}_{e}-N)$, the s-wave Cooper-pairing operator $\Delta_{+}=\sum_{j}c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}$, $\Delta_{-}=\Delta_{+}^{\dagger}$, and the CDW operator $\Delta_{z}=\sum_{j\sigma}e^{-i\bm{Q}\cdot j}n_{j\sigma}$ with the wave-vector $\bm{Q}$.

To study the collective modes, let us consider the following response function

[TABLE]

Since the Hamiltonian is time-independent, we have $D_{+}(t,t^{\prime})=D_{+}(t-t^{\prime})=-\frac{i}{N}\Theta(t-t^{\prime})\langle[J_{+}(t-t^{\prime}),\Delta_{z}]\rangle$. From Eq.(6) in the main text, the time dependence of pseudospin operators in the Heisenberg representation can be determined explicitly, in particular, $J_{+}(t-t^{\prime})=e^{i(\lambda+U)(t-t^{\prime})}J_{+}$. Using the commutator in Eq.(25), Eq.(26) can be thus calculated exactly, i.e., $D_{+}(t-t^{\prime})=2i\rho\Theta(t-t^{\prime})e^{i(\lambda+U)(t-t^{\prime})}$, where $\rho=\frac{\langle\Delta_{+}\rangle}{N}$ denotes the superconducting order parameter.

Information about collective modes are now embedded in the analytic properties of the Fourier transform of Eq.(26), say

[TABLE]

where $\delta$ is a positive infinitesimal. Similarly, we have

[TABLE]

and

[TABLE]

Therefore, as can be seen from Eqs.(27), (28), and (29), if the ground state of our spin-orbit coupled Hubbard model is in the superconducting phase (i.e., $\rho\neq 0$), then there must exist a triplet of collective modes with energies $-(\lambda+U)$, $\lambda+U$, and 0, respectively. Note that the collective mode with energy 0 in Eq.(29) is simply the gapless Goldstone mode resulting from the spontaneous $U(1)$ symmetry breaking of superconductivity. In addition, we find that there are a pair of new collective modes with energies $\pm(\lambda+U)$, which is also a consequence of superconductivity.

Appendix F The $SU(2)$ pseudospin symmetry of Eq.(7) in the main text

As shown in the main text, the $\eta$-pairing operator can be an eigenoperator of a class of Hubbard models with SOC. In fact, the $\eta$-pairing operator can also be an eigenoperator of an extended Hubbard model with SOC, such as the interacting Dirac-semimetal model Eq.(7) we studied in the main text.

To see this, we calculate the commutator of the Hamiltonian (7) and the $\eta$-pairing operator $\eta_{\bm{\pi}}^{\dagger}$, which reads

[TABLE]

Therefore, the $\eta$-pairing operator $\eta_{\bm{\pi}}^{\dagger}$ becomes an eigenoperator of $H^{\prime}$ and $(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle$ become exact eigenstates with eigenvalues $m(U-2\mu-4P)$, when $B=-t,V=-\frac{P}{2}$. Moreover, if $B=-t,V=-\frac{P}{2}$, and $\mu=\frac{U}{2}+4V$, $H^{\prime}$ commutes with $\eta_{\bm{\pi}}^{\dagger}$ and hence possesses the $SU(2)$ pseudospin symmetry. Accordingly, within the regime of the PDW phase where $(\eta_{\bm{\pi}}^{\dagger})^{m}|0\rangle$ are the exact ground states (see Eq.(9) in the main text), the corresponding Hamiltonian (7) exhibits an exact pseudospin symmetry which gets spontaneously broken.

Appendix G Derivation of Eq.(12) in the main text

Eq.(12) describes the common phase boundary shared by the three exactly solvable regions (9), (10), and (11), which respectively host the PDW, uniform superconductor, and CDW phases. To obtain Eq.(12), we need first to work out the boundary of each region.

From Eq.(9), we see that its boundary $\partial_{\texttt{PDW}}$ consists of two disjoint sets

[TABLE]

Similarly, the boundary of Eq.(10) reads

[TABLE]

Lastly, the boundary of Eq.(11) reads

[TABLE]

Note that we have omitted $B=-t$ in all of the above equations.

From these equations, we find that any two of the boundaries $\partial_{\texttt{PDW}}$, $\partial_{\texttt{SC}}$, and $\partial_{\texttt{CDW}}$ have equal intersections

[TABLE]

As a consequence, Eqs.(9),(10),(11) in the main text share a common phase boundary

[TABLE]

which is a multicritical line between these distinct symmetry-breaking phases.