# Swarming bottom feeders: Flocking at solid-liquid interfaces

Niladri Sarkar<EMAIL_ADDRESS>Instituut-Lorentz, Leiden University,
P.O. Box 9506, 2300 RA Leiden, The Netherlands  Abhik Basu
<EMAIL_ADDRESS><EMAIL_ADDRESS>Condensed Matter Physics Division,
Saha Institute of Nuclear Physics, Calcutta 700064, West Bengal, India  John
Toner<EMAIL_ADDRESS>Department of Physics and Institute of Theoretical
Science, University of Oregon, Eugene, Oregon 97403, USA

###### Abstract

We present the hydrodynamic theory of coherent collective motion (“flocking”)
at a solid-liquid interface, and many of its predictions for experiment. We
find that such systems are stable, and have long-range orientational order,
over a wide range of parameters. When stable, these systems exhibit “giant
number fluctuations”, which grow as the 3/4th power of the mean number. Stable
systems also exhibit anomalous rapid diffusion of tagged particles suspended
in the passive fluid along any directions in a plane parallel to the solid-
liquid interface, whereas the diffusivity along the direction perpendicular to
the plane is not anomalous. In the remaining parameter space, the system
becomes unstable.

Many “active” systems consist of macroscopically large numbers of self-
propelled particles that align their directions of motion. This occurs both in
living kruse04 ; kruse05 ; goldstein13 ; saintillan08 ; hatwalne04  and
synthetic saha14 ; cates15 ; narayan07 ; lubensky09 ; marchetti08  systems.
Such “active orientationally ordered phases” exhibit many phenomena impossible
in their equilibrium analogs (e.g., nematics deGennes ), including spontaneous
breaking of continuous symmetries in two dimensionsvicsek95 ; tonertu95 ;
toner98 ; toner05 , instability in the extreme Stokesian limitsimha2002 , and
giant number fluctuations Chate+Giann ; toner2019giant ; ramaswamy03 .

“Dry” active systems - i.e., those lacking momentum conservation due to, e.g.,
friction with a substrate wolgemuth2002 ; toner98 ; ramaswamy03 \- behave
quite differently from “wet” active fluids (i.e., those with momentum
conservation) lushi2014 .

In this paper, we present the first theory of a natural hybrid of these two
cases: polar active particles at a solid-liquid interface (see figure (1)). We
are motivated by experiments schaller13  in which highly concentrated actin
filaments on a solid-fluid interface are propelled by motor proteins, and
those of Bricard et al bricard2013 ; Geyer17 , who studied the emergence of
macroscopically directed motion in “Quincke rollers”. The latter are motile
colloids, spontaneously rolling on a solid substrate when a sufficiently
strong electric field is applied.

These systems differ from both dry and wet active matter, as defined above, by
having both friction from the underlying solid substrate and the long range
hydrodynamic interactions due to the overlying bulk passive fluid.

The geometry we consider here, as in Ref. schaller13 ; bricard2013 , places a
collection of polar, self-propelled particles at the flat interface (the
$x$-$y$ plane of our coordinate system) between a solid substrate and a semi-
infinite bulk isotropic and incompressible passive liquid. as illustrated in
Fig. 1. We consider the extreme Stokesian limit, in which inertial forces are
completely negligible compared to viscous forces.

Figure 1: (Color online) Schematic diagram of our system: a layer of active
polar particles moving on a solid substrate with a passive ambient (“bulk”)
fluid above.

The most surprising result of our work is that, even in the presence of noise,
this system can be in a stable, long-range ordered polar state, in sharp
contrast to “wet” active systems, which are generically unstable simha2002  at
low Reynolds number, and equilibrium systems, which cannot display long range
orientational order in two dimensions at finite temperatureMW ; xtalfoot ;
2dxtal ; teth .

Remarkably, this ordered state is predicted even by a linear theory.
Furthermore, this linear theory provides an asymptotically exact long
wavelength description, in contrast to dry polar active systems, which can
only be correctly described by a non-linear theory. Indeed, dry polar active
systems can only exhibit long range order due to non-linear effects vicsek95 ;
tonertu95 ; toner98 ; toner05 .

Concomitant with the long-range polar order, the density fluctuations are
giant: the standard deviation $\sqrt{\langle(N-\langle N\rangle)^{2}\rangle}$
of the number $N$ of the active particles contained in a fixed open area
scales with its average $\langle N\rangle$ according to

$\sqrt{\langle(N-\langle N\rangle)^{2}\rangle}\propto\langle
N\rangle^{3/4}\,.$ (1)

This agrees very well with the experiments of schaller13 , which found
$\sqrt{\langle(N-\langle N\rangle)^{2}\rangle}\propto\langle N\rangle^{0.8}$.
Note that our prediction should not be confused with qualitatively similar
predictions for dry active matter Chate+Giann ; GNF  and active nematics AN ,
for which the exponent is different, because they belong to different
universality classes.

We also find that the fluctuations in the active fluid layer stir the bulk
fluid above it, making the diffusion of a passive tagged particle parallel to
the active fluid layer anomalous: specifically, the mean squared displacement
grows with time $t$ as $t\ln t$, whereas the diffusive motion perpendicular
to the active fluid layer remains conventional, i.e., the mean squared
displacement scales like $t$.

To understand the physics of this system, we have constructed a theory which,
when linearized for small fluctuations about a uniform reference state, is
asymptotically exact in the long wavelength limit, and gives the above
results. We define $\hat{{\bf p}}(\mathbf{r}_{{}_{\parallel}},t)$ as the
coarse grained polarization of the active particles, and
$\rho(\mathbf{r}_{{}_{\parallel}},t)$ as the conserved areal density of the
active polar particles on the surface. Taking our uniform reference state to
be $\hat{{\bf p}}(\mathbf{r}_{{}_{\parallel}},t)=\hat{{\bf x}}$ (see Fig. 1),
and $\rho=\rho_{0}$, one hydrodynamic variable is the transverse fluctuations
$p_{y}$ of $\hat{{\bf p}}(\mathbf{r}_{{}_{\parallel}},t)$, which we take to
have unit magnitude, i.e., $|\hat{{\bf p}}|^{2}=1$. This is a non-conserved
broken symmetry - i.e., “Goldstone” - mode. Our second hydrodynamic variable
is the fluctuations
$\delta\rho(\mathbf{r}_{{}_{\parallel}},t)\equiv\rho(\mathbf{r}_{{}_{\parallel}},t)-\rho_{0}$
of the density from its mean value.

These variables couple to the bulk passive fluid velocity
$\mathbf{v}(\mathbf{r}_{{}_{\parallel}},z,t)$ via an active boundary condition
given below in (16). Eliminating $\mathbf{v}(\mathbf{r}_{{}_{\parallel}},z,t)$
by solving the Stokes equation for the bulk fluid subject to this active
boundary condition gives the equations of motion for the spatially Fourier
transformed fields $p_{y}(\mathbf{q},t)$ and $\delta\rho(\mathbf{q},t)$:

$\displaystyle\partial_{t}\delta\rho(\mathbf{q},t)=-iv_{\rho}[q_{x}\delta\rho(\mathbf{q},t)+\rho_{c}q_{y}p_{y}]+i\mathbf{q}\cdot{\bf
f}_{\rho}(\mathbf{q},t)\,,$ (2) $\displaystyle\partial_{t}p_{y}(\mathbf{q},t)$
$\displaystyle=$ $\displaystyle-
iv_{p}q_{x}p_{y}(\mathbf{q},t)-\gamma\left({q^{2}+q_{y}^{2}\over
q}\right)p_{y}(\mathbf{q},t)-\left({\gamma_{\rho}\over\rho_{c}}\right)\left({q_{x}q_{y}\over
q}\right)\delta\rho(\mathbf{q},t)-i\sigma_{t}q_{y}\delta\rho(\mathbf{q},t)+f_{y}(\mathbf{q},t)\,,$
(3)

where $v_{\rho}$, $v_{p}$, $\gamma$, $\gamma_{\rho}$, $\rho_{c}$, and
$\sigma_{t}$ are parameters of our model. Note the non-analytic character of
the damping $\gamma$ and $\gamma_{\rho}$ terms in (3); due to long-ranged
hydrodynamic interactions mediated by the bulk passive fluid.

In (2) and (3), ${\bf f}_{\rho}$ and $f_{y}$ are zero-mean Gaussian white
noises whose variances are parameters of our model.

For stability, fluctuations must decay for all directions of $\mathbf{q}$. We
show in the associated long paper (ALP) alp  that this condition is satisfied
provided that the analogs of the bulk compressibility and the shear and bulk
viscosities in our system are all positive, and that the coupling of the
density of the active particles to their self-propelled speeds is not too
strong.

Thus, in contrast to “wet” active matter in the “Stokesian” limit simha2002 ;
toner05 , our “mixed” system can be generically stable. Indeed, the
requirements for stability are almost as easily met for these systems as for
an equilibrium fluid. Furthermore, when the stability conditions are met,
fluctuations about the uniform ordered state in this model decay with a rate
that scales linearly with $q$, quite different from the linear theory of dry
active matter. The also propagate nondispersively with a wavespeed independent
of $q$.

This unusual damping in this linear theory is responsible for many novel
phenomena: most strikingly, it makes $\langle p_{y}^{2}({\bf
r}_{\perp},t)\rangle$ asymptotically independent of the lateral size of the
system, a tell-tale signature of orientational long-range order. It also leads
to giant number fluctuations of the active particles given by (1), as
mentioned earlier.

In the ordered state, the active particles “stir” the passive fluid above
them. The mean squared components $\langle
v_{x}^{2}(\mathbf{r}_{{}_{\parallel}},z,t)\rangle$, and $\langle
v_{y}^{2}(\mathbf{r}_{{}_{\parallel}},z,t)\rangle$ of the passive fluid
velocity field $\mathbf{v}(\mathbf{r}_{{}_{\parallel}},z,t)$ thereby induced
are inversely proportional to the distance $z$ from the solid-fluid interface.

The unequal-time correlations $\langle
v_{x,y}(\mathbf{r}_{{}_{\parallel}},z,t)v_{x,y}(\mathbf{r}_{{}_{\parallel}},z,0)\rangle$
of the in-plane velocity fluctuations of the passive fluid also exhibit long
temporal correlations, which decay as $1/t$, whereas the correlation $\langle
v_{z}(\mathbf{r}_{{}_{\parallel}},z,t)v_{z}(\mathbf{r}_{{}_{\parallel}},z,0)\rangle$
of the the bulk fluid velocity perpendicular to the surface decays as
$1/t^{3}$.

The correlations of the in-plane velocity in turn lead to anomalous diffusion
of neutrally buoyant passive particles in the $x$\- and $y$-direction, with
variances of the displacements growing faster with time than the linear
dependence found for simple brownian particles. Specifically, we find, for a
particle that is initially a height $z_{0}$ above the solid-liquid interface,:

