Spin texture in a bilayer high-temperature cuprate superconductor

To describe the locally non-centrosymmetric bilayer high-$T_{c}$ cuprate, we use the following tight-binding bilayer Rashba model Gotlieb et al. (2018); Raines et al. (2019); Hitomi and Yanase (2019); Harrison et al. (2015):

The non-interacting part includes three terms,

where $c_{l\mathbf{k}\sigma}(c_{l\mathbf{k}\sigma}^{\dagger})$ is the annihilation (creation) operator of an electron on the $l$th layer ($l=1,2$) with spin $\sigma={\uparrow},\downarrow$ and wave vector $\mathbf{k}$. The dispersion relation on a square lattice is $\varepsilon(\mathbf{k})=-2t(\cos k_{x}+\cos k_{y})+4t^{\prime}\cos k_{x}\cos k_{y}-2t^{\prime\prime}(\cos 2k_{x}+\cos 2k_{y})-\mu$, in which the hopping terms up to third nearest neighbor ($t$, $t^{\prime}$, and $t^{\prime\prime}$) and the chemical potential $\mu$ are included. $\mathbf{g}_{\mathbf{k}}=(-\sin k_{y},\sin k_{x},0)$ defines the antisymmetric SOC of Rashba type and $\bm{\sigma}$ is the vector of Pauli matrices. Due to the global inversion symmetry, the SOC in layers 1 ad 2 are opposite in sign: $\lambda_{1}=-\lambda_{2}$. The interlayer hoping is $t_{\perp}(\mathbf{k})=t_{z}(\cos k_{x}-\cos k_{y})^{2}$, which causes the bilayer splitting in high-$T_{c}$ cuprate. For Bi2212, the nearest neighbor hopping is $t=360$ meV Markiewicz et al. (2005); Drozdov et al. (2018). In the remainder of this paper, we set $t=1$ as the energy unit, and choose the other tight-binding parameters to be $t^{\prime}=-0.3$, $t^{\prime\prime}=0.1$ and $t_{z}=0.08$ Drozdov et al. (2018). Finally, we set $\lambda_{1}=-\lambda_{2}=0.06$.

The interacting part of Hamiltonian reads

in which $n_{\mathbf{r},l,\sigma}$ is the number of electrons of spin $\sigma$ at lattice site $\mathbf{r}$ of layer $l$. Only on-site interactions are considered.

In order to reveal the spin texture arising in model (1), we use cluster dynamical mean-field theory (CDMFT) Lichtenstein and Katsnelson (2000); Kotliar et al. (2001); Liebsch et al. (2008); Sénéchal (2015) with an exact diagonalization solver at zero temperature (or ED-CDMFT). In this approach, the infinite lattice is tiled into identical units (or supercells) defining a superlattice. The supercell is made of one or more clusters, each of which coupled to a bath of uncorrelated, auxiliary orbitals. The parameters describing this bath (energy levels, hybridization, etc.) are then found by imposing a self-consistency condition.

In this work the supercell is made of two superimposed, four-site plaquettes (one per layer), each of which coupled to a bath of eight uncorrelated orbitals. The cluster-bath system, or impurity model, is illustrated on Fig. 1 and defined by the following Anderson impurity model (AIM):

where $H_{c}$ is the Hamiltonian (1), but restricted to a single cluster, and $c_{\mu}$ and $a_{\alpha}$ destroy electrons on the cluster sites and the bath orbitals, respectively. Probing superconductivity forces us to use the Nambu formalism, in which each degree of freedom is occurring in particle and hole form in a multiplet. Thus, the index $\mu$ is a composite index comprising cluster site $i$, spin and Nambu indices: $c_{i}=(c^{\phantom{\dagger}}_{i{\uparrow}},c^{\phantom{\dagger}}_{i\downarrow},c^{\dagger}_{i{\uparrow}},c^{\dagger}_{i\downarrow})$. This index takes $4\times 4=16$ values in the particular AIM that we use. Likewise, the index $\alpha$ comprises bath orbital index $r$, spin and Nambu indices and takes $4\times 8=32$ values: $a_{r}=(a^{\phantom{\dagger}}_{r{\uparrow}},a^{\phantom{\dagger}}_{r\downarrow},a^{\dagger}_{r{\uparrow}},a^{\dagger}_{r\downarrow})$. $\theta_{\mu\alpha}$ is a complex-valued, $16\times 32$ hybridization matrix between cluster and bath orbitals, whereas $\epsilon_{\alpha\beta}$ is a $32\times 32$ matrix of one-body terms within the bath, including possible superconducting pairing. In principle, the matrix $\epsilon_{\alpha\beta}$ could be diagonalized (this would change the values of the hybridizations $\theta_{\mu\alpha}$), but we find it convenient and intuitive to allow pairing operators between bath orbitals.

