Phase diagram of strongly-coupled Rashba systems

We first consider the noncentrosymmetric monolayer, in point group $C_{4v}$, which is described by the Hamiltonian

	$\displaystyle H=$	$\displaystyle H_{\text{int}}+H_{0},$		(1)
	$\displaystyle H_{0}=$	$\displaystyle\sum_{\langle i,j\rangle}\sum_{s,s^{\prime}}-tc^{\dagger}_{i,s}(\sigma^{0})_{ss^{\prime}}c_{j,s^{\prime}}+i\lambda c^{\dagger}_{i,s}\hat{r}_{ij}\cdot\left(\sigma^{x}\hat{y}-\sigma^{y}\hat{x}\right)_{s,s^{\prime}}c_{j,s^{\prime}}$		(2)
	$\displaystyle H_{\text{int}}=$	$\displaystyle U\sum_{i}n_{i,\uparrow}n_{i,\downarrow}.$		(3)

The normal-state Hamiltonian $H_{0}$ is characterized by a nearest-neighbour hopping $t$ on a square lattice and a Rashba SOC of strength $\lambda$, with the latter allowed by the broken inversion symmetry. In (2) the $\sigma^{\mu}$ matrices denote the Pauli spin matrices, $\hat{r}_{ij}$ is the unit vector between sites $i$ and $j$, and $\langle i,j\rangle$ denotes summation over nearest neighbours only. The interaction term is a Hubbard on-site repulsion of strength $U$. Here $c_{j,s}$ and $c^{\dagger}_{j,s}$ have the usual meaning of annihilation and creation operators for spin-$s$ electrons on site $j$, and $n_{j,s}$ is the corresponding number operator.

We analyze the Hamiltonian (1) in the strong coupling limit, where $U$ far exceeds the bandwidth of the normal-state Hamiltonian. We adopt the approximation that this excludes doubly-occupied sites except as virtual processes. Projected onto the subspace excluding double occupancy, we have the effective Hamiltonian

$\displaystyle H_{\text{eff}}$

$\displaystyle=\tilde{H_{0}}+H_{J}$

(4)

where $\tilde{H_{0}}$ is as in (2) with the electron operators replaced by their projections in the no-double-occupation space and $H_{J}$ is obtained via canonical transformation [17]:

$\displaystyle\begin{split}H_{J}=\sum_{\langle i,j\rangle}&J\left[\mathbf{S}_{i}\cdot\mathbf{S}_{j}-\frac{n_{i}n_{j}}{4}\right]+D(\hat{\mathbf{z}}\times\mathbf{\hat{e}_{ij}})\cdot(\mathbf{S}_{i}\times\mathbf{S}_{j})\\ +&J^{\prime}\left[\alpha_{|j-i|}\left(S_{i}^{x}S_{j}^{x}-S_{i}^{y}S_{j}^{y}\right)-S_{i}^{z}S_{j}^{z}-\frac{n_{i}n_{j}}{4}\right]\end{split}$

(5)

where $J=4t^{2}/U$ is the Heisenberg coupling, $J^{\prime}=-4\lambda^{2}/U$ the compass interaction and $D=8t\lambda/U$ the Dzyaloshinskii-Moriya (DM) coupling. Note that the compass interaction includes the direction-dependent factor $\alpha_{|j-i|}$ where $\alpha_{x}=-1$ and $\alpha_{y}=1$. The DM interaction depends non-trivially on the hopping direction, $\hat{\mathbf{e}}_{ij}$.

We use the auxiliary boson method [18, 19] to develop a mean-field (MF) theory for the effective Hamiltonian. The electron operators are replaced as $c^{\dagger}_{i,\sigma}\rightarrow f^{\dagger}_{i,\sigma}b_{i}$, where the $f_{\sigma}$ and $b$ are annihilation operators for the fermionic spinons and bosonic holons, respectively. At the MF level, the no-double-occupancy condition is enforced on the average with

$\displaystyle 1=\sum_{i,\sigma}\langle f^{\dagger}_{i\sigma}f_{i\sigma}\rangle+\langle b^{\dagger}_{i}b_{i}\rangle=\sum_{i,\sigma}\langle f^{\dagger}_{i\sigma}f_{i\sigma}\rangle+\delta.$

