# Resonant inelastic x-ray scattering of the $J_{eff}$ = 1/2 Mott insulator
Sr2IrO4 from the density-functional theory

V.N. Antonov  G. V. Kurdyumov Institute for Metal Physics of the N.A.S. of
Ukraine, 36 Academician Vernadsky Boulevard, UA-03142 Kyiv, Ukraine  Max-
Planck-Institute for Solid State Research, Heisenbergstrasse 1, 70569
Stuttgart, Germany  D.A. Kukusta  G. V. Kurdyumov Institute for Metal Physics
of the N.A.S. of Ukraine, 36 Academician Vernadsky Boulevard, UA-03142 Kyiv,
Ukraine  L.V. Bekenov  G. V. Kurdyumov Institute for Metal Physics of the
N.A.S. of Ukraine, 36 Academician Vernadsky Boulevard, UA-03142 Kyiv, Ukraine

###### Abstract

We have investigated the electronic structure of Sr2IrO4 within the density-
functional theory (DFT) using the generalized gradient approximation while
taking into account strong Coulomb correlations (GGA+$U$) in the framework of
the fully relativistic spin-polarized Dirac linear muffin-tin orbital band-
structure method. We have investigated the x-ray absorption spectra (XAS),
x-ray magnetic circular dichroism (XMCD), and resonant inelastic x-ray
scattering (RIXS) spectra at the Ir $L_{3}$ and O $K$ edges. The calculated
results are in good agreement with the experimental data. The RIXS spectrum of
Sr2IrO4 at the Ir $L_{3}$ edge in addition to the elastic scattering peak at 0
eV possesses a sharp feature below 1.5 eV corresponding to transitions within
the Ir $t_{2g}$ levels. The excitation located from 2 eV to 5 eV is due to
$t_{2g}$ $\rightarrow$ $e_{g}$ transitions. The third wide structure situated
at 5$-$12 eV appears due to transitions between the Ir 5$d_{\rm{O}}$ states
derived from the ”tails” of oxygen 2$p$ states and $e_{g}$ and $t_{2g}$
states. The RIXS spectrum of Sr2IrO4 at the O $K$ edge consists of three major
inelastic excitations at 0.7 eV, 3.5 eV, and around 6.2 eV. We have found that
the first low energy feature is due to interband transitions between occupied
and empty O${}_{t_{2g}}$ transitions, which appear due to the strong
hybridization between oxygen 2$p$ and Ir $t_{2g}$ states in the close vicinity
of the Fermi level. The next two peaks at around 3.5 and 6.2 eV reflect the
iterband transitions from the occupied O 2$p$ states and the empty oxygen
states which arise from the hybridization with Ir $t_{2g}$ and $e_{g}$ states,
respectively. We have found that the theory reproduces well the shape and
energy position of the low energy feature, but the second and the third peaks
are shifted towards smaller energy in comparison with the experimental
measurements. To reproduce the correct energy position of the oxygen 2$p$ band
we have used a self-interaction-like correction procedure. We have found that
the dependence of the RIXS spectrum at the oxygen $K$ edge on the incident
photon energy and the momentum transfer vector Q is much stronger in
comparison with the correspondent dependence at the Ir $L_{3}$ edge.

###### pacs:

75.50.Cc, 71.20.Lp, 71.15.Rf

##  I Introduction

In 5$d$ transition metal compounds the energy scale of the spin-orbit coupling
(SOC) is comparable to the on-site Coulomb interaction and the crystal-field
energy. Due to the strong competition between these interactions fascinating
electronic states can arise. The SOC in such systems splits the $t_{2g}$
orbitals into a quartet ($j_{eff}$ = 3/2) and a doublet ($j_{eff}$ = 1/2)
Jackeli and Khaliullin (2009); Chen et al. (2010); Witczak-Krempa et al.
(2014). In 5$d^{5}$ (Ir4+) iridium oxides, the quartet $j_{eff}$ = 3/2 is
fully occupied, and the relatively narrow $j_{eff}$ = 1/2 doublet occupied by
one electron can be splitted by a moderate Hubbard $U_{eff}$ with opening a
small band gap called the relativistic Mott gap Kim et al. (2008); Martins et
al. (2011); Antonov et al. (2018). Iridates have been at the center of an
intensive search in recent years for novel phenomena, such as topological
insulators Qi and Zhang (2010); Ando (2013); Wehling et al. (2014); Bansi et
al. (2016), Mott insulators Kim et al. (2008, 2009); Jackeli and Khaliullin
(2009); Watanabe et al. (2010); Martins et al. (2011), Weyl semimetals
Witczak-Krempa and Kim (2012); Go et al. (2012); Sushkov et al. (2015), and
quantum spin liquids Jackeli and Khaliullin (2009); Kimchi et al. (2014).

Among iridium compounds, Sr2IrO4, a single-layer member of the Ruddlesden-
Popper series iridates, is of special interest. It has a quasi-two-dimensional
(2D) square-lattice perovskite structure and was the first discovered spin-
orbit $j_{eff}$ = 1/2 Mott insulator Kim et al. (2008). Besides, it has the
structural and physical similarity with La2CuO4, a parent compound to
high-$T_{c}$ cuprates, such as the presence of a pseudogap Kim et al. (2015);
Yan et al. (2015); Battisti et al. (2017), similar Fermi surfaces and Fermi
arcs (in electron- and hole-doped compounds) Damascelli et al. (2003); Kim et
al. (2014), $d$-wave symmetry Kim et al. (2016); Terashima et al. (2017),
electron-boson coupling Hu et al. (2019), and similarities in the magnetic
ordering and magnetic excitations Kim et al. (2012a); Clancy et al. (2019).
However, the superconductivity in Sr2IrO4 has not been found yet Bertinshaw et
al. (2019).

In this work we focus our attention on the RIXS properties in Sr2IrO4. Since
the first publication by Kao et al. on NiO Kao et al. (1996), the resonant
inelastic X-ray scattering method has shown remarkable progress in condensed
matter physics research as a spectroscopic technique to record the momentum
and energy dependence of inelastically scattered photons in complex materials
Ament et al. (2011). RIXS nowadays has rapidly become the forefront of
experimental photon science. It provides a direct probe of spin and orbital
states and dynamics. RIXS has a number of unique features in comparison with
other spectroscopic techniques. It covers a large scattering phase space and
requires only small sample volumes. It also is bulk sensitive, polarization
dependent, as well as element and orbital specific Ament et al. (2011).

Depending on the x-ray resonant energy, RIXS can be divided into two classes:
soft x-ray and hard x-ray Ament et al. (2011). For high atomic number
transition metal elements, such as 5$d$ transition metal compounds, the
$L$-edge resonant energies are in the hard x-ray region. For such spectra high
quality single crystals are needed as the key optical elements. The RIXS
resolution crucially depends on the availability of a high quality single
crystal with a Bragg diffraction peak close to back-scattering at the energy
of the $L$-edge of the targeted element. This requirement severely limits the
application of RIXS, and thus by far the majority of hard x-ray RIXS studies
have been focused on 5$d$ iridates Liu et al. (2012); Hozoi et al. (2014);
Clancy et al. (2016); Nag et al. (2018); Takayama et al. (2019); Aczel et al.
(2022) and osmates Calder et al. (2016); Taylor et al. (2017); Calder et al.
(2017). In a soft RIXS setup, the x-ray energy range is usually below about 2
keV Liu et al. (2015). The $L$-edges of the 3$d$ transition metal elements all
fall below this energy scale. The energy resolution in the soft x-ray region
is relatively high. For example, the combined energy resolution was 150 meV at
the Ni $L_{3}$ edge ($\sim$850 eV) in Ta2NiSe5 Monney et al. (2020). The best
energy resolution currently achieved for RIXS at the oxygen $K$-edge
($\sim$530 eV) of the common ligand atoms is about 45-50 meV Lee et al.
(2013); Lu et al. (2018), which is much better than the majority of 5$d$
elements probed using $L$-edge RIXS in the hard x-ray region up to date. We
should mention, however, that in recent years the experimentalists have
achieved a remarkable progress in increasing the resolution for hard RIXS
spectra. For example, Kim et al. Kim et al. (2023) have been able to obtain
the total resolution of 34.2 meV at the Ir $L_{3}$ edge in Sr2IrO4. Such
resolution permits direct measurements of single-magnon excitations as well as
other many-body excitations in strongly correlated systems.

In the x-ray absorption (XA), x-ray magnetic circular dichroism (XMCD), and
RIXS processes at the O $K$ edge, the 1$s$ core level is involved. The
exchange splitting of the 1$s$-core state is extremely small and SOC is absent
in the O 1$s$ orbitals, therefore, only the exchange and spin-orbit splitting
of the 2$p$ states is responsible for the observed spectra at the oxygen $K$
edge. On the other hand, the oxygen valence 2$p$-states of the surrounding
ligand atoms are sensitive to the electronic states at neighboring sites
because of their delocalized nature. They strongly hybridize with the 5$d$
orbitals. Due to such hybridization combined with high SOC at the 5$d$ ion,
information on the elementary excitations can be extracted using an indirect
RIXS process at the O $K$ edge Lu et al. (2018). Although O $K$ RIXS has a
much smaller penetration depth ($\sim$100 nm) than 5$d$ $L$ RIXS, a comparison
between O $K$ and Ir $L_{3}$ spectra measured on Sr2IrO4 suggests that they
have a comparable counting efficiency Lu et al. (2018). The lower penetration
depth of soft x-rays has its own advantages providing high sensitivity to
ultrathin samples such as films. Soft x-ray RIXS at the O $K$ edge is a
promising method for studying the electronic and magnetic excitations in 5$d$
compounds. There are several experimental investigations of the RIXS spectra
at the oxygen $K$ edge Liu et al. (2015); Lu et al. (2018); Paris et al.
(2020); Kim et al. (2023) in Sr2IrO4. The Ir $L_{3}$ RIXS spectra in this
oxide are investigated in Refs. Ishii et al. (2011); Kim et al. (2012a);
Lupascu et al. (2014); Paris et al. (2020); Bertinshaw et al. (2020); Kim et
al. (2023); Clancy et al. (2023).

We carry out here a detailed study of the electronic structure, the XA, XMCD,
and RIXS spectra of Sr2IrO4 in terms of the density functional theory. Our
study sheds light on the important role of band structure effects and
transition metal 5$d$ $-$ oxygen 2$p$ hybridization in the spectral properties
in 5$d$ oxides. The energy band structure, the XA, XMCD, and RIXS spectra of
Sr2IrO4 are investigated in the ab initio approach using the fully
relativistic spin-polarized Dirac linear muffin-tin orbital band-structure
method. We use both the generalized gradient approximation (GGA) and the
GGA+$U$ approach to assess the sensitivity of the RIXS results to different
treatment of the correlated electrons.

The paper is organized as follows. The crystal structure of Sr2IrO4 and
computational details are presented in Sec. II. Sec. III presents the
electronic and magnetic structures of Sr2IrO4. In Sec. IV the theoretical
investigations of the XA, XMCD, and RIXS spectra of Sr2IrO4 at the Ir
$L_{2,3}$ edges are presented, the theoretical results are compared with the
experimental measurements. In Sec. V we present the theoretical investigations
of the XA and RIXS spectra at the O $K$ edge. Finally, the results are
summarized in Sec. VI.

##  II Computational details

### X-ray magnetic circular dichroism.

Magneto-optical (MO) effects refer to various changes in the polarization
state of light upon interaction with materials possessing a net magnetic
moment, including rotation of the plane of linearly polarized light (Faraday,
Kerr rotation), and the complementary differential absorption of left and
right circularly polarized light (circular dichroism). In the near visible
spectral range these effects result from excitation of electrons in the
conduction band. Near x-ray absorption edges, or resonances, magneto-optical
effects can be enhanced by transitions from well-defined atomic core levels to
transition symmetry selected valence states.

Within the one-particle approximation, the absorption coefficient
$\mu^{\lambda}_{j}(\omega)$ for incident x-ray polarization $\lambda$ and
photon energy $\hbar\omega$ can be determined as the probability of electronic
transitions from initial core states with the total angular momentum $j$ to
final unoccupied Bloch states

$\displaystyle\mu_{j}^{\lambda}(\omega)$ $\displaystyle=$
$\displaystyle\sum_{m_{j}}\sum_{n\bf k}|\langle\Psi_{n\bf
k}|\Pi_{\lambda}|\Psi_{jm_{j}}\rangle|^{2}\delta(E_{n\bf
k}-E_{jm_{j}}-\hbar\omega)$ (1) $\displaystyle\times\theta(E_{n\bf
k}-E_{F})\,,$

where $\Psi_{jm_{j}}$ and $E_{jm_{j}}$ are the wave function and the energy of
a core state with the projection of the total angular momentum $m_{j}$;
$\Psi_{n\bf k}$ and $E_{n\bf k}$ are the wave function and the energy of a
valence state in the $n$-th band with the wave vector k; $E_{F}$ is the Fermi
energy.

$\Pi_{\lambda}$ is the electron-photon interaction operator in the dipole
approximation

$\Pi_{\lambda}=-e\mbox{\boldmath$\alpha$}\bf{a_{\lambda}},$ (2)

where $\bm{\alpha}$ are the Dirac matrices and $\bf{a_{\lambda}}$ is the
$\lambda$ polarization unit vector of the photon vector potential, with
$a_{\pm}=1/\sqrt{2}(1,\pm i,0)$, $a_{\parallel}=(0,0,1)$. Here, $+$ and $-$
denote, respectively, left and right circular photon polarizations with
respect to the magnetization direction in the solid. Then, x-ray magnetic
circular and linear dichroisms are given by $\mu_{+}-\mu_{-}$ and
$\mu_{\parallel}-(\mu_{+}+\mu_{-})/2$, respectively. More detailed expressions
of the matrix elements in the electric dipole approximation may be found in
Refs. Guo et al. (1994); Antonov et al. (2004); Arola et al. (2004). The
matrix elements due to magnetic dipole and electric quadrupole corrections are
presented in Ref. Arola et al. (2004).

### Resonant inelastic x-ray scattering.

In the direct RIXS process Ament et al. (2011) the incoming photon with energy
$\hbar\omega_{\mathbf{k}}$, momentum $\hbar\mathbf{k}$, and polarization
$\bm{\epsilon}$ excites the solid from the ground state $|{\rm g}\rangle$ with
energy $E_{\rm g}$ to the intermediate state $|{\rm I}\rangle$ with energy
$E_{\rm I}$. During relaxation an outgoing photon with energy
$\hbar\omega_{\mathbf{k}^{\prime}}$, momentum $\hbar\mathbf{k}^{\prime}$ and
polarization $\bm{\epsilon}^{\prime}$ is emitted, and the solid is in state
$|f\rangle$ with energy $E_{\rm f}$. As a result, an excitation with energy
$\hbar\omega=\hbar\omega_{\mathbf{k}}-\hbar\omega_{\mathbf{k}^{\prime}}$ and
momentum $\hbar\mathbf{q}$ = $\hbar\mathbf{k}-\hbar\mathbf{k}^{\prime}$ is
created. Our implementation of the code for the calculation of the RIXS
intensity uses Dirac four-component basis functions Nemoshkalenko et al.
(1983) in the perturbative approach Arola et al. (1997a, b). RIXS is the
second-order process, and its intensity is given by

$\displaystyle
I(\omega,\mathbf{k},\mathbf{k}^{\prime},\bm{\epsilon},\bm{\epsilon}^{\prime})$
$\displaystyle\propto$ $\displaystyle\sum_{\rm f}\left|\sum_{\rm
I}{\langle{\rm
f}|\hat{H}^{\prime}_{\mathbf{k}^{\prime}\bm{\epsilon}^{\prime}}|{\rm
I}\rangle\langle{\rm I}|\hat{H}^{\prime}_{\mathbf{k}\bm{\epsilon}}|{\rm
g}\rangle\over E_{\rm g}-E_{\rm I}}\right|^{2}$ (3)
$\displaystyle\times\delta(E_{\rm f}-E_{\rm g}-\hbar\omega),$

where the RIXS perturbation operator in the dipole approximation is given by
the lattice sum
$\hat{H}^{\prime}_{\mathbf{k}\bm{\epsilon}}=\sum_{\mathbf{R}}\hat{\bm{\alpha}}\bm{\epsilon}\exp(-{\rm
i}\mathbf{k}\mathbf{R})$, where $\bm{\alpha}$ are the Dirac matrices. The sum
over the intermediate states $|{\rm I}\rangle$ includes the contributions from
different spin-split core states at the given absorption edge. The matrix
elements of the RIXS process in the frame of fully relativistic Dirac LMTO
method were presented in Ref. Antonov et al. (2022a).

### Crystal structure.

The powder-neutron-diffraction measurements show that Sr2IrO4 possesses the
tetragonal $I4_{1}/acd$ perovskite structure (group number 142) [Fig. 1(a)]
Shimura et al. (1995). The IrO6 octahedra in Sr2IrO4 are rigidly aligned, just
as the CuO6 octahedra in cuprates, rotated by $\sim$11∘ about the $c$ axis in
the $a-b$ plane [Fig. 1(b)], and have a local distortion of 4.5% axial
elongation.

Figure 1: (Color online) (a) The schematic representation of the body centered
tetragonal $I4_{1}/acd$ (group number 142) Sr2IrO4 crystal structure Shimura
et al. (1995); (b) the positions of ions in the IrO2 plane perpendicular to
the $c$ axis.

Atomic positions of Sr2IrO4 at 10 K (the lattice constants $a$ = 5.48164 Å,
$c$ = 25.80019 Å) for Sr, Ir, O1, and O2 are (0, 1/4, $z_{\rm{Sr}}$), (0, 1/4,
3/8), (0, 1/4, $z_{\rm{O}}$), and ($x$, $x$ \+ 1/4, $z_{\rm{O}}$),
respectively, with $x$ = 0.1996, $z_{\rm{Sr}}$ = 0.5506, and $z_{\rm{O}}$ =
0.4548 Shimura et al. (1995). The oxygen atoms surrounding the Ir sites
provide an octahedral environment. The Ir$-$O1 and Ir$-$O2 interatomic
distances are equal to 2.05886 Å and 1.97704 Å, respectively. Around each Ir
atom there are eight Sr atoms with the Ir-Sr distance $d_{\rm{Ir-Sr}}$ =
3.34615 Å. The Ir-Ir distance $d_{\rm{Ir-Sr}}$ = 3.87610 Å.

Note that in our electronic structure calculations, we rely on experimentally
measured internal parameters $x$, $z_{\rm{Sr}}$, $z_{\rm{O}}$ and lattice
constants because they are well established for this material and are probably
still more accurate than can be obtained from DFT.

### Calculation details

The details of the computational method are described in our previous papers
Antonov et al. (2006, 2007, 2010); Antonov et al. (2022a) and here we only
mention several aspects. The band structure calculations were performed using
the fully relativistic linear muffin-tin orbital (LMTO) method Andersen
(1975); Antonov et al. (2004). This implementation of the LMTO method uses
four-component basis functions constructed by solving the Dirac equation
inside an atomic sphere Nemoshkalenko et al. (1983). The exchange-correlation
functional of the generalized gradient approximation (GGA)-type was used in
the version of Perdew, Burke and Ernzerhof Perdew et al. (1996). The Brillouin
zone integration was performed using the improved tetrahedron method Blöchl et
al. (1994). The basis consisted of Ir and Sr $s$, $p$, $d$, and $f$; and O
$s$, $p$, and $d$ LMTO’s.

To take into account the electron-electron correlation effects we used in this
work the ”relativistic” generalization of the rotationally invariant version
of the LSDA+$U$ method Yaresko et al. (2003) which takes into account that in
the presence of spin-orbit coupling the occupation matrix of localized
electrons becomes non-diagonal in spin indexes. Hubbard $U$ was considered as
an external parameter and varied from 0.65 eV to 3.65 eV. We used in our
calculations the value of exchange Hund coupling $J_{H}$=0.65 eV obtained from
constrained LSDA calculations Dederichs et al. (1984); Pickett et al. (1998).
Thus, the parameter $U_{eff}=U-J_{H}$, which roughly determines the splitting
between the lower and upper Hubbard bands, varied between 0 eV and 3.0 eV. We
adjusted the value of $U$ to achieve the best agreement with the experiment.

In the RIXS process an electron is promoted from a core level to an
intermediate state, leaving a core hole. As a result, the electronic structure
of this state differs from that of the ground state. In order to reproduce the
experimental spectrum the self-consistent calculations should be carried out
including a core hole. Usually, the core-hole effect has no impact on the
shape of XAS at the $L_{2,3}$ edges of 5$d$ systems and just a minor effect on
the XMCD spectra at these edges Antonov et al. (2004). However, the core hole
has a strong effect on the RIXS spectra in transition metal compounds Antonov
et al. (2022a, b), therefore we take it into account.

##  III Electronic and magnetic structures

Figure 2: (Color online) The AFM ordering with Ir moments parallel to $c$ axis
used in GGA+SO+$U$ calculations.

We performed GGA, GGA+SO, and GGA+SO+$U$ calculations of the electronic and
magnetic structures of Sr2IrO4 for the experimental crystal structure Shimura
et al. (1995) and AFM ordering along the $c$ direction (Fig. 2). Although
Sr2IrO4 has a canted AFM structure Kim et al. (2009), test calculations for a
non-collinear AFM structure gave almost identical optical, XA, XMCD, and RIXS
spectra. Fig. 3 (the upper panel) presents the energy band structure of
Sr2IrO4 in the energy range of the Ir $t_{2g}$ manifold from $-$2.1 to 1.0 eV,
calculated in the fully relativistic Dirac GGA+SO approximation. The GGA+SO
bands are presented by circles proportional in size to their orbital character
projected onto the basis set of Ir $d_{3/2}$ (blue) and $d_{5/2}$ (red)
states. The strong SOC splits the $t_{2g}$ manifold into a lower $J_{eff}$ =
3/2 quartet and an upper $J_{eff}$ = 1/2 doublet. The functions of the
$J_{eff}$ = 3/2 quartet are dominated by $d_{3/2}$ states with some weight of
$d_{5/2}$ ones, which is determined by the relative strength of SOC and
crystal-field splitting. The $J_{eff}$ = 1/2 functions are almost completely
given by the linear combinations of $d_{5/2}$ states. This allows one to
identify the bands with pure $d_{5/2}$ character as originating from the
$J_{eff}$ = 1/2 states. The GGA+SO approach produces a metallic state in
Sr2IrO4. The GGA+SO+$U$ approximation shifts the occupied and empty $t_{2g}$
bands downward and upward, respectively, by $U_{eff}$/2 producing a dielectric
ground state (the lower panel of Fig. 3). The energy gap is increased with
increasing Hubbard $U$.

Figure 3: (Color online) The $t_{2g}$ energy band structure of Sr2IrO4
calculated in the fully relativistic Dirac GGA+SO approximation (the upper
panel). The bands crossing the Fermi level which have almost pure $d_{5/2}$
character (open red circles) are formed by $t_{2g}$ states with $j_{eff}$ =
1/2. The lower panel presents the $t_{2g}$ energy bands calculated in the
GGA+SO+$U$ approximation with $U_{eff}$ = 1.2 eV.  Figure 4: (Color online)
The experimentally measured real part of the optical conductivity (open
magenta circles), $\sigma_{1xx}$, Pröpper et al. (2016) (in 103 $\Omega^{-1}$
cm-1) in Sr2IrO4 of the in-plane response compared with the theoretical
spectra calculated in the GGA+SO+$U$ approximation for different $U_{eff}$
values.

Figure 4 presents the experimentally measured real part of the optical
conductivity (open magenta circles), $\sigma_{1xx}$, Pröpper et al. (2016) for
the energy below 2 eV in Sr2IrO4 compared with the theoretical spectra
calculated in the GGA+SO+$U$ approximation for different $U_{eff}$ values. The
experimental optical absorption consists of two peaks at around 0.5 and 1.0
eV. We found that the low energy peak is derived from transitions between
initial and final bands formed by pure $J_{eff}$ = 1/2 states near band
boundaries, e.g., around the X point or the P-N high-symmetry line. The
antiferromagnetic ordering of Ir moments within the $ab$ plane stabilized by
the on-site Coulomb repulsion $U$ causes a gap opening near the zone boundary
between two pairs of bands which show nearly parallel dispersion what insures
high joint density of states for the interband transitions responsible for the
low energy peak. This is in line with previous theoretical calculations Zhang
et al. (2013); Pröpper et al. (2013) and experimental photoemission results
Wang et al. (2013). The high energy peak located around 1 eV is dominated by a
contribution from transitions with $J_{eff}$ = 3/2 initial states. Our
calculations give the lower absorption peak about twice as strong as the
higher energy one, while in the experimental spectra the strength is
approximately the same for both. A similar trend was observed also by Pröpper
et al. Pröpper et al. (2013) and Kim et al. Kim et al. (2012b). The later
authors relate this to an interband mixing of $J_{eff}$ = 3/2 and $J_{eff}$ =
1/2 states, which reflects the itinerancy of the system, i.e., the
hybridization of Ir 5$d$ states via neighboring oxygen 2$p$ states. The best
agreement between the calculated and experimentally measured energy positions
of the optical absorption peaks can be achieved for $U_{eff}$ = 1.2 eV.

