Electronic structure and optical properties of Sr$_2$IrO$_4$ under epitaxial strain
Churna Bhandari Zoran S. Popović S. Satpathy
   





Department of Physics & Astronomy, University of Missouri,
Columbia, MO 65211, United States of America










  
We study the modification of the electronic structure in the strong spin-orbit coupled   Sr_2IrO_4 by epitaxial strain  
using density functional methods. 
Structural optimization shows that strain changes the  internal structural parameters such as the Ir-O-Ir bond angle, which 
has an important effect on  the band structure.
An interesting prediction is the Γ -X crossover of the valence band maximum with strain, 
while the conduction minimum  at M remains unchanged.
This in turn suggests strong   strain dependence of the transport properties for the hole doped system,
but not when the system is electron doped.  
Taking the measured value of the Γ-X separation for the unstrained case, we predict the 
Γ -X crossover of the valence band maximum to occur for the tensile epitaxial strain  e_xx≈ 3%.
A minimal tight-binding model within the J_ eff = 1/2 subspace is developed to describe the main features of the band structure.
The optical absorption spectra under epitaxial strain are computed using density-functional theory, which explains the  observed anisotropy in the optical spectra with the polarization of the incident light.
We show that the optical transitions between the Ir (d) states, which are dipole forbidden, can be explained in terms of the admixture 
of Ir (p) orbitals with  the Ir (d) bands.















Electronic structure and optical properties of Sr_2IrO_4 under epitaxial strain
    Churna Bhandari, Zoran S. Popović and S. Satpathy
    December 30, 2023
===============================================================================










§ INTRODUCTION


The   5d oxides such as Sr_2IrO_4 (SIO) are of considerable current interest due to the presence of a large spin-orbit coupling (SOC)  which leads to many  novel features such as the  spin-orbit assisted J_ eff = 1/2 Mott insulator<cit.> and spin-orbital entangled electron states.

It has been suggested that the spin-orbital entanglement induced by the strong SOC in these structures could 
make these materials hosts for several unconventional features such as  the Kitaev model<cit.>, 
quantum spin Hall effect at room temperature,<cit.> or unconventional superconductivity<cit.>. 




Strain is an important parameter for probing the nature of the  electron states, and it can be induced by pressure or by epitaxial growth on  lattice-matched substrates. 
Tuning of the band gap is an important aspect for functional manipulation for potential device applications.  

There have been several studies of SIO under epitaxial strain condition, both from theory and experiments.
Samples of SIO have been grown epitaxially on a number of substrates such as SrTiO_3, LaAlO_3, GdScO_3, etc.<cit.>.
Resistivity and optical absorption measurements on these structures have shown that the  Mott-Hubbard gap is preserved under epitaxial strain, but its magnitude can be tuned by varying the strain. The changes in the electronic structure show up in the optical properties
as red or blue shift of the optical absorption under strain condition<cit.>.  

Several  theoretical works have also addressed the electronic structure of SIO under strain<cit.>. 
The optical properties  were studied by 
Zhang et al.<cit.> and Kim et al.<cit.>. 
These calculations were limited to the low energy range (0 - 2 eV) and, furthermore, the polarization dependence
 of the optical absorption
has not been studied theoretically,
even though experiments show a strong anisotropy for absorption with  E ∥ plane vs. E ∥ẑ .
Thus, a full understanding of the electronic and optical properties is still missing.

In this paper, we study the effect of the epitaxial strain on the electronic structure and optical absorption of SIO from density-functional theory. We find an interesting Γ-X crossover of the valence band top under strain, which we explain from a tight-binding model,
 and also find the anisotropy in the absorption spectrum for light polarized along the plane vs. normal to the plane, in agreement with the experiments. The dipole-forbidden d-d transitions are explained in terms of admixture of the Ir p orbitals with the d bands.



 

The organization of the paper is as follows. In Section <ref>,
we describe the crystal structure, the method of calculation, and the results of the density-functional structural optimization 
of the crystal structure under strain.
Density-functional results for the band structure and magnetic moments under strain are discussed in Section <ref>. 
Section <ref> discusses the optical absorption spectrum under strain, and  
the results are  summarized in Section <ref>.
  




§ CRYSTAL STRUCTURE AND METHOD OF CALCULATION
 

Bulk Crystal Structure –
The crystal structure of SIO with space group (142) I4_1/acd consists of eight formula units in the √(2) a_0 × c unit cell<cit.> as shown in Fig. <ref>, with the lattice constants a = b = 5.497 Å  and c = 25.798 Å. 
(An equivalent unit cell with the body-centered tetragonal lattice and four formula units in the basis may also be used.) 
Note that a = √(2) a_0.
 The structure shows a staggered rotation of the IrO_6 octahedra about the c-axis by the angle θ = 11.5^∘. The unit cell has four IrO_6 planes stacked along the c-axis, each plane consisting of two octahedra,
with staggered rotation angles and antiferromagnetic Ir moments. 
The electronic structure is essentially controlled by the Ir atoms with the 5d^5 configuration placed in the crystal field of the 
oxygen octahedra, with the Sr atoms playing a passive role of donating electrons to the system, described by the nominal chemical formula Sr_2^ 2-  Ir^4+ O_4^ 2-.
The structure can be thought of as IrO_2 layers separated from each other by two intervening SrO planes, making the structure quasi-two dimensional (2D). This means that the basic electronic structure can be modeled by a single plane of the Ir atoms in a minimal model, which we discuss later.











Density-Functional Methods –
Density functional theory (DFT)  with linearized full potential muffin-tin orbital method (FP-LMTO)<cit.> was used to solve the Kohn-Sham equations within the local spin-density approximation for the exchange and correlational
 functional<cit.>.
  The LMTO basis set consisted of spdf orbitals for Ir and Sr and spd orbitals for O, 
which were augmented outside the muffin-tin spheres by two Hankel functions of two different decay lengths.
 The semicore states Ir 5p and Sr 3p were  treated as valence electrons.  
 The band calculations were carried out for the optimized structures within the local spin density approximation including both the SOC and the Hubbard U terms  (LSDA + SO + U). Following earlier authors, we used U = 2.7 eV. We also employed Vienna ab initio simulation package (VASP) in projector augmented wave (PAW) formalism including Hubbard U=2 eV and spin-orbit interaction<cit.> for computing magnetic moments. 


 Structural optimization –

