BCS-BEC crossover in a $(t_{2g})^{4}$ Excitonic Magnet

Excitonic condensation has attracted considerable attention for over half a century, since its early theoretical prediction [1]. An exciton is a bound state of an electron-hole pair. This composite particle resembles the Cooper pair of superconductors and follows hard-core bosonic statistics. Early work involving semiconductors showed that in the weak-coupling limit (small $U$), near the semimetal to semiconductor transition and depending on the excitonic binding energy and the band gap, the system can become unstable against the formation of multiple excitons [2, 3] and this can lead to a condensation into a BCS-like macroscopic state called Excitonic Insulator. This BCS wavefunction smoothly transforms into a BEC state, if the gap between the bands increases at fixed Hubbard repulsion $U$. The phenomenon of excitonic condensation in strongly correlated systems can also lead to a BEC-like Excitonic Insulator in the strong coupling limit (large $U$). This regime also attracted considerable theoretical investigations  [4, 5, 6] and has been studied using perturbation theory at large $U$, where the hopping amplitude $t$ is the small parameter.

The extended Falicov-Kimball model has been used often as a minimal theoretical model to investigate the above discussed physics [7, 8, 9]. Only recently, the more realistic, but also more difficult, three-orbital Hubbard models were explored to study excitonic condensation [10, 11, 12, 13]. Due to the theoretical progress, on the experimental side there have been some studies showing reliable signatures of the excitonic condensate in real materials such as in a transition metal dichalcogenide [14], Ta₂NiSe₅ [15], and in bilayer systems [16, 17]. To increase our understanding of this interesting Excitonic Insulator it is important to find additional candidate materials and additional theoretical models where the excitonic condensation occurs and can be studied in detail. Recently, it was recognized that materials with strong spin-orbit coupling can also provide a new avenue for the study of excitons and excitonic condensation [18, 19]. For example, Sr₂IrO₄, with Ir⁴⁺ ions and filling $n=5$ electrons per Iridium ion ($t_{2g}^{5}$), is a celebrated material due to the presence of long-range antiferromagnetic ordering in quasi two-dimensional layers, as in the case of superconducting cuprates [20]. Recent resonant inelastic X-ray scattering (RIXS) experiments on Sr₂IrO₄ have reported the presence of excitons as an excitation at approximately $0.5-0.6$ eV [21, 22, 23]. RIXS experiments on one-dimensional stripes of Sr₂IrO₄ have also shown the presence of excitons at nearly the same energy [24]. These excitons are present in spin-orbit entangled states, hence known as spin-orbit excitons.

Another promising avenue to study excitonic condensates induced by spin-orbit coupling involves the $t_{2g}^{4}$ case, which is the focus of the present paper. Theoretically, it has been predicted that the three-orbital Hubbard model with a degenerate $t_{2g}$ space and in the $LS$ coupling limit ($U\gg t,\lambda$) can lead to the BEC state of triplons (singlet-triplet excitations) [25, 26, 27]. We will show that these triplon excitations are low-energy manifestations of spin-orbit excitons. Cluster Dynamical Mean Field Theory (DMFT) calculations have also supported similar findings [10, 11], although it is difficult to distinguish between BCS and BEC states in the small clusters used. A recent computational study of the ground state of a one-dimensional spin-orbit coupled three-orbital Hubbard model in a non-degenerate (tetragonal) $t_{2g}$ space using the numerically-exact density matrix renormalization group (DMRG) also reported a phase with staggered spin-orbit excitonic correlations [12]. All the above mentioned studies reveal an antiferromagnetic ordering accompanying the excitonic condensate. The $t_{2g}^{4}$ case is relevant for materials with Ir⁵⁺ ions and other $4d/5d$ transition metal oxides with the same doping $n=4$. The presence of triplon condensation was initially discussed for double perovskite materials, such as Sr₂YIrO₆ and Ba₂YIrO₆ with a $5d^{4}$ configuration [28, 29, 30, 31, 32, 33], but later on RIXS experiments showed that the bandwidth of triplon excitations is not sufficiently large in comparison to $\lambda$ to develop condensation [34]. It should be remarked that recent RIXS experiments on Ca₂RuO₄ suggests that this compound could be a candidate material for excitonic magnetism [35], and $ab-initio$ calculations have reached the same conclusion [36].