(6)

where the density of holons is equal to the doping, i.e. $\langle b^{\dagger}_{i}b_{i}\rangle=\delta$ [20]. Following the usual approach, we assume that the holons are condensed at low temperatures [21, 22], and so we can use the MF approximation $\langle b_{i}\rangle=\sqrt{\delta}$. This auxiliary boson technique arises as a natural MF treatment of strongly coupled systems, which are well described by a resonating valence bond (RVB) state background [23, 24, 25].

We proceed to apply the MF approximation to the projected Hamiltonian as follows: we replace the electron annihilation and creation operators by the spinon and holon operators, with the holon operators further replaced by their expectation value, e.g. $b_{i}^{\dagger}b_{j}\approx\langle b_{i}^{\dagger}\rangle\langle b_{j}\rangle=\delta$. We further decouple the nearest-neighbour interactions in the particle-particle and particle-hole channels, introducing the MF amplitudes

	$\displaystyle\chi^{\mu}_{a}$	$\displaystyle=\langle f_{is_{1}}^{\dagger}(\sigma^{\mu})_{s_{1}s_{2}}f_{i+as_{2}}\rangle$		(7)
	$\displaystyle\Delta^{\mu}_{a}$	$\displaystyle=\langle f_{i+as_{1}}(i\sigma^{\mu}\sigma^{y})_{s_{1}s_{2}}f_{is_{2}}\rangle$		(8)

The normal-state amplitudes $\chi$ and $\chi^{\text{SOC}}$ are defined as

	$\displaystyle\chi$	$\displaystyle=\chi_{\pm x}^{0}=\chi_{\pm y}^{0}$		(9)
	$\displaystyle\chi^{\text{SOC}}$	$\displaystyle=\pm\chi^{y}_{\pm x}=\mp\chi^{x}_{\pm y}.$		(10)

The renormalization or bond parameters $\chi$ and $\chi^{\text{SOC}}$ are a measure of the probability of spin-preserving and spin-flipping hopping between neighbouring sites, respectively. This can alternately be thought of as singlet bonds and opposite-spin triplet bonds; e.g., an electron undergoing a spin-preserving hopping requires an opposite spin partner in the neighbouring site, effectively creating a singlet bond. The amplitudes $\chi$ and $\chi^{\text{SOC}}$ renormalize the nearest-neighbour hopping and the Rashba SOC in the noninteracting Hamiltonian, respectively. Doping away from the half-filled RVB background state and adding SOC allows spin mobility and spin-flip processes, and hence singlet/triplet mixed-state superconductivity and multiple bond parameters.

For now we restrict our attention for the gap amplitudes $\Delta^{\mu}_{a}$ to the $A_{1}$ and $B_{1}$ irreducible representations (irreps) of the point group $C_{4v}$, as these contain singlet pairing states which are favoured by the nearest-neighbour interaction in (5); the pairing amplitudes for other irreps are small or vanishing. For the $A_{1}$ irrep, the gap amplitudes are

	$\displaystyle\Delta^{s}$	$\displaystyle=\Delta^{0}_{\pm x}=\Delta^{0}_{\pm y}$		(11)
	$\displaystyle\Delta^{t}$	$\displaystyle=\pm\Delta^{y}_{\pm x}=\mp\Delta^{x}_{\pm y}$		(12)

while for the $B_{1}$ irrep we have

	$\displaystyle\Delta^{s}$	$\displaystyle=\Delta^{0}_{\pm x}=-\Delta^{0}_{\pm y}$		(13)
	$\displaystyle\Delta^{t}$	$\displaystyle=\pm\Delta^{y}_{\pm x}=\pm\Delta^{x}_{\pm y}.$		(14)

The pairing potentials in momentum space are listed in table 1. Note that the singlet $A_{1}$ component corresponds to an extended $s$-wave pairing, whereas the singlet $B_{1}$ is $d_{x^{2}-y^{2}}$-wave. The MF amplitudes are determined by solving the self-consistency equations, and the stable superconducting state is determined from the free energy. In our numerical calculations we utilize a lattice of $200\times 200$ ${\bf k}$-points, and the temperature is taken to be $T=0.001J$.