$\displaystyle\langle(r_{i}(t)-r_{i}(0))^{2}\rangle=\left\\{\begin{array}[]{ll}2D_{i}t\left[\ln\left({v_{0}t\over
z_{0}}\right)+O(1)\right]\,\,,\,\,\,\,\,\,t\ll{z_{0}^{2}\over D_{z}}\,,\\\ \\\
D_{i}t\left[\ln\left({v_{0}^{2}t\over
D_{z}}\right)+O(1)\right]\,\,,\,\,\,\,\,\,t\gg{z_{0}^{2}\over
D_{z}}\,,\end{array}\right.$ (7) (8)

where $i=x,y$, $v_{0}$ is a system-dependent characteristic speed (roughly
speaking, the self-propulsion speed of the active particles), $z_{0}$ is the
initial distance from the surface, and $D_{x,y,z}$ are diffusion constants
which are independent of $z_{0}$. Note that the mean square displacements
depend on the initial height $z_{0}$ for short times $t\ll z_{0}^{2}/D_{z}$,
but not for long times $t\gg z_{0}^{2}/D_{z}$. The crossover between these
limits is the time $t=z_{0}^{2}/D_{z}$ it takes for a neutrally buyoant
particle to diffuse a distance $z_{0}$ in the $z$-direction.

Diffusion in the $z$-direction remains conventional, controlled by a
$z$-independent diffusivity.

This set of predictions could also be tested experimentally by particle
tracking of neutrally buoyant tracer particles in the passive fluid.

Particles denser than the passive fluid, which therefore sediment, will also
be affected by this activity induced flow. We find that particles sedimenting
at a speed $v_{\rm sed}\ll v_{0}$ from an initial height $z_{0}$ will, when
they reach the surface, be spread out over a region of RMS dimensions
$\sqrt{\langle(x(z=0)-x(z=z_{0}))^{2}\rangle}$ and
$\sqrt{\langle(y(z=0)-y(z=z_{0}))^{2}\rangle}$ in the $x$ and $y$ directions,
respectively, with

$\langle(r_{i}(t)-r_{i}(0))^{2}\rangle=2D_{i}\left({z_{0}\over v_{\rm
sed}}\right)\ln\left({v_{0}\over v_{\rm sed}}\right)\ \ \ ,\ \ \ v_{\rm
sed}\ll v_{0}\,,$ (9)

where $v_{0}$ is roughly the mean speed of the active particles, and $v_{\rm
sed}$ is the speed at which the sedimenting particles sink.

Once again, these predictions should be readily testable in particle tracking
experiments.

We find that the polarization $\hat{{\bf p}}$, has a simple scaling form for
its spatio-temporally Fourier transformed correlation function:

$C_{pp}({\bf
q},\omega)\equiv\langle|p_{y}(\mathbf{q},\omega)|^{2}\rangle=\left({1\over
q^{2}}\right)F_{pp}\bigg{(}\left({\omega\over
q}\right),\theta_{\mathbf{q}}\bigg{)}\,,$ (10)

where the scaling function $F_{pp}(u,\theta_{\mathbf{q}})$ is given in the
ALP; and $\theta_{\bf q}\equiv\tan({q_{y}/q_{x}})$ is the angle between
$\mathbf{q}$ and the direction $\hat{{\bf x}}$ of the mean polarization. The
positions of the peaks in $C_{pp}(\mathbf{q},\omega)$ versus $\omega$ unreal
(but most definitely not their widths), are precisely those found for dry
active matter in tonertu95 ; toner98 ; toner05 ; i.e., $\omega_{\rm
peak}=c_{\pm}(\theta_{\mathbf{q}})q$, where $c_{\pm}(\theta_{\mathbf{q}})$ is
given by (19) and plotted in Figure (2).

Figure 2: (Color online) Polar plot of the sound speeds; the polarization
points directly to the right. That is, the distance along a straight line line
drawn from the origin and making an angle $\theta$ with the $x$-axis to its
intersection with the curve is proportional to the sound speed of a mode
propagating at the same angle $\theta$ to the mean polarization direction
$\hat{{\bf x}}$. There are two intersections for each such line, corresponding
to the two roots given in equation (19) for the sound speeds. Here we have
taken $v_{\rho}=1$, $v_{p}=c_{0}=2$, and $\gamma=.3$ (all in arbitrary units).

These peak positions agree with those found in the experiments of Geyer17  on
Quinke rollers.

The density-density correlation function
$C_{\rho\rho}(\mathbf{q},\omega)\equiv\langle|\delta\rho(\mathbf{q},\omega)|^{2}\rangle$,
and the density-polarization cross-correlation
$C_{p\rho}(\mathbf{q},\omega)\equiv\langle
p_{y}(\mathbf{q},\omega)\delta\rho(-\mathbf{q},-\omega)\rangle$, both obey
similar scaling laws, which are given in detail in the ALP.

Integrating these spatio-temporally Fourier-transformed correlation functions
over all frequencies $\omega$ shows that the equal time correlation functions
$C_{pp}(\mathbf{q})\equiv\langle|p_{y}(\mathbf{q},t)|^{2}\rangle$,
$C_{\rho\rho}(\mathbf{q})\equiv\langle|\delta\rho(\mathbf{q},t)|^{2}\rangle$,
and $C_{p\rho}(\mathbf{q})\equiv\langle
p_{y}(\mathbf{q},t)\delta\rho(-\mathbf{q},t)\rangle$ all scale like $1/q$.
Their dependence on the direction of $\mathbf{q}$ is given explicitly in the
ALP.

Fourier transforming these in space shows that the real space, equal-time
correlation functions $C_{pp}(\mathbf{r})=\langle
p_{y}(\mathbf{r}+\mathbf{R},t)p_{y}(\mathbf{R},t)\rangle$,
$C_{\rho\rho}(\mathbf{r})\equiv\langle\delta\rho(\mathbf{r}+\mathbf{R},t)\delta\rho(\mathbf{R},t)\rangle$,
and $C_{p\rho}(\mathbf{r})\equiv\langle
p_{y}(\mathbf{r}+\mathbf{R},t)\delta\rho(\mathbf{R},t)\rangle$ all scale like
$1/r$, and depend on the direction of $\mathbf{r}$. Explicit expressions for
this direction-dependence are given in the ALP.

These predictions could also be tested experimentally in systems in which the
active particles can be imaged, like those of schaller13 ; bricard2013 .
Although the anisotropy of the system ensures that all the correlators are
anisotropic functions of distance $\bf r$, nonetheless, their spatial scaling
remains isotropic. That is, the anisotropy exponent $\zeta$ that determines
the relative scaling between $x$ and $y$ is $\zeta=1$, in contrast to the
Toner-Tu model toner98 .

The correlator $C_{\rho\rho}(\mathbf{r}-\mathbf{r}^{\prime})$ can be used to
obtain the result (1) for the giant number fluctuations. The bulk velocity can
be obtained from $p_{y}(\mathbf{r},t)$ and $\delta\rho(\mathbf{r},t)$ through
the aforementioned solution of the Stokes equation subject to the active
boundary condition. This in turn allows us to derive the anomalous diffusion
(8); see the ALPalp  for detailed derivations.

We will now provide an outline of how we obtained these results. Details can
be found in the ALP.

In the presence of friction from the substrate, there is no momentum
conservation on the surface, so the only conserved variable on the surface is
the active particle number. We also include the bulk fluid velocity
$\mathbf{v}(\mathbf{r},t)$, which is defined throughout the semi-infinite
three dimensional (3D) space above the surface, since in that space momentum
(which is equivalent to velocity in the limit of an incompressible bulk fluid)
is conserved. However, we work in the Stokesian limit, in which viscous forces
dominate inertial ones.

We formulate the hydrodynamic equations for these variables by expanding their
equations of motion phenomenologically in powers of fluctuations of both
fields $\hat{{\bf p}}$ and $\rho$ from their mean values, and in spatio-
temporal gradients. In doing so, we respect all symmetries and conservation
laws of the underlying dynamics. In our non-equilibrium system, additional
equilibrium constraints like detailed balance do not apply. Our system has
underlying rotational invariance in the plane of the surface, which is
spontaneously broken by the active particles when they align their
polarizations.

Conservation of the active particles implies that
$\rho(\mathbf{r}_{{}_{\parallel}},t)$ obeys a continuity equation:

$\displaystyle\partial_{t}\rho+{\bm{\nabla}}_{s}\cdot{\bf J}_{\rho}$
$\displaystyle=$ $\displaystyle 0\,,$ (11)

where ${\bm{\nabla}}_{s}\equiv{\hat{\bf x}}\partial/\partial x+{\hat{\bf
y}}\partial/\partial y$ is the 2D gradient operator, with ${\hat{\bf x}}$ and
${\hat{\bf y}}$ the unit vectors along the $x$ and $y$ axis respectively. We
phenomenologically expand the active particle current ${\bf J}_{\rho}$ to
leading order in powers of the bulk velocity evaluated at the surface ${\bf
v}(\mathbf{r}_{{}_{\parallel}},z=0)$, and gradients, while respecting rotation
invariance. In practice, this means we can make the vector ${\bf J}_{\rho}$
only out of vectors the system itself chooses, i.e., out of gradients, the
surface velocity $\mathbf{v}_{s}(\mathbf{r}_{{}_{\parallel}},t)\equiv{\bf
v}(\mathbf{r}_{{}_{\parallel}},z=0,t)$, and the polarization $\hat{{\bf
p}}(\mathbf{r}_{{}_{\parallel}},t)$. These constraints force ${\bf J}_{\rho}$
to take the form:

$\displaystyle{\bf J_{\rho}}(\mathbf{r}_{{}_{\parallel}})$ $\displaystyle=$
$\displaystyle\rho_{e}(\rho,|\mathbf{v}_{s}|){\bf
v}_{s}(x,y)+\kappa(\rho,|\mathbf{v}_{s}|)\hat{{\bf p}}$ (12)

to leading order in gradients. The factor $\kappa(\rho,|\mathbf{v}_{s}|)$ is
an active parameter reflecting the self-propulsion of the particles through
interaction with the solid substrate, while the $\rho_{e}$ term reflects
convection of the active particles by the passive fluid above them. The
parameter $\rho_{e}\neq\rho$ in general due to drag between the active
particles and the substrate.

In calculating the bulk velocity
$\mathbf{v}(\mathbf{r}_{{}_{\parallel}},z,t)$, we assume the bulk fluid is in
the extreme “Stokesian” limit, in which inertia is negligible relative to
viscous drag. This should be appropriate for most systems in which the active
particles are microscopic, since the Reynolds’ number will be extremely low
for such particles. It is, however, certainly not valid for bottom-feeding
fish, so the title of this paper takes some poetic license!

In this limit, the three-dimensional (3D) incompressible bulk velocity field
${\bf v}=(v_{i},v_{z}),\,i=x,y$ satisfies the 3D Stokes’ equation

$\eta\nabla^{2}_{3}v_{\alpha}(\mathbf{r}_{{}_{\parallel}},z)=\partial_{\alpha}\Pi(\mathbf{r}_{{}_{\parallel}},z),$
(13)

where $\eta$ is the bulk viscosity of the fluid, together with the
incompressibility constraint ${\bm{\nabla}}_{3}\cdot{\bf v}=0$. Here
${\bm{\nabla}}_{3}\equiv{\hat{\bf x}}\partial/\partial x+{\hat{\bf
y}}\partial/\partial y+{\hat{\bf z}}\partial/\partial z$ is the full three-
dimensional gradient operator, with ${\hat{\bf x}}$, ${\hat{\bf y}}$, and
${\hat{\bf z}}$ as the unit vectors along the $x$, $y$, and $z$ axes
respectively, and $\Pi$ is the bulk pressure which enforces the
incompressibility constraint.

This equation (13) can be solved exactly for the bulk velocity
$\mathbf{v}(\mathbf{r}_{{}_{\parallel}},z,t)$ in terms of the surface velocity
$\mathbf{v}_{s}(\mathbf{r}_{{}_{\parallel}},t)$. If we Fourier expand the
surface velocity:

$\mathbf{v}_{s}(\mathbf{r}_{{}_{\parallel}},t)={1\over\sqrt{L_{x}L_{y}}}\sum_{\mathbf{q}}\mathbf{v}_{s}(\mathbf{q},t)e^{i\mathbf{q}\cdot\mathbf{r}_{{}_{\parallel}}}$
(14)

where $(L_{x},L_{y})$ are the linear dimensions of our (presumed rectangular)
surface, then, as we show in the ALP, the bulk velocity
$\mathbf{v}(\mathbf{r}_{{}_{\parallel}},z,t)$ is given by

$\mathbf{v}(\mathbf{r}_{{}_{\parallel}},z,t)={1\over\sqrt{L_{x}L_{y}}}\sum_{\mathbf{q}}[\mathbf{v}_{s}(\mathbf{q},t)-z(\mathbf{q}\cdot\mathbf{v}_{s})({\hat{\mathbf{q}}}+i{\hat{\mathbf{z}}})]e^{-qz+i\mathbf{q}\cdot\mathbf{r}_{\perp}}\,.$
(15)

The last ingredient in our theory is the boundary condition on the bulk fluid
velocity at the interface. The active particles at the solid-liquid interface
generate active forces, which change the boundary condition from the familiar
partial-slip boundary condition to:

$v_{si}(\mathbf{r}_{{}_{\parallel}},t)=v_{a}(\rho)p_{i}(\mathbf{r}_{{}_{\parallel}},t)+\zeta_{1}(\rho)\hat{{\bf
p}}\cdot\nabla_{S}p_{i}+\zeta_{2}(\rho)p_{i}\nabla_{S}\cdot\hat{{\bf
p}}+p_{i}\hat{{\bf p}}\cdot\nabla_{S}\zeta(\rho)+\mu\eta\bigg{(}\frac{\partial
v_{i}(\mathbf{r}_{{}_{\parallel}},z,t)}{\partial
z}\bigg{)}_{z=0}-\partial_{i}P_{s}(\rho)\,,$ (16)

where $v_{a}(\rho)$ is the spontaneous self-propulsion speed of the active
particles relative to the solid substrate, $\zeta_{1,2}$ and $\zeta$ are
coefficients of the active stresses permitted by symmetry, and $P_{s}(\rho)$
is a surface osmotic pressure. As before $i=(x,y)$. For a system in thermal
equilibrium, $v_{a}=0=\zeta_{1,2}(\rho)=\zeta(\rho)$, and (16) reduces to the
well-known equilibrium partial slip boundary condition partial-slip .

We now turn to the equation of motion for $\hat{{\bf p}}$. As the active
particles are polar, the system lacks $\hat{{\bf p}}\rightarrow-\hat{{\bf p}}$
symmetry. This allows $\partial_{t}\hat{{\bf p}}$ to contain terms even in
$\hat{{\bf p}}$. The most general equation of motion for $p_{k}$ allowed by
symmetry, neglecting “irrelevant” terms,

$\displaystyle\partial_{t}p_{k}=T_{ki}\bigg{(}\alpha
v_{si}-\lambda_{pv}(\mathbf{v}_{s}\cdot\nabla_{s})p_{i}+\left({\nu_{1}-1\over
2}\right)p_{j}\partial_{i}v_{sj}+\left({\nu_{1}+1\over 2}\right)(\hat{{\bf
p}}\cdot\nabla_{s})v_{si}-\lambda(\hat{{\bf
p}}\cdot\nabla_{s})p_{i}-\partial_{i}P_{p}(\rho)+f_{i}\bigg{)},$ (17)

where the projection operator $T_{ki}\equiv\delta^{s}_{ki}-p_{k}p_{i}$ insures
that the fixed length condition $|\hat{{\bf p}}|=1$ on $\hat{{\bf p}}$ is
preserved. It is the breaking of Galilean invariance by the solid substrate
that allows $\lambda_{pv}$ to differ from $1$, and the presence of the “self-
advection” term $\alpha$ in (17). The terms proportional to $\nu_{1}$ are
“flow alignment terms”, identical in form to those found in nematic liquid
crystals martin1972 . The term with coefficient $\lambda$ is allowed by the
polar symmetry of the particles, and can be interpreted as self advection of
the particle polarity in its own direction. The function $P_{p}(\rho)$ is a
density dependent “surface polarization pressure” independent of the “osmotic
pressure” $P_{s}(\rho)$ introduced earlier. We have also added to the equation
of motion (17) a white noise ${\bf f}$ with statistics

$\langle f_{i}({\bf r}_{{}_{\perp}},t)f_{j}({\bf
r}_{{}_{\perp}}^{\prime},t^{\prime})\rangle=2D_{p}\delta_{ij}\delta({\bf
r}_{{}_{\perp}}-{\bf r}_{{}_{\perp}}^{\prime})\delta(t-t^{\prime})\,.$ (18)

Our hydrodynamic model, then, is summarized by the equations of motion (11),
(12), and (17) for $\rho$ and $\hat{{\bf p}}$, respectively, and the solution
(15) of the Stokes equation (13) for the bulk velocity field
$\mathbf{v}(x,y,z,t)$ obtained with the boundary condition (16). Fluctuations
also involve the noise correlations (18) .

These equations of motion and boundary conditions have an obvious spatially
uniform, steady state solution:
$\rho(\mathbf{r}_{{}_{\parallel}},t)=\rho_{0}\,,\hat{{\bf
p}}(\mathbf{r}_{{}_{\parallel}},t)=\hat{{\bf x}}$, where we have defined
$v_{0}\equiv v_{a}(\rho_{0})$ and have chosen the $\hat{{\bf x}}$ axis of our
coordinate system to be along the (spontaneously chosen) direction of
polarization, as illustrated in figure (1).

To study fluctuations about this steady state, we expand the equations of
motion (11), (12), and (17) for $\rho$ and $\hat{{\bf p}}$, and the boundary
condition (16), to linear order in $\delta\rho$ and $p_{y}$. We obtain the
bulk velocity $\mathbf{v}(\mathbf{r}_{{}_{\parallel}},z,t)$ from the surface
velocity $\mathbf{v}_{s}(\mathbf{r}_{{}_{\parallel}},t)$ using our solution
(15) of the Stokes equation. This ultimately produces Eqs. (2) and (3) , where
the phenomenological hydrodynamic parameters $v_{\rho}$, $v_{p}$, $\gamma$,
$\gamma$, $\rho_{c}$, and $\sigma_{t}$ are all related to the expansion
coefficients of the various parameters introduced above when expanded in
powers of the small fluctuations $\delta\rho$ and $p_{y}$. The rather involved
details of this calculation are given in the ALP.

The correlation functions can be straightforwardly determined from these
equations of motion, and shown to have peaks at $\omega_{\rm
peak}=c_{\pm}(\theta_{\mathbf{q}})q$, where $c_{\pm}(\theta_{\mathbf{q}})$ is
given by

$\displaystyle c_{\pm}\left(\theta_{\mathbf{q}}\right)$ $\displaystyle=$
$\displaystyle\pm\sqrt{{1\over
4}\left(v_{\rho}-v_{p}\right)^{2}\cos^{2}\theta_{\mathbf{q}}+c^{2}_{0}\sin^{2}\theta_{\mathbf{q}}}$
(19) $\displaystyle+\left({v_{\rho}+v_{p}\over
2}\right)\cos\theta_{\mathbf{q}}\quad\,.$

We have presented a comprehensive hydrodynamic theory of flocking at a solid-
liquid interface. This theory makes quantitative , experimentally testable
predictions about orientational long range order, spatio-temporal scaling of
fluctuations, giant number fluctuations and anomalous diffusion along
directions parallel to the solid-liquid interface. These predictions are exact
in the asymptotic long wavelength limit, as will be shown in the ALP  using
renormalization group arguments. One simple variant on our system would be to
replace the bulk isotropic fluid of our system with a nematic.

Acknowledgements: One of us (AB) thanks the SERB, DST (India) for partial
financial support through the MATRICS scheme [file no.: MTR/2020/000406]. NS
is partially supported by Netherlands Organization for Scientific Research
(NWO), through the Vidi grant No. 2016/N/00075794. We thank S. Ramaswamy for
sharing reference maitra2018  with us. NS thanks Institut Curie and MPIPKS for
their support through postdoctoral fellowships while some of this work was
being done. AB thanks the MPIPKS, Dresden for their hospitality, and their
support through their Visitors’ Program, while a portion of this work was
underway. JT likewise thanks the MPIPKS for their hospitality, and their
support through the Martin Gutzwiller Fellowship, and the Higgs Center of the
University of Edinburgh for their support with a Higgs Fellowship.

## References

  * (1) K. Kruse, J.F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto, Asters, vortices, and rotating spirals in active gels of polar filaments, Phys. Rev. Lett. 92, 078101(2004).
  * (2) K. Kruse, J.F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto, Generic theory of active polar gels: a paradigm for cytoskeletal dynamics, Eur. Phys. J. E, 16, 5, (2005).
  * (3) H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler, and R.E. Goldstein, Confinement stabilizes a bacterial suspension into a spiral vortex, Phys. Rev. Lett., 110, 268102, (2013).
  * (4) D. Saintillan, and M. J. Shelley, Instabilities, pattern formation, and mixing in active suspensions, Phys. Fluids 20, 123304, (2008).
  * (5) Y. Hatwalne, S. Ramaswamy, M. Rao, S. Madan and A. Simha, Rheology of active-particle suspensions, Phys. Rev. Lett. 92, 118101 (2004).
  * (6) S. Saha, R. Golestanian, and S. Ramaswamy, Clusters, asters, and collective oscillations in chemotactic colloids, Phys. Rev. E, 89, 062316, (2014).
  * (7) B. Liebchen, D. Marenduzzo, I. Pagonabarraga, and M.E. Cates, Clustering and Pattern Formation in Chemorepulsive Active Colloids, Phys. Rev. Lett. 115, 258301(2015).
  * (8) V. Narayan, S. Ramaswamy, and N. Menon, Long-Lived Giant Number Fluctuations in a Swarming Granular Nematic, Science 317, 105, (2007).
  * (9) L.J. Daniels, Y. Park, T. C. Lubensky, and D. J. Durian, Dynamics of gas-fluidized granular rods, Phys. Rev. E79, 041301(2009).
  * (10) A. Baskaran, and M.C. Marchetti, Hydrodynamics of self-propelled hard rods, Phys. Rev. E 77, 011920 (2008).
  * (11) P.G. de Gennes and J. Prost, _The Physics of Liquid Crystals_ (Oxford University Press, Oxford, 1995)
  * (12) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75, 1226, (1995).
  * (13) J. Toner and Y. Tu, Long-Range Order in a Two-Dimensional Dynamical $\mathrm{XY}$ Model: How Birds Fly Together, Phys. Rev. Lett. 75, 4326 (1995).
  * (14) J. Toner and Y. Tu, Flocks, herds, and schools: A quantitative theory of flocking, Phys. Rev. E 58, (1998).
  * (15) J. Toner, Y. Tu, and S. Ramaswamy, Hydrodynamics and phases of flocks, Annals of Physics 318, 170 (2005).
  * (16) R.A. Simha and S. Ramaswamy, Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles, Phys. Rev. Lett. 89, 058101(2002).
  * (17) H. Chaté, F. Ginelli, G. Gregoire and F. Raynaud, Collective motion of self-propelled particles interacting without cohesion, Phys Rev E 77, 046113 (2008); F. Ginelli, The Physics of the Vicsek model, Eur. Phys. J. Special Topics 225, 2099 (2016).
  * (18) J. Toner, Giant number fluctuations in dry active polar fluids: A shocking analogy with lightning rods, J. Chem. Phys. 150, 154120 (2019).
  * (19) S. Ramaswamy, R.A. Simha, and J. Toner, Active nematics on a substrate: Giant number fluctuations and long-time tails, Europhys. Lett. 62, 196 (2003).
  * (20) C. Wolgemuth, E. Hoiczyk, D. Kaiser, and G. Oster, How myxobacteria glide, Current Biology, 12, 369 (2002).
  * (21) E. Lushi, and H. Wioland, and R. Goldstein, Fluid flows created by swimming bacteria drive self-organization in confined suspensions, Proceedings of the National Academy of Sciences 111, 9733 (2014).
  * (22) V. Schaller, and A.R. Bausch, Topological defects and density fluctuations in collectively moving systems, Proceedings of the National Academy of Sciences 110, 4488 (2013).
  * (23) A. Bricard,J.-B. Caussin, N. Desreumaux, O. Dauchot, and D. Bartolo, Emergence of macroscopic directed motion in populations of motile colloids, Nature 503, 95 (2013).
  * (24) D. Geyer, A. Morin and D. Bartolo, Sounds and hydrodynamics of polar active fluids, Nature Materials 17, 789 (2018).
  * (25) N. D. Mermin and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models, Phys. Rev. Lett. 17, 1133 (1966); P. C. Hohenberg, Existence of Long-Range Order in One and Two Dimensions, Phys. Rev. 158, 383 (1967); N. D. Mermin, Absence of Ordering in Certain Classical Systems, J. Math. Phys. 8, 1061 (1967).
  * (26) Exceptions to this result are two dimensional crystals (see, e.g., 2dxtal ) and fluctuating tethered membranes (see, e.g., teth ).
  * (27) B. I. Halperin and D. R. Nelson, Theory of Two-Dimensional Melting, Phys. Rev. Lett. 41, 121 (1978).
  * (28) Y. Kantor, M. Kardar, and D. R. Nelson, Statistical Mechanics of Tethered Surfaces, Phys. Rev. Lett. 57, 791 (1986).
  * (29) S. Shankar, S. Ramaswamy, and M. C. Marchetti, Low-noise phase of a two-dimensional active nematic system, Phys. Rev. E 97, 012707 (2018); S. Mishra, A. Baskaran, and M. C. Marchetti, Fluctuations and pattern formation in self-propelled particles, Phys. Rev. E 81, 061916 (2010).
  * (30) S. Ramaswamy, R. A. Simha, and J. Toner, Active nematics on a substrate: Giant number fluctuations and long-time tails, Europhys. Lett. 62, 196 (2003).
  * (31) N. Sarkar, A. Basu and J. Toner, Associated long paper.
  * (32) These are not actually the sound speeds one would obtain from the real part of the eigenfrequencies $\omega$.
  * (33) Y. Zhu and S. Granick, No-Slip Boundary Condition Switches to Partial Slip When Fluid Contains Surfactant, Langmuir 18, 10058 (2002).
  * (34) P.C. Martin, O. Parodi, and P.S. Pershan, Unified hydrodynamic theory for crystals, liquid crystals, and normal fluids, Phys. Rev. A 6, 2401 (1972).
  * (35) A. Maitra, P. Srivastava, M.C. Marchetti, J.S. Lintuvuori, S. Ramaswamy, and M. Lenz, A nonequilibrium force can stabilize 2D active nematics, Proceedings of the National Academy of Sciences 115, 6934 (2018).


# Learn to Intervene: An Adaptive Learning Policy for Restless Bandits in
Application to Preventive Healthcare

Arpita Biswas1111The research was conducted when A. Biswas worked at Google.
Gaurav Aggarwal2 Pradeep Varakantham2&Milind Tambe2
1Harvard University
2 Google<EMAIL_ADDRESS>{gauravaggarwal, pvarakantham,
<EMAIL_ADDRESS>
###### Abstract

In many public health settings, it is important for patients to adhere to
health programs, such as taking medications and periodic health checks.
Unfortunately, beneficiaries may gradually disengage from such programs, which
is detrimental to their health. A concrete example of gradual disengagement
has been observed by an organization that carries out a free automated call-
based program for spreading preventive care information among pregnant women.
Many women stop picking up calls after being enrolled for a few months. To
avoid such disengagements, it is important to provide timely interventions.
Such interventions are often expensive and can be provided to only a small
fraction of the beneficiaries. We model this scenario as a restless multi-
armed bandit (RMAB) problem, where each beneficiary is assumed to transition
from one state to another depending on the intervention. Moreover, since the
transition probabilities are unknown a priori, we propose a Whittle index
based Q-Learning mechanism and show that it converges to the optimal solution.
Our method improves over existing learning-based methods for RMABs on multiple
benchmarks from literature and also on the maternal healthcare dataset.

##  1 Introduction

Preventive timely intervention (e.g., to check adherence or to provide
medicines) can be used to significantly alleviate many public health issues
such as diabetes Newman et al. (2018), hypertension Brownstein et al. (2007),
tuberculosis Chang et al. (2013); Rahedi Ong’ang’o et al. (2014), HIV
Christopher et al. (2011); Kenya et al. (2013), depression Löwe et al. (2004);
Mundorf et al. (2018), etc. In this paper, we focus on the area of maternal
health (well-being of women during pregnancy, childbirth and the post-natal
period) where preventive interventions can impact millions of women. A key
challenge in such preventive intervention problems is the limited amount of
resources for intervention (e.g., availability of health care workers).
Furthermore, the human behavior (w.r.t. taking medicines or adhering to a
protocol) changes over time and to interventions, thereby requiring a
strategic assignment of the limited resources to the beneficiaries most in
need.

We are specifically motivated towards improving maternal health among low-
income communities in developing countries, where maternal deaths remain
unacceptably high due to not having timely preventive care information
Thaddeus and Maine (1994). We work with a non-profit organization, ARMMAN
2015, that provides a free call-based program for around $2.3$ million
pregnant women. This is similar to other programs, such as Mom-connect
(https://tinyurl.com/mom-connect-jnj). Each enrolled woman receives $140$
automated calls to equip them with critical life-saving healthcare information
across $80$ weeks (pregnancy and up to one year after childbirth).
Unfortunately, the engagement behavior (overall time spent listening to
automated calls) changes and, for most women, the overall engagement
decreases. This can have serious implications on their health. We ask the
question: how do we systematically choose whom to provide interventions
(personal visit by a healthcare worker) in order to maximize the overall
engagement of the beneficiaries?

Preventive intervention problems of interest in this paper are challenging,
owing to multiple key reasons: i) number of interventions are budgeted and
much smaller than the total number of beneficiaries; ii) beneficiary’s
engagement is uncertain and may change after a few months; iii) post-
intervention improvement of beneficiary’s engagement is uncertain; iv)
intervention decisions at a time step have an impact on the state of
beneficiaries and decisions to be taken at the next step. A relevant model for
this setting is restless multi-armed bandits (RMABs). RMABs with prior
knowledge of uncertainty model have been studied for health interventions Lee
et al. (2019); Mate et al. (2020, 2021); Killian et al. (2021); Bhattacharya
(2018), sensor monitoring tasks Iannello et al. (2012); Glazebrook et al.
(2006), anti-poaching patrols Qian et al. (2016), and uplift modeling in
eCommerce platforms Gubela et al. (2019). Due to the unpredictability of human
beneficiaries, it is unrealistic to know the uncertainty model a priori. Our
key contribution is in designing an intervention scheme under limited budget,
without relying on a priori knowledge of the uncertainty model.