The bath parameters are assumed to be spin independent, since we are not looking for magnetic ordering. In order to probe superconductivity, we include singlet and triplet pairing operators within the bath. Given two bath orbitals labeled by $r$ and $s$, the following pairing operators are defined:

In terms of the bath orbital numbering scheme defined on Fig. 1b, the pairing terms in $H_{\rm imp}$ are

These terms are included in the one-body matrix $\epsilon_{\alpha\beta}$.

The AIM is characterized by 10 variational parameters, all illustrated on Fig. 1b: bath energy levels $\epsilon_{1,2}$ (diagonal elements of $\epsilon_{\alpha\beta}$), hybridization amplitudes $\theta_{1,2}$, singlet pairing amplitudes $\Delta_{1,2}$ and triplet pairing amplitudes $p_{x1}$, $p_{x2}$, $p_{y1}$, $p_{y2}$. It turns out that, owing to the rather small value of $\lambda_{1,2}$ in Eq. (3), the triplet bath parameters are too small to have an observable effect and can be neglected. This reduces the number of independent bath parameters to six. The two clusters forming the supercell also happen to have the same bath parameter values in the converged solutions, which is expected from symmetry.

For a given set of bath parameters, the AIM (6) may be solved and the electron Green function computed. The latter may be expressed as

where $\mathbf{t}_{c}$ is the one-body matrix in the cluster part of the impurity Hamiltonian $H_{\rm imp}$, $\bm{\Sigma}_{c}(\omega)$ is the associated self-energy, and $\bm{\Gamma}_{c}(\omega)$ is the bath hybridization matrix:

where $\bm{\theta}$ is the $16\times 32$ matrix with components $\theta_{\mu\alpha}$ and $\bm{\epsilon}$ the $32\times 32$ matrix with components $\epsilon_{\alpha\beta}$. Equation (10) allows us to extract the cluster self-energy $\bm{\Sigma}_{c}(\omega)$ from computed quantities. The fundamental approximation of CDMFT is to replace the full self-energy of the problem with the local self-energy $\bm{\Sigma}$. More precisely, when the supercell contains more than one cluster, as is the case here, the supercell self-energy is the direct sum of the self-energies of the different clusters: $\bm{\Sigma}(\omega)=\bigoplus_{c}\bm{\Sigma}_{c}(\omega)$. The lattice Green function is then approximated as

where $\mathbf{\tilde{k}}$ is a wave vector of the reduced Brillouin zone (associated with the superlattice) and $\mathbf{t}(\mathbf{\tilde{k}})$ is the one-body Hamiltonian (2)-(4) expressed in that mixed basis of reduced wave vectors and supercell orbitals. In our system, the matrix $\mathbf{G}(\mathbf{\tilde{k}},\omega)$ has dimension $2\times 16=32$, because of the two clusters forming the supercell.

Let us finally summarize the self-consistent procedure used to set the bath parameters, as proposed initially in Caffarel and Krauth (1994): (i) trial values of the bath parameters are chosen on the first iteration. (ii) For each iteration, the AIM (6) is solved, i.e., the cluster Green functions $\mathbf{G}_{c}(\omega)$ are computed using the Lanczos method, for each cluster. (iii) The bath parameters are updated, by minimizing the distance function:

where $\bar{\mathbf{G}}_{c}(\omega)$ is the restriction to cluster $c$ of the projected Green function $\bar{\mathbf{G}}$, defined as

(iv) We go back to step (ii) and iterate until the bath parameters or the bath hybridization functions $\bm{\Gamma}_{c}(\omega)$ stop varying within some preset tolerance.

Ideally, $\bar{\mathbf{G}}_{c}(\omega)$ should coincide with the impurity Green function $\mathbf{G}_{c}(\omega)$, but the finite number of bath parameters does not allow for this correspondence at all frequencies. This is why a distance function $d(\bm{\epsilon},\bm{\theta})$ is defined, with emphasis on low frequencies along the imaginary axis. The weight function $W(i\omega_{n})$ is where the method has some arbitrariness; in this work $W(i\omega_{n})$ is taken to be a constant for all Matsubara frequencies lower than a cutoff $\omega_{c}=2t$, with a fictitious temperature $\beta^{-1}=t/50$.

