Evolution of electronic and magnetic properties of Sr2IrO4 under strain

Motivated by properties-controlling potential of the strain, we investigate strain dependence of structure, electronic, and magnetic properties of Sr2IrO4 using complementary theoretical tools: ab-initio calculations, analytical approaches (rigid octahedra picture, Slater-Koster integrals), and extended t−J\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t-{{{\mathcal{J}}}}$$\end{document} model. We find that strain affects both Ir-Ir distance and Ir-O-Ir angle, and the rigid octahedra picture is not relevant. Second, we find fundamentally different behavior for compressive and tensile strain. One remarkable feature is the formation of two subsets of bond- and orbital-dependent carriers, a compass-like model, under compression. This originates from the strain-induced renormalization of the Ir-O-Ir superexchange and O on-site energy. We also show that under compressive (tensile) strain, Fermi surface becomes highly dispersive (relatively flat). Already at a tensile strain of 1.5%, we observe spectral weight redistribution, with the low-energy band acquiring almost purely singlet character. These results can be directly compared with future experiments.


INTRODUCTION
Exploring the physics of quasi-two-dimensional (2D) spin-orbit Mott insulators can help to understand high-temperature superconductivity as well as the general interplay of spin-orbit coupling, Hund's, and Coulomb interactions. In particular, a lot of studies have been devoted to the quasi-2D iridates Sr 2 IrO 4 and Ba 2 IrO 4 1-3 . Iridates show eminent similarities to the cuprate family of hightemperature superconductors, both in structure and low-energy physics, and were expected to become superconducting upon doping. However, so far no superconductivity has been reported in iridates.
In general, Sr 2 IrO 4 behavior often deviates from theoretical predictions. For example, Mott insulators normally become metallic at high enough pressure as the unit cell becomes smaller and the bands broaden. This is also true for spin-orbit coupled Mott insulators, such as ruthenates 3,4 . In Sr 2 IrO 4 , resistance indeed decreases until the pressure of around 25-30 GPa (which according to ref. 5 , corresponds approximately to a strain of −4%) 5 , or, according to a very recent study 6 , 32-38 GPa (−5.1% strain). Then, however, resistance starts to increase, showing a peculiar U-shaped dependency and persisting insulating behavior up to at least 185 GPa 6 . So far, no metallization in Sr 2 IrO 4 or other iridates (Sr 3 Ir 2 O 7 7,8 , BaIrO 3 3,9 , etc) has been observed at pressures up to 40-185 GPa 3,6 . Moreover, there is also surprisingly little correlation between the insulating behavior and magnetism 10 as the latter disappears at around 20 GPa (roughly −2.9% strain) in Sr 2 IrO 4 5,11 and 14.4 GPa (roughly −2.1% strain) in Sr 3 Ir 2 O 7 8 , without the onset of a metal-insulator transition.
Furthermore, iridates emerge as a good functional playground for manipulation of the magnetic and electronic properties, which is an exciting goal both fundamentally and practically 3,12 . Iridiumbased heterostructures and superlattices have therefore emerged as a whole new field very recently [13][14][15][16][17][18][19][20][21] . Strain and pressure in particular are powerful tools on hand to control the magnetic properties of the material. It has been shown that misfit strain can directly control dispersion of magnetic excitations in Sr 2 IrO 4 22-24 , as well as transport properties 25 . A shift of the two-magnon Raman peak to higher energies was observed under tensile strain 24 , albeit much weaker than the shift observed in the canonical Mott-Hubbard insulator K 2 NiF 4 and cuprates like Bi 1.98 Sr 2.06 Y 0.68 Cu 2 O 8+δ 26 . In ref. 23 , the authors used resonant inelastic scattering (RIXS) to show that magnetic dispersion in Sr 2 IrO 4 is strongly affected by strain. In particular, the contribution of the second and third nearest-neighbor (NN) exchange was suppressed (enhanced) upon tensile (compressive) strain. The tensile strain was shown to drive the system closer to a shorter-range first-NN only Heisenberg limit, with only little magnon branch softening left at (π/2, π/2) already upon the tensile strain of +2%. Upon compressive strain, the energy of (π, 0) magnon was shown to increase 22,23 .
A clear understanding of the electronic and magnetic properties of iridates and their evolution with strain is, therefore, of interest not only from a fundamental point of view but also for applications 3,12 . Unveiling the details of the interplay of lattice, magnetic, and other degrees of freedom in Sr 2 IrO 4 is needed to understand the recently observed electrical control of octahedra rotation 27,28 and the much-debated strong magnetoelastic coupling [29][30][31][32] . Currently, a clear understanding of neither how exactly nor by which mechanism do superexchange and hopping parameters in Sr 2 IrO 4 change with strain is available, not even on a phenomenological level. One of the interesting questions is whether the change in electronic and magnetic properties upon the strain is mostly associated with bond length change, as argued in e.g., ref. 22 , or the change of the in-plane rotation angles θ of the oxygen octahedra (see Fig. 1) 33 .
Studying the behavior of iridates under strain and pressure is a demanding task not only experimentally, but also theoretically. On one hand, iridates are strongly correlated Mott insulators 1 , so one needs to resort to theoretical methods where correlations are 1 treated non-perturbatively, employing effective descriptions like Hubbard or Heisenberg models. On the other hand, microscopic changes of orbitals, their overlap, and structural changes are essential to understand the behavior of a crystal under strain [34][35][36][37][38] , so ab-initio methods are demanded. Another difficulty is that as one eventually approaches a possible metal-insulator transition at high pressure and/or strain, effective models, such as the Heisenberg superexchange model, fail.
In this paper, we focus on the effect of strain and combine various complementary theoretical tools to provide a comprehensive analysis of how the magnetic properties are affected by strain. For different (compressive and tensile) strain values, we use density functional theory (DFT) based ab-initio calculations to access microscopic changes in the crystal structure, and study the corresponding changes in the electronic properties through Wannierization of the scalar-relativistic DFT bandstructure obtained within the generalized gradient approximation 39 . Subsequently, we solve an extended t À J model within the self-consistent Bohr approximation (SCBA) to obtain the angleresolved photoemission spectra (ARPES) and study the straincontrolled evolution of the Fermi surface. Realistic values of the input parameters for these calculations were used: the hopping parameters were obtained from the DFT calculations, while the extended-range exchange couplings were obtained by direct comparison to the magnon dispersion measured with RIXS. In this way, the presented analysis contains no free parameters apart from an overall constant energy shift (chemical potential) in the SCBA calculations.