Atomic positions were relaxed in all our calculations within the local spin-density approximation and the force convergence was obtained with a tolerance of 10^-3 Ryd/Bohr. For the bulk, unstrained structure, we took the experimental unit cell and optimized
the positions of the atoms. An important structural parameter is the IrO_6 octahedral rotation angle θ, which we found to be 
13.5 ^∘ as compared to the measured angle of 11.5 ^∘. 


For the epitaxially strained structure, the in-plane lattice constant was changed according to the strain parameter e_xx≡ (a - a_ bulk)/ a_ bulk, while the out-of-plane lattice constant
was determined from the condition that the cell volume is preserved under strain, so that e_zz≡ (c -c_ bulk)/ c_ bulk = - 2 e_xx. Structural parameters were then optimized for each strain condition, with the fixed lattice constants. We varied the in-plane strain e_xx by ± 2 %, which is the same order of magnitude as the strains present in the experimental structures, e. g., e_xx≈ 2% for SIO grown on the GdScO_3  substrate
and -2 % for the NdGaO_3 substrate<cit.>. The computed structural parameters are shown in Table <ref> and
Fig. <ref>. The trend in the structural parameters shown here is similar to the results obtained from an earlier work
using the dynamical mean-field theory<cit.>.

For test purposes, we also performed a second set of calculations, where we took the in-plane lattice constant to be the 
same as the substrate as usual, but used the measured c/a ratio from the experiments<cit.>.
No significant changes in the bond angles and bond lengths were found.  



§ DENSITY FUNCTIONAL BAND STRUCTURE
  



 §.§ Electronic structure and Γ-X crossover
      


The electron bands for the optimized structures under three different strain conditions are shown in Fig. <ref>,
calculated within the local spin-density approximation with Coulomb interaction and spin-orbit coupling included (LSDA+SO+U). 
The Fermi surfaces for the doped SIO are shown in Fig. <ref>.

The results show systematic changes of the  band structure features in the gap region, some of which are shown in Fig. <ref>. The magnitude of the fundamental gap Δ_g increases with tensile strain, and so does the direct gap Δ_X at the X point, while the direct gap at the Γ point does not change very much. 
This trend in the gap values may be expected, since tensile strain would reduce the band widths of the LHB and the UHB, while the Coulomb U is relatively unchanged. The opposite happens for the compressive strain. The Γ point beats the trend because the valence top at
Γ has a considerable contribution from the oxygen states. 

The approximately linear variation of the gap with strain is clearly visible in the partial density-of-states shown in Fig. <ref>.
This trend of the gap variation with strain has already been noted in earlier works<cit.>.


 Strain induced Γ - X crossover –  An interesting feature of the band structure is the drastic change of the valence band maximum under strain, which would have a significant effect on the transport properties of the hole-doped structure. As can be seen from the band structure, Fig. <ref>, the conduction band minimum does not change under strain, always occurring at the M point, while it changes with strain. In the unstrained structure, DFT results show that the energy of the valence maximum at X and Γ are nearly the same. 
 From our calculations, the Γ - X separation is E_Γ X≡ E_Γ - E_X = 0.01 eV, which compares well with the earlier DFT calculations in the literature, where E_Γ X≈ 0.01 - 0.07 eV<cit.>.
Thus according to the DFT results, Γ is slightly below X, which is opposite to the ARPES measurement <cit.>, where the energy at X point is found to be above Γ by about 2 eV, i. e., E_Γ X≈ 0.2 eV <cit.>.
For the tensile strain, the valence maximum at Γ increases in energy, while that of X goes down; For the compressive strain, the reverse happens, leading to a Γ-X crossover with strain. 


The Γ - X crossover may be described within a tight-binding model involving the |J_ eff = 1/2 ⟩ orbitals on the square lattice of Ir atoms
and by including the dependence of the TB hopping integrals of the octahedral rotation angles which vary with strain.
This TB model is developed in <ref>.
The TB model shows the correct trend with the result:
  Δ E_Γ X (TB)  ≈ 1.65   eV×  e_xx (Eq. <ref>) as compared to 
  Δ E_Γ X ( DFT)  ≈ 4.0   eV× e_xx,  obtained from the DFT results shown in Fig. <ref>.
  The TB analysis indicates that both the angle and distance changes with strain are important for the description of the Γ-X crossover. 

  Based on the computed slope of E_Γ X from DFT and the experimental Γ-X separation of 0.2 eV, we would need 
an estimated epitaxial compressive strain of e_xx≥ 3% in order to switch the valence band maximum from X to Γ.
 The Γ-X crossover is clearly seen from the Fermi surface plots shown in Fig. <ref>, where we have shown the
 Fermi surfaces for 5% dopant concentration. 

For the case of electron doping, we find that the conduction band minimum always occurs at the M point in the 
Brillouin zone and the electron pocket is elliptical in shape. 


 







 



 





 §.§ Magnetic moment under strain



The calculated spin and orbital magnetic moments with and without strain are listed in Table <ref>. 
We note that the total energy obtained with spin moments constrained along the planar direction (x̂) has a lower energy than
the spin moments constrained along ẑ, consistent with the x̂ orientation of the magnetic moments.
Table <ref> shows the calculated magnetic moments using VASP with the LSDA+U+SO functional.

There has been much interest in the magnetic moments in SIO, in particular on the ratio  μ_l / μ_s <cit.>,  
since a ratio of two indicates the spin-orbital entanglement of the wave function and a deviation from this value is indicative of a mixture 
between J_ eff = 1/2 and 3/2 sectors, as suggested from Table <ref>.
For the unstrained material (e_xx = 0), the calculated total moment μ_l + μ_s  = 0.38 μ_B is in good agreement with values obtained from  the magnetic susceptibility measurements (0.5 μ_B) <cit.> as well as  from previous calculations (0.36 μ_B)<cit.>. 