The lattice Green function (12) can be used to compute the average of any one-body operator defined on the lattice. In addition, we can go back to a fully wave vector-dependent representation $\mathcal{G}(\mathbf{k},\omega)$, where $\mathbf{k}$ now belongs to the original Brillouin zone and $\mathcal{G}$ is a smaller, $8\times 8$ matrix, by a procedure called periodization. The simplest periodization scheme is to Fourier transform $\mathbf{G}$ directly from the supercell to the original Brillouin zone, as follows Sénéchal et al. (2000):

where $i,j$ are composite spin, Nambu and layer indices, and the difference between $\mathbf{k}$ and $\mathbf{\tilde{k}}$ is an element $\mathbf{K}$ of the reciprocal superlattice: $\mathbf{k}=\mathbf{\tilde{k}}+\mathbf{K}$. Note that $\mathbf{G}(\mathbf{k},\omega)=\mathbf{G}(\mathbf{\tilde{k}},\omega)$ since $\mathbf{G}$ is by construction a periodic function of the reduced Brillouin zone.

Figure 2 shows the $d$-wave order parameter, computed from the Green function (12), as a function of hole doping, for $U=8$. This is the ground state average of the following operator:

where $\mathbf{x}$ and $\mathbf{y}$ denote the nearest-neighbor vectors on the square lattice. The blue curve is obtained in a single-layer model, without spin orbit coupling. The black squares are obtained in the current bilayer model, and differ very little from the single layer values, because of the small value of both $t_{z}$ and $\lambda_{1,2}$. The green diamonds are the average of the following triplet operator:

Note the factor of 100 in the scale. A similar operator defined along the $y$ axis with the $x$ component of the triplet $d$-vector has equal expectation values.

The periodized Green function (15) give us access to quantities observable by SARPES, such as the spectral function:

where $\operatorname{Tr}_{N}$ means a trace over spin and layer indices that excludes the Nambu sector. One can also extract spin information by projecting the Green function $\bm{\mathcal{G}}$ on various spin directions. We thus define the following spin spectral functions:

One can also define the corresponding layer-resolved quantities.

Figure 3 shows the spin texture measured in three of the superconducting solutions obtained, for three different values of density $n$ (underdoped, optimally doped and overdoped), on the first layer. The left panel shows a color plot of the spectral function (18) at the Fermi level (the non-interacting Fermi surface is indicated by a dashed line). The right panel shows the magnitude of the projection on the $x$-$y$ plane of the spin spectral function (19), on which we superimposed a (diluted) vector plot indicating the direction of the spin in the $x$-$y$ plane, as given by the vector $\mathbf{S}_{\bot}(\mathbf{k})=(S_{x}(\mathbf{k},0),S_{y}(\mathbf{k},0))$. On this plot the arrows only provide a direction, not a magnitude. The color scale is the same for all three values of electron density in the figure. The structure of the spin texture is similar for all three cases illustrated: it is made of a clockwise rotating pattern around $\mathbf{k}=(0,0)$ and around $\mathbf{k}=(\pi,\pi)$, on the first layer. The corresponding quantities on the second layer are obtained simply by reversing the arrows.

As expected, the amplitude of the spin texture is maximal in the vicinity of the non-interacting Fermi surface, with the important proviso that the pseudogap phenomenon suppresses the spectral weight (and the amplitude of the spin texture) away from the diagonals as we go in the underdoped regime. Exactly on the Fermi surface, at least near the diagonals, the spin texture vanishes as it must change direction, leading to split maxima along the diagonals on the right panels of Fig. 3. Figure 4 shows the drop of the spectral weight (indeed its maximum value over the Brillouin zone at the Fermi level) as doping $x$ decreases (blue circles). This suppression is even more pronounced for the maximum value of $|\mathbf{S}_{\bot}|$ (red squares).

Overall, the known physics of the two-dimensional, one-band Hubbard model is not affected by the presence of the spin-orbit coupling between the two layers or by the inter-layer coupling. The $d$-wave order parameter (Fig. 2) is not affected in any visible way and the triplet component of the order parameter is two orders of magnitude smaller than the singlet component. The spin-orbit coupling and the associated spin texture do not interfere with the pseudogap physics illustrated by the concentration of the spectral weight along the diagonal as one gets closer to half-filling. The spin-orbit coupling is too small to make the system topological Nakosai et al. (2012), and we checked that the Chern number, computed by adding the real part of the zero-frequency self-energy to the non-interacting Hamiltonian Wang and Zhang (2012a, b), is indeed zero in all the cases studied. The amplitude of the spin texture, however, is more suppressed in the underdoped region than the spectral function is (roughly twice as much).