Phase diagram of the square-lattice Hubbard model with Rashba-type antisymmetric spin-orbit coupling

We clarify the ground-state phase diagram of the half-filled square-lattice Hubbard model with Rashba spin-orbit coupling (SOC) characterized by the spin-split energy bands due to broken inversion symmetry. Although the Rashba metals and insulating magnets have been studied well, the intermediate interaction strength of the system remained elusive due to the lack of appropriate theoretical tools to unbiasedly describe the large-scale magnetic structures. We complementarily apply four different methods (sine-square deformed mean-field theory, random phase approximation, Luttinger-Tisza method, and density matrix embedding theory) and succeed in capturing the incommensurate spin-density-wave (SDW) phases with very long spatial periods which were previously overlooked. The transition to the SDW phases from the metallic phase is driven by an unprecedented instability that nests the two parts of the Fermi surface carrying opposite spins. For large SOC, the spiral, stipe, and vortex phases are obtained, when the four Dirac points exist near the Fermi level and their whole linear dispersions nest by a wavelength $\pi$, opening a band gap. These two types of transition provide Fermiology that distinguishes the antisymmetric SOC systems, generating a variety of magnetic phases that start from the relatively weak correlation regime and continue to the strongly interacting limit.

Figure 1

(a) 2D square lattice with Rashba SOC in Eq. (3), where ${\mathit{n}}_{\mu }$ ($\mu =x,y$) denotes the unit vector that rotates the spins in hopping to its neighbors by the SU(2) gauge field shown schematically in the right panel. (b) Representative noninteracting energy band structure ${\varepsilon }_{\pm }\left(\mathit{k}\right)$ of ${\widehat{\mathcal{H}}}_0$ for $\theta =2\pi /3$. The left panel is the one along the symmetric lines. (c) Schematic illustration of the four Dirac cones with the same velocities. The broken arrows represent the $\mathit{q}=(0,\pi ),(\pi ,0)$ vector that “nests” the whole Dirac cones in forming the stripe and vortex phases.

Figure 2

(a) Ground-state magnetic phase diagram obtained by the SSDMF with $L=28$. There are three incommensurate spin-density-wave (SDW) orders labeled by SDWI, SDWII, and SDWIII. For large $\theta$, the system hosts the spiral, stripe, and vortex orders. (b) Spin structure factor (upper panels) and the spin-density distribution in real space (lower panels) given as density plots, normalizing the minimum to maximum values as $[0:1]$ for the SDW phases with single and double $\mathit{Q}$ (denoted as $Q$ and $2Q$). (c) Spin structure factor and the spin density as a function of $U$ for spiral, stripe, and vortex phases.

Figure 3

(a) Fermi surfaces and (b) the density plots of the largest eigenvalues of the bare magnetic susceptibility ${\lambda }_3\left(\mathit{q}\right)$ as functions of $\mathit{q}$ for several choices of $\theta$. The purple arrows in the upper panels represent the nesting vector $\mathit{Q}$, and the purple circles in the lower panels represent the peak position of ${\lambda }_3\left(\mathit{q}\right)$ within $q_x\ge q_y\ge 0$. For $\theta \ne 0$, there are two zero-energy Dirac points at $(\pi ,0)$ and $(0,\pi )$, and the other two Dirac points reach the zero energy at (0,0) and $(\pi ,\pi )$ for $\theta =\pi$, represented by the green circles. The density plot of ${\lambda }_3\left(\mathit{q}\right)$ is normalized to ${min}_{\mathit{q}}{\lambda }_3\left(\mathit{q}\right)=0$ and ${max}_{\mathit{q}}{\lambda }_3\left(\mathit{q}\right)=1$ for clarification. (c) Peak position of ${\lambda }_3\left(\mathit{q}\right)$ in the region $q_x\ge q_y\ge 0$ as a function of $\theta$. There are three phases separated by the boundaries at ${\theta }_{c1}^{\left(\mathrm{RPA}\right)}\sim \pi /4$ and ${\theta }_{c2}^{\left(\mathrm{RPA}\right)}\sim 2\pi /3$. The purple dashed lines are the wave numbers from Eqs. (19, 20, 21), which are in good agreement with those obtained by the RPA.

Figure 4

Comparison of the two vectors (a) ${\mathit{Q}}_{\mathrm{SDWI}}$ and (b) ${\mathit{Q}}_{\pi }$, regarding the contribution to ${\chi }_0\left(\mathit{Q}\right)$. The left panels show the spin orientation (black arrows) on the Fermi surfaces and how the nesting vectors connect the Fermi surfaces. The right panels show the density plot of the contributions to ${\lambda }_3\left(\mathit{Q}\right)$ at $\mathit{k}$ which interpolates to $\mathit{k}+\mathit{Q}$ in Eq. (12). Broken line in the right panels is the location of the Fermi surface for the guide to the eye.