As already mentioned, the d^5 configuration of Ir can be thought of as a single hole in the t_2g manifold, leading to the cubic-field  values,  μ_l = 2/3  μ_B and μ_s = 1/3  μ_B,  as seen from Table <ref>.
 When a tetragonal field is present as in the case of SIO, the | J_ eff , m ⟩ states get mixed among themselves, and the 
 magnetic moments can be substantially altered from the cubic-field value. 
 
 From our calculations, we find μ_l / μ_s ≈ 2.6 (Table <ref>) for the unstrained structure, in general agreement with earlier calculations<cit.> as well as with a recent measurement 
 using non-resonant magnetic x-ray diffraction which obtained the value μ_l / μ_s ≈ 2.5. 
 An earlier x-ray absorption measurement<cit.>  yielded the  ratio
μ_l / μ_s ≈ 1.1; the reason for the discrepancy between the two measurements is unclear. 

The magnetic moment ratio as well as its variation with strain can be approximately described by invoking a tetragonal crystal field for 
the single ion in the presence of the spin-orbit coupling. For more accurate description,  
a renormalized spin-orbit coupling has been invoked<cit.>.
The results for the single ion in the tetragonal field are given in <ref>. For the spin moment aligned along x̂, which is the case for SIO, we have 
μ_l / μ_s = 2 + 4 ξ / 3, where ξ = ε / λ is the ratio of the tetragonal field to the spin-orbit coupling strength
λ≈ 0.4 eV. For the unstrained case, ε≈ 0.14 eV <cit.>,
so that the ratio μ_l / μ_s = 2.47,
in reasonable agreement with the DFT result stated in Table <ref>.     
If we use the results of Ref. <cit.> which suggest that the magnitude of ε increases (decreases) by about 0.1 eV for
tensile (compressive) epitaxial strain of 2 %, then μ_l / μ_s = 2.8 for the tensile case and 2.1 for the compressive case, which more or less explains the DFT calculated trend shown in Table <ref>. 

 
 
 

 
















§ OPTICAL ABSORPTION


The change in the band structure with strain is reflected in the optical absorption spectrum.
The basic quantity to compute is the  imaginary part of the 
dimensionless
dielectric constant ε_2(ω), from which the real part ε_1 (ω) and the 
refractive index n (ω) are
computed using the Kramers-Krönig relation. 
In the dipole-approximation, the optical absorption coefficient α (ω)  is  given by


    α (ω) =    ω  c^-1  n(ω)^-1×ε_2(ω) ,



    ε_2(ω) =   4π^2e^2/m^2ω^2∑_v, c∫_BZd^3 k/( 2π)^3      
    | ê· M_cv ( k)|^2       δ(ħω-ϵ_c k+ϵ_v k),


where 
M_cv ( k)      = ⟨ψ_c k|   p |ψ_v  k⟩
is the momentum matrix element between the conduction and the valence states (defined as unoccupied and occupied states, respectively), 
and ê is the
light polarization vector. A closely related subsidiary function, useful to the discussion of the optical absorption, is the
joint density of states


    JDOS (ω) =    ∑_v, c∫_BZd^3 k/( 2π)^3δ(ħω-ϵ_c k+ϵ_v k).




It is important to note that even though the JDOS (ω) scales up linearly with the size of the unit cell chosen,
the optical absorption coefficient α (ω) and
ε_2(ω) are both independent of the size of the
unit cell used in the band calculation, as they must be. This is because a large number of transitions in the larger unit cell are simply forbidden.

This is easily seen by realizing that
 if a larger unit cell is used in the calculation, then a large number of vertical transitions in the folded Brillouin zone (larger unit cell)
is disallowed because the matrix element M_cv ( k) becomes zero from Bloch symmetry, 
even though the conduction and valence states both have nominally the same momentum k.
Alternatively, one can start with the allowed optical transitions in the Brillouin zone corresponding to the smallest unit cell and then
fold them into the smaller Brillouin zone corresponding to the larger unit cell and explicitly see that not all vertical transitions are allowed in the smaller Brillouin zone. The joint density of states JDOS (ω), in contrast, scales with the size of the unit cell, which can be seen by computing its integral over energy by taking advantage of the δ function in the definition.
Note  that the dipole approximation used in the expressions omits the local field and excitonic effects. 

Results are presented below for plane-polarized light with  polarization direction in the plane  
or normal to the plane. 














 §.§ Unstrained bulk 
 

The calculated absorption spectra for bulk SIO without strain are shown in Fig. (<ref>) for light polarization along the plane   (E ∥  plane), which is also compared to the measured data<cit.>.
The spectra show three distinct peaks which can be understood from the transitions between the 
J_ eff = 1/2 and J_ eff = 3/2 states as shown in Fig. (<ref>). 
The J_ eff = 3/2 states are split into two as seen from Table <ref> 
due to the tetragonal crystal field V_t = 2ε / 3.

In the experiment, only two distinct peaks are seen (these are labeled α and β), while the weaker γ peak is missing, possibly
due to instrumental broadening. 

We did not find any significant asymmetry in α(ω) for polarization directions within the plane. 
However, there is a large difference whether  E ∥  plane or E ∥  z, as discussed later . 




 §.§  Dipole selection rules and the d-d optical transitions
 


 


The optical absorption in the low energy range occurs due to transitions within the Ir t_2g manifold,
transitions that are however not 
dipole allowed due to the optical selection rules (δ L = ± 1 for plane-polarized light).
In this subsection, we discuss the optical transition matrix elements
and  conclude that they become dipole allowed due to the admixture of the Ir p states into the Ir d bands.
To get the dipole matrix elements, we need to estimate two things, as indicated from Eq. (<ref>) below: 
(i) The amount of the Ir p - d admixture η, which we estimate from perturbation theory as well as
from density-functional results and (ii) The optical matrix element for dipole-allowed transition between the Ir p and d orbitals
⟨  p |  ê· p |d ⟩, which we estimate from an effective hydrogenic model.


The d-d transitions are 

 non-zero in the crystal only because there is deviation from the spherical symmetry, the so called crystal field,
 which also produces splitting between the t_2g and the e_g states<cit.>. 
 The crystal field not only splits the Ir (5d) orbitals, but it also mixes  the Ir (5p) and (4f) orbitals into the 5d states, so that the d-d optical transitions have non-zero matrix element.
 Alternatively, such transitions may be equivalently described to be due to the transfer between Ir d atomic orbitals to the adjacent atoms because of the wave function overlap, e. g., within the linear combination of atomic orbitals (LCAO) model.