RESULTS AND DISCUSSION
Evolution of hopping parameters under strain Sr 2 IrO 4 shows an in-plane staggered octahedra rotation characterized by a single parameter: ∡O − Ir − Ir angle denoted by θ in Fig. 1. Under ambient conditions, the octahedral rotation is found to be θ ≈ 13. 6 ∘ for the relaxed structure, which is close to the reported experimentally value of 11. 8 ∘40 . The epitaxial strain on iridates then affects not only the distance between the Ir atoms but also the Ir-O-Ir bond angle, as can be seen in Fig. 2a. The inplane octahedra rotation angle θ, obtained using DFT (see Methods for details), monotonically increases (decreases) upon compressive (tensile) strain in the studied range of −7.5% to 7.5%, where negative strains correspond to compression.
To ascertain the influence of structural changes on the electronic properties, we study the evolution of Wannier tight-binding model hoppings derived from DFT (see Methods for details), as a function of strain (Fig. 2b). The notations for the hoppings are shown in Fig. 1: the intraorbital hoppings between xy orbitals along a 0 or b 0 axes is denoted as t 1 , between xz(yz) along a 0 (b 0 ) axis as t 2 , and between xz(yz) along b 0 (a 0 ) as t 3 . The interorbital hopping between yz and xz orbitals is denoted as t 4 , all other interorbital hoppings are negligible. Further neighbor interorbital hoppings are denoted as t 0 and are shown in Fig. 1.
Upon compression, direction-dependent hopping parameter t 2 is increasing, but surprisingly, t 1 is decreasing (Fig. 2b). This emerging anisotropy in hopping parameters is interesting, as t 2 hopping describes the propagation of an electron with xz (yz) orbital character along only one axis, a 0 (b 0 ), whereas t 1 allows an xy electron to hop in both directions. We thus see that upon compressive strain, the system favors the separation of the entire Fermi sea into two Fermi seas with bond-dependent propagation (xz carriers which can only propagate along a 0 , and yz carriers which can only propagate along b 0 ) and suppression of the bondindependent and thus truly two-dimensional xy carriers. This compass-model-like 41 propagation is quite unusual and could cause the formation of charge density wave.
Upon tensile strain, t 1 is nearly independent of the strain value and is the dominant hopping, while t 2 decreases steadily (Fig. 2b). Different behavior of t 1 upon compressive and tensile strain reflects the change of Fermi surface topology between compressive and tensile strain.
It is also interesting to note that the smallest hopping parameter t 3 , describing the hopping between almost parallel dorbitals with very small overlap goes to zero around −3%, which corresponds to compression of~20 GPa, not too far from the value of resistivity minimum under pressure 5,6 .
To disentangle the contribution of inter-atomic distance d and the octahedral rotation θ to the hopping parameter trends with strain, we employ the analytical approaches of Glazer and Slater-Koster. The Glazer picture 42 is often used in rigid octahedra approximation whereby the main effect of the modest strain is assumed to be the change of the in-plane rotation angle θ. However, as detailed in Supplementary Fig. 2, the trends obtained within the Glazer picture disagree with the DFT results in Fig. 2b, and even contradict them in rigid octahedra approximation. Therefore, the Glazer picture has limited applicability for iridates, and rigid octahedra approximation is improper.
We then proceed with a more specific orbital-resolved Slater-Koster-integrals-based approach 43,44 . Slater-Koster integrals are hybridization matrix elements E between atomic d-states on neighboring atoms obtained via integrating over relevant spherical harmonics. The resulting interatomic matrix elements E are proportional to the d-wave functions overlap and can be expressed via cubic harmonic matrix elements V ddσ , V ddπ , V ddδ for a known bond direction l, m, n as tabulated in Slater-Koster tables 43,44 . In Sr 2 IrO 4 , we also need to account for the rotation of the d orbitals within the t 2g sector due to the in-plane octahedral rotation 45 . Therefore, we decompose the rotated d orbital on the basis of non-rotated d orbitals before evaluating the Slater-Koster matrix elements. For example, the hybridization matrix elementẼ between the two rotated NN xy orbitals can be obtained as a superposition of hybridization matrix elements E of non-rotated xy and x−y 2 orbitals obtained as: where θ is the in-plane rotation of the IrO 6 octahedra (see Fig. 1). Similarly, for the overlap between the rotated xz(yz) orbitals along the a 0 direction, we get E xz;xzðyz;yzÞ ¼ cos 2 ðθÞE xz;xzðyz;yzÞ À sin 2 ðθÞE yz;yzðxz;xzÞ ; and for interatomic interorbital overlap along x: Figure 3 shows the resulting hybridization matrix elementsẼ as a function of the in-plane octahedral rotation θ, the Ir-Ir distance d, as well as both the parameters (see Fig. 2a). We find that at least in the Slater-Koster approximation, accounting for the change of the distance d alone (Fig. 3b) can provide a better approximation to a full dependency of matrix elementsẼ on strain (Fig. 3c) then  accounting for the change of bond angle θ. This is also consistent with the quantum chemistry study 22 .
However, not all trends obtained from the DFT calculations are well reproduced: the hopping parameter E § †, § † is increasing under compressive strain (Fig. 3c), unlike the t 1 hopping extracted from DFT (Fig. 2b). To address this, we also consider the O-mediated indirect Ir−O−Ir hoppings.
The indirect oxygen-mediated overlap between the two rotated NN xy orbitals can be calculated as a sum of the hopping integrals between two Ir atoms via α = p x , p y orbitals of the oxygen, E xy;O;xy ¼ P α¼p x ;p y E xy;α;xy . The hopping integral is calculated as E xy;α;xy ¼ ðcosð2θÞE Àl;m;n α;xy þ sinð2θÞE Àl;m;n α;x 2 Ày 2 Þ ðcosð2θÞE l;m;n α;xy À sinð2θÞE l;m;n α;x 2 Ày 2 Þ=Δ pd ; where l ¼ cos θ, m ¼ sin θ, n = 0 are the directional cosines of the vector from the oxygen O to the Ir atom 43 along the a 0 -axis in the units of Ir-O distance d IrÀO ¼ 0:5d IrÀIr 0 = cos θ (d IrÀIr 0 is the distance between the NN Ir ions under ambient conditions), and E α,xy , E α,xy are the p-d Slater-Koster integrals 43 . We note that l, m, n indexes were omitted for Eqs. (1-3), because for the d−d overlap, Slater-Koster integrals are quadratic in directional cosines 43 , and simply {l, m, n} = {1, 0, 0} for a pair of Ir atoms along the a 0 bond. For the indirect oxygen-mediated overlap, however, one has to account for the sign of the directional cosines. Moreover, for the indirect hopping, the Ir-O-Ir hopping has to be renormalized by the charge transfer energy Δ pd ¼ E onÀsite xy À E onÀsite α , the energy difference between corresponding Ir-d and O-p orbitals (see Supplementary  Fig. 3a). Surprisingly, the charge transfer energy Δ pd has a strongly non-linear dependency on the strain (see Supplementary Fig. 3a).
We plot the resulting indirect superexchange matrix elementsẼ as a function of the in-plane octahedral rotation θ, Ir-Ir distance d in Fig. 4. Indeed, the indirect Ir-Ir overlap decreases drastically under compressive strain, unlike the direct d−d overlap. Interestingly, this behavior is directly linked to both the change in the distance between the atoms as well as the angle θ describing octahedra rotation (see Fig. 4a,b). We also note that taking into account the strain dependence of the iridium and oxygen on-site energies is crucial to obtain correct trends. In fact, nonlinear behavior of charge transfer energy Δ pd seems to be directly responsible for the non-linear strain dependency of the indirect oxygen-mediated hoppingẼ. This suggests that the role of oxygens in the low-energy physics of strained iridates and other transitional-metal oxides might be underestimated and requires further investigation.
Accordingly, the contribution of the indirect orbital overlap should be small for xz and yz orbitals. Indeed, the indirect yz-yz orbitals overlap along a 0 -axis is zero. The xz orbitals hybridize with O-p z orbital, however, this hybridization decreases with strain much slower than for xy orbitals, explaining the different behavior of xy and xz orbitals under compressive strain.
The fact that relative Ir-O hybridization is directly responsible for the resulting suppression of t 1 (xy-xy hopping) under compressive strain suggests an electronic state crossover as the role of xy orbitals in the composite J eff = 1/2 is decreasing. Notably, a pressure-induced phase transition was also suggested in a recent X-ray powder diffraction study 46 at pressures around 20 GPa, which should correspond to~−3% strain and is in good agreement with our findings.
Overlap of the spin-orbit coupled J eff states We now estimate the overlap between the J eff = 1/2 states for NN, 2NN, and 3NN (denoted τ, τ 0 and τ″, correspondingly), which can be calculated from t 2g orbitals overlap using the Clebsh-Gordon coefficients 2,47 . In Fig. 5, we show the change of the overlap of NN J eff = 1/2 states calculated from the DFT values obtained here. As experimental values of J eff = 1/2 states overlap τ are hard to measure, one can try to compare hopping parameters τ with available experimental estimates of magnetic exchange interactions J (see Table 1 in ref. 23 ), which in ambient conditions are assumed to scale with τ 22 . However, in a recent RIXS study on strained Sr 2 IrO 4 23 , the authors suggest that the simple J / τ 2 =U relationship fails for strained Sr 2 IrO 4 due to the polaronic renormalization of the charge excitations. In particular, the firstneighbor exchange interaction J 1 was shown to decrease slightly upon the tensile strain, while J 2 and J 3 decreased much faster, based on a fit of the Heisenberg model to the measured magnon dispersion 23 . An earlier RIXS study also suggested that magnetic exchange interaction J increases upon the compressive strain 22 .
Consistent with both RIXS studies 22,23 , calculated here values of NN, 2NN, and 3NN hopping parameters are all decreasing upon tensile strain (Fig. 5). As discussed in ref. 23 , this trend for τ's is significantly slower than that observed for superexchange interaction J , indicating that J / τ 2 =U relationship indeed fails for strained Sr 2 IrO 4 . Our DFT strain trends are also consistent with the modest increase of magnetic exchange interaction J 1 under compressive strain reported in the two-magnon Raman study 24 .
It is interesting to compare the trends observed in Sr 2 IrO 4 to those in 3d transition metal oxides-cuprates. In ref. 48 authors used XAS at Cu L 3 -edge of La 2 CuO 4 together with analytical and DFT theoretical approaches to show that both bandwidth and electron-electron correlations were increasing upon the compressive strain. As a result, magnetic exchange J was shown to increase (decrease) almost linearly upon the compressive (tensile) strain. In comparison, while the orbital-dependent hoppings in Sr 2 IrO 4 behave very differently from this (Fig. 2b), showing the surprising decrease in t 1 with compressive strain, effective J eff = 1/ 2 orbitals in Sr 2 IrO 4 have strain dependence (Fig. 5) somewhat similar to those of cuprates, echoing the famous parallel between Sr 2 IrO 4 and La 2 CuO 4 .