To investigate the spin-orbit excitonic condensation, this paper uses a simple degenerate three-orbital Hubbard model. Using numerically exact DMRG simulations on one-dimensional chains, we show that even in the intermediate Hubbard repulsion regime ($U/W\approx 1$) an excitonic condensation is induced accompanied by antiferromagnetic ordering. This regime is stabilized by increasing sufficiently $\lambda$ starting at the incommensurate spin-density wave metallic region of $\lambda=0$. This numerical result of a spin-orbit excitonic condensate at intermediate $U/W$ cannot be understood using large $U$ perturbation theories. Moreover, in two-dimensional ($2d$) lattices, by using the Hartree-Fock Approximation (HFA), we have also found a similar excitonic insulator phase in both the weak and intermediate $U/W$ regimes. As our main result, we show that there is a BCS-BEC crossover inside the excitonic condensate phase, both in the $1d$ and $2d$ lattices we studied. The BCS limit of the excitonic insulator occurs at intermediate $U/W$ (and also for weak $U/W$ in $2d$), and by increasing $\lambda$ (at fixed $U/W$) the BEC state is reached. The previously known BEC state present at large $U/W$, due to the condensation of triplons, is here shown to be smoothly connected to the BCS excitonic insulator of intermediate and weak $U$. At strong coupling, in the non-magnetic phase at large $\lambda$, higher than needed for the excitonic magnetic phase, these electron-hole pair excitations transform into low-energy triplons ($\Delta J_{\textrm{eff}}=1$) and higher-energy quintuplons ($\Delta J_{\textrm{eff}}=2$) excitations and those triplons condense on decreasing $\lambda$. We provide numerical evidence for this triplon condensation using DMRG by calculating the excitonic pair-pair susceptibility. We also provide the full $\lambda$ vs $U$ phase diagrams for $1d$ and $2d$ lattices using the many-body techniques discussed above.

The organization of this paper is as follows. In Sec. II, the model used and the computational methodology are presented. We discuss the atomic limit of the Hamiltonian in Sec. III. In Sec. IV, the main results are shown. In particular, in Sec. IVA, the results on $1d$ lattices using the DMRG technique are presented. In Sec. IVB, the results on $2d$ lattices using HFA are displayed to further support our main proposal of a BCS-BEC crossover in the model studied. In Sec. V, we discuss the overall results and present our conclusions.

For the study presented in this publication, we used a degenerate three-orbital Hubbard model. The Hamiltonian has three primary terms: the tight-binding kinetic energy, the standard on-site multiorbital Hubbard interaction, and the on-site spin-orbit coupling (SOC): $H=H_{K}+H_{\mathrm{int}}+H_{SOC}$. The tight-binding term is written as

The hopping amplitudes $t_{\gamma\gamma^{\prime}}$ connect only nearest-neighbor lattice sites (in both the $1d$ chain and $2d$ square lattices). As discussed earlier, here we use degenerate orbitals, that is, $t_{\gamma\gamma^{\prime}}=\delta_{\gamma\gamma^{\prime}}t$. In some cases we fixed $t=0.5$ for simplicity but most of our results are expressed in terms of dimensionless ratios, such as $U/W$ or $\lambda/t$. The total bandwidth is $W=4.0|t|$ and $8.0|t|$ for the $1d$ and $2d$ lattices, respectively. The above mentioned orbitals – labeled as 0, 1, and 2 – could be associated respectively to the $d_{yz}$, $d_{xz}$, and $d_{xy}$ orbitals, namely the $t_{2g}$ sector. The on-site multiorbital Hubbard interaction term in the Hamiltonian consists of the standard several components