We now estimate the dipole matrix elements, which are non-zero due to the p-d and f-d mixing for the Ir atom. 
We first consider the p-d mixing, which as we shall see makes the larger contribution.
Symbolically, this is
given by


    M_dd   =   ⟨ψ_c  k|  ê· p |ψ_v  k⟩
       =   ⟨ d + η p |  ê· p |d + η p⟩≈      2 η×⟨  p |  ê· p |d ⟩.


The admixture η of the Ir p states into the Ir d bands, although small due to the large p-d energy difference, is nevertheless significant enough and is largely responsible for the optical transition in  SIO. 


The magnitude of η can be estimated from the partial density of states (PDOS) in the band calculations in the energy range of the Ir d bands. With the definition for the PDOS,
ρ_d (ε) = ∑_i | ⟨ψ_i | d ⟩ |^2 δ (ε - ε_i) and  
ρ_p (ε) = ∑_i | ⟨ψ_i | p ⟩ |^2 δ (ε - ε_i), 
and with the wave functions 
ψ≈ d + η p, one can estimate the average value of the admixture from the two PDOS in the energy range of interest. 
The PDOS, computed using the LAPW method,<cit.> are shown in Fig. <ref> and we have obtained the value of η from the ratio
of the integral of the two PDOS  in the energy window of the Ir t_2g bands, viz., from -2 to 1 eV. 
The result is: η^2 ≈ρ_p / ρ_d ≈ 0.19 × 10^-2 or η≈ 4%. We find a similar admixture of the 4f states, ρ_f/ρ_d ≈ 0.2 × 10^-2, indicating a nearly equal 4f admixture η_4f≈ 4% into the Ir (5d) bands.

Alternatively, the strength of the admixture η can be estimated from perturbation theory. From the second order  perturbation theory, we have


    η = 3 ×⟨ p | V_cf | d ⟩  / ( ε_d - ε_p ),


where the factor of three is from the degeneracy of the p states. 
If we further approximate  ⟨ p | V_cf | d ⟩∼⟨ d | V_cf | d ⟩≈Δ_cf, admittedly a crude approximation, and take the  t_2g -  e_g splitting Δ_cf∼ 2 eV, and use the result  ε_d - ε_p ∼ 43 eV for Ir from the standard Atomic tables<cit.>, we get the perturbation theory result of  η≈ 14 %, which is a factor of three too high as compared to the DFT result. However, given the crudeness of the approximation we made for the matrix elements in the perturbation theory, the order of magnitude agreement is reasonable. In the following, we shall use the DFT value for the p-d admixture  η≈ 4 %.


The second part is the estimation of the 
dipole matrix element M_pd≡⟨  p |  ê· p |d ⟩ appearing in Eq. (<ref>), which we do  
from an effective hydrogenic  model<cit.>.
In this model, the atomic wave functions are described as hydrogenic wave functions, but with an orbital dependent effective nuclear charge Z.  Thus


    Ψ_nlm   =    R_nl (r) Y_lm (Ω), 
    
    R_nl (ρ)    =     N_nl (Z) ρ^l e^-ρ/2 L_n+l^2l +1 (ρ) ,


where ρ = (2Z/na_0) r. Note that Z here is an effective atomic number, which takes into account the 
screening of the core electrons and depends on the principal quantum number n also. We take the value<cit.> Z = 18.7 appropriate for Ir
n=5 orbitals (5p and 5d).
 The integration can be performed analytically to yield the result M_pd for the plane polarized light:


    M_pd = ⟨Ψ_510 |  p_z |Ψ_520⟩ = - 0.25 i


in atomic units (viz., ħ = 1, Bohr radius a_0 =1, m_e = 1/2, Energy unit = 1 Ryd). Note that for atoms, the direction of polarization doesn't matter due to spherical symmetry. Plugging in the estimated magnitudes of M_pd and η  in Eq. (<ref>),
we find the d-d matrix element due to the admixture of Ir (5p) into the Ir (5d) states


    |M_dd^p|  = 2η× |M_pd| = 2 × 4%× 0.25 = 2 ×10^-2.



A similar calculation for the matrix element due to the admixture with the Ir (4f) orbitals yields a number, which is an order of magnitude
smaller than |M_dd^p|, viz.,

    |M_dd^f| = 2 η_f × |M_fd| = 2 × 4%× 0.064 = 5 ×10^-3,
 

where we estimated the 4f admixture η_f from Fig. (<ref>) and used the effective atomic number Z = 38.3 for the Ir 4f orbitals<cit.>.  
Thus the Ir (4f) orbitals contribute a much smaller amount 
as compared to the Ir (5p) orbitals. 


The estimate  given by  Eq. (<ref>), M_dd∼ 10^-2,  is very comparable to the average matrix element M_cv∼ 10^-2  obtained from the band calculations (Fig. <ref>), indicating that the admixture of the Ir (5p) orbitals is indeed responsible for the optical absorption. 


 












 §.§ Optical absorption under strain


Uniform strain –

The calculated optical absorption for uniform compression (uniform pressure) is compared to the same for bulk SIO without any compression in Fig. (<ref>). Apart from a shift in the peak positions to lower energies, caused by a reduction in band gap due to the
increase of the band width due to compression, we note that the overall absorption coefficient is increased for both light polarizations.
This can be explained from our above argument of Ir 5p - 5d admixture, because a larger crystal field upon compression mixes the Ir 5p more
into the Ir 5d bands, making the optical matrix element larger. 
Indeed, following the same logic, as the lattice constant is increased, there is less and less admixture of the Ir (5p) orbitals, and in the limit of infinite lattice constant, the admixture vanishes (perfect spherical symmetry), so that optical absorption would be zero as would be expected for the d → d  transition due to the dipole selection rules.

Epitaxial strain –

The change of the optical absorption with epitaxial strain is shown in Fig. (<ref>). 
With compressive epitaxial strain (e_xx < 0 and e_zz > 0), the distances in the plane are reduced, leading  to a larger optical absorption for E ∥ plane  as argued for the uniform compression case.
The same effect leads to an overall reduction of the optical absorption for E ∥ẑ. 
For the absorption with E ∥ plane,
all three peaks α, β, and γ are still there, 
but the peaks are red (blue) shifted with compressive (tensile) epitaxial strain, consistent
with the band structure changes under strain shown in Fig. (<ref>).


   


  