Figure 5

Transverse and longitudinal susceptibilities (a) $U=0.1$ (metal) and (b) $U=3.0$ (SDWI) at $\theta =\pi /6$ and $k_{\mathrm{B}}T=0.05$. For the paramagnetic metallic phase, the peak position locates at $\mathit{q}=(\pi ,\pi )$ for the transverse susceptibilities, and $\mathit{q}=(\pm \pi ,\pm q^{*}),(\pm q^{*},\pm \pi )$ ($q^{*}\ne \pi$) for the longitudinal one. For the SDWI phase, both the transverse and longitudinal susceptibilities have the peak at $\mathit{q}=(\pi ,\pi )$. The density plot is normalized for clarification.

Figure 6

(a) Heisenberg exchange interaction $J$, DM interaction $D$, and bond-dependent Ising-type exchange interaction $\sqrt{J^2+D^2}-J$ as functions of $\theta$. (b) Density plot of the lowest eigenvalue of $F\left(\mathit{q}\right)$ as a function of $\mathit{q}$ for $\theta =\pi /6$, $\pi /2$, and $5\pi /6$, with its values normalized to $[0:1]$. (c) Ordering wave vector ${\mathit{Q}}_{\mathrm{LT}}$ at $q_x\ge q_y\ge 0$ as a function of $\theta$ obtained by the Luttinger-Tisza method (solid line): those in the spiral-$Q$ region $\theta <{\theta }_c^{\left(\mathrm{LT}\right)}$ have ${\mathit{m}}_{\mathit{r}}<1$. Data points are obtained as the mean-field solution, and show the spiral $2Q$ at $\pi /2<\theta <{\theta }_c^{\left(\mathrm{LT}\right)}$: a small plateau and other regions have different $Q$ but they are both the spiral $2Q$. (d) Lowest-order quantum energy correction $\delta E_{\mathrm{MF}}$ from the linear spin-wave theory. (e) Energy difference between the stripe ($\zeta =0$) and vortex ($\zeta =\pi /4$) phases.

Figure 7

(a) Schematic illustration of the model with the SSD. (b) System-size dependence of the average magnetization for the AFM and metallic phase at $\theta =\pi /3$ and $U=0.5$. (c) Structure factor as the function of $\mathit{q}$ obtained from Eq. (45) together with the PBC result along the symmetric line in the Brillouin zone. (d) System-size dependence of the peak position at $\theta =\pi /6$ and $U=3.0$. (e) Energy difference between the SDWII and metallic phase of the same $\theta$ as a function of $U$ (left) and $U-U_c$ (right panel) for several choices of $\theta$, where $U_c$ is the phase boundary.

Figure 8

Schematic illustration of the algorithm of the DMET. (b) How the one-body potential develops during the DMET iteration for $\theta =\pi /6$ and $U=2.0$.

Figure 9

(a) Gauge-invariant quantities, $W_{\mathrm{f}}$, ${\mathrm{\Phi }}_{\mathrm{f}}$, $W_{\mathrm{R}}$, ${\mathrm{\Phi }}_{\mathrm{R}}$, as functions of $\theta$, together with the magnetic phase diagram obtained by the SSDMF at $U\sim 8$. The red line indicates $\theta$ where $W_{\mathrm{f}}=0$ or $W_{\mathrm{R}}=0$. The inset shows the four site cluster we consider. (b) Schematic illustration of the local SU(2) gauge transformation defined on four sublattices A (red), B (yellow), C (green), and D (blue). The local SU(2) gauge transformation rotates the spin quantization axis by $\pm \pi /2$ about the gray arrows, followed by the rotation about the $z$ axis.

Figure 10

(a) Relationships between the approximations used in this study. The Luttinger-Tisza method is verified at around $U=\infty$, and the RPA at $U\lesssim U_c$, which give consistent results with the SSDMF results. (b) Summary of the results of the previous calculations on the $\pi$-flux Hubbard model which is identical to the $\theta =\pi$ case of our model after we perform the local SU(2) gauge transformation. The AFM(market with $*$) of the $\pi$-flux Hubbard model is the vortex phase of our model. The large-scale QMC results of larger sizes all predict the absence of the QSL phase ($U_{c1}$ being absent) [108, 109, 110, 111], showing that Chang et al. [112] using $L\le 14$ suffer a numerical finite-size effect.