# Surface recombination and out of plane diffusivity of free excitons in
hexagonal boron nitride

Sébastien Roux  Laboratoire d’Etude des Microstructures, ONERA-CNRS,
Université Paris-Saclay, BP 72, 92322 Châtillon Cedex, France  Université
Paris-Saclay, UVSQ, CNRS, GEMaC, 78000, Versailles, France  Christophe Arnold
Université Paris-Saclay, UVSQ, CNRS, GEMaC, 78000, Versailles, France  Etienne
Carré  Laboratoire d’Etude des Microstructures, ONERA-CNRS, Université Paris-
Saclay, BP 72, 92322 Châtillon Cedex, France  Université Paris-Saclay, UVSQ,
CNRS, GEMaC, 78000, Versailles, France  Eli Janzen  Tim Taylor Department of
Chemical Engineering, Kansas State University Manhattan, KS 66506, USA  James
H. Edgar  Tim Taylor Department of Chemical Engineering, Kansas State
University Manhattan, KS 66506, USA  Camille Maestre  Laboratoire des
Multimatériaux et Interfaces, UMR CNRS 5615, Univ Lyon
Université Claude Bernard Lyon 1, F-69622 Villeurbanne, France  Bérangère
Toury  Laboratoire des Multimatériaux et Interfaces, UMR CNRS 5615, Univ Lyon
Université Claude Bernard Lyon 1, F-69622 Villeurbanne, France  Catherine
Journet  Laboratoire des Multimatériaux et Interfaces, UMR CNRS 5615, Univ
Lyon
Université Claude Bernard Lyon 1, F-69622 Villeurbanne, France  Vincent
Garnier  Laboratoire MATEIS, UMR CNRS 5510, Univ Lyon, INSA Lyon, F-69621
Villeurbanne, France  Philippe Steyer  Laboratoire MATEIS, UMR CNRS 5510, Univ
Lyon, INSA Lyon, F-69621 Villeurbanne, France  Takashi Taniguchi  Research
Center for Materials Nanoarchitectonics, National Institute for Materials
Science, 1-1 Namiki, Tsukuba 305-0044, Japan  Kenji Watanabe  Research Center
for Electronic and Optical Materials, National Institute for Materials
Science, 1-1 Namiki, Tsukuba 305-0044, Japan  Annick Loiseau
<EMAIL_ADDRESS>Laboratoire d’Etude des Microstructures, ONERA-CNRS,
Université Paris-Saclay, BP 72, 92322 Châtillon Cedex, France  Julien Barjon
<EMAIL_ADDRESS>Université Paris-Saclay, UVSQ, CNRS, GEMaC, 78000,
Versailles, France

###### Abstract

We present a novel experimental protocol using cathodoluminescence
measurements as a function of the electron incident energy to study both
exciton diffusion in a directional way and surface exciton recombination. Our
approach overcomes the challenges of anisotropic diffusion and the limited
applicability of existing methods to the bulk counterparts of 2D materials.
The protocol is then applied at room and at cryogenic temperatures to four
bulk hexagonal boron nitride crystals grown by different synthesis routes. The
exciton diffusivity depends on the sample quality but not on the temperature,
indicating it is limited by defect scattering even in the best quality
crystals. The lower limit for the diffusivity by phonon scattering is 0.2
cm2.s-1. Diffusion lengths were as much as 570 nm. Finally, the surface
recombination velocity exceeds 105 cm2.s-1, at a level similar to silicon or
diamond. This result reveals that surface recombination could strongly limit
light-emitting devices based on 2D materials.

##  .1 I - Introduction

In most 2D semiconductors, a surprising increase in the luminescence
efficiency appears when reducing the thickness at the atomic level, especially
between bilayer and monolayer crystals where the bandgap undergoes its
indirect to direct transition [1, 2, 3]. However, in the specific case of
hexagonal boron nitride (hBN), the bulk crystal exhibits a high radiative
efficiency up to 50% in cathodoluminescence (CL) [4], with a short radiative
lifetime of 30 ns [5], which arises from the high compacity of its indirect
excitons. Conversely, the luminescence intensity for hBN thin films [6] is
dramatically decreased : there is no CL signal for less than 6 layers. The
luminescence of the monolayer hBN free exciton has only recently been
observed, thanks to extremely long photoluminescence experiments [7, 8]. The
underlying reasons for the decrease in intensity with the thickness of the hBN
layers are yet to be fully understood and elucidated. It will be explored in
this work in terms of surface effect and exciton diffusion.

Often ignored in 2D materials, surfaces play a well-known role in 3D
semiconductors in limiting light emission efficiency in standard
optoelectronic devices. The rate of surface recombination of free charge
carriers is governed by the physical properties of the defect states at the
surface, so that it could be reduced by using passivation strategies [9, 10,
11]. It is quantitatively characterised by a parameter called the surface
recombination velocity, $s$, which typically ranges between 1 and 106 cm.s-1
depending on the material and its surface termination. [9, 12, 11, 13]. In
hBN, due to the sp2 orbital hybridisation, all chemical bonds are satisfied
between atoms on the same atomic plane, so there are no dangling bonds at the
surface. Considering this, the surface recombination rate for free carriers
should be low in hBN compared to sp3 hybridised semiconductors with dandling
bounds [9]. However, strong surface effects are suggested by the low
luminescence efficiency of hBN crystals of a few atomic layers [6]. In hBN,
and more generally in 2D materials, surface recombinations are still
unexplored and remain an open question.

Surface recombination and the associated luminescence losses are promoted by
the diffusion of charge carriers injected into the volume of a crystal towards
its surfaces. In the context of lamellar materials, such as the 2Ds, the out-
of-plane diffusivity is the key parameter for understanding the surface
recombination properties. Furthermore, understanding exciton diffusion in 2D
semiconductors is of strategic importance for the development of the new
”excitronic” devices based on the control of exciton flux, which has attracted
attention with the advent of 2D heterostructures [14, 15, 16, 17]. However,
experiments and devices have only focused on the in-plane exciton diffusion in
2D monolayers of transition metal dichalcogenides [18, 19, 20].

The study of exciton diffusion can be carried out by PL experiments using
various pump-probe techniques [21, 22], or by CL experiments, which benefit
from a better control of the excitation volume, down to tens of $nm^{3}$, by
analysing the transfer of excitons from a localised excitation to a physical
probe [23, 24, 25, 26, 27, 24, 28]. To properly analyse such experiments, one
must decouple the combined effects of (i) in-plane diffusion, (ii) out-of-
plane diffusion, (iii) surface recombinations, and (iv) inhomogeneous in-depth
energy deposition. Due to the anisotropy and the still unexplored surface
recombination properties of 2D materials, the use of traditional protocols for
their out-of-plane exciton diffusion poses significant challenges.

This work presents a novel experimental protocol in CL that provides a
powerful tool for investigating the out-of-plane diffusion and the surface
recombination properties of a bulk crystal using a 1D excitation geometry. The
method is applicable to all types of semiconductors, including those with
anisotropic diffusion, and does not require prior knowledge of the surface
recombination properties. After introducing the principle of the experiment, a
theoretical analysis of diffusion is given, providing an analytical
description of experimental results. The method is then applied to four bulk
hBN samples at room and cryogenic temperature, allowing the evaluation of out-
of-plane diffusion length, diffusion constant, surface recombination velocity,
and their dependence on the crystal quality and temperature.

##  .2 II - Analytical model for a 1D diffusion experiment in CL

We first present the principle of the CL experiment. As illustrated in Figure
1, a 1D diffusion geometry is obtained by defocusing significantly the
incident electron beam. The CL intensity is measured as a function of the
excitation depth, which is modulated by varying the acceleration voltage of
the incident electrons. At low excitation depths, surface recombination
effects lead to a reduction of the CL intensity and the key point is to use
the recombinant surface as a probe for exciton diffusion.

Figure 1: Scheme of the CL experiment in a 1D geometry to analyse the out-of-
plane exciton diffusion in the bulk crystal and the recombination properties
of the surface.

The quantitative analysis of such an experiment requires the consideration of
two essential aspects: (i) the development of a simple and quantitative
description of the excitation depth as a function of the acceleration voltage,
and (ii) the establishment of an analytical model that explains the decrease
in CL intensity at low excitation depths, taking into account the surface
effect and in relation to the relevant physical parameters. We first address
theoretically these two issues.

The free exciton 111Assuming strongly bounded excitons, i.e. negligible free
carrier concentrations. concentration $n(z)$ that is established under steady-
state generation $\mathcal{G}$ in a semiconductor crystal results from the
diffusion equation studied here in a 1D symmetry :

$D\frac{\partial^{2}n}{\partial z^{2}}-\frac{n}{\tau}+\mathcal{G}=0$ (1)

$D$ is the exciton diffusivity along the $z$ axis, and its determination is
one of the goals of the CL experiments. $\tau$ is the exciton lifetime in the
bulk of the crystal, i.e. far from the surface. It is measured by time-
resolved CL (TRCL) experiments at high voltage as detailed in a previous work
[5].

We consider here a semi-infinite crystal surrounded by vacuum for $z<0$. The
surface recombinations at $z=0$ are described as a boundary condition of the
diffusion equation, with the surface recombination velocity $s$ :

$D\left[\frac{\partial n}{\partial z}\right]_{z=0}=s\>n(0)$ (2)

Dealing with the generation $\mathcal{G}$ due to an excitation spread in depth
is a rather difficult mathematical problem. The solution for a surface source
located in a plane at the depth $z=z^{\prime}$ was formulated in 1955 by van
Roosbroeck [30], where the depth distribution of particles in the crystal is
defined separately on each side of the generation plane:

$f_{z^{\prime},a}(z)=\left\\{\begin{array}[]{ll}\frac{1}{2L}e^{-\frac{z^{\prime}}{L}}\left(e^{\frac{z}{L}}+ae^{-\frac{z}{L}}\right)\text{
if }z\leqslant z^{\prime}\\\
\frac{1}{2L}e^{-\frac{z}{L}}\left(e^{\frac{z^{\prime}}{L}}+ae^{-\frac{z^{\prime}}{L}}\right)\text{
if }z\geqslant z^{\prime}\end{array}\right.$ (3)

where $L=\sqrt{D\tau}$ is the diffusion length along the $z$ axis. An
interpretation of this result with virtual image sources is as follows: the
recombinant surface at $z=0$ might be viewed as a semi-reflecting mirror for
the real sources in the crystal with a reflection coefficient
$a=\frac{1-S}{1+S}$. In this expression $S=\frac{s}{v_{D}}$ is the reduced
recombination velocity, i.e. normalised by the diffusion velocity
$v_{D}=\frac{L}{\tau}=\sqrt{\frac{D}{\tau}}$. The limit cases are:
$S$=+$\infty$ ($a$=-1) for a surface with infinite recombination velocity and
$S$=0 ($a$=1) assimilated to a total reflection.

Since the diffusion equation is linear, an arbitrary in-depth profile of
exciton generation $g(z)$ normalised such as $\int_{0}^{+\infty}g(z)dz=1$
results in an exciton concentration distribution that is written here
222$g(z)=\delta(z-z^{\prime})$ for a plane source in $z=z^{\prime}$:

$n(z)=\frac{G\tau}{\mathcal{A}}\int_{0}^{+\infty}f_{z^{\prime},a}(z)g(z^{\prime})dz^{\prime}$
(4)

This expression indicates how excitons are distributed in the semiconductor
crystal by the combined effects of diffusion from a given excitation profile
and recombinations at the crystal surface. The pre-factor contains the
experimental parameters of a CL experiment. $\mathcal{A}$ is the area of the
electron beam impinging on a 3D crystal. Experimentally one should ensure that
$\mathcal{A}\gg L^{2}$ to assume a 1D diffusion. For instance, the electron
beam spot should be sufficiently defocused to be much larger than the in-plane
diffusion length (see Figure 1). $G$ (s-1) is the electron-hole generation
rate in the semiconductor. $G=\frac{(1-f)Vi}{\langle E_{eh}\rangle}$ is a
simple function of the acceleration voltage, $V$, the beam current, $i$, the
retrodiffusion factor, $f$, and the average electron-hole energy, $\langle
E_{eh}\rangle$ required for the formation of an electron-hole pair. $\langle
E_{eh}\rangle\approx 3E_{g}$ is a commonly used relation as a function of the
bandgap energy $E_{g}$ when $\langle E_{eh}\rangle$ is not known
experimentally for the material of interest [32].

The recombinations at the crystal surface at $z=0$ are responsible for exciton
losses and in fine for a decrease of their luminescence intensity $I$. Instead
of integrating Eq. 4 over the crystal volume to evaluate the total number of
steady-state excitons in the crystal, the quantity of lost excitons
$N^{\prime}$ is calculated more directly [33] considering the diffusion flux
at $z=0$ :

$N^{\prime}=D\left[\frac{\partial n}{\partial z}\right]_{z=0}\mathcal{A}\tau$
(5)

which gives the following result after calculation:

$\frac{N^{\prime}}{N}=\frac{(1-a)}{2}\int_{0}^{+\infty}e^{-\frac{z^{\prime}}{L}}g(z^{\prime})dz^{\prime}$
(6)

where $N=G\tau$ is the steady-state exciton population in a crystal without
surface losses ($s=0$). The luminescence intensity $I$ is proportional to the
total number of excitons present in the crystal under steady-state excitation,
via the constant radiative rate of excitons. Noting $I_{0}$ the luminescence
intensity when surface recombination is negligible, measured experimentally in
the high voltages limit, we obtain $\frac{I}{I_{0}}=1-\frac{N^{\prime}}{N}$.
Finally, for a constant generation rate $G$, the CL intensity evolves as:

$\frac{I}{I_{0}}=1-\frac{(1-a)}{2}\int_{0}^{+\infty}e^{-\frac{z^{\prime}}{L}}g(z^{\prime})dz^{\prime}$
(7)

Various analytical functions have been proposed to describe $g(z)$, the in-
depth profile of excitation by an electron beam as function of its
acceleration voltage [34, 35]. Their success was limited for quantitative
analysis of cathodoluminescence experiments because their validity is
restricted to limited ranges of acceleration voltage and material density.
Today, Monte Carlo (MC) algorithms are preferred since they provide
simulations of electron energy losses in a simple and relatively accurate way
[36]. As shown in the Supplementary Materials, the in-depth distribution of
the generation rate can be assessed numerically with Monte-Carlo simulations
as well as the acceleration voltage dependence of the mean excitation depth
$\langle z_{e}\rangle=\int_{0}^{+\infty}zg(z)dz$. An example of calculated
profile for a 5kV acceleration voltage is shown in Figure 2 (a) together with
the corresponding value of $\langle z_{e}\rangle$.

Considering its shape, this profile can be approximated by a square profile
also shown in Figure 2(a) and defined by a depth distribution
$g(z)=\sqcap(z)$, with a constant generation rate down to an electron stopping
depth $R=2\langle z_{e}\rangle$, i.e. :

$\sqcap(z)=\left\\{\begin{array}[]{ll}1/R&\mbox{if }0\leqslant z\leqslant R\\\
0&\mbox{if }z>R\end{array}\right.$

Using this simpler depth distribution is of key interest since its provides an
analytical expression of $\frac{I}{I_{0}}$ involving three parameters namely
$a$, $L$ and $\langle z_{e}\rangle$ :

$\frac{I}{I_{0}}=1-\frac{(1-a)}{2}\;\frac{1-e^{-\frac{2\langle
z_{e}\rangle}{L}}}{\frac{2\langle z_{e}\rangle}{L}}$ (8)

Figure 2(b) depicts the variations of $\frac{I}{I_{0}}$ as a function of
$\langle z_{e}\rangle$ for different values of $a$ and $L$ as generated from
the analytical equation Eq. 8. They are compared with numerical calculations
of $\frac{I}{I_{0}}$ obtained by discretising the integral term in Eq. 7 to
inject the generation profiles $g(z)$ obtained from MC simulations.

Figure 2:  (a) Generation profile for a 5 kV incident electron beam in hBN
calculated by MC simulations (in black), and its square approximation (in
red). Both profiles share the same mean excitation depth $\langle
z_{e}\rangle$ calculated by MC simulation. (b) $I/I_{0}$ as a function of
$\langle z_{e}\rangle$ for a steady excitation inside the crystal. The results
are obtained either numerically with discrete generation profiles from MC
simulations (cross) or analytically with Eq. 8 considering the square profile.

This comparison attests that an excellent agreement is found between
analytical and numerical variations. This remarkable result shows that the
precise shape of the generation profile $g(z)$ has a negligible effect on the
luminescence intensity, while the mean excitation depth $\langle z_{e}\rangle$
turns out to be the critical parameter for the interpretation of these CL
experiments. This observation holds true even for GaN, which has a
significantly higher density ($6.15g.cm^{-3}$) than hBN ($2.18g.cm^{-3}$),
resulting in drastically different excitation profiles (for further details
see the Supplementary Materials). Thus, it turns out that the analytical
description of Eq. 8 can be applied to a wide range of semiconductor
materials, as soon as the CL reabsorption by the crystal remains negligible,
which is the case in indirect band gap semiconductors for instance. 333The
present model is inspired by the early works of Wittry et al. in 1967. [38]
Despite the title of their paper, “Measurement of Diffusion Lengths in Direct
Gap Semiconductors by Electron Beam Excitation”, the present work rather
applies to indirect semiconductors, since the efficient luminescence
reabsorption in direct bandgap semiconductors requires careful consideration
when evaluating luminescence intensity.

Figure 2(b) also reveals that the parameters $a$ and $L$ have distinct
influences on the $\frac{I}{I_{0}}(\langle z_{e}\rangle)$ curve. In
particular, the parameter $a$ determines the shape of the curve at low
excitation depths, while the parameter $L$ governs the position of the curve
along the horizontal axis. This observation provides the guidance for the
analysis of experimental measurements. Provided to vary the acceleration
voltage over a sufficiently large range, the analytical description of Eq. 8
has the ability to accurately fit experimental measurements of the CL
intensity independently with reliable values of the diffusion length $L$ and
the parameter $a$ which quantifies the level of surface recombinations.
Further, combining these experiments with measurements of the bulk exciton
lifetime $\tau$ using TRCL at high voltage as introduced above gives access to
the evaluation of the diffusivity $D$ and the surface recombination velocity
$s$ which are the key physical quantities of interest.

##  .3 III - Experimental results for free excitons in hBN : analysis and
discussion

We now turn to the experiments. Four bulk hBN crystals were investigated. One
was grown at the Kansas State University using an Atmospheric Pressure High
Temperature process [39, 40, 41] from Ni/Cr solvent and natural boron isotop
content (APHT sample). Two crystals were synthesised at the NIMS laboratory
using a High Pressure High Temperature route [42, 43] (samples HPHT1 and
HPHT2). And the last crystal was grown in the LMI laboratory using a Polymer
Derived Ceramic method [44, 45] (PDC sample). The HPHT samples are widely
recognised in the 2D material scientific community as reference hBN materials.
Prior to analysis, the samples were cleaved to obtain a clean, uncontaminated
surface, and then mounted on conducting doped silicon wafers. A conductive
bond was established on one edge of the crystal using Ag paint. This mounting
allows an efficient charge evacuation during electron beam irradiation. The
thickness of the crystals was measured to be greater than 10 $\mu$m using a
mechanical profilometer.

Cathodoluminescence spectra were collected out using a Scanning Electron
Microscope (SEM) JEOL7001F equipped with a Horiba Jobin-Yvon CL system
optimised for UV spectroscopy, as described in detail in reference [4]. To
avoid non-linear effects such as exciton-exciton annihilation [46], a low
excitation power of about 1 $\mu$W was employed. To achieve a 1D geometry
configuration (shown in Figure 1), the electron beam was significantly
defocused, resulting in an excitation spot size of 3.5 $\mu$m. The excitation
depth was controlled by varying the electron acceleration voltage from 500 V
to 20 kV. For each acceleration voltage, the CL intensity $I$ measurements
were normalised by the excitation power $P=iV\cdot(1-f)$, taking into account
slight variations ($\pm$10%) in the incident electron beam power ($iV$) at
each acceleration voltage $V$, as well as the voltage dependence of the
electron backscattering factor $f$ (see Supplementary Materials). The
uncertainty in the CL intensity was estimated to be 20%. Finally, time-
resolved CL was performed using a custom-built beam blanker installed in the
SEM column, achieving a temporal resolution of 100 ps (details of the setup
can be found in [5]).

Figure 3: (a) CL intensity $I/I_{0}$ as a function of the acceleration voltage
and the corresponding mean excitation depth $\langle z_{e}\rangle$ (see
Supplementary Materials). The analytical equation (8) is depicted for a
diffusion length $L$ of 92 nm and various values of the parameter $a$. (b)
Temporal decay of the free exciton luminescence measured in the same area
after the interruption of a 15 kV electron beam. The bulk lifetime of the free
exciton $\tau$ is determined by fitting the data with an asymptote.

Figure 3 presents the first application of the protocol to the APHT sample. In
Figure 3(a), the CL intensity $I/I_{0}$ is plotted as a function of the mean
excitation depth $\langle z_{e}\rangle$ as extracted from Monte Carlo
simulations. At high excitation depths, the intensity plateaus indicate that
excitons are generated too far from the top surface to cause significant
surface effects. The decrease in luminescence intensity at low excitation
depths is attributed to surface recombinations enhanced by exciton diffusion
towards the top surface. Fit of the experimental data by the equation 8 is
discussed from Figure 3(a). Taking into account the uncertainties, the 95%
confidence interval for $L$ is between 84 nm and 100 nm, and the best fit is
with $L=92$ nm as shown in the Figure 3(a). This Figure also illustrates that
only a few points at low excitation depths contribute to the determination of
the parameter $a$. Equation 8 is displayed for $L=92$ nm and $a=-1$, $a=-0.9$,
and $a=-0.8$. Reasonable agreement is found for $a<-0.9$, while higher values
lead to a poorer fit, establishing an upper limit of $-0.9$ for $a$.

These results make it possible to evaluate and discuss the diffusivity $D$ and
the surface recombination velocity $s$ by starting with the latter. The
surface recombination velocity $s=\frac{L}{\tau}\frac{1-a}{1+a}$ is known from
$L$, $\tau$, and $a$. The exciton lifetime $\tau$ is measured independently by
time-resolved CL with an uncertainty of 10%, following the procedure described
in [5] and shown in Figure 3(b). Knowing the upper limit of $a$, a lower limit
for the surface recombination velocity $s$ is deduced, yielding $s>7\times
10^{4}$ cm/s for the APHT sample.

The high value of $s$ is very robust as it could be confirmed from the
analysis of data obtained from the different samples and at various
temperatures (as shown in Tables 1 and 2), leading to refinement of its lower
limit to $10^{5}$ cm/s at 6 and 300 K. This is definetely a surprinsingly high
value, approaching the highest levels reported around $10^{6}$ cm/s for 3D
semiconductors such as in GaAs [12] and Diamond [13]. A possible explanation
for this result could be the presence of defects at the hBN surface. While it
is challenging to completely rule out this possibility, note that the surfaces
studied in this work were clean and freshly cleaved, minimising potential
contamination. In addition, there is no change in the shape of the spectra
with the depth of excitation, suggesting a defect concentration homogeneous in
depth and a surface free of contamination, as shown in the spectra in the
appendices. Another possible explanation relates to the surface exciton states
that have been identified theoretically [47] and experimentally [4] in hBN,
and which are approximately 100 meV below the bulk hBN exciton. These surface
states suggest that the hBN surface acts as a sink for excitons, with a
sufficient energy depth to capture them efficiently at room temperature,
consistently with the high surface recombination rate observed in hBN.

Whatever its origin, the importance of surface recombination in hBN allows us
to understand the difficulty of measuring the luminescence of its atomic
layers. In the presence of such a high surface recombination rate, the
luminescence of hBN crystals made of a few atomic layers is inevitably very
low, as experimentally measured by L. Schué et al. [6]. Contrary to initial
assumptions, the surfaces of 2D materials can also limit the efficiency of
their luminescence, emphasising the need to passivate them for use in
optoelectronic devices. Recent studies on atomic layers of transition metal
dichalcogenides (TMDs) have demonstrated that the luminescence quantum yields
only reach unity when the material is properly passivated [48, 49]. The
increase in quantum efficiency following passivation indicates that non-
radiative recombinations occur at the TMD-air surface. Similar passivation
techniques should be considered to mitigate non-radiative recombinations at
the hBN surface, therby enabling enhanced luminescence intensity and improving
the potential of hBN for optoelectronic applications.

We now focus on the diffusion properties. The measured diffusion lengths $L$
and the exciton lifetimes $\tau$ on the four samples at 300 K are depicted in
Figure 4 (a-b). In the case of the PDC sample, the exciton lifetime is too
short to be measured directly by TRCL. However, using the linearity
relationship between $\tau$ and the CL efficiency given in reference [5], the
exciton lifetime $\tau$ on this sample was estimated as 0.03 $\pm$ 0.015 ns
from the measurement of a 0.1% CL efficiency (details of the measurement in
reference [4]). For all other samples, $\tau$ is measured directly by TRCL
with a 10% uncertainty. Considering the uncertainty of the CL intensity
measurement, the 95% confidence intervals of $L$ are always found within $\pm$
15%. Finally, the diffusivity $D$ is extracted considering the relation
$D=L^{2}/\tau$ and taking into account the uncertainty on $L$ and $\tau$. The
results are reported in the Table 1.

Figure 4: (a) $I/I_{0}$ as function of the mean excitation depth $\langle
z_{e}\rangle$, and (b) temporal decays of the free exciton luminescence after
interruption of a 15 kV steady excitation at $t$=0, recorded on HPHT1, HPHT2,
APHT and PDC samples. The out-of-plane diffusion length $L$ is deduced by
fitting the equation (8) in (a), while the bulk exciton lifetime $\tau$ is
measured by TRCL in (b). (c) Diffusivity as function of the exciton lifetime
$\tau$ for the 4 samples. Note that $\tau$ is used to quantitatively compare
the global quality of the samples.

As explained in reference [5], the higher the quality of the sample, the
longer the exciton lifetime is. In the Figure 4 (c), the exciton lifetime
$\tau$ is used to quantitavely compare the globale quality of the 4 samples.
The figure exhibits a variation of the diffusivity $D$ of more than one order
of magnitude between the different samples and clearly demonstrates that the
diffusivity is related to the defects present in the samples.

Table 1: Results of the CL experiments on the four bulk hBN samples at 300 K : the out-of-plane exciton diffusion length $L$, the parameter $a$ and the exciton lifetime $\tau$. An estimate of the diffusivity $D=L^{2}/\tau$ and a lower bound of the surface recombination velocity $s=\frac{L}{\tau}\frac{1-a}{1+a}$ are derived. |  $L$ (nm) | $a$ |  $\tau$ (ns) |  $s$ (cm/s) |  $D$ (cm2/s)
---|---|---|---|---|---
HPHT1 | 570 |  $<$ -0.98 | 16.5 |  $>$3*105 | 0.20 $\pm$ 0.07
HPHT2 | 246 |  $<$ -0.96 | 8.4 |  $>$1*105 | 0.07 $\pm$ 0.03
APHT | 92 |  $<$ -0.9 | 2.2 |  $>$7*104 | 0.04 $\pm$ 0.01
PDC |  $<$ 5 | - | 0.03 | - |  $<$ 0.02

Figure 5 presents the measurements of the diffusion length $L$ and exciton
lifetime $\tau$ on the best sample (HPHT1) at both 6 K and 300 K using our
experimental protocol. The corresponding results are summarised in Table 2.
The diffusivity is 0.25 $\pm$ 0.08 cm2/s at 6 K and 0.20 $\pm$ 0.07 cm2/s at
the room temperature. The difference is not significant given the
uncertainties of the measurement. In the case of intrinsic diffusion limited
by phonon scattering, the diffusion constant at cryogenic temperatures is
expected to be several orders of magnitude higher than at room temperature
[50]. Our results indicate that even in the highest-quality HPHT1 sample,
exciton diffusion is not limited by phonon scattering but rather by defects
scattering. This is consistent with impurity concentrations in oxygen and
carbon in the order of 1018 cm-3 or below measured in the reference HPHT
samples [43, 51] a value which remains high compared to ultrapure
semiconductors (1011-1013 cm-3 typically) where the intrinsic diffusivity is
reached [50].

Figure 5: (a) $I/I_{0}$ as function of the mean excitation depth $\langle z_{e}\rangle$ during continuous excitation and (b) temporal decays of the free exciton luminescence after interruption of the 15 kV electron beam, recorded on the HPHT1 sample at 6 and 300 K.  Table 2: Results of the CL experiments at 6 K and 300 K on the HPHT1 sample. T° (K) |  $L$ (nm) | $a$ |  $\tau$ (ns) |  $s$ (cm/s) |  $D$ (cm2/s)
---|---|---|---|---|---
300 K | 570 |  $<$ -0.98 | 16.5 |  $>$3*105 | 0.20 $\pm$ 0.07
6 K | 275 |  $<$ -0.95 | 3 |  $>$3*105 | 0.24 $\pm$ 0.08

An estimate of the lower limit of the intrinsic diffusivity in hBN is provided
by the highest value measured on our sample : $D$ $>$ 0.2 cm2/s. Few values of
the exciton diffusion constant are available in the literature for comparison,
since excitons in standard semiconductors are weakly bounded and generally
dissociated at room temperature. Among them, the diffusion constant is 5 cm2/s
in diamond at 300 K [50]. In organic semiconductors where the diffusion is
limited by exciton hoping from site to site, $D$ is found in the range 10-6 \-
0.01 cm2.s-1 [52]. The lower bound found for $D$ in hBN indicates that the
hypothesis of a hoping diffusion mechanism can be safely discarded in this
material, which has a diffusion constant typical of a standard semiconductor.

##  .4 IV - Conclusion

In conclusion, we have developed a new CL protocol thanks to which the out-of-
plane exciton diffusion length $L$, the diffusivity $D$ and the surface
recombination velocity $s$ have been recorded for the first time in a 2D
semiconductor. It consists in recording the emitted light as function of the
acceleration voltage over a large range and using a widely defocused beam to
be in 1D geometry. By a detailed analysis of the 1D diffusion equation
applying to these CL experiments, we show that the emitted light in the
presence of surface recombinations and diffusion is not affected by the in-
depth profile of the exciton generation function, but by the mean excitation
depth $\langle z_{e}\rangle$ identified as the key parameter for our
experiment. A simple analytical function is derived to describe the
experimental results and to estimate the physical quantities in turn. The
application of the protocol to four bulk hBN samples reveals diffusion lengths
$L$ up to 570 nm at 300 K in the best crystal. From temperature dependent
experiments, the exciton diffusivity is shown to be limited by defects
scattering. Its lower limit for the intrinsic diffusivity is 0.2 cm2/s,
similar to the value expected for a standard semiconductor. Finally, the hBN
surface recombination velocity is extremely high, at the level of silicon or
diamond. This surprising result explains the low intensity of hBN thin films,
and reveals that exciton recombination at surfaces could be a strong
limitation for light-emitting devices based on 2D materials.

###### Acknowledgements.

The research leading to these results has received funding European Union’s
Horizon 2020 research and innovation program under grant agreements No 785219
(Graphene Core 2) and No 881603 (Graphene Core 3). Support for the APHT hBN
crystal growth comes from the Office of Naval Research, Award No.
N00014-20-1-2474. K.W. and T.T. acknowledge support from the JSPS KAKENHI
(Grant Numbers 21H05233 and 23H02052) and World Premier International Research
Center Initiative (WPI), MEXT, Japan.

## References

  * Mak  _et al._ [2010] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Atomically thin MoS2: a new direct-gap semiconductor, Physical Review Letters 105, 136805 (2010).
  * Tonndorf  _et al._ [2013] P. Tonndorf, R. Schmidt, P. Böttger, X. Zhang, J. Börner, A. Liebig, M. Albrecht, C. Kloc, O. Gordan, D. R. T. Zahn, S. M. de Vasconcellos, and R. Bratschitsch, Photoluminescence emission and raman response of monolayer MoS2, MoSe2, and WSe2, Optics Express 21, 4908 (2013).
  * Zhao  _et al._ [2012] W. Zhao, Z. Ghorannevis, L. Chu, M. Toh, C. Kloc, P.-H. Tan, and G. Eda, Evolution of electronic structure in atomically thin sheets of WS2 and WSe2, ACS Nano 7, 791 (2012).
  * Schué  _et al._ [2019] L. Schué, L. Sponza, A. Plaud, H. Bensalah, K. Watanabe, T. Taniguchi, F. Ducastelle, A. Loiseau, and J. Barjon, Bright luminescence from indirect and strongly bound excitons in hBN, Physical Review Letters 122, 067401 (2019).
  * Roux  _et al._ [2021] S. Roux, C. Arnold, F. Paleari, L. Sponza, E. Janzen, J. H. Edgar, B. Toury, C. Journet, V. Garnier, P. Steyer, T. Taniguchi, K. Watanabe, F. Ducastelle, A. Loiseau, and J. Barjon, Radiative lifetime of free excitons in hexagonal boron nitride, Physical Review B 104, L161203 (2021).
  * Schué  _et al._ [2016] L. Schué, B. Berini, A. C. Betz, B. Plaçais, F. Ducastelle, J. Barjon, and A. Loiseau, Dimensionality effects on the luminescence properties of hBN, Nanoscale 8, 6986 (2016).
  * Elias  _et al._ [2019] C. Elias, P. Valvin, T. Pelini, A. Summerfield, C. Mellor, T. Cheng, L. Eaves, C. Foxon, P. Beton, S. Novikov, _et al._ , Direct band-gap crossover in epitaxial monolayer boron nitride, Nature communications 10, 2639 (2019).
  * Rousseau  _et al._ [2021] A. Rousseau, L. Ren, A. Durand, P. Valvin, B. Gil, K. Watanabe, T. Taniguchi, B. Urbaszek, X. Marie, C. Robert, _et al._ , Monolayer boron nitride: Hyperspectral imaging in the deep ultraviolet, Nano Letters 21, 10133 (2021).
  * Nolte [1990] D. D. Nolte, Surface recombination, free-carrier saturation, and dangling bonds in InP and GaAs, Solid-State Electronics 33, 295 (1990).
  * Hsu  _et al._ [1992] J. W. P. Hsu, C. C. Bahr, A. vom Felde, S. W. Downey, G. S. Higashi, and M. J. Cardillo, Si(111)(7×7) dangling bond contribution to surface recombination, Journal of Vacuum Science and Technology A 10 (1992).
  * Fitzgerald and Grove [1968] D. Fitzgerald and A. Grove, Surface recombination in semiconductors, Surface Science 9, 347 (1968).
  * Jastrzebski  _et al._ [1975] L. Jastrzebski, J. Lagowski, and H. C. Gatos, Application of scanning electron microscopy to determination of surface recombination velocity: GaAs, Applied Physics Letters 27, 537 (1975).
  * Kozák  _et al._ [2012] M. Kozák, F. Trojánek, and P. Malý, Large prolongation of free-exciton photoluminescence decay in diamond by two-photon excitation, Optics Letters 37, 2049 (2012).
  * High  _et al._ [2008] A. A. High, E. E. Novitskaya, L. V. Butov, M. Hanson, and A. C. Gossard, Control of exciton fluxes in an excitonic integrated circuit, Science 321, 229 (2008).
  * Butov [2017] L. V. Butov, Excitonic devices, Superlattices and Microstructures 108, 2 (2017).
  * Jauregui and Kim [2017] L. A. Jauregui and P. Kim, Curved paths of electron–hole pairs, Nature Materials 16, 1169 (2017).
  * Unuchek  _et al._ [2018] D. Unuchek, A. Ciarrocchi, A. Avsar, K. Watanabe, T. Taniguchi, and A. Kis, Room-temperature electrical control of exciton flux in a van der waals heterostructure, Nature 560, 340 (2018).
  * Cadiz  _et al._ [2018] F. Cadiz, C. Robert, E. Courtade, M. Manca, L. Martinelli, T. Taniguchi, K. Watanabe, T. Amand, A. C. Rowe, D. Paget, B. Urbaszek, and X. Marie, Exciton diffusion in WSe2 monolayers embedded in a van der waals heterostructure, Applied Physics Letters 112, 152106 (2018).
  * Kumar  _et al._ [2014] N. Kumar, Q. Cui, F. Ceballos, D. He, Y. Wang, and H. Zhao, Exciton diffusion in monolayer and bulk MoSe2, Nanoscale 6, 4915 (2014).
  * Uddin  _et al._ [2020] S. Z. Uddin, H. Kim, M. Lorenzon, M. Yeh, D.-H. Lien, E. S. Barnard, H. Htoon, A. Weber-Bargioni, and A. Javey, Neutral exciton diffusion in monolayer MoS2, ACS Nano 14, 13433 (2020).
  * Heller  _et al._ [1996] W. Heller, A. Filoramo, P. Roussignol, and U. Bockelmann, Direct measurement of exciton diffusion in quantum wells, Solid-State Electronics 40, 725 (1996).
  * Lawson  _et al._ [1982] C. M. Lawson, R. C. Powell, and W. K. Zwicker, Transient grating investigation of exciton diffusion and fluorescence quenching in NdxTa1-xPO crystals, Physical Review B 26, 4839 (1982).
  * Zarem  _et al._ [1989] H. A. Zarem, P. C. Sercel, J. A. Lebens, L. E. Eng, A. Yariv, and K. J. Vahala, Direct determination of the ambipolar diffusion length in GaAs/AlGaAs heterostructures by cathodoluminescence, Applied Physics Letters 55, 1647 (1989).
  * Ino and Yamamoto [2008] N. Ino and N. Yamamoto, Low temperature diffusion length of excitons in gallium nitride measured by cathodoluminescence technique, Applied Physics Letters 93, 232103 (2008).
  * Hocker  _et al._ [2016] M. Hocker, P. Maier, L. Jerg, I. Tischer, G. Neusser, C. Kranz, M. Pristovsek, C. J. Humphreys, R. A. R. Leute, D. Heinz, O. Rettig, F. Scholz, and K. Thonke, Determination of axial and lateral exciton diffusion length in GaN by electron energy dependent cathodoluminescence, Journal of Applied Physics 120, 85703 (2016).
  * Kenkre [1981] V. M. Kenkre, Determination of the exciton diffusion constant from variation of quantum yield with penetration length, Chemical Physics Letters 82, 301 (1981).
  * Barjon  _et al._ [2012] J. Barjon, F. Jomard, A. Tallaire, J. Achard, and F. Silva, Determination of exciton diffusion lengths in isotopically engineered diamond junctions, Applied Physics Letters 100, 122107 (2012).
  * Nogues  _et al._ [2014] G. Nogues, T. Auzelle, M. D. Hertog, B. Gayral, and B. Daudin, Cathodoluminescence of stacking fault bound excitons for local probing of the exciton diffusion length in single GaN nanowires, Applied Physics Letters 104, 102102 (2014).
  * Note [1] Assuming strongly bounded excitons, i.e. negligible free carrier concentrations.
  * Van Roosbroeck [1955] W. Van Roosbroeck, Injected current carrier transport in a semi‐infinite semiconductor and the determination of lifetimes and surface recombination velocities, Journal of Applied Physics 26, 380 (1955).
  * Note [2] $g(z)=\delta(z-z^{\prime})$ for a plane source in $z=z^{\prime}$.
  * Klein [1968] C. A. Klein, Bandgap dependence and related features of radiation ionization energies in semiconductors, Journal of Applied Physics 39, 2029 (1968).
  * Yacobi and Holt [2013] B. G. Yacobi and D. B. Holt,  _Cathodoluminescence microscopy of inorganic solids_ (Springer Science & Business Media, 2013).
  * Kanaya and Okayama [1972] K. Kanaya and S. Okayama, Penetration and energy-loss theory of electrons in solid targets, Journal of Physics D: Applied Physics 5, 43 (1972).
  * Sieber [2012] B. Sieber, Cathodoluminescence - principes physiques et systèmes de détection, Techniques de l’ingénieur, Analyses de surface et de materiaux , Réf : P3792 v1 (2012).
  * Drouin  _et al._ [2007] D. Drouin, A. R. Couture, D. Joly, X. Tastet, V. Aimez, and R. Gauvin, CASINO V2.42—a fast and easy-to-use modeling tool for scanning electron microscopy and microanalysis users, Scanning 29, 92 (2007), https://onlinelibrary.wiley.com/doi/pdf/10.1002/sca.20000 .
  * Note [3] The present model is inspired by the early works of Witry et al. in 1967. [38] Despite the title of their paper, “Measurement of Diffusion Lengths in Direct Gap Semiconductors by Electron Beam Excitation”, the present work rather applies to indirect semiconductors, since the efficient luminescence reabsorption in direct bandgap semiconductors requires careful consideration when evaluating luminescence intensity.
  * Wittry and Kyser [1967] D. B. Wittry and D. F. Kyser, Measurement of diffusion lengths in direct‐gap semiconductors by electron‐beam excitation, Journal of Applied Physics 38, 375 (1967).
  * Liu  _et al._ [2018] S. Liu, R. He, L. Xue, J. Li, B. Liu, and J. H. Edgar, Single crystal growth of millimeter-sized monoisotopic hexagonal boron nitride, Chemistry of Materials 30, 6222 (2018).
  * Li  _et al._ [2020a] J. Li, C. Elias, G. Ye, D. Evans, S. Liu, R. He, G. Cassabois, B. Gil, P. Valvin, B. Liu, _et al._ , Single crystal growth of monoisotopic hexagonal boron nitride from a fe–cr flux, Journal of Materials Chemistry C 8, 9931 (2020a).
  * Kubota  _et al._ [2007] Y. Kubota, K. Watanabe, O. Tsuda, and T. Taniguchi, Deep ultraviolet light-emitting hexagonal boron nitride synthesized at atmospheric pressure, Science 2007, 932 (2007).
  * Watanabe  _et al._ [2004] K. Watanabe, T. Taniguchi, and H. Kanda, Direct-bandgap properties and evidence for ultraviolet lasing of hexagonal boron nitride single crystal, Nature Materials 3, 404 (2004).
  * Taniguchi and Watanabe [2007] T. Taniguchi and K. Watanabe, Synthesis of high-purity boron nitride single crystals under high pressure by using Ba–BN solvent, Journal of Crystal Growth 303, 525 (2007).
  * Li  _et al._ [2018] Y. Li, V. Garnier, C. Journet, J. Barjon, A. Loiseau, I. Stenger, A. Plaud, B. Toury, and P. Steyer, Advanced synthesis of highly crystallized hexagonal boron nitride by coupling polymer-derived ceramics and spark plasma sintering processes—influence of the crystallization promoter and sintering temperature, Nanotechnology 30, 035604 (2018).
  * Li  _et al._ [2020b] Y. Li, V. Garnier, P. Steyer, C. Journet, and B. Toury, Millimeter-scale hexagonal boron nitride single crystals for nanosheet generation, ACS Applied Nano Materials 3, 1508 (2020b).
  * Plaud  _et al._ [2019] A. Plaud, L. Schué, K. Watanabe, T. Taniguchi, A. Loiseau, and J. Barjon, Exciton-exciton annihilation in hBN, Applied Physics Letters 114, 232103 (2019).
  * Paleari  _et al._ [2018] F. Paleari, T. Galvani, H. Amara, F. Ducastelle, A. Molina-Sánchez, and L. Wirtz, Excitons in few-layer hexagonal boron nitride: Davydov splitting and surface localization, 2D Materials 5, 045019 (2018).
  * Lien  _et al._ [2019] D.-H. Lien, S. Z. Uddin, M. Yeh, M. Amani, H. Kim, J. W. Ager III, E. Yablonovitch, and A. Javey, Electrical suppression of all nonradiative recombination pathways in monolayer semiconductors, Science 364, 468 (2019).
  * Amani  _et al._ [2016] M. Amani, P. Taheri, R. Addou, G. H. Ahn, D. Kiriya, D.-H. Lien, J. W. Ager III, R. M. Wallace, and A. Javey, Recombination kinetics and effects of superacid treatment in sulfur-and selenium-based transition metal dichalcogenides, Nano letters 16, 2786 (2016).
  * Morimoto  _et al._ [2015] H. Morimoto, Y. Hazama, K. Tanaka, and N. Naka, Ultrahigh exciton diffusion in intrinsic diamond, Physical Review B 92, 201202(R) (2015).
  * Onodera  _et al._ [2019] M. Onodera, K. Watanabe, M. Isayama, M. Arai, S. Masubuchi, R. Moriya, T. Taniguchi, and T. Machida, Carbon-rich domain in hexagonal boron nitride: Carrier mobility degradation and anomalous bending of the landau fan diagram in adjacent graphene, Nano Letters 19, 7282 (2019).
  * Mikhnenko  _et al._ [2015] O. V. Mikhnenko, P. W. M. Blom, and T.-Q. Nguyen, Exciton diffusion in organic semiconductors, Energy and Environmental Science 8, 1867 (2015).


# Root Cause Identification for Collective Anomalies in Time Series given an
Acyclic Summary Causal Graph with Loops

Charles K. Assaad
EasyVista &Imad Ez-zejjari
EasyVista & Lei Zan
EasyVista
Univ Grenoble Alpes,
CNRS, Grenoble INP, LIG

###### Abstract

This paper presents an approach for identifying the root causes of collective
anomalies given observational time series and an acyclic summary causal graph
which depicts an abstraction of causal relations present in a dynamic system
at its normal regime. The paper first shows how the problem of root cause
identification can be divided into many independent subproblems by grouping
related anomalies using d-separation. Further, it shows how, under this
setting, some root causes can be found directly from the graph and from the
time of appearance of anomalies. Finally, it shows, how the rest of the root
causes can be found by comparing direct effects in the normal and in the
anomalous regime. To this end, an adjustment set for identifying direct
effects is introduced. Extensive experiments conducted on both simulated and
real-world datasets demonstrate the effectiveness of the proposed method.

##  1 INTRODUCTION

The need for high availability of information systems requires efficient
monitoring tools and a new generation of AIOps software to automate the
identification of actionable root causes of anomalies in an IT monitoring
system that can be used to eliminate the anomalies. In recent years, many
approaches have been developed to identify the root causes of anomalies in
multivariate time series. The most common direction that is explored considers
discovering the causal graph (Pearl, 2000; Spirtes et al., 2000) that
represents the anomalous regime of the dynamic system (Wang et al., 2018; Meng
et al., 2020) using observational time series. Causal discovery methods are
known to rely on strong assumptions which imply that they necessitate
validation by an expert especially when there is no guarantee that these
assumptions are satisfied. In addition, causal discovery methods usually need
large sample sizes (Malinsky and Danks, 2018; Glymour et al., 2019; Assaad et
al., 2022a). However, in many domains, the size of anomalous data depends on
the sampling rate of the system. Thus systems with low sampling rates will
collect a small size of anomalous data compared to systems with high sampling
rates. In consequence, sometimes, it is difficult to have enough data for
causal discovery methods. And even if a sufficient amount of data is
collected, the process of validation of the causal graph by an expert is time
consuming and can delay the elimination of anomalies.

To tackle this issue, we follow a different approach for root cause analysis
which consists of discovering and reasoning about the summary causal graphs
(Assaad et al., 2022a) which depicts an abstraction of causal relations
present in a dynamic system at its normal regime. Usually, the size of data
collected in the normal regime is significantly greater than the size of data
collected in the anomalous regime since anomalies are supposed to be rare. In
addition, system experts can have sufficient time to validate the graph long
before the appearance of anomalies. Note that in this work, we do not adress
the problem of causal discovery of the summary causal graph of the normal
regime and we assume that the graph is already learned and validated by a
system expert.

This paper presents a new method for root cause identification, which we call
EasyRCA, which consists of using a summary causal graph of the normal regime
in order to divide the problem of root cause identification into many
independent subproblems by grouping related anomalies using d-separation. Then
for each group, EasyRCA finds the root causes either directly from the graph
and the time of appearance of anomalies or by comparing direct effects in the
normal and anomalous regime. To this end, an adjustment set for identifying
direct effects is introduced.

The remainder of the paper is organized as follows: Section 2 introduces some
terminology and formalizes the problem. Section 3 describes related work.
Section 4 presents our method EasyRCA which is evaluated on simulated and real
datasets in Section 5. Finally, Section 6 concludes the paper.

##  2 Problem setup

In this section, we first introduce some terminology, tools, and assumptions
which are standard for the major part. Then, we formalize the problem we are
going to solve.

Suppose that a dynamic system can be represented by a structural causal model
(SCM) (Pearl, 2000) in which each point in a time series is given by a
function (so-called causal mechanism) of its parents and an unobserved noise:

$Y_{t}:=f^{y}_{t}(Parents(Y_{t}),\xi^{y}_{t})$ (1)

where the noise variables are jointly statistically independent so that there
are no hidden confounding , i.e., causal sufficiency is satisfied (Spirtes et
al., 2000). The qualitative causal relations induced by such SCM can be
represented by a causal graph in which, under the causal Markov condition
(Spirtes et al., 2000) each vertex is independent of all other vertices given
its parents, except for its descendants. In dynamic systems, these causal
graphs are referred to as full-time causal graphs. An example of such a graph
is presented in Figure 1(a). The main difficulty in working with this type of
graph is that it is infinite and so in practice, inferring it is unfeasible.
However, it is very likely that causal relations between two time series will
hold throughout time as such relations are generally associated with
underlying physical processes. Thus we can assume consistency throughout time.

######  Assumption 1 (Consistency throughout time, (Assaad et al., 2022a)).

A full time causal graph is said to be _consistent throughout time_ if all the
causal relationships remain constant in direction throughout time.

Under this assumption, the full time causal graph can be contracted to give a
finite graph which is called a window causal graph. It is a representation of
the causal relations through a time window, the size of which depends on the
maximum lag between a cause and an effect in the full time causal graph. An
example of a window causal graph is given in Figure1(b). This said, it is
usually difficult for an expert to validate, analyze let alone provide a
window causal graph because it is difficult to determine the temporal lag
between a cause and an effect. Thus, experts usually rely on the so-called
summary causal graph which is a compact version of the window causal graph
that represents the causal relations between time series without giving any
information about the temporal lags of these relations. In this work, we
assume that the summary causal graph is acyclic but loops are allowed to
represent temporal dependencies within the same time series. An example of
such graph is given in Figure 1(c) and formally it is defined as follows:

######  Definition 1 (Acyclic summary causal graph with loops).

Consider $\mathcal{G}=(\mathcal{V},\mathcal{E})$ is a _summary causal graph_.
The set of vertices in that graph consists of the set of time series. The arcs
$\mathcal{E}$ of the graph are defined as follows: $\forall
X,Y\in\mathcal{V}$, $X$ causes $Y$ if and only if there exists some time lag
$\gamma$ such that $X_{t-\gamma}$ causes $Y_{t}$ such that $\gamma\geq 0$ for
$X\neq Y$ and $\gamma>0$ for $X=Y$. If $\mathcal{G}=(\mathcal{V},\mathcal{E})$
has no directed cycles other than the edges going from one vertex to itself,
then $\mathcal{G}=(\mathcal{V},\mathcal{E})$ is said to be an _acyclic summary
causal graph with loops_ (ASCGL).

In this work, we suppose that the ASCGL is given either using experts
knowledge or learned directly from observational time series (Peters et al.,
2013; Assaad et al., 2021, 2022b) or by first discovering a window causal
graph111In our framework, if one decides using a causal discovery algorithm
that learns a window causal graph, then one needs to carefully incorporate to
these algorithms the constraint that the summary causal graphs should be
acyclic. (Runge et al., 2019; Runge, 2020) and then deduce the ASCGL from it.
The correctness of such learned graphs usually rests on untestable assumptions
and depends on the quality of the data so it is important to validate it by an
expert or to simplify the problem of causal discovery by providing some
background knowledge.

$Z_{t-2}$$Z_{t-1}$$Z_{t}$$X_{t-2}$$X_{t-1}$$X_{t}$$W_{t-2}$$W_{t-1}$$W_{t}$$Y_{t-2}$$Y_{t-1}$$Y_{t}$
(a) Full time causal graph

$Z_{t-1}$$Z_{t}$$X_{t-1}$$X_{t}$$W_{t-1}$$W_{t}$$Y_{t-1}$$Y_{t}$ (b) Window
causal graph

$X$$Z$$W$$Y$ (c) Summary causal graph

Figure 1: Different causal graphs that one can infer from three time series:
full time causal graph (1(a)), window causal graph (1(b)) and summary causal
graph (1(c)). Note that the first one gives more information but cannot be
inferred in practice, the second one is a schematic viewpoint of the full
behavior, whereas the last one is an abstraction and can be deduced from the
window causal graph.

Now, we turn our focus on anomalies. In this work, we assume that anomalies
are collective.

######  Definition 2 (Collective anomaly, Chandola et al. (2009)).

In time series, a collective anomaly is a sequence of data instances that is
anomalous with respect to the entire time series.

Point anomalies are disregarded because we are interested in finding
actionable root causes that can eliminate anomalies. If an anomaly appears for
one time instant and then disappears this means it was eliminated on its own
and do not require any action to resolve it. However, we also consider that
those collective anomalies have a limited size since information systems are
expected to be highly available (without anomalies). In addition, we assume
that each anomaly can be eliminated by removing the intervention that caused
it directly or via a causal path. This is given by the following assumption:

######  Assumption 2.

All anomalies are propagated from an external intervention through the
structural causal model.

Given the definition of ASCGL and collective anomalies, we define root
causes222The ”root causes” are relative to the set of observed time series. as
follows:

######  Definition 3 (Root causes).

Given an ASCGL and a set of anomalous vertices $\mathcal{A}$, the set of root
causes $\mathcal{C}$ of $\mathcal{A}$ is a set of vertices that were affected
by an external intervention which led to marginal distribution change in
$\mathcal{A}$.

In the literature, there exists two sorts of interventions (Eberhardt and
Scheines, 2007) and both are crucial to root cause analysis. The first is
known as parametric intervention and it is defined as follows:

######  Definition 4 (Parametric intervention).

Consider an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$. An intervention on
a vertex $Y\in\mathcal{V}$ is parametric if the causal mechanism before
intervention is different than after the intervention but
$Parents(Y,\mathcal{G})$ remains unchanged.

The second type of intervention is known as structural intervention333Our
definition of structural intervention is less restrictive than the classical
definition which states that the intervention alone completely determines the
probability distribution of the variable that underwent the intervention,
i.e., this variable becomes independent of all of its parents. and it is
defined as follows:

######  Definition 5 (Structural intervention).

Consider an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$. An intervention on
a vertex $Y\in\mathcal{V}$ is structural if $\exists X\in\mathcal{V}$ such
that $X\in Parents(Y,\mathcal{G})$ before the intervention and $X\not\in
Parents(Y,\mathcal{G})$ after the intervention.

Structural interventions can be regarded as a special case of parametric
interventions. But we distinguish between them because an expert might want to
differentiate between interventions that provoke a disruption in the system
from the ones that do not. As it will be shown in Section 4, to find these two
types of interventions we will estimate the direct effects in the normal and
anomalous regime. In consequence, we assume that causal mechanisms are fixed
throughout time within the same regime (e.g., in Equation 1, $f^{y}(.)$ is
fixed for all $t$) and we assume the minimality condition444The minimality
condition is usually assumed by causal discovery methods either directly or by
assuming a stronger assumption called faithfulness which implies the
minimality condition (Glymour et al., 2019). (Spirtes et al., 2000) which
implies that adjacent vertices in the ASCGL are statistically dependent in the
normal regime. Finally, to simplify the problem we assume linear SCMs.

Now that we have introduced the needed tools and assumptions, the problem we
are trying to solve is formalized as follows:

######  Problem.

Given an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$, a set of anomalous
vertices $\mathcal{A}\subset\mathcal{V}$, the distribution of the time series
in the normal regime $\mathcal{{N}}$ and in the anomalous regime
$\mathcal{\bar{N}}$, and the maximal lag between a cause and an effect
$\gamma_{max}$, we want to find the smallest set of root causes $\mathcal{C}$
of $\mathcal{A}$.

##  3 RELATED WORKS

Recently, there has been an increase in the popularity of automating the
process of root cause analysis. Among the most popular unsupervised methods
that deal with time series is CloudRanger (Wang et al., 2018) which is
decomposed into two steps. First, it discovers the summary causal graph
between anomalous time series using the PC algorithm (Spirtes et al., 2000)
which was introduced for non-temporal data. Then, identifies root causes
through random walk based on a transition matrix computed using the
correlation between time series. As one might expect, the main limitation of
this method is that it uses a non-temporal algorithm that does not take into
account temporal lags that might exist between two time series. In addition,
correlation does not necessarily represent the causal effect of one variable
on another. To fix these issues, Meng et al. (2020) proposed a similar method,
called MicroCause, where the causal discovery is done using the PCMCI555There
exists two versions of PCMCI, one that allows for instantaneous relations
(Runge, 2020) and one that does not (Runge et al., 2019). In the
experimentation section, we use the version that allows instantaneous
relations. (Runge et al., 2019; Runge, 2020) algorithm, an extension of the PC
algorithm for time series which infers a window causal graph. Then, MicroCause
deduces from the inferred window causal graph a summary causal graph.
Furthermore, to compute the transition matrix for the random walk, MicroCause
estimates the partial correlation between each causally related time series
given their parents in the graph. Conditioning on the parents is a sufficient
condition to eliminate all spurious correlations when there are no hidden
common causes. Note that PC and PCMCI use conditional independencies to infer
the causal graph and such methods are known to be correct when the
faithfulness condition is satisfied. Closer to our proposal, Budhathoki et al.
(2021) introduced a formal method, that we will denote as WhyMDC, to detect
the root cause of a change in a marginal distribution from non temporal data.
WhyMDC considers that a directed acyclic causal graph is given and as far as
we know, it is the first method to identify root causes by searching for
changes in causal mechanisms.

There exist other root cause identification methods which are beyond the scope
of this paper. For example, Budhathoki et al. (2022) proposed a root cause
analysis framework to detect the root cause of a point anomaly using non-
temporal structural causal models and Zhang et al. (2022) proposed a
supervised learning approach to find root causes.

##  4 ROOT CAUSE IDENTIFICATION USING ASCGLs

###  4.1 Grouping related anomalies

We first give an extension of the concept of d-seperation to ASCGL and then
show how it can be used to divide the root cause identification problem into
many independent subproblems.

A path is said to be _blocked_ by a set of vertices
$\mathcal{Z}\in\mathcal{V}$ if it contains a chain $X\rightarrow W\rightarrow
Y$ or a fork $X\leftarrow W\rightarrow Y$ and $W\in\mathcal{Z}$, or it
contains a collider $X\rightarrow W\leftarrow Y$ such that no descendant of
$W$ is in $\mathcal{Z}$. A path is said to be _active_ if it is not blocked.
Using blocked paths, the notion of d-separation is defined as follows:

######  Definition 6 (d-separation, Pearl (2000)).

Given a DAG $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and disjoint sets
$\mathcal{X},\mathcal{Y},\mathcal{Z}\subseteq\mathcal{V}$, $\mathcal{X}$ and
$\mathcal{Y}$ are _d-separated_ by $\mathcal{Z}$ if every path between a
vertex in $\mathcal{X}$ and a vertex in $\mathcal{Y}$ is blocked by
$\mathcal{Z}$. We write d-separated as
$\mathcal{X}\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{Y}\mid\mathcal{Z}$.

Note that d-separation was introduced for directed acyclic graph (DAG), so it
is directly applicable for full time graphs and window causal graphs but not
for summary causal graphs. However, it turned out that the extention to ASCGL
is simple. If there are no loops, d-separation in an ASCGL is equivalent to
the one in Definition 6. For example, in Figure 1(c), if we omit the loops
then it is obvious that
$X\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}W\mid Z$. If there are
loops, at first glance, Definition 6 might seem to fail. However, The
following proposition shows how Definition 6 can still be used for ASCGL.

######  Proposition 1.

Given an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and disjoint sets
$\mathcal{X},\mathcal{Y},\mathcal{Z}\subseteq\mathcal{V}$,
$\mathcal{X}\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{Y}\mid\mathcal{Z}$
if $\mathcal{Z}=Parents(\mathcal{X},\mathcal{G})\cup
Parents(\mathcal{Y},\mathcal{G})$ and
$\mathcal{X}\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}^{\prime}}\mathcal{Y}\mid\mathcal{Z}^{\prime}$
such that $\mathcal{G}^{\prime}$ is a DAG identical to $\mathcal{G}$ but loops
are omitted and
$\mathcal{Z}^{\prime}=Parents(\mathcal{X},\mathcal{G}^{\prime})\cup
Parents(\mathcal{Y},\mathcal{G}^{\prime})$.

###### Proof.

Consider a DAG $\mathcal{G}^{\prime}$ without loops. Suppose
$\mathcal{X},\mathcal{Y}\subseteq\mathcal{V}^{\prime}$,
$\mathcal{Z}^{x}=Parents(\mathcal{X},\mathcal{G}^{\prime})$ and
$\mathcal{Z}^{y}=Parents(\mathcal{Y},\mathcal{G}^{\prime})$ such that
$\mathcal{X}\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}^{\prime}}\mathcal{Y}\mid\mathcal{Z}^{x}\cup\mathcal{Z}^{y}$.
Now consider an identifical graph $\mathcal{G}$ but with loops. $\forall
X_{t-\gamma_{xy}},Y_{t}$ such that $\gamma_{xy}\in\mathds{N}$, the set
$\mathcal{Z}^{x,t}=\mathcal{Z}^{x}_{t-\gamma_{xy}}\cup\cdots\cup\mathcal{Z}^{x}_{t-\gamma_{xy}-\gamma_{max}}$
contains all parents of $X_{t-\gamma}$ in $\mathcal{Z}^{x}$ and none of its
decedants, thus $\mathcal{Z}^{x,t}$ blocks all active paths going into
$X_{t-\gamma_{xy}}$ that does not pass by the past $X_{t-\gamma_{xy}}$. It
follows that $\mathcal{Z}^{x,t}\cup
X_{t-\gamma_{xy}-1},\cdots,X_{t-\gamma_{xy}-\gamma_{max}}$ blocks all active
paths going into $X_{t-\gamma_{xy}}$. Similarly,
$\mathcal{Z}^{y}_{t}\cup\cdots\cup\mathcal{Z}^{y}_{t-\gamma_{max}}\cup
Y_{t-1},\cdots,Y_{t-\gamma_{max}}$ blocks all active paths going into $Y_{t}$.
Therefore,
$\mathcal{X}\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{Y}\mid\mathcal{Z}^{x}\cup\mathcal{Z}^{y}$
in $\mathcal{G}$. ∎

Note that we focused on parents and excluded ancestors to avoid separation
sets of infinite size. For example, in Figure 1(c), we can explain why
$X\not\\!\\!\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}W\mid Z$ by
looking at the compatible full-time causal graph in Figure 1(a) where
$X_{t}\not\\!\\!\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}W_{t}\mid
Z_{t-1},X_{t-1},W_{t-1}$.

Assumption 2 and Definition 3 imply that root causes of anomalous vertices
$\mathcal{A}$ is a subset of $\mathcal{A}$. So in the following, we consider
that $\mathcal{A}=\mathcal{C}\cup\mathcal{\bar{C}}$. Such that $\mathcal{C}$
represents root causes and $\mathcal{\bar{C}}$ represents non root causes.

######  Proposition 2.

Given an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and anomalous vertices
$\mathcal{A}\subseteq\mathcal{V}$ such that
$\mathcal{A}=\mathcal{C}\cup\mathcal{\bar{C}}$,
$\forall\mathcal{S}\subseteq\mathcal{V}\backslash\mathcal{A}$,
$\mathcal{C}\not\\!\\!\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{\bar{C}}\mid\mathcal{S}$.

###### Proof.

Consider anomalous vertices $\mathcal{A}=\mathcal{C}\cup\mathcal{\bar{C}}$,
such that $\mathcal{C}$ is the set of root causes of $\mathcal{\bar{C}}$. If
$\exists\mathcal{S}\subseteq\mathcal{V}\backslash\mathcal{A}$ such that
$\mathcal{C}\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{\bar{C}}\mid\mathcal{S}$
then by definition of d-separation all paths between $\mathcal{C}$ and
$\mathcal{\bar{C}}$ are blocked given $\mathcal{S}$ which means all directed
paths from $\mathcal{C}$ to $\mathcal{\bar{C}}$ are blocked given
$\mathcal{S}$. In this case, $\forall X\in\mathcal{C},\forall
Y\in\mathcal{\bar{C}}$, there exists no directed path $\pi$ from $X$ to $Y$
such that each vertex on $\pi$ belongs to $\mathcal{A}$. It follows that
$\mathcal{\bar{C}}$ is not propagated from $\mathcal{C}$ which contradicts
Assumption 2. Hence it must be the case that
$\forall\mathcal{S}\subseteq\mathcal{V}\backslash\mathcal{A}$,
$\mathcal{C}\not\\!\\!\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{\bar{C}}\mid\mathcal{S}$.
∎

######  Definition 7 (Linked anomalous graph).

Given an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and a set of anomalous
vertices $\mathcal{A}\subset\mathcal{V}$.
$\mathcal{L}=\\{\mathcal{L}^{1},\cdots,\mathcal{L}^{m}\\}$ is a set of linked
anomalous graphs if $\forall i\in\\{1,\cdots,m\\}$
$\mathcal{L}^{i}=(\mathcal{A}^{i},\mathcal{E}^{i})$ is a subgraph of
$\mathcal{G}$ such that $\mathcal{A}^{i}\subset\mathcal{A}$ and there exists a
set of vertices $\mathcal{S}\subset\mathcal{V}\backslash\mathcal{A}$ such that
$\mathcal{A}^{i}\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{A}\backslash\mathcal{A}^{i}\mid\mathcal{S}$.

For example, consider the ASCGL in Figure 2: $B,C,D,W,X,Y,Z$ are anomalous
vertices and $A$ is a normal vertex. Since
$B,C,D\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}W,X,Y,Z\mid A$ and
$B,C,D$ and $W,X,Y,Z$ are respectively d-connected given $A$, then $B,C,D$ and
$W,X,Y,Z$ form two linked anomalous graphs.

######  Proposition 3.

Given an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$ if the set of linked
anomalous graphs is
$\mathcal{L}=\\{\mathcal{L}^{1},\cdots,\mathcal{L}^{m}\\}$, then $\forall
i,j\in\\{1,\cdots,m\\},\mathcal{L}^{i}\cap\mathcal{L}^{j}=\emptyset$.

###### Proof.

Consider two different linked anomalous graphs
$\mathcal{L}^{1}=(\mathcal{A}^{1},\mathcal{E}^{1})$ and
$\mathcal{L}^{2}=(\mathcal{A}^{2},\mathcal{E}^{2})$ such that
$\mathcal{A}^{1}\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{A}\backslash\mathcal{A}^{1}\mid\mathcal{S}^{1}$
and
$\mathcal{A}^{2}\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{A}\backslash\mathcal{A}^{2}\mid\mathcal{S}^{2}$
such that
$\mathcal{S}^{1},\mathcal{S}^{2}\subseteq\mathcal{V}\backslash\mathcal{A}$. It
follows that if $\mathcal{L}^{1}\cap\mathcal{L}^{2}\neq\emptyset$ then
$\exists X\in\mathcal{V}$ such that $X\in\mathcal{A}^{1}$ and
$X\in\mathcal{A}^{2}$. In consequence,
$\not\exists\mathcal{S}\subseteq\mathcal{V}\backslash\mathcal{A}$ such that
$\mathcal{A}^{1}~{}\not\\!\\!\rotatebox[origin={c}]{90.0}{$\models$}_{\mathcal{G}}\mathcal{A}^{2}\mid\mathcal{S}$.
Which means according to Definition 7, $\mathcal{A}^{1}$ and $\mathcal{A}^{2}$
belong to the same linked anomalous graph. ∎

######  Proposition 4.

Given an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$ if the set of linked
anomalous graphs is
$\mathcal{L}=\\{\mathcal{L}_{1},\cdots,\mathcal{L}_{m}\\}$, then $\forall
i,j\in\\{1,\cdots,m\\},\mathcal{C}_{i}\cap\mathcal{C}_{j}=\emptyset$ such that
$\mathcal{C}_{i}$ is the set of root causes of $\mathcal{L}_{i}$ and
$\mathcal{C}_{j}$ is the set of root causes of $\mathcal{L}_{j}$.

###### Proof.

This follows from Proposition 3. ∎

Propositions 3 and 4 suggest that linked anomalous graphs are modular with
respect to each other, which implies that the set of root causes of each
linked anomalous graph can be identified independently of the rest of the
anomalies in the graph.

$X$$Z$$W$$Y$$A$$B$$C$$D$ Figure 2: An ASCGL with two linked anomalous graphs.
White vertices represents normal vertices and orange vertices represents
anomalous vertices.

Next, we will show how to detect a subset of the root causes uniquely by
looking at the graph and the time of the first appearance of anomalies on each
vertex.

###  4.2 Identifying root causes from the graph

######  Definition 8 (Sub-root vertex).

A sub-root vertex is root vertex in a linked anomalous graph.

######  Definition 9 (Time defying vertex).

Consider a linked anomalous graph
$\mathcal{L}^{i}=(\mathcal{A}^{i},\mathcal{E}^{i})$. $Y$ is a time defying
vertex if and only if $\forall X\in Parents(Y,\mathcal{L}^{i})$ the appearance
time of the anomaly on $Y$ precedes the appearance time of the anomaly on $X$.

######  Proposition 5.

Given a linked anomalous graph
$\mathcal{L}^{i}=(\mathcal{A}^{i},\mathcal{E}^{i})$, its set of sub-root
vertices $\mathcal{R}^{i}$, and its set of time defying vertices
$\mathcal{T}^{i}$, then
$\mathcal{R}^{i}\cup\mathcal{T}^{i}\subseteq\mathcal{C}^{i}$ such that
$\mathcal{C}^{i}$ is the true set of root causes in $\mathcal{L}^{i}$.

###### Proof.

Consider a linked anomalous graph
$\mathcal{L}^{i}=(\mathcal{A}^{i},\mathcal{E}^{i})$. 1) By Definition 8, a
sub-root vertex in $\mathcal{L}^{i}$ does not have any anomalous parent, it
follows from Assumption 2 and Propositions 3,4 that the anomaly in this vertex
cannot be propagated from other vertices, which implies that it was itself
directly affected by an external intervention which means it is a root cause.
2) Consider two anomalous vertices $X,Y\in\mathcal{A}^{i}$ such that
$Parents(Y,\mathcal{L}^{i})=\\{X\\}$. According to Assumption 2, the anomaly
in $X$ which appeared at time $t$ was propagated to $Y$ according to the lag
$\gamma\geq 0$ between $X$ and $Y$ of the SCM. It follows that if the
appearance time of the anomaly on $Y$ is $t^{\prime}$ such that
$t^{\prime}<t\leq t+\gamma$ then the anomaly on $Y$ was not propagated from
$X$, which implies that $Y$ was itself directly affected by an external
intervention which means it is a root cause. Using induction, suppose this is
true for $Parents(Y,\mathcal{L}^{i})=\\{X^{1},\cdots,X^{p}\\}$. If
$Parents(Y,\mathcal{L}^{i})=\\{X^{1},\cdots,X^{p+1}\\}$ such that $t^{\prime}$
precedes the appearance time of the anomalies on $\\{X^{1},\cdots,X^{p}\\}$
then if $t^{\prime}$ does not precede the appearance time on $X^{p+1}$, then
the anomaly on $Y$ could have been propagated from $X^{p+1}$, otherwise, we
conclude that $Y$ was directly affected by an external intervention which
means it is a root cause. ∎

Proposition 5 states that if a vertex is a sub-root or a time-defying vertex
in a linked anomalous graph then it belongs to the set of root causes.
However, this does not mean that every element in the set of root causes is
necessarily a sub-root or a time-defying vertex. Therefore, there might be a
vertex
$X\in\mathcal{A}^{i}\backslash\\{\mathcal{R}^{i}\cup\mathcal{T}^{i}\\}$, such
that $X\in\mathcal{C}^{i}$.

###  4.3 Identifying root causes from data

To find root causes that are neither sub-roots nor time-defying vertices, we
search for changes in the causal mechanisms. In some cases, we can find these
changes, for each $X\rightarrow Y$ in a given linked anomalous graph, by
estimating the total effect666The _do_ operator represents an external
intervention. (Pearl, 2000) of $X$ on $Y$ in the normal regime, defined as

$\displaystyle TE^{{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}=$ $\displaystyle
E_{N}[Y_{t}|do(X_{t-\gamma_{xy}}=x)]$ (2) $\displaystyle-
E_{N}[Y_{t}|do(X_{t-\gamma_{xy}}=x^{\prime})]$

and the total effect of $X$ on $Y$ in the anomalous regime, defined as

$\displaystyle TE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}=$
$\displaystyle E_{\bar{N}}[Y_{t}|do(X_{t-\gamma_{xy}}=x)]$ (3) $\displaystyle-
E_{\bar{N}}[Y_{t}|do(X_{t-\gamma_{xy}}=x^{\prime})],$

where $E_{N}$ and $E_{\bar{N}}$ are respectivly the expectations in the normal
regime and in the anomalous regime and the temporal lag $\gamma_{xy}$ is
between $\bar{\gamma}_{xy}$ and $\gamma_{max}$ such that $\bar{\gamma}_{xy}$
can be found by the substracting the time of appearance of anomalies on $Y$
from the time of appearance of anomalies on $X$. If there is no directed path
from $X$ to $Y$ (other than $X\rightarrow Y$) or if directed paths between $X$
and $Y$ exist but we know that all vertices (other than $X$ and $Y$) on these
paths cannot be root causes, then if $TE^{{N}}_{X_{t-\gamma_{xy}}\rightarrow
Y_{t}}\neq TE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}$ we can conclude
that there is a change in the causal mechanism of $Y$ provoked by an external
intervention on $Y$.

In order to estimate $TE^{N}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}$ and
$TE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}$ from observational data we
need to eliminate the _do_ from the total effect expression. This can be
achieved, when there is no hidden common causes, by the back-door criterion
(Pearl, 2000) which searches for an adjustment set of vertices, called back-
door set, that eliminates all spurious correlations between $X$ and $Y$.
However, the back-door criterion cannot directly be applied to ASCGL because
of loops. In the following we present an adjustment set for identifying total
effects777A similar result was presented in (Eichler and Didelez, 2007) for
summary causal graphs assuming there is no instantaneous relations and
allowing for sets of infinite size. from ASCGLs:

######  Definition 10 (An adjustment set for total effects in an ASCGL).

Consider an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$, a maximal lag
$\gamma_{max}$, two vertices $X$ and $Y$ such that $X\rightarrow Y$ in
$\mathcal{G}$, and the temporal lag $\gamma_{xy}$ between $X$ and $Y$. The
adjustment set for identifying the total effect relative to
($X_{t-\gamma_{xy}},Y_{t}$) is
$\mathcal{B}_{t-\gamma_{xy}}\cup\cdots\cup\mathcal{B}_{t-\gamma_{xy}-\gamma_{max}}\cup\mathcal{X}$
such that:  1. $\mathcal{B}=Parents(X,\mathcal{G})\backslash\\{X\\}$;  2.
$\mathcal{X}=\\{X_{t-\gamma_{xy}-1},\cdots,X_{t-\gamma_{xy}-\gamma_{max}}\\}$if
there exists a loop on $X$ in $\mathcal{G}$, otherwise
$\mathcal{X}=\\{\emptyset\\}$.

######  Proposition 6.

Given an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and a maximal lag
$\gamma_{max}$, if
$\mathcal{B}_{t-\gamma_{xy}}\cup\cdots\cup\mathcal{B}_{t-\gamma_{xy}-\gamma_{max}}\cup\mathcal{X}$
satisfies Definition 10 in $\mathcal{G}$ relative to
($X_{t-\gamma_{xy}},Y_{t}$) then
$\mathcal{B}_{t-\gamma_{xy}}\cup\cdots\cup\mathcal{B}_{t-\gamma_{xy}-\gamma_{max}}\cup\mathcal{X}$
blocks all activated paths between $X_{t-\gamma_{xy}}$ and $Y_{t}$ going into
$X_{t-\gamma_{xy}}$ in every window causal graph associated with
$\mathcal{G}$.

###### Proof Sketch.

If $X$ and $Y$ have no loops in $\mathcal{G}$, then
$\mathcal{B}_{t-\gamma_{xy}}\cup\cdots\cup\mathcal{B}_{t-\gamma_{xy}-\gamma_{max}}$
is sufficient to block all paths between $X_{t-\gamma_{xy}}$ and $Y_{t}$ going
into $X_{t-\gamma_{xy}}$ since all possible parents of $X_{t-\gamma_{xy}}$ are
in
$\mathcal{B}_{t-\gamma_{xy}}\cup\cdots\cup\mathcal{B}_{t-\gamma_{xy}-\gamma_{max}}$.
Given that $\mathcal{G}$ is acyclic,
$\mathcal{B}_{t-\gamma_{xy}}\cup\cdots\cup\mathcal{B}_{t-\gamma_{xy}-\gamma_{max}}$
cannot block any directed path nor create any new activated path between
$X_{t-\gamma_{xy}}$ and $Y_{t}$ as there cannot be any descendant of
$X_{t-\gamma_{xy}}$ in
$\mathcal{B}_{t-\gamma_{xy}}\cup\cdots\cup\mathcal{B}_{t-\gamma_{xy}-\gamma_{max}}$.
If $X$ and $Y$ have loops then adjusting on $\mathcal{B}$ up to $\gamma_{max}$
cannot block all back-door paths because there will always be an activated
path passing by $\mathcal{B}_{t-\gamma_{xy}-\gamma_{max}-i}$ such that $i>0$.
The only way to block this path is to add the past of $X_{t-\gamma_{xy}}$ to
the adjustment set. ∎

When there exist directed path between $X$ and $Y$ the total effect would no
longer be reliable to detect external interventions. For example, in Figure 2,
if there’s an external intervention on $B$ then
$TE^{N}_{C_{t-\gamma_{cd}}\rightarrow D_{t}}\neq
TE^{\bar{N}}_{C_{t-\gamma_{cd}}\rightarrow D_{t}}$ due to the change in the
causal mechanism of the mediator $B$ of $C$ and $D$. To avoid such cases, for
each $X\rightarrow Y$ in a given linked anomalous graph, we need to estimate
the direct effect (Pearl, 2000) of $X$ on $Y$ in the normal regime, defined as

$\displaystyle DE^{{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}=$ $\displaystyle
E_{N}[Y_{t}|do(X_{t-\gamma_{xy}}=x,\mathcal{W}=\textit{w})]$ (4)
$\displaystyle-
E_{N}[Y_{t}|do(X_{t-\gamma_{xy}}=x^{\prime},\mathcal{W}=\textit{w})]$

and the direct effect of $X$ on $Y$ in the anomalous regime, defined as

$\displaystyle DE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}=$
$\displaystyle
E_{\bar{N}}[Y_{t}|do(X_{t-\gamma_{xy}}=x,\mathcal{W}=\textit{w})]$ (5)
$\displaystyle-
E_{\bar{N}}[Y_{t}|do(X_{t-\gamma_{xy}}=x^{\prime},\mathcal{W}=\textit{w})],$

where
$\mathcal{W}=\mathcal{V}_{t-\gamma_{max}}\cup\cdots\cup\mathcal{V}_{t-\gamma_{xy}}\backslash\\{X_{t-\gamma_{xy}}\\}\cup\cdots\cup\mathcal{V}_{t}\backslash\\{Y_{t}\\}$.

Assuming linearity, the do from the direct effect expression can be eliminated
(i.e., direct effect can be identifyed) using any adjustment set given by the
single-door criterion (Pearl, 2000) which is not applicable in ASCGLs. In the
following, we present an adjustment set for direct effects in ASCGLs:

######  Definition 11 (An adjustment set for direct effects in an ASCGL).

Consider an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$, a maximal lag
$\gamma_{max}$, two vertices $X$ and $Y$ such that $X\rightarrow Y$ in
$\mathcal{G}$, and the temporal lag $\gamma_{xy}$ between $X$ and $Y$. An
adjustment set for identifying the direct effect relative to
($X_{t-\gamma_{xy}},Y_{t}$) is
$\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}\cup\mathcal{X}\cup\mathcal{Y}$
such that:  1. $\mathcal{B}=Parents(Y,\mathcal{G})\backslash\\{X,Y\\}$;  2.
$\mathcal{X}=\\{X_{t},\cdots,X_{t-\gamma_{max}}\\}\backslash\\{X_{t-\gamma_{xy}}\\}$and
$\mathcal{Y}=\\{Y_{t-1},\cdots,Y_{t-\gamma_{max}}\\}$ if there exists a loop
on $Y$ in $\mathcal{G}$, otherwise $\mathcal{Y}=\\{\emptyset\\}$.

######  Proposition 7.

Given an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and a maximal lag
$\gamma_{max}$, if
$\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}\cup\mathcal{X}\cup\mathcal{Y}$
satisfies Definition 11 in $\mathcal{G}$ relative to
($X_{t-\gamma_{xy}},Y_{t}$) then
$\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}\cup\mathcal{X}\cup\mathcal{Y}$
blocks all activated paths between $X_{t-\gamma_{xy}}$ and $Y_{t}$ in every
window causal graph associated with $\mathcal{G}$ except the direct path
$X_{t-\gamma_{xy}}\rightarrow Y_{t}$.

###### Proof Sketch.

If $X$ and $Y$ have no loops in $\mathcal{G}$, then
$\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}$ is sufficient to
block all paths from $X_{t-\gamma_{xy}}$ to $Y_{t}$ except
$X_{t-\gamma_{xy}}\rightarrow Y_{t}$ since and all possible parents of $Y_{t}$
are in $\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}$. Given that
$\mathcal{G}$ is acyclic
$\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}$ cannot create any
new activated path that is not blocked by
$\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}$ since no vertex in
$\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}$ is a descendant of
$Y_{t}$. If $X$ and $Y$ have loops then adjusting on $\mathcal{B}$ up to
$\gamma_{max}$ cannot block all paths because there can be an activated path
passing by the future of $X_{t-\gamma_{xy}}$ or by the past of $Y_{t}$. The
only way to block these types of paths is to add $\mathcal{X}\cup\mathcal{Y}$
to the adjustment set. ∎

###  4.4 An algorithm for root cause identification

Here, we describe our main method called EasyRCA888Code available at
https://github.com/ckassaad/
EasyRCA, in which the pseudocode is provided in Algorithm 1. The algorithm
starts by finding linked anomalous graphs (line 1). Then for each linked
anomalous graph, it searches for the sub-roots and time-defying vertices (line
3). Finally, it searches for the rest of the root causes by comparing direct
effects in the normal regime with direct effects in the anomalous regime
(lines 4-16). The conditions in lines 12 and 14 need the minimality condition
because if $X$ and $Y$ are statistically independent given the adjustment set
in the normal regime then an intervention on $Y$ might not imply any change to
the statistical distribution thus one cannot conclude on the presence of
interventions. The for-loop in line 2 can be parallelized since as showed in
Proposition 3 and 4, linked anomalous graphs are modular.

######  Theorem 1.

Given an ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$, a set anomalous
vertices $\mathcal{A}\subseteq\mathcal{V}$, the distribution of the time
series in the normal regime $\mathcal{{N}}$ and in the anomalous regime
$\mathcal{\bar{N}}$, and the maximal lag between a cause and an effect
$\gamma_{max}$, under Assumption 2 and the minimality condition, EasyRCA is
capable of identifying the set of root cause $\mathcal{C}$ of $\mathcal{A}$.

###### Proof Sketch.

It follows from Propositions 3, 4, 5, 7. ∎

Note that we can also distinguish between parametric and structural
interventions. Given that $DE^{N}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}\neq
DE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}$, if
$DE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}=0$, we conclude that the
intervention on $Y$ is structural, otherwise, we conclude that it is
parametric.

Algorithm 1 EasyRCA

0: ASCGL $\mathcal{G}=(\mathcal{V},\mathcal{E})$, distribution of the time
series in the normal regime $\mathcal{{N}}$ and in the anomalous regime
$\mathcal{\bar{N}}$, maximal lag $\gamma_{max}$, Anomalies $\mathcal{A}$

1: $\mathcal{L}^{1},\cdots,\mathcal{L}^{m}=$ list of linked anomalous graphs
as in Definition 7

2: for $i\in\\{1,\cdots,m\\}$ do

3: Identify sub-root vertices $\mathcal{R}^{i}$ and time defying vertices
$\mathcal{T}^{i}$ using Definition 8 and 9

4: Let $\mathcal{D}^{i}=[]$

5: Let $\mathcal{A}^{i}$ be the set of vertices in $\mathcal{L}^{i}$

6: for $Y$ in
$\mathcal{A}^{i}\backslash\\{\mathcal{R}^{i}\cup\mathcal{T}^{i}\\}$ do

7: for $X$ in $Parents(Y,\mathcal{G})$ do

8: $\bar{\gamma}_{xy}$: anomaly lag between $X$ and $Y$

9: for $\gamma_{xy}$ in $\\{\bar{\gamma}_{xy},\cdots,\gamma_{max}\\}$ do

10: Identify
$\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}\cup\mathcal{X}\cup\mathcal{Y}$
using Definition 11

11: Estimate $DE^{N}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}$

12: if $DE^{N}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}\neq 0$ then

13: Estimate $DE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}$

14: if $DE^{N}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}\neq
DE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}$ then

15: $\mathcal{D}^{i}=[\mathcal{D}^{i},Y]$

16: Break

17: Return $\mathcal{R}$, $\mathcal{T}$, $\mathcal{D}$

##  5 EXPERIMENTS

We propose first an extensive analysis on simulated data, generated from
random causal graphs; then we perform an analysis on a real word dataset.

###  5.1 Experimental Setup

In practice, to test if $DE^{N}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}\neq
DE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}$, we fit 11 multiple linear
regressions:

$Y_{t}=\hat{a}_{x}X_{t-\gamma_{xy}}+\sum_{B_{t-\gamma_{by}}\in\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}\cup\mathcal{X}\cup\mathcal{Y}}\hat{a}_{b}B_{t-\gamma_{by}}+\epsilon^{y}_{t},$

such that
$\mathcal{B}_{t}\cup\cdots\cup\mathcal{B}_{t-\gamma_{max}}\cup\mathcal{X}\cup\mathcal{Y}$
is identified using Definition 11. One of them fitted on the anomalous data
and 10 fitted on different chunks of the normal data. Then using the Grubbs-
test (Grubbs, 1950) we check if the coefficient $\hat{a}_{x}$ of the anomalous
data is significantly different than the 10 others $\hat{a}_{x}$ from the
normal data. To test if $TE^{\bar{N}}_{X_{t-\gamma_{xy}}\rightarrow Y_{t}}=0$,
we use a t-test on the coefficient $\hat{a}_{x}$ of the anomalous data.

Baselines: We compare EasyRCA with three other methods999We implemented
CloudRanger and MicroCause and we adapted WhyMDC based on the DoWhy package.:
CloudRanger, MicroCause and a naive adaptation of WhyMDC to time series for
which we provide a window causal graph. Since CloudRanger and MicroCause try
to solve a harder problem compared to EasyRCA by inferring the summary causal
graph from anomalous data (while EasyRCA considers that the summary causal
graph is given), we also consider another version of EasyRCA, denoted as
EasyRCA∗, where we suppose that the summary causal graph is not given. In the
first step of EasyRCA∗, we infer the window causal graph from normal data
using PCMCI (Runge, 2020), the same causal discovery algorithm used by
MicroCause, and then we deduce the summary causal graph from it. Note that the
summary causal graph obtained from the window graph that is inferred by PCMCI
can be cyclic even if the true summary causal graph is acyclic (due to
estimation errors). In such cases, we consider that EasyRCA∗ does not identify
any root cause.

Hyper-parameters: For EasyRCA, EasyRCA∗ and MicroCause, we set the maximal lag
$\gamma_{max}$ to $3$ and for all methods (even though the true $\gamma_{max}$
is smaller in our simulation study), we set the significance threshold to
$0.01$. For EasyRCA∗, CloudRanger and MicroCause we use a Fisher-z-test
(Kalisch and Bühlmann, 2007), which is commonly used for causal discovery when
linearity and gaussianity are satisfied. Furthermore, for CloudRanger and
MicroCause we set the walk length to $1000$ and the backward step threshold to
$0.1$. Lastly, other hyper-parameters in WhyMDC were set to the default values
in the DoWhy package.

Evaluation: To assess the quality of identifying root causes, we use the
F1-score. Since by construction EasyRCA identifies sub-root and time defying
vertices as root causes we do not evaluate their detection. This ensures a
fair comparison with other methods.

###  5.2 Simulated Data

For simulated data, we start by randomly generating $30$ different ASCGLs such
that each graph contains $6$ vertices, has a maximal degree between $4$ and
$5$, and has one root vertex. We consider that all lags in the window causal
graph associated with any of the generated ASCGL are equal to $1$. So the
generative process (the SCM) is the following:

$Y_{t}=\sum_{X_{t-1}\in
Parents(Y_{t},\mathcal{G}_{w})}aX_{t-1}+0.1\xi^{y}_{t}$

where $a\sim U\\{0.1,1\\}$, $\xi^{y}_{t}\sim\mathcal{N}(0,1)$, $Y_{t}$ denotes
the value of the vertex at time $t$, $Parents(Y_{t})$ denotes the direct
parents of $Y_{t}$ in the window causal graph.

For each ASCGL, we choose two root causes (two vertices that will undergo an
intervention): the first is the root of the ASCGL and the second is a randomly
chosen vertex among the non-root vertices. We propagate the effect of each
intervention according to the generating process toward all the descendants of
the vertex which underwent an intervention. We set the starting time of each
intervention according to the generative process. For example, if the root
vertex $X$ has an intervention at time $t$ and a vertex $Y$ randomly selected
to undergo an intervention is a child of $X$ then the starting time of the
intervention on $Y$ is $t+1$. In general, the starting time of the
intervention does not need to respect the generative process, but we chose to
respect it to avoid any time-defying vertices which would give an advantage to
our method.

In our experiments, we consider the two types of interventions separately and
we vary the anomaly size between $100$ and $2000$. In the case of structural
interventions, values of the root vertex and values of a randomly chosen non-
root vertex in the anomalous interval are replaced by data drawn form the
distribution $Exp(2)$. In the case of parametric interventions, values of the
root are set similarly to structural interventions. Then values of the non-
root vertex are regenerated with new coefficients from $U(0.1,1)$.

Results: In Figure 3, we report the performance of each method at detecting
structural interventions with respect to the anomaly size. As one can see,
EasyRCA and EasyRCA∗ clearly outperform other methods, and their performances
increases (and their variance decreases) significantly between the anomaly of
size $100$ and the anomaly of size $2000$ reaching an F1-score of $1$ for both
EasyRCA and EasyRCA∗. The small difference in the performance of EasyRCA and
EasyRCA∗ shows that our method is robust with respect to small errors in the
ASCGL. WhyMDC and MicroCause have similar results and outperform CloudRanger.
In Figure 4, we report the performance of each method at detecting parametric
interventions with respect to the anomaly size. As before, EasyRCA and
EasyRCA∗ outperforms other methods but now the difference between EasyRCA and
EasyRCA∗ is more visible and all other methods suffer. However, it is worth
noting, that unlike CloudRanger, WhyMDC and MicroCause were able to detect the
root vertex of each graph as a root cause.

$100$$200$$500$$1000$$2000$$0$$0.5$$1$Anomaly
sizeF1EasyRCAEasyRCA∗WhyMDCMicroCauseCloudRanger Figure 3: Mean and variance
of F1-scores with respect to structural interventions over 30 graphs
containing one linked anomalous graph with one sub-root vertex and one
structural intervention. $100$$200$$500$$1000$$2000$$0$$0.5$$1$Anomaly
sizeF1EasyRCAEasyRCA∗WhyMDCMicroCauseCloudRanger Figure 4: Mean and variance
of F1-scores with respect to parametric interventions over 30 graphs
containing one linked anomalous graph with one sub-root vertex and one
parametric intervention.

###  5.3 Real Data

For real data, we consider a dataset101010The real IT monitoring data is
available at https://easyvista2015-my.sharepoint.com/personal/aait-
bachir_easyvista_com/_layouts/15/onedrive.aspx?id=%2Fpersonal%2Faait%2Dbachir%5Feasyvista%5Fcom%2FDocuments%2FLab%2FPublicData&ga=1
which consists of eight time series collected from an IT monitoring system
with a one-minute sampling rate provided by
EasyVista111111https://www.easyvista.com/fr/produits/ev-observe such that each
of these time series is considered anomalous and all collective anomalies are
considered to have the same time of appearance and of size $100$. The
corresponding ASCGL is provided in Figure 5 where PMDB represents the
extraction of some information about the messages received by the Storm
ingestion system; MDB refers to an activity of a process that orient messages
to other process with respect to different types of messages; CMB represents
the activity of extraction of metrics from messages; MB represents the
activity of insertion of data in a database; LMB reflects the updates the last
values of metrics in Cassandra; RTMB represents the activity of searching to
merge of data with information coming from the check message bolt; GSIB
represents the activity of insertion of historical status in database. ESB
represents the activity of writing data in Elasticsearch. According to
EasyVista’s system experts, PMDB and ESB are expected to be the root causes of
these anomalies.

EasyRCA inferred $3$ roots causes, PMDB as a root vertex, in addition to RTMB
and ESB as structural interventions. EasyRCA∗ inferred $5$ roots causes, PMDB,
GSIB and MB as root vertices, in addition to RTMB and ESB as structural
interventions. MicroCause inferred that PMDB and MB are the root causes of the
anomalies. CloudRanger inferred that GSIB and MDB are the root causes. We did
not apply WhyMDC in this real world application because the true window causal
graph is unknown. In terms of the trade-off between false positives and false
negatives, EasyRCA gives the best result.

PMDBMDBCMBMBLMBRTMBGSIBESB Figure 5: ASCGL of the normal regime of an IT
monitoring system. All vertices are anomalous in the anomalous regime.
According to EasyVista’s system experts, PMDB and ESB are expected to be the
root causes of these anomalies.

##  6 Conclusion

We adressed the problem of identifying root causes of collective anomalies
using observational time series and an ASCGL of the normal regime of a given
system. We showed that the problem can be divided into many independent
subproblems and that all root causes can be identified using the graph and the
data. For future work, it would be interesting to extend this method for
cyclic summary causal graphs, for nonlinear SCMs and to allow for hidden
common causes.

#### Acknowledgements

We thank Ali Aït-Bachir, Christophe de Bignicourt and Hosein Mohanna from
EasyVista for providing the IT monitoring data along with the underlying
causal graph and for localizing anomalies in the data.

## References

  * Assaad et al. (2021) Charles K. Assaad, Emilie Devijver, Eric Gaussier, and Ali Ait-Bachir.  A mixed noise and constraint-based approach to causal inference in time series.  In _Machine Learning and Knowledge Discovery in Databases. Research Track_ , pages 453–468, Cham, 2021. Springer International Publishing.  ISBN 978-3-030-86486-6.
  * Assaad et al. (2022a) Charles K. Assaad, Emilie Devijver, and Eric Gaussier.  Survey and evaluation of causal discovery methods for time series.  _J. Artif. Int. Res._ , 73, apr 2022a.  doi: 10.1613/jair.1.13428.
  * Assaad et al. (2022b) Charles K. Assaad, Emilie Devijver, and Eric Gaussier.  Entropy-based discovery of summary causal graphs in time series.  _Entropy_ , 24(8), 2022b.  ISSN 1099-4300.  doi: 10.3390/e24081156.
  * Budhathoki et al. (2021) Kailash Budhathoki, Dominik Janzing, Patrick Bloebaum, and Hoiyi Ng.  Why did the distribution change?  In Arindam Banerjee and Kenji Fukumizu, editors, _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , pages 1666–1674. PMLR, 13–15 Apr 2021.
  * Budhathoki et al. (2022) Kailash Budhathoki, Lenon Minorics, Patrick Bloebaum, and Dominik Janzing.  Causal structure-based root cause analysis of outliers.  In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, _Proceedings of the 39th International Conference on Machine Learning_ , volume 162 of _Proceedings of Machine Learning Research_ , pages 2357–2369. PMLR, 17–23 Jul 2022.
  * Chandola et al. (2009) Varun Chandola, Arindam Banerjee, and Vipin Kumar.  Anomaly detection: A survey.  _ACM Comput. Surv._ , 41:15:1–15:58, 2009.
  * Eberhardt and Scheines (2007) Frederick Eberhardt and Richard Scheines.  Interventions and causal inference.  _Philosophy of Science_ , 74(5):981–995, 2007\.  ISSN 00318248, 1539767X.
  * Eichler and Didelez (2007) Michael Eichler and Vanessa Didelez.  Causal reasoning in graphical time series models.  In _Proceedings of the Twenty-Third Conference on Uncertainty in Artificial Intelligence_ , UAI’07, page 109–116, Arlington, Virginia, USA, 2007\. AUAI Press.  ISBN 0974903930.
  * Glymour et al. (2019) Clark Glymour, Kun Zhang, and Peter Spirtes.  Review of causal discovery methods based on graphical models.  _Frontiers in Genetics_ , 10, 2019.  ISSN 1664-8021.  doi: 10.3389/fgene.2019.00524.
  * Grubbs (1950) Frank E. Grubbs.  Sample Criteria for Testing Outlying Observations.  _The Annals of Mathematical Statistics_ , 21(1):27 – 58, 1950.  doi: 10.1214/aoms/1177729885.
  * Kalisch and Bühlmann (2007) Markus Kalisch and Peter Bühlmann.  Estimating high-dimensional directed acyclic graphs with the pc-algorithm.  _J. Mach. Learn. Res._ , 8:613–636, may 2007.  ISSN 1532-4435.
  * Malinsky and Danks (2018) Daniel Malinsky and David Danks.  Causal discovery algorithms: A practical guide.  _Philosophy Compass_ , 13(1):e12470, 2018.  doi: https://doi.org/10.1111/phc3.12470.
  * Meng et al. (2020) Yuan Meng, Shenglin Zhang, Yongqian Sun, Ruru Zhang, Zhilong Hu, Yiyin Zhang, Chenyang Jia, Zhaogang Wang, and Dan Pei.  Localizing failure root causes in a microservice through causality inference.  In _2020 IEEE/ACM 28th International Symposium on Quality of Service (IWQoS)_ , pages 1–10, 2020.  doi: 10.1109/IWQoS49365.2020.9213058.
  * Pearl (2000) Judea Pearl.  _Causality: Models, Reasoning, and Inference_.  Cambridge University Press, New York, NY, USA, 2000.  ISBN 0-521-77362-8.
  * Peters et al. (2013) Jonas Peters, D. Janzing, and B. Schölkopf.  Causal inference on time series using restricted structural equation models.  In _Advances in Neural Information Processing Systems 26_ , pages 154–162, 2013.
  * Runge (2020) Jakob Runge.  Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets.  In Jonas Peters and David Sontag, editors, _Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI)_ , volume 124 of _Proceedings of Machine Learning Research_ , pages 1388–1397. PMLR, 03–06 Aug 2020.
  * Runge et al. (2019) Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic.  Detecting and quantifying causal associations in large nonlinear time series datasets.  _Science Advances_ , 5(11), 2019.  doi: 10.1126/sciadv.aau4996.
  * Spirtes et al. (2000) Peter Spirtes, Clark Glymour, and Richard Scheines.  _Causation, Prediction, and Search_.  MIT press, 2nd edition, 2000.
  * Wang et al. (2018) Ping Wang, Jingmin Xu, Meng Ma, Weilan Lin, Disheng Pan, Yuan Wang, and Pengfei Chen.  Cloudranger: Root cause identification for cloud native systems.  In _2018 18th IEEE/ACM International Symposium on Cluster, Cloud and Grid Computing (CCGRID)_ , pages 492–502, 2018.  doi: 10.1109/CCGRID.2018.00076.
  * Zhang et al. (2022) Chaoli Zhang, Zhiqiang Zhou, Yingying Zhang, Linxiao Yang, Kai He, Qingsong Wen, and Liang Sun.  Netrca: An effective network fault cause localization algorithm, 2022\.