Optical anisotropy – 
The polarization dependence of the optical absorption has been reported in the literature<cit.>, 
where a significant anisotropy is found in the optical spectra. The absorption is significantly reduced for the  polarization
E ∥ẑ as compared to E ∥ plane in the low energy region 
corresponding to the transition within the Ir t_2g manifold, while for the higher energy region, the difference in the relative strength is not as drastic.
This is reproduced quite well from our calculations as seen from Fig. <ref>, where we also compare with the existing experiments. 
The calculated spectra is for the unstrained structure, while the experimental spectra is for a system with non-uniform strain,
the only case
for which experimental results are available for both polarizations.

The anisotropy is directly attributable to the differences in the matrix element | ê· M_cv ( k)| (Eq. <ref>),
since the variation of the refractive index n (ω) with energy or polarization is relatively weak, as indicated from Fig. <ref>.





§ SUMMARY
      

In summary, we studied the electronic properties and the optical absorption spectra of Sr_2IrO_4 under epitaxial strain condition using density-functional methods. Systematic structural changes with strain were found including the staggered rotation angle θ, which was found to have important effect on the electronic structure. An interesting result is the Γ - X crossover of the J_ eff = 1/2 valence band maximum with strain, allowing for engineering of the hole pocket in the hole-doped material, with potentially drastic changes in the transport properties. 
A minimal tight-binding Hamiltonian was developed for the J_ eff = 1/2 sector, which is capable of qualitatively describing the
important features of the band structure under strain, including the Γ- X crossover. 


We calculated the optical spectra under strain, interpreting the results in terms of small admixtures of the Ir (5p) and (4f) states with the Ir (5d) bands, without which the optical transitions would be dipole forbidden. The calculated spectra were compared to the experiments, where available, and the observed anisotropy in the optical absorption, which is very strongly anisotropic for the low energy region but 
less strongly anisotropic for the high energy region (see Fig. <ref>), was correctly  
explained from the calculated results.

Our work opens up the possibility of exploring the strain manipulation of the transport properties of epitaxially grown spin-orbit coupled Mott systems.




We thank Jamshid Moradi Kurdestany for valuable discussions 
and the U.S. Department of Energy, 
Office of Basic Energy Sciences, Division of Materials Sciences and Engineering (Grant  No. DE-FG02-00ER45818) for financial support. Computational resources were provided by the National Energy Research Scientific Computing Center, a User Facility also supported by the U.S. Department of Energy. 







§ TIGHT-BINDING MODEL AND THE Γ - X CROSSOVER 
 

In this Appendix, we construct a minimal tight-binding model to describe the band structure within the J_ eff = 1/2 sector on the square lattice
of Ir atoms, appropriate for SIO, in order to understand the Γ - X crossover of the valence band with epitaxial strain.




We are primarily interested in the J_ eff = 1/2 sector, viz., the lower and the upper Hubbard bands (LHB / UHB)  that form the valence and the 
conduction bands in the gap region. 
A tight-binding description can follow two different paths: One is to keep the t_2g orbitals (in total six orbitals in the basis per Ir atom including spin) in the Hamiltonian and then fit the the two Hubbard bands in the DFT band structure and the other is to keep only the two J_ eff = 1/2  orbitals per Ir atom (the minimal model). 
For the fitting with the t_2g orbitals, the NN TB parameters that best fit the two Hubbard bands are (in eVs): V_π = -0.20, V_δ = 0.04, U = 0.78, λ = 0.4. 
In our discussions below, we employ the second approach, where we use the minimal  TB model with just the two orbitals per Ir atom to describe the J_ eff  = 1/2 sector.



We consider the Hubbard model on a square lattice with anti-ferromagnetic order as appropriate for SIO, 
keeping the two spin-orbital entangled J_ eff  = 1/2 orbitals on each site, 
which we call e_1 and e_2, defined with respect to the 
local octahedral axes, with the staggered rotations as indicated in Fig. (<ref>), viz.,


    |e_1⟩   ≡    |1/2, -1/2⟩ =
    (|xy↑⟩ + |yz↓⟩ + i|xz↓⟩)/√(3)
    
    |e_2 ⟩   ≡    |1/2, 1/2⟩ =
    (|yz↑⟩ - i|xz↑⟩ - |xy↓⟩)/√(3).


The TB Hamiltonian is  given by


    H = ∑_⟨ ij⟩α t_ij c^†_iαc_jα+h.c. + U/2∑_iα n_iα n_iα̅,


where c^†_iα creates an electron at the site i (which may be in sublattice A or B) and in the orbital   e_α, t_ij is the hopping integral, which is non-zero only for hopping between the same type of orbitals e_1 or e_2, U is the on-site Coulomb interaction, and the summation ⟨ ij⟩
indicates sum over distinct pairs of bonds.  Note that in general the hopping integral t_ij are complex numbers as discussed below.
We will retain hoppings up to the fourth NN as indicated  in Fig. (<ref>). 
Furthermore, we find that at least three (and sometimes four) NN hoppings need to be kept  for an accurate description of the band structure in this minimal
model for the  J_ eff  = 1/2 sector.  


To obtain the hopping integrals between these orbitals, it is convenient to first obtain the integrals in the unrotated d basis from standard Tables<cit.>
and then rotate the basis.
Under rotation the angular momentum functions transform among one another, and in the present case, we have a 
 site-dependent 5 × 5 rotation matrix appropriate for L=2.
 
We denote the unrotated basis as |α⟩ (xy, yz, zx, x^2 - y^2, and 3z^2 -1, in that order), rotated basis by 
|α^'⟩ on sublattice A and 
|α^''⟩ on sublattice B, 
and the corresponding rotation matrices that
transform one basis into another  by 
R^' and R^'',  i. e., |α^'⟩ = R^' |α⟩ and
|α^''⟩ = R^'' |α⟩.

The hopping integrals in the rotated basis are then given by  

H_α^'β^''≡⟨α^' | H | β^''⟩ = 
⟨α | R^' T HR^'' |β⟩, or


    H =  R^' T  H R^''.