##### Contributions.

First, we represent the preventive intervention problem as a Restless Multi-
Arm Bandit (RMAB), where the uncertainty model associated with beneficiary
behaviors with and without intervention are not known a priori, thereby
precluding the direct application of Whittle Index based methods Whittle
(1988). Second, we develop a model-free learning method based on Q-learning,
referred to as Whittle Index based Q-Learning (WIQL) that executes actions
based on the difference between Q-values of active (intervention) and passive
actions. We show that WIQL converges to the optimal solution asymptotically.
Finally, to show that WIQL is a general approach and applicable to multiple
domains, we evaluate WIQL on various examples. We then simulate the
intervention problem in the context of maternal healthcare. We demonstrate
that the model assumptions we make capture the problem of concern and the
intervention scheme employed by WIQL significantly improves the engagement
among beneficiaries compared to the existing benchmarks.

##  2 Related Work

The restless multi-armed bandit (RMAB) problem was introduced by Whittle 1988.
The main result involves formulating a relaxation of the problem and solving
it optimally using a heuristic called Whittle Index policy. This policy is
optimal when the underlying Markov Decision Processes satisfy indexability,
which is computationally intensive to verify. Moreover, Papadimitriou and
Tsitsiklis 1994 established that solving RMAB is PSPACE-hard.