In the equation above, $\mathbf{S}_{{i}\gamma}={{1}\over{2}}\sum_{\alpha,\beta}c_{{i}\alpha\gamma}^{\dagger}\sigma_{\alpha\beta}c^{\phantom{\dagger}}_{{i}\beta\gamma}$ represents the spin at orbital $\gamma$ and lattice site ${i}$, while $n_{{i}\gamma}$ is the electronic density operator also at orbital $\gamma$ and site $i$. The first two terms are the intra- and inter-orbital electronic repulsion, respectively. The Hund coupling is contained in the third term, which favors the ferromagnetic alignment of the spins at different orbitals and the same site. Finally, the $P_{{i}\gamma}=c_{{i}\downarrow\gamma}c_{{i}\uparrow\gamma}$ operator in the fourth term (pair-hopping) is the pair anhilation operator that arises from the matrix elements of the “1/r” Coulomb repulsion as in the early studies of Kanamori. We used the standard relation $U^{\prime}=U-2J_{H}$ due to rotational invariance, and we fix $J_{H}=U/4$ for all the calculations as employed widely in previous efforts [12, 13, 37, 38, 39].

The SOC term is

where the coupling $\lambda$ is the SOC strength.

Using the model described above, we performed calculations on one-dimensional chains employing the DMRG technique [40, 41] for various system lengths, such as $L$ = 16, 24, and 32 sites. We have used up to 1600 states for the DMRG process and have maintained a truncation error below $10^{-8}$ throughout the finite algorithm sweeps. We used the corrected single-site DMRG algorithm [42] with correction $a$ = 0.001-0.008, and performed 35 to 40 finite sweeps to gain proper convergence depending on the system size. After this convergence, we calculated the spin-structure factor $S(q)$, orbital-structure factor $L(q)$, excitonic momentum distribution function $\Delta_{m=1/2}(q)$, and coherence length $r_{coh}$. With the information gathered from all these observables, we constructed the phase diagram. To calculate spectral functions with the DMRG, we have used the correction vector method [43], with the Krylov-space approach [44]. We obtain the single-particle spectral function $A(q,\omega)$ and the excitonic pair-pair susceptibility $\Delta_{m=1/2}(q,\omega)$. These frequency-dependent observables require heavy computational time, and multiple compute nodes. In our DMRG calculations, we target the total $J_{\textrm{eff}}^{z}$ of the system to reduce the computational cost [12]. Details of the Hartree-Fock calculations are discussed in Sec. IVB below.

Before describing our results for chains and planes, we will discuss in detail the atomic limit of our Hamiltonian to introduce to the readers particular aspects of the $t_{2g}^{4}$ system. The magnetic properties in the $n=4$ case with a finite spin-orbit coupling are fascinating because in the atomic limit the ground state is a singlet having $J_{\textrm{eff}}=0$ for any finite value of $\lambda/U$. Only at $\lambda=0$, the $J_{\textrm{eff}}=0,1,2$ states are degenerate in the ground state manifold. Figure 1(a) shows the energies of the excited states relative to the $J_{\textrm{eff}}=0$ ground state. The evolution of magnetic moments and occupations in the single-particle spin-orbit ($j_{\textrm{eff}},m$) states is displayed in Fig. 1(b). Note that $j_{\textrm{eff}}$ is the effective total angular momentum of the electron and $m$ is the projection along the $z$-axis. Sometimes the quantum number $j_{\textrm{eff}}$ will be denoted as $j$, if it appears in subindex.

The $jj$ coupling limit ($\lambda/U>>1$) can be understood by diagonalizing the spin-orbit coupling term. The transformation between the $t_{2g}$ orbitals and the $j_{\textrm{eff}}$ basis is given by (site $i$ index dropped)

where $s$ is $1(-1)$ when $\sigma$ is $\uparrow(\downarrow)$ and $\bar{\sigma}=-\sigma$. The $H_{SOC}$ term in the $j_{\textrm{eff}}$ basis becomes

The expression above clearly shows that for $n=4$ the ground state has all $j_{\textrm{eff}}=3/2$ states fully occupied, and all $j_{\textrm{eff}}=1/2$ states empty (i.e. $|GS\rangle_{jj}=a_{\frac{3}{2},\frac{-3}{2}}^{\dagger}a_{\frac{3}{2},\frac{3}{2}}^{\dagger}a_{\frac{3}{2},\frac{-1}{2}}^{\dagger}a_{\frac{3}{2},\frac{1}{2}}^{\dagger}|0\rangle$) leading to a total $J_{\textrm{eff}}=0$. The first excited state is located with a gap $3\lambda/2$, as shown in Fig. 1(a). Also it can be checked that ${}_{jj}<GS|\mathbb{L}^{2}|GS>_{jj}=_{jj}<GS|\mathbb{S}^{2}|GS>_{jj}=4/3$. On the other hand, in the $LS$ coupling limit only one-third fraction of the ($j_{\textrm{eff}}=1/2,m$) states is occupied, with both the magnitudes of S and L being 1.