The evolution of the Fermi surface under strain
Structural response to strain may also be accompanied by changes in the Fermi surface. Thus, strong changes in Fermi surface upon uniaxial pressure have recently been reported in Rubased compound Sr 2 RuO 4 49,50 , along with the more than double increase of the superconducting transition temperature 51 . A recent work on Sr 2 IrO 4 employing a tight-binding model has pointed out that out-of-plane tilting of the oxygen octahedra can induce shrinking of the Fermi surface and suppress nesting and the d-wave superconductivity 52 . There is no out-of-plane tilting in Sr 2 IrO 4 under pristine conditions or modest strain-it was very recently shown to appear only under the pressure of as much as 40 GPa 6 . However, in-plane octahedra rotation is strongly affected by the modest strain already, and it is important to understand if and how the Fermi surface is affected, particularly on the beyondmean-field level.
To study the evolution of the Fermi surface under strain we calculate photoemission spectral functions of strained Sr 2 IrO 4 , using extended t À J model formalism developed in ref. 53 . The extended t À J model used in the calculation (see Methods) depends on two sets of parameters: the magnetic exchange parameters J 1 , J 2 , J 3 , Ising anisotropy coefficient Δ 54 , and the hopping parameters t i describing overlap of the t 2g orbitals. We obtain the set of t i 's for each strain value from DFT calculations as discussed in detail above. Using this Wannier Hamiltonian as a starting point describing single-particle hopping processes, we consider all possible many-body hopping processes to derive the hopping part of the extended t À J model 53 .
To properly account for the changes in the electronic structure, we need to account for the evolution of the magnetic exchange parameters with the strain. It has been obtained from the published 23 fits to the RIXS spectra on strained samples. As experimental data is available for small strain values range only, we restrict calculations of photoemission spectra to that range and show the photoemission (and inverse photoemission) spectra in Fig. 6 for two substrates: (LaAlO 3 ) 0 . 3 (Sr 2 AlTaO 6 ) 0 . 7 (100) (LSAT) and GdScO 3 (110) (GSO), providing a strain of −0.52% and +1.53%, correspondingly.
Calculated photoemission spectra show one conductance band at positive energies, and two valence bands at negative energies: a sharp singlet band around −0.25 eV and a more incoherent triplet band at −0.5 eV. (see Fig. 7 for a detailed discussion on the band character). We see that the strain-induced changes of the photoemission spectra are quite prominent for samples with a strain difference of 2%. First, for tensile strain, as compared to compressive strain, the Mott gap increases (Fig. 6b), suggesting stronger polaron binding of the photoinduced hole to the magnetic background 23 . Second, upon compressive strain, the photoemission spectra of Sr 2 IrO 4 show a highly dispersive singlet band (Fig. 6a), while upon tensile strain, both singlet and triplet bands are much less dispersive, and the Fermi surface of Sr 2 IrO 4 becomes relatively flat (Fig. 6b). It is important to note that the relative flattening of the Fermi sheet upon tensile strain is a many-body effect distinct from the anisotropic compass-like hoppings under compressive strain. One should be able to observe such significant renormalization of the spectral weight in the ARPES data even for small, realistic values of strain.
The conduction band is only weakly affected by strain. Figure 6 shows a marginal flattening of the conduction band upon the tensile strain. We, therefore, expect a minimal effect of epitaxial strain on possible superconductivity.
To explore the effect of strain on the ARPES spectra in more detail, we plot in Fig. 7 separate contributions of the singlet (J = 0) and triplet (J = 1) charge excitations to the full spectra. As one can see, under tensile strain, J = 0 contributes the most at (π,0), whereas under compressive strain, its contribution is slightly more widespread. On the contrary, the J = 1 spectral weight is shifted to (π/2, π/2) at the tensile strain and to (π, 0) at the compressive strain. In particular, J = 1 contribution to the "lower energy" band of the photoemission spectra is strongly reduced upon the tensile strain. We thus observe strain-controlled spectral weight redistribution between the charge carriers of singlet and triplet characters. Already moderate tensile strain is sufficient to make the lower energy band of almost purely singlet character. In summary, we predict a dramatic strain dependence of the electronic properties of Sr 2 IrO 4 for compressive v/s tensile strain. The most remarkable feature is the appearance of the compassmodel-like contribution of electron propagation due to the separation of the Fermi sea in Sr 2 IrO 4 into two subsets of bondand orbital-dependent carriers under compressive strain. This enables the formation of charge density wave due to nesting and could be connected to the puzzling metalization avoidance in Sr 2 IrO 4 upon pressure. The Fermi sea separation originates from strain dependency of relative Ir-O hybridization, as well as on-site O energy, suggesting an important role of oxygens in a lowenergy physics of strained iridates and other transitional-metal oxides.
Despite the suppression of the bond-independent xy hopping t 1 under compressive strain, the hopping amplitude of the composite J eff = 1/2 state still increases under compressive strain owing to the contribution from hoppings between xz/yz orbitals (t 2 )-the largest of two in-plane nearest-neighbor directiondependent hoppings. The obtained trends for J eff = 1/2 hopping, τ, are in good agreement with available experimental data.
We also calculated the photoemission spectra of Sr 2 IrO 4 upon compressive and tensile strain (for samples grown on LSAT and GSO substrates, respectively). We find that under compressive (tensile) strain, the singlet band becomes significantly more (less) dispersive, and both the singlet and triplet bands shift up (down) in energy. We also show that the electronic properties of the lowenergy model can be controlled by strain, since the already moderate tensile strain is sufficient to make the lower energy band of almost purely singlet character, and shift the triplet spectral weight to (π/2, π/2) point. These features can be readily  In-plane structure of Sr 2 IrO 4 and the notation for orbital overlap parameters. For nearest neighbors (NN), interorbital overlap parameters are denoted t 1 , t 2 , t 3 , and intraorbital t 4 , and for second nearest neighbors (2NN) we consider interorbital overlap only: t 0 , t 0 , t 0 . Orbital overlap between third neighbors is defined in the same way, not shown. Due to the symmetry considerations, the overlap of out-ofplane t 2g orbitals is anisotropic: t 2 describes hopping between xz(yz) along a 0 (b 0 ) axis. The O-Ir-Ir angle θ characterizes the in-plane octahedra rotation as shown.
E.M. Pärschke et al. observed in the future ARPES measurements-a smoking gun test of our findings.