A monolayer model similar to the one considered here was studied in [26] using the Gutzwiller approximation. These authors considered relatively weak spin-orbit coupling, and our results in this limit are in agreement with theirs.

In figure 1 we show the variation of the MF parameters as a function of hole doping and SOC strength. Our main result here is that there is a transition between the $B_{1}$ and $A_{1}$ superconducting states with increasing SOC strength, and the transition line moves to higher SOC strengths with increasing doping. In both superconducting states the singlet pairing amplitude dominates over the triplet.

The transition from the $d$-wave-dominated $B_{1}$ to extended-$s$-wave-dominated $A_{1}$ pairing state can be understood in terms of the lifting of the spin-degeneracy by the SOC, which results in two Fermi surfaces (FS) in the normal state (NS) as illustrated in figure 2. In the absence of SOC and weak doping, the FS passes close to the gap maxima of the $B_{1}$ state; in contrast, the FS lies close to the nodal line of the $A_{1}$ extended-$s$-wave state. Thus, the average gap opened by the $B_{1}$ state over the FS is much larger than that opened by the $A_{1}$, which energetically favours the former. Upon switching on the SOC, however, the spin-split FSs will tend to move away from the nodal line of the $A_{1}$ state, thus enhancing its stability relative to the $B_{1}$ state, and eventually ensuring that it has the lower free energy. As we increase the doping, the Fermi surface moves off the nodal line. Although the extended $s$-wave state can now open a full gap, it remains less stable than the $d$-wave. Moreover, it requires a larger SOC to stabilize, since the gap on one of the spin-split Fermi surfaces initially decreases as the SOC is switched on; only after it crosses the nodal line can the SOC stabilize the extended $s$-wave.

As for the subdominant triplet gaps, figure 1 shows that the triplet amplitude is immediately enhanced upon the transition into the $A_{1}$ state. This likely reflects the fact that the $A_{1}$ triplet state is insensitive to the SOC, whereas the SOC is pair-breaking for the $B_{1}$ state, i.e. the $A_{1}$ triplet state only pairs quasiparticles in the same spin-split (“helicity”) bands. Although a $B_{1}$ triplet state with same-helicity pairing is possible, this requires $f$-wave pair amplitude and thus interactions beyond nearest neighbour [27, 28].

As discussed in [26], the monolayer NCS hs nontrivial topological properties. We have not considered the topological character of the phases in our system, but this is potentially a fruitful direction for further work.

The bilayer system can be considered as two copies of the monolayer with opposite sign of the SOC, which are coupled by an interlayer hopping (ILH) term $t_{\perp}$ [7, 4, 14, 12]. The layer degree of freedom restores the inversion symmetry and takes the system to the point group $D_{4h}$, and we encode this degree of freedom in the $\eta^{\mu}$ Pauli matrices. The non-interacting part of the Hamiltonian is thus

	$\displaystyle H_{0}=$	$\displaystyle\sum_{\langle i,j\rangle}\sum_{s,s^{\prime},l,l^{\prime}}-tc^{\dagger}_{i,s,l}\left(\eta^{0}_{ll^{\prime}}\sigma^{0}_{ss^{\prime}}\right)c_{j,s^{\prime},l^{\prime}}+i\lambda c^{\dagger}_{i,s,l}\eta^{z}_{ll^{\prime}}\hat{r}_{ij}\cdot\left(\sigma^{x}_{ss^{\prime}}\hat{y}-\sigma^{y}_{ss^{\prime}}\hat{x}\right)c_{j,s^{\prime},l^{\prime}}$
		$\displaystyle+t_{\perp}\sum_{j}\sum_{s,s^{\prime},l,l^{\prime}}c^{\dagger}_{j,s,l}\left(\eta^{x}_{ll^{\prime}}\sigma^{0}_{ss^{\prime}}\right)c_{j,s^{\prime},l^{\prime}}$		(15)

where $c_{j,s,l}$ is the annihilation operator for a spin-$s$ electron on layer $l=1,2$ of unit cell $j$. As in the monolayer case we include an on-site Hubbard interaction, which we assume to be the dominant energy scale. Accordingly, we develop an effective theory excluding double-occupancy at each site of the bilayer. This gives effective nearest-neighbour interactions: in each layer this has the same form as (5), but with opposite sign of the DM interaction. The ILH term gives rise to an additional Heisenberg interaction with exchange constant $J_{\perp}=4t_{\perp}^{2}/U$ between the spins in each layer of the same unit cell.