The rotation matrix  R (θ) for L =2 with rotation θ about the ẑ  axis is well known<cit.>

 
    R (θ) =       [  cos2θ      0      0  sin2θ      0;      0   cosθ   sinθ      0      0;      0  -sinθ   cosθ      0      0; -sin2θ      0      0  cos2θ      0;      0      0      0      0      1 ].
 

Due to the staggered rotations, the hopping integrals between orbitals on the same sublattice and opposite sublattices are given using Eq. (<ref>) as follows:
H_AA =  R (θ)^ T  H_AA R (θ), H_AB =  R (θ)^ T  H_AB R (-θ), H_BB =  R (-θ)^ T  H_BB R (-θ), and H_BA =  R (-θ)^ T  H_BA R (θ). 
Denoting the two rotated J_ eff = 1/2 states corresponding to Eq. (<ref>) by |e_1A⟩, |e_1B⟩, |e_2A⟩, and |e_2B⟩, one can then obtain the hopping integrals between these set of orbitals from the standard d-d hopping integrals from Harrison's Tables<cit.> and using the above rotation matrices. 
With the standard notations for the direction cosines nlm for the distance vector joining the first atom to the second, these matrix elements are readily obtained. The hopping amplitude between the opposite sublattice is 


    ⟨e_iA| H | e_j B⟩ = h_i δ_ij,  
           h_1 = t + i t^',   h_2 = h_1^*,


    3t       =[  3l^2m^2V_σ+(1-4l^2m^2)V_π  ]  cos^22θ
          
     - [ 3/4(l^2-m^2)^2V_σ+ (1-(l^2-m^2)^2 )  V_π ]   sin^2   2θ
            + V_πcos 2θ,


    t^' = -    V_π/3sin2θ.


Notice that the hopping matrix is diagonal with complex elements for non-zero θ, which however can be made real by making a 
gauge transformation. For hopping between the same sublattice, the hopping matrix is diagonal, but with real elements this time:


    ⟨e_iA| H | e_jA⟩ = t δ_ij,  
           3t =  V_π + [ V_π -(l^2-m^2)^2   
    (V_π- 3    V_σ / 4) ]   sin^2   2θ
           + [ V_π-l^2m^2(4V_π-3V_σ) ]  cos^2 2θ
           lm (l^2-m^2) (2V_π- 3 V_σ / 2)sin4θ.


We make a further simplification by taking V_σ = -(3/2 )V_π following Harrison's scaling<cit.>  in order to construct a minimal model.



The same expression Eq. (<ref>) is valid for hopping between B sublattice atoms ⟨e_iB| H | e_jB⟩, except that 
the sign of θ is changed due to the staggered rotations. 

Sometimes, in the literature, a simplistic angle dependence for the hopping h (θ) = h_0 cosθ is used. 
However, as the two equations above demonstrate, the angle dependence is more complicated and cannot be written down as cosθ even to
the lowest order in the angle.

In our TB model, we have retained up to four NN hoppings, as indicated in Fig. <ref>. 
 



Gauge transformation – As seen from Eq. (<ref>), the hopping integrals are in general complex,

and of course there is no problem working with the complex hoppings. 
However, it is convenient to make the NN integral real by a gauge transformation<cit.>. 

In the transformation, one simply adds a multiplicative phase factor to the orbital definitions, viz.,


    |e_1j⟩→          |e_1j⟩       e^i ε_j ϕ /2 
           |e_2j⟩→   |e_2j⟩         e^-i ε_j ϕ /2,


where  ε_j = ± 1 for j = A/ B sublattices. With the choice ϕ = tan^-1 (t^' /t), where t, t^' are 
given in  Eq. (<ref>)
for the 1NN hopping,
the new hopping integral becomes real, viz.,  t +i t^'→t̃, where t̃ = (t^2 + t^' 2)^1/2. 
This transformation leaves  all same-sublattice hoppings real, while the opposite-sublattice hoppings including t_4 and beyond continue to remain complex. One can choose to work with such complex hoppings, or else, simply ignore the imaginary parts, since hoppings between far-away neighbors are small anyway. We denote the gauge-transformed basis (<ref>) as
|e_1⟩ and |e_2⟩ in Eq. (<ref>)  and the 
 hopping amplitudes (real)   as t_i as indicated in Fig. (<ref>).

The antiferromagnetic lattice structure together with the fact that |e_1⟩ and |  e_2⟩ subspaces do not mix leads
to the 2 × 2 TB Hamiltonian in the momentum space


    H ( k) = ([ -Δ + h_11      h_12;    h_12^*  Δ + h_11 ]),

where h_11 = 4t_2  cos k_x cos k_y + 2 t_3 (cos 2 k_x + cos 2 k_y), 
h_12 = 2 t_1 (  cos k_x + cos k_y)  + 4 t_4 (cos 2 k_x  cos  k_y + cos  2 k_y  cos   k_x), a_0 = 1 here, and Δ = U/2 is the staggered field. Diagonalization readily yields the energies of the upper and lower J_ eff = 1/2 Hubbard bands to be


    ε_± ( k) = h_11  ( k) ±√(Δ^2 + h_12^2  ( k) ).




The TB parameters obtained from fitting Eq. (<ref>) to the density-functional bands for the optimized structure with no strain
are, in units of eV: 

U = 0.65, t_1 = -0.095, t_2 =0.015, t_3 = 0.035, and t_4 = 0.01. 
The TB fit to the DFT bands are shown in Fig. (<ref>) as 
dotted lines. An important feature of the band structure is the occurrence of the conduction minimum at the M point of the Brillouin zone, which the electrons would occupy in the electron doped system. It can be easily shown that for this to happen, the
TB parameters must satisfy the condition<cit.>

    t_3 > t_2/2 > -t_1^2/U,


which is clearly satisfied by our parameters given above.
Furthermore, the eccentricity of the elliptical energy contours around the  M point in the conduction band,
as seen in Fig. <ref>, is given by<cit.>
  e  ≡ (1 - r^2)^1/2, where the axis ratio 
r = [ (2 t_3 - t_2) /  (t_2 + 2 t_3 + 4 (t_1-2t_4)^2 / U)   ]^ 1/2≈ 0.6 for SIO.