DFT calculations
The DFT results are calculated by the Quantum Espresso package 55 using Perdew--Burke-Ernzerhof exchange-correlation functional 39 and projector-augmented-wave pseudopotentials 56 with the unit cell containing 16 Sr, 8 Ir, and 32 O atoms. In the calculations of structure relaxation and electronic structure, spin-orbit interaction is not included. The unstrained structure (a = 5.5098 Å, c = 26.1522 Å) is acquired through relaxing both the lattice geometry and the atomic positions. For biaxialstrained structures, the aand b-axis are fixed as per the strain value while the c-axis and the atomic positions are simultaneously relaxed. The plane wave cutoff energy is set to be 40 Ry, and a Γ-centered 10 × 10 × 2k-points in the full Brillouin zone are sampled. The structural and electronic convergence criteria were set to be 10 −4 Ry/a 0 and 10 −6 Ry, respectively.
With the DFT eigenvalues and eigenvectors, we then use Wannier90 package 57 and implement a disentanglement procedure to obtain the Wannier functions and the hopping parameters. During the disentanglement, 120 bands were considered, which accounted for contributions from 96 (32 × 3) O-2p orbitals and 24 (8 × 3) Ir-t 2g orbitals. Four iridium layers (i.e., 8 Ir atoms) are involved to construct 24 wannierized t 2g orbitals in one unit cell (6 Wannier functions per Ir layer).