The additional layer degree of freedom also increases the number of MFs, with both intralayer

	$\displaystyle\chi^{\mu}_{a,l}$	$\displaystyle=\langle f_{i,s_{1},l}^{\dagger}(\sigma^{\mu})_{s_{1}s_{2}}f_{i+a,s_{2},l}\rangle$		(16)
	$\displaystyle\Delta^{\mu}_{a,l}$	$\displaystyle=\langle f_{i+a,s_{1},l}(i\sigma^{\mu}\sigma^{y})_{s_{1}s_{2}}f_{i,s_{2},l}\rangle$		(17)

and interlayer amplitudes

	$\displaystyle\chi_{\perp}$	$\displaystyle=\langle f^{\dagger}_{isl}(\eta^{x}\otimes\sigma^{0})_{ss^{\prime}ll^{\prime}}f_{is^{\prime}l^{\prime}}\rangle$		(18)
	$\displaystyle\Delta_{\perp}$	$\displaystyle=\langle f_{isl}(\eta^{x}\otimes\sigma^{0}i\sigma^{y})_{ss^{\prime}ll^{\prime}}f_{is^{\prime}l^{\prime}}\rangle$		(19)

The normal-state intralayer amplitudes have the in-plane variation as defined in (9) and (10), but $\chi^{SOC}$ changes sign across the layers. For the pairing amplitudes, we focus on the irreps of the $D_{4h}$ point group which have nearest-neighbour spin-singlet pairing, namely the $A_{1g}$, $A_{2u}$, $B_{1g}$, and $B_{2u}$ irreps. The two intralayer gap amplitudes in the $A_{1g}$ irrep are

	$\displaystyle\Delta^{s}$	$\displaystyle=\Delta^{0}_{\pm x,l=1,2}=\Delta^{0}_{\pm y,l=1,2}$		(20)
	$\displaystyle\Delta^{t}$	$\displaystyle=\pm\Delta^{y}_{\pm x,1}=\mp\Delta^{x}_{\pm y,1}=\mp\Delta^{y}_{\pm x,2}=\pm\Delta^{x}_{\pm y,2}$		(21)

The intralayer amplitudes are essentially the same as those defining the $A_{1}$ irrep of the monolayer, and with $s$-wave-like dominant pairing and with the subdominant triplet amplitudes $\Delta^{t}<\Delta^{s}$ reversing sign between each layer. There is also an interlayer gap amplitude $\Delta_{\perp}$, which is vanishing in the other irreps we consider here. For the $A_{2u}$ irrep, the singlet gap amplitude reverses sign across the layers whilst the triplet gap amplitude does not. Explicitly, the intralayer amplitudes are defined

	$\displaystyle\Delta^{s}$	$\displaystyle=\Delta^{0}_{\pm x,1}=\Delta^{0}_{\pm y,1}=-\Delta^{0}_{\pm x,2}=-\Delta^{0}_{\pm y,2}$		(22)
	$\displaystyle\Delta^{t}$	$\displaystyle=\pm\Delta^{y}_{\pm x,l=1,2}=\mp\Delta^{x}_{\pm y,l=1,2}.$		(23)

Similarly, the $B_{1g}$ and $B_{2u}$ states have intralayer amplitudes which are the same as in the $B_{1}$ irrep of the monolayer, but where the $d$-wave-like dominant singlet and subdominant triplet components reverse sign between the layers, respectively. That is, for the $B_{1g}$ we have

	$\displaystyle\Delta^{s}$	$\displaystyle=\Delta^{0}_{\pm x,l=1,2}=-\Delta^{0}_{\pm y,l=1,2}$		(24)
	$\displaystyle\Delta^{t}$	$\displaystyle=\pm\Delta^{y}_{\pm x,1}=\pm\Delta^{x}_{\pm y,1}=\mp\Delta^{y}_{\pm x,2}=\mp\Delta^{x}_{\pm y,2}$		(25)