We determine the TB parameters under strain condition with the following ansatz. 
Using the calculated angle θ≈ 13^∘ and the hopping expressions Eqs. (<ref>) and (<ref>), there is a one to one correspondence between t_i and V_π^i for the i-th neighbor.
 We obtain the values of V_π^i for different neighbors under no strain condition. In order to compute the hopping integrals under strain conditions, 
 we back substitute  V_π^i and the DFT optimized angles under strain into Eqs. (<ref>) and (<ref>).

In addition to the change of the angles, distances between atoms also change.  
Taking the variation of V_π to follow Harrison's
R^-5 scaling with distance R  and including the effect of θ from Eqs. (<ref>)  and (<ref>),  we can compute the TB hopping parameters t_i under the epitaxial strain conditions. 
We assume that strain enters the hopping parameters only via change of the hopping distances and the rotation angle θ.


With this simple TB model, we find that the important trends of the 
band structure under strain are described qualitatively correctly,
though not quantitatively with a factor of two to three discrepancy as compared to the DFT results. 


In particular, the model can qualitatively describe the Γ - X 
crossover of the valence band top under strain (hole pocket in the doped structure), predicted from the DFT calculations. 
From the TB energy expression Eq. (<ref>), we readily find the Γ - X  energy difference for the valence 
band

    E_Γ X = ε_- (Γ) - ε_- (X) = 8 t_2 + Δ - [ Δ^2 + 16 (t_1+2 t_4)^2 ]^1/2.


For our TB parameters, this is almost zero consistent with the DFT results. 
The strain dependence of  E_Γ X is obtained from the TB expression (<ref>), viz.,


    Δ E_Γ X = ∂   E_Γ X/∂ RΔ R +   
    ∂   E_Γ X/∂θΔθ≈ 1.65   eV×  e_xx,


where we have used the fact that the distance change between the atoms in the plane is simply
Δ R / R = e_xx,
and  also the DFT result for the change in angle with strain, viz., 
  Δθ≈ (-1.5   rad )  e_xx, obtained from Table <ref> assuming a linear strain dependence,
  since the strain is small. It turns out that the two terms in Eq. (<ref>) contribute nearly half each to the final result, so that both the angle and distance changes are important for the description of the Γ-X crossover. 




§ SINGLE ION IN A TETRAGONAL FIELD WITH SPIN-ORBIT COUPLING
 


It is instructive to examine the eigenstates of the single site Hamiltonian for the d orbitals in the presence of the SOC and a tetragonal crystal field, which this Appendix deals with.

We assume a large cubic crystal field splitting Δ_cf→∞, so that the e_g and the t_2g sectors don't mix. Within the  t_2g sector, the Hamiltonian is given by


    H_SOC=
         |xz↑⟩    |yz↓⟩    |xy↓⟩    |xz↓⟩    |yz↑⟩    |xy↑⟩⟨ xz↑|     0     -iλ/2     iλ/2     0     0     0⟨ yz↓|     iλ/2     0       -λ/2     0     0     0 ⟨ xy↓|     -iλ/2      -λ/2       -ϵ    0     0     0⟨ xz↓|    0     0     0     0     iλ/2     iλ/2⟨ yz↑|     0     0     0    -iλ/2     0    λ/2⟨ xy↑|   0     0     0    -iλ/2    λ/2    -ϵ.




 


By diagonalizing the Hamiltonian and keeping the tetragonal field to the lowest order, we obtain the eigenvalues and eigenfuctions, which are summarized in Table <ref> including some relevant expectation values.
 Note that the wave functions | J_ eff , m ⟩  listed in Table <ref>   are not necessarily eigenstates of J^2 and J_z, because they were not obtained by diagonalizing the  full 10 × 10 λL⃗·S⃗ matrix, but rather just the 6 × 6 matrix in the t_2g sector. 
 

The basic electronic structure of the iridates is determined by the Ir^4+ ions with the 5d^5 configuration
and a large crystal field due to the oxygen octahedra. The  crystal field combined with the SOC results in the spin-orbital entangled states as summarized in Table <ref>.
The  cubic  crystal field splits the 5d states into  e_g plus t_2g states. With the SOC included,  
the six-fold degenerate t_2g states (including spin) split into a two-fold 
J_ eff = 1/2 (higher energy)
and a four-fold   J_ eff = 3/2  (lower energy)
 state.  With the 5d^5 configuration, four electrons fill the J_ eff = 3/2 state, while the remaining electron
 occupies the lower Hubbard band separated from the upper Hubbard band by the Coulomb interaction U thus resulting in a SOC-induced Mott insulator. 
  The  tetragonal  field  has been  modeled in Table <ref> by  
asymmetric on-site energies for the Ir d orbitals, viz., 
 ε =  ε (xz/yz) - ε (xy)  and  δ = ε (3 z^2 - 1) - ε (x^2-y^2).
 The d^5 configuration is therefore equivalent to a single hole t_2g^1 with the |1/2, 1/2 ⟩ configuration.
 

 In Table <ref>, we have also listed the  angular momentum expectation values,
 which will be useful for the interpretation of the computed magnetic moments 
 discussed in Section <ref>. 
 The Table lists
   ⟨ L_z ⟩ and ⟨ 2S_z ⟩ for the spin state along ẑ;
    It is straightforward to show that  ⟨ L_x  ⟩ = ⟨ L_y ⟩ = ⟨ S_x  ⟩ = ⟨ S_y ⟩ = 0 in this case.  
Also listed are the expectation values 
 ⟨ L_x ⟩ and ⟨ 2S_x ⟩ 
  for the spin state along x̂ direction (which happens to be the case for the ground state of SIO); 
Again, ⟨ L_y  ⟩ = ⟨ L_z ⟩ = ⟨ S_y  ⟩ = ⟨ S_z ⟩ = 0 for this case. 
One point to note here is that the tetragonal field ξ  mixes up the J_ eff = 1/2 and 3/2 wave functions. 

 For spin magnetization along ẑ, the expectation values of the magnetic moments are 


    μ_l^z =  μ_B ⟨ψ_1/2   |L_z |ψ_1/2⟩ =  2/3+8ξ/27,


    μ_s^z= μ_B  ⟨ψ_1/2   |2S_z |ψ_1/2⟩=1/3+16ξ/27,