Figures 5 and 6 present the energy band structure and partial density of
states (DOS), respectively, in Sr2IrO4 calculated in the GGA+SO+$U$
approximation with $U_{eff}$ = 1.2 eV. Five electrons occupy the $t_{2g}$-type
low energy band (LEB) manifold in the energy interval from $-$1.5 eV to
$E_{F}$ in Sr2IrO4. The empty $t_{2g}$ states [the upper energy band (UEB)]
consist of one peak and occupy the energy range from 0.41 eV to 0.82 eV (see
Fig. 6). The $e_{g}$-type states of Ir are distributed far above the Fermi
level from 2.2 eV to 5.1 eV. The 4$d$ states of Sr ions are mostly situated
above the Fermi level from 5.3 to 10.8 eV.

Figure 5: (Color online) The energy band structure and total DOS [in
states/(cell eV)] of Sr2IrO4 calculated with taking into account Coulomb
correlations in the GGA+SO+$U$ approximation ($U_{eff}$ = 1.2 eV).  Figure 6:
(Color online) The partial DOSs for Sr2IrO4 calculated in the GGA+SO+$U$
($U_{eff}$= 1.2 eV) approximation.

The electronic structure of apical O1 and in-plane O2 ions significantly
differ from each other. The apical O1 2$s$ states consist of two very narrow
peaks situated at $-$15.1 and $-$14.4 eV. The in-plane O2 2$s$ states possess
a relatively wider two peak structure from $-$17.9 to $-$16.5 eV. The O1 2$p$
states are situated just below Ir LEB between $-$3.5 and $-$2.1 eV. There is
also a narrow peak at $-$5.2 eV. The in-plane O2 2$p$ states occupy a
relatively wide energy interval from $-$8.7 to $-$3.6 eV. The small peaks in
the close vicinity of the Fermi level from $-$1.5 eV to $E_{F}$ and from 0.41
eV to 0.82 eV are due to the strong hybridization between O 2$p$ and Ir
$t_{2g}$ LEB and UEB, respectively.

The occupation number of 5$d$ electrons in the Ir atomic sphere in Sr2IrO4 is
equal to 6.3, which is much larger than the expected value of five $t_{2g}$
electrons. The excessive charge is provided by the ”tails” of oxygen 2$p$
states. These 5$d_{\rm{O}}$ states are located at the bottom of oxygen 2$p$
states from $-$8.7 eV to $-$3.1 eV and play an essential role in the RIXS
spectrum at the Ir $L_{3}$ edge (see Section IV).

The theoretically calculated spin $M_{s}$, orbital $M_{s}$, and total
$M_{total}$ magnetic moments using the GGA+SO+$U$ approach ($U_{eff}$ = 1.2
eV) for the AFM solution are equal to 0.2647 $\mu_{\text{B}}$, 0.4447
$\mu_{\text{B}}$, and 0.7094 $\mu_{\text{B}}$, respectively. The spin and
orbital magnetic moments at the Sr site are relatively small ($M_{s}$ = 0.0007
$\mu_{\text{B}}$ and $M_{l}$ = 0.0015 $\mu_{\text{B}}$). The magnetic moments
for apical O1 ions are equal to $M_{s}$ = 0.0294 $\mu_{\text{B}}$, $M_{l}$ =
0.0254 $\mu_{\text{B}}$. For in-plane O2 ions the magnetic moments almost
vanish.

##  IV Ir XMCD and RIXS spectra

Figure 7 presents the experimentally measured x-ray absorption spectra (the
upper panel) and XMCD spectra (the lower panel) at the Ir $L_{2,3}$ edges for
Sr2IrO4 Haskel et al. (2012) (open circles) compared with the theoretically
calculated ones in the GGA+SO+$U$ ($U_{eff}$ = 1.2 eV) approximation (full
blue curves). The theoretically calculated Ir $L_{2,3}$ XA and XMCD spectra
are in good agreement with the experiment. The isotropic XA spectra are
dominated by the empty $e_{g}$ states with a smaller contribution from the
empty $t_{2g}$ orbitals at lower energy. The XMCD spectra, however, mainly
come from the $t_{2g}$ orbitals ($J_{eff}$ = 1/2). This results in a shift
between the maxima of the XA and XMCD spectra.

Figure 7: (Color online) The experimental x-ray absorption (upper panels) and
XMCD spectra (lower panels) at the Ir $L_{2,3}$ edges in the Sr2IrO4 thin film
(magenta circles) Haskel et al. (2012) measured at 6 K under a 0.8 T magnetic
field compared with the theoretically calculated spectra in the GGA+SO+$U$
approximation (full blue curves). The dotted black curves in the upper panels
show the background scattering intensity.

Due to the importance of SOC effects in iridates, it is natural to quantify
the strength of the SO interactions in these compounds. One method of
accomplishing this is provided by the x-ray absorption spectroscopy. Van der
Laan and Thole showed that the so-called branching ratio BR =
$I_{L_{3}}/I_{L_{2}}$ ($I_{L_{2,3}}$ is the integrated intensity of the
isotropic XAS at the $L_{2,3}$ edges) is an important quantity in the study of
5$d$ oxides related to the SO interaction van der Laan and Thole (1988). The
BR is directly related to the ground-state expectation value of the angular
part of the spin-orbit coupling $<{\bf L\cdot S}>$ through BR = $(2+r)/(1-r)$,
with $r$= $<{\bf L\cdot S}>/n_{h}$ and $n_{h}$ is the number of holes in 5$d$
states van der Laan and Thole (1988). As a result, XAS provides a direct probe
of SO interactions, which is complementary to other techniques such as the
magnetic susceptibility, electron paramagnetic resonance, and Mossbauer
spectroscopy (which probe SOC through the value of the Lande g-factor). In the
limit of negligible SOC effects the statistical branching ratio BR = 2, and
the $L_{3}$ white line is twice the size of the $L_{2}$ feature van der Laan
and Thole (1988). The measured BR in Sr2IrO4 is close to 4.1 Haskel et al.
(2012), which differs significantly from the statistical BR = 2 in the absence
of orbital magnetization in 5$d$ states. A strong deviation from 2 indicates a
strong coupling between the local orbital and spin moments. Our DFT
calculations produce BR = 3.56 for the GGA+SO+$U$ ($U_{eff}$ = 1.2 eV)
approximation which is rather close to the experimental data of Haskel et al.
Haskel et al. (2012). It should be mentioned that the effect of Coulomb
correlations changes the energy band structure of transition metal compounds
in two ways. First, $d$ occupied states are shifted downward by $U_{eff}$/2
and empty $d$ states are shifted upward by this value relative to the Fermi
energy. Second, the Coulomb correlations enhance an effective spin-orbit
coupling constant Liu et al. (2008). The relative influence of such effect is
increased in a row of 5$d\rightarrow 4d\rightarrow 3d$ transition metal
compounds due to an increase of Hubbard $U$ and a decrease of the atomic SO
coupling constant $\lambda_{SO}$.

Figure 8: (Color online) The experimental RIXS spectrum (open green circles)
measured by Ishii et al. Ishii et al. (2011) at the Ir $L_{3}$ edge in Sr2IrO4
compared with the theoretically calculated one in the GGA+SO+$U$ approximation
($U_{eff}$=1.2 eV). Inset: The experimental RIXS spectrum (open magenta
circles) measured by Clancy et al. Clancy et al. (2023) at the Ir $L_{3}$ edge
for $t_{2g}$ $\rightarrow$ $t_{2g}$ transitions in Sr2IrO4 compared with the
theoretically calculated one in the GGA+SO+$U$ approach ($U_{eff}$=1.2 eV).

The RIXS spectra at Ir $L_{2,3}$ occur from a local excitation between the
filled and empty 5$d$ states. More precisely, the incoming photon excites a
2$p_{1/2}$ core electron ($L_{2}$ spectrum) or 2$p_{3/2}$ one ($L_{3}$
spectrum) into an empty 5$d$ state what is followed by the de-excitation from
the occupied 5$d$ state into the core level. Because of the dipole selection
rules, apart from 6$s_{1/2}$-states (which have a small contribution to RIXS
due to relatively small 2$p$ $\rightarrow$ 6$s$ matrix elements Antonov et al.
(2004)) only 5$d_{3/2}$-states occur for $L_{2}$ RIXS, whereas for $L_{3}$
RIXS 5$d_{5/2}$-states also contribute. Although the 2$p_{3/2}$ $\rightarrow$
5$d_{3/2}$ radial matrix elements are only slightly smaller than the
2$p_{3/2}$ $\rightarrow$ 5$d_{5/2}$ ones, the angular matrix elements strongly
suppress the 2$p_{3/2}$ $\rightarrow$ 5$d_{3/2}$ contribution Antonov et al.
(2004). Therefore, the RIXS spectrum at the Ir $L_{3}$ edge can be viewed as
interband transitions between 5$d_{5/2}$ states.

Figure 8 shows the experimental RIXS spectrum (open green circles) measured by
Ishii et al. Ishii et al. (2011) at the Ir $L_{3}$ edge in Sr2IrO4 compared
with the theoretically calculated one in the GGA+SO+$U$ approximation
($U_{eff}$=1.2 eV). The Ir $L_{3}$ RIXS spectrum consists of two peaks below 5
eV. We found that the low energy peak corresponds to intra-$t_{2g}$
excitations. This fine structure has a two peak structure in our calculations
but the measurements of Ishii et al. Ishii et al. (2011) show only one peak.
However, Clancy et al. Clancy et al. (2023) using higher resolution were able
to distinguish two peaks (see the inset in Fig. 8). The low energy peak at 0.5
eV is due to interband transitions between occupied and empty Ir $J_{eff}$ =
1/2 states (the red curves in the lower panel of Fig. 3). These transitions
also contribute to the second high energy peak at around 0.7 eV together with
$J_{3/2}$ $\rightarrow$ $J_{1/2}$ transitions (the green curve in the inset of
Fig. 8). The intensive peak at around 3.4 eV (the red curve in Fig. 8) is due
to $t_{2g}\rightarrow e_{g}$ transitions. The next fine structure from 4.5 eV
to 12 eV (the magenta curve) is due to 5$d_{\rm{O}}$ $\rightarrow$ $t_{2g}$,
$e_{g}$ transitions.

Figure 9: (Color online) The RIXS spectra as a function of incident photon
energy $E_{i}$ calculated at the Ir $L_{3}$ edge in Sr2IrO4 with the momentum
transfer vector Q = (0, 0, 33).

Figure 9 shows the Ir $L_{3}$ RIXS spectrum as a function of incident photon
energy $E_{i}$ above the corresponding edge with the momentum transfer vector
Q = (0, 0, 33). We found that the low energy fine structure corresponding to
intra-$t_{2g}$ excitations is slightly decreased when the incident photon
energy changes in the interval of 2 eV above the Ir $L_{3}$ edge, while the
high energy peak corresponding to the $t_{2g}\rightarrow e_{g}$ transitions is
monotonically increased (the upper panel of Fig. 9). A similar trend was
observed in the measurements of Ishii et al. Ishii et al. (2011). They
measured the RIXS spectra of Sr2IrO4 at two representative incident photon
energies in the interval of 3 eV and discovered that with increasing $E_{i}$
the low energy peak is decreased but the high energy peak is increased [for
the same momentum transfer vector Q = (0, 0, 33)]. The lower panel of Fig. 9
shows the Ir $L_{3}$ RIXS spectrum as a function of incident photon energy in
the larger energy interval of 18 eV. With increasing $E_{i}$ the low energy
peak is steadily decreased, but the high energy peak shows more complex
behavior. First, the intensity of the peak is increased, but then decreased
with changing the relative intensity of the two peaks at 3 and 5 eV .

Figure 10: (Color online) The RIXS spectra at the Ir $L_{3}$ edge in Sr2IrO4
calculated as a function of $Q_{y}$ (the upper panel) and $Q_{z}$ (the lower
panel) with the momentum transfer vector Q = (0, Qy, Qz) for incident photon
energy $\hbar\omega_{in}$ = 11210 eV.

It is widely believed that $d-d$ excitations show only small momentum transfer
vector Q dependence in 5$d$ transition metal compounds Liu et al. (2012);
Krajewska et al. (2020). In particular, Sr2IrO4 has a layered-perovskite
structure, therefore, the momentum dependence along the $c$ axis is expected
to be small, as in high-$T_{c}$ cuprates Ishii et al. (2005). Indeed, as we
see in the lower panel of Fig. 10, the RIXS spectra are almost identical for
the transfer vectors Q = (0, 0, 33) and (0, 0, 25). Similar dependence was
observed also in the measurements of Ishii et al. Ishii et al. (2011). The
upper panel of Fig. 10 shows the RIXS spectra at the Ir $L_{3}$ edge in
Sr2IrO4 calculated as a function of $Q_{y}$ with the momentum transfer vector
Q = (0, Qy, 33) for incident photon energy $\hbar\omega_{in}$ = 11210 eV. We
found that with increasing Qy the first low energy peak is increased and the
high energy fine structure is decreased. Analyzing Fig. 10 we can conclude
that the momentum dependence of the excitations in Sr2IrO4 is rather small as
it was earlier observed in other iridates such as Sr3CuIrO6 Liu et al. (2012)
or In2Ir2O7 Krajewska et al. (2020).

##  V O XAS and RIXS spectra

Figure 11: (Color online) (The upper panel) The experimental RIXS spectrum
(open magenta circles) measured by Liu et al. Liu et al. (2015) at the O $K$
edge in Sr2IrO4 compared with the theoretically calculated ones in the
GGA+SO+$U$+SIC approximation ($U_{eff}$ = 1.2 eV, $V_{O_{2p}}$ = $-$3 eV).
(The lower panel) The experimental RIXS spectrum (open magenta circles)
measured by Liu et al. Liu et al. (2015) at the O $K$ edge in Sr2IrO4 compared
with the theoretically calculated ones in the GGA+SO+$U$+SIC approach
($U_{eff}$ = 1.2 eV) for different parameters $V_{l}$.

The RIXS spectra at the O $K$ edge in Sr2IrO4 were measured by Liu et al. Liu
et al. (2015), Lu et al. Lu et al. (2018), Paris et al. Paris et al. (2020),
and Kim et al. Kim et al. (2023). The last three investigations concentrate on
the analysis of low energy excitations below 1.2 eV. Liu et al. Liu et al.
(2015) present the RIXS spectrum up to 12 eV using circular and $\pi$
polarizations of the incident beam. The O $K$ RIXS spectrum consists of a peak
centered at zero energy loss, which comprises the elastic line and other low-
energy features such as phonons, magnons, etc., and three major inelastic
excitations at 0.7 eV, 3.5 eV, and around 6.2. We found that the first low
energy feature is due to the interband transitions between occupied and empty
O${}_{t_{2g}}$ states, which appear as a result of the strong hybridization
between oxygen 2$p$ states with Ir $t_{2g}$ LEB and UEB in the close vicinity
of the Fermi level (see Fig. 6), therefore, the oxygen $K$ RIXS spectroscopy
can be used for the estimation of the energy band gap and positions of Ir 5$d$
Hubbard bands. The next two peaks at around 3.5 and 6.2 eV reflect the
interband transitions from the occupied O 2$p$ states and the empty oxygen
states which originate from the hybridization with Ir $t_{2g}$ and $e_{g}$
states, respectively. We found that the theory reproduces well the shape and
energy position of the low energy feature, but the second and the third peaks
are shifted towards smaller energy in comparison with the experimental
measurements. It means that the DFT calculations cannot produce the correct
energy position of the oxygen 2$p$ bands. These bands are almost fully
occupied in Sr2IrO4, therefore, they cannot be described by the GGA+$U$
method. To reproduce the correct energy position of the oxygen 2$p$ band in
Sr2IrO4 we used a self-interaction-like correction procedure as proposed by
Kaneko et al. Kaneko et al. (2013), where the valence bands are shifted
downwards by adding a SIC-like orbital-dependent potential $V_{l}$ into the
Hamiltonian. We used $V_{l}$ as a parameter an adjusted it to produce the
correct energy position of the oxygen 2$p$ bands. We found that the best
agreement with the experiment can be achieved for $V_{O_{2p}}$ = $-$3.0 eV
(see the lower panel of Fig. 11).

Figure 12: (Color online) The valence band photoemission spectrum of Sr2IrO4
King et al. (2013) compared with total DOS calculated in the GGA+SO+$U$+SIC
approach ($U_{eff}$ = 1.2 eV) for different parameters $V_{l}$.

Figure 12 presents the valence band photoemission spectrum of Sr2IrO4 King et
al. (2013) compared with total DOS calculated in the GGA+SO+$U$+SIC approach
($U_{eff}$ = 1.2 eV). Two dominant peaks at around $-$1 and $-$3.2 eV are
observed, which might be attributed to the photoemmision from Ir $t_{2g}$ and
oxygen 2$p$ states, respectively King et al. (2013). Like in the case of O $K$
RIXS spectrum, the GGA+SO+$U$ approach cannot reproduce the correct energy
position of the peak at $-$3.2 eV. However, the SIC-like approach with
$V_{O_{2p}}$ = $-$3 eV improves the situation.

Figure 13: (Color online) The RIXS spectra as a function of incident photon
energy calculated at the O $K$ edge in Sr2IrO4 with circular polarization.
Figure 14: (Color online) The RIXS spectra at the O $K$ edge in Sr2IrO4
calculated as a function of $Q_{y}$ (the upper panel) and $Q_{z}$ (the lower
panel) with the momentum transfer vector Q = (0, Qy, Qz) for incident photon
energy $\hbar\omega_{in}$ = 529 eV.

Figure 13 presents the RIXS spectra as a function of incident photon energy
$E_{i}$ calculated at the O $K$ edge in Sr2IrO4 with circular polarization. We
found much stronger dependence on the incident photon energy in the case of
the O $K$ RIXS spectrum in comparison with the corresponding dependence at the
Ir $L_{3}$ edge (compare Figs. 9 and 13). With increasing the incident photon
energy both peaks at 0.7 eV and 3.5 eV are increased, the later one is
increased dramatically. And this occurs in a small energy interval for $E_{i}$
of 0.6 eV.

Figure 14 shows the RIXS spectra at the O $K$ edge in Sr2IrO4 calculated as a
function of $Q_{y}$ (the upper panel) and $Q_{z}$ (the lower panel) with the
momentum transfer vector Q = (0, Qy, Qz). With decreasing parameters Qx and Qz
the intensity of the major peak at 3.5 eV is decreased but the low energy peak
at 0.7 eV is increased. There is also a strong change in the shape of the low
energy peak at 0.7 eV and the third peak at around 6.2 eV with the change of
the parameter Qz.

Figure 15: (Color online) The experimental O $K$ polarization dependent XA
spectra (open magenta circles) Kim et al. (2008) in Sr2IrO4 compared with the
theoretically calculated ones in the GGA+SO+$U$+SIC approach ($U_{eff}$ = 1.2
eV, $V_{O_{2p}}$ = $-$3 eV).

Figure 15 presents the experimental O $K$ polarization dependent x-ray
absorption spectra (open magenta circles) Kim et al. (2008) compared with the
theoretically calculated ones in the GGA+SO+$U$+SIC approach ($U_{eff}$ = 1.2
eV, $V_{O_{2p}}$ = $-$3 eV). Due to the quasi-two-dimensional structure of
Sr2IrO4 there is strong anisotropy in the x-ray absorption spectra. There are
two small peaks $a$ and $b$ at 529.1 and 529.9 eV and larger peak $c$ at
around 531.8 eV for the E $\perp$ $c$ polarization and only two peaks $A$ and
$B$ for E $\parallel$ $c$ at 529.9 and 531.3 eV. We found that the low energy
peak $a$ for E $\perp$ $c$ and the large peak $B$ for E $\parallel$ $c$ are
derived from the apical oxygens O1. The peaks $b$ and $A$ are due to the 1$s$
$\rightarrow$ 2$p$ x-ray absorption on the in-plane O2 oxygens. The theory
relatively well reproduces the experimentally measured XA spectra.

##  VI Conclusions

The electronic and magnetic structures of Sr2IrO4 were investigated
theoretically in the frame of the fully relativistic spin-polarized Dirac LMTO
band-structure method in order to understand the importance of Coulomb
interaction and spin-orbit coupling. We also present comprehensive theoretical
calculations of the XA, XMCD, and RIXS spectra at the Ir $L_{2,3}$ and oxygen
$K$ edges.

The strong SOC splits the $t_{2g}$ manifold into a lower $J_{eff}$ = 3/2
quartet and an upper $J_{eff}$ = 1/2 doublet. The functions of the $J_{eff}$ =
3/2 quartet are dominated by $d_{3/2}$ states with some weight of $d_{5/2}$
ones, which is determined by the relative strength of the SOC and crystal-
field splitting. The $J_{eff}$ = 1/2 functions are almost completely given by
linear combinations of $d_{5/2}$ states. This allows one to identify the bands
with pure $d_{5/2}$ character as originating from the $J_{eff}$ = 1/2 states.
The GGA+SO approach produces a metallic state in Sr2IrO4. The GGA+SO+$U$
approximation shifts the occupied and empty $t_{2g}$ bands downward and
upward, respectively, by $U_{eff}$/2 producing a dielectric ground state in
Sr2IrO4.

We found that the best agreement between the calculated and experimentally
measured optical conductivity spectrum can be achieved for $U_{eff}$ = 1.2 eV.
The experimental optical absorption consists of two peaks at around 0.5 and
1.0 eV. The low energy peak is derived from transitions between initial and
final bands formed by pure $J_{eff}$ = 1/2 states near band boundaries. The
high energy peak located around 1 eV is dominated by a contribution from
transitions with $J_{eff}$ = 3/2 initial states.

The theoretically calculated Ir $L_{2,3}$ XAS and XMCD spectra are in good
agreement with the experiment. The isotropic XA spectra are dominated by the
empty $e_{g}$ states with a smaller contribution from the empty $t_{2g}$
orbitals at lower energy. The XMCD spectra, however, mainly come from the
$t_{2g}$ orbitals ($J_{eff}$ = 1/2). This results in a shift between the
maxima of the XA and XMCD spectra. The ratio BR = $I_{L_{3}}/I_{L_{2}}$ is an
important quantity in the study of the SO interaction in 5$d$ oxides. It is
directly related to the ground-state expectation value of the angular part of
the spin-orbit coupling. Our DFT calculations produce BR = 3.56 which is
rather close to the experimental value of 4.1 Haskel et al. (2012). These
values differ significantly from the statistical BR = 2 in the absence of
orbital magnetization in 5$d$ states. A strong deviation from 2 indicates a
strong coupling between the local orbital and spin moments.

The Ir $L_{3}$ RIXS spectrum consists of two peaks below 5 eV. We found that
the low energy peak corresponds to intra-$t_{2g}$ excitations. This fine
structure has a two peak structure. The low energy peak at 0.5 eV is due to
interband transitions between occupied and empty Ir $J_{eff}$ = 1/2 states.
These transitions also contribute to the second high energy peak at around 0.7
eV together with $J_{3/2}$ $\rightarrow$ $J_{1/2}$ transitions. The intensive
peak at around 3.4 eV is due to $t_{2g}\rightarrow e_{g}$ transitions. The
next fine structure from 4.5 eV to 12 eV is due to 5$d_{\rm{O}}$ $\rightarrow$
$t_{2g}$, $e_{g}$ transitions.

We investigated theoretically the influence of the incident photon energy
$E_{i}$ and the momentum transfer vector Q on the shape of the Ir $L_{3}$ RIXS
spectrum. We found that with inreasing of $E_{i}$ in the interval of 2 eV
above the Ir $L_{3}$ edge the low energy fine structure corresponding to
intra-$t_{2g}$ excitations is decreased and the high energy peak corresponding
to the $t_{2g}\rightarrow e_{g}$ transitions is monotonically increased in
agreement with the measurements of Ishii et al. Ishii et al. (2011). The
momentum dependence of Ir $L_{3}$ RIXS was found to be relatively small. With
increasing Qy the first low energy peak is increased and the high energy fine
structure is slightly decreased. The variation of parameter Qz almost does not
influence the RIXS spectrum at the Ir $L_{3}$ edge.

The RIXS spectrum of Sr2IrO4 at the O $K$ edge consists of three major
inelastic excitations at 0.7 eV, 3.5 eV, and around 6.2 eV. We found that the
first low energy feature is due to the interband transitions between occupied
and empty O${}_{t_{2g}}$ states, which appear as a result of the strong
hybridization between oxygen 2$p$ states with Ir $t_{2g}$ LEB and UEB in the
close vicinity of the Fermi level. The next two peaks at around 3.5 and 6.2 eV
reflect the interband transitions from the occupied O 2$p$ states to the empty
oxygen states which originate from the hybridization with Ir $t_{2g}$ and
$e_{g}$ states, respectively. We found that the theory reproduces well the
shape and energy position of the low energy feature, but the second and the
third peaks are shifted towards smaller energy in comparison with the
experimental measurements. It means that the DFT calculations cannot produce
the correct position of oxygen 2$p$ bands. To reproduce the correct energy
position of the oxygen 2$p$ band in Sr2IrO4 we used a self-interaction-like
correction procedure. We added a SIC-like orbital-dependent potential $V_{l}$
into the Hamiltonian and found that the best agreement with the experiment can
be achieved for $V_{O_{2p}}$ = $-$3.0 eV. We found that the dependence of the
RIXS spectrum at the oxygen $K$ edge on the incident photon energy and the
momentum transfer vector Q is much stronger in comparison with the
correspondent dependency at the Ir $L_{3}$ edge.

Due to the quasi-two-dimensional structure of Sr2IrO4 there is strong
anisotropy in the x-ray absorption spectra at the oxygen $K$ edge. There are
two small peaks $a$ and $b$ at 529.1 and 529.9 eV and a larger peak $c$ at
around 531.8 eV for the E $\perp$ $c$ polarization and only two peaks $A$ and
$B$ for E $\parallel$ $c$ at 529.9 and 531.3 eV. We found that the low energy
peak $a$ for E $\perp$ $c$ and the large peak $B$ for E $\parallel$ $c$ are
derived from the apical oxygens O1. The peaks $b$ and $A$ are due to the 1$s$
$\rightarrow$ 2$p$ x-ray absorption on the in-plane O2 oxygens. The theory
relatively well reproduces the experimentally measured XA spectra.

## Acknowledgments

We are thankful to Dr. Alexander Yaresko from the Max-Planck-Institute FKF in
Stuttgart for helpful discussions.

## References

  * Jackeli and Khaliullin (2009) G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 017205 (2009).
  * Chen et al. (2010) G. Chen, R. Pereira, and L. Balents, Phys. Rev. B 82, 174440 (2010).
  * Witczak-Krempa et al. (2014) W. Witczak-Krempa, G. Chen, Y. B. Kim, and L. Balents, Annu. Rev. Condens. Matter Phys. 5, 57 (2014).
  * Kim et al. (2008) B. J. Kim, H. Jin, S. J. Moon, J.-Y. Kim, B.-G. Park, C. S. Leem, J. Yu, T. W. Noh, C. Kim, S.-J. Oh, et al., Phys. Rev. Lett. 101, 076402 (2008).
  * Martins et al. (2011) C. Martins, M. Aichhorn, L. Vaugier, and S. Biermann, Phys. Rev. Lett. 107, 266404 (2011).
  * Antonov et al. (2018) V. N. Antonov, S. Uba, and L. Uba, Phys. Rev. B 98, 245113 (2018).
  * Qi and Zhang (2010) X.-L. Qi and S.-C. Zhang, Physics Today 63, 33 (2010).
  * Ando (2013) Y. Ando, J. Phys. Soc. Jpn. 82, 102001 (2013).
  * Wehling et al. (2014) T. O. Wehling, A. Black-Schafferc, and A. Balatsky, Adv. Phys. 63, 1 (2014).
  * Bansi et al. (2016) A. Bansi, L. H, and T. Das, Rev. Mod. Phys. 88, 021004 (2016).
  * Kim et al. (2009) B. J. Kim, H. Ohsumi, T. Komesu, S. Sakai, T. Morita, H. Takagi, and T. Arima, Science 323, 1329 (2009).
  * Watanabe et al. (2010) H. Watanabe, T. Shirakawa, and S. Yunoki, Phys. Rev. Lett. 105, 216410 (2010).
  * Witczak-Krempa and Kim (2012) W. Witczak-Krempa and Y. B. Kim, Phys. Rev. B 85, 045124 (2012).
  * Go et al. (2012) A. Go, W. Witczak-Krempa, G. S. Jeon, K. Park, and Y. B. Kim, Phys. Rev. Lett. 109, 066401 (2012).
  * Sushkov et al. (2015) A. B. Sushkov, J. B. Hofmann, G. S. Jenkins, J. Ishikawa, S. Nakatsuji, S. DasSarma, and H. D. Drew, Phys. Rev. B 92, 241108 (2015).
  * Kimchi et al. (2014) I. Kimchi, J. G. Analytis, and A. Vishwanath, Phys. Rev. B 90, 205126 (2014).
  * Kim et al. (2015) B. J. Kim, H. Jin, S. J. Moon, J.-Y. Kim, B.-G. Park, C. S. Leem, J. Yu, T. W. Noh, C. Kim, S.-J.-H. Park, et al., Phys. Rev. B 115, 176402 (2015).
  * Yan et al. (2015) Y. Yan, M. Ren, H. Xu, B. Xie, R. Tao, H. Choi, N. Lee, Y. Choi, T. Zhang, and D. Feng, Phys. Rev. X 5, 041018 (2015).
  * Battisti et al. (2017) I. Battisti, K. M. Bastiaans, V. Fedoseev, A. de la Torre, N. Iliopoulos, A. Tamai, E. C. Hunter, R. S. Perry, J. Zaanen, F. Baumberger, et al., Nat. Phys. 13, 21 (2017).
  * Damascelli et al. (2003) A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod. Phys. 75, 473 (2003).
  * Kim et al. (2014) Y. K. Kim, O. Krupin, J. D. Denlinger, A. Bostwick, E. Rotenberg, Q. Zhao, J. F. Mitchell, J. W. Allen, and B. J. Kim, Science 345, 187 (2014).
  * Kim et al. (2016) Y. K. Kim, N. H. Sung, J. D. Denlinger, and B. J. Kim, Nat. Phys. 12, 37 (2016).
  * Terashima et al. (2017) K. Terashima, M. Sunagawa, H. Fujiwara, T. Fukura, M. Fujii, K. Okada, K. Horigane, K. Kobayashi, R. Horie, J. Akimitsu, et al., Phys. Rev. B 96, 041106(R) (2017).
  * Hu et al. (2019) Y. Hu, X. Chen, S.-T. Peng, C. Lane, M. Matzelle, Z.-L. Sun, M. Hashimoto, D.-H. Lu, E. Schwier, M. Arita, et al., Phys. Rev. Lett. 123, 216402 (2019).
  * Kim et al. (2012a) J. Kim, D. Casa, M. H. Upton, T. Gog, Y.-J. Kim, J. F. Mitchell, M. van Veenendaal, M. Daghofer, J. van den Brink, G. Khaliullin, et al., Phys. Rev. Lett. 108, 177003 (2012a).
  * Clancy et al. (2019) J. P. Clancy, H. Gretarsson, M. H. Upton, J. Kim, G. Cao, and Y.-J. Kim, Phys. Rev. B 100, 104414 (2019).
  * Bertinshaw et al. (2019) J. Bertinshaw, Y. Kim, G. Khaliullin, and B. Kim, Annu. Rev. Condens. Matter Phys. 7, 315 (2019).
  * Kao et al. (1996) C.-C. Kao, W. A. L. Caliebe, J. B. Hastings, and J.-M. Gillet, Phys. Rev. B 54, 16361 (1996).
  * Ament et al. (2011) L. J. P. Ament, M. van Veenendaal, T. P. Devereaux, J. P. Hill, and J. van den Brink, Rev. Mod. Phys. 83, 705 (2011).
  * Liu et al. (2012) X. Liu, V. M. Katukuri, L. Hozoi, W.-G. Yin, M. P. M. Dean, M. H. Upton, J. Kim, D. Casa, A. Said, T. Gog, et al., Phys. Rev. Lett. 109, 157401 (2012).
  * Hozoi et al. (2014) L. Hozoi, H. Gretarsson, J. P. Clancy, B.-G. Jeon, B. Lee, K. H. Kim, V. Yushankhai, P. Fulde, D. Casa, T. Gog, et al., Phys. Rev. B 89, 115111 (2014).
  * Clancy et al. (2016) J. P. Clancy, H. Gretarsson, E. K. H. Lee, D. Tian, J. Kim, M. H. Upton, D. Casa, T. Gog, Z. Islam, B.-G. Jeon, et al., Phys. Rev. B 94, 024408 (2016).
  * Nag et al. (2018) A. Nag, S. Bhowal, A. Chakraborty, M. M. Sala, A. Efimenko, F. Bert, P. K. Biswas, A. D. Hillier, M. Itoh, S. D. Kaushik, et al., Phys. Rev. B 98, 014431 (2018).
  * Takayama et al. (2019) T. Takayama, A. Krajewska, A. S. Gibbs, A. N. Yaresko, H. Ishii, H. Yamaoka, K. Ishii, N. Hiraoka, N. P. Funnell, C. L. Bull, et al., Phys. Rev. B 99, 125127 (2019).
  * Aczel et al. (2022) A. A. Aczel, Q. Chen, J. P. Clancy, C. dela Cruz, D. R. i Plessis, G. J. MacDougall, C. J. Pollock, M. H. Upton, T. J. Williams, N. LaManna, et al., Phys. Rev. Mater. 6, 094409 (2022).
  * Calder et al. (2016) S. Calder, J. Vale, N. Bogdanov, X. Liu, C. Donnerer, M. Upton, D. Casa, A. Said, M. Lumsden, Z. Zhao, et al., Nature Commun. 7, 11651 (2016).
  * Taylor et al. (2017) A. E. Taylor, S. Calder, R. Morrow, H. L. Feng, M. H. Upton, M. D. Lumsden, K. Yamaura, P. M. Woodward, and A. D. Christianson, Phys. Rev. Lett. 118, 207202 (2017).
  * Calder et al. (2017) S. Calder, J. G. Vale, N. Bogdanov, C. Donnerer, D. Pincini, M. M. Sala, X. Liu, M. H. Upton, D. Casa, Y. G. Shi, et al., Phys. Rev. B 95, 020413(R) (2017).
  * Liu et al. (2015) X. Liu, M. P. M. Dean, J. Liu, S. G. Chiuzbaian, N. Jaouen, A. Nicolaou, W. G. Yin, C. R. Serrao, R. Ramesh, H. Ding, et al., J. Phys.: Condens. Matter 27, 202202 (2015).
  * Monney et al. (2020) C. Monney, M. Herzog, A. Pulkkinen, Y. Huang, J. Pelliciari, P. Olalde-Velasco, N. Katayama, M. Nohara, H. Takagi, T. Schmitt, et al., Phys. Rev. B 102, 085148 (2020).
  * Lee et al. (2013) W. S. Lee, S. Johnston, B. Moritz, J. Lee, M. Yi, K. J. Zhou, T. Schmitt, L. Patthey, V. Strocov, K. Kudo, et al., Phys. Rev. Lett. 110, 265502 (2013).
  * Lu et al. (2018) X. Lu, P. Olalde-Velasco, Y. Huang, V. Bisogni, J. Pelliciari, S. Fatale, M. Dantz, J. G. Vale, E. C. Hunter, J. Chang, et al., Phys. Rev. B 97, 041102(R) (2018).
  * Kim et al. (2023) J.-K. Kim, C. Dietl, H.-W. J. Kim, S.-H. Ha, J. Kim, A. H. Said, J. Kim, and B. J. Kim, J. Synchrotron Rad. 30, 643 (2023).
  * Paris et al. (2020) E. Paris, Y. Tseng, E. M. Pärschke, W. Zhang, M. H. Upton, A. Efimenko, K. Rolfs, , D. E. McNally, L. Maurel, et al., Proceed. Nat. Acad. Sciences USA 117, 24764 (2020).
  * Ishii et al. (2011) K. Ishii, I. Jarrige, M. Yoshida, K. Ikeuchi, J. Mizuki, K. Ohashi, T. Takayama, J. Matsuno, and H. Takagi, Phys. Rev. B 83, 115121 (2011).
  * Lupascu et al. (2014) A. Lupascu, J. P. Clancy, H. Gretarsson, Z. Nie, J. Nichols, J. Terzic, G. Cao, S. S. A. Seo, Z. Islam, and M. H. U. et al., Phys. Rev. Lett. 112, 147201 (2014).
  * Bertinshaw et al. (2020) J. Bertinshaw, J. K. Kim, J. Porras, K. Ueda, N. H. Sung, A. Efimenko, A. Bombardi, J. Kim, B. Keimer, and B. J. Kim, Phys. Rev. B 101, 094428 (2020).
  * Clancy et al. (2023) J. P. Clancy, H. Gretarsson, A. Lupascu, J. A. Sears, Z. Nie, M. H. Upton, J. Kim, Z. Islam, M. Uchida, D. G. Schlom, et al., Phys. Rev. B 107, 054423 (2023).
  * Guo et al. (1994) G. Y. Guo, H. Ebert, W. M. Temmerman, and P. J. Durham, Phys. Rev. B 50, 3861 (1994).
  * Antonov et al. (2004) V. Antonov, B. Harmon, and A. Yaresko, _Electronic Structure and Magneto-Optical Properties of Solids_ (Kluwer, Dordrecht, 2004).
  * Arola et al. (2004) E. Arola, M. Horne, P. Strange, H. Winter, Z. Szotek, and W. M. Temmerman, Phys. Rev. B 70, 235127 (2004).
  * Nemoshkalenko et al. (1983) V. V. Nemoshkalenko, A. E. Krasovskii, V. N. Antonov, V. N. Antonov, U. Fleck, H. Wonn, and P. Ziesche, Phys. status solidi B 120, 283 (1983).
  * Arola et al. (1997a) E. Arola, P. Strange, and B. L. Gyorffy, Phys. Rev. B 55, 472 (1997a).
  * Arola et al. (1997b) E. Arola, P. Strange, and B. L. Gyorffy, Phys. Rev. B 55, 472 (1997b).
  * Antonov et al. (2022a) V. N. Antonov, D. A. Kukusta, and L. V. Bekenov, Phys. Rev. B 105, 155144 (2022a).
  * Shimura et al. (1995) T. Shimura, Y. Inaguma, T. Nakamura, M. Itoh, and Y. Morii, Phys. Rev. B 52, 9143 (1995).
  * Antonov et al. (2006) V. N. Antonov, O. Jepsen, A. N. Yaresko, and A. P. Shpak, J. Appl. Phys. 100, 043711 (2006).
  * Antonov et al. (2007) V. N. Antonov, B. N. Harmon, A. N. Yaresko, and A. P. Shpak, Phys. Rev. B 75, 184422 (2007).
  * Antonov et al. (2010) V. N. Antonov, A. N. Yaresko, and O. Jepsen, Phys. Rev. B 81, 075209 (2010).
  * Andersen (1975) O. K. Andersen, Phys. Rev. B 12, 3060 (1975).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl et al. (1994) P. E. Blöchl, O. Jepsen, and O. K. Andersen, Phys. Rev. B 49, 16223 (1994).
  * Yaresko et al. (2003) A. N. Yaresko, V. N. Antonov, and P. Fulde, Phys. Rev. B 67, 155103 (2003).
  * Dederichs et al. (1984) P. H. Dederichs, S. Blügel, R. Zeller, and H. Akai, Phys. Rev. Lett. 53, 2512 (1984).
  * Pickett et al. (1998) W. E. Pickett, S. C. Erwin, and E. C. Ethridge, Phys. Rev. B 58, 1201 (1998).
  * Antonov et al. (2022b) V. N. Antonov, D. A. Kukusta, and L. V. Bekenov, Phys. Rev. B 105, 155145 (2022b).
  * Pröpper et al. (2016) D. Pröpper, A. N. Yaresko, M. Höppner, Y. Matiks, Y.-L. Mathis, T. Takayama, A. Matsumoto, H. Takagi, B. Keimer, and A. V. Boris, Phys. Rev. B 94, 035158 (2016).
  * Zhang et al. (2013) H. Zhang, K. Haule, and D. Vanderbilt, Phys. Rev. Lett. 111, 246402 (2013).
  * Pröpper et al. (2013) D. Pröpper, A. N. Yaresko, M. Höppner, Y. Matiks, Y.-L. Mathis, T. Takayama, A. Matsumoto, H. Takagi, B. Keimer, A. V. Boris, et al., Phys. Rev. B 88, 085125 (2013).
  * Wang et al. (2013) Q. Wang, Y. Cao, J. A. Waugh, S. R. Park, T. F. Qi, O. B. Korneta, G. Cao, and D. S. Dessau, Phys. Rev. B 87, 245109 (2013).
  * Kim et al. (2012b) B. H. Kim, G. Khaliullin, and B. I. Min, Phys. Rev. Lett. 109, 167205 (2012b).
  * Haskel et al. (2012) D. Haskel, G. Fabbris, M. Zhernenkov, P. P. Kong, C. Q. Jin, G. Cao, and M. van Veenendaal, Phys. Rev. Lett. 109, 027204 (2012).
  * van der Laan and Thole (1988) G. van der Laan and B. T. Thole, Phys. Rev. Lett. 60, 1977 (1988).
  * Liu et al. (2008) G.-Q. Liu, V. N. Antonov, O. Jepsen, and O. K. Andersen, Phys. Rev. Lett. 101, 026408 (2008).
  * Krajewska et al. (2020) A. Krajewska, T. Takayama, R. Dinnebier, A. Yaresko, K. Ishii, M. Isobe, and H. Takagi, Phys. Rev. B 101, 121101(R) (2020).
  * Ishii et al. (2005) K. Ishii, K. Tsutsui, Y. Endoh, T. Tohyama, S. Maekawa, M. Hoesch, K. Kuzushita, M. Tsubota, T. Inami, J. Mizuki, et al., Phys. Rev. Lett. 94, 207003 (2005).
  * Kaneko et al. (2013) T. Kaneko, T. Toriyama, T. Konishi, and Y. Ohta, Phys. Rev. B 87, 035121 (2013).
  * King et al. (2013) P. D. C. King, T. Takayama, A. Tamai, E. Rozbicki, S. M. Walker, M. Shi, L. Patthey, R. G. Moore, D. Lu, K. M. Shen, et al., Phys. Rev. B 87, 241106(R) (2013).


# Highly Efficient Structural Learning of Sparse Staged Trees

Manuele Leonelli<EMAIL_ADDRESS>
School of Science and Technology, IE University, Madrid, Spain  Gherardo
Varando<EMAIL_ADDRESS>
Image Processing Laboratory, Universitat de València, València, Spain

###### Abstract

Several structural learning algorithms for staged tree models, an asymmetric
extension of Bayesian networks, have been defined. However, they do not scale
efficiently as the number of variables considered increases. Here we introduce
the first scalable structural learning algorithm for staged trees, which
searches over a space of models where only a small number of dependencies can
be imposed. A simulation study as well as a real-world application illustrate
our routines and the practical use of such data-learned staged trees.

Keywords: Asymmetric conditional independence; Bayesian networks;
Probabilistic graphical models; Staged trees; Structural learning.

##  1 Introduction

Probabilistic graphical models, and in particular Bayesian networks (BNs), are
nowadays widely used in machine learning to conveniently represent the
relationships existing between the components of a random vector. The directed
acyclic graph (DAG) associated to a BN represents graphically (symmetric)
conditional independence statements, which can be assessed using the
d-separation criterion (Pearl, 1988). Although the underlying DAG can be
expert-elicited, this is often learned from data using algorithms that explore
the space of all possible DAGs (see e.g. Scutari et al., 2019).

For quite some time it has been noticed that the strict assumption of
symmetric conditional independence may be too restrictive to fully represent
the relationship between variables in a dataset (Boutilier et al., 1996;
Chickering et al., 1997; Friedman and Goldszmidt, 1996). However, the
development and use in practice of probabilistic graphical models embedding
asymmetric conditional independence has been limited (see Hyttinen et al.,
2018; Nicolussi and Cazzaro, 2021; Talvitie et al., 2019, for some recent
proposals). Possible reasons behind the limited use of such models could be:
(i) the lack of widely available software; (ii) the complexity of the learning
routines; and (iii) the less intuitive visualization of the associated
independences which are not explicitly represented by a single graph.

Staged trees (Collazo et al., 2018; Smith and Anderson, 2008) are
probabilistic graphical models which, starting from an event tree, represent
any type of asymmetric conditional independence by a partitioning/coloring of
the vertices of the tree. All the model information and the conditional
independences can be read directly from the tree and, in particular, from the
coloring of the vertices. The stagedtrees R package (Carli et al., 2022)
provides a user-friendly implementation of a wide array of structural learning
and inferential routines to fit staged trees to data. Therefore, two of the
main limitations to the use of asymmetric probabilistic graphical models do
not apply to staged trees. Furthermore, a wide toolkit of procedures to work
with staged trees have been developed, including handling missing data,
sensitivity analysis and exploration of equivalence classes, among others (see
Collazo et al., 2018, for details).

On the other hand, learning staged trees from data is complex. Although
efficient structural learning algorithms for such models have been implemented
(e.g. Freeman and Smith, 2011; Leonelli and Varando, 2021; Silander and Leong,
2013), they can only work with a limited number of variables. The main reason
behind this is the explosion of the size of the model search space as the
number of variables considered increases. As an illustration, the number of
DAGs over 6 binary variables is 3781503, whilst there are $1.20019\times
10^{44}$ staged trees under the same conditions (Duarte and Solus, 2021).
Furthermore, the number of DAGs remains constant if variables have more than
two levels, whilst the number of staged trees would further increase
dramatically.

Recent proposals for efficient structural learning of staged trees look at
sub-classes of staged trees, with the aim of reducing the size of the model
space. Carli et al. (2020) defined naive staged trees which have the same
complexity of a naive BN over the same variables. Leonelli and Varando (2022)
considered simple staged trees which have a constrained type of partitioning
of the vertices. Duarte and Solus (2021) defined CStrees which only embed
symmetric and context-specific types of independence, and not others (Pensar
et al., 2016).

One of the first solutions to make structural learning of BNs scalable was to
limit the number of parents each variable can have (Friedman et al., 1999;
Tsamardinos et al., 2006). This was imposed not only to restrict the model
space of possible DAGs, but it also made sense from an applied point of view
since most often only a limited number of variables can be expected to have a
direct influence on another. The option of setting a maximum number of parents
is also available in the standard bnlearn software (Scutari, 2010).

Here, we define a sub-class of staged trees embedding the same idea of
limiting the number of variables that can have a direct influence to another.
As we formalize below, this means that the BN representation associated to
such a staged tree is sparse, meaning that it has a small number of edges. A
structural learning algorithm for this class of staged trees is introduced and
its features are explored in an extensive simulation study.

##  2 Bayesian Networks and Staged Trees

Before introducing staged trees, we give a formal definition of BNs. We then
describe their relationships with staged tree models.

###  2.1 Bayesian Networks

Let $G=([p],E)$ be a DAG with vertex set $[p]=\\{1,\dots,p\\}$ and edge set
$E$. Let $\bm{X}=(X_{i})_{i\in[p]}$ be categorical random variables with joint
mass function $P$ and sample space
$\mathbb{X}=\times_{i\in[p]}\mathbb{X}_{i}$. For $A\subset[p]$, we let
$\bm{X}_{A}=(X_{i})_{i\in A}$ and $\bm{x}_{A}=(x_{i})_{i\in A}$ where
$\bm{x}_{A}\in\mathbb{X}_{A}=\times_{i\in A}\mathbb{X}_{i}$. We say that $P$
is Markov to $G$ if, for $\bm{x}\in\mathbb{X}$,

$P(\bm{x})=\prod_{k\in[p]}P(x_{k}\mid\bm{x}_{\Pi_{k}}),$

where $\Pi_{k}$ is the parent set of $k$ in $G$ and
$P(x_{k}|\bm{x}_{\Pi_{k}})$ is a shorthand for
$P(X_{k}=x_{k}|\bm{X}_{\Pi_{k}}=\bm{x}_{\Pi_{k}})$. It is customary to label
the vertices of a BN so to respect the topological order of $G$ and we
henceforth assume that $1,2,\dots,p$ is a topological order of $G$.

The ordered Markov condition implies conditional independences of the form

$X_{i}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\bm{X}_{[i-1]}\,|\,\bm{X}_{\Pi_{i}}.$

Henceforth, $P$ is assumed to be strictly positive. Let $G$ be a DAG and $P$
Markov to $G$. The _Bayesian network_ model (associated to $G$) is

$\mathcal{M}_{G}=\\{P\in\Delta^{\circ}_{|\mathbb{X}|-1}\,|\,P\mbox{ is Markov
to }G\\}.$

where $\Delta^{\circ}_{|\mathbb{X}|-1}$ is the ($|\mathbb{X}|-1$)-dimensional
open probability simplex.

###  2.2 Staged Trees

Differently to BNs, whose graphical representation is a DAG, staged trees
visualize conditional independence by means of a colored tree. Let $(V,E)$ be
a directed, finite, rooted tree with vertex set $V$, root node $v_{0}$ and
edge set $E$. For each $v\in V$, let $E(v)=\\{(v,w)\in E\\}$ be the set of
edges emanating from $v$ and $\mathcal{C}$ be a set of labels.

An $\bf X$-compatible staged tree is a triple $T=(V,E,\theta)$, where $(V,E)$
is a rooted directed tree and:

  1. 1.

$V={v_{0}}\cup\bigcup_{i\in[p]}\mathbb{X}_{[i]}$;

  2. 2.

For all $v,w\in V$, $(v,w)\in E$ if and only if
$w=\bm{x}_{[i]}\in\mathbb{X}_{[i]}$ and $v=\bm{x}_{[i-1]}$, or $v=v_{0}$ and
$w=x_{1}$ for some $x_{1}\in\mathbb{X}_{1}$;

  3. 3.

$\theta:E\rightarrow\mathcal{L}=\mathcal{C}\times\cup_{i\in[p]}\mathbb{X}_{i}$
is a labelling of the edges such that
$\theta(v,\bm{x}_{[i]})=(\kappa(v),x_{i})$ for some function
$\kappa:V\to\mathcal{C}$. The function $k$ is called the colouring of the
staged tree $T$.

If $\theta(E(v))=\theta(E(w))$ then $v$ and $w$ are said to be in the same
_stage_. Therefore, the equivalence classes induced by $\theta(E(v))$ form a
partition of the internal vertices of the tree in _stages_.

Points 1 and 2 above construct a rooted tree where each root-to-leaf path, or
equivalently each leaf, is associated to an element of the sample space
$\mathbb{X}$. Then a labeling of the edges of such a tree is defined where
labels are pairs with one element from a set $\mathcal{C}$ and the other from
the sample space $\mathbb{X}_{i}$ of the corresponding variable $X_{i}$ in the
tree. By construction, $\bf X$-compatible staged trees are such that two
vertices can be in the same stage if and only if they correspond to the same
sample space.

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$$v_{9}$$v_{10}$$v_{11}$$v_{12}$$v_{13}$$v_{14}$$X_{1}=0$$X_{1}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$

Figure 1: Example of an $\bm{X}$-compatible staged tree over four binary
random variables.

Figure 1 reports an $(X_{1},X_{2},X_{3},X_{4})$-compatible staged tree over
four binary variables. The coloring given by the function $\kappa$ is shown in
the vertices and each edge $(\cdot,(x_{1},\ldots,x_{i}))$ is labeled with
$X_{i}=x_{i}$. The edge labeling $\theta$ can be read from the graph combining
the text label and the color of the emanating vertex. The staging of the
staged tree in Figure 1 is given by the partition $\\{v_{0}\\}$,
$\\{v_{1}\\}$, $\\{v_{2}\\}$, $\\{v_{3},v_{4}\\}$, $\\{v_{5},v_{6}\\}$,
$\\{v_{7},v_{11}\\}$, $\\{v_{8},v_{12}\\}$, $\\{v_{9},v_{13}\\}$ and
$\\{v_{10},v_{14}\\}$.

The parameter space associated to an $\bf X$-compatible staged tree
$T=(V,E,\theta)$ with labeling $\theta:E\rightarrow\mathcal{L}$ is defined as

$\Theta_{T}=\Big{\\{}\bm{y}\in\mathbb{R}^{|\theta(E)|}\;|\;\forall~{}e\in
E,y_{\theta(e)}\in(0,1)\textnormal{ and }\sum_{e\in
E(v)}y_{\theta(e)}=1\Big{\\}}.$ (1)

Equation (1) defines a class of probability mass functions over the edges
emanating from any internal vertex coinciding with conditional distributions
$P(x_{i}|\bm{x}_{[i-1]})$, $\bm{x}\in\mathbb{X}$ and $i\in[p]$. In the staged
tree in Figure 1 the staging $\\{v_{3},v_{4}\\}$ implies that the conditional
distribution of $X_{3}$ given $X_{1}=0$ and $X_{2}=0$, represented by the
edges emanating from $v_{3}$, is equal to the conditional distribution of
$X_{3}$ given $X_{1}=0$ and $X_{2}=1$. A similar interpretation holds for the
staging $\\{v_{5},v_{6}\\}$. This in turn implies that
$X_{3}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{2}|X_{1}$, thus illustrating that the
staging of a tree is associated to conditional independence statements.

Let $\bm{l}_{T}$ denote the leaves of a staged tree $T$. Given a vertex $v\in
V$, there is a unique path in $T$ from the root $v_{0}$ to $v$, denoted as
$\lambda(v)$. The _depth_ of a vertex $v\in V$ equals the number of edges in
$\lambda(v)$. For any path $\lambda$ in $T$, let $E(\lambda)=\\{e\in
E:e\in\lambda\\}$ denote the set of edges in the path $\lambda$.

The _staged tree model_ $\mathcal{M}_{T}$ associated to the $\bf X$-compatible
staged tree $(V,E,\theta)$ is the image of the map

$\begin{array}[]{llll}\phi_{T}&:&\Theta_{T}&\to\Delta_{|\bm{l}_{T}|-1}^{\circ}\\\
&&\bm{y}&\mapsto\Big{(}\prod_{e\in
E(\lambda(l))}y_{\theta(e)}\Big{)}_{l\in\bm{l}_{T}}\end{array}$ (2)

Therefore, staged trees models are such that atomic probabilities are equal to
the product of the edge labels in root-to-leaf paths and coincide with the
usual factorization of mass functions via recursive conditioning.

###  2.3 Staged Trees and Bayesian Networks

Although the relationship between BNs and staged trees was already formalized
by Smith and Anderson (2008), a formal procedure to represent a BN as a staged
tree has been only recently introduced in Duarte and Solus (2020) and Varando
et al. (2021).

Assume $\bm{X}$ is topologically ordered with respect to a DAG $G$ and
consider an $\bf X$-compatible staged tree with vertex set $V$, edge set $E$
and labeling $\theta$ defined via the coloring
$\kappa(\bm{x}_{[i]})=\bm{x}_{\Pi_{i}}$ of the vertices. The staged tree
$T_{G}$, with vertex set $V$, edge set $E$ and labeling $\theta$ so
constructed, is called _the staged tree model of $G$_. Importantly,
$\mathcal{M}_{G}=\mathcal{M}_{T_{G}}$, i.e. the two models are exactly the
same, since they entail exactly the same factorization of the joint
probability. Clearly, the staging of $T_{G}$ represents the Markov conditions
associated to the graph $G$.

Varando et al. (2021) approached the reverse problem of transforming a staged
tree into a BN. Of course, since staged trees represent more general
asymmetric conditional independences, given a staged tree $T$ most often there
is no BN with DAG $G_{T}$ such that $\mathcal{M}_{T}=\mathcal{M}_{G_{T}}$.
However, Varando et al. (2021) introduced an algorithm that, given an
$\bm{X}$-compatible staged tree $T$, finds the minimal DAG $G_{T}$ such that
$\mathcal{M}_{T}\subseteq\mathcal{M}_{G_{T}}$. Minimal means that such a DAG
$G_{T}$ embeds all symmetric conditional independences that are in
$\mathcal{M}_{T}$ and that there are no DAGs with less edges than $G_{T}$
embedding the same conditional independences.

$1$$2$$3$$4$ Figure 2: The DAG $G$ such that $T_{G}$ is the staged tree in
Figure 1.

As an illustration, the staged tree in Figure 1 can be constructed as the
$T_{G}$ from the BN with DAG in Figure 2, embedding the conditional
independences $X_{3}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{2}\;|\;X_{1}$ and
$X_{4}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{1}\;|\;X_{2},X_{3}$. Conversely, consider
the staged tree $T$ in Figure 3, which differs from the one in Figure 1 only
on a different coloring of the vertices $v_{6}$ and $v_{14}$. Such a staged
tree does not embed any symmetric conditional independence, only non-symmetric
ones, and therefore there is no DAG $G_{T}$ such that
$\mathcal{M}_{G_{T}}=\mathcal{M}_{T}$. Furthermore, the minimal DAG $G_{T}$
such that $\mathcal{M}_{T}\subseteq\mathcal{M}_{G_{T}}$ is the complete one
since the staging of the tree implies direct dependence between every pair of
variables.

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$$v_{9}$$v_{10}$$v_{11}$$v_{12}$$v_{13}$$v_{14}$$X_{1}=0$$X_{1}=1$$X_{2}=0$$X_{2}=1$$X_{2}=0$$X_{2}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{3}=0$$X_{3}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$$X_{4}=0$$X_{4}=1$

Figure 3: Example of an $\bm{X}$-compatible staged tree over four binary
random variables.

###  2.4 Structural Learning Algorithms Inducing Sparsity

Since the model search space of staged trees is huge, we consider here
structural learning for a subclass of staged trees that we define next.

###### Definition 1

A staged tree $T$ is in the class of $k$-parents staged trees if the maximum
in-degree in $G_{T}$ is less or equal to $k$.

For instance, the staged tree in Figure 1 is in the class of 2-parents staged
trees, whilst the one in Figure 3 is not, since its associated minimal DAG is
such that $X_{4}$ has three parents.

Of course, the class of $k$-parents staged trees is much smaller than the one
of $\bm{X}$-compatible ones, for small values of $k$, and therefore structural
learning is expected to be quicker. Here we define a structural learning
algorithm to learn a staged tree in the class of $k$-parents staged trees
which consists of the following steps: (i) learn a BN with DAG $G$ having at
most $k$ parents (for instance using bnlearn); (ii) construct the equivalent
staged tree $T_{G}$; (iii) run the backward hill-climbing algorithm of Carli
et al. (2022) which only joins stages together (no splitting of stages) based
on the minimization of the model BIC (Görgen et al., 2022). Call the resulting
staged tree $T$. It can be easily proven that $G_{T}$ has at most $k$ parents
and $\mathcal{M}_{G_{T}}\subseteq\mathcal{M}_{G}$.

Although the idea of using the staged tree equivalent to a BN as starting
point of a structural learning algorithm (or at least using a partial ordering
associated to such a BN) is not new (see e.g. Barclay et al., 2013), here we
specifically use such a strategy to limit the complexity of the learned staged
tree. This has two major advantages: (i) the speed of the algorithms increases
greatly; (ii) the non-symmetric conditional independences can be easily
visualized even when a large number of variables are present, as illustrated
in Section 4.

##  3 Experiments

We perform simulation experiments to evaluate the proposed learning strategy
for $k$-parents staged trees. Moreover we it compare to standard learning of
staged trees and DAGs.

In all the simulated experiments we generated data from random $k$-parents
staged trees, which are obtained as follow: (1) A random DAG $G$ with fixed
topological order $X_{1},\ldots,X_{p}$ is obtained by randomly selecting up to
$k$ parents uniformly from $\\{X_{1},\ldots,X_{j-1}\\}$ for each $X_{j}$; (2)
The equivalent staged tree $T_{G}$ is obtained; (3) Stages in $T_{G}$ are
randomly merged with probability $0.5$. The obtained staged tree $T$ is such
that $G_{T}$ is a sub-graph of $G$ and thus it is a $k$-parents staged tree;
(4) Lastly, we assign random probabilities (uniformly from the simplex) to
each stage of the staged tree $T$. Once we generate the random $k$-parents
staged tree $T$, we can easily sample observations of $X_{1},\ldots,X_{p}$
from it via sequential sampling (as implemented in the stagedtrees package).
For each fixed parameters ($k$, $p$ and number of observation sampled) we
repeat the experiments $20$ times and report averages and standard errors.

###  3.1 Oracle DAG

Figure 4: Average computation times for learning the stagedtrees (solid lines)
with the standard BHC approach (blue) and the proposed $k$-parents BHC method
(black). Additionally the time spent creating the initial model are shown
(dashed lines). Results for simulated data from $p=3,\ldots,20$ binary
variables, BHC results are obtained only up to $10$ variables. In different
columns results for the maximum number of parents $2,3,4$.

We first consider an ideal scenario where we evaluate the performance of the
proposed method when the starting DAG is the graph $G$ used to generate the
true staged tree $T$. We thus run a standard backward hill-climbing (BHC)
procedure, as implemented in the stagedtrees package, starting from both the
full saturated tree model and starting from the model $T_{G}$. All heuristic
searches optimize the BIC score.

We plot in Figure 4 the average computation time as a function of the system
size $p=3,\ldots,20$ and varying the number of maximum parents $k=2,3,4$. For
each one of the two approaches we split the computation time in: build time,
the time spent building the starting tree; and search time, the time spent
running the search algorithm. We can observe that starting the search from the
DAG-equivalent tree allows us to scale the algorithm easily up to $20$
variables, while the standard approach starting from the full tree become
quickly infeasible after $10$ variables.

As a sanity check we also compute the normalized hamming distance (the sum
across the depth of the tree of the average number of nodes for which coloring
needs to change to obtain the same staged tree) and the context intervention
distance (Leonelli and Varando, 2021) between the true and learned models. As
expected the $k$-parents trees obtain better results, which we do not report
here for lack of space.

###  3.2 Learned DAG

We perform now a simulation study similar to the oracle setting but, as in a
more realistic scenario, we do not assume knowledge of the DAG $G$. We thus,
first estimate a DAG $\hat{G}$ from data using the hill-climbing approach in
the bnlearn package (Scutari, 2010), and then we apply the BHC learning
algorithm starting from the staged tree $T_{\hat{G}}$ (bhcdag). We compare the
obtained tree to $T_{\hat{G}}$ (the tree equivalent to the learned DAG, dag)
and the output of the BHC algorithm starting from the full saturated model
(bhc). We run $20$ replications of the experiment for different system sizes
($p=6,10,20$), for different number of maximum parents in the true DAG
($k=2,3,4$) and with sample sizes ranging from $100$ to $10000$.

Figure 5: Average normalized hamming distance between true and estimated
staged trees as a function of sample size, for different methods, system size
($p$, rows) and maximum parents ($k$, columns).

In Figure 5 we plot the average (across repetitions) normalized hamming
distance between the estimated tree and the true one. We can observe that the
$k$-parents staged trees, obtained by the BHC algorithm starting from
$T_{\hat{G}}$ (bhcdag), are closer to the true data-generating models, with
respect to both the output of BHC starting from the saturated model and the
tree $T_{\hat{G}}$ obtained from the estimated DAG.

##  4 COVID-19 Drivers and Country Risks

We next extend the analysis of Qazi and Simsekler (2022) who developed a BN to
investigate how various country risks and risks associated to the COVID-19
epidemics relate to each other. In particular here we focus on how various
types of risks affect the overall country risk associated to COVID-19.

For this purpose, as in Qazi and Simsekler (2022), the dataset used in the
analysis comes from the combination of two sources. Country-level exposure to
COVID-19 risks are retrieved from INFORM (INFORM, 2022). The data comprises of
a score between zero and ten for 191 countries for three drivers of COVID-19
risks, namely hazard and exposure, vulnerability, and lack of coping capacity.
An overall COVID-19 risk index, again between zero and ten, is constructed
from these three drivers. Country-level exposure to various socioeconomic risk
factors are collected from Euler Hermes (Euler Hermes, 2022). The ratings for
five drivers of country risk, namely economic, political, financing,
commercial and business environment are collected for 188 countries (the
indexes are integer-valued between one and four or six). The combined dataset
comprises 181 countries. Each variable is discretized into two levels using
the clustering method from the arules package (Hahsler et al., 2005).

Figure 6: BN learned for the COVID-19 drivers and country risks.

A BN is learned using the hc function of the bnlearn package with the
constraint that the overall COVID-19 risk must be a leaf of the network and is
reported in Figure 6. Without specifying it, the learned BN is such that each
vertex has at most two parents. The DAG suggests that COVID-19 risk is
conditionally independent of all other drivers given the lack of coping
capacity and hazard & exposure.

Starting from this BN, a staged tree in the class of 2-parents staged trees is
learned using the algorithm of Section 2.4. This staged tree provides a better
representation of the data since it has a BIC of 1521.965, compared to the BIC
of 1547.634 of the BN. The learned staged tree embeds the same set of
symmetric conditional independences as in the BN of Figure 2, but also non-
symmetric ones. Of course, the full tree cannot be easily visualized since,
for instance, there are $2^{9}=512$ vertices with depth nine. However, since
it is known that COVID-19 risk only depends on the lack of coping capacity and
hazard & exposure we can construct the staged tree over these three variables
only and easily visualize further non-symmetric dependences. This is reported
in Figure 7, which shows the presence of the context-specific independence
between COVID-19 risk and hazard & exposure for lack of coping capacity equal
to low. Similar interpretations could be drawn by constructing the “partial”
staged trees associated to other variables.

Of course a generic staged tree would provide a better representation of the
data. For instance, one learned with the backward hill-climbing of Carli et
al. (2022) starting from the saturated model has a BIC of 1453.587. However,
its complete visualization is again unfeasible and plots as the one of Figure
7 are in general not viable since there are no constraints on the number of
parents of $G_{T}$. Indeed, whilst the DAG $G_{T}$ for the tree in Figure 7
has 13 edges (as the DAG in Figure 6), the DAG $G_{T}$ from the generic staged
tree is complete and consisting of 36 edges, meaning that all variables are
directly related to one another.

$v_{0}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$H = highH = lowC = highC =
lowC = highC = lowR = highR = lowR = highR = lowR = highR = lowR = highR = low

Figure 7: 2-parents staged tree for the COVID-19 risk constructed over hazard
& exposure (H), lack of coping capacity (C) and COVID-19 risk (R).

##  5 Discussion

We defined a novel sub-class of staged tree models borrowing the idea of
limiting the number of parents in BNs. A structural learning algorithm for
such a class has been introduced and its properties illustrated both in
simulation experiments and in a real-world application.

The number of parents of a variable in the staged tree is limited by those of
the learned BN. For instance, in our application, although no limit on the
number of parents was set, there were at most two parents. However, non-
symmetric dependences may have been missed by the BN learning algorithm which
is specifically designed to account for symmetric ones. One possibility could
be to add edges to the learned BN, for instance between variables such that
their conditional mutual information is large, and then run a backward hill-
climbing algorithm of Carli et al. (2022) over that BN. The feasibility of
such an algorithm is the focus of current research.

### Acknowledgements

Gherardo Varando’s work was funded by the European Research Council (ERC)
Synergy Grant “Understanding and Modelling the Earth System with Machine
Learning (USMILE)” under Grant Agreement No 855187.

## References

  * Barclay et al. (2013) L. M. Barclay, J. L. Hutton, and J. Q. Smith.  Refining a Bayesian network using a chain event graph.  _International Journal of Approximate Reasoning_ , 54:1300–1309, 2013.
  * Boutilier et al. (1996) C. Boutilier, N. Friedman, M. Goldszmidt, and D. Koller.  Context-specific independence in Bayesian networks.  In _Proceedings of the Twelfth Conference on Uncertainty in Artificial Intelligence_ , pages 115–123, 1996.
  * Carli et al. (2020) F. Carli, M. Leonelli, and G. Varando.  A new class of generative classifiers based on staged tree models.  _arXiv:2012.13798_ , 2020.
  * Carli et al. (2022) F. Carli, M. Leonelli, E. Riccomagno, and G. Varando.  The R package stagedtrees for structural learning of stratified staged trees.  _Journal of Statistical Software_ , 102(6):1–30, 2022.
  * Chickering et al. (1997) D. M. Chickering, D. Heckerman, and C. Meek.  A Bayesian approach to learning Bayesian networks with local structure.  In _Proceedings of 13th Conference on Uncertainty in Artificial Intelligence_ , pages 80–89, 1997.
  * Collazo et al. (2018) R. Collazo, C. Görgen, and J. Smith.  _Chain event graphs_.  Chapmann & Hall, 2018.
  * Duarte and Solus (2020) E. Duarte and L. Solus.  Algebraic geometry of discrete interventional models.  _arXiv:2012.03593_ , 2020.
  * Duarte and Solus (2021) E. Duarte and L. Solus.  Representation of context-specific causal models with observational and interventional data.  _arXiv:2101.09271_ , 2021.
  * Euler Hermes (2022) Euler Hermes.  Country risk reports.  Retrieved from: https://www.eulerherm es.com/en_global/economic-research/country-reports.html, 2022.
  * Freeman and Smith (2011) G. Freeman and J. Q. Smith.  Bayesian MAP model selection of chain event graphs.  _Journal of Multivariate Analysis_ , 102(7):1152–1165, 2011.
  * Friedman and Goldszmidt (1996) N. Friedman and M. Goldszmidt.  Learning Bayesian networks with local structure.  In _Proceedings of the 12th Conference on Uncertainty in Artificial Intelligence_ , pages 252–262, 1996.
  * Friedman et al. (1999) N. Friedman, I. Nachman, and D. Pe’er.  Learning Bayesian network structure from massive datasets: The “sparse candidate” algorithm.  In _Proceedings of the 15th Conference on Uncertainty in Artificial Intelligence_ , pages 206–215, 1999.
  * Görgen et al. (2022) C. Görgen, M. Leonelli, and O. Marigliano.  The curved exponential family of a staged tree.  _Electronic Journal of Statistics_ , 16(1):2607–2620, 2022.
  * Hahsler et al. (2005) M. Hahsler, B. Gruen, and K. Hornik.  arules – A computational environment for mining association rules and frequent item sets.  _Journal of Statistical Software_ , 14(15):1–25, 2005.
  * Hyttinen et al. (2018) A. Hyttinen, J. Pensar, J. Kontinen, and J. Corander.  Structure learning for Bayesian networks over labeled DAGs.  In _International Conference on Probabilistic Graphical Models_ , pages 133–144, 2018.
  * INFORM (2022) INFORM.  COVID-19 risk index.  Retrieved from: https://drmkc.jrc. ec.europa.eu/inform-index/INFORM-Covid-19, 2022.
  * Leonelli and Varando (2021) M. Leonelli and G. Varando.  Context-specific causal discovery for categorical data using staged trees.  _arXiv:2106.04416_ , 2021.
  * Leonelli and Varando (2022) M. Leonelli and G. Varando.  Structural learning of simple staged trees.  _arXiv:2203.04390_ , 2022.
  * Nicolussi and Cazzaro (2021) F. Nicolussi and M. Cazzaro.  Context-specific independencies in stratified chain regression graphical models.  _Bernoulli_ , 27(3):2091–2116, 2021.
  * Pearl (1988) J. Pearl.  _Probabilistic reasoning in intelligent systems: networks of plausible inference_.  Morgan Kaufmann, 1988.
  * Pensar et al. (2016) J. Pensar, H. Nyman, J. Lintusaari, and J. Corander.  The role of local partial independence in learning of Bayesian networks.  _International Journal of Approximate Reasoning_ , 69:91–105, 2016.
  * Qazi and Simsekler (2022) A. Qazi and M. C. E. Simsekler.  Nexus between drivers of COVID-19 and country risks.  _Socio-Economic Planning Sciences_ , page 101276, 2022.
  * Scutari (2010) M. Scutari.  Learning Bayesian networks with the bnlearn R package.  _Journal of Statistical Software_ , 35(3):1–22, 2010.
  * Scutari et al. (2019) M. Scutari, C. E. Graafland, and J. M. Gutiérrez.  Who learns better Bayesian network structures: Accuracy and speed of structure learning algorithms.  _International Journal of Approximate Reasoning_ , 115:235–253, 2019.
  * Silander and Leong (2013) T. Silander and T. Leong.  A dynamic programming algorithm for learning chain event graphs.  In _Proceedings of the International Conference on Discovery Science_ , pages 201–216, 2013.
  * Smith and Anderson (2008) J. Smith and P. Anderson.  Conditional independence and chain event graphs.  _Artificial Intelligence_ , 172(1):42 – 68, 2008\.
  * Talvitie et al. (2019) T. Talvitie, R. Eggeling, and M. Koivisto.  Learning Bayesian networks with local structure, mixed variables, and exact algorithms.  _International Journal of Approximate Reasoning_ , 115:69–95, 2019.
  * Tsamardinos et al. (2006) I. Tsamardinos, L. E. Brown, and C. F. Aliferis.  The max-min hill-climbing Bayesian network structure learning algorithm.  _Machine Learning_ , 65(1):31–78, 2006.
  * Varando et al. (2021) G. Varando, F. Carli, and M. Leonelli.  Staged trees and asymmetry-labeled DAGs.  _arXiv:2108.01994_ , 2021.