Calculation of the photoemission spectra
The ARPES spectral function is obtained as imaginary part of the Green function G(k, ω) describing propagation of the photohole in the ground state, dressed in the low-energy magnetic excitations (magnons): Gðk; ωÞ ¼ Tr AF h jh k 1 ωÀHtÀJ þiδ h y k AF j i. Here, photohole h k is a vector in a full spin-orbital Hilbert space, obtained by projecting three t 2g orbitals of Ir atom onto the spin-orbit coupled basis 53,58 . The motion of the hole is governed by t À J Hamiltonian H tÀJ ¼ H t þ H J . The magnetic H J part includes NN, 2NN, and 3NN Heisenberg interactions J 1 , J 2 , J 3 and an anisotropic Δ term, arising from non-negligible Hund's coupling 54,59 : x iS The hopping H t part is derived by projecting multiorbital Hubbard model employing orbital-dependent hopping parameters t i onto spin-orbit coupled basis 53,58 . The motion of a charge excitation in the new spinorbit coupled basis is then expressed analytically in terms of these t i which are obtained directly from DFT 53,58 . We evaluate the Green function G(k, ω) using the SCBA 60 . SCBA is a diagrammatic approach that evaluates Green function of a quasiparticle dressed with bosons (here, photohole dressed with magnons) that form diagrams of rainbow type 60 . The spectral functions are calculated numerically for a 16 × 16 cluster.

DATA AVAILABILITY
All the data that support the findings of this study are available from the corresponding author (E.M.P.) upon reasonable request.