whereas for the $B_{2u}$

	$\displaystyle\Delta^{s}$	$\displaystyle=\Delta^{0}_{\pm x,1}=-\Delta^{0}_{\pm y,1}=-\Delta^{0}_{\pm x,2}=\Delta^{0}_{\pm y,2}$		(26)
	$\displaystyle\Delta^{t}$	$\displaystyle=\pm\Delta^{y}_{\pm x,l=1,2}=\pm\Delta^{x}_{\pm y,l=1,2}$		(27)

Table 2 shows the momentum-space form of the gap amplitudes for each irrep of the bilayer system. Note that since the interlayer pairing interaction acts only between sites in the same unit cell, the interlayer SC parameter is restricted to $s$-wave.

The phase diagram in $\lambda-t_{\perp}$ space can be seen in figure 3. The odd states were not found to be favoured at any point in the chosen parameter ranges. The intralayer renormalization parameters were almost independent of $t_{\perp}$ and the ILH renormalization was found to be almost independent of $\lambda$, all with almost imperceptible change across the $B_{1g}$ to $A_{1g}$ transition. As such only the SC gaps are shown.

The stabilizing effect of the interlayer pairing for the $A_{1g}$ irrep is evident for non-zero $t_{\perp}$. Increasing the ILH has the effect of shifting the $A_{1g}-B_{1g}$ transition line to lower values of $\lambda$. This appears to be driven primarily by the interlayer singlet in the $A_{1g}$, which reaches a comparable amplitude to the dominant intralayer singlet pairing above about $t_{\perp}>0.8t$ and intermediate SOC strength. For ILH strengths $t_{\perp}>t$, the $B_{1g}$ is completely suppressed in favour of the $A_{1g}$, as the ILH also pushes the FSs outward – similar to the effect the SOC had on the FSs which favoured the extended $s$-wave state in the monolayer.

Figure 4 shows the phase diagram in $\lambda-\delta$ space at fixed $t_{\perp}=0.2t$. Only the variation of the interlayer MFs are shown, since the intralayer MFs varying negligibly from the monolayer results. The transition between the $d$-wave-like ($B_{1g}$ and $B_{2u}$) and the extended-$s$-wave-like ($A_{1g}$ and $A_{2u}$) states is very similar to the monolayer $B_{1}$ to $A_{1}$ transition. At low doping the odd-parity $A_{2u}$ and $B_{2u}$ states are more stable than their even-parity counterparts. Since singlet pairing is dominant in our model, this is equivalent a transition from zero to $\pi$ phase difference between the layers. This is remarkable, as the weak coupling of the layers should favour the zero phase difference, in analogy to a Josephson junction.

To investigate how the odd states are stabilized, the momentum space-resolved free energy difference was obtained to elucidate where the largest contributions to the stability of the odd states lie in momentum space. This is shown in figure 5, where it is clear that the $A_{2u}$ is most stabilized about the $X$-point, whilst for the $B_{2u}$ this occurs about the point $k_{x}=k_{y}=\frac{\pi}{2}$.

For a $\hat{z}$-aligned external magnetic field $H_{Z}$, the Zeeman Hamiltonian is

$\displaystyle H_{Z}=\frac{-g\mu_{B}}{2}\sum_{kss^{\prime}l}H_{z}\sigma_{z}c^{\dagger}_{ksl}c_{ks^{\prime}l}.$

(38)

Applying this field adds an $A_{2g}$ symmetry distortion in the bilayer, the same irrep as the $z$-rotation. This lowers the symmetry of the system, taking it from $D_{4h}$ to $C_{4h}$. This allows the $A_{1u}$, $A_{2g}$, $B_{1u}$, $B_{2g}$ to mix with the $A_{2u}$, $A_{1g}$, $B_{2u}$, $B_{1g}$ irreps respectively (Table 2). The $\eta^{y}\sigma^{z}$ interlayer triplet in $A_{u}$ is sensitive only to the AFM interlayer coupling, with respect to which it is repulsive; as such it was found to have vanishing magnitude. Again, only the uniform extended $s$-wave dominant state ($A_{g}$) hosts an interlayer pairing. Of all new intralayer MFs introduced in table 2 only the additional triplet state in $A_{g}$ and $A_{u}$ was found to be nonzero. Nevertheless it is ignored in the following since the amplitude was negligible compared to the other MFs.