The information presented above in the atomic limit is relevant for cases with $U/W\ll 1$ and a finite spin-orbit coupling, to find out in which limit the system is located. For example, recently compounds with isolated RuCl₆ octahedra exhibited a single-ion physics with a non-magnetic $J_{\textrm{eff}}=0$ state [45]. Next, we will discuss our results for the $1d$ and $2d$ lattices where the presence of kinetic energy can drive the system away from this atomic limit, leading to magnetic ordering.

This section presents the DMRG numerical results for one-dimensional chains, and discuss their implications. Figure 2 visually summarizes our conclusions. Within our numerical accuracy, we observed that the excitonic condensation starts at intermediate Hubbard correlation strength $U\approx\mathcal{O}(W)$. This condensation of excitons opens a gap in the $j_{\textrm{eff}}=3/2$ and $j_{\textrm{eff}}=1/2$ bands at momentum $q=0$ and $q=\pi$, respectively, as shown in Fig. 2(a). A similar perspective for the gap opening near the Fermi level by excitonic condensation was discussed in early research for semiconductors [1]. At intermediate $U$, we noticed the excitons condense at finite momentum $\pi$ and in the triplet channel (for details see Supplementary [46]). We also noticed that increasing $\lambda$ (concomitantly adjusting $U$ to remain inside the excitonic condensation region) decreases the coherence length ($r_{coh}$) of electron-hole pairs from approximately one lattice spacing ($a$) to a much smaller number $r_{coh}<<a$, resembling the BCS-BEC crossover. Although in the extreme BCS limit $r_{coh}$ can be as large as hundreds of lattice spacings, outside our resolution, in our finite and one dimensional system we only found a robust indication for a maximum $r_{coh}$ of approximately $1.0a$ which definitely is different from the BEC limit $r_{coh}<<a$. Confirming that indeed we found a BCS-BEC crossover at intermediate $U$, we performed mean-field calculations on $2d$ lattices (Sec. V), where we found $r_{coh}$ as large as $\mathcal{O}(10a)$ in the BCS limit.

We have also found clear distinctions between the exciton-exciton correlation decay in the BCS and BEC limits in one-dimensional chains. It can be shown that the excitonic operator $\Delta_{1/2,m}^{\dagger 3/2,m^{{}^{\prime}}}$ in the single-atom $LS$ coupling limit can be written in terms of triplon and quintuplon excitations (see Supplementary [46]), corresponding, respectively, to $\mathbb{T}_{n}^{\dagger}|J_{\textrm{eff}}^{z}=0,\mathbb{J_{\textrm{eff}}}=0\rangle=|n,1\rangle$, and $\mathbb{Q}_{l}^{\dagger}|J_{\textrm{eff}}^{z}=0,\mathbb{J_{\textrm{eff}}}=0\rangle=|l,2\rangle$, where $n\in\{\pm 1,0\}$ and $l\in\{\pm 2,\pm 1,0\}$. By calculating the pair-pair susceptibity of excitons $\Delta(q,\omega)$ for one-dimensional chains, in the non-magnetic phase present in the strong coupling limit, we found bands of triplon and quintuplon excitations, with minima at $q=\pi$ but both bands are gapped. We noticed that decreasing $\lambda$ drives the system to the BEC state, where $\Delta(q,\omega)$ has clear signatures of gapped triplon and quintuplon excitations only near $q=0$ but mainly a continuum of excitations at other momenta above the Goldstone-like modes emerging from $q=\pi$, as sketched in Fig. 2(c).