There are three main threads of relevant research. The first category focuses
on specific classes of RMABs. Akbarzade et al. 2019 provide a class of bandits
with “controlled restarts” and state-independent policies, which possess the
indexability property and a Whittle index policy. Mate et al. 2020 consider
two-state bandits to model a health intervention problem, where taking an
intervention collapses the uncertainty about their current state. Bhattacharya
2018 models the problem of maximizing the coverage and spread of health
information as an RMAB problem and proposes a hierarchical policy. Lee et al.
2019 study the problem of screening patients to maximize early-stage cancer
detection under limited resource, by formulating it as a subclass of RMAB.
Similarly, Glazebrook et al. 2006, Hsu 2018, Sombabu et al. 2020, Liu and Zhao
2010 give Whittle indexability results for different subclasses of (hidden)
Markov bandits. This category of research assumes that the transition and
observation models are known beforehand. Instead, our focus is on providing
learning methods when the transition model is unknown a priori.

The second category contains different learning methods for RMABs. Fu et al.
2019 provide a Q-learning method where the Q value is defined based on the
Whittle indices, states, and actions. However, they do not provide proof of
convergence to optimal solution and experimentally, do not learn
(near-)optimal policies. Along similar lines, Avrachenkov and Borkar 2020
provide a fundamental change to the Q-value definition for computing optimal
whittle index policy. However, their convergence proof requires all
homogeneous arms with same underlying MDPs. We provide a learning method that
is not only shown to theoretically converge but also empirically outperforms
the above mentioned methods on benchmark instances and also on a real
(heterogeneous arms) problem setting.

The third relevant line of work is to predict adherence to a health program
and effects of interventions, by formulating these as supervised learning
problems. Killian et al. 2019 use the data from 99DOTS Cross et al. (2019) and
train a deep learning model to target patients at high risk of not adhering to
the health program. On similar lines, Nisthala et al. 2020 use the engagement
data of a health program and train deep learning model to predict patients who
are at high risk of dropping out of the program. Son et al. 2010 use Support
Vector Machine to predict adherence to medication among heart failure
patients. There are many other papers Howes et al. (2012); Lauffenburger et
al. (2018) that train models on historical data for predicting adherence and
provide interventions to patients who have low adherence probability. These
works assume the training data to be available beforehand. In contrast, we
consider the online nature of the problem where feedback is received after an
intervention is provided, which in turn is used in making future decisions.

##  3 Preliminaries

RMAB models various stochastic scheduling problems, where an instance is a
3-tuple ($N$, $M$, $\\{MDP_{i}\\}_{i\in N}$) with $N$ denoting the set of
arms, $M$ is the budget restriction denoting how many arms can be pulled at a
given time, and MDPi is the Markov Decision Process for each arm $i$. An MDP
consists of a set of states $\mathcal{S}$, a set of actions $\mathcal{A}$,
transition probabilities
$\mathcal{P}:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\mapsto[0,1]$, and
reward function $R:\mathcal{S}\times\mathcal{A}\mapsto\mathbb{R}$. The action
set $\mathcal{A}$ of each MDP consists of two actions: an active action ($1$)
and a passive action ($0$). At each time step $t$, an action $A_{i}(t)$ is
taken on an arm $i$, such that $\sum_{i}A_{i}(t)=M$. Then, each arm $i$
transitions to a new state and observes a reward, according to the underlying
MDPi. Let $X_{i}(t)\in\mathcal{S}$ and $R_{i}^{X_{i}(t)}(A_{i}(t))$ denote the
current state and reward obtained at time $t$ respectively. Now, policy per
step $\pi:X_{1}(t)\times\ldots\times X_{N}(t)\mapsto\\{A_{i}(t)\\}_{i\in N}$
can be defined as a mapping from the current states of all beneficiaries to
the actions to be taken on each arm. Thus, given a policy $\pi$, the action on
an arm $i$ is denoted as:

$A_{i}^{\pi}(t)=\begin{cases}1&\quad\text{if }i\text{ is selected by policy
}\pi\text{ at time }t\\\ 0&\quad\text{if }i\text{ is not selected by policy
}\pi\text{ at time }t\end{cases}$

The goal is to find a policy $\pi^{*}$ that maximizes the expected reward
until $T$ time steps, subject to the budget.

$\displaystyle\underset{\pi}{\max}$ $\displaystyle\frac{1}{T}\
\mathbb{E}\left[\sum_{i\in
N}\sum_{t=1}^{T}R_{i}^{X_{i}(t)}(A^{\pi}_{i}(t))\right]$ (1) s.t.
$\displaystyle\sum_{i\in N}A^{\pi}_{i}(t)=M\qquad\mbox{ for all
}t=\\{1,\ldots,T\\}$

To deal with the computational hardness of solving this problem, an index-
based heuristic policy based on the Lagrangian relaxation of the RMAB problem
was proposed by Whittle 1988—at each time step $t$, an index is computed for
each arm depending on the current state of the arm, transition probabilities,
and reward function of its MDP. Then, the top $M$ arms with highest index
values are selected.

##### Whittle Index-based policy.

Whittle’s relaxation is to replace the budget constraint of $M$ on the number
of arms to a time-averaged constraint, i.e.,

$\frac{1}{T}\ \mathbb{E}\left[\displaystyle\sum_{t=1}^{T}\sum_{i\in
N}A^{\pi}_{i}(t)\right]=M$

Further, using Lagrange’s relaxation (with $\lambda$ as the Lagrange’s
multipliers) and dropping other constants, the objective function can be
rewritten as:

$\underset{\pi}{\max}\
\frac{1}{T}\mathbb{E}\left[\displaystyle\sum_{t=1}^{T}\sum_{i\in
N}\\!\left(\\!R_{i}^{X_{i}(t)}(A^{\pi}_{i}(t))+\lambda\cdot(1\\!\\!-\\!\\!A^{\pi}_{i}(t))\right)\right]$
(2)

Whittle showed that this problem can be decoupled and solved for each arm by
computing the index $\lambda_{i}(Z)$ which acts like a subsidy that needs to
be given to an arm $i$ at state $Z$, so that taking the action $1$ is as
beneficial as taking the action $0$. Assuming indexability, choosing arms with
higher subsidy $\lambda_{i}(Z)$ leads to the optimal solution for Equation 2.

######  Definition 1 (Indexability)

Let $\Phi(\lambda)$ be the set of states for which it is optimal to take
action $0$ when taking an active action costs $\lambda$. An arm is indexable
if $\Phi(\lambda)$ monotonically increases from $\emptyset$ to $\mathcal{S}$
when $\lambda$ increases from $-\infty$ to $+\infty$. An RMAB problem is
indexable if all the arms are indexable.

##  4 The Model

We formulate the problem of selecting $M$ out of $N$ state-transitioning arms
at each time step, as an RMAB, with beneficiaries being the arms in preventive
healthcare intervention scenarios. We represent the engagement pattern as an
MDP (Figure 1) with three (abstract) “behavioral” state: (i) self motivated
($S$): in this state, beneficiary shows high engagement and there is no need
to intervene in this state. (ii) persuadable ($P$): in this state, the
beneficiary engages less frequently with a possibility of increasing
engagement when intervened, which makes this the most interesting state, and
(iii) lost Cause ($L$): The engagement is very low in this state and very
likely to remain low irrespective of intervention.

These three states capture different levels of engagement as well as
differences in terms of the benefit obtained by interventions, which is not
possible to represent using other existing intervention models, such as the
model described by Mate et al. 2020. Note that the more the states, slower is
the convergence of any online algorithm. Thus, for short-term intervention
programs (a year or two), we recommend a three-states model. In Section 7.2,
we provide a mechanism to obtain the states from real-data. Let
$A_{i}(t)\in\\{0,1\\}$ be the action taken on beneficiary $i$ at time $t$; $1$
denotes an intervention and $0$ otherwise. Also, $\sum_{i}A_{i}(t)=M$ for each
time slot $t$. Depending on the action taken, each beneficiary changes its
states according to the transition probabilities.

Figure 1: MDP for modeling beneficiary $i$’s engagement. Blue and red figures
represent transitions $p_{i,a}^{YZ}$ from a state $Y$ to $Z$, with
intervention ($a\\!=\\!1$) and without ($a\\!=\\!0$), respectively.

Though our model is inspired by a preventive intervention in healthcare, this
formulation captures the intervention problem in many other domains, such as
sensor monitoring, anti-poaching patrols, and uplift modeling in eCommerce.
Depending on the domain, a persuadable state is where the benefit of an
intervention is maximum. When a beneficiary is in the persuadable state and
intervention is provided, an arm is more likely to transition to the self-
motivated state whereas, when no intervention is provided, it is more likely
to transition to the lost-cause state. These transitions are denoted as
$p^{PS}_{i,1}$ and $p^{PL}_{i,0}$ for action $1$ and $0$, respectively. Every
transition generates a state-dependent reward. Since better engagements have
higher rewards, we assume $R_{i}^{S}>R_{i}^{P}>R_{i}^{L}$. The transition
probabilities $p_{i,a}^{YZ}$ may be different for each beneficiary and thus,
the effect of intervention would vary from one beneficiary to another. If the
transition probabilities were known beforehand, one could compute the Whittle-
index policy and select accordingly. However, here, transition probabilities
are unknown, and hence, we consider the problem of learning the Whittle index
while simultaneously selecting a set of best arms depending on the estimated
Whittle Index.

##  5 Whittle Index based Q-learning (WIQL)

Q-Learning Watkins and Dayan (1992) is a well-studied reinforcement learning
algorithm for estimating $Q^{*}(Z,a)$ for each state-action pair ($Z,a$) of an
MDP.

$Q^{*}(Z,a):=R^{Z}+\displaystyle\sum_{Z^{\prime}\in\\{S,P,L\\}}\mathcal{P}(Z,a,Z^{\prime})\cdot
V^{*}(Z^{\prime}),$

where the optimal expected value of a state is given by

$V^{*}(Z):=\underset{a\in\\{0,1\\}}{\max}\left(R^{Z}+\displaystyle\sum_{Z^{\prime}\in\\{S,P,L\\}}\mathcal{P}(Z,a,Z^{\prime})\cdot
V^{*}(Z^{\prime})\right)$

Q-Learning estimates $Q^{*}$ using point samples—at each time $t$, an agent
(policy maker) takes an action $a$ using estimated $Q$ values at the current
state $Z$, a reward $R$ is observed, a new state $Z$ is reached, and $Q$
values are updated according to the following update rule:

$\displaystyle Q^{t+1}(Z,a)$ $\displaystyle\leftarrow$
$\displaystyle(1-\alpha_{t}(Z,a))\cdot Q^{t}(Z,a)+$ (3)
$\displaystyle\alpha_{t}(Z,a)\\!\cdot\\!\left(R^{Z}(t)+\underset{a^{\prime}\in\\{0,1\\}}{\max}Q^{t}(Z^{\prime},a^{\prime})\right)\
$

Here, $\alpha_{t}(\cdot)\in[0,1]$ is the learning parameter. When $\alpha=0$,
the agent does not learn anything new and retains the value obtained at
$t^{th}$ step, while $\alpha=1$ stores only the most recent information and
overwrites all previously obtained rewards. Setting $0<\alpha<1$ strikes a
balance between the new values and the old ones. With mild assumptions on
$\alpha$, the convergence of Q-Learning to the optimal $Q^{*}$ values has been
established Watkins and Dayan (1992); Jaakkola et al. (1994); Borkar and Meyn
(2000). We build on the results to provide a Q-Learning approach for learning
Whittle-index policy.

We adopt Q-Learning for RMABs and store Q values separately for each arm. Q
values of state-action pairs are typically used for selecting the best action
for an arm at each state; however, for RMAB, the problem is to select $M$
arms. Our action-selection method ensures that an arm whose estimated benefit
is higher is more likely to get selected. Algorithm 1 describes the action-
selection and the update rule. This algorithm is not specific to any
particular state representation and can be used for any finite state RMAB
instance. However, when the number of states are large (10 states), the
convergence is typically slow and not suitable for short horizon.

Input: $N$, $M$, $\alpha(c)$, and initial states $X_{i}(0)\in\mathcal{S}\
\forall\ i\in[N]$.

1 Initialize: $Q_{i}^{0}(Z,a)\leftarrow 0$ and $\lambda_{i}^{0}(Z)\leftarrow
0$ for each state $Z\in\mathcal{S}_{i}$, each action $a\in\\{0,1\\}$, each arm
$i\in[N]$.

2for  _$t$ in $1,\ldots,T$_ do

// Select $M$ arms using $\epsilon$-decay policy

3 $\epsilon\leftarrow\frac{N}{N+t}$.

4 With probability $\epsilon$, select $M$ arms uniformly at random. Otherwise,
select top $M$ arms according to their values $\lambda_{i}^{t}(X_{i}(t))$.
Store the selected arms in $\Psi$.

// Take suitable actions on the arms

5 for  _ $i$ in $1:N$_ do

6 if  _$i\in\Psi$_ then

7 Take action $A_{i}(t)=1$ on arm $i$.

8 else

9 Take action $A_{i}(t)=0$ on arm $i$.

10

11 Observe reward $r$ and next state $X_{i}(t+1).$

12

// Update Q, $\lambda$ using $Z=X_{i}(t)$, $a=A_{i}(t)$

13 $c^{h}_{i,Z,a}=\displaystyle\sum_{h=0}^{t}\mathbb{I}\\{X_{i}(h)=Z\mbox{ \&
}A_{i}(h)=a\\}$.

14 for  _ $i$ in $1:N$_ do

15 $Q_{i}^{t+1}(Z,a)\leftarrow(1-\alpha(c^{h}_{i,Z,a}))\cdot
Q_{i}^{t}(Z,a)+\alpha(c^{h}_{i,Z,a})\cdot\left(r+\underset{a^{\prime}\in\\{0,1\\}}{\max}Q_{i}^{t}(X_{i}(t+1),a^{\prime})\right)$.

16 $\lambda_{i}^{t+1}(X_{i}(t)=Z)\leftarrow Q_{i}^{t+1}(Z,1)-Q_{i}^{t+1}(Z,0)$

17

18

Algorithm 1 Whittle Index based Q-learning(WIQL)

We propose Whittle Index based Q-learning (WIQL), that uses an
$\epsilon$-decay policy to select $M$ arms at each time step $t$. During early
steps, arms are likely to be selected uniformly at random. As time proceeds,
arms with higher values of their estimated $\lambda_{i}(X_{i}(t))$ gets more
priority. The selected set of $M$ arms (who receive interventions) is called
$\Psi$. Each arm is restless, i.e., each arm transitions to a new state and
observes a reward, with or without interventions. These observations are then
used for updating the Q values in Step $13$ (Eq. 3). While updating
$Q_{i}(Z,a)$, we use a learning parameter $\alpha(c_{i,Z,a})$ that decreases
with increase in $c^{t}_{i,Z,a}$ (number of times the arm $i$ received action
$a$ at state $Z$); eg, $\alpha(c^{t}_{i,Z,a})=1/(c^{t}_{i,Z,a}+1)$ satisfies
this criteria. These Q values are then used to estimate the Whittle index
$\lambda(X_{i}(t))$.

##  6 Theoretical Results

In this section, we show that the WIQL does not alter the optimality
guarantees of Q-learning. First, we show that taking intervention decisions
based on the benefit, (i.e., difference in Q values of active and passive
actions) is equivalent to optimizing joint Q over all arms subject to the
budget on intervention actions.

###### Theorem 1

Taking action $1$ on the top $M$ arms according to
$(Q^{*}_{i}(s_{i},1)-Q^{*}_{i}(s_{i},0))$ is equivalent to maximizing
$\sum_{i}Q_{i}^{*}(s_{i},a_{i})$ over all possible action profiles
$(a_{1},\ldots,a_{N})\in\\{0,1\\}^{N}$ such that $\sum_{i}a_{i}=M$.

Proof Sketch. For ease of explanation, we prove this for M=1. Let $i^{*}$ be
the arm that maximizes the benefit of taking an intervention action
($Q^{*}_{i^{*}}(s_{i^{*}},1)-Q^{*}_{i^{*}}(s_{i^{*}},0)$) at its current state
$s_{i^{*}}$. Then, for any $j\in\\{1,\ldots,N\\}\setminus\\{i^{*}\\}$.

$\displaystyle Q^{*}_{i^{*}}(s_{i^{*}},1)-Q^{*}_{i^{*}}(s_{i^{*}},0)\geq
Q^{*}_{j}(s_{j},1)-Q^{*}_{j}(s_{j},0)$ $\displaystyle
Q^{*}_{i^{*}}(s_{i^{*}},1)+Q^{*}_{j}(s_{j},0)\geq
Q^{*}_{j}(s_{j},1)+Q^{*}_{i^{*}}(s_{i^{*}},0)$ (4) Adding $\sum_{i\neq
i*,i\neq j}Q_{i}^{*}(s_{i},0)$ on both sides $\displaystyle
Q^{*}_{i^{*}}(s_{i^{*}},1)\\!\\!+\\!\\!\displaystyle\sum_{i\neq
i^{*}}Q^{*}_{i}(s_{i}0)\geq
Q^{*}_{j}(s_{j},1)\\!\\!+\\!\\!\displaystyle\sum_{i\neq j}Q^{*}_{i}(s_{i}0).$
(5)

Eq. (5) shows that taking intervention action on $i^{*}$ and passive actions
on all other arms would maximize $\sum_{i}Q^{*}(s_{i},a_{i})$ when $M$=$1$.
This argument holds true for any $M\geq 1$ (complete proof is in the
Appendix).$\hfill\square$

###### Theorem 2

WIQL converges to the optimal with probability $1$ when
$\sum_{h}\alpha(c^{h}_{i,Z,a})=\infty\mbox{ and
}\sum_{h}\alpha(c^{h}_{i,Z,a})^{2}<\infty.$

Proof Sketch. This proof follows from (1) the convergence guarantee of
Q-Learning algorithm Watkins and Dayan (1992), (2) $\epsilon$-decay selection
policy, and (3) theorem 1. It has been established in Watkins and Dayan (1992)
that the update rule of Q-Learning converges to $Q^{*}$ whenever
$\sum_{t}\alpha_{t}(Z,a)=\infty$ and $\sum_{t}\alpha_{t}(Z,a)^{2}<\infty.$
These assumptions require that all state-action pairs be visited infinitely
often, which is guaranteed by the $\epsilon$-decay selection process, where
each arm has a non-zero probability of being selected uniformly at random.
Thus, $Q^{t}_{i}(Z,a)$ converges to $Q^{*}_{i}(Z,a)$ which implies that
$Q^{t}_{i}(Z,1)-Q^{t}_{i}(Z,0)$ converges to $Q^{*}_{i}(Z,1)-Q^{*}_{i}(Z,0)$
(using series convergence operation). Finally, using Theorem 1, we claim that
selecting arms based on highest values of $Q^{*}_{i}(Z,1)-Q^{*}_{i}(Z,0)$
would lead to an optimal solution problem. This completes the proof that WIQL
converges to the optimal. $\hfill\square$

##  7 Experimental Evaluation

We compare WIQL against five benchmarks: (1) OPT: assumes full knowledge of
the underlying transition probabilities and has access to the optimal Whittle
Indices, (2) AB Avrachenkov and Borkar (2020), (3) Fu Fu et al. (2019), (4)
Greedy: greedily chooses the top $M$ arms with the highest difference in their
observed average rewards between actions $1$ and $0$ at their current states,
and (5) Random: chooses $M$ arms uniformly at random at each step.

We consider a numerical example and a maternal healthcare application to
simulate RMAB instances using beneficiaries’ behavioral pattern from the call-
based program. For each experiment, we plot the total reward averaged over
$30$ trials, to reduce the effect of randomness in the action-selection
policy.

###  7.1 Numerical Example: Circulant Dynamics

This example has been studied in the existing literature on learning Whittle
Indices Avrachenkov and Borkar (2020); Fu et al. (2019). Each arm has four
states $\mathcal{S}=\\{0,1,2,3\\}$ two actions $\mathcal{A}=\\{0,1\\}$. The
rewards are $R^{0}=-1$, $R^{1}=R^{2}=0$, and $R^{3}=1$ for $a\in\\{0,1\\}$.
The transition probabilities for each action $a\in\\{0,1\\}$ are represented
as a $|\mathcal{S}|\times|\mathcal{S}|$ matrix:

$\mathcal{P}^{1}\\!\\!=\\!\\!\begin{pmatrix}0.5&0.5&0&0\\\ 0&0.5&0.5&0\\\
0&0&0.5&0.5\\\
0.5&0&0&0.5\end{pmatrix},\mathcal{P}^{0}\\!\\!=\\!\\!\begin{pmatrix}0.5&0&0&0.5\\\
0.5&0.5&0&0\\\ 0&0.5&0.5&0\\\ 0&0&0.5&0.5\end{pmatrix}$

The optimal Whittle Indices for each state are as follows:
$\lambda^{*}(0)=-0.5$, $\lambda^{*}(1)=0.5$, $\lambda^{*}(2)=1$, and
$\lambda^{*}(3)=-1$. Whittle Index policy would prefer taking action $1$ on
arms who are currently at state $2$, followed by those at state $1$, then
those at state $0$, and lastly those at state $3$.

Figure 2 demonstrates that WIQL gradually increases towards the OPT policy,
for two sets of experiments—(1) $N$=$5$ and $M$=$1$ and (2) $N$=$100$ and
$M$=$20$. We observe that AB converges towards an average reward of zero. This
happens because it prioritizes arms who are currently at state $1$ over other
states. Since, the expected reward of taking action $1$ at state $1$ is
$\mathcal{P}(1,1,1)\cdot R^{1}+\mathcal{P}(1,1,2)\cdot R^{2}=0$, the total
average reward also tends to zero. The result obtained by the algorithm Fu is
the same as what is shown in Figure $2$ of their paper Fu et al. (2019) where
the total average reward converges to a value of $0.08$. As expected, Greedy
and Random, one being too myopic and the other being too exploratory in
nature, are unable to converge to the optimal value.

These observations show that WIQL outperforms the existing algorithms on the
example that was considered in the earlier papers. Note that, while
implementing the algorithms AB and Fu, we fixed the hyperparameters to the
values specified for this example. However for the real-world application,
that we consider next, it is not obvious how to obtain the best set of
hyperparameters for their algorithms. Thus, we do not compare these algorithms
for the maternal healthcare application. Next, we compare the performance of
WIQL algorithm with Greedy, Random and a Myopic policy (defined in the
subsequent paragraph).

(a) Average reward: $1$ out of $5$ arms chosen.

(b) Average reward: $20$ out of $100$ arms chosen.

Figure 2: Circulant Dynamics: Results averaged over $30$ trials.

###  7.2 Real-world Application: Maternal Healthcare

Figure 3: Evaluating WIQL on maternal healthcare application with static and
dynamic behavior change (shown in first and second row respectively). The
graphs show the total reward accumulated on a particular week, averaged over
30 iterations.

We now focus on the maternal healthcare intervention problem where only a
small subset of beneficiaries can be selected for providing interventions
every week. We use the data, obtained from the call-based preventive care
program, that contains call-records of enrolled beneficiaries—how long they
listened to the call, whether an intervention was given, and when. The data
also contain the ID of a healthcare worker dedicated to provide personalized
intervention to each enrolled beneficiary. The data was collected for an
experimental study towards building up a robust intervention program (which is
the focus of this work). During the experimental study, only one intervention
(in-person visit by a health worker) was provided to a selected set of
beneficiaries ($1559$ chosen from among $3031$ beneficiaries) who were more
likely to drop-out of the program. We call this the Myopic intervention and
use it as a benchmark to compare our approach.

We use the three-state MDP model (Figure 1) to simulate the behavior of the
beneficiaries. During a particular week, the beneficiaries listening to
$>50\%$ of the content of the automated calls are in state $S$ (self-
motivated), those listening to $5-50\%$ are in state $P$ (persuadable), and
those listening to $<5\%$ are in state $L$ (lost cause). A reward of $2$ is
obtained when a beneficiary is in state $S$, a reward of $1$ is obtained in
state $P$, and a reward of $0$ is obtained in state $L$. Thus, a high total
reward accumulated per week implies that a large number of beneficiaries are
at either state $S$ or $P$.

Some other observations on the call-records are as follows. The beneficiaries
who were in state $P$ on a particular week never transitioned to state $S$ in
the immediate next week, unless they received an intervention. On the other
hand, a few beneficiaries who were in state $L$ transitioned to state $P$ even
without intervention. Moreover, the fraction of times a transition from state
$L$ to $P$ occurred is almost the same with and without the intervention at
state $L$. Even though these transition probabilities at the level of all
users are known based on this data, it is difficult to know the actual
transitions for a given beneficiary a priori and hence, we simulate this
behavior. We conduct two set of simulations, namely static (transition
probabilities of each beneficiary remain same throughout their enrollment),
and dynamic (transition probabilities change a few weeks after their
enrollment).

##### Static State-Transition Model.

We assume three categories of arms—(A) high chance of improvement: highly
likely to improve their engagement on receiving an intervention, and
deteriorate in absence of an intervention when they are at state $P$, i.e.,
$p^{PS}=0.8$ and $p^{PL}=0.8$, (B) medium chance of improvement: $p^{PS}=0.4$
and $p^{PL}=0.6$, and (C) low change improvement: $p^{PS}=0.1$ and
$p^{PL}=0.6$. We assume $10$ arms belong to category-A, $10$ arms belong to
category-B and $30$ arms belong to category-C. This assumption helps us
determine the efficacy of any learning algorithm; in particular, the most
efficient algorithm would quickly learn to intervene the $10$ arms of
category-A whenever they are at state $P$. We compare WIQL with greedy,
random, and Myopic algorithm.

##### Dynamic State-Transition Model.

For the dynamic setting, we simulate the first $28$ weeks as described in
Section 7.2.1. Further, we assume that the arms in category (A) change to
“medium” chance of improvement after they are enrolled for around $28$ weeks,
those in category (B) change to “low” chance of improvement. Also, there is a
new set of $10$ arms with “high” chance of improvement. Note that, in reality,
the behavior change would happen at any arbitrary time; however, to check if
WIQL adapts to dynamic change of transition probabilities, we set this
simulation environment.

We run each simulation for $80$ weeks (the length of the program). In Figure 3
we provide results obtained by considering various values of
$M\in\\{500,1000,1500\\}$, where a value of $M$ represents the total number of
personalized visits made by the $100$ health-care workers on a week. We
observe that, for $M=500$ (each health-care worker visits only $5$
beneficiaries per week), WIQL performs (only) marginally better than the other
benchmarks. In contrast, when $M\in\\{1000,1500\\}$, the reward obtained by
WIQL is higher than Greedy and significantly more than Myopic and Random.
Comparing WIQL and Greedy based on their total reward per week, we observe
that the per-week engagement outperforms Greedy by a significant margin.
Observe that, the convergence of WIQL to a total reward of $5000$ is quicker
when $M$ is higher. This is because more sample points are observed per week.
Additionally, we observe that the myopic algorithm leads to overall low
engagement among the beneficiaries, even under the static setting. These
results show that WIQL is able to learn which arms should be intervened at
which state without any prior knowledge about the transition probabilities,
and also adapts better to the underlying transition dynamics.

##  8 Conclusion and Discussion

We focus on a limited-resource sequential decision problem and formulate it as
an RMAB setting. We provide a mechanism to systematically learn and decide on
whom to intervene, and hence, improve the overall benefit of intervention. Our
method possesses the capacity of impacting and improving the well-being of
millions, for example, in the maternal healthcare domain, as demonstrated in
this paper.

WIQL is a general solution for learning RMABs, and we demonstrate this using
examples from other domains, such as Circulant Dynamics. Additionally, WIQL
can be used in a more general setting where new arms arrive over time. For the
departing arms, we can assume that each arm ends in a new state “dropped” and
never transitions to the earlier states. Moreover, our experiments show that
WIQL adapts to the dynamic behavioral change. In practice, however, WIQL may
not be directly applicable for domains where there are beneficiaries with
extreme health risk. The non-zero probability of being selected randomly may
come at a cost of a patient in critical need of intervention. One way to
mitigate this issue is to assign priorities to beneficiaries depending on
their comorbidities, possible complication during pregnancy and after
childbirth, etc. If the high-risk class is small, we can target the
intervention for all of them, and run WIQL on the remaining beneficiaries.
This constraint may hamper the convergence guarantee of WIQL; however, it
would benefit the enrolled women at large.

Going ahead, it would be interesting to solve the RMAB learning problem
considering a large number of states and more than two actions, each with a
different cost of operation.

## References

  * Akbarzadeh and Mahajan [2019] N. Akbarzadeh and A. Mahajan.  Restless bandits with controlled restarts: Indexability and computation of whittle index.  In 2019 IEEE Conference on Decision and Control. IEEE, 2019.
  * ARMMAN [2015] ARMMAN.  Title.  Leveraging technology to create scalable solutions empowering mothers and enabling healthy children, 2015.  Accessed: 2021-05-18.
  * Avrachenkov and Borkar [2020] Konstantin Avrachenkov and Vivek S Borkar.  Whittle index based q-learning for restless bandits with average reward.  arXiv preprint arXiv:2004.14427, 2020.
  * Bhattacharya [2018] Biswarup Bhattacharya.  Restless bandits visiting villages: A preliminary study on distributing public health services.  In Proceedings of the 1st ACM SIGCAS Conference on Computing and Sustainable Societies, pages 1–8, 2018.
  * Borkar and Meyn [2000] Vivek S Borkar and Sean P Meyn.  The ode method for convergence of stochastic approximation and reinforcement learning.  SIAM Journal on Control and Optimization, 38(2):447–469, 2000.
  * Brownstein et al. [2007] J Nell Brownstein, Farah M Chowdhury, Susan L Norris, Tanya Horsley, Leonard Jack Jr, Xuanping Zhang, and Dawn Satterfield.  Effectiveness of community health workers in the care of people with hypertension.  American journal of preventive medicine, 32(5):435–447, 2007.
  * Chang et al. [2013] Alicia H Chang, Andrea Polesky, and Gulshan Bhatia.  House calls by community health workers and public health nurses to improve adherence to isoniazid monotherapy for latent tuberculosis infection: a retrospective study.  BMC public health, 13(1):894, 2013.
  * Christopher et al. [2011] Jason B Christopher, Alex Le May, Simon Lewin, and David A Ross.  Thirty years after alma-ata: a systematic review of the impact of community health workers delivering curative interventions against malaria, pneumonia and diarrhoea on child mortality and morbidity in sub-saharan africa.  Human resources for health, 9(1):27, 2011.
  * Cross et al. [2019] Andrew Cross, Nakull Gupta, Brandon Liu, Vineet Nair, Abhishek Kumar, Reena Kuttan, Priyanka Ivatury, Amy Chen, Kshama Lakshman, Rashmi Rodrigues, et al.  99dots: a low-cost approach to monitoring and improving medication adherence.  In Proceedings of the Tenth International Conference on Information and Communication Technologies and Development, pages 1–12, 2019\.
  * Fu et al. [2019] Jing Fu, Yoni Nazarathy, Sarat Moka, and Peter G Taylor.  Towards q-learning the whittle index for restless bandits.  In 2019 Australian & New Zealand Control Conference (ANZCC), pages 249–254. IEEE, 2019.
  * Glazebrook et al. [2006] K.D. Glazebrook, D. Ruiz-Hernandex, and C. Kirkbride.  Some indexable families of restless bandit problems.  Adv. Appl. Probab, pages 643–672, 2006.
  * Gubela et al. [2019] Robin Gubela, Artem Bequé, Stefan Lessmann, and Fabian Gebert.  Conversion uplift in e-commerce: A systematic benchmark of modeling strategies.  International Journal of Information Technology & Decision Making, 18(03):747–791, 2019.
  * Howes et al. [2012] Christine Howes, Matthew Purver, Rose McCabe, Patrick Healey, and Mary Lavelle.  Predicting adherence to treatment for schizophrenia from dialogue transcripts.  In Proceedings of the 13th Annual Meeting of the Special Interest Group on Discourse and Dialogue, pages 79–83, 2012.
  * Hsu [2018] Y. Hsu.  Age of information: Whittle index for scheduling stochastic arrivals.  In 2018 IEEE International Symposium on Information Theory. IEEE, 2018.
  * Iannello et al. [2012] F. Iannello, O. Simeone, and U. Spagnolini.  Optimality of myopic scheduling and whittle indexability for energy harvesting sensors.  In Conference on Information Sciences and Systems (CISS). IEEE, 2012\.
  * Jaakkola et al. [1994] Tommi Jaakkola, Michael I Jordan, and Satinder P Singh.  Convergence of stochastic iterative dynamic programming algorithms.  In Advances in neural information processing systems, pages 703–710, 1994.
  * Kenya et al. [2013] Sonjia Kenya, Jamal Jones, Kristopher Arheart, Erin Kobetz, Natasha Chida, Shelly Baer, Alexis Powell, Stephen Symes, Tai Hunte, Anne Monroe, et al.  Using community health workers to improve clinical outcomes among people living with hiv: a randomized controlled trial.  AIDS and Behavior, 17(9):2927–2934, 2013.
  * Killian et al. [2019] J. A. Killian, B. Wilder, A. Sharma, V. Choudhary, B. Dilkina, and M. Tambe.  Learning to prescribe interventions for tuberculosis patients using digital adherence data.  In KDD. ACM, 2019.
  * Killian et al. [2021] Jackson A Killian, Andrew Perrault, and Milind Tambe.  Beyond” to act or not to act”: Fast lagrangian approaches to general multi-action restless bandits.  In Proceedings of the 20th International Conference on Autonomous Agents and MultiAgent Systems, pages 710–718, 2021.
  * Lauffenburger et al. [2018] Julie C Lauffenburger, Jessica M Franklin, Alexis A Krumme, William H Shrank, Olga S Matlin, Claire M Spettell, Gregory Brill, and Niteesh K Choudhry.  Predicting adherence to chronic disease medications in patients with long-term initial medication fills using indicators of clinical events and health behaviors.  Journal of managed care & specialty pharmacy, 24(5):469–477, 2018\.
  * Lee et al. [2019] Elliot Lee, Mariel S Lavieri, and Michael Volk.  Optimal screening for hepatocellular carcinoma: A restless bandit model.  Manufacturing & Service Operations Management, 21(1):198–212, 2019\.
  * Liu and Zhao [2010] K. Liu and Q. Zhao.  Indexability of restless bandit problems and optimality of whittle index for dynamic multichannel access.  IEEE Transactions on Information Theory, pages 5547–5567, 2010\.
  * Löwe et al. [2004] Bernd Löwe, Jürgen Unützer, Christopher M Callahan, Anthony J Perkins, and Kurt Kroenke.  Monitoring depression treatment outcomes with the patient health questionnaire-9.  Medical care, pages 1194–1201, 2004.
  * Mate et al. [2020] Aditya Mate, Jackson Killian, Haifeng Xu, Andrew Perrault, and Milind Tambe.  Collapsing bandits and their application to public health interventions.  In Neural Information Processing Systems, NeurIPS, 2020.
  * Mate et al. [2021] Aditya Mate, Andrew Perrault, and Milind Tambe.  Risk-aware interventions in public health: Planning with restless multi-armed bandits.  In Proceedings of the 20th International Conference on Autonomous Agents and MultiAgent Systems, pages 880–888, 2021.
  * Mundorf et al. [2018] Christopher Mundorf, Arti Shankar, Tracy Moran, Sherry Heller, Anna Hassan, Emily Harville, and Maureen Lichtveld.  Reducing the risk of postpartum depression in a low-income community through a community health worker intervention.  Maternal and child health journal, 22(4):520–528, 2018.
  * Newman et al. [2018] Patrick M Newman, Molly F Franke, Jafet Arrieta, Hector Carrasco, Patrick Elliott, Hugo Flores, Alexandra Friedman, Sophia Graham, Luis Martinez, Lindsay Palazuelos, et al.  Community health workers improve disease control and medication adherence among patients with diabetes and/or hypertension in chiapas, mexico: an observational stepped-wedge study.  BMJ Global Health, 2018.
  * Nishtala et al. [2020] Siddharth Nishtala, Harshavardhan Kamarthi, Divy Thakkar, Dhyanesh Narayanan, Anirudh Grama, Ramesh Padmanabhan, Neha Madhiwalla, Suresh Chaudhary, Balaraman Ravindra, and Milind Tambe.  Missed calls, automated calls and health support: Using ai to improve maternal health outcomes by increasing program engagement.  arXiv preprint arXiv:2006.07590, 2020.
  * Papadimitriou and Tsitsiklis [1994] Christos H Papadimitriou and John N Tsitsiklis.  The complexity of optimal queueing network control.  In Proceedings of IEEE 9th Annual Conference on Structure in Complexity Theory, pages 318–322. IEEE, 1994.
  * Qian et al. [2016] Y. Qian, C. Zhang, B. Krishnamachari, and B. Tambe.  Restless poachers: Handling exploration- exploitation tradeoffs in security domains.  In International Joint Conference on Autonomous Agents and Multi-Agent Systems, AAMAS. IFAAMAS, 2016.
  * Rahedi Ong’ang’o et al. [2014] Jane Rahedi Ong’ang’o, Christina Mwachari, Hillary Kipruto, and Simon Karanja.  The effects on tuberculosis treatment adherence from utilising community health workers: a comparison of selected rural and urban settings in kenya.  PLoS One, 9(2):e88937, 2014.
  * Sombabu et al. [2020] B. Sombabu, A. Mate, D. Manjunath, and S. Moharir.  Whittle index for aoi-aware scheduling.  In IEEE International Conference on Communication Systems & Networks (COSMSNETS). IEEE, 2020.
  * Son et al. [2010] Youn-Jung Son, Hong-Gee Kim, Eung-Hee Kim, Sangsup Choi, and Soo-Kyoung Lee.  Application of support vector machine for prediction of medication adherence in heart failure patients.  Healthcare informatics research, 16(4):253–259, 2010.
  * Thaddeus and Maine [1994] Sereen Thaddeus and Deborah Maine.  Too far to walk: maternal mortality in context.  Social science & medicine, 38(8):1091–1110, 1994.
  * Watkins and Dayan [1992] Christopher JCH Watkins and Peter Dayan.  Q-learning.  Machine learning, 8(3-4):279–292, 1992.
  * Whittle [1988] Peter Whittle.  Restless bandits: Activity allocation in a changing world.  Journal of applied probability, pages 287–298, 1988.

##  Appendix A Proof for Whittle Index Q-Learning

###### Theorem 3

Taking action $1$ on the top $M$ arms according to
$(Q^{*}_{i}(s_{i},1)-Q^{*}_{i}(s_{i},0))$ is equivalent to maximizing
$\sum_{i}Q_{i}^{*}(s_{i},a_{i})$ over all possible action profiles
$(a_{1},\ldots,a_{N})\in\\{0,1\\}^{N}$ such that $\sum_{i}a_{i}=M$.

Proof Let $I^{*}$ be the set containing the top $M$ beneficiaries (i.e.,
beneficiaries with highest values of $Q_{i}(s_{i},1)-Q_{i}(s_{i},0)$), and
$I^{-*}$ be the set contains the non top M beneficiaries. Let $I_{M}$
represent any set of k beneficiaries. We then have:

$\displaystyle\sum_{i^{*}\in
I^{*}}[Q^{*}_{i^{*}}(s_{i^{*}},1)-Q^{*}_{i^{*}}(s_{i^{*}},0)]$
$\displaystyle\quad\qquad\qquad\geq\sum_{j\in
I_{M}}[Q^{*}_{j}(s_{j},1)-Q^{*}_{j}(s_{j},0)]$ Adding $\sum_{i\notin
I^{*},i\notin I_{M}}Q_{i}^{*}(s_{i},0)$ on both sides, we obtain:
$\displaystyle\sum_{i^{*}\in I^{*}}Q_{i*}^{*}(s_{i^{*}},1)+\sum_{i\in
I^{-*}}Q_{i}^{*}(s_{i},0)$ $\displaystyle\quad\qquad\qquad\geq\sum_{i\in
I_{M}}Q_{i}(s_{i},1)+\sum_{i\notin I_{M}}Q_{i}(s_{i},0)$ (7)

Equation 7 shows that providing intervention to beneficiaries in $I^{*}$ and
taking passive actions on all other arms would maximize
$\sum_{i}Q^{*}(s_{i},a_{i})$.

Now the rest of the proof follows from (1) the convergence guarantee of
Q-Learning algorithm Watkins and Dayan [1992], (2) $\epsilon$-decay selection
policy, and (3) the indexability assumption of RMABs Whittle [1988]. It has
been established in Watkins and Dayan [1992] that the update rule of
Q-Learning converges to $Q^{*}$ whenever $\sum_{t}\alpha_{t}(Z,a)=\infty$ and
$\sum_{t}\alpha_{t}(Z,a)^{2}<\infty.$ These assumptions require that all
state-action pairs be visited infinitely often, that is, the underlying MDP do
not have a state $Z$ such that after some finite time $t$, $Z$ is not
reachable via any sequence of actions. Similarly, the assumptions
$\displaystyle\sum_{h}\alpha(c^{h}_{i,Z,a})=\infty\mbox{ and
}\sum_{h}\alpha(c^{h}_{i,Z,a})^{2}<\infty$ require that each arm, state and
action pair would be visited infinitely often. Thus, in addition to the
reachability assumption on MDP, we need to show that each arm is visited
infinitely often. This is guaranteed by the $\epsilon$-decay selection
process, where each arm has a non-zero probability of being selected uniformly
at random. Thus, $Q^{t}_{i}(Z,a)$ converges to $Q^{*}_{i}(Z,a)$ which implies
that $Q^{t}_{i}(Z,1)-Q^{t}_{i}(Z,0)$ converges to
$Q^{*}_{i}(Z,1)-Q^{*}_{i}(Z,0)$ (using series convergence operation). Note
that the value $Q^{*}_{i}(Z,1)-Q^{*}_{i}(Z,0)$ is the Whittle Index. Using
these values for selecting the best arms would lead to an optimal solution for
the relaxed Lagrangian problem, because of the indexability assumption
(Definition 1). $\hfill\square$

##  Appendix B Experimental Evaluation

Here we provide the details of the benchmark algorithms that we use for
evaluation. Also, we provide examples from two other domains, one of them was
considered in Avrachenkov and Borkar [2020] and the other was considered in Fu
et al. [2019].

###  B.1 Benchmark Algorithms

We compare WIQL with five other algorithms:

  * •

OPT: This algorithm assumes full knowledge of the underlying transition
probabilities and has access to the optimal Whittle Indices
$\lambda_{i}^{*}(Z)$ for each arm $i$ and each state $Z$. At each time step
$t$, the top $M$ arms are chosen according to the $\lambda^{*}$ values of
their current states. OPT is our benchmark to visualize the convergence of
other algorithms.

  * •

AB: Proposed by Avrachenkov and Borkar [2020] this method is based on a
fundamental change in the definition of Q values that aims at converging to
the optimal Whittle Index policy. Arms are assumed to be homogeneous and they
use a shared Q update step for all the arms. The shared update helps in
collecting the point samples of Q values very quickly and results in fast
convergence to the optimal algorithm. However, we consider scenarios where the
arms may have different transition probabilities and we need to differentiate
among arms, even when they are in the same state. Thus, for our experiments,
we store Q values separately for each arm.

  * •

Fu: This method has been proposed by Fu et al. [2019] to update Q values
separately for each arm, each $\lambda\in\Lambda$, and each state-action pair.
They assume $\Lambda$ to be a set of input parameters. At each time step,
$\lambda_{min}\in\Lambda$ is computed which represents the subsidy that
minimizes the gap between Q values for taking action $1$ and that of action
$0$. For the experiments, we set $\Lambda$ and other hyperparameters the same
as what they considered in their paper. This method is a heuristic with no
convergence guarantee.

  * •

Greedy: This method maintains a running average reward for each arm, state,
and action, depending on the history of observations. At each time step, it
greedily chooses the top $M$ arms with the highest difference in their average
rewards between actions $1$ and $0$ at their current states.

  * •

Random: This method chooses $M$ arms uniformly at random at each step. Note
that this algorithm is the same as setting $\epsilon=1$ in Step $4$ of the
WIQL algorithm.

###  B.2 Numerical Examples

In this section, we provide two additional numerical examples for more general
RMAB instances.

####  B.2.1 Example with restart

This RMAB problem is considered by Avrachenkov and Borkar [2020]. Here, each
arm is assumed to be in one of the five states $\mathcal{S}=\\{0,\ldots,4\\}$
at any point of time. There are two actions—active ($a=1$) and passive
($a=0$). The active action forces an arm to restart from the first state. At
each time step, action $a=1$ can be taken only on $M$ arms. The reward from an
arm $i$ at a state $Z\in\mathcal{S}$ is assumed to be $R^{Z}_{i}=0.9^{Z}$ when
passive action is taken. The transition probabilities for each action
$a\in\\{0,1\\}$ are represented as a $|\mathcal{S}|\times|\mathcal{S}|$
matrix, where $p^{1}=1.0$, $q^{1}=0.0$, $p^{0}=0.1$, and $q^{0}=0.9$.

$\mathcal{P}^{a}=\begin{pmatrix}p^{a}&q^{a}&0&0&0\\\ p^{a}&0&q^{a}&0&0\\\
p^{a}&0&0&q^{a}&0\\\ p^{a}&0&0&0&q^{a}\\\ p^{a}&0&0&0&q^{a}\end{pmatrix}$

.

(a) Average reward accumulated when only $1$ arm out of $5$ arms can be chosen
at each iteration.

(b) Average reward accumulated when only $20$ arms out of $100$ arms can be
chosen at each iteration.

Figure 4: Comparing various algorithms for the example with Restart. The
results are averaged over $30$ trials.

(a) Average reward accumulated when only $1$ arm out of $5$ arms can be chosen
at each iteration.

(b) Average reward accumulated when only $20$ arms out of $100$ arms can be
chosen at each iteration.

Figure 5: Comparing various algorithms for Mentoring Instructions example. The
results are averaged over $30$ trials.

##### Results.

Figure 4 compares the performance of our proposed method with the benchmarks.
Similar to the Circulant Dynamics example, we observe that WIQL performs
better in terms of the average total reward, compared to the other benchmark
algorithms for the Example with Restart problem where the number of states is
higher.

####  B.2.2 Mentoring Instructions

This RMAB problem is considered by Fu et al. [2019]. Here, each arm,
representing a student is assumed to be in one of the ten states
$\mathcal{S}=\\{0,\ldots,9\\}$ as any point of time. At each time step, only
$M$ mentors are available and thus, action $a=1$ can be taken only on $M$
arms, which would possibly improve the state of the student. The higher index
of the state implies better reward (state $9$ represents the best study level
and state $1$ represents the worst study level). In particular, the reward at
a state $Z\in\mathcal{S}$ is assumed to be $R^{Z}=\sqrt{\frac{Z}{10}}$. The
transition probabilities for each action $a\in\\{0,1\\}$ are represented as a
$|\mathcal{S}|\times|\mathcal{S}|$ matrix, where $p^{1}=0.7$, $q^{1}=0.3$,
$p^{0}=0.7$, and $q^{0}=0.3$.

$\mathcal{P}^{a}=\begin{pmatrix}q^{a}&p^{a}&&&\\\ q^{a}&0&\ddots&&\\\
&q^{a}&\ddots&p^{a}&\\\ &&\ddots&0&p^{a}\\\ &&&q^{a}&p^{a}\end{pmatrix}$

##### Results.

Figure 5 compares the performance of our proposed method with the benchmarks.
We observe that WIQL performs marginally better in terms of the average total
reward, compared to the other benchmark algorithms for the Mentoring
Instructions problem. However, unlike the previous example, the margin of
improvement in performance is very low. A possible reason is the number of
states being much more than that of the earlier example, which takes more time
for the algorithms to learn the best policy. However, this example clearly
shows that Random policy may perform even worse when the number of states is
as high as 10 states.