where |ψ_1/2⟩ = | 1/2, 1/2 ⟩ , the Bohr magneton μ_B = 1, and other components of the magnetic moments are zero,
leading to the ratio of the magnetic moments to be

    μ_l^z  /   μ_s^z   =  2-8ξ /3 .



For Sr_2IrO_4, the magnetization is along  x̂. 
With the spin wave function along x̂, | ψ̅_1/2⟩ = 2^-1/2  (| 1/2, 1/2 ⟩ + | 1/2, -1/2 ⟩ ), we then have


    μ_l^x  =  μ_B ⟨ψ̅_1/2   |L_x| ψ̅_1/2⟩= 2/3-4ξ/27,



    μ_s^x      = μ_B  ⟨ψ̅_1/2   |2S_x| ψ̅_1/2⟩ = 1/3-8ξ/27,

and the ratio

    μ_l^x /   μ_s^x   =  2+ 4ξ /3 .




  
  
  


§ REFERENCES

23

KimPRL08Kim B J, Jin H, Moon S J, Kim J -Y, Park B -G, Leem C S, Yu J, Noh T W, Kim C, Oh S -J, Park J -H, Durairaj V, Cao G and Rotenberg E 2008 Phys. Rev. Lett. 101 076402





JackeliPRL09Jackeli G and Khaliullin G 2009
Phys. Rev. Lett. 102 017205

Balentsnat10Balents L 2010 Nature 464 199 

ShitadePRL09Shitade A, Katsura H, Kuneś H J, Qi X L, Zhang S -C and Nagaosa N 2009 
Phys. Rev. Lett. 102 256403

WangPRL011Wang F and Senthil T 2011
Phys. Rev. Lett. 106 136402

WantanbePRL13Watanabe H, Shirakawa T and Yunoki S 2013
Phys. Rev. Lett. 110 027002


RyanPRB13Serrao C R, Liu J, Heron J T, Singh-Bhalla G, Yadav A, Suresha S J, Paull R J, Yi D, Chu J -H, Trassin M, Vishwanath A, Arenholz E, Frontera C, Zelezny J, Jungwirth T, Marti X and Ramesh R,
2013 Phys. Rev. B 87 085121

NicholsAPL13Nichols J, Terzic J J, Bittle E G, Korneta O B, De Long L E, Brill J W, Cao G and  Seo S S A, 
2013 Appl. Phys. Lett. 102 141908


LupascuLupascu A et al. 2014 Phys. Rev. Letts. 112 147201

Nichols2-APL13 Nichols J, Korneta O B, Terzic J, De Long L E, Brill J W, Cao G and SeO S S A 2013 
Appl. Phys. Lett. 103 131910



ZhangPRL13 Zhang H, Haule K and Vanderbilt D 2013
Phys. Rev. Lett. 111 246402


LadoPRB015Lado J L and Pardo V 2015
Phys. Rev. B 92 155151 



Kimsr016Kim B H, Kim K and Min B I 2016 Sci. Rep. 6 27095


CrawfordCrawford M K, Subramanian M A, Harlow R L, Fernandez-Baca J A, Wang Z R and Johnston D C
1994 Phys. Rev. B 49 9198

 
 
 MethfesselMethfessel M, Schilfgaarde M van and Casali R A, 
A Full-Potential LMTO Method Based on Smooth Hankel Functions, Electronic Structure and Physical Properties of Solids. The Use of the LMTO Method, Lecture Notes in Physics 535, 114 (2000)

 
Kotani10 Kotani  T and Schilfgaarde M van 2010 Phys. Rev. B 81 125117

lmsuite See, for details of the method: https://www.questaal.org


 vonBarthHedin Barth U von and Hedin L 1972 
 J. Phys. C: Solid State Phys. 5 1629
 
Kohn-Sham Kohn W and Sham L J, 
1965 Phys. Rev. 140  A1133 

KressePRB99 Kresse G and Joubert D 1999 Phys. Rev. B 59 1758

 
 
 
PeitaoPRM18Peitao L, Bongjae K, Xing-Qiu C, Sarma D D, Georg K and Cesare F 2018 Phys. Rev. Mat. 2 075003

 
YueNatcom16 Cao Y et al. 2016 Nat. Communs. 7 11367


Haskel Haskel  D et al. 2012 Phys. Rev. Letts. 109  027204 

Fujiyama Fujuyama S et al. 2014
 Phys. Rev. Letts. 112 016405
 
 Franchini Kim B, Liu P and Franchini C 2017  Phys. Rev. B 95, 115111
 
CaoCao G, Bolivar J, McCall S, Crow J E and Guertin R P 1998 Phys. Rev. B 57 R11039


KimNAT14 Kim J, Daghofer M, Said A H, Gog T,  van den Brink J, Khaliullin G and Kim B J 2014 Nat. Commun. 5 4453


BhandariJPCSBhandari C and Satpathy S 2017 J. Phys. Chem. Solids  https://doi.org/10.1016/j.jpcs.2017.10.047

LAPWSchwarz K, Blaha P and Madsen G K H 2002 Comp. Phys. Commun. 147 71

HermanHerman F and Skillman S, Atomic Structure Calculations (Prentice-Hall, New Jersey, 1963)


ZZClementi E, Raimondi D L and Reinhardt W P 1967 J. Chem. Phys. 47 1300


PaulingPauling L and Wilson E B, Jr.
Introduction to Quantum Mechanics with Applications to Chemistry
(Dover Publications, New York, 1985)

Harrison Harrison W A, 
Electronic Structure and the Properties of Solids: The Physics of the
Chemical Bond (Dover Publications, New York, 1989)

Avinash-osmateMohapatra S, Bhandari C, Satpathy S and A. Singh A  2018
Phys. Rev. B, 97 155154  https://doi.org/10.1103/PhysRevB.97.155154



Avinash-2017Mohapatra S, Brink J van den and Singh A 
2017 Phys. Rev. B 95 094435



Sayantika Bhowal S, Moradi Kurdestany J M and Satpathy S “Stability of the Antiferromagnetic State in the Electron Doped Iridates" (unpublished)

Bhandari18Bhandari C and Satpathy S 2018 Phys. Rev. B 98, 041303(R) https://doi.org/10.1103/PhysRevB.98.041303