In this section, we will present and discuss the results obtained in two-dimensional lattices by performing mean-field calculations in real-space. All the four-fermionic terms in the Hubbard interaction Eq.(2) are treated under the Hartree-Fock approximation, where the single-particle density matrix expectation values $\langle c_{i\alpha\sigma}^{\dagger}c_{i\beta\sigma^{\prime}}\rangle$ serve as order parameters at each site $i$ ($\alpha,\beta$ are orbitals and $\sigma,\sigma^{\prime}$ are spin projection indexes). We reach self-consistency iteratively in the order parameters while tuning the chemical potential $\mu$ accordingly to attain the required electronic density. Given the many order parameters, converged results often require using 10-20 initial random initial configurations and inspecting the lowest energy achieved after the iterative process, searching for patterns that are then uniformized and test for their energy. Often also converged results at other couplings are used as seeds at new couplings, until a consistent phase diagram emerges. The modified Broyden’s method was used to gain fast convergence [50].

In Fig. 8(a), we show the full $\lambda$ vs $U$ phase diagram for the two-dimensional lattice. To identify the various phases, we calculated the spin structure factor and orbital structure factor on 16$\times$16 clusters with periodic boundary conditions for all the points explicitly shown in the phase diagram. Our mean-field calculations in two dimensions capture almost all the phases found in our numerical exact one-dimensional results, the main difference being having shifted boundaries which are to be expected considering the different dimensionality and different many-body approximations. The only notable difference is that our mean-field calculations do not capture the non-magnetic phase in the strong coupling limit, because the lattice non-magnetic state can be written as a direct product of $J_{\textrm{eff}}=0$ at each site i.e. … $|J_{\textrm{eff}=0}\rangle_{i}\otimes|J_{\textrm{eff}=0}\rangle_{i+1}\otimes|J_{\textrm{eff}=0}\rangle_{i+2}...$, where each $J_{\textrm{eff}}=0$ state in the $LS$ coupling limit is a sum of Slater determinants and hence cannot be reproduced by the Hartree-Fock approximation that relies on single determinants. However, the excitonic condensate phase, the focus of the present paper, is properly captured by the Hartree-Fock calculations and it is present even in a larger region of the phase diagram than in one dimension, hence giving us a good opportunity to discuss the presence of the BCS-BEC crossover in two-dimensional lattices.

To proceed with our discussion, we fix $U/W=0.5$ ($W=8t$) where the excitonic condensation is present in a narrow but finite range of spin-orbit coupling, while for smaller $\lambda$’s the IC-SDW phase is present. Similar to our one-dimensional DMRG calculations, in two-dimensional lattices we found that inside the excitonic condensate region (at a fixed weak or intermediate $U/W$ values), $r_{coh}$ decreases on increasing $\lambda$, as shown in Fig. 8(d) depicting the BCS-BEC crossover. We have calculated $r_{coh}$ for system sizes 8$\times$8, 16$\times$16, and 24$\times$24. We found that in the BCS limit, $r_{coh}$ increases as the system size increases and $r_{coh}$ can reach values as high as $\approx 15.0a$ for the 24$\times$24 lattice. This clearly supports our claim for the presence of the BCS state above the IC-SDW region, as in our DMRG chain calculations (but in one-dimension perhaps due to size-effects or lack of resolution, the largest $r_{coh}$ obtained was only of order of one lattice spacing). On the other hand, we found that in the BEC limit, located below the relativistic band insulator in the phase diagram, $r_{coh}$ is $\mathcal{O}(0.1a)$, as shown in Fig. 8(d). Figure 8(e) displays $S(\pi,\pi)$ and $L(\pi,\pi)$, for $U/W=0.5$, to show that only for a finite range of $\lambda$ the antiferromagnetic ordering develops. We also show the momentum distribution function of excitons $\Delta_{1/2}(\mathbb{q})$ at $\lambda/W=0.53$ and $\lambda/W=0.592$ in the BCS and BEC limits, respectively. We noticed that in the BEC limit $\Delta(\mathbb{q})$ is much sharper near $\mathbb{q}=(\pi,\pi)$ than in the BCS limit, as expected because in BEC a larger ratio of excitons is expected at the condensation momentum than other momenta. To further investigate the above claims we calculated the ratio of excitons at wavevector $\mathbb{q}=(\pi,\pi)$ and at other wavevectors using $N(\pi,\pi)=\Delta_{1/2}(\pi,\pi)/\langle\Delta_{1/2}(\mathbb{q}\neq(\pi,\pi))\rangle$, where $\langle\Delta_{1/2}(\mathbb{q}\neq(\pi,\pi))\rangle=\frac{1}{L^{2}-1}\sum_{\mathbb{q}\neq(\pi,\pi)}\Delta_{1/2}(\mathbb{q})$. It is evident from Fig. 8(e) that $N(\pi,\pi)$ increases as we transition from the BCS to the BEC limits.

The single-particle density of states (DOS) provides further evidence for the existence of the excitonic condensate. Figure  9(a) and (b) show the DOS $\rho_{j}(\omega-\mu)$, at fixed $U/W=0.5$. We chose $\lambda/W=0.02$ and $\lambda/W=0.55$ corresponding to the IC-SDW and excitonic condensate phases, respectively. As shown in Fig. 9(b), at $\lambda/W=0.02$ a finite density-of-states at the chemical potential is present for both bands $j_{\textrm{eff}}=3/2$ and $j_{\textrm{eff}}=1/2$ indicating that the system is metallic in the IC-SDW phase. In the excitonic condensate phase, a small gap is present near the chemical potential in both the bands due to the condensation of excitons. The continous distribution in eigenenergies was attained by solving 24$\times$24 lattices with 48$\times$48 twisted boundary conditions [51]. Two peaks arising from the small hole pocket of the $j_{\textrm{eff}}=1/2$ states and electron pocket ot the $j_{\textrm{eff}}=3/2$ states are also visible. The inset shows the magnified $\rho_{j}(\omega-\mu)$ near the chemical potential. To discuss the evolution of these hole and electron pockets, as the system crossovers from the BCS to BEC limit, we show the momentum distribution function of electrons $n_{jm}({\bf k})$ in Fig. 9(c). The hole pocket is present at $\mathbb{k}=(0,0)$ in the $j_{\textrm{eff}}=3/2$ band and the electron pocket is at $\mathbb{k}=(\pi,\pi)$ in the $j_{\textrm{eff}}=1/2$ band. The nesting vector between these pockets also explains the ordering momentum of the excitons. As $\lambda$ increases, we noticed a gradual decrease in both electron and hole pockets which indicates that the number of carriers decreases as the system crossovers from the BCS to BEC limits.

In this publication, we studied the degenerate $t_{2g}^{4}$ multiorbital Hubbard model in the presence of spin-orbit coupling, using one-dimensional chains and numerically exact DMRG and also using two-dimensional clusters within the Hartree-Fock approximation. In both calculations, we provide evidence for a BCS-BEC crossover in the spin-orbit excitonic condensate, in the regime where the Hund coupling $J_{H}$ is fixed at a robust value $J_{H}=U/4$. Within our accuracy, we established that in this model and at intermediate $U/W$, the system transits from an IC-SDW metallic phase to the BCS limit of an antiferromagnetic excitonic condensate, and on further increasing $\lambda$ the coherence length of electron-hole pairs decreases rapidly as the system crossovers to the BEC regime. This BEC regime ends as eventually the system transits to the relativistic band insulator on increasing further $\lambda$. Our work provides the first indications of a BCS-BEC crossover in the excitonic magnetic state at intermediate $U/W$, a region of couplings that cannot be explored within approximations developed in the large $U/W$ regime.

We hope our study encourages further theoretical and experimental investigations on $t_{2g}^{4}$ coumpounds with robust spin-orbit coupling. Although our study is performed using degenerate orbitals, we believe that our findings could be generic and relevant for materials showing magnetic excitonic condensation at intermediate values of the Hubbard repulsion and spin-orbit coupling.

N.K., R.S., and E.D. were supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), Materials Sciences and Engineering Division. A.N. was supported by the Canada First Research Excellence Fund. G.A. was supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. DOE, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences, Division of Materials Sciences and Engineering. Part of this work was conducted at the Center for Nanophase Materials Sciences, sponsored by the Scientific User Facilities Division (SUFD), BES, DOE, under contract with UT-Battelle.