Progression, Detection and Remission: Evolution of Chronic Myeloid Leukemia
using a three-stage probabilistic model
Sonjoy Pan †, Siddhartha P. Chakrabarty ∗ and Soumyendu Raha ¶
† Department of Mathematics, Swami Vivekananda University, Kolkata-700121,
India

∗ Department of Mathematics, Indian Institute of Technology Guwahati,
Guwahati-781039, India

¶ Department of Computational and Data Sciences, Indian Institute of Science,
Bangalore-560012, India

†<EMAIL_ADDRESS>∗<EMAIL_ADDRESS>and ¶<EMAIL_ADDRESS>
###### Abstract

We present a three-stage probabilistic model for the progression of Chronic
Myeloid Leukemia (CML), as manifested by the leukemic stem cells, progenitor
cells and mature leukemic cells. This progression is captured through the
process of cell division and cell mutation, with probabilities of occurrence
being assigned to both of them. The key contributions of this study include,
the determination of the expected number of the leukemic stem cells,
progenitor cells, mature leukemic cells, as well as total number of these
cells (in terms of probabilities, and contingent on the initial cell count),
expected time to reach a threshold level of total and injurious leukemic
cells, as well as the critical time when the disease changes its phases, the
probability of extinction of CML, and the dynamics of CML evolution consequent
to primary therapy. Finally, the results obtained are demonstrated with
numerical illustrations.

Keywords: CML; Disease progression; Extinction probability; Primary therapy

##  1 Introduction

The occurrence of leukemia is triggered by a disruption of the highly
regulated and complex process of hematopoiesis, resulting from the generation
of a handful of mutated blood cells [1]. Leukemia can be acute or chronic and
myeloid or lymphocytic, contingent on the maturity state and type of cells,
respectively. Chronic Myeloid Leukemia (CML), a commonly occurring type of
leukemia, is typically attributed to the translocation between chromosomes 9
and 22, resulting in a longer chromosome 9, and a shorter chromosome 22, than
the typical chromosomes [2]. CML, which is characterized by elevated levels of
white blood cells, accounts for a fifth of the total leukemia incidences, of
which 90 percent cases are associated with the Philadelphia (Ph) chromosome,
resulting from the translocation of chromosomes 9 and 22. The typical
progression of CML occurs in three stages [3, 4], namely, the chronic phase (a
period of relatively low cell count, accompanied by a high level of cellular
differentiation), the accelerated phase (during which higher cell levels is
observed, along with decline in the extent of differentiation), and the blast
phase (where an effectively uncontrolled cell growth, but low extent of
differentiation, is observed). In absence of intervention, the disease can
progress to the accelerated phase over a span of 7-10 years [4], which then
progresses to the blast phase.

Michor et al. [5], in a basic model, proposed to capture the dynamics and
progression of CML, and considered an ordinary differential equation (ODE)
driven setup for normal and leukemic (non-resistant and resistant) cells, for
a four-compartment model of stem cells, progenitor cells, differentiated cells
and terminally differentiated cells. The model evaluated the success of
imatinib during the molecularly targeted therapy in case of CML, and concluded
that, while imatinib achieves success in case of inhibition of differentiated
leukemic cells, it does not cause a decline in the levels of leukemic stem
cells. Further, the authors also make an estimation of the likelihood of
development of resistance to imatinib, as well as the time-point at which the
resistance is detected. It is believed that the proliferation of cancerous
stem cells, and progenitor cells brings forth, and abets the cellular growth
in CML [6]. In this paper [6], a stochastic model is considered, with
primitive and proliferating, as well as quiescent CML cell populations, to
explain the clinically observed data exhibiting, bi-phasic, one-phasic and
reverse bi-phasic decline pattern, as also the emanation of quiescent stem
cells and treatment resistance. An agent based model (ABM) to capture the
“age-structure” of CML stem cells, is proposed in [7] through the usage of
normal hematopoiesis, leukemia genesis (for normal and leukemia cells) and
leukemia treatment (using imatinib). A mathematical model considered in [8],
included the leukemic and normal cells, both stem-like and differentiated,
leading to a four-dimensional coupled delay differential equation (DDE)
system. An ODE driven model to analyze the dynamics of imatinib therapy is
developed in [9], with the focus being on patients exhibiting a non-monotonic
BCR-ABL ratio, during the course of the therapy, by considering mature cells
(quiescent, cycling, progenitor and mature), along with the concentration of
immune cells. A simplified version of the model in [9] was considered in [10],
based on the argument that a three-compartment model (two for leukemic cells
and one for immune cells) can capture the dynamics of the CML-immune response.

Population of hematopoietic stem cells and differentiated cells, both normal
and cancerous, were considered in a compartment model in [11], in the paradigm
of an optimal control problem, with the objective of minimizing the treatment
cost and load of the cancerous cells. In another control theoretic problem
[12], a two population model, involving short-term hematopoietic stem cells
and differentiated cells, possessing the ability to self-renew, and lacking
the same, respectively, is considered, in presence of imatinib. In a recent
work on optimization of combination therapy of leukemia, of imatinib and
interferon-$\alpha$, populations of CML cells and effector T cells cytotoxic
to CML, are incorporated into the model [13].

A stochastic approach to address the problem of acquired drug resistance to
imatinib, and the consequent disease progression was considered in [14],
wherein a multi-type branching process was used to encapsulate the evolution
of several clones of CML cells, resistant to imatinib. The model for
hematopoiesis in CML, emphasized the stochastic aspects applicable in case of
abnormalities in hematopoiesis, that is typical in cases of blood cancers. To
this end, stochastic compartments and stochastic competition were taken into
account in the hybrid discrete-continuous modeling framework. A causal
Bayesian network model for explaining the progression of CML from the chronic
phase to blast crisis is presented in [4]. In a recent article, the authors
presented deterministic and stochastic models for a triple treatment scenario
of Wnt/$\beta$-catenin signaling, combined with Tyrosine Kinase Inhibitors
(TKI) and IFN-$\alpha$ [15].

We also briefly review here, some of the literature on usage of branching
process in epidemiological and cancer modeling. An epidemiological model for
determination of the number of generations until extinction (for a branching
process), in the paradigm of infectious disease modeling was discussed in
[16]. Accordingly, several processes were used to obtain the generation
distributions, and reconciled in terms of application, to several infectious
disease settings. In [17], the modeling was carried out for a population of
cancerous cells, exhibiting exponential growth, progressing to treatment
resistance, as a result of one mutation, followed by development of disease
phenotype, after subsequent mutations. From a stochastic perspective, the
distribution of the arrival time of the $k$-th mutation, as well as the growth
progression of the count of type $k$ cells, was established. A random fitness
increments approach, to model the progression of tumors was studied in [18]. A
multi-type branching process was analyzed for accumulated mutations in cells,
chosen from some distribution. An ovarian cancer model, with the objective of
ascertaining an opportune window for screening was discussed in [19]. The
model is driven by a branching process, to describe the growth and progression
of three cellular populations, namely, primary, peritoneal, and metastatic,
with the model parameters being reconciled with the clinical data.

From a therapeutic perspective, TKIs, such as imatinib can have curative
impact on CML, albeit the longer periods of treatment being required to
achieve a desired result [20]. A mathematical model focused on the combination
treatment of imatinib (for molecular targeted therapy) and interferon alfa-2a
(IFN-$\alpha$) (immunotherapy) presented in [21] incorporated time-varying
delays in treatment terms, and demonstrated how the dual therapy led to
improved outcomes in case of CML patients.

In this paper, we study on the dynamics of CML evolution, considering various
stages of disease progression and different types of cells involved in
leukemic cancer. We also determine the expected time to reach a possible state
of CML and the chance of disease extinction. The rest of the paper is
organized as follows. In Section 2, the probabilistic modeling of CML
progression is described, using the theory of multi-type branching processes.
Some analytical results on cells count, estimated time, and associated
probabilities for occurrence of various kinds of leukemic cells and subsequent
conditions are established in Section 3. Section 4 is devoted to several
numerical simulations. Finally, the conclusion is summarized in Section 5.

##  2 Model

In CML patients, the leukemic cells develop through various cell
differentiation and cell mutations. Initially, the leukemic stem cells are
produced in the bone marrow, and then these cells get transformed to
progenitor cells, which finally differentiate to mature leukemic cells, which
are responsible for leukemic cancer. As observed earlier, the CML typically
progresses in three phases, namely, chronic phase, accelerated phase and blast
phase. The chronic phase is the initial stage of CML, where the myeloblasts
capture up to $10\%$ of blood cells or bone marrow, and these increase to
$10\%-19\%$ in the accelerated phase. The blast phase (acute phase or blast
crisis) is the severe stage where the blood cells or bone marrow contain more
than $20\%$ myeloblasts. Using the theory of multi-type branching processes,
we analyze the three stages of CML progression based on three different types
of leukemic cells, namely, stem cells, progenitor cells and mature leukemic
cells. By calculating the percentage of total leukemic cells in blood, this
study also helps in estimating the time when CML progresses from one phase to
another (chronic, accelerated or blast).

We categorize the cells produced in the three stages of CML to be denoted as
Type-$i$ cells, for $i=0,1,2$. Let us assume that a cell of these three types,
generate new cells, by following the Poisson distribution, with parameter
$\lambda$, which are then called the offsprings of the parent cell. The
offspring mean denoted by $\lambda$ is the same for all types of cells.
Further, assume that each Type-$i$ cell ($i=0,1$) generates both Type-$i$
cells (through cell division) and Type-($i+1$) cells (through cell mutation),
but a Type-$2$ cell generates only Type-$2$ cells (through cell division). Let
$p$ be the probability of the birth of the new cell through mutation, and
hence $(1-p)$ is the probability of cell division for Type-$i$ cells, $i=0,1$.
Further, the assumption is that each of the three types of cells either
generates new cells (through cell division or cell mutation) with probability
$a$ after time $\Delta t$, or dies within the interval $\Delta t$ with
probability $(1-a)$. Thus $a$ is the survival probability for each cell. When
a cell generates progenies, it bursts and hence after generating some new
cells, it dies. We take $\Delta t$ to be a unit time in which this happens.
Let $Z_{i}(t)$ be the number of Type-$i$ cells at time $t$. Suppose that the
development of CML can be traced back to the initial source of $s$ number of
Type-$0$ cells. We denote $X_{i}^{(k)}(t)$ and $Y_{i}^{(k)}(t)$ to be the
number of cells produced by the cell division and cell mutation, respectively,
from the $k$-th individual of Type-$i$ cells, at time $t$. Accordingly, we
obtain,

$\displaystyle Z_{0}(0)$ $\displaystyle=$ $\displaystyle s,$ $\displaystyle
Z_{0}(t+1)$ $\displaystyle=$
$\displaystyle\sum\limits_{k=1}^{Z_{0}(t)}X_{0}^{(k)}(t),$ $\displaystyle
Z_{1}(t+1)$ $\displaystyle=$
$\displaystyle\sum\limits_{k=1}^{Z_{0}(t)}Y_{0}^{(k)}(t)+\sum\limits_{k=1}^{Z_{1}(t)}X_{1}^{(k)}(t),$
$\displaystyle Z_{2}(t+1)$ $\displaystyle=$
$\displaystyle\sum\limits_{k=1}^{Z_{1}(t)}Y_{1}^{(k)}(t)+\sum\limits_{k=1}^{Z_{2}(t)}X_{2}^{(k)}(t).$

We assume that each individual of Type-$0$ cells gives birth to a random
number of Type-$1$ cells denoted by $X$. Hence $X$ is a random variable which
follows the Poisson distribution, with parameter $\lambda$, and this is called
offspring distribution, with $\lambda$ being the offspring mean. Suppose
$X_{i}^{(k)}$ and $Y_{i}^{(k)}$ are independent identically distributed
(i.i.d.) random variable with the same probability distribution as $X$. If
$k$-th individual of Type-$i$ cells has $U_{i}^{(k)}$ offspring at time $t$,
where $U_{i}^{(k)}$’s are i.i.d. random variable, then after $(t+1)$ time,
Type-$i$ cells have total of $\displaystyle{\sum\limits_{k}U_{i}^{(k)}}$
offspring, which is also a Poisson random variable with parameter
($\displaystyle{\lambda\times}$ the number of Type-$i$ cells). The schematic
representation of the model is presented in Figure 1.

$0$$0$$1$$0$$1$$1$$2$$0$$1$$1$$2$$1$$2$$2$ Figure 1: Schematic Representation
of the Model.

##  3 Model Analysis

In this section, we carry out the model analysis in detail. We begin with the
estimation of the total number of offspring (total progeny) after development
of leukemic cancer cells and calculate the time duration of the process. In
the discussion, we consider $\lambda a>1$, as the number of CML cells
initially increases over time.

######  Theorem 1.

The expected number of total leukemic cells at time $t$ is given by,

$E\left[Z(t)|Z(0)=s\right]=s\lambda^{t}a^{t},\leavevmode\nobreak\ t\geq 0.$

###### Proof.

We begin with $Z_{0}(0)=s$ cells. Now, recall that the probability of survival
and the offspring mean for each cell, in a unit interval of time, are given by
$a$ and $\lambda$, respectively. Hence, each of the $s$ cells produces
$\lambda a$ offsprings, on an average, after a unit interval of time.
Therefore, at time $t=1$, the expected number of CML cells (both Type-$0$ and
Type-$1$ cells) becomes $s\lambda a$. Using a similar argument, at time $t=2$,
the expected number of CML cells (Type-$0$, Type-$1$ and Type-$2$ cells) is
given by $s\lambda^{2}a^{2}$. Proceeding in the same manner, we can conclude
that the expected number of total leukemic cells (all three types of cells) at
time $t$ will be $s\lambda^{t}a^{t}$. ∎

######  Theorem 2.

The expected number of Type-$0$ cells at time $t$ is given by,

$E[Z_{0}(t)|Z(0)=s]=s\lambda^{t}a^{t}(1-p)^{t},\leavevmode\nobreak\ t\geq 0.$

###### Proof.

We begin with $Z_{0}(0)=s$ cells. Further, $a$ is the probability of survival
and $\lambda$ is the offspring mean for a cell, in a unit interval of time,
with $(1-p)$ being the probability of cell division. Hence, each of Type-$0$
cells produces $\lambda a(1-p)$ off-springs on an average, after a unit
interval of time. Hence, at time $t=1$, the expected number of Type-$0$ cells
comes out to be $s\lambda a(1-p)$. Similarly, at time $t=2$, the expected
number of Type-$0$ cells is given by $s\lambda^{2}a^{2}(1-p)^{2}$. Proceeding
in the same way, we can conclude that the expected number of Type-$0$ cells at
time $t$ will be $s\lambda^{t}a^{t}(1-p)^{t}$. ∎

######  Theorem 3.

The expected number of Type-$1$ cells at time $t$ is given by,

$E[Z_{1}(t)|Z(0)=s]=ts\lambda^{t}a^{t}p(1-p)^{t-1},\leavevmode\nobreak\ t\geq
1.$

###### Proof.

We use the principle of mathematical induction in order to prove the theorem.
Initially, at time $t=0$, there are no Type-$1$ cell. At time $t=1$, the
expected number of Type-$1$ cells is given by $s\lambda ap$, which comes only
from Type-$0$ cells, through the process of cell mutation. Then Type-$1$ cells
will be produced through both the processes, namely, the cell mutation from
Type-$0$ cells, and the cell division from Type-$1$ cells. Therefore, the
number of Type-$1$ cells at time $t=2$ depends on the number of Type-$0$
cells, as well as the number of Type-$1$ cells at $t=1$. At time $t=2$, the
number of Type-$1$ cells produced from Type-$0$ cells is
$s\lambda^{2}a^{2}p(1-p)$, which is the same as the number of Type-$1$ cells
produced from Type-$1$ cells. Therefore, the expected number of Type-$1$ cells
at time $t=2$ is equal to $2s\lambda^{2}a^{2}p(1-p)$.

Let us now assume that the theorem holds for $t=m$. Hence,

$E[Z_{1}(m)|Z(0)=s]=ms\lambda^{m}a^{m}p(1-p)^{m-1}.$

Now, the value of $\displaystyle{E[Z_{1}(m+1)|Z(0)=s]}$ is dependent on
$\displaystyle{E[Z_{0}(m)|Z(0)=s]}$ and $\displaystyle{E[Z_{1}(m)|Z(0)=s]}$.
Therefore, we have,

$E[Z_{0}(m)|Z(0)=s]=s\lambda^{m}a^{m}(1-p)^{m}.$

The number of Type-$1$ cells produced from
$\displaystyle{s\lambda^{m}a^{m}(1-p)^{m}}$ Type-$0$ cells, through the cell
mutation is $\displaystyle{s\lambda^{m+1}a^{m+1}p(1-p)^{m}}$. Further, the
number of Type-$1$ cells produced from
$\displaystyle{ms\lambda^{m}a^{m}p(1-p)^{m-1}}$ Type-$1$ cells through the
cell division is
$\displaystyle{ms\lambda^{m+1}a^{m+1}p(1-p)^{m}}$. Therefore, we finally
obtain,

$E[Z_{1}(m+1)|Z(0)=s]=(m+1)s\lambda^{m+1}a^{m+1}p(1-p)^{m}.$

Hence, by the principle of mathematical induction, the Theorem holds. ∎

######  Theorem 4.

The expected number of Type-$2$ cells at time $t$ is given by,

$E[Z_{2}(t)|Z(0)=s]=s\lambda^{t}a^{t}p^{2}\left[1+2(1-p)+3(1-p)^{2}+\dots+(t-1)(1-p)^{t-2}\right],\leavevmode\nobreak\
t\geq 2,$

or,

$E[Z_{2}(t)|Z(0)=s]=s\lambda^{t}a^{t}\left[1-(1-p)^{t}-tp(1-p)^{t-1}\right],\leavevmode\nobreak\
t\geq 2.$

###### Proof.

We use the principle of mathematical induction to prove the theorem. At $t=0$
and $t=1$, there is no Type-$2$ cell. At $t=2$, the expected number of
Type-$2$ cells is given by $\displaystyle{s\lambda^{2}a^{2}p^{2}}$, which
comes only from $\displaystyle{s\lambda a(1-p)}$ Type-$1$ cells through the
process of cell mutation. Then, the Type-$2$ cells will be produced through
both the processes, namely, the cell mutation from Type-$1$ cells, and the
cell division from Type-$2$ cells. Therefore, the number of Type-$2$ cells at
time $t=3$ depends on the number of Type-$1$ cells as well as Type-$2$ cells
at time $t=2$. At time $t=3$, the number of Type-$2$ cells that are produced
from $2s\lambda^{2}a^{2}p(1-p)$ Type-$1$ cells is
$2s\lambda^{3}a^{3}p^{2}(1-p)$ and those which are produced from
$\displaystyle{s\lambda^{2}a^{2}p^{2}}$ Type-$2$ cells is
$\displaystyle{s\lambda^{3}a^{3}p^{2}}$. Therefore, the expected number of
Type-$2$ cells at time $t=3$ is equal to
$\displaystyle{s\lambda^{2}a^{2}p^{2}+s\lambda^{3}a^{3}p^{2},\leavevmode\nobreak\
\textit{i.e.,}\leavevmode\nobreak\ s\lambda^{3}a^{3}p^{2}[1+2(1-p)]}$. Let us
consider that the theorem holds for $t=m$. Hence,

$E[Z_{2}(m)|Z(0)=s]=s\lambda^{m}a^{m}p^{2}\left[1+2(1-p)+3(1-p)^{2}+\dots+(t-1)(1-p)^{m-2}\right].$

Now, the value of $\displaystyle{E[Z_{2}(m+1)|Z(0)=s]}$ is dependent on
$\displaystyle{E[Z_{1}(m)|Z(0)=s]}$ and $\displaystyle{E[Z_{2}(m)|Z(0)=s]}$.
We have $\displaystyle{E[Z_{1}(m)|Z(0)=s]=ms\lambda^{m}a^{m}p(1-p)^{m-1}}$.
The number of Type-$2$ cells produced from
$\displaystyle{ms\lambda^{m}a^{m}p(1-p)^{m-1}}$ Type-$1$ cells through the
cell mutation is $\displaystyle{ms\lambda^{m+1}a^{m+1}p^{2}(1-p)^{m-1}}$.
Further, the number of Type-$2$ cells produced from
$\displaystyle{s\lambda^{m}a^{m}p^{2}\left[1+2(1-p)+3(1-p)^{2}+\dots+(t-1)(1-p)^{m-2}\right]}$
Type-$2$ cells through the cell division is
$\displaystyle{s\lambda^{m+1}a^{m+1}p^{2}\left[1+2(1-p)+3(1-p)^{2}+\dots+(t-1)(1-p)^{m-2}\right]}$.
Therefore,

$\displaystyle E[Z_{1}(m+1)|Z(0)=s]$ $\displaystyle=$ $\displaystyle
ms\lambda^{m+1}a^{m+1}p^{2}(1-p)^{m-1}$ $\displaystyle+$ $\displaystyle
s\lambda^{m+1}a^{m+1}p^{2}\left[1+2(1-p)+3(1-p)^{2}+\dots+(t-1)(1-p)^{m-2}\right]$
$\displaystyle=$ $\displaystyle
s\lambda^{m+1}a^{m+1}p^{2}\left[1+2(1-p)+3(1-p)^{2}+\dots+(m-1)(1-p)^{m-2}+m(1-p)^{m-1}\right]$

Hence, by the principle of mathematical induction, the Theorem holds.
Furthermore, we obtained in Theorem 1 that
$E[Z(t)|Z(0)=s]=E[Z_{0}(t)|Z(0)=s]+E[Z_{1}(t)|Z(0)=s]+E[Z_{2}(t)|Z(0)=s]=s\lambda^{t}a^{t}.$
Hence,

$E[Z_{2}(t)|Z(0)=s]=s\lambda^{t}a^{t}-E[Z_{0}(t)|Z(0)=s]-E[Z_{1}(t)|Z(0)=s]=s\lambda^{t}a^{t}\left[1-(1-p)^{t}-tp(1-p)^{t-1}\right],\leavevmode\nobreak\
t\geq 2.$

∎

The results presented in Theorems 1 to 4 are presented in a summarized form in
Table 1.

$i$ | Expected number of Type-$i$ cells
---|---
$i=0$ |  $E[Z_{0}(t)|Z(0)=s]=s\lambda^{t}a^{t}(1-p)^{t},\leavevmode\nobreak\ t\geq 0$.
$i=1$ |  $E[Z_{1}(t)|Z(0)=s]=ts\lambda^{t}a^{t}p(1-p)^{t-1},\leavevmode\nobreak\ t\geq 1$.
$i=2$ |  $E[Z_{2}(t)|Z(0)=s]=s\lambda^{t}a^{t}\left[1-(1-p)^{t}-tp(1-p)^{t-1}\right],\leavevmode\nobreak\ t\geq 2$.
Table 1:  Summary of Theorems 1 to 4 .

######  Theorem 5.

The time to reach the leukemia of $M$ cells on average is given by,

$T=\frac{\ln\left(M/s\right)}{\ln(\lambda a)},\leavevmode\nobreak\ \lambda
a>1.$

###### Proof.

If $T$ be the time required to reach $M$ number of leukemic cells on average,
then we can write, from Theorem 1,

$s\lambda^{T}a^{T}=M,$

which gives,

$T=\frac{\ln\left(M/s\right)}{\ln(\lambda a)},\leavevmode\nobreak\ \lambda
a>1.$

∎

######  Theorem 6.

If $T_{12}$ be the time required to reach injurious leukemic cells of average
size $M_{12}$, the following relation holds,

$s\lambda^{T_{12}}a^{T_{12}}\left[1-(1-p)^{T_{12}}\right]=M_{12}.$

###### Proof.

Let $T_{12}$ be the time required to reach $M_{12}$ number of injurious
leukemic cells. Then, from Theorem 3 and Theorem 4, we arrive at the
following,

$\displaystyle M_{12}$ $\displaystyle=$ $\displaystyle
T_{12}s\lambda^{T_{12}}a^{T_{12}}p(1-p)^{T_{12}-1}+s\lambda^{T_{12}}a^{T_{12}}\left[1-(1-p)^{T_{12}}-T_{12}p(1-p)^{T_{12}-1}\right]$
$\displaystyle=$ $\displaystyle
s\lambda^{T_{12}}a^{T_{12}}\left[1-(1-p)^{T_{12}}\right].$

∎

Finally, we determine the probability of CML extinction, under the model
considerations. The event $\\{Z(t)=0,\leavevmode\nobreak\ \text{for
some}\leavevmode\nobreak\ t>0\leavevmode\nobreak\ |\leavevmode\nobreak\
Z(0)=s\\}$ is called the disease extinction. Then we have the following
Theorem.

######  Theorem 7.

If $P(t)$ be the probability of extinction of CML at time $t$, then,

$P(t)\geq(1-a)^{s\lambda^{t-1}a^{t-1}},\leavevmode\nobreak\ t\geq 1.$

###### Proof.

According to Theorem 1, the expected number of total leukemic cells at time
$t$ is $s\lambda^{t}a^{t}$. The disease CML becomes extinct at time $t$ if all
leukemic cells produced till time $t-1$ die before $t-1$ or during $t-1$ to
$t$. Hence, the probability of extinction of CML at time $t$ is given by,

$P(t)\geq(1-a)^{s\lambda^{t-1}a^{t-1}},\leavevmode\nobreak\ t\geq 1.$

∎

We now discuss about the consequences of the primary therapy on CML patients.
The primary therapy for CML prevents the expansion of Type-0 cells, but it
does not have a direct effect on Type-1 and Type-2 cells, since these cells
are resistant to primary therapy. The effect on Type-0 cells, due to primary
therapy results in changes in the development of Type-1 and Type-2 cells.
Accordingly, we let the effect of the primary therapy be manifested through
change in the probability of cell survival and cell mutation for Type-0 cells
to $\bar{a}$ and $\bar{p}$, respectively, whereas those for Type-1 and Type-2
cells remains the same as prior to therapy. Then the results can be written in
the following Theorems.

######  Theorem 8.

If the primary CML therapy is initiated, the expected number of total leukemic
cells at time $t$ is given by,

$E[\bar{Z}(t)|Z(0)=s]=s\lambda^{t}\bar{a}\left[\bar{p}a^{t-1}+\bar{p}a^{t-2}\\{\bar{a}(1-\bar{p})\\}+\bar{p}a^{t-3}\\{\bar{a}(1-\bar{p})\\}^{2}+\cdots+\bar{p}a\\{\bar{a}(1-\bar{p})\\}^{t-2}+\\{\bar{a}(1-\bar{p})\\}^{t-1}\right].$

###### Proof.

In order to prove this Theorem, we consider the following three facts:

  * •

Type-0 cells produce the same cells at a rate $\lambda$, with probability
$1-\bar{p}$ (cell division), where the survival probability of Type-0 cells is
$\bar{a}$.

  * •

Type-0 cells produce Type-1 cells at a rate $\lambda$, with probability
$\bar{p}$ (cell mutation), where the survival probability of Type-0 cells is
$\bar{a}$.

  * •

Both Type-1 and Type-2 cells produce new leukemic cells of any type at a rate
$\lambda$, with probability $1$ (cell division and mutation together), where
the survival probability of Type-1 or Type-2 cells is $a$.

Therefore, the number of total leukemic cells at $t=n$ is dependent on the
following three evaluations:

  * •

Number of Type-0 cells produced from Type-0 cells, present at $t=n-1$, through
cell division.

  * •

Number of Type-1 cells produced from Type-0 cells, present at $t=n-1$, through
cell mutation.

  * •

Total number of Type-1 and Type-2 cells produced from both of Type-1 and
Type-2 cells, present at $t=n-1$, through cell division and mutation together.

  1. (A)

Calculation of total leukemic cells at $t=1$,

     * •

The number of Type-0 cells produced from Type-0 cells, present at $t=0$,
through cell division is $s\lambda\bar{a}(1-\bar{p})$.

     * •

The number of Type-1 cells produced from Type-0 cells, present at $t=0$,
through cell mutation is $s\lambda\bar{a}\bar{p}$.

     * •

There is no such type of leukemic cells which is produced from Type-1 or
Type-2 cells, present at $t=0$.

Therefore, at $t=1$, the number of total leukemic cells becomes
$s\lambda\bar{a}(1-\bar{p})+s\lambda\bar{a}\bar{p}=s\lambda\bar{a}$.

  2. (B)

Calculation of total leukemic cells at $t=2$,

     * •

The number of Type-0 cells produced from Type-0 cells, present at $t=1$,
through cell division is $s\lambda^{2}\bar{a}^{2}(1-\bar{p})^{2}$.

     * •

The number of Type-1 cells produced from Type-0 cells, present at $t=1$,
through cell mutation is $s\lambda^{2}\bar{a}^{2}\bar{p}(1-\bar{p})$.

     * •

Total number of Type-1 and Type-2 cells produced from both of Type-1 and
Type-2 cells, present at $t=1$, through cell division and mutation together is
$s\lambda^{2}a\bar{a}\bar{p}$.

Therefore, at $t=2$, the number of total leukemic cells is,

$s\lambda^{2}\bar{a}^{2}(1-\bar{p})^{2}+s\lambda^{2}\bar{a}^{2}\bar{p}(1-\bar{p})+s\lambda^{2}a\bar{a}\bar{p}=s\lambda^{2}\bar{a}\left[\bar{p}a+\\{\bar{a}(1-\bar{p})\\}\right].$

  3. (C)

Calculation of total leukemic cells at $t=3$,

     * •

The number of Type-0 cells produced from Type-0 cells, present at $t=2$,
through cell division is $s\lambda^{3}\bar{a}^{3}(1-\bar{p})^{3}$.

     * •

The number of Type-1 cells produced from Type-0 cells, present at $t=2$,
through cell mutation is $s\lambda^{3}\bar{a}^{3}\bar{p}(1-\bar{p})^{2}$.

     * •

Total number of Type-1 and Type-2 cells produced from both of Type-1 and
Type-2 cells, present at $t=2$, through cell division and mutation together is
$s\lambda^{3}a\bar{a}^{2}\bar{p}(1-\bar{p})+s\lambda^{3}a^{2}\bar{a}\bar{p}$.

Therefore, at $t=3$, the number of total leukemic cells is,

$\displaystyle
s\lambda^{3}\bar{a}^{3}(1-\bar{p})^{3}+s\lambda^{3}\bar{a}^{3}\bar{p}(1-\bar{p})^{2}+s\lambda^{3}a\bar{a}^{2}\bar{p}(1-\bar{p})+s\lambda^{3}a^{2}\bar{a}\bar{p}$
$\displaystyle=$ $\displaystyle
s\lambda^{3}\bar{a}\Big{[}\bar{p}a^{2}+\bar{p}a\\{\bar{a}(1-\bar{p})\\}+\\{\bar{a}^{2}(1-\bar{p})^{2}\\}\Big{]}.$

Proceeding in the same manner, one can obtain that the expected number of
total leukemic cells at time $t$ is equal to,

$s\lambda^{t}\bar{a}\left[\bar{p}a^{t-1}+\bar{p}a^{t-2}\\{\bar{a}(1-\bar{p})\\}+\bar{p}a^{t-3}\\{\bar{a}(1-\bar{p})\\}^{2}+\cdots+\bar{p}a\\{\bar{a}(1-\bar{p})\\}^{t-2}+\\{\bar{a}(1-\bar{p})\\}^{t-1}\right].$

∎

######  Theorem 9.

If the primary CML therapy is initiated, the expected number of Type-$0$ cells
at time $t$ is given by,

$E[\bar{Z}_{0}(t)|Z(0)=s]=s\lambda^{t}{\bar{a}}^{t}(1-{\bar{p}})^{t}.$

###### Proof.

The proof is similar to Theorem 2. ∎

######  Theorem 10.

If the primary CML therapy is initiated, the expected number of Type-$1$ cells
at time $t$ is given by

$\displaystyle E[\bar{Z}_{1}(t)|Z(0)=s]$ $\displaystyle=$ $\displaystyle
s\lambda^{t}\bar{a}\bar{p}\left[\\{a(1-p)\\}^{t-1}+\\{a(1-p)\\}^{t-2}\\{\bar{a}(1-\bar{p})\\}+\\{a(1-p)\\}^{t-3}\\{\bar{a}(1-\bar{p})\\}^{2}+\cdots\right.$
$\displaystyle\leavevmode\nobreak\ +\leavevmode\nobreak\
\left.\\{a(1-p)\\}\\{\bar{a}(1-\bar{p})\\}^{t-2}+\\{\bar{a}(1-\bar{p})\\}^{t-1}\right].$

###### Proof.

In order to prove this theorem, we consider the following two facts:

  * •

Type-0 cells produce Type-1 cells at a rate $\lambda$ with probability
$\bar{p}$ (cell mutation), where the survival probability of Type-0 cells is
$\bar{a}$.

  * •

Type-1 cells produce the same cells at a rate $\lambda$ with probability $1-p$
(cell division), where the survival probability of Type-0 cells is $a$.

Therefore, the number of Type-1 cells at $t=n$ is dependent on the following
two evaluations:

  * •

Number of Type-1 cells produced from Type-0 cells, present at $t=n-1$, through
cell mutation.

  * •

Number of Type-1 cells produced from Type-1 cells, present at $t=n-1$, through
cell division.

  1. (A)

Calculation of Type-1 cells at $t=1$,

     * •

The number of Type-1 cells produced from Type-0 cells, present at $t=0$,
through cell mutation is $s\lambda\bar{a}\bar{p}$.

     * •

There is no such Type-1 cells which is produced from Type-1 cells present at
$t=0$.

Therefore, at $t=1$, the number of Type-1 cells becomes
$s\lambda\bar{a}\bar{p}$.

  2. (B)

Calculation of Type-1 cells at $t=2$,

     * •

The number of Type-1 cells, produced from Type-0 cells present at $t=1$,
through cell mutation is $s\lambda^{2}\bar{a}^{2}\bar{p}(1-\bar{p})$.

     * •

The number of Type-1 cells produced from Type-1 cells, present at $t=1$,
through cell division is $s\lambda^{2}a\bar{a}\bar{p}(1-p)$.

Therefore, at $t=2$, the number of Type-1 cells becomes
$s\lambda^{2}\bar{a}^{2}\bar{p}(1-\bar{p})+s\lambda^{2}a\bar{a}\bar{p}(1-p)=s\lambda^{2}\bar{a}\bar{p}\Big{[}a(1-p)+\bar{a}(1-\bar{p})\Big{]}$.

  3. (C)

Calculation of Type-1 cells at $t=3$,

     * •

The number of Type-1 cells produced from Type-0 cells, present at $t=2$,
through cell mutation is $s\lambda^{3}\bar{a}^{3}\bar{p}(1-\bar{p})^{2}$.

     * •

The number of Type-1 cells produced from Type-1 cells, present at $t=2$,
through cell division is
$s\lambda^{3}a\bar{a}^{2}\bar{p}(1-p)(1-\bar{p})+s\lambda^{3}a^{2}\bar{a}\bar{p}(1-p)^{2}$.

Therefore, at $t=3$, the number of Type-1 cells is

$\displaystyle
s\lambda^{3}\bar{a}^{3}\bar{p}(1-\bar{p})^{2}+s\lambda^{3}a\bar{a}^{2}\bar{p}(1-p)(1-\bar{p})+s\lambda^{3}a^{2}\bar{a}\bar{p}(1-p)^{2}$
$\displaystyle=$ $\displaystyle
s\lambda^{3}\bar{a}\bar{p}\Big{[}\\{a(1-p)\\}^{2}+\\{a(1-p)\\}\\{\bar{a}(1-\bar{p})\\}+\\{\bar{a}(1-\bar{p})\\}^{2}\Big{]}.$

Proceeding in the same manner, one can obtain that the expected number of
Type-1 cells at time $t$ equals,

$s\lambda^{t}\bar{a}\bar{p}\Big{[}\\{a(1-p)\\}^{t-1}+\\{a(1-p)\\}^{t-2}\\{\bar{a}(1-\bar{p})\\}+\\{a(1-p)\\}^{t-3}\\{\bar{a}(1-\bar{p})\\}^{2}+\cdots+\\{a(1-p)\\}\\{\bar{a}(1-\bar{p})\\}^{t-2}+\\{\bar{a}(1-\bar{p})\\}^{t-1}\Big{]}.$

∎

######  Theorem 11.

If the primary CML therapy is initiated, the expected number of Type-$2$ cells
at time $t$ can be calculated using the following relation,

$E[\bar{Z}_{2}(t)|Z(0)=s]=E[\bar{Z}(t)|Z(0)=s]-E[\bar{Z}_{0}(t)|Z(0)=s]-E[\bar{Z}_{1}(t)|Z(0)=s].$

###### Proof.

The proof follows from Theorem 8, Theorem 9 and Theorem 10. ∎

##  4 Results

In this penultimate Section, we present various illustrative numerical
results, pertaining to the derivations made in Section 3. The levels of the
three types of CML cells, as well as the total number of cells, as they
progress over time are numerically illustrated. We also demonstrate how the
number of CML cells can be impacted by the values of the parameters, such as
the initial levels of the leukemic stem cells, the probability of survival of
the CML cells, cell division or mutation rates, and the offspring mean. The
simulation is carried out in order to observe the probability of extinction of
the disease, with respect to the variation in the values of the parameters.
Finally, we determine the optimal time, at which the progression of the
disease changes from one state to another. For the purpose of some
illustrative results, we choose the values of the parameters to be, $s=1$,
$p=0.3$, $a=0.5$, $\lambda=3$, with the time window under consideration being
$T=50$.

We begin with the results for the progression of the cells over time.
Accordingly, Figure 2 shows the changes in the levels of the three types of
CML cells over time. We observe from the figure, that from the initial time of
$t=0$, until the time of $t=40$, the leukemic stem cells and the progenitor
cells, gradually increase with time, but there is no noticeable increment in
the levels of the mature leukemic cells. However, beyond this point, there is
a rapid and accelerated growth of the mature leukemic cells, thereby resulting
in the host being predominantly invaded by these cells.

Figure 2: Progression of the leukemic stem cells ($Z_{0}$), the progenitor
cells ($Z_{1}$) and the mature leukemic cells ($Z_{2}$)

In Figure 3, we plot the percentage of injurious cells, defined as
$\displaystyle{\frac{Z_{1}+Z_{2}}{Z_{0}+Z_{1}+Z_{2}}}\times 100$, as a
function of time. It captures the progression of percentage of the progenitor
cells and the mature leukemic cells, as a percentage of the total number of
the leukemic cells, over time. As can be seen from Figure 3, the injurious
cells, for all practical purpose have taken over the host after about time
$t=12$, since the detection.

Figure 3: Percentage of injurious cells as a function of $t$

Now, we present the three dimensional results for the progression of the three
types of cells, against the ranges of $a\in[0,0.5]$ and $\lambda\in[0,3]$ in
Figure 4. It can be observed that until one reaches fairly close to the
maximum values in the ranges of $a$ and $\lambda$, the progression of the
disease is reasonably under control, beyond which the cell levels dramatically
reach fatal levels in the host.

((a)) $Z_{0}$ as a function of $a$ and $\lambda$

((b)) $Z_{1}$ as a function of $a$ and $\lambda$

((c)) $Z_{2}$ as a function of $a$ and $\lambda$

Figure 4: Progression of the leukemic stem cells ($Z_{0}$), the progenitor
cells ($Z_{1}$) and the mature leukemic cells ($Z_{2}$) as a function of $a$
and $\lambda$

We next look at the probability of extinction for various parameter values in
Figure 5. For the first figure, we choose the ranges of values of $\lambda$ to
be $[0,3]$. The extinction is almost certain upto nearly the value of
$\lambda=2$. However, beyond this point, there is a very small window of
$\lambda$ (among the values of $\lambda a>1$) for which there is a positive
probability of extinction, and then the probability is virtually nonexistent
(demonstrating extremely poor prognosis).

Figure 5: Probability of extinction for various values of $\lambda$ and $t$

Finally, we look at the illustration for the results in Theorems 5 and 6. For
the implementation of the result in Theorem 5, we consider the values of
$0\leq M\leq 10^{8}$, and plot the corresponding expected times $T$ to reach a
level of $M$ leukemic cells can be seen in Figure 6. This result helps in
determining the expected time at which the CML disease changes from being the
less injurious chronic phase, to the more serious and fatal condition based on
the percentage of leukemic cells in the patient’s blood. For the result in
Theorem 6, we consider the ranges $0\leq T_{12}\leq 45$ and $10^{2}\leq
M_{12}\leq 10^{8}$, and present the implicit plot for the same in Figure 7.

Figure 6: Expected times $T$ to reach a level of $M$ leukemic cells Figure 7:
Implicit plot of expected time $T_{12}$ to reach $M_{12}$ injurious leukemic
cells

##  5 Conclusion

In this study, the evolution of CML is discussed through various stages of
disease progression with three types of leukemic cells. We have analytically
calculated the expected number of total leukemic cells, as well as three types
of leukemic cells, before and after initiation of the therapy. The expected
times to reach a certain number of total leukemic cells and injurious leukemic
cells are estimated. The probability of extinction of the CML is also
analytically determined. The progression of different leukemic cells are
graphically plotted which exhibited a sudden and large increment in Type-$2$
cells, as compared to other two types of cells. The numerical findings also
suggests that after a while, the leukemic cells are dominated by the injurious
cells. The changes in the number of leukemic cells, contingent on the values
of survival probability and offspring mean are shown graphically. The chance
of disease extinction is numerically shown for different values of offspring
mean, which depicts that the chance of CML extinction is almost impossible
without therapy, when the offspring mean exceeds $2$. Thus, the findings
obtained can help in determining the percentage of leukemic cells in a
patient’s blood at a specific time and to diagnose the stage which the patient
is going through. Also, it can help in estimating the time needed to go to the
next stage of CML, and to predict the condition of the patient in advance, so
that the medical practitioner can advice on the appropriate therapeutic
intervention.

## Acknowledgment

SPC was supported by Grant. No. MTR/2019/000225 from the Science and
Engineering Research Board, India. SR’s work was supported partially by Grant
No. MTR/2020/000186 from the Science and Engineering Research Board, India.

## References

  * [1] Clapp, G.D., 2016. Applying mathematical models to study the role of the immune system in chronic myelogenous leukemia (Doctoral dissertation).
  * [2] Pujo-Menjouet, L. and Mackey, M.C., 2004. Contribution to the study of periodic chronic myelogenous leukemia. Comptes Rendus Biologies, 327(3), pp.235-244.
  * [3] Komarova, N.L., Katouli, A.A. and Wodarz, D., 2009. Combination of two but not three current targeted drugs can improve therapy of chronic myeloid leukemia. PLoS One, 4(2), p.e4423.
  * [4] Koch, D., Eisinger, R.S. and Gebharter, A., 2017. A causal Bayesian network model of disease progression mechanisms in chronic myeloid leukemia. Journal of Theoretical Biology, 433, pp.94-105.
  * [5] Michor, F., Hughes, T.P., Iwasa, Y., Branford, S., Shah, N.P., Sawyers, C.L. and Nowak, M.A., 2005. Dynamics of chronic myeloid leukaemia. Nature, 435(7046), pp.1267-1270.
  * [6] Komarova, N.L. and Wodarz, D., 2007. Effect of cellular quiescence on the success of targeted CML therapy. PLoS One, 2(10), p.e990.
  * [7] Roeder, I., Herberg, M. and Horn, M., 2009. An “age” structured model of hematopoietic stem cell organization with application to chronic myeloid leukemia. Bulletin of Mathematical Biology, 71(3), pp.602-626.
  * [8] Radulescu, I.R., Candea, D. and Halanay, A., 2014. A study on stability and medical implications for a complex delay model for CML with cell competition and treatment. Journal of Theoretical Biology, 363, pp.30-40.
  * [9] Clapp, G.D., Lepoutre, T., El Cheikh, R., Bernard, S., Ruby, J., Labussière-Wallet, H., Nicolini, F.E. and Levy, D., 2015. Implication of the autologous immune system in BCR–ABL transcript variations in chronic myelogenous leukemia patients treated with Imatinib. Cancer Research, 75(19), pp.4053-4062.
  * [10] Besse, A., Clapp, G.D., Bernard, S., Nicolini, F.E., Levy, D. and Lepoutre, T., 2018. Stability analysis of a model of interaction between the immune system and cancer cells in chronic myelogenous leukemia. Bulletin of Mathematical Biology, 80(5), pp.1084-1110.
  * [11] Ainseba, B.E. and Benosman, C., 2010. Optimal control for resistance and suboptimal response in CML. Mathematical Biosciences, 227(2), pp.81-93.
  * [12] Radulescu, I.R., Candea, D. and Halanay, A., 2016. Optimal control analysis of a leukemia model under imatinib treatment. Mathematics and Computers in Simulation, 121, pp.1-11.
  * [13] Bunimovich-Mendrazitsky, S. and Shklyar, B., 2017. Optimization of Combined Leukemia Therapy by Finite-Dimensional Optimal Control Modeling. Journal of Optimization Theory and Applications, 175(1), pp.218-235.
  * [14] Leder, K., Foo, J., Skaggs, B., Gorre, M., Sawyers, C.L. and Michor, F., 2011. Fitness conferred by BCR-ABL kinase domain mutations determines the risk of pre-existing resistance in chronic myeloid leukemia. PloS One, 6(11), p.e27682.
  * [15] Pan, S., Raha, S. and Chakrabarty, S.P., 2020. A quantitative study on the role of TKI combined with Wnt/$\beta$-catenin signaling and IFN-$\alpha$ in the treatment of CML through deterministic and stochastic approaches. Chaos, Solitons & Fractals, 133, p.109627.
  * [16] Farrington, C.P. and Grant, A.D., 1999. The distribution of time to extinction in subcritical branching processes: applications to outbreaks of infectious disease. Journal of Applied Probability, pp.771-779.
  * [17] Durrett, R. and Moseley, S., 2010. Evolution of resistance and progression to disease during clonal expansion of cancer. Theoretical Population Biology, 77(1), pp.42-48.
  * [18] Durrett, R., Foo, J., Leder, K., Mayberry, J. and Michor, F., 2010. Evolutionary dynamics of tumor progression with random fitness values. Theoretical Population Biology, 78(1), pp.54-66.
  * [19] Danesh, K., Durrett, R., Havrilesky, L.J. and Myers, E., 2012. A branching process model of ovarian cancer. Journal of Theoretical Biology, 314, pp.10-15.
  * [20] Lenaerts, T., Pacheco, J.M., Traulsen, A. and Dingli, D., 2010. Tyrosine kinase inhibitor therapy can cure chronic myeloid leukemia without hitting leukemic stem cells. Haematologica, 95(6), p.900.
  * [21] Berezansky, L., Bunimovich-Mendrazitsky, S. and Domoshnitsky, A., 2012. A mathematical model with time-varying delays in the combined treatment of chronic myeloid leukemia. Advances in Difference Equations, 2012(1), pp.1-13.


11institutetext:  Department of Physics, University of Fribourg, 1700
Fribourg, Switzerland
Department of Physics and Astronomy, University of Tennessee, Knoxville,
Tennessee 37996-1200, USA
Department of Physics, University of Erlangen-Nuremberg, 91058 Erlangen,
Germany

# Nonequilibrium-DMFT based RIXS investigation of the two-orbital Hubbard
model

Philipp Werner  Steven Johnston 11 and Martin Eckstein 2233112233

###### Abstract

Resonant inelastic X-ray scattering (RIXS) detects various types of high- and
low-energy elementary excitations in correlated solids, and this tool will
play an increasingly important role in investigations of time-dependent
phenomena in photo-excited systems. While theoretical frameworks for the
computation of equilibrium RIXS spectra are well established, the development
of appropriate methods for nonequilibrium simulations are an active research
field. Here, we apply a recently developed nonequilibrium dynamical mean field
theory (DMFT) based approach to compute the RIXS response of photo-excited
two-orbital Mott insulators. The results demonstrate the feasibility of multi-
orbital nonequilibrium RIXS calculations and the sensitivity of the quasi-
elastic fluorescence-like features and $d$-$d$ excitation peaks on the
nonequilibrium population of the Hubbard bands.

###### pacs:

71.10.Fd

##  1 Introduction

Resonant inelastic X-ray scattering (RIXS) [1] is a powerful experimental tool
for probing the charge [2], spin [3, 4], orbital [5, 6] and lattice [7, 8, 9]
degrees of freedom in correlated solids. It is a photon-in, photon-out
technique, where the incoming photon with energy $\omega_{\text{in}}$ and
momentum $\vec{q}_{\text{in}}$ excites an electron from a core level into the
conduction band. What is measured is the photon emitted with energy
$\omega_{\text{out}}$ and momentum $\vec{q}_{\text{out}}$, when the excited
electron or another one fills the core hole. A recent development are time-
resolved RIXS experiments [10, 11, 12], which track the RIXS signal in a
material driven out of its equilibrium state by a laser pulse. Thanks to a
high time resolution, these measurements can reveal the interplay between
spin, orbital, and lattice degrees of freedom during the laser driving and
relaxation dynamics.

The equilibrium RIXS amplitude $I$ can be obtained from cluster
diagonalization and the Kramers-Heisenberg formula [13, 14]
$I_{\vec{q}_{\text{in}},\vec{q}_{\text{out}}}(\omega_{\text{in}},\omega_{\text{out}})=\frac{1}{\pi}\text{Im}\langle\psi|(H_{\text{cl}}-E_{0}-(\omega_{\text{in}}-\omega_{\text{out}})-i0_{+})^{-1}|\psi\rangle$,
with
$|\psi\rangle=\sum_{m}e^{i(\vec{q}_{\text{in}}-\vec{q}_{\text{out}})\cdot\vec{r}_{m}}D^{\dagger}_{m}(H_{\text{cl}}-E_{0}-\omega_{\text{in}}-i\Gamma)^{-1}D_{m}|\phi_{0}\rangle.$
Here, $H_{\text{cl}}$ is the Hamiltonian of the cluster (including core
levels), $E_{0}$ is the energy of the ground state $|\phi_{0}\rangle$, $D_{m}$
is the dipole transition operator for lattice site $m$, which can also depend
on the photon polarization and electron spin, and $\Gamma$ is a parameter
determining the lifetime of the core hole. These calculations capture the
atomic multiplet structure very well, but provide only a crude description of
the itinerant nature of the conduction electrons.

The nonequilibrium generalization of this cluster based RIXS formalism has
recently been formulated by Devereaux, Freericks and co-workers in Ref. [14].
It is based on the measurement of the four-point correlation function

$\displaystyle
S_{mn}(t_{1}^{\phantom{\prime}},t_{2}^{\phantom{\prime}},t_{1}^{\prime},t_{2}^{\prime})=\langle
U(-\infty,t_{1}^{\prime})D_{n}(t_{1}^{\prime})U(t_{1}^{\prime},t_{2}^{\prime})D^{\dagger}_{n}(t_{2}^{\prime})$
$\displaystyle\hskip 14.22636pt\times
U(t_{2}^{\prime},t_{2}^{\phantom{\prime}})D_{m}(t_{2})U(t_{2},t_{1})D^{\dagger}_{m}(t_{1})U(t_{1},-\infty)\rangle,$
(1)

with $U$ the time evolution operator, and has initially been demonstrated for
a noninteracting problem, where this complicated correlation function can be
explicitly calculated. The application of the method to interacting systems is
numerically challenging, because of the difficulties of measuring and storing
$S_{mn}(t_{1}^{\phantom{\prime}},t_{2}^{\phantom{\prime}},t_{1}^{\prime},t_{2}^{\prime})$.
It has been successfully used to study a single-band Hubbard model with
additional core-hole potential on a small cluster [19], but extensions of this
approach to multi-orbital problems will be extremely costly. Moreover, while a
small cluster may be sufficient for capturing the dynamics during the RIXS
process, finding a finite cluster that can accurately represent the non-
equilibrium dynamics over a longer period of time after a pump laser
excitation is challenging.

In a separate effort, a dynamical mean field theory (DMFT) [15] based
equilibrium RIXS formalism has been developed by Hariki, Kunes, and co-workers
[16, 17]. In this scheme, the orbitals near the Fermi level are replaced by an
impurity problem with a hybridization function derived from a DMFT
calculation. The advantage of this approach is that it can treat the itinerant
nature of the conduction electrons via the DMFT construction. It thus allows
one to compute a fluorescence-like RIXS signal with a linear dependence of
$\omega_{\text{loss}}=\omega_{\text{in}}-\omega_{\text{out}}$ on
$\omega_{\text{in}}$. A current limitation of this method is that it does not
provide information on the $\vec{q}$-dependence of the RIXS spectra.

In Ref. [20] we showed how the DMFT based RIXS approach can be implemented on
the real-time axis, and thereby extended to nonequilibrium problems. The basic
idea is to simulate the RIXS probe pulse explicitly using nonequilibrium DMFT
[21], and to measure a time-dependent correlation function from which the
number of emitted photons can be obtained. We will use here an implementation
without rotating wave approximation, and with the noncrossing approximation
(NCA) [22, 23] as impurity solver (for details, see Ref. [20]), to investigate
the equilibrium and nonequilibrium RIXS spectra of the half-filled an quarter-
filled two-orbital Hubbard model in the Mott insulating regime. In particular,
we will show how these spectra are modified by a short but strong photo-doping
pulse.

##  2 Model and method

We consider a two-orbital Hubbard model with a coupling to core levels. The
local (atomic) Hamiltonian reads

$\displaystyle H_{\text{loc}}=$
$\displaystyle\sum_{\alpha}Un_{\alpha\uparrow}n_{\alpha\downarrow}$
$\displaystyle+\sum_{\sigma}[(U-2J)n_{1\sigma}n_{2\bar{\sigma}}+(U-3J)n_{1\sigma}n_{2\sigma}]$
$\displaystyle+U_{cd}n_{c}n_{d}+\tfrac{\Delta}{2}(n_{d}-2n_{c})-\mu(n_{d}+n_{c}),$
(2)

where $\alpha=1,2$ denotes the orbitals, $\sigma=\uparrow,\downarrow$ denotes
spin, $d_{\alpha\sigma}$ and $n_{\alpha\sigma}$ are the corresponding
annihilation and density operators, $c_{\sigma}$ and $n_{c\sigma}$ are the
annihilation and density operators for the core level, $U$ is the intra-
orbital interaction, $J$ is the Hund coupling, $U_{cd}$ is the interaction
with the core level, and $\mu$ is the chemical potential. We use the notations
$n_{d}=\sum_{\alpha,\sigma}n_{\alpha\sigma}$, $n_{c}=\sum_{\sigma}n_{c\sigma}$
and introduce a parameter $\Delta$ to shift the core energy relative to the
$d$ levels. To mimic the probe pulse, we add a dipolar coupling between the
core level and the first $d$-orbital:
$H_{\text{probe}}=E_{\text{probe}}(t)\sum_{\sigma}(d^{\dagger}_{1\sigma}c^{\phantom{\dagger}}_{\sigma}+h.c.)$,
with a field of the form
$E_{\text{probe}}(t)=E_{\text{probe}}f_{\text{probe}}(t-t_{\text{probe}})\sin(\omega_{\text{in}}(t-t_{\text{probe}}))$
whose amplitude $E_{\text{probe}}$ is weak enough that the RIXS signal is
quadratic in $E_{\text{probe}}$, $\omega_{\text{in}}$ the incoming frequency
of the order of the $c$–$d$ energy splitting, and $f(t-t_{\text{probe}})$ an
envelope function centered at time $t_{\text{probe}}$. The $d$-electrons hop
between lattice sites with the kinetic term $H_{\text{kin}}=\sum_{\langle
i,j\rangle,\alpha,\sigma}v(t)(d^{\dagger}_{i\alpha\sigma}d^{\phantom{\dagger}}_{j\alpha\sigma}+h.c.)$.
A time-dependent modulation of the hopping parameter $v$ with frequency
$\omega_{\text{pump}}$ of the order of $U$, amplitude $v_{\text{pump}}$ and
envelope $f_{\text{pump}}$ will be used to “photo-excite” the system:
$v(t)=v_{0}+v_{\text{pump}}f_{\text{pump}}(t-t_{\text{pump}})\sin(\omega_{\text{pump}}(t-t_{\text{pump}}))$.
Furthermore, to mimic the finite core-hole lifetime (parameter $\Gamma$ in the
Kramers-Heisenberg approach), we couple a fermionic bath with a box-shaped
density of states (DOS) to the core level.

We consider an infinite-dimensional Bethe lattice with appropriately
renormalized hopping amplitude, so that the DMFT self-consistency gives a
direct relation between the impurity hybridization function
$\Lambda_{\nu,\mu,\sigma}$ and the local Green’s functions
$G_{\mu,\nu,\sigma}$ (with $\nu,\mu=(1,2,c)$) [15]:
${\bf\Lambda_{\sigma}}(t,t^{\prime})={\bf V}(t){\bf
G_{\sigma}}(t,t^{\prime}){\bf V}(t^{\prime})$. Here, the bold symbols
represent $3\times 3$ matrices. We choose $V_{11}(t)=V_{22}(t)=v(t)$ and
$V_{cc}(t)=0$ (localized core level). The hybridization function
${\bf\Lambda_{\sigma}}$ of the (equilibrium or photo-excited) lattice system
represents the environment of the probed site and is not modified in the
simulation of the RIXS process. The RIXS probe pulse, however, introduces
coherence between the $c$ and $d_{1}$ orbitals, and produces nonzero off-
diagonal elements in ${\bf G_{\sigma}}$.

The RIXS signal can be computed from the correlation functions
$D_{\alpha,\sigma,\sigma^{\prime}}(t,t^{\prime})=-i\langle
T_{\mathcal{C}}P^{\phantom{\dagger}}_{\alpha\sigma}(t)P_{\alpha\sigma^{\prime}}^{\dagger}(t^{\prime})\rangle$
with
$P_{\alpha\sigma}=c^{\dagger}_{\sigma}d^{\phantom{\dagger}}_{\alpha\sigma}$,
see Ref. [20] for details. In the following, we will focus on the signal
emitted by orbital $\alpha=1$, and use $v_{0}=1$ as the unit of energy
($\hbar/v_{0}$ as the unit of time), which corresponds to a noninteracting
lattice DOS of bandwidth $4$.

##  3 Results

Half-filled system. We first consider a half-filled two-orbital Hubbard model
with $U=10$, $J=2$, and $U_{cd}=3$. Because of the Hund coupling, the dominant
atomic configurations in equilibrium are the half-filled high-spin doublon
states with one electron in each orbital and parallel spins. The splitting
between the Hubbard bands is $U+J$. We choose the parameters $\mu$ and
$\Delta$ such that (i) the energies for $d$-electron addition and removal
(upper and lower Hubbard bands) are symmetric at $\pm 6$, and (ii) the core
removal energy is $E_{\text{core}}=-20$ [24]. For model (2), condition (i)
implies that the electron addition energy is
$E(3,1)-E(2,2)=2U-2J+2U_{cd}+\tfrac{\Delta}{2}-\mu\equiv 6$, while (ii)
implies $E(2,1)-E(2,2)=2U_{cd}-\Delta-\mu\equiv-20$. Here, $E(n_{d},n_{c})$ is
the energy of the corresponding state with filling $n_{d}$ and $n_{c}$ for the
$d$ and $c$ orbitals, respectively. These two conditions fix $\Delta=6.667$
and $\mu=19.33$. The corresponding equilibrium spectral functions are shown in
the left panel of Fig. 1. Apart from the main Hubbard bands, we recognize
shoulder and satellite structures at higher energies, which are associated
with local spin excitations (triplon insertion plus hopping to a neighboring
site, which leaves behind a low-spin doublon [25]).

Figure 1:  Left panel: equilibrium spectral functions of the half-filled model
with $U=10$, $J=2$, $U_{cd}=3$, $E_{\text{core}}=-20$. The coupling to a
Fermion bath with box-shaped DOS from $-25$ to $-15$ (black lines) and
coupling strength $v_{\text{bath}}=1$ [20] leads to a broadening of the
$c$-DOS. The main Hubbard bands in the $d$-DOS are at $\pm 6$. The middle
panel shows the local equilibrium RIXS response, normalized by
$E_{\text{probe}}^{2}$, and the right panel the normalized local RIXS signal
of the photo-doped system with $\omega_{\text{pump}}=12$ and
$t_{\text{probe}}-t_{\text{pump}}\approx 2$.

The middle panel shows the equilibrium RIXS signal $I_{\text{RIXS}}$ extracted
from the $D_{1,\sigma,\sigma^{\prime}}$ correlation functions, summed over the
spin channels, and normalized by $E_{\text{probe}}^{2}$. The probe pulse has
an envelope $f_{\text{probe}}(t)=\exp(-t^{2}/2)$, corresponding to a width of
$\Delta t\approx 2$ (and hence to an energy resolution of about $\Delta
E\approx\tfrac{1}{4}$), and an amplitude $E_{\text{probe}}=0.05$, which is
well within the perturbative regime where $I_{\text{RIXS}}\propto
E_{\text{probe}}^{2}$. The signal exhibits a single strong peak at
$\omega_{\text{in}}-|E_{\text{core}}|\approx 3$, which is associated with
processes where a core electron is excited into the upper Hubbard band and
decays back. Without the core-valence interaction $U_{cd}$, the corresponding
energy would be given by the electron insertion into the dominant subband of
the upper Hubbard band at
$\omega_{\text{in}}-|E_{\text{core}}|\approx\tfrac{U+J}{2}=6$, but because of
the core hole interaction, it is shifted to $\approx\tfrac{U+J}{2}-U_{cd}=3$.
Although it is not very clear on the scale of the plot, the signal peaks at a
value of $\omega_{\text{out}}$ which is almost independent of
$\omega_{\text{in}}$ within an energy range approximately equal to the width
of the Hubbard band. This implies that
$\omega_{\text{loss}}=\omega_{\text{in}}-\omega_{\text{out}}$ is approximately
proportional to $\omega_{\text{in}}$ in the corresponding energy range, i.e.
kinetic energy is transmitted to other electrons during the RIXS process
(fluorescent line). In addition, we see a very weak signal at
$\omega_{\text{out}}-|E_{\text{core}}|\approx-9$, corresponding to
$\omega_{\text{loss}}\approx U+J$. This signal originates from a RIXS process
that leaves a singlon-triplon pair in the system after the core hole is filled
($d$-$d$ excitation). The energy of $U+J$ corresponds to the energy difference
$E(3,2)+E(1,2)-2E(2,2)$ between two sites with a singlon-triplon pair and two
sites in their ground state.

The right panel of Fig. 1 shows $I_{\text{RIXS}}/E_{\text{probe}}^{2}$ of the
photo-doped system. Here, a few-cycle hopping modulation is applied to
generate a substantial density of singlons and triplons, before the RIXS
spectrum is probed. Figure 2 shows the evolution of the triplon probability,
$p(\text{triplon})=\langle
n_{1\uparrow}n_{1\downarrow}n_{2}+n_{1}n_{2\uparrow}n_{2\downarrow}-4n_{1\uparrow}n_{1\downarrow}n_{2\uparrow}n_{2\downarrow}\rangle$,
with $n_{\alpha}=n_{\alpha\uparrow}+n_{\alpha\downarrow}$, which is initially
very low, because the system is predominantly in half-filled high-spin states
with one electron per orbital. During the hopping modulation pulse with
frequency $\omega_{\text{pulse}}=U+J$ (blue line), the triplon population
increases even beyond the infinite temperature value of $\tfrac{1}{4}$, while
the probability of high-spin doublons decreases from $0.99$ to $0.17$. Due to
particle-hole symmetry, there is a corresponding increase in the density of
singlons. The RIXS probing pulse, which partially overlaps with the pump
pulse, measures the system immediately after the pump excitation (time delay
$t_{\text{probe}}-t_{\text{pump}}\approx 2$). We have also measured the RIXS
spectrum after a longer time delay $t_{\text{probe}}-t_{\text{pump}}\approx
4$, and obtained an almost identical result. This can be explained by the fact
that the transfer of kinetic energy to local spin excitations occurs quickly
(during the pump pulse) [26], while the life-time of the singlons and triplons
in this large-gap insulator is orders of magnitude longer than our simulation
times [27].

Figure 2:  Evolution of the density of triplons (red line) induced by a short
hopping modulation pulse (blue line). The state of the system immediately
after the pulse is probed by a RIXS pulse (green line, here for
$\omega_{\text{in}}=|E_{\text{core}}|$). The violet curve is proportional to
the photon count (for $\omega_{\text{out}}=|E_{\text{core}}|$).

As we can see in the right panel of Fig. 1, the photodoping generates several
prominent new peaks in $I_{\text{RIXS}}$ with $\omega_{\text{loss}}\approx 0$
(diagonal), and several weaker peaks corresponding to $d$-$d$ excitations with
$\omega_{\text{loss}}>0$ (lower-right triangle) and to $d$-$d$ de-excitations
with $\omega_{\text{loss}}<0$ (upper-left triangle). The presence of high- and
low-spin doublons in the photo-doped system now leads to three quasi-elastic
features associated with doublon-to-triplon excitations, at energies
$E(3,1)-E(2,2)=3$, $3-J=1$, $3-3J=-3$. The new quasi-elastic peaks at
$\omega_{\text{in,out}}-|E_{\text{core}}|\approx-6$, $-4$ and $0$ correspond
to RIXS processes where an electron is excited from the core level to a
singlon, and eventually decays back. The peak near $-6$ can be associated with
an intermediate high-spin doublon state (e. g.
$|\\!\uparrow,\uparrow\rangle$), the peak near $-4$ with an intermediate low-
spin doublon state with one electron in each orbital (e. g.
$|\\!\uparrow,\downarrow\rangle$), and the peak near $0$ with an intermediate
doublon state with a doubly occupied and an empty orbital (e. g.
$|0,\uparrow\downarrow\rangle$). Indeed, since the $U_{cd}$\- and
$\mu$-contributions to the local energy do not change in these processes, the
corresponding excitation energies $E(2,1)-E(1,2)$ are
$U-3J+\tfrac{3}{2}\Delta=14$, $U-2J+\tfrac{3}{2}\Delta=16$ and
$U+\tfrac{3}{2}\Delta=20$, which upon subtraction of $|E_{\text{core}}|$
yields the observed peak positions. The excitation of a core electron to a
triplon state produces a change in the $U_{cd}$ contribution, and the
corresponding excitation energy $E(4,1)-E(3,2)$ becomes
$3U-5J-2U_{cd}+\tfrac{3}{2}\Delta=24$. Hence, these excitations explain the
prominent peak at $\omega_{\text{in,out}}-|E_{\text{core}}|\approx 4$, which
partially overlaps with the peak at $\approx 3$ corresponding to excitations
to doublon states. The latter is substantially weaker in the photo-doped
system than in equilibrium, because the high-spin doublon states no longer
dominate. There is also a small density of empty sites in the photo-doped
state, which results in a weak quasi-elastic feature near $E(1,1)-E(0,2)=13$,
i.e. $\omega_{\text{in,out}}-|E_{\text{core}}|\approx-7$. We summarize the
different quasi-elastic features in the half-filled equilibrium and photo-
doped systems in Tab. 1.

relevant transitions | half-filled | quarter-filled
---|---|---
equilibrium: |  |
$(2^{h},2)\leftrightarrow(3,1)$ | 3 |
$(1,2)\leftrightarrow({\bf 2},1)$ |  | $-1,1,5$
photo-doped: |  |
$(0,2)\leftrightarrow(1,1)$ | $-7$ | $-2$
$(1,2)\leftrightarrow({\bf 2},1)$ | $-6,-4,0$ | $-1,1,5$
$({\bf 2},2)\leftrightarrow(3,1)$ | $-3,1,3$ | $2,6,8$
$(3,2)\leftrightarrow(4,1)$ | $4$ | $9$
Table 1: Energies $\omega_{\text{in,out}}-|E_{\text{core}}|$ of the quasi-
elastic features in the half-filled and quarter-filled system. The left
(right) number in the brackets indicates the filling of the $d$ ($c$) orbital.
$2^{h}$ denotes the high-spin doublon state, and ${\bf 2}$ the three different
types of doublons.

Figure 3:  Left panel: equilibrium spectral functions of the quarter-filled
model with $U=10$, $J=2$, $U_{cd}=3$, $E_{\text{core}}=-20$. The coupling to a
Fermion bath with box-shaped DOS from $-25$ to $-15$ (black lines) and
coupling strength $v_{\text{bath}}=1$ [20] leads to a broadening of the
$c$-DOS. The lower Hubbard band (singly occupied sites) in the $d$-DOS is at
$-2$, while the upper Hubbard band (doubly occupied sites) has three subbands
at energies $2$, $4$ and $6$. The middle panel shows the local equilibrium
RIXS response, normalized by $E_{\text{probe}}^{2}$, and the right panel the
normalized local RIXS signal of the photo-doped system with
$\omega_{\text{pump}}=10$ and $t_{\text{probe}}-t_{\text{pump}}\approx 2$.

The weaker RIXS features on the $\omega_{\text{loss}}>0$ side (lower-right
triangle) can be associated with processes involving the neighboring sites.
For example the small peak at $\omega_{\text{in}}-|E_{\text{core}}|\approx-4$,
$\omega_{\text{out}}-|E_{\text{core}}|\approx-6$ can be explained by a process
where a singlon is converted into a low-spin doublon with interaction energy
$U-2J$, which is then (via a second-order hopping process) flipped into a
high-spin doublon with energy $U-3J$, before an electron decays to fill the
core hole:
$(|0,\downarrow\rangle|\\!\uparrow\downarrow\rangle)$$\xrightarrow{\omega_{\text{in}}}(|\\!\\!\uparrow,\downarrow\rangle|\\!\\!\downarrow\rangle)\xrightarrow{\text{ex}}(|\\!\\!\uparrow,\uparrow\rangle|\\!\\!\downarrow\rangle)\xrightarrow{\omega_{\text{out}}}(|0,\uparrow\rangle|\\!\\!\uparrow\downarrow\rangle)$,
where the configurations $(|d\rangle|c\rangle)$ indicate the valence
$(|d\rangle)$ and core $(|c\rangle)$ configurations on the probe site. The
energy loss $\omega_{\text{loss}}\approx J$ is either converted into a local
spin excitation on the neighboring site, or converted into kinetic energy.
Similarly, the feature at $\omega_{\text{in}}-|E_{\text{core}}|\approx 0$,
$\omega_{\text{out}}-|E_{\text{core}}|\approx-6$ corresponds to an
intermediate state with energy $U$ (doublon with both electrons in the same
orbital). The much weaker feature at
$\omega_{\text{in}}-|E_{\text{core}}|\approx 0$,
$\omega_{\text{out}}-|E_{\text{core}}|\approx-4$ indicates that the states
involved in the corresponding second-order processes are less populated or
characterized by weaker exchange couplings.

$d$-$d$ excitation peaks at $\omega_{\text{in}}-|E_{\text{core}}|\approx 3$
are related to excitations of a core electron to a high-spin doublon state
(see Tab. 1). The resulting triplon can hop to a neighboring site, leaving
behind doublons in various configurations. The decay to a singlon plus filled
core hole then produces photons at
$\omega_{\text{out}}-|E_{\text{core}}|\approx 0$, $-4$, $-6$. In the photo-
doped state, the intermediate step (hopping of triplon) is aided by the
presence of low-spin doublons on neighboring sites, which reduces the amount
of kinetic energy that needs to be absorbed or emitted. Alternatively, the
triplon can change its configuration in a second-order process and then decay
to a low-spin doublon plus a filled core hole. In this case, the emitted
energy peaks around $\omega_{\text{out}}-|E_{\text{core}}|\approx 3-J=1$ and
$\approx 3-3J=-3$, while the intermediate step is again aided by the presence
of neighboring low-spin doublons. The superposition of these two series of
peaks explains the observed loss features. Note that the peak at
$\omega_{\text{out}}-|E_{\text{core}}|=-9$, which appeared in the equilibrium
spectrum, becomes very weak in the photo-doped state. This is because pairs of
neighboring high-spin doublons are unlikely in the photo-doped system.

We next turn our attention to the features in the upper left triangle,
corresponding to $\omega_{\text{loss}}<0$. In a photo-doped system, energy can
be gained in processes involving singlon-triplon recombination or local state
transitions. The gain features at
$\omega_{\text{in}}-|E_{\text{core}}|\approx-6$ and
$\omega_{\text{out}}-|E_{\text{core}}|\approx-4$, $0$, as well as
$\omega_{\text{in}}-|E_{\text{core}}|\approx-4$ and
$\omega_{\text{out}}-|E_{\text{core}}|\approx 0$ can be explained as energy
absorptions of $J$, $2J$ and $3J$ from local state transitions in second order
hoppings between neighboring doublons. These become possible because low-spin
doublons are already present before the RIXS pulse. The strong feature near
$\omega_{\text{out}}-|E_{\text{core}}|\approx 3$ and
$\omega_{\text{in}}-|E_{\text{core}}|$ between $-4$ and $-3$ originates from
emission processes that leave behind a high-spin doublon state. The initial
state is either a doublon with local energy $U$
($\omega_{\text{in}}-|E_{\text{core}}|\approx-3$) or a singlon. In the latter
case, the doublon state after the RIXS excitation has local energy $U-2J$
(hence $\omega_{\text{in}}-|E_{\text{core}}|\approx-4$), and before the de-
excitation, a triplon hops in from a neighboring site:
$(|0,\downarrow\rangle|\\!\\!\uparrow\downarrow\rangle)\xrightarrow{\omega_{\text{in}}}(|\\!\\!\uparrow,\downarrow\rangle|\\!\\!\downarrow\rangle)\xrightarrow{\text{hop}}(|\\!\\!\uparrow\downarrow,\uparrow\rangle|\\!\\!\downarrow\rangle)\xrightarrow{\omega_{\text{out}}}(|\\!\\!\downarrow,\downarrow\rangle|\\!\\!\uparrow\downarrow\rangle)$.
Again, such processes become possible because singlons and triplons are
present with high probability in the photo-doped system.

Quarter-filled system. As a second example, we consider the quarter-filled
two-band model with $U=10$, $J=2$, $U_{cd}=3$, $E_{\text{core}}=-20$, which
has a very different equilibrium DOS, as illustrated in the left panel of Fig.
3. Here, in order to place the core level at $-20$, and the highest energy
subband of the upper Hubbard band at $8$, we choose the electron removal and
addition energies as $-\mu+U_{cd}-\Delta\equiv-20$ and
$-\mu+2U_{cd}+U+\tfrac{\Delta}{2}\equiv 8$, which gives $\Delta=10$ and
$\mu=13$. With this choice, the lower Hubbard band (representing singly
occupied states) is centered at $\approx-2$, and the subbands of the upper
Hubbard band (corresponding to high-spin and low-spin doublon states) at
$\approx 2$, $4$, and $8$, see left panel of Fig. 3. In the atomic limit, the
energy of creating a high-spin doublon becomes
$E(2,1)-E(1,2)=U-3J+\tfrac{3}{2}\Delta=19$, and that of creating a low-spin
doublon with electrons in different orbitals (the same orbital)
$U-2J+\tfrac{3}{2}\Delta=21$ ($U+\tfrac{3}{2}\Delta=25$). Hence, in the
equilibrium RIXS spectrum we expect dominant features around
$\omega_{\text{in/out}}-|E_{\text{core}}|\approx-1$, $1$, $5$, which are
associated with the creation and decay of these different types of doublons.
This is indeed approximately the case as shown in the middle panel of Fig. 3.
Again, the horizontal elongation of these dominant features indicates energy
loss through scattering within a given Hubbard-subband, and hence a
fluorescence-like behavior.

Because the energy gaps are smaller than in the half-filled case, prominent
$d$-$d$ excitation peaks are visible already in the equilibrium RIXS spectrum.
These energy loss features appear near
$\omega_{\text{out}}-|E_{\text{core}}|\approx-2$ and $-1$. The signal at $-2$
can be associated with the transition from a singly occupied $d$ site with
core hole to an empty $d$ site with filled core hole. This type of emission
occurs in processes that leave behind a doublon-hole pair in the singly
occupied background:
$(1,2)(1,2)\xrightarrow{\omega_{\text{in}}}(2,1)(1,2)\xrightarrow{\text{hop}}(1,1)(2,2)\xrightarrow{\omega_{\text{out}}}(0,2)(2,2)$.
Here, the two brackets correspond to neighboring sites, the first entry to the
$d$ occupation and the second entry to the $c$ occupation. The first arrow
represents the excitation by the probe pulse, the last arrow the de-
excitation, and the middle arrow a hopping process, which costs an energy
$\Delta E=U_{cd}-nJ$, where $n$ depends on the specific doublon configuration.
In the case of the signal at $-2$, this energy is absorbed from the system (e.
g. from the kinetic energy of the carriers) and does not show up in
$\omega_{\text{out}}$. If on the other hand, we combine the last two processes
above into a single “hopping plus emission” process, the energy $\Delta E$
must be subtracted from $\omega_{\text{out}}$, which explains the weaker
signals at $\approx-2-U_{cd}=-5$ and $-2-(U_{cd}-J)=-3$. Photons with
$\omega_{\text{out}}-|E_{\text{core}}|\approx-1$ are emitted by the de-
excitation of a high-spin doublon. For example, the peak around
$\omega_{\text{in}}-|E_{\text{core}}|\approx 5$ and
$\omega_{\text{out}}-|E_{\text{core}}|\approx-1$ originates from a RIXS
process where the incoming pulse creates a low-spin doublon with interaction
energy $U$, which in an intermediate step is converted into a high-spin
doublon. This intermediate step releases an energy $\Delta E=3J=6$, which is
converted into kinetic energy or additional doublon-holon excitations via the
Hund-excitation analogue of impact ionization [28].

The right hand panel of Fig. 3 shows the RIXS spectrum after a photo-doping
pulse. Here, we use the same set-up as illustrated in Fig. 2, but with a pump
pulse frequency $\omega_{\text{pump}}=U$. As in the case of the half-filled
system, the photo-doping results in additional peaks on the diagonal (quasi-
elastic processes), as well as new gain and loss features in the upper left
and lower-right triangle. The most prominent new peak appears at
$\omega_{\text{in,out}}\approx-2$. This feature corresponds to the excitation
of a core electron to an initially empty site, and the decay of this electron
back to the core. Indeed, the corresponding excitation energy is
$E(1,1)-E(0,2)=U_{cd}+\tfrac{3}{2}\Delta=18$, so that
$\omega_{\text{in,out}}-|E_{\text{core}}|=-2$. There are also doublons in the
initial photo-doped state. Electron addition to a high-spin doublon yields
$E(3,1)-E(2,2)=2U-2J-U_{cd}+\tfrac{3}{2}\Delta=28$. Hence, three features on
the diagonal at energies $\omega_{\text{in,out}}-|E_{\text{core}}|\approx 8$,
$8-2J=6$, $8-3J=2$ can be associated with the creation and decay of triplons
(Tab. 1). The prominent peak near
$\omega_{\text{in,out}}-|E_{\text{core}}|\approx 9$ has a different origin: it
can be associated with the creation and decay of quadruplons (addition of a
core electron to a triplon), since
$E(4,2)-E(3,2)=3U-5J-2U_{cd}+\tfrac{3}{2}\Delta=29$. Indeed, our photo-doping
pulse is strong enough that it not only creates doublons, but also a non-
negligible density of triplons, see the top panel of Fig. 4.

Figure 4:  Top panel: Time evolution of the doublon and triplon probabilities
in the photo-doped quarter-filled system. Here, we label the type of doublon
by the corresponding local energy. Bottom panels: Time dependent change in
$I_{\text{RIXS}}/E_{\text{probe}}^{2}$. The left panel shows the difference
for probe pulses centered at $t_{\text{probe}}=4.5$ and $3.5$ and the right
panel the difference for probe pulses centered at $t_{\text{probe}}=5.5$ and
$4.5$.

The features associated with $d$-$d$ excitations ($\omega_{\text{loss}}>0$) in
the equilibrium spectrum get weaker in the photo-doped state, because the
corresponding initial states (singlons) are less probable, and the production
of doublon-hole pairs is suppressed if there is already a high density of
doublons and holons in the system. A closer inspection further reveals that
the dominant loss peak near $\omega_{\text{in}}-|E_{\text{core}}|\approx 1$
seems to have shifted to a lower energy. We interpret this peak as a holon-
assisted feature. The RIXS pulse creates a low-energy doublon, which is
converted into a singlon with the help of a neighboring holon. In this case,
the intermediate step releases $\Delta E=U_{cd}-(U-3J)=-3$, which can be
converted into kinetic energy or Hund excitations. The filling of the core
hole from a singlon then emits light with
$\omega_{\text{out}}-|E_{\text{core}}|\approx-2$.

Gain features with $\omega_{\text{loss}}<0$ can be associated with doublon or
triplon assisted intermediate processes. For example, the peak at
$\omega_{\text{in}}-|E_{\text{core}}|\approx-2$ and
$\omega_{\text{out}}-|E_{\text{core}}|\approx 5$ corresponds to the creation
of a singlon on an empty site by the RIXS pulse. This singlon is then
converted into a low-spin same-orbital doublon through electron exchange with
a neighboring doublon. This doublon subsequently decays into a singlon by
filling the core hole. Energy is gained in this process from the recombination
of a pump-induced doublon-hole pair. There are additional signals at the same
$\omega_{\text{out}}$ from processes with
$\omega_{\text{in}}-|E_{\text{core}}|\approx-1$ and $1$ (one of them
overlapping with the former peak). Here, a high-spin or low-spin (two-orbital)
doublon is created by the RIXS pulse and converted into a same-orbital doublon
via a second-order hopping process with a neighboring doublon. In this case,
energy is gained from local state transitions (Hund de-excitations) of photo-
induced doublons. The gain features with
$\omega_{\text{out}}-|E_{\text{core}}|\approx 8$ on the one hand originate
from triplon assisted processes. In this case, a doublon is created by the
RIXS pulse near $\omega_{\text{in}}\approx-1$, $1$, $5$ and converted into a
triplon by the hopping of a neighboring triplon. On the other hand, RIXS
pulses can produce triplons from low-spin doublons at
$\omega_{\text{in}}\approx 2$, $6$. The decay of the triplon into a high-spin
doublon plus filled core hole emits radiation at
$\omega_{\text{out}}-|E_{\text{core}}|\approx 8$. Note that there are no
apparent gain features associated with the peak at
$\omega_{\text{out}}-|E_{\text{core}}|\approx 9$. Such features would require
the generation of quadruplons as intermediate states,
$(2,2)\xrightarrow{\omega_{\text{in}}}(3,1)\xrightarrow{\text{hop}}(4,1)\xrightarrow{\omega_{\text{out}}}(3,2)$,
which is very unlikely if the density of photo-doped quadruplons is small.

Finally, we illustrate in the bottom panels of Fig. 4 the ultrafast evolution
of the RIXS signal after the photo-doping pulse in the quarter filled system.
The left panel shows the difference between the signals measured with probe
pulses centered at $t_{\text{probe}}=4.5$ and $3.5$, and the right panel the
difference of the signals for $t_{\text{probe}}=5.5$ and $4.5$. We see that
the quasi-elastic signals at $\omega_{\text{in,out}}\approx-2$ and $9$
increase, which indicates that the probability of holes and triplons in the
system continues to increase after the pump at least until $t\approx 4.5$.
This is directly confirmed in the top panel, which plots the evolution of the
triplon population (the population of holes similarly increases). Also the
quasi-elastic features associated with high-spin ($U-3J$) and low-spin both
orbital ($U-2J$) doublons at $\omega_{\text{in,out}}\approx 8$ and $6$ gain
some weight. This observation is again consistent with the direct measurement
of the corresponding doublon populations in the upper panel. On the other
hand, the quasi-elastic signal and the gain and loss features near
$\omega_{\text{in}}\approx 2$ are suppressed, which indicates that RIXS
processes involving low-spin same-orbital ($U$) doublons become less likely
with increasing pump-probe delay. Indeed, these are the doublons created
predominantly by the resonant pump pulse with $\omega_{\text{pump}}=U$, and
they decay into the lower-energy doublons through scattering events (see red
line in the upper panel). This example shows how the analysis of the time-
resolved RIXS signal allows to extract information on the evolution of the
different local state populations during and after the photo-excitation.

##  4 Conclusions

We have calculated the equilibrium and nonequilibrium RIXS spectra of the
half-filled and quarter-filled two-orbital Hubbard model using a
nonequilibrium-DMFT based formalism. The main idea behind this scheme is to
avoid the computation and storage of four-time correlation functions by
explicitly simulating the RIXS excitation pulse and measuring the self-energy
of the emitted photons. This paper provides a proof-of-principle and
demonstrates the feasibility of this approach in a multi-orbital context. We
showed that in equilibrium our method captures $d$-$d$ excitations as well as
fluorescent lines, similar to the DMFT based studies employing a
configuration-interaction solver and the Kramers-Heisenberg formalism [16,
17]. After a strong photo-excitation, we observed additional peaks in the RIXS
signal associated with the presence of photo-doped charge carriers, as well as
new or modified gain and loss features connected to intermediate states and
transitions involving these photo-carriers. We furthermore demonstrated how
time-dependent changes in the RIXS spectrum can be tracked by shifting the
probe pulse relative to the pump and how they reveal the evolution of
different local states.

Our current implementation lacks momentum resolution and involves a certain
number of approximations. For example, the core-hole lifetime is controlled by
an electron bath coupled to the core level. As shown in Ref. [20], this is
essentially equivalent to the use of the life-time parameter $\Gamma$ in
previous equilibrium [16, 17] and nonequilibrium [14, 19] RIXS studies.
Potentially more severe is the use of the NCA in the solution of the DMFT
impurity model and the calculation of the photon self-energy. While
qualitatively correct results can be expected in the Mott regime[20],
systematic tests with beyond-NCA approaches will be needed to fully understand
the implications of this approximation on the RIXS signal. It is also
important to point out that while we have used here an implementation with
explicit core level and classical light, there also exist alternative
nonequilibrium DMFT implementations where the core level is shifted to
$E_{\text{core}}=-\infty$, or where the outgoing photons are treated quantum
mechanically [20]. Testing these alternative implementations on multi-orbital
Hubbard models will be an interesting topic for future investigations.

###### Acknowledgements.

We thank C. Monney, S. Johnson and H. Strand for helpful discussions. P. W.
acknowledges support from ERC Consolidator Grant 724103 and M. E. from ERC
Starting Grant No. 716648. S. J. acknowledges support from the National
Science Foundation under Grant No. DMR-1842056. The calculations were run on
the beo04 cluster at the University of Fribourg.

## References

  * [1] Ament L., van Veenendaal M., Devereaux T., Hill J., and van den Brink J., Resonant inelastic x-ray scattering studies of elementary excitations, Rev. Mod. Phys. 83, 705 (2011).
  * [2] Hill J., Kao C.-C., Caliebe W., Matsubara M., Kotani A., Peng J., and Greene P., Resonant Inelastic X-Ray Scattering in Nd2CuO4, Phys. Rev. Lett. 80, 4967 (1998).
  * [3] Braicovich L., Ament L. J. P. , Bisogni V., Forte F., Aruta C., Balestrino G., Brookes N. B., De Luca G. M., Medaglia P. G., Miletto Granozio F., Radovic M., Salluzzo M., van den Brink J., and Ghiringhelli G., Dispersion of Magnetic Excitations in the Cuprate La2CuO4 and CaCuO2 Compounds Measured Using Resonant X-Ray Scattering, Phys. Rev. Lett. 102, 167401 (2009).
  * [4] Schlappa, J., Kumar, U., Zhou, K. J., Singh, S., Mourigal, M., Strocov, V. N., Revcolevschi, A., Patthey, L., Rønnow, H. M. , Johnston, S. and T. Schmitt, T., Probing multi-spinon excitations outside of the two-spinon continuum in the antiferromagnetic spin chain cuprate Sr2CuO3. Nat. Commun. 9, 5394 (2018).
  * [5] Chen, C. C., Moritz, B., Vernay, F., Hancock, J. N., Johnston, S., Jia, C. J., Chabot-Couture, G., Greven, M., Elfimov, I., Sawatzky, G. A., and Devereaux, T. P., Unraveling the Nature of Charge Excitations in La2CuO4 with Momentum-Resolved Cu $K$-Edge Resonant Inelastic X-Ray Scattering. Phys. Rev. Lett. 105, 177401 (2010).
  * [6] Schlappa J., Wohlfeld, K., Zhou, K. J., Mourigal, M., Haverkort, M. W., Strocov, V. N., Hozoi, L., Monney, C., Nishimoto, S., Singh, S., Revcolevschi, A., Caux, J.-S., Patthey, L., Rønnow, H. M., van den Brink, J. and T. Schmitt, T., Spin-orbital separation in the quasi-one-dimensional Mott insulator Sr2CuO3, Nature 485, 82 (2012).
  * [7] Ament L. J. P., van Veenendaal M., and van den Brink J., Determining the electron-phonon coupling strength from Resonant Inelastic X-ray Scattering at transition metal $L$-edges, Europhys. Lett. 95, 27008 (2011).
  * [8] Johnston, S., Monney, C., Bisogni, V., Zhou, K.-J., Kraus, R., Behr, G., Strocov, V. N., Málek, J., Drechsler, S.-L., Geck, J., Schmitt, T. and van den Brink, J. Electron-lattice interactions strongly renormalize the charge-transfer energy in the spin-chain cuprate Li2CuO2, Nat. Commun. 7, 10563 (2016).
  * [9] Chaix L., Ghiringhelli, G., Peng, Y. Y., Hashimoto, M., Moritz, B., Kummer, K., Brookes, N. B., He, Y., Chen, S., Ishida, S., Yoshida, Y., Eisaki, H., Salluzzo, M., Braicovich, L., Shen, Z.-X., Devereaux, T. P. and Lee, W.-S., Dispersive charge density wave excitations in Bi2Sr2CaCu2O8+δ, Nat. Phys. 13, 952 (2017).
  * [10] Dean M. et al., Ultrafast energy- and momentum-resolved dynamics of magnetic correlations in the photo-doped Mott insulator Sr2IrO4, Nature Mat. 15, 601 (2016).
  * [11] Mitrano M. et al., Ultrafast time-resolved x-ray scattering reveals diffusive charge order dynamics in La2-xBaxCuO4, Science Advances 5, eaax3346 (2019).
  * [12] Mitrano, M. and Wang, Y., Light-Driven Quantum Materials with Ultrafast Resonant Inelastic X-Ray Scattering. Commun. Phys. 3, 184 (2020).
  * [13] Kramers H. and Heisenberg W., Über die Streuung von Strahlung durch Atome, Z. Phys. 31, 681 (1925).
  * [14] Chen Y., Wang Y., Jia C., Moritz B., Shvaika A. M., Freericks J., and Devereaux T., Theory for time-resolved resonant inelastic x-ray scattering, Phys. Rev. B 99, 104306 (2019).
  * [15] Georges A., Kotliar G., Krauth W., and Rozenberg M. J., Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions. Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys., 130 68 (1996).
  * [16] Hariki A., Winder M., and Kunes J., Continuum Charge Excitations in High-Valence Transition-Metal Oxides Revealed by Resonant Inelastic X-Ray Scattering, Phys. Rev. Lett. 121, 126403 (2018).
  * [17] Hariki A., Winder M., Uozumi T., and Kunes J., LDA+DMFT approach to resonant inelastic x-ray scattering in correlated materials, Phys. Rev. B 101, 115130 (2020).
  * [18] Ghiringhelli G., Piazzalunga A., Dallera C., Schmitt T., Strocov V., Schlappa J., Patthey L., Wang X., Berger H., and Grioni M., Observation of Two Nondispersive Magnetic Excitations in NiO by Resonant Inelastic Soft-X-Ray Scattering, Phys. Rev. Lett. 102, 027401 (2009).
  * [19] Wang Y., Chen Y., Jia C., Moritz B., and Deveraux T., Time-Resolved Resonant Inelastic X-Ray Scattering in a Pumped Mott Insulator, Phys. Rev. B 101, 165126 (2020).
  * [20] M. Eckstein and P. Werner, Simulation of time-dependent resonant inelastic X-ray scattering with non-equilibrium dynamical mean-field theory, arxiv (2020).
  * [21] Aoki H., Tsuji N., Eckstein M., Kollar M., Oka T., and Werner P., Nonequilibrium dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014).
  * [22] Keiter H. and Kimball J. C., A revised diagrammatic technique for the degenerate Anderson model, Int. J. Magn. 1, 233 (1971).
  * [23] Eckstein M. and Werner P., Nonequilibrium dynamical mean-field calculations based on the noncrossing approximation and its generalizations, Phys. Rev. B 82, 115115 (2010).
  * [24] A realistic core level would be at much lower energy. Here, we choose a shallow core level to enable calculations with a time step $\Delta t=0.01$ that nevertheless resolve the short timescale $\hbar/|E_{\text{core}}|$. For larger $E_{\text{core}}$, one should write the core-valence Hamiltonian in the rotating wave approximation [20]. However, already for $E_{\text{core}}=20$, counter-rotating terms give only a small contribution to the result.
  * [25] Lysne M., Murakami Y. and Werner P., Signatures of bosonic excitations in high-harmonic spectra of Mott insulators, Phys. Rev. B 101, 195139 (2020).
  * [26] Strand H., Golez D., Eckstein M., and Werner P., Hund’s coupling driven photocarrier relaxation in the two-band Mott insulator, Phys. Rev. B 96, 165104 (2017).
  * [27] Eckstein M. and Werner P., Thermalization of a pump-excited Mott insulator, Phys. Rev. B 84, 035122 (2011).
  * [28] Petocchi F., Beck S., Ederer C., and Werner P., Hund excitations and the efficiency of Mott solar cells, Phys. Rev. B 100, 075147 (2019).