# Dynamical Triplet Unravelling: A quantum Monte Carlo algorithm for
reversible dynamics

Romain Chessex<EMAIL_ADDRESS>Massimo Borrelli  Hans Christian Öttinger
Polymer Physics, Department of Materials, ETH Zürich, CH-8093 Zürich,
Switzerland

###### Abstract

We introduce a quantum Monte Carlo method to simulate the reversible dynamics
of correlated many-body systems. Our method is based on the Laplace transform
of the time-evolution operator which, as opposed to most quantum Monte Carlo
methods, makes it possible to access the dynamics at longer times. The Monte
Carlo trajectories are realised through a piece-wise stochastic-deterministic
reversible evolution where free dynamics is interspersed with two-process
quantum jumps. The dynamical sign problem is bypassed via the so-called
deadweight approximation, which stabilizes the many-body phases at longer
times. We benchmark our method by simulating spin excitation propagation in
the XXZ model and dynamical confinement in the quantum Ising chain, and show
how to extract dynamical information from the Laplace representation.

††preprint: APS/123-QED

##  I Introduction

Recent years have witnessed an ever-growing interest in the many-body dynamics
of strongly correlated quantum systems [1, 2]. Technical advancements in
quantum state preparation, control and measurement have made it possible to
observe the real-time evolution of typical systems in condensed matter physics
at very low temperatures. Several experimental platforms, ranging from cold
atoms in optical lattices [3, 4, 5, 6, 7], to trapped ions [8, 9] and
superconducting circuits [10, 11], have emerged as quantum simulators and
different many-body Hamiltonians have been realized. The plethora of phenomena
that have been investigated is very rich and includes many-body equilibration
versus localization [12, 13, 14], many-body dissipation, dynamics of quantum
correlations [15, 16, 17], quenched many-body dynamics, and dynamical phase
transitions [18, 19, 20].

Parallel to the abundance of experimental investigations, numerical methods
[21, 22, 19] have been developed to simulate condensed matter systems
efficiently, including molecular and atomic gases in optical lattices. At the
intersection between many-body physics and quantum information theory, tensor
networks have been extensively explored in 1d systems, and nowadays, they are
considered as a universal representation of low-energy, weakly entangled
quantum states. Within this framework, time-dependent density matrix
renormalization group [23, 24, 25, 26] has proven very effective in handling
large particle systems and has established itself as a reliable standard tool.
For higher dimensional systems and longer simulation times, nonequilibrium
dynamical mean field theory [27, 28] has been successfully applied to study
spectroscopy [29], nonlinear optics [30] and transport in solid-state physics
[31] as well as quenches [32] and relaxation dynamics in cold atoms in optical
lattices [33, 34]. Finally, other computationally efficient methods have found
a multitude of applications in quantum chemistry, with prominent examples
being time-dependent Hartree-Fock [35], multiconfigurational time-dependent
Hartree-Fock [36], in which some electron-electron correlations are also
included, and density functional theory [37].

Most of the above techniques rely explicitly on symmetries in the Hamiltonian
to cope with the large size of the Hilbert space. As opposed to deterministic
techniques, quantum Monte Carlo methods (QMC) sample stochastic trajectories
in a configuration space, onto which the original Hilbert space is directly
mapped, and have proven very reliable in a number of different applications
[38, 39, 40, 41, 42, 43]. Although QMC methods were originally formulated to
study ground state and equilibrium properties, dynamical adaptations have also
been put forward that are based on stochastic sampling of diagrammatic
expansions of many-body Green’s functions defined along a Keldysh contour [40,
44, 45]. All these methods, however, have been limited in accessing long-time
behaviour by the well-known dynamical sign problem. Notable exceptions are the
inchworm algorithm presented in [46] and real-time full configuration
interaction QMC algorithm [39]. The first technique can access medium-to-long-
time $t\sim 10$ dynamical properties, due to corrections of short-time
diagrams, whereas the second can allow for simulation times up to $t\sim 40$,
although at the expense of probability conservation.

Here, we present a new QMC algorithm that allows us to simulate the reversible
dynamics of a quantum many-body system at long times. In some technical
aspects, this method draws from the fixed point quantum Monte Carlo method
(FPQMC), which was first introduced in [47] for ground state calculations. The
reversible dynamics, as dictated by the von Neumann equation, is reformulated
in the Laplace representation and unraveled by stochastic trajectories. These
are generated by a piece-wise deterministic-stochastic evolution where the
free part of the many-body Hamiltonian is solved exactly, whilst particle-
particle interactions are simulated using two-process quantum jumps. The QMC
walkers sampling the configuration space are called triplets and consist of
two state vectors, and a statistical weight. Depending on this weight, a
triplet might spawn a new triplet via a two-process quantum jump or be removed
from the simulation altogether. Two key ingredients are introduced. First, the
deadweight approximation, which stabilizes the phases originating from fast
oscillating exponentials that are responsible for the dynamical sign problem.
Second, as previously formulated in FPQMC and relying on dynamic norm,
importance sampling guarantees fast convergence of the time-dependent
quantities under scrutiny. Near-exact trace conservation of the density matrix
is guaranteed for the whole duration of the simulation. We benchmark our
algorithm against two case-studies on quantum lattices i) dynamics of a spin
excitation following a quench in the Heisenberg XXZ model and ii) excitation
confinement in the quantum Ising model [48, 49, 50]. In both cases we find
results in excellent agreement with previous established literature.

This manuscript is organized as follows; in Sec. II, we present the
theoretical background for our numerical method and introduce some physical
quantities of interest. In Sec. III, we break down the algorithm itself and
explain in detail all its major components. In Sec. IV, we test our method for
the case studies mentioned earlier and analyze the results. Finally, we draw
some conclusion and illustrate some open perspectives in Sec. V.

##  II Theoretical basis

###  II.1 Dynamics in the Laplace space

We focus on many-body systems whose Hamiltonian $H$ can be split into a free
part and an interacting part, that is $H=H^{\text{free}}+H^{\text{int}}$. The
adjectives free and interacting do not necessarily refer to the actual kinetic
energy versus interactions within the Hamiltonian, but rather to whether the
exact eigenstates are known a priori, which is assumed to be always the case
for $H^{\text{free}}$. The time evolution is dictated by the von Neumann
equation

$\frac{d}{dt}\rho=\mathcal{L}\rho\equiv-i\left[H,\rho\right],$ (1)

where $[\cdot,\cdot]$ denotes the commutator and $\rho$ the density matrix of
the many-body system. If the total Hamiltonian is time independent, the formal
solution of Eq. (1) reads

$\rho_{t}=\mathcal{E}_{t}\rho_{0}=e^{t\mathcal{L}}\rho_{0},$ (2)

for some initial state $\rho_{0}$. To use the framework developed in [47], we
transform the time-evolution superoperator $\mathcal{E}_{t}$ in the Laplace
domain, that is

$\mathcal{R}_{s}=\int_{0}^{\infty}\mathcal{E}_{t}e^{-st}\text{d}t.$ (3)

We define the superoperators

$\mathcal{L}^{\text{free}}\rho=-i[H^{\text{free}},\rho],\qquad\mathcal{L}^{\text{int}}\rho=-i[H^{\text{int}},\rho],$
(4)

generating the dynamics described by the superoperators
$\mathcal{R}^{\text{free}}_{s}$ and $\mathcal{R}^{\text{int}}_{s}$,
respectively. By first Laplace transforming the formal solution of Eq. (1) and
then applying the geometric series (see [51] for details), one arrives at the
following expression for the Laplace transform of the total evolution
superoperator $\mathcal{R}_{s}$

$\mathcal{R}_{s}=\sum_{m=0}^{\infty}r^{m}\left[\mathcal{R}_{s+r}^{\text{free}}\left(1+\frac{\mathcal{L}^{\text{int}}}{r}\right)\right]^{m}\mathcal{R}_{s+r}^{\text{free}},$
(5)

with $r>0$. We refer to Eq. (5) as the magical formula because it unifies a
perturbative expansion and a numerical scheme in a single equation. If $r=0$,
the perturbative expansion is recovered, since the $m^{\text{th}}$ term in the
sum corresponds to the $m^{\text{th}}$ order in perturbation theory.
Conversely, if $r>0$, the sum can be interpreted as the time-evolution
operator of a numerical integration scheme with a $1/r$ time-step. Hence,
truncating the magical formula at some $M_{\text{trunc}}$ term sets a natural
time limit for the numerical integration, which is
$t_{\text{max}}\sim\frac{1}{s_{\text{min}}}\sim\frac{M_{\text{trunc}}}{r}$.

The Laplace transform of the time-evolution superoperator (3) does not
guarantee trace conservation by itself. However, noting that

$\text{Tr}(\mathcal{R}_{s}\rho_{0})=\int_{0}^{\infty}\text{Tr}(\rho_{t})e^{-st}\text{d}t=\frac{1}{s},$
(6)

one recovers the correct normalization for the quantity $\tilde{\rho}_{s}$
defined by

$\tilde{\rho}_{s}=s\mathcal{R}_{s}\rho_{0}.$ (7)

Even though $\tilde{\rho}_{s}$ possess all properties of a density matrix, we
have to be careful when interpreting it as a physical density matrix because
it is constructed from an integration over the whole time domain. In the limit
$s\to\infty$, the density matrix is associated to a physical state, however
for any finite value $s$, $\tilde{\rho}_{s}$ can be interpreted only as a
formal density matrix.

Once the density matrix $\tilde{\rho}_{s}$ is obtained, relevant quantities of
interest can be calculated. For instance, the $s$-dependent correlation
functions read

$C_{s}^{AB}=\text{Tr}\left(AB\tilde{\rho}_{s}\right)=s\text{Tr}\left(AB\mathcal{R}_{s}\rho_{0}\right),$
(8)

where $A$ and $B$ are generic self-adjoint operators.

###  II.2 Two-process stochastic unravelling

As mentioned earlier, the magical formula (5) can be interpreted as an
integration scheme for any $r>0$. In this respect, the superoperator

$T_{r}(s)=r\left[\mathcal{R}_{s+r}^{\text{free}}\left(1+\frac{\mathcal{L}^{\text{int}}}{r}\right)\right],$
(9)

should be regarded as the fundamental propagator in the Laplace domain. As
anticipated earlier, a truncation order $M_{\text{trunc}}$ is set in the
numerical implementation, which leads to the following approximation of the
evolution superoperator

$\mathcal{R}_{s}\approx\sum_{m=0}^{M_{\textrm{trunc}}}\left[T_{r}(s)\right]^{m}\mathcal{R}_{s+r}^{\text{free}}.$
(10)

From this approximation, we develop a stochastic process to unravel the von
Neumann equation. We take inspiration from the method originally introduced as
triplet unravelling in [47]. Our method tracks the evolution of triplets of
the form $(c_{m},|\phi_{m}\rangle,|\psi_{m}\rangle)$ in the many-body Hilbert
space where the subscript ${m}$ refers to the $m$-th order iterative term in
the evolution described by the propagator (9). The evolution of the piece-
wise-deterministic stochastic processes $|\phi_{m}\rangle$, $|\psi_{m}\rangle$
alternates between an exact continuous evolution, governed by
$H^{\text{free}}$, and stochastic state jumps accounting for $H^{\text{int}}$.
We construct these processes as to reproduce the application of the magical
formula Eq. (5) for infinite order of iteration $M_{\text{trunc}}\to\infty$
via

$\tilde{\rho}_{s}=s\sum_{m=0}^{\infty}\mathbb{E}[c_{m}|\phi_{m}\rangle\langle\psi_{m}|],$
(11)

where $\mathbb{E}[\cdot]$ represents the expectation value of the stochastic
processes. This construction ensures that the solution of the von Neumann
equation in the Laplace domain is recovered in the limit of a large order of
iteration. Because eigenstates of $H^{\text{free}}$ are assumed to be known
exactly, we use these states as the computational basis allowing us to compute
the free evolution exactly. In order to minimize the statistical correlations
between the processes $|\phi_{m}\rangle$ and $|\psi_{m}\rangle$, we use two-
sided jump processes and model the effect of $H^{\text{int}}$ as follows

$|\phi_{m}\rangle\langle\psi_{m}|\mapsto|\phi_{m}\rangle\langle\psi_{m}|-\frac{i}{r}(H^{\text{int}}|\phi_{m}\rangle\langle\psi_{m}|-|\phi_{m}\rangle\langle\psi_{m}|H^{\text{int}}),$
(12)

which coincides with a stochastic application of the superoperator
$\left(1+\frac{\mathcal{L}^{\text{int}}}{r}\right)$.

##  III Dynamic algorithm

The present algorithm is the dynamical adaptation of the FPQMC method
introduced in [47], of which it preserves the overall structure. The latter
consists of a main loop, whose application is iterated over an ensemble of
walkers until a solution emerges. Two distinct steps compose the main loop,
the interaction and the free evolution that are implemented by the
superoperators $\left(1+\frac{\mathcal{L}^{\text{int}}}{r}\right)$ and
$\mathcal{R}_{s+r}^{\text{free}}$, respectively.

###  III.1 The main loop

The algorithm considers the evolution of an ensemble of Monte Carlo walkers,
called triplets, defined as $\\{(w_{n},|i_{n}\rangle,|j_{n}\rangle)\\}_{n}$,
where $w_{n}$ is the (complex) weight factor and $|i_{n}\rangle,|j_{n}\rangle$
are the free Hamiltonian eigenstates. Here the subscript $n$ refers to the
index in the triplet ensemble. In what follows, we will use the lighter
notation $|i\rangle\equiv i$. As anticipated earlier, during the evolution,
the triplets repeatedly experience stochastic spawning events, realized by
state jumps associated with $H^{\textrm{int}}$, interspersed with continuous
weight updates coming from the exact, free evolution $H^{\textrm{free}}$. The
rate of the spawning events is $r$ and it is fixed prior to the simulation.

In order to control the evolution of the ensemble size, we introduce the
following parameters: i) the walker’s unit of weight $w_{u}$, which also fixes
the initial walker’s population of the statistical ensemble, ii) the dead
weight $u_{\text{dw}}$ which is used throughout the simulation to identify
statistically unimportant triplets, and iii) the spring constant $\kappa$ of
the importance sampling procedure. The dead weight and spring constant are
approximation and bias parameters, respectively. As such, they are not rooted
in the theoretical framework introduced in Sec. II, but they are extra
computational features that are necessary for numerical stabilization and
convergence. Both deadweight approximation and importance sampling are
discussed in details in subsection III.2. However, as they are crucial for the
execution of the main loop described in this subsection, we briefly recall why
they are needed. The deadweight approximation aims to reduce the effect of the
dynamical sign problem which is known to cause numerical instabilities
connected to fast oscillating exponentials, especially at longer times.
Importance sampling reduces the number of triplets needed to obtain a
convergent simulation by preventing unwanted extra exploration of the Hilbert
space. During the simulation, the walker’s population is controlled by both
$\kappa$ and $u_{\text{dw}}$. In general, if the spring constant is too large,
the triplet population decreases fast. The same being true for the deadweight
parameter, both parameters should be optimized. In principle, one should first
find the smallest $u_{\text{dw}}$ that overcomes the dynamical sign problem
and, subsequently, decrease the spring constant until a target value of the
triplet’s population is achieved.

We now discuss in detail the main loop of our algorithm, depicted in Fig. 1.
The compression and decompression steps are thoroughly described in Appendix
A.

  1. 1.

Spawnings

    1. (a)

Preparation to interaction.

      1. i.

Deactivation of unimportant triplets $(w_{n},i_{n},j_{n})$. If
$|w_{n}|<u_{\text{dw}}$, the triplet becomes inactive and survives with
probability $|w_{n}|/u_{\text{dw}}$. The weight of surviving triplets is
updated to $u_{\text{dw}}w_{n}/|w_{n}|$. Surviving inactive triplets do not
experience any collision and are therefore evolved only freely.

      2. ii.

Pre-spawning decompression. The active triplets are split into
$N_{c}=\lfloor|w_{n}|/w_{u}\rfloor$ child triplets with weight $w_{n}/N_{c}$.

    2. (b)

Spawning. For each child triplet, a side is randomly chosen (quantum
mechanically, either $i_{n}$ or $j_{n}$) and a new state $k_{n}$ is selected
with equal probability among the $n_{t}$ possible spawning transitions. For
instance, if the ket is chosen, the following triplet is added to the ensemble

$(i\frac{H^{\text{int}}_{i_{n}k_{n}}}{r}\frac{w_{n}}{N_{c}}2n_{t}T_{b},k_{n},j_{n}),$
(13)

where $T_{b}=e^{\kappa(n_{ij}^{2}-n_{kj}^{2})/2}$ is the transition bias, with
$n_{ij}$ the minimum number of application of $H^{\text{int}}$ to transition
from state $i$ to $j$ (see next subsection for details). If $j_{n}$ is
selected, an equivalent triplet is spawned on the bra side.

  2. 2.

Free evolution

    1. (a)

Full compression. A class containing all the triplets associated to a specific
pair $(i_{n},j_{n})$ is replaced by a single triplet whose weight is the sum
of all the weights in the class.

    2. (b)

Weight update. For each class $(w_{n},i_{n},j_{n})$, a weight update is
performed, accounting for the free evolution

$w_{n}\mapsto\frac{r}{s+r+i(H_{i_{n}i_{n}}^{\text{free}}-H_{j_{n}j_{n}}^{\text{free}})}w_{n}.$
(14)

In a typical simulation, correlation functions of the form of Eq. (8) are
calculated via a procedure involving an initialization, an iteration on the
main loop, and a computation of correlation functions. After state
initialization, which matches the triplet ensemble to the desired initial
state $\rho_{0}$, the main loop is executed $M_{\text{trunc}}$ times, with the
cutoff $M_{\text{trunc}}$ being set prior to the simulation. At the end of
each loop $m\leq M_{\text{trunc}}$, the quantity

$\text{Tr}\left(AB\tilde{\rho}_{s}^{(m)}\right)$ (15)

is calculated and stored, with $\tilde{\rho}_{s}^{(m)}$ being the matrix
representation of the output of the $m-1$ loop. Finally, the correlation
function can be computed by summing over all the contribution from all loops,
i.e.

$C_{s}^{AB}\approx\sum_{m=1}^{M_{\text{trunc}}}\text{Tr}\left(AB\tilde{\rho}_{s}^{(m)}\right).$
(16)

Due to the importance sampling procedure, the ensemble is biased by the spring
constant $\kappa$, and triplets’ weight have to be compensated to obtain
unbiased physical results. For a general operator $X$, physical expectation
value is computed using the physical weights
$e^{\frac{\kappa}{2}n_{ij}^{2}}w_{n}$ as follows

$\text{Tr}\left(X\tilde{\rho}\right)\approx\sum_{n}w_{n}e^{\frac{\kappa}{2}n_{ij}^{2}}X_{i_{n}j_{n}},$
(17)

where $\tilde{\rho}$ is represented by the ensemble
$\\{w_{n},i_{n},j_{n}\\}_{n}$ and $X_{i_{n}j_{n}}=\langle
i_{n}|X|j_{n}\rangle$.

Figure 1: Flowchart of the algorithm describing the application of the
superoperator $\mathcal{L}^{\text{int}}/r$. The green diamond represents the
deadweight approximation.

###  III.2 Numerical stabilization and convergence

The size of a typical many-body Hilbert space imposes a strong demand on the
computational resources and limits simulations to small systems. Furthermore,
the dynamical sign problem represents a non-negligible technical issue in
terms of numerical stability. In order to tackle these two problems, we
complement the two-process unravelling loop described above with two extra
procedures, importance sampling and the deadweight approximation.

The importance sampling procedure relies on two main concepts, originally
introduced and explained in details in Sec. III.C of [47]. The first is the
introduction of a dynamic norm $n_{ij}$, a distance between states in the
Hilbert space, defined as the minimum number of applications of
$H^{\text{int}}$ to transition from state $i$ to $j$. The second concept
associates to each triplet two weights, a physical weight $c$ and an ensemble
weight $w$. The physical weight $c$ is used to compute expectation values,
e.g. $\text{Tr}(A\tilde{\rho}_{s})$, whereas the ensemble weight $w$ reflects
the number of jumps performed by the triplet. The latter weight allows to rate
the statistical relevance of the triplet within its ensemble without altering
the physical averages. These two are related via the equation $c=bw$, by a
norm-dependent bias $b\equiv b(n_{ij})$. A standard choice for $b$ is
$b=\exp\left(\frac{\kappa}{2}n_{ij}^{2}\right)$, which is a harmonic
interaction with spring constant $\kappa$ between the states $i$ and $j$ that
aims at reducing the number of triplets with a large dynamic norm.

The dynamical sign problem manifests itself as numerical instabilities
occurring at medium times. In order to avoid divergences, we need to prevent
statistically irrelevant triplet from over-contributing to the dynamics. The
basic idea of the deadweight approximation is to forbid all the triplets with
an ensemble weight $w$ below a threshold $u_{\text{dw}}$ from experiencing
jumps while still allowing them to undergo free evolution. The triplets
undergoing the jumps are called active triplets whereas the ones below the
threshold $u_{\text{dw}}$ as referred to as inactive. Intuitively, one can
imagine these inactive triplets as an effective environment whose free
dynamics is needed in order to guarantee an ergodic exploration of the Hilbert
space by the remaining active triplets.

An extra computational feature which greatly improves the efficiency of the
algorithm concerns the parametric nature of $s$ in the magical formula (5).
Unlike time evolution, where typically each iteration at a given physical time
$t^{\prime}=t+\Delta t$ produces an updated density matrix
$\rho_{t^{\prime}}$, the evolution dictated by the magical formula yields
results for a density matrix $\tilde{\rho}_{s}$ that represents an integral
over the whole time range, as reflected in the use of the Laplace transform
(3). At a first glance, one might be skeptical about the efficiency of our
method as one would need to produce an independent simulation for each $s$
value to reconstruct $\tilde{\rho}_{s}$. However, by noting that the $s$
variable only appears in the exact free part of the total evolution, one can
produce results for the whole $s$ range in a single simulation at a minimal
additional computational cost. For example, two triplets
$(c^{(1)}_{m},|\phi^{(1)}_{m}\rangle,|\psi^{(1)}_{m}\rangle)$ and
$(c^{(2)}_{m},|\phi^{(2)}_{m}\rangle,|\psi^{(2)}_{m}\rangle)$ having the very
same trajectory for two different values $s_{1}\neq s_{2}$ differ from one
another only in the application of the free evolution operator
$\mathcal{R}_{s+r}^{\textrm{free}}$. The idea is then to evolve a single copy
of these triplets and reweight each application of the free evolution operator
when computing observables simultaneously for $s_{1}$ and $s_{2}$. In
practice, we select the range of $s$ values we are interested in before the
simulation, and at each iterative order $m$, we reweight the triplets for the
whole range of $s$ values during the free evolution part.

##  IV Results

In this section, we test our algorithm on two 1d models in quantum magnetism,
namely the Heisenberg XXZ and the Ising model, and look at the dynamics of
excitations following a quench. We stress that, being a QMC method, it can be
applied to 2d and 3d systems as well, with the only extra difficulty being the
increasing size of the Hilbert space. Although the use of the Laplace domain
indicates that the long-term quantities can be naturally computed, we will
show that certain oscillation frequencies can also be extracted from functions
computed on the Laplace domain.

###  IV.1 Heisenberg XXZ model: analysis of the method and quenched dynamics

In order to assess the efficiency of our method, we first present a
quantitative analysis of our numerical method. In particular we focus on the
walker’s population control via the parameters $w_{u}$ and $u_{\text{dw}}$,
the dependence of the end result on the importance sampling parameter $\kappa$
and the trace preserving properties of the algorithm.

Like in any Monte Carlo method, the walker population is the main quantity
responsible for the size of the statistical error bars. Within this method,
the population for a given loop refers to the number of spawning attempts
during that loop, and due to the high number of spawnings, it naturally
increases at an exponential rate. It is therefore vital to slow down any
possible overgrowth and limit the exploration of the Hilbert space to
statistically important triplets only. In practice, whenever the target value
of the triplets’ population is reached, the deadweight approximation is
enabled by setting a non-vanishing, fine-tuned threshold $u_{\text{dw}}$. This
drastically reduces exponential increase in population and further keeps it at
a nearly constant value.

To be more concrete we consider the well-known Heisenberg XXZ model with $L$
spins and open boundary conditions. Its Hamiltonian reads

$H=J_{xy}\sum_{i=1}^{L-1}(\sigma_{i}^{x}\sigma_{i+1}^{x}+\sigma_{i}^{y}\sigma_{i+1}^{y})+J_{z}\sum_{i=1}^{L-1}\sigma_{i}^{z}\sigma_{i+1}^{z},$
(18)

where $\sigma^{x},\sigma^{z}$ are the standard Pauli matrices and the sum runs
over the spins in the chain. Following the recipe introduced earlier, the
Hamiltonian (18) is split into a free and an interacting part,
$H=H^{\text{free}}+H^{\text{int}}$, as follows

$H^{\text{free}}=J_{z}\sum_{i}\sigma_{i}^{z}\sigma_{i+1}^{z},\;H^{\text{int}}=J_{xy}\sum_{i}\sigma_{i}^{+}\sigma^{-}_{i+1}+\sigma_{i}^{-}\sigma_{i+1}^{+},$
(19)

where $\sigma^{\pm}=\sigma^{x}\pm i\sigma^{y}$. This particular splitting of
the Hamiltonian sets the eigenbasis of $H^{\text{free}}$ as the computational
basis, that is

$|e_{s_{1},\dots,s_{L}}\rangle=|s_{1}\rangle\otimes|s_{2}\rangle\otimes\dots\otimes|s_{L}\rangle,$
(20)

where $s_{a}=s^{z}_{a}=\pm 1$. The single-particle operators $\sigma^{z}$ and
$\sigma^{\pm}$ act of the basis states according to the standard algebra of
Pauli matrices

$\sigma^{z}|\pm 1\rangle=\pm|\pm 1\rangle,\qquad\sigma^{\pm}|\mp 1\rangle=|\pm
1\rangle,\qquad\sigma^{\pm}|\pm 1\rangle=0.$ (21)

The initial state is set to

$\rho_{\text{ini}}=|\text{ini}\rangle\langle\text{ini}|,$ (22)

with

$|\text{ini}\rangle=\bigotimes_{i=1}^{L/2}|+1\rangle_{i}\bigotimes_{i=L/2+1}^{L}|-1\rangle_{i},$
(23)

in which strongly localized spin excitations serve as an initial source of
energy. This initial state will be useful later to understand the propagation
of spin excitations across the chain.

Figure 2 shows three examples of the evolution of the population and the
accuracy of the numerical integration for these populations in the $L=10$ XXZ
model, for $J_{z}=0.9$, $J_{xy}=1$, $s=5\cdot 10^{-2}$, $r=30$ and $\kappa=2$.
In all the examples, the population control is performed by enabling the
deadweight approximation at $m/r=8$, represented by the vertical line in the
upper panel. The value of the threshold $u_{\text{dw}}$ then controls the
second phase of the exponential increase of the population. In passing, we
note that $w_{u}$ can be responsible for some secondary effects later during
the evolution, e.g., the decreasing distance between the dotted red and dashed
orange lines in the first stage in Fig. 2.

Figure 2: Upper panel: Evolution of the triplet population for various values
of the deadweight threshold $u_{\text{dw}}$ and triplet’s unit weight $w_{u}$
for the $L=10$ XXZ model with $J_{xy}=1$, $J_{z}=0.9$, $s=1/20$, $r=30$ and
$\kappa=2$. The error bars were calculated over 30 independent simulations.
Lower panel: evolution of the quantity
$\text{Tr}(\sigma_{L/2}^{z}\tilde{\rho}^{(m)}_{s})$ defined in Eq. (15) for
the same parameters as above.

The introduction of the importance sampling procedure was motivated by the
need to reduce the number of triplets in the ensemble to only the
statistically important ones. Thus, the main effect of the $\kappa$ parameter
is a reduction of the number triplets with a large dynamic norm. A larger
$\kappa$ leads to a smaller exponential increase of the population. However,
if the spring constant is chosen too large, it may result in an underestimate
of the final values of the correlation functions (8). This gap originates from
breaking the ergodicity of Hilbert space exploration, with the triplets being
confined to a small region of the Hilbert space around their initial state.
Figure 3 illustrates the underestimation of the quantity
$C_{s}^{\sigma_{L/2}^{z}}$ in the $L=8$ XXZ model with $J_{z}=1.5$, $J_{xy}=1$
for different values of the spring constant. We note that the importance
sampling induces a systematic error for values $\kappa\leq 4$, a critical
value below which the population is large enough to avoid ergodicity breaking.
On a similar note, we remark that the deadweight approximation can also affect
the final results and a $u_{\text{dw}}$ analysis should be carried out,
similarly to $\kappa$. In order to obtain the optimal threshold, it is more
efficient to start from a larger value of $u_{\text{dw}}$ and further decrease
it until the final result becomes constant within the acceptable statistical
error.

Figure 3: Dependence of the spring constant $\kappa$ on the final result
$C_{s}^{\sigma^{z}_{L/2}}(\kappa)$ defined in Eq. (8) for a $L=8$ XXZ model
with $J_{z}=1.5$, $J_{xy}=1$ and $s=1/10$. The errors were calculated over 100
independent simulations with $u_{\text{dw}}=1.125w_{u}$ and $r=30$.

The trace conserving character of the evolution operator in the Laplace domain
$\mathcal{R}_{s}$ can be checked using the following condition [51]

$\lim_{s\to 0}\text{Tr}\left(\tilde{\rho}_{s}\right)=\text{Tr}(\rho_{0}).$
(24)

If the algorithm violates the unitary character of the time-evolution,
instabilities will occur and grow exponentially with the simulation length.
Figure 4 shows trace conservation according to Eq. (24) with $\rho_{0}$ chosen
as the ground state of the XXZ Hamiltonian with $L=20$. The exponential growth
of the error bars for small $s$ signals the occurring of the dynamical sign
problem. The severity the dynamical sign problem generally depends on the
details of the system as well as the observables being computed.

Figure 4: Trace conservation of the magical formula Eq. (5) for the $L=20$ XXZ
model with $J_{z}=J_{xy}=1$. The population was about $2\cdot 10^{4}$ triplets
with $r=100$, $\kappa=4$, $u_{\text{dw}}=1$ at $m/r=4$ and the error bars were
calculated over 30 independent simulations.

We now focus on the dynamical behaviour of spin excitations following an
instantaneous quench. For $J_{z}<J_{xy}$ these travel ballistically across the
chain regardless of its length, whilst, on the contrary, for $J_{z}>J_{xy}$
propagation is inhibited by the strong $\sigma_{z}$ coupling [24]. The
propagation of these spin excitations can be quantified by analyzing the
evolution of the magnetization profile. Starting from $\rho_{\text{ini}}$,
this profile is calculated as

$C_{s}^{\sigma_{i}^{z}}\equiv\;\text{Tr}(\sigma_{i}^{z}\tilde{\rho}_{s}).$
(25)

Figure 5 illustrates the evolution of the absolute value of the magnetization
profile (25) for $J_{xy}=1$ in both the ballistic and strongly interacting
regimes at $J_{z}=0.6$ (top) and $J_{z}=1.5$ (bottom), respectively. In the
upper panel, a power-law decay can be observed. This is consistent with
ballistic propagation where the spin excitations travel within a light-cone
type of spatial and temporal region. On the contrary, in the lower panel at
$J_{z}=1.5$ magnetization propagation is practically suppressed as a result of
the localizing longitudinal interactions.

Figure 5: Evolution in the Laplace domain of the magnetization profile
$|C_{s}^{\sigma^{z}_{i}}|$ for the $L=40$ XXZ chain for $J_{xy}=1$. Ballistic
behaviour is clearly visible in the upper panel, whereas transport is almost
completely suppressed in the lower panel. The values plotted were averaged
over 100 independent simulations with a population of about $10^{5}$ triplets
and $r=100$, $\kappa=4$, $u_{\text{dw}}=1$ at $m/r=12$.

To measure this sharp transition between ballistic and suppressed propagation,
we consider the Loschmidt echo
$L(t)=|\langle\psi_{0}|e^{-iHt}|\psi_{0}\rangle|^{2}$, a measure of the
disturbance induced on a quantum system by an external perturbation. In the
Laplace domain, this reads,

$L(s)=s\;\text{Tr}(\rho_{0}\mathcal{R}_{s}\rho_{0})\equiv\text{Tr}(\rho_{0}\tilde{\rho}_{s}),$
(26)

where $\rho_{0}=|\psi_{0}\rangle\langle\psi_{0}|$ is the initial state. Figure
6 illustrates the evolution of $L(s)$ for $J_{xy}=1$ at various $J_{z}$ for a
spin chain of length $L=40$ from the initial state Eq. (22). Increasing
$J_{z}$ in the Hamiltonian (18) suppresses the propagation of magnetization.
The quantum state of the spin chain essentially freezes up as a consequence of
strong localization, and this results in a weak decay to nearly constant
Loschmidt echo. This observation is supported by previous results obtained via
exact diagonalization [52].

Figure 6: Evolution of the Loschmidt echo in the Laplace domain
$L(s)=s\;\text{Tr}(\rho_{\text{ini}}\mathcal{R}_{s}\rho_{\text{ini}})$ for the
$L=40$ XXZ chain with $J_{yx}=1$ and various values of $J_{z}$. The errors
were estimated over 30 independent simulations, with $r=30$, $\kappa=1$ and
$u_{\text{dw}}=1.2$ at $m/r=12$. The population at the end of each simulation
was about $2\cdot 10^{6}$ triplets. The continuous curve represents the exact
solution for $J_{z}=0$.

###  IV.2 Confinement in the Quantum Ising model

For our second case study, we consider a ferromagnetic quantum Ising chain of
length $L$ coupled to the transverse and longitudinal magnetic fields $h_{x}$
and $h_{z}$, respectively. The Hamiltonian reads

$H=-J\sum_{i=1}^{L-1}\sigma_{i}^{z}\sigma_{i+1}^{z}-h_{z}\sum_{i=1}^{L}\sigma_{i}^{z}-h_{x}\sum_{i=1}^{L}\sigma_{i}^{x},$
(27)

where $\sigma^{x},\sigma^{z}$ are the standard Pauli matrices and the sum runs
over the spins in the chain. In the absence of a longitudinal field, the model
is exactly solvable and has a phase transition at $h_{x}=J$. In the ordered
phase $h_{x}<J$, the ground state is degenerate due to a spin-flip symmetry
breaking and corresponds to domain walls of various lengths. For instance, for
$h_{x}=h_{z}=0$, the domain walls’ lengths are maximized, and the degenerate
ground states are
$|\Psi_{\textup{u}}\rangle=\bigotimes_{i=1}^{L}|\uparrow\rangle_{i}$ and
$|\Psi_{\textup{d}}\rangle=\bigotimes_{i=1}^{L}|\downarrow\rangle_{i}$. A non-
zero value of the longitudinal field $h_{z}$ creates an energy gap between
these states by increasing the energy of the spin domains along the field
$h_{z}$. Hence, the latter field acts as an attracting potential between the
two walls delimiting a domain.

Using our algorithm, we now show how spin dynamics can be suppressed in an
Ising chain in the above conditions. Similarly to the XXZ model studied in the
previous subsection, we initialize the evolution of the spin chain from the
domain wall state
$|\text{ini}\rangle=|\uparrow\dots\uparrow\downarrow\dots\downarrow\rangle$.
We quench the longitudinal field from $h_{z}=0$ to a non-zero value and study
the resulting dynamics. The initial state $|\text{ini}\rangle$ is
characterized by a spin kink exactly in the middle of the chain. Hence, a non-
zero longitudinal field will induce an energy imbalance between the two spin
domains. Figure 7 shows the evolution in the Laplace domain of the energy
density profile $\text{Tr}(H_{i}\tilde{\rho}_{s})$, with

$H_{i}=-J\sigma_{i}^{z}\sigma_{i+1}^{z}-\frac{h_{x}}{2}(\sigma_{i}^{x}+\sigma_{i+1}^{x})-\frac{h_{z}}{2}(\sigma_{i}^{z}+\sigma_{i+1}^{z}),$
(28)

starting from the initial state Eq. (22) for an $L=40$ chain with $J=1$,
$h_{x}=0.2$, and $h_{z}=1.2$. Neither energy exchange nor spin excitation
propagation occurs between the two halves. This can be explained by
considering the mid kink as a quasi-particle, whose motion is triggered by a
non-zero longitudinal field $h_{z}$. The kinetic energy gain, which is of
order $\sim h_{z}$, allows the kink to move within the potential, but due to
energy conservation it has to periodically bounce back, leading to a confined
dynamics (see [50] for details). These oscillations are centered around the
central bond connecting $L/2$ and $L/2+1$ and they can be measured by
considering the evolution of $\langle\sigma^{z}_{L/2}\rangle$. This is
illustrated in the upper panel of Fig. 8 in the Laplace domain for a $L=20$
chain with $J=1$, $h_{x}=0.2$ and several confining potentials. Because the
result is shown in the Laplace domain, the oscillations take the form of
Lorezian curves. The mean frequency of the confining oscillations can be
extracted as the inflection point of these curves. This is shown on the lower
panel of Fig. 8, which illustrates the logarithmic derivative of the signal in
the upper panel. The derivatives were calculated using a cubic spline fitting,
and the oscillations frequency is represented by the location of the peak.
Furthermore, the oscillation amplitudes can be estimated from the distance
between the values $C_{s=0}^{\sigma_{L/2}^{z}}$ and
$C_{s=\infty}^{\sigma_{L/2}^{z}}$. These amplitudes and frequencies scale as
$\sim h_{x}/h_{z}$ and $\sim h_{z}$, respectively, which is in agreement with
the quasi-particle interpretation of the kink dynamics and with the findings
reported in [50].

Figure 7: Evolution in the Laplace domain of the energy density profile
$C_{s}^{H_{i}}$ for a $L=40$ quantum Ising chain with $J=1$, $h_{x}=0.2$, and
$h_{z}=1.2$ starting from a state with a single kink in the middle. The
profile was computed over 30 independent simulations for $r=30$, $\kappa=3$
and $u_{\text{dw}}=1$ at $m/r=6$. Figure 8: Upper panel: evolution in the
Laplace domain of the mean magnetization of the spin on the left of the
central opposite pair $C_{s}^{\sigma_{L/2}^{z}}$ for the $L=20$ quantum Ising
chain with $J=1$, $h_{x}=0.2$ and the parameters $r=30$, $\kappa=2$, and
$w_{u}=4\cdot 10^{-6}$. The dead weight approximation was enable at $m/r=6$
with $u_{\text{dw}}=1.6,~{}1.5,~{}1.8$ for $h_{z}=1.2,~{}0.6,~{}0.3$,
respectively. The choice of specific thresholds $u_{\text{dw}}$ originates
from the optmization procedure. Lower panel: logarithmic derivative of the
magnetization $C_{s}^{\sigma_{L/2}^{z}}$. The derivative was calculated using
a cubic spline fitting.

This is further supported by Fig. 9, which displays the time evolution of the
mean value of the magnetization $\langle\sigma_{L/2}^{z}\rangle_{t}$ obtained
via an inverse Laplace transform of the data in the upper panel in Fig. 8,
using the Zakian method. The inverting procedure relies on a rational
polynomial fit of the original data in order to extrapolate the whole Laplace-
transformed function. It starts by fitting the signal in Fig. 8 by a function
of the type $P_{2}/P^{\prime}_{2}$ where $P_{2}$, $P^{\prime}_{2}$ are
polynomials of degree two. The optimal parameters of the first fit are then
used as an initial guess for a fit of higher order polynomials, say
$P_{3}/P^{\prime}_{3}$. This iterative scheme continues until the parameters
start diverging. As illustrated in Fig. 9, the frequency of the oscillations
depends linearly on $h_{z}$ while the amplitude is inversely proportional. The
numerical inverse Laplace method induces numerical amplitude damping, due to
the use of finite of the $s$ range. The sharp increase in the magnetization
for $h_{z}=0.6$ at about $t\approx 8$ comes from the dynamical sign problem as
well as an interplay between the numerical inverse method and the
extrapolation over the whole $s>0$ domain, suggesting that the results in the
real time domain can be trusted up to $t\approx 5$. It is however worthy of
noting that the Laplace signal does not contain any signs of divergences and
can hence be trusted at least up to $s\approx 10$.

Figure 9: Time evolution of the magnetization
$\langle\sigma_{L/2}^{z}\rangle_{t}$ calculated from an inverse Laplace
transform from the data of Fig. 8.

##  V Conclusions and open perspectives

In this work, we have introduced an alternate Monte Carlo method to study the
reversible dynamics of many-body quantum systems. The validity and the
accuracy of the method was benchmarked against excitation propagation in the
Heisenberg XXZ model and dynamical confinement in the quantum Ising chain. In
both cases, we were able to resolve long-time dynamical properties.

Similarly to FPQMC, this method is based on a piece-wise stochastic
deterministic two-process unravelling to solve the von Neumann equation. The
introduction of an importance sampling procedure allows us to limit the
exploration of the Hilbert space to statistically important states.
Furthermore, large times can be reached with the deadweight approximation by
preventing statistically unimportant triplets from spawning. Due to the
unravelling and the form of the evolution operator in the Laplace domain,
these triplets can be kept in the simulation by only participating to the free
evolution, which greatly reduces the statistical errors on the final result
and delays the appearance of the divergences due to the dynamical sign
problem. The trace of the density matrix is conserved at all times during the
simulation. The maximum simulation time accessible is highly dependent on the
observable, as it has been shown between the Loschmidt echo and the
oscillations in the quantum Ising model. Unfortunately, we doubt that larger
time scales can be reached.

Regarding future directions, we believe that the introduction of a small
dissipation parameter can reduce the severity of the sign problem, an idea
that is already under current investigation. Thus, we are confident that our
method can be extended to dissipative many-body quantum systems.

##  Appendix A Compression and decompression steps

Prior to the execution of the loop the ensemble is modified as to improve the
statistics without influencing directly the averages. This modification is
carried out through compressions or decompressions. In a compression, classes
of triplets are formed by grouping together all the triplets associated to a
fixed pair of states, for instance, $(i,j)$. These are then replaced by a
single triplet whose weight is equal to the sum of the weights of all the
members of the class. Decompression is applied on a compressed ensemble. A
single class of triplets $(i,j)$ is split into triplets with weight $w_{u}>0$
in absolute value, $(\text{sgn}(w_{n})w_{u},i_{n},j_{n})$, and a single rest
triplet $(w_{r},i_{n},j_{n})$, with
$w_{r}=\text{sgn}(w_{n})(|w_{n}/w_{u}|-\lfloor|w_{n}/w_{u}|\rfloor)$
($\lfloor\cdot\rfloor$ is the floor function). The rest triplet is then
removed from the simulation with probability $1-|w_{r}|$; otherwise its weight
is updated to $\text{sgn}(w_{n})w_{u}$. That way, the total statistical weight
is conserved on average.

## References

  * Eisert  _et al._ [2015] J. Eisert, M. Friesdorf, and C. Gogolin, Quantum many-body systems out of equilibrium, Nature Physics 11, 124 (2015).
  * Lake  _et al._ [2010] B. Lake, A. M. Tsvelik, S. Notbohm, D. Alan Tennant, T. G. Perring, M. Reehuis, C. Sekar, G. Krabbes, and B. Büchner, Confinement of fractional quantum number particles in a condensed-matter system, Nature Physics 6, 50 (2010).
  * Hild  _et al._ [2014] S. Hild, T. Fukuhara, P. Schauß, J. Zeiher, M. Knap, E. Demler, I. Bloch, and C. Gross, Far-from-equilibrium spin transport in Heisenberg quantum magnets, Phys. Rev. Lett. 113, 147205 (2014).
  * Ronzheimer  _et al._ [2013] J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, S. Langer, I. P. McCulloch, F. Heidrich-Meisner, I. Bloch, and U. Schneider, Expansion dynamics of interacting bosons in homogeneous lattices in one and two dimensions, Phys. Rev. Lett. 110, 205301 (2013).
  * Bloch  _et al._ [2012] I. Bloch, J. Dalibard, and S. Nascimbène, Quantum simulations with ultracold quantum gases, Nature Physics 8, 267 (2012).
  * Gross and Bloch [2017] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357, 995 (2017).
  * Schäfer  _et al._ [2020] F. Schäfer, T. Fukuhara, S. Sugawa, Y. Takasu, and Y. Takahashi, Tools for quantum simulation with ultracold atoms in optical lattices, Nature Reviews Physics 2, 411 (2020).
  * Blatt and Roos [2012] R. Blatt and C. F. Roos, Quantum simulations with trapped ions, Nature Physics 8, 277 (2012).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021).
  * Romero  _et al._ [2017] G. Romero, E. Solano, and L. Lamata,  _Quantum Simulations with Circuit Quantum Electrodynamics_, edited by D. G. Angelakis (Springer International Publishing, 2017).
  * Petrescu  _et al._ [2018] A. Petrescu, H. E. Türeci, A. V. Ustinov, and I. M. Pop, Fluxon-based quantum simulation in circuit QED, Phys. Rev. B 98, 174505 (2018).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Alet and Laflorencie [2018] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics, Comptes Rendus Physique 19, 498 (2018).
  * Deutsch [2018] J. M. Deutsch, Eigenstate thermalization hypothesis, Reports on Progress in Physics 81, 082001 (2018).
  * Cheneau  _et al._ [2012] M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P. Schauß, T. Fukuhara, C. Gross, I. Bloch, C. Kollath, and S. Kuhr, Light-cone-like spreading of correlations in a quantum many-body system, Nature 481, 484 (2012).
  * Geiger  _et al._ [2014] R. Geiger, T. Langen, I. E. Mazets, and J. Schmiedmayer, Local relaxation and light-cone-like propagation of correlations in a trapped one-dimensional bose gas, New Journal of Physics 16, 053034 (2014).
  * Faupin  _et al._ [2022] J. Faupin, M. Lemm, and I. M. Sigal, Maximal speed for macroscopic particle transport in the Bose-Hubbard model, Phys. Rev. Lett. 128, 150602 (2022).
  * Canovi  _et al._ [2014] E. Canovi, P. Werner, and M. Eckstein, First-order dynamical phase transitions, Phys. Rev. Lett. 113, 265702 (2014).
  * Heyl [2019] M. Heyl, Dynamical quantum phase transitions: A brief survey, EPL (Europhysics Letters) 125, 26001 (2019).
  * Link and Strunz [2020] V. Link and W. T. Strunz, Dynamical phase transitions in dissipative quantum dynamics with quantum optical realization, Phys. Rev. Lett. 125, 143602 (2020).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Bertini  _et al._ [2021] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Žnidarič, Finite-temperature transport in one-dimensional quantum lattice models, Rev. Mod. Phys. 93, 025003 (2021).
  * Gobert  _et al._ [2005] D. Gobert, C. Kollath, U. Schollwöck, and G. Schütz, Real-time dynamics in spin-$\frac{1}{2}$ chains with adaptive time-dependent density matrix renormalization group, Phys. Rev. E 71, 036102 (2005).
  * Sirker  _et al._ [2009] J. Sirker, R. G. Pereira, and I. Affleck, Diffusion and ballistic transport in one-dimensional quantum systems, Phys. Rev. Lett. 103, 216602 (2009).
  * Al-Hassanieh  _et al._ [2006] K. A. Al-Hassanieh, A. E. Feiguin, J. A. Riera, C. A. Büsser, and E. Dagotto, Adaptive time-dependent density-matrix renormalization-group technique for calculating the conductance of strongly correlated nanostructures, Phys. Rev. B 73, 195304 (2006).
  * Aoki  _et al._ [2014] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Nonequilibrium dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014).
  * Lacroix and Ayik [2014] D. Lacroix and S. Ayik, Stochastic quantum dynamics beyond mean field, The European Physical Journal A 50, 95 (2014).
  * Dong  _et al._ [2017] Q. Dong, I. Krivenko, J. Kleinhenz, A. E. Antipov, G. Cohen, and E. Gull, Quantum Monte Carlo solution of the dynamical mean field equations in real time, Phys. Rev. B 96, 155126 (2017).
  * Mizuno  _et al._ [2000] Y. Mizuno, K. Tsutsui, T. Tohyama, and S. Maekawa, Nonlinear optical response and spin-charge separation in one-dimensional mott insulators, Phys. Rev. B 62, R4769 (2000).
  * Amaricci  _et al._ [2012] A. Amaricci, C. Weber, M. Capone, and G. Kotliar, Approach to a stationary state in a driven Hubbard model coupled to a thermostat, Phys. Rev. B 86, 085110 (2012).
  * Eckstein  _et al._ [2009] M. Eckstein, M. Kollar, and P. Werner, Thermalization after an interaction quench in the Hubbard model, Phys. Rev. Lett. 103, 056403 (2009).
  * Eckstein and Kollar [2008] M. Eckstein and M. Kollar, Nonthermal steady states after an interaction quench in the Falicov-Kimball model, Phys. Rev. Lett. 100, 120404 (2008).
  * Eckstein and Kollar [2010] M. Eckstein and M. Kollar, Near-adiabatic parameter changes in correlated systems: influence of the ramp protocol on the excitation energy, New Journal of Physics 12, 055012 (2010).
  * McLachlan and Ball [1964] A. D. McLachlan and M. A. Ball, Time-dependent Hartree-Fock theory for molecules, Rev. Mod. Phys. 36, 844 (1964).
  * Nest  _et al._ [2005] M. Nest, T. Klamroth, and P. Saalfrank, The multiconfiguration time-dependent Hartree–Fock method for quantum chemical calculations, The Journal of Chemical Physics 122, 124102 (2005).
  * Jones [2015] R. O. Jones, Density functional theory: Its origins, rise to prominence, and future, Rev. Mod. Phys. 87, 897 (2015).
  * Mühlbacher and Rabani [2008] L. Mühlbacher and E. Rabani, Real-time path integral approach to nonequilibrium many-body quantum systems, Phys. Rev. Lett. 100, 176403 (2008).
  * Guther  _et al._ [2018] K. Guther, W. Dobrautz, O. Gunnarsson, and A. Alavi, Time propagation and spectroscopy of fermionic systems using a stochastic technique, Phys. Rev. Lett. 121, 056401 (2018).
  * Werner  _et al._ [2009] P. Werner, T. Oka, and A. J. Millis, Diagrammatic Monte Carlo simulation of nonequilibrium systems, Phys. Rev. B 79, 035320 (2009).
  * Church and Rubenstein [2021] M. S. Church and B. M. Rubenstein, Real-time dynamics of strongly correlated fermions using auxiliary field quantum Monte Carlo, The Journal of Chemical Physics 154, 184103 (2021).
  * Goth and Assaad [2012] F. Goth and F. F. Assaad, Time and spatially resolved quench of the fermionic Hubbard model showing restricted equilibration, Phys. Rev. B 85, 085129 (2012).
  * Alexandru  _et al._ [2016] A. Alexandru, G. Başar, P. F. Bedaque, S. Vartak, and N. C. Warrington, Monte Carlo study of real time dynamics on the lattice, Phys. Rev. Lett. 117, 081602 (2016).
  * Schiró and Fabrizio [2009] M. Schiró and M. Fabrizio, Real-time diagrammatic Monte Carlo for nonequilibrium quantum transport, Phys. Rev. B 79, 153302 (2009).
  * Gull  _et al._ [2011] E. Gull, D. R. Reichman, and A. J. Millis, Numerically exact long-time behavior of nonequilibrium quantum impurity models, Phys. Rev. B 84, 085134 (2011).
  * Cohen  _et al._ [2015] G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, Taming the dynamical sign problem in real-time evolution of quantum many-body problems, Phys. Rev. Lett. 115, 266802 (2015).
  * Chessex  _et al._ [2022] R. Chessex, M. Borrelli, and H. C. Öttinger, Fixed-point quantum monte carlo method: A combination of density-matrix quantum monte carlo method and stochastic unravellings, Phys. Rev. A 105, 062803 (2022).
  * Kormos  _et al._ [2017] M. Kormos, M. Collura, G. Takács, and P. Calabrese, Real-time confinement following a quantum quench to a non-integrable model, Nature Physics 13, 246 (2017).
  * Verdel  _et al._ [2020] R. Verdel, F. Liu, S. Whitsitt, A. V. Gorshkov, and M. Heyl, Real-time dynamics of string breaking in quantum spin chains, Phys. Rev. B 102, 014308 (2020).
  * Mazza  _et al._ [2019] P. P. Mazza, G. Perfetto, A. Lerose, M. Collura, and A. Gambassi, Suppression of transport in nondisordered quantum spin chains due to confined excitations, Phys. Rev. B 99, 180302 (2019).
  * Öttinger [2017] H. C. Öttinger,  _A Philosophical Approach to Quantum Field Theory_ (Cambridge University Press, Cambridge UK, 2017).
  * Viti  _et al._ [2016] J. Viti, J.-M. Stéphan, J. Dubail, and M. Haque, Inhomogeneous quenches in a free fermionic chain: Exact results, EPL (Europhysics Letters) 115, 40011 (2016).


capbtabboxtable[][] ruled listingtblst listingListing

#  Improving Health Information Access in the World’s Largest Maternal Mobile
Health Program via Bandit Algorithms

Arshika Lalan1, Shresth Verma1∗, Paula Rodriguez Diaz2, Panayiotis Danassis2,
Amrita Mahale3, Kumar Madhu Sudan3, Aparna Hegde3, Milind Tambe1, Aparna
Taneja1 Joint first author

###### Abstract

Harnessing the wide-spread availability of cell phones, many nonprofits have
launched mobile health (mHealth) programs to deliver information via voice or
text to beneficiaries in underserved communities, with maternal and infant
health being a key area of such mHealth programs. Unfortunately, dwindling
listenership is a major challenge, requiring targeted interventions using
limited resources. This paper focuses on Kilkari, the world’s largest mHealth
program for maternal and child care – with over 3 million active subscribers
at a time – launched by India’s Ministry of Health and Family Welfare (MoHFW)
and run by the non-profit ARMMAN. We present a system called CHAHAK that aims
to reduce automated dropouts as well as boost engagement with the program
through the strategic allocation of interventions to beneficiaries. Past work
in a similar domain has focused on a much smaller scale mHealth program and
used markovian restless multiarmed bandits to optimize a single limited
intervention resource. However this paper demonstrates the challenges in
adopting a markovian approach in Kilkari; therefore CHAHAK instead relies on
non-markovian time-series restless bandits, and optimizes multiple
interventions to improve listenership. We use real Kilkari data from the
Odisha state in India to show CHAHAK’s effectiveness in harnessing multiple
interventions to boost listenership, benefiting marginalized communities. When
deployed CHAHAK will assist the largest maternal mHealth program to date.

## Introduction

Ensuring the well-being and health of pregnant women and newborn children is a
significant challenge in the developing world (United Nations Women 2023).
According to WHO (World Health Organization 2023), almost 800 women died
everyday in 2020 from pregnancy and childbirth related issues. Many of these
deaths are preventable with access to critical health information in a timely
manner. Mobile Health (mHealth) programs make use of mobile phones and other
wireless communication technologies to improve access and raise awareness on
healthcare information. ARMMAN (https://armman.org/) is a non-profit
organization in India that conducts several such mHealth programs to improve
the access of pregnant women, mothers of infants, and their families to
critical preventive care information (ARMMAN 2023). One program that ARMMAN
runs is _Kilkari_ (ARMMAN 2019), the _world’s largest maternal mHealth
messaging program, and it is the focus of this paper_. Kilkari was launched in
2016 by India’s Ministry of Health & Family Welfare (MoHFW) as a free mHealth
education service sending women preventive care information though _Automated
Voice Messages (AVMs)_ during pregnancy and child infancy.

A major challenge for mHealth programs (including Kilkari) is that of
dwindling engagement over time (Amagai et al. 2022; Mate et al. 2022), i.e.,
beneficiaries stop listening to the AVMs. In addressing this challenge, SAHELI
(Verma et al. 2023) is a key previous deployed system that has focused on
minimizing such disengagement by optimizing limited interventions (phone calls
by health workers). SAHELI was the first system to deploy bandit algorithms,
specifically markovian restless multi-armed bandits, for allocating limited
interventions in ARMMAN’s smaller scale maternal mHealth program, mMitra
(ARMMAN 2023). mMitra serves approximately 200,000 beneficiaries within
Mumbai’s urban population. In contrast, Kilkari serves around _3.2 million_
beneficiaries at a time across 18 states and union territories in India.
Despite the similarities in both being maternal mHealth programs, mMitra and
Kilkari have several fundamental differences. Most importantly, the lack of
demographic and personal information in Kilkari necessitates the use of novel
_non-markovian_ Time-Series Bandits (TSB) to plan interventions – i.e., past
SOTA techniques on restless multi-armed bandits are not applicable in this
context. Additionally, there is the need to infer time-slot preferences for
sending AVMs to beneficiaries, requiring the use of stochastic (not restless)
bandits.

We propose to deploy a system called Campaign for Health-resource Allocation
using Historical data in AI based system for Kilkari or CHAHAK (Hindi for
“chirping of birds”) to boost the engagement of beneficiaries with the mHealth
program, and, thus, ultimately positively impact behavioural and health
outcomes. To that end, we address two critical optimization challenges. First,
we optimize the appropriate time-slot to send AVMs to each beneficiary.
Practical problems such as limited access to phones due to shared family
phones for many women, working hours, house chore responsibilities etc.,
significantly affect the likelihood of engagement in a given time-slot.
Second, we optimize the limited interventions (e.g., support calls, or home
visits from health workers) so as to maximize engagement, given the history of
the beneficiaries’ behaviour. Thus, our contributions are:

  1. 1.

AVM time-slot optimization. We develop a UCB-based (Auer, Cesa-Bianchi, and
Fischer 2002) approach that improves the timing of AVMs according to
beneficiaries’ listenership patterns. This work is the _first_ to optimize
such time slot selection.

  2. 2.

Intervention planning. In absence of rich demographic features in Kilkari, we
leverage beneficiaries’ listenership history and model the problem as a non-
markovian Time-Series Bandit (TSB) to determine which beneficiaries should
receive which intervention (such as a visit by an ASHA worker or an automated
voice call reminder), given our limited intervention budget. This work is the
first to model multiple inteventions via TSB.

  3. 3.

Secondary evaluation with data from 4000 real-world beneficiaries. We use
anonymized real-world Kilkari data from Odisha state in India to show CHAHAK’s
effectiveness in identifying beneficiaries’ preferred time slot, and in
harnessing multiple interventions to boost listenership.

When deployed, CHAHAK will benefit millions of underserved and marginalized
women and children in India. While we focus on Kilkari, the proposed
methodological contributions have the potential to _perform well on any other
multi-action, non-Markovian bandit setting_.

## Discussion and Related Work

Adherence to healthcare programs is a well-studied problem in a plethora of
domains such as cardiac problems (Son et al. 2010; Corotto et al. 2013),
tuberculosis (Killian et al. 2019; Pilote et al. 1996), HIV (Tuldrà et al.
1999), elderly care (Pollack et al. 2002), lifestyle changes (Liao et al.
2020), and more. These works either perform one-shot prediction of high drop-
out risk beneficiaries, or involve sequential decision making where one has to
select a subset of beneficiaries to intervene at each time-step given limited
resources (number of interventions), in order to maximize adherence.
Specifically on the latter, state-of-the-art approaches tend to model the
sequential allocation problem as a Restless Multi-Armed Bandit (RMAB) problem
(Nishtala et al. 2020; Mate et al. 2020, 2022).

SAHELI (Verma et al. 2023) is the first deployed application of RMAB in public
health, and is in continuous use by ARMMAN for mMitra, the maternal mHealth
program operating in the Indian city of Mumbai. Despite a similar setting,
there are some notable limitations of SAHELI. First, SAHELI assumes that
beneficiaries follow the Markov property, which might not be the case in many
applications. In fact, in a section ahead we show that beneficiaries’
behaviour in Kilkari is co-related to their past listenership patterns, which
violates the Markov assumption. Additionally, SAHELI requires access to a rich
set of demographic features which are key to transferring knowledge from
historic data to newly enrolled beneficiaries. Instead, CHAHAK only has access
to listenership logs. Hence utilizing these is crucial in predicting
beneficiaries’ behaviour. Moreover, SAHELI relied on live service calls as the
only intervention, however there are multiple interventions possible in
CHAHAK, with varying impact and budget. Inspired by (Danassis et al. 2023), we
introduce a novel Multi-action, Non-Markovian Time Series Bandit Model
formulation for the problem.

Finally, another key difference between CHAHAK and past work (e.g., (Verma et
al. 2023; Danassis et al. 2023)) is that CHAHAK attempts to reduce the number
of _automatic dropouts_ , by optimizing for each beneficiary’s preferred time
slots for receiving AVMs. AVMs are not currently optimized as per beneficiary
time-slot preferences due to lack of knowledge of phone ownership, working
hours, or any other indicator on time-slot preference (Bashingwa et al. 2023).
This aggravates the problem of drop-offs whereby if a beneficiary has less
than 25% engagement in the past 6 weeks, they are automatically dropped off
from the system. This automatic dropoff is imposed due to the scale at which
Kilkari operates; serving millions of beneficiaries imposes infrastructural
and logistic limitations, and thus critical bandwidth can not be spent on AVMs
to incorrectly recorded numbers or ones no longer in service.

## CHAHAK Overview

Kilkari sends 72 Automated Voice Messages (AVMs), one per week, on topics such
as maternal and child health, immunization, family planning, etc. to pregnant
women and new mothers to share critical health information starting from the
second trimester of pregnancy until the child is one year old. To increase the
reach of these messages, a beneficiary receives up to a maximum of 9 call
attempts across 4 days (until they pick up) to maximise the probability of
receiving the message for the intended week.

Despite multiple call attempts, our secondary analysis (see Fig. 1 (a)) shows
that approximately 23% of beneficiaries are not reached on average. Moreover,
there is evidence of beneficiaries having preferences over the call time-
slots, which could contribute to reduced listenership and adherence (Bashingwa
et al. 2021). These are crucial challenges to address since beneficiaries who
listen to less than 25% of the AVM for 6 weeks in a row, are automatically
dropped out of the program – losing access to critical health information.
Secondary analysis also shows that more than 50% of beneficiaries enrolled in
Kilkari in January 2022 in Odisha were dropped off within 6 months (Fig. 1
(b)). It is important to note that new beneficiaries are constantly being
enrolled in the system, and hence AVM delivery bandwidth does not go
unutilized. This motivates the need for optimizing preferred time slots,
especially for those whose listenership is mainly impacted due to their
availability to answer calls. Besides AVM time slot optimization,
interventions continue to assist and retain beneficiaries in the system, and
are of particular importance to beneficiaries who need encouragement and
support. CHAHAK ensures that interventions are directed towards beneficiaries
who stand to gain the most, by optimizing both AVM time slots and intervention
planning.

Figure 1: (a) Around 23% beneficiaries are unreached despite 9 attempts (b)
Number of beneficiaries (from the set of those enrolled in January 2022)
remaining in the program after each month

### AVM Time Slot Optimization

Optimizing call slots for AVMs is essential to improve engagement. Currently,
Kilkari uses a random calling approach, attempting to reach beneficiaries at
seven different timeslots throughout the day. Instead, CHAHAK optimizes AVM
calls to maximize the likelihood of pick up based on past behaviours.

### Intervention Planning

CHAHAK has two possible interventions available - ASHA visit, and an automated
phone call reminder. Given the past history of listenership, beneficiaries are
selected for an intervention based on the urgency, i.e, how long they are
predicted to remain in an engaging state, to avoid drop-offs as well as boost
listenership.

ASHA visit scheduling: An ASHA (Accredited Social Health Activist) is a
trained female community health worker who raises awareness about healthcare
and bridges the gap between the community and the public health system. Among
their many duties, ASHAs also help bring attention to the Kilkari program and
register pregnant women or young mothers in their area into the program. In
cases when a beneficiary is predicted to be at a risk of dropout, an ASHA can
visit the beneficiary’s home to discuss the importance of Kilkari in their
pregnancy journey and ways to resolve any potential hindrances to their
listenership. However, due to the limited bandwidth of the relatively few ASHA
workers per district, there is a limit on how many beneficiaries can be
reached. In order to test our system in a demanding, low-resource scenario, we
opted for a $1\%$ (of the total number of beneficiaries) budget for the number
of ASHA interventions. This number was selected as a realistic testcase in
deliberation with our partner, ARMMAN.

Automated Voice Call reminders: The automated intervention calls highlight the
importance of Kilkari and remind the beneficiaries to keep listening. Similar
to the ASHA interventions, we also selected an $1\%$ (of the total number of
beneficiaries) budget for the calls.

## Data Description

ARMMAN stores anonymized call logs of beneficiaries’ listenership of Kilkari
AVM’s. These logs contain information such as duration of call listened to,
number of attempts to reach the beneficiary and the date and time of call
attempts.

### Pickup Rates

To model the time slot planning algorithm, we estimate the pick-up rate from
anonymized call logs data of $4000$ beneficiaries who were registered in
Odisha and who had received a call in the first week of 2022. We consider $5$
time slots of two hours each from 8AM to 6PM plus one additional slot each for
before 8AM and after 6PM, resulting in 7 time slots total. The pickup rate for
time slot x is then estimated as the ratio between the number of calls picked
up, divided by the number of calls made (in that timeslot).

### Intervention Campaigns

Using anonymised call logs, we perform secondary analysis of awareness
campaigns conducted by ARMMAN in June 2023 in 6 blocks in the Indian state of
Odisha. In 3 blocks, the ASHA health workers were informed to visit randomly
selected beneficiaries and encourage them to listen to the AVMs. To confirm
that the visit happened, beneficiaries would send a missed call to a
predefined phone number, which helped identify the date of the ASHA visit. In
3 other blocks, beneficiaries were sent an automated voice call reminder to
encourage listening to the scheduled AVMs. Our secondary analysis also
considers 6 other blocks with similar demographics but without any campaigns
conducted. For the ASHA campaign, there were $148$ missed calls, while $1519$
Automated Voice Call Reminders were answered. The automated call logs provide
the listenership of these beneficiaries before and after the intervention.
These can be compared with the listenership trend of $3078$ beneficiaries from
the blocks which didn’t receive any interventions.

We observe a $16\%$ increase in listenership through ASHA visits and $5\%$
increase in listenership from CALL campaigns (across all intervened
beneficiaries). ASHA visits are, as expected, more effective because of the
personal interaction with the beneficiaries. We use this real data to train
the proposed system to learn an optimal allocation of intervention resources,
and to evaluate it in simulation to show its potential.

## Time Slot Scheduling

As mentioned, currently Kilkari employs a random calling approach, reaching
out to beneficiaries at one of seven different timeslots throughout the day.
Under this "retry" scheme, the system makes a maximum of $R$ call attempts per
week until the beneficiary picks up, where currently $R=9$. However, the
program does not remember which timeslot worked best for reaching the
beneficiary. Consequently, there is no regular anticipation of the call
timing, leading to potential wastage of a large number of retries.

To address this issue and reduce the number of retries, we propose
implementing a multi-armed bandit model – specifically Upper Confidence Bound
(UCB) – for each beneficiary. Here, each "arm" represents a specific timeslot.
"Pulling the arm" corresponding to timeslot $t$, implies sending the AVM at
timeslot $t$. After each AVM call attempt, we record whether the call was
answered or not. CHAHAK maintains independently its own estimates on the
empirical mean and upper confidence bounds of each arm (beneficiary). Let
$\mu_{i,j}$ denote the empirical mean payoff of slot $j$ for beneficiary $i$.
Let $\nu_{i,j}$ denote the number of times slot $j$ was selected by the
beneficiary $i$ so far. Since CHAHAK will run the algorithm independently for
each beneficiary, for simplicity we can drop the subscript $i$. The number of
pulls (or AVMs to beneficiary) each week is determined by the number of
retries $R$, pulling up to $R$ slots until a call is successfully answered.
The timeslot ($a_{\tau}$ at timestep $\tau$) is chosen based on the following
UCB equation:

$a_{\tau}=\arg\max_{j}\mu_{j}+\sqrt{2(\log\tau)/\nu_{j}}$

## Intervention Planning

In a RMAB setting, a planner can act (intervene) on $k$ out of $N$ arms
(beneficiaries) at each timestep (due to budget constraints). Each arm is
associated with a state, which transitions according to a Markovian transition
function at each time-step. Several RMAB based solutions have been
successfully developed (Mate et al. 2022; Verma et al. 2023; Wang* et al.
2023; Hsu 2018; Kadota et al. 2016), all of which, assume that transition
probabilities follow the Markov property. Recent work, though, has
demonstrated the importance of incorporating history (Danassis et al. 2023),
to capture _temporal dependencies_ in human behaviour. Unlike (Verma et al.
2023), Kilkari also lacks the rich demographic features and hence CHAHAK
relies on the use of call history data. Moreover, the section below validates
that beneficiaries’ behaviour in the Kilkari program is indeed history-
dependent (non-Markovian). Therefore, we propose a non-markovian Time-Series
Bandit model for multi-action setting to allocate health resources in Kilkari.

### Non-Markovian Behaviour in Kilkari

In this section we demonstrate the benefit for a non-Markovian (or higher
order Markov) model in explaining the beneficiaries’ behaviour in Kilkari. For
an $h$-order Markov process, the next state depends on the preceding $h$
states. More formally, let $x_{1},\dots,x_{t}$ be the elements of the process.
For a Markov process of order $h$, $(X_{t})_{t=1}^{\infty}$, it holds that
$\forall t$, $\Pr[X_{t}=x_{t}\mid
x_{t-1},x_{t-2},\ldots,x_{1}]=\Pr[X_{t}=x_{t}\mid
x_{t-1},x_{t-2},\ldots,x_{t-h}]$. Hence $h=1$ represents a traditional MDP,
while for $h>1$, an $h$-order Markov process can be viewed as a first order
Markov process on an expanded state space $s^{\prime}=\times_{i=1}^{h}s_{i}$.

As a measure of how closely the model can explain the beneficiaries’
behaviour, we compute the log-liklihood of observing the historical
trajectories of each beneficiary on Kilkari assuming an underlying MDP of
order $h$. This is easily done by calculating empirical transition
probabilities (assuming MDP of order $h$), and then computing the log-
likelihood of each trajectory. $l(h)\triangleq-\log\mathcal{L}(h\mid
x)\triangleq-\log\Pr(X=x\mid\text{model of order $h$})$. Finally, we average
over all trajectories. Figure 2 shows the relative improvement
$-\left(\frac{l(h)-l(h=1)}{l(h=1)}\right)$ in negative log-likelihood for
order $h$ Markov processes, with up to $15\%$ improvement for $h=6$.

Figure 2: Relative (with respect to $h=1$) improvement in log likelihood. The
probability of observing the historical trajectories in our dataset increases
as we increase the order $h$ of the underlying Markov model. This suggests
non-Markovian behaviour for the beneficiaries in Kilkari.

### Multi-Action Time Series Bandits

We model the intervention scheduling as a Multi-action TSB problem
$\\{\mathcal{N},\mathcal{K},(X^{n\in N}_{t})_{t=1}^{H},R\\}$, where
$\mathcal{N}$ is the set of _independent_ arms, $\mathcal{K}$ is a set of
budget constraints – one per type of action $i$ –
$(X^{n\in\mathcal{N}}_{t})_{t=1}^{H}$ is the associated transition process for
arm $n$ for time horizon $H$, and $R:(X_{t})_{t=1}^{H}\rightarrow\Re$ is the
reward function. The goal is to maximize the cumulative reward over the time
horizon $H$. In our problem setting, the reward for a beneficiary at every
timestep is the duration of call listened to in a week.

Each arm $n$ follows a history-dependent (i.e., non-Markovian) transition
process. Let $s_{n,t}\in\mathcal{S}$, and $a_{i,n,t}\in\mathcal{A}$ denote the
state and action taken on arm $n$, respectively, in timestep $t$. States are
continuous in $[0,1]$, and represent the ‘level of engagement’ of a
beneficiary, with higher numbers representing a higher level of listenership
in seconds. States are fully observable. The action set consists of three
actions: (i) a visit by an ASHA worker, (ii) an automated voice call or (iii)
a no-cost passive action.

### MA-TARI: A New Arm Ranking Index

We follow a two stage process for planning which arms to choose for a
particular action (intervention). First, we learn a model of beneficiaries’
future listenership based on past behaviour through a Time-Series Forecasting
(TSF) model. Next for each arm $n$ independently, we calculate two quantities
recursively using our TSF model for each intervention. (i) The number of weeks
$u_{n,i}$ it will take for beneficiary $n$ to fall below a threshold of
listenership if we act once (using one of our intervention actions) on arm $n$
using intervention $i$ and never act again; (ii) The number of weeks $v_{n}$
it will take for beneficiary $n$ to fall below a threshold of listenership if
we never act (passive action at each timestep). Then the Multi-Action Time-
series Arm Ranking Index (MA-TARI) for arm $n$ and intervention $i$ is given
by:

$m_{n,i}=u_{n,i}/v_{n}$

Since we independently compute one index per intervention, hereafter we will
use TARI (instead of MA-TARI) when we refer to the index of a specific
intervention. To optimize the allocation of interventions to beneficiaries
given the budget per intervention, we propose and compare four policies for
the multi-action setting:

  1. 1.

MA-TARI Integer Linear Programming : We solve an ILP to find the optimal
action $i$ for every arm $n$, such that we maximise the sum of MA-TAR indices
$m_{n,i}$ (i.e., cumulative expected benefit from the interventions) while
respecting the budget (first constraint) and making sure each beneficiary will
only receive one intervention (second constraint). Intuitively, maximizing the
sum of MA-TARI, optimizes the number of weeks the beneficiaries remain
engaging with the program (based on our TSF model predictions), and thus the
total listenership as well. We solve the following program (where $y_{n,i}$
denote the selection variables):

$\displaystyle\max_{y_{n,i}\in\\{0,1\\}}\ \ $
$\displaystyle\sum_{n,i}y_{n,i}m_{n,i}$ $\displaystyle s.t.\ \ $
$\displaystyle\sum_{n}y_{n,i}\leq k_{i},\forall i$
$\displaystyle\sum_{i}y_{n,i}\leq 1,\forall n$

  2. 2.

MA-TARI Greedy: Solving an LP for a large number of arms can be
computationally challenging. As such, we also propose a greedy approach, which
is more computationally tractable, especially for a low resource NGO. In the
greedy policy, for each intervention $i$, we select the top $k_{i}$ arms with
the highest TAR index without considering other interventions. If an arm is
selected by two separate interventions, then we favor the one with the higher
index. As demonstrated in the Experiments and Results Section, in this
context, the computationally tractable greedy approach performs on par with
the ILP.

  3. 3.

Random : Randomly choosing $k_{i}$ beneficiaries per intervention $i$.

  4. 4.

Control : No interventions were given to this group.

Finally, note that while we focus on Kilkari as our application domain, _the
proposed MA-TARI is potentially applicable to any multi-action, non-Markovian
RMAB problem_.

## Experiments and Results

### Time Slot Planning

Figure 3 evaluates the effectiveness of UCB in learning an optimised time slot
for each beneficiary while adhering to the retry scheme $R=9$. We average over
30 runs for each beneficiary, and over 4000 beneficiaries.

Figure 3 [left] illustrates how UCB effectively reduces the number of calls
needed to reach a beneficiary each week. A uniformly random algorithm that has
no update on the timeslots’ pickup rates – denoted Uniform (no update)) –
requires an average of 4.86 calls per week for a successful call, while UCB
accomplishes the same results by cutting the number of retries in almost
_half_ , using as few as 2.5 calls per week by week 6.

To further illustrate the effectiveness of UCB, we also compare it with
Uniform (update), where we give the uniform random algorithm the additional
ability to update its mean of call success rate per slot. Figure 3 [right]
demonstrates the effectiveness of using UCB in learning the true pickup rate
of the best timeslot. Within 4 weeks, UCB achieves a pickup rate difference of
less than 0.15, compared to the true pickup rate, for 75% of the population.
In contrast, the Uniform (update) approach takes 8 weeks to reach a similar
level of accuracy.

Figure 3: [Left] Average number of calls performed each week for each
beneficiary until a call is successfully picked up. [Right] Fraction of
beneficiaries for which the learned pickup rate of the best slot has converged
(within a difference of $\leq 0.15$ from the true mean).

### Intervention Planning

The call logs provide us with beneficiaries listenership behaviour before and
after the intervention (ASHA visit or CALL), and thereby each beneficiary’s
listenership trajectory as well as the action at every timestep.

#### Time Series Forecasting Models

We split the time series trajectories of beneficiaries into a 70-30 train-test
split. We compared a range of time series forecasting models, such as LSTMs,
BiLSTMS, Transformers, and adding time-based vector representations and
attention layers to predict future listenership. Specifically, all these
models take as input past states (duration of call listened to) and actions
(ASHA intervention, CALL intervention, or no intervention) to predict
listenership in the next timestep (output).

Comparing the performance of these models revealed that all models have a
similar Mean Absolute Error. We conducted this experiment over a range of
context lengths, and observed similar results. On average, the maximum
deviation of MAE of any model from the MAE of the LSTM model was less than
0.01. Thus we opted to use the simple LSTM architecture for computational
efficiency as that would be preferred by our NGO partner which operates in a
low resource environment. Our chosen model achieved MAE of $0.23$ on average
for one step prediction on real data. We used a history of $h=4$ weeks for
training and evaluation.

#### Multi-Action TAR

We use the TSF model to compute the MA-TAR index. In order to test its
efficacy, we use a ground-truth behaviour simulator which simulates the true-
listenership given the action suggested by the policy and past listenership.
This ground-truth simulator is also a Time-Series Forecast model but is fitted
on the entire dataset.

In Fig. 4, we show how the two MA-TARI policies and the Random policy perform
with respect to the control policy. Specifically, we measure the cumulative
increase in listenership (in hours) every week, over the control group, for up
to 8 weeks. In the simulation, we consider 4000 beneficiaries and allow 1%
ASHA interventions and 1% CALL interventions every week. Fig. 4 shows that the
interventions significantly help in improving listenership over time. In
particular, the MA-TARI policy results in 49.6h and 49.47h of additional
content listened by the beneficiaries by the 8th week (for the LP and Greedy
variants, respectively). Compared to the 31.5h increase with the random
policy, this corresponds to an almost 57% increase in exposure to content.
Moreover, Fig. 5(a) shows the distribution of the MA-TAR indices for the ASHA
visit and call actions, highlighting the potential to delay the poor
listenership state by 3x weeks in the best case. As expected, the ASHA
intervention is more effective than the call reminder.

Another important metric is identifying and proactively preventing drop-outs
from the program. Fig. 5(b) depicts the percentage of drop-outs prevented
(cumulatively), compared to the control group. Again, MA-TARI – both the
greedy and LP variants – achieve a significant improvement in preventing drop-
outs (33% at the last week).

These results directly translate to real-world impact, as listening to content
improves health literacy, which ultimately leads to better health outcomes for
the mothers and their children. Importantly, note also that the greedy variant
performs on par to the cumbersome ILP formulation, thus providing a
computationally efficient, and scalable policy for a low resource organization
such as our NGO partner.

Figure 4: Comparison of different policies for planning multi-action
interventions as compared to Control group.

(a)

(b)

Figure 5: [Left] Distribution of Multi-Action TAR Indices for beneficiaries
that receive interventions. [Right] Prevention of dropouts compared to
control.

## Data Usage and Ethics

Acknowledging the responsibility associated with real-world AI systems for
undeserved communities, we have closely coordinated with domain experts from
our partner NGO, ARMMAN, throughout our analysis. This study falls into the
category of secondary analysis of the aforementioned anonymized dataset. There
is no demographic information available. We use the previously collected
engagement trajectories of different beneficiaries participating in the
Kilkari program to train the predictive model and evaluate the performance.
All the data collected through the program is owned by the NGO and only the
NGO is allowed to share data. The research group only accesses an anonymized
version of the data.

Bias and fairness There is no demographic data available for Kilkari
beneficiaries. Nonetheless, prior studies such as (Mohan et al. 2021) point
out that exposure to Kilkari helps improve health behaviors among the most
marginalised, thus helping to close the gap across some inequities in the
population. They also indicate that the more marginalised population benefits
from higher number of retries in Kilkari calls. CHAHAK potentially helps
reduce inequities by reducing the number of retries required, and by improving
listening of Kilkari messages particularly amongst low listeners.

## Path to Deployment

CHAHAK is intended to be deployed at a national scale in India. To that end,
the next step will involve a cohort study in one state in India to validate
CHAHAK’s benefits in the field. Following a trial deployment, the program will
then be gradually scaled up. Of course we expect such an endeavour to present
many new challenges, such as regional differences in listenership patterns, or
local constraints on the number and types of interventions. We aim to draw
from our experience in deploying SAHELI (Verma et al. 2023) – the system for
ARMMAN’s smaller scale maternal mHealth program, mMitra – in order to tackle
the challenges that will arise. Most importantly, though, just as with
CHAHAK’s development, all of the steps will be done in _close collaboration
with our partner ARMMAN_ ; with ARMMAN ultimately in charge of the actual
deployment.

## Conclusion

In this work, we present CHAHAK, a novel AI system designed to improve
listenership and adherence in Kilkari, the largest maternal mHealth program in
the world, with 3 million active beneficiaries at any given time. Key
innovations in CHAHAK include a novel multi-action, non-Markovian Time Series
Bandit model to optimise the allocation of interventions across beneficiaries
to improve engagement. We also presented a bandit algorithm for identifying
preferred time slots of beneficiaries. With real world anonymized call data
from the Indian state of Odisha, we show the benefits of our novel AI
techniques in terms of improved listenership and identification of preferred
beneficiary timeslots.

## Acknowledgments

This material is based upon work supported by the National Science Foundation
under Grant No. IIS-2229881. Any opinions, findings, and conclusions or
recommendations expressed in this material are those of the author(s) and do
not necessarily reflect the views of the National Science Foundation.

## References

  * Amagai et al. (2022) Amagai, S.; Pila, S.; Kaat, A. J.; Nowinski, C. J.; and Gershon, R. C. 2022.  Challenges in participant engagement and retention using mobile health apps: literature review.  _Journal of medical Internet research_ , 24(4): e35120.
  * ARMMAN (2019) ARMMAN. 2019.  Kilkari.  https://armman.org/kilkari/.  Accessed: 2023-4-5.
  * ARMMAN (2023) ARMMAN. 2023.  mMitra.  https://armman.org/mmitra/.  Accessed: 2023-12-06.
  * Auer, Cesa-Bianchi, and Fischer (2002) Auer, P.; Cesa-Bianchi, N.; and Fischer, P. 2002.  Finite-time analysis of the multiarmed bandit problem.  _Machine learning_ , 47: 235–256.
  * Bashingwa et al. (2021) Bashingwa, J. J. H.; Mohan, D.; Chamberlain, S.; Arora, S.; Mendiratta, J.; Rahul, S.; Chauhan, V.; Scott, K.; Shah, N.; Ummer, O.; et al. 2021.  Assessing exposure to Kilkari: a big data analysis of a large maternal mobile messaging service across 13 states in India.  _BMJ Global Health_ , 6(Suppl 5): e005213.
  * Bashingwa et al. (2023) Bashingwa, J. J. H.; Mohan, D.; Chamberlain, S.; Scott, K.; Ummer, O.; Godfrey, A.; Mulder, N.; Moodley, D.; and LeFevre, A. E. 2023.  Can we design the next generation of digital health communication programs by leveraging the power of artificial intelligence to segment target audiences, bolster impact and deliver differentiated services? A machine learning analysis of survey data from rural India.  _BMJ open_ , 13(3): e063354.
  * Corotto et al. (2013) Corotto, P. S.; McCarey, M. M.; Adams, S.; Khazanie, P.; and Whellan, D. J. 2013.  Heart failure patient adherence: epidemiology, cause, and treatment.  _Heart failure clinics_ , 9(1): 49–58.
  * Danassis et al. (2023) Danassis, P.; Verma, S.; Killian, J. A.; Taneja, A.; and Tambe, M. 2023.  Limited Resource Allocation in a Non-Markovian World: The Case of Maternal and Child Healthcare.  In Elkind, E., ed., _Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence, IJCAI-23_ , 5950–5958. International Joint Conferences on Artificial Intelligence Organization.  AI for Good. https://doi.org/10.24963/ijcai.2023/660.
  * Hsu (2018) Hsu, Y.-P. 2018.  Age of information: Whittle index for scheduling stochastic arrivals.  In _2018 IEEE International Symposium on Information Theory (ISIT)_ , 2634–2638. IEEE.
  * Kadota et al. (2016) Kadota, I.; Uysal-Biyikoglu, E.; Singh, R.; and Modiano, E. 2016.  Minimizing the age of information in broadcast wireless networks.  In _2016 54th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , 844–851. IEEE.
  * Killian et al. (2019) Killian, J. A.; Wilder, B.; Sharma, A.; Choudhary, V.; Dilkina, B.; and Tambe, M. 2019.  Learning to Prescribe Interventions for Tuberculosis Patients Using Digital Adherence Data.  _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_.
  * Liao et al. (2020) Liao, P.; Greenewald, K.; Klasnja, P.; and Murphy, S. 2020.  Personalized heartsteps: A reinforcement learning algorithm for optimizing physical activity.  _Proceedings of the ACM on Interactive, Mobile, Wearable and Ubiquitous Technologies_ , 4(1): 1–22.
  * Mate et al. (2020) Mate, A.; Killian, J. A.; Xu, H.; Perrault, A.; and Tambe, M. 2020.  Collapsing Bandits and Their Application to Public Health Interventions.  In _Advances in Neural and Information Processing Systems (NeurIPS)_ , volume 34.
  * Mate et al. (2022) Mate, A.; Madaan, L.; Taneja, A.; Madhiwalla, N.; Verma, S.; Singh, G.; Hegde, A.; Varakantham, P.; and Tambe, M. 2022.  Field study in deploying restless multi-armed bandits: Assisting non-profits in improving maternal and child health.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, 12017–12025.
  * Mohan et al. (2021) Mohan, D.; Scott, K.; Shah, N.; Bashingwa, J. J. H.; Chakraborty, A.; Ummer, O.; Godfrey, A.; Dutt, P.; Chamberlain, S.; and LeFevre, A. E. 2021.  Can health information through mobile phones close the divide in health behaviours among the marginalised? An equity analysis of Kilkari in Madhya Pradesh, India.  _BMJ Global Health_ , 6(Suppl 5): e005512.
  * Nishtala et al. (2020) Nishtala, S.; Kamarthi, H.; Thakkar, D.; Narayanan, D.; Grama, A.; Hegde, A.; Padmanabhan, R.; Madhiwalla, N.; Chaudhary, S.; Ravindran, B.; et al. 2020.  Missed calls, Automated Calls and Health Support: Using AI to improve maternal health outcomes by increasing program engagement.  _arXiv preprint arXiv:2006.07590_.
  * Pilote et al. (1996) Pilote, L.; Tulsky, J. P.; Zolopa, A. R.; Hahn, J. A.; Schecter, G. F.; and Moss, A. R. 1996.  Tuberculosis Prophylaxis in the Homeless: A Trial to Improve Adherence to Referral.  _Archives of Internal Medicine_ , 156(2): 161–165.
  * Pollack et al. (2002) Pollack, M. E.; Brown, L.; Colbry, D.; Orosz, C.; Peintner, B.; Ramakrishnan, S.; Engberg, S.; Matthews, J. T.; Dunbar-Jacob, J.; McCarthy, C. E.; et al. 2002.  Pearl: A mobile robotic assistant for the elderly.  In _AAAI workshop on automation as eldercare_ , volume 2002. AAAI, 2002, Edmonton, Alberta, Canada.
  * Son et al. (2010) Son, Y.-J.; Kim, H.-G.; Kim, E.-H.; Choi, S.; and Lee, S.-K. 2010.  Application of support vector machine for prediction of medication adherence in heart failure patients.  _Healthcare informatics research_ , 16(4): 253–259.
  * Tuldrà et al. (1999) Tuldrà, A.; Ferrer, M. J.; Fumaz, C. R.; Bayés, R.; Paredes, R.; Burger, D. M.; and Clotet, B. 1999.  Monitoring Adherence to HIV Therapy.  _Archives of Internal Medicine_ , 159(12): 1376–1377.
  * United Nations Women (2023) United Nations Women. 2023.  SDG-3.  Accessed: March 31, 2023.
  * Verma et al. (2023) Verma, S.; Singh, G.; Mate, A.; Verma, P.; Gorantala, S.; Madhiwalla, N.; Hegde, A.; Thakkar, D.; Jain, M.; Tambe, M.; and Taneja, A. 2023.  Increasing Impact of Mobile Health Programs: SAHELI for Maternal and ChildCare.  In _Innovative Applications of Artificial Intelligence (IAAI)_.
  * Wang* et al. (2023) Wang*, K.; Verma*, S.; Mate, A.; Shah, S.; Taneja, A.; Madhiwalla, N.; Hegde, A.; and Tambe, M. 2023.  Scalable Decision-Focused Learning in Restless Multi-Armed Bandits with Application to Maternal and Child Health.  In _Proceedings of the AAAI Conference on Artificial Intelligence (AAAI)_.
  * World Health Organization (2023) World Health Organization. 2023.  Maternal Mortality.  Accessed: March 31, 2023.