Theoretical studies of a three-band Hubbard model with a strong spin-orbit coupling for 5d transition metal oxide Sr2IrO4

Motivated by recent experiments on Sr2IrO4, we study the ground state properties of a two-dimensional three-band Hubbard model with a strong relativistic spin-orbit coupling. Using the exact diagonalization technique, the dynamical magnetic structure factor M(q, ω) is calculated to examine the low-energy magnetic excitations. We find that the low-energy excitations in M(q, ω) are well described by an effective Heisenberg model composed of a local Kramers doublet of an effective total angular momentum Jeff = |s⃗ − ℓ⃗| = 1/2. The antiferromagnetic exchange interaction estimated from M(q, ω) is as large as ~ 80 meV, which is in good quantitative agreement with experiments. To study a possible long-range ordered state in the thermodynamic limit, we use the variational cluster approximation based on the self-energy functional theory, which is parallelized to accelerate the calculations. We find the ground state where the local Kramers doublet is in-plane antiferromagnetic ally ordered.

1. Introduction 5d transition metal oxides in a layered perovskite structure such as Sr 2 IrO 4 [1] and Ba 2 IrO 4 [2] have attracted much attention because of their unique properties caused by a strong relativistic spin-orbit coupling (SOC) of 5d transition element. It has been observed experimentally that these materials behave like an effective total angular momentum J eff = | s − ℓ | = 1/2 Mott insulator at low temperatures [2,3,4], where s ( ℓ) is the spin (orbital) angular momentum operator. Moreover, very recent experiments using muon spin rotation (µSR), magnetic susceptibility, and resonant inelastic x-ray scattering (RIXS) have revealed that the low-energy magnetic excitations can be described by an "isospin"-1/2 Heisenberg model with an effective exchange interaction as large as ∼ 60-100 meV [2,5,6], which is almost comparable to the ones in the parent compounds of high-T c cuprate superconductors [7].
Motivated by these experiments, we study theoretically the nature of magnetic excitations for Sr 2 IrO 4 using a three-band Hubbard model with the SOC. For this purpose, here the numerically exact diagonalization technique is employed. We find that the low-energy magnetic excitations are described by an antiferromagnetic Heisenberg model with the nearest-neighbor exchange coupling as large as ∼ 80 meV. Furthermore, to study the stability of this antiferromagnetically ordered state, we use the variational cluster approximation (VCA) [8] based on the self-energy where Here, c † r,α,σ (c r,α,σ ) represents the creation (annihilation) operator of an electron with spin σ = (↑, ↓) and orbital α = (xy, yz, zx) on site r, n r,α,σ = c † r,α,σ c r,α,σ , c † k,α,σ is the Fourier transform of c † r,α,σ , andσ indicates the opposite spin of σ. The number N e (N h ) of electrons (holes) is set to be 5N (N ), where N is the number of sites. ε α k is the energy dispersion of orbital α defined by where ∆ is the energy level difference between d xy and the others [10,11]. The relativistic SOC term is described by H so in (3), where φ rασ is the Wannier orbital with spin σ and orbital α at site r. We assume the matrix elements φ rασ | ℓ · s φ r ′ βσ ′ are finite only for r = r ′ . H so is then written in the matrix form where s ↑ = +1 (s ↓ = −1). Equation (4) describes the Coulomb interaction term which includes intra-orbital Coulomb repulsion U , inter-orbital Coulomb repulsion U ′ , Hund's coupling J, and pair hopping term J ′ with U = U ′ + 2J and J = J ′ [12].  [11] is determined to reproduce the energy band dispersion calculated by the first-principles electronic structure calculations based on the density functional theory [10]. In this paper, we set U/t 1 = 6.111 and J/U = 0.1. These values are similar to the ones estimated by the constrained random phase approximation [13,14].

Exact diagonalization study on magnetic excitations
In this section, the magnetic excitations are studied using the numerically exact diagonalization technique. Because of the exponential increase of the Hilbert space with N , the system sizes that can be treated are limited. However, this method provides the unbiased results, which are valuable and complementary to other approximate calculations such as VCA.

Method
We employ the numerically exact diagonalization technique for a cluster of √ 8× √ 8 with periodic boundary conditions. To reduce the Hilbert space dimensions, we fully use the symmetry of Hamiltonian, i.e., translational symmetry and spin inversion symmetry. We first use the Lanczos algorithm [15] to calculate the ground state |ψ 0 and the ground state energy E 0 . Then, the continued fraction expansion [15] is used to calculate the dynamical correlation function for an operator O, which is defined by where η(> 0) is a broadening factor. Note that the static correlation function is related to the dynamical correlation function by S O = dωC O (ω). Following the notation in (9), the dynamical magnetic structure factor is given by is the Fourier transform of the ξ (= x, y, z) component of the local spin (orbital) angular momentum operator S ξ r (L ξ r ), andσ ξ is the ξ component of the Pauli matrix.
quasi-particle operator, we construct a pseudo-spin operator J ξ θ (r): Then, we calculate the dynamical correlation function for J ξ θ (q), i.e., J ξ is the Fourier transform of J ξ θ (r). The results for J ξ θ (q, ω) are shown in the right panels of figure 1. Comparing these two spectra M x (q, ω) and J x θ (q, ω) in figure  1, we find that the position of the low-lying peak at each q in M x (q, ω) are exactly the same as the one in J x θ (q, ω) [18]. This is the direct numerical evidence that the low-energy magnetic excitations correspond nicely to the excitations of the local Kramers doublet.
Assuming that the low-energy excitation dispersions are expressed by the magnon dispersion ω k of an isotropic "isospin"-1/2 Heisenberg model, i.e., ω k = 2J 1 − (cos k x + cos k y ) 2 /4, we can estimate an effective exchange interaction J. As shown in figure 1, we find that J is as large as ∼ 0.22t 1 (= 79 meV). This value is in good agreement with experimental observations [2,6,5].

Variational Cluster Approximation (VCA) study on the long-range ordered ground state
In the previous section, we have found that the low-energy magnetic excitation is well described by the effective antiferromagnetic Heisenberg model composed of the local Kramers doublet. However, the method used in Section 3 is not suitable to study properties in the thermodynamic limit. To access the thermodynamic limit and study a possible long-rage ordered state, here we employ an approximation technique called VCA.

Method
The VCA [8] based on the SFT [9] is one of cluster approaches to study interacting fermion systems on lattices. In this method, we first divide the original system into small clusters, such that each of them can be treated by the exact diagonalization technique. This set of small clusters is called a "reference system" below. In order to apply the VCA, the interacting part in the reference system must be the same as that in the original system, which is certainly satisfied in (1). However, the non-interacting part such as the kinetic energy term and the SOC term can be different. Here, we shall denote by t ′ the parameter set of the non-interacting part of the reference system. Within the VCA, the grand potential Ω of the original system per site is written by the following equation: where Ω t ′ is the exact ground potential of the reference system andĜ t ′ (z) is the Green's function of a decoupled cluster of the reference system. The matrix elements ofĜ t ′ (z) are simply given by (15) where a is the shorthand notation of (r s , α, σ) and r s is the site position in a decoupled cluster of the reference system. H t ′ is the Hamiltonian in a decoupled cluster of the reference system with the ground state |ψ t ′ and the ground state energy E t ′ . The chemical potential for the reference system is denoted by µ ′ .Ĝ t ′ (z) can be calculated using the continued fraction method. The matrix elements ofV (K) is given by where N c is the number of clusters, A is the shorthand notation of (r c + r s , α, σ), r c is the position of the cluster, and t A,A ′ (t ′ A,A ′ ) indicates the coefficient of the term c † r,α,σ c r ′ ,α ′ ,σ ′ in the non-interacting part of the Hamiltonian for the original (reference) system, which includes the hopping parameter t i (i = 1-5), the SOC λ, and the potential difference ∆. K is the momentum for the superlattice of clusters defined by r c .
The choice of the parameter set t ′ is arbitrary and, in principle, the best result should be obtained by searching for the optimized t ′ which satisfies the stationary condition for Ω. However, computing directly Ω is impractical. For this reason, it is more convenient to start with an appropriate Weiss field. In previous section, we found that the low-energy magnetic excitations is similar to the ones of the effective Heisenberg model composed of the effective local Kramers doublet. Therefore, it is natural to introduce the following Weiss field for the reference system: with Q = (π, π). The variational parameters in t ′ treated in this study are thus h ′ , θ ′ , and φ ′ , in addition to chemical potential µ ′ . For the optimization of multi-variables (h ′ ,θ ′ , φ ′ ,µ ′ ), we use the direction set method by optimizing one parameter at a time (with the remaining parameters fixed) iteratively until the full optimization is achieved for all parameters. The Brent's method is used for the optimization [19]. The most time consuming part in the VCA is the calculation for all the matrix elements of G(z). Because this calculation is independent, it is readily parallelized by using MPI without any trouble. We also parallelize the integral on x and summation over K in (14). Figure 2 (a) shows the total elapsed time (t e ) versus the total number N core of cores used to optimize two parameters, µ ′ and θ ′ , with a fixed system size. It is clearly seen in figure 2 (a) that t e monotonically decrease with increasing N core up to N core ∼ 300 without saturation. We define the following quantity r s to evaluate the efficiency on the strong scaling of the parallelization: where N * core is a reference number of cores and t * e is the total elapsed time when N * core number of cores are used. Thus, r s = 1 corresponds to 100% efficiency and usually r s is less than one. As shown in fig. 2 (b), we find that our numerical scheme keeps r s > 0.8 up to N core ∼ 300. The reference number of cores is N * core = 8. Figure 3 shows the variational parameter h ′ dependence of the ground state energy E(h ′ ) obtained by the VCA using a 4-site cluster. The optimal h ′ is found to be h ′ = 0.1892t 1 . A finite value of the optimal h ′ means that the system prefers to be symmetry broken. For all calculations, we find φ ′ = 0, which indicates the effective Kramers doublet prefers to order antiferromagnetically with easy xy-plane anisotropy. This anisotropy is caused by the presence of the SOC and the Hund's coupling J. This finding is consistent with other theoretical results obtained by a strong-coupling theory [16], a 2-site exact diagonalization study [20], and a variational Monte Carlo study [11]. We also find in the VCA calculation that θ determined by minimizing the hole density n rθσ is θ ∼ 0.35π, which is close the the value for the exact J eff = 1/2 limit, i.e., θ = arctan √ 2 ∼ 0.304π.

Conclusion
Using the numerically exact diagonalization technique, we have studied the magnetic excitations of the three-band Hubbard model with the strong SOC for Sr 2 IrO 4 . We have found that the lowenergy magnetic excitations are well described by those of the J eff = 1/2 Kramers doublet, and that the lowest excitation dispersions can be reproduced by the antiferromagnetic Heisenberg model. We have estimated the effective exchange interaction as large as ∼ 80 meV, which is in good agreement with recent experimental observations. We have also discussed the stability of the antiferromagnetically ordered ground state by using the VCA. We have found the ground state where the local Kramers doublet is in-plane antiferromagnetically ordered. These results strongly support the SOC induced J eff = 1/2 Mott insulator for Sr 2 IrO 4 .