# Hydrodynamic effects on the energy transfer from dipoles to metal slab

Daniel Brown  Hai-Yao Deng<EMAIL_ADDRESS>School of Physics and
Astronomy, Cardiff University, 5 The Parade, Cardiff CF24 3AA, Wales, United
Kingdom

###### Abstract

A systematic study of nonlocal and size effects on the energy transfer of a
dipole (e.g. a molecule or a quantum dot) induced by the proximity of a metal
slab is presented. Nonlocal effects are accounted for using the hydrodynamic
model (HDM). We derive a general relation that connects the energy transfer
rate to the linear charge density-density response function of the slab. This
function is explicitly evaluated for the HDM and the local Drude model. We
show that a thin metal slab can support a series of higher-frequency surface
plasma wave (SPW) modes in addition to the normal SPW modes thanks to the
nonlocal effects. These modes markedly alter the response and the energy
transfer process, as revealed in the structure of the energy transfer rate in
the parameter space. Our findings are important for applications such as the
recently developed metal-induced energy transfer imaging, which relies on
accurate modeling of the energy transfer rate.

##  I Introduction

As well known since the 1970s, the radiation properties of a dipole, which may
be represented by a molecule or a quantum dot (QD), can be strongly affected
by a nearby metal morawitz1969 ; kuhn1970 ; barton1970 ; chance1974 ;
morzwitz1974 ; chance1975 ; philpott1975 ; moskovits1978 . For example
pineda1985 ; nustovit2012 , fluorescence can be quenched as the dipole
channels most of its energy to the metal in the form of charge density waves
propagating along the metal surface, known as surface plasma waves (SPW).
These waves can also modify the dipole-dipole interaction and enable coherent
energy transmission between dipoles far apart philpott1975 ; stuart1998 .
Using conducting structures to obtain desired properties of an adjacent object
has become an versatile tool in a plethora of research areas kimble2008 ;
curto2010 ; dey2018 ; grob2018 ; chandra2020 ; lu2020 .

In the past few years, a powerful sub-nanometer imaging technique, called
Metal-Induced Energy Transfer Imaging (MIETI) was developed chizhik2014 ;
karedla2014 ; chizhik2017 ; baronsky2017 ; isbaner2018 ; ghosh2019 , utilizing
the sensitive dependence of the fluorescence time (rather than the radiative
power) on the dipole-metal distance. MIETI relies on accurate modeling of the
dipole-to-metal energy transfer rate. However, most theoretical studies
morawitz1969 ; kuhn1970 ; barton1970 ; chance1974 ; morzwitz1974 ; chance1975
; philpott1975 ; moskovits1978 ; cao2021  treat the metal as a simple
dielectric by essentially the Drude model and non-local effects have been
ignored. In studies where those effects were included vagov2016 ; Deng2017a ,
only a semi-infinite metal was considered rather than a metal slab, which is
usually the geometry in reality wu2009 .

In the present work, we systematically explore the non-local effects, by means
of a hydrodynamic model (HDM) barton1979 ; pendry2013 ; Deng2019 ; Deng2020 ,
on the energy transfer rate as a function of the dipole-metal distance and the
slab thickness. We find that hydrodynamic effects – together with size effects
– can significantly impact the energy transfer process, entailing a very
different picture than ignoring the effects. In this investigation, we also
provide a general relation between the transfer rate and the electrodynamic
response function of the metal, which is applicable beyond the HDM.

In the next section, we describe the system and the formalism, explaining the
connection between the energy transfer rate and the charge density-density
response function of the metal and also giving explicit expressions for the
energy transfer rate in the HDM and the LDM. In Sec. III, we present
comprehensive numerical results, revealing the existence of a sequence of
higher-frequency SPW branches and demonstrating their significance in the
structure of the energy transfer rate in the parameter space spanned by the
intrinsic dipole frequency, the slab thickness and the dipole-slab distance.
We summarize the paper in Sec. IV, wherein we also discuss experimental
relevance of the results.

##  II System and formalism

The system under study is depicted in Fig. 1, where a dipole exemplified by a
two-level QD is placed at a distance $a$ over the surface of a metal slab of
thickness $L$. The two levels of the QD are denoted by $|\uparrow\rangle$ and
$|\downarrow\rangle$, respectively, which are separated in energy by
$\hbar\omega$, where $\hbar$ is the reduced Planck constant and $\Omega$ is
the associated frequency. The Hamiltonian of the system can be written as a
sum of three parts,
$H_{QD}=\frac{\hbar\omega}{2}\left(|\uparrow\rangle\langle\uparrow|-|\downarrow\rangle\langle\downarrow|\right)$,
$H_{\text{Metal}}=\sum_{M}\varepsilon_{M}|M\rangle\langle M|$ and
$H^{\prime}$, which describe the QD, the metal and the coupling between them,
respectively. Here $M$ labels the eigenstates, $|M\rangle$, and eigenvalues,
$\varepsilon_{M}$ of $H_{\text{Metal}}$. We shall assume that the distance $a$
is not too small so that electronic tunneling can be ignored between the QD
and the metal. We also neglect retardation effects. As such, the coupling is
of electrostatic character and can be written as

$H^{\prime}=e\int
d^{3}\mathbf{x}d^{3}\mathbf{y}\frac{\hat{\rho}(\mathbf{x})\hat{n}(\mathbf{y})}{\mathinner{\\!\left\lvert\mathbf{x}-\mathbf{y}\right\rvert}},$
(1)

where $e$ denotes the charge of an electron, $\hat{n}$ is the electron density
operator subtracted by the mean density of electrons in the QD and
$\hat{\rho}$ is the charge density operator for the metal.

Figure 1: Schematic of the system: a dipole, represented by a QD is placed a
distance $a$ over a metal slab of thickness $L$. The QD has two energy levels
$|\uparrow\rangle$ and $\downarrow\rangle$, apart by $\hbar\omega$ in energy.
The two surfaces of the slab are located at $z=0$ and $z=L>0$, respectively.

Now we imagine the QD being excited to its upper state $|\uparrow\rangle$ by
absorbing a photon and set to calculate the rate at which the absorbed energy
is transferred to the metal by the coupling Hamiltonian $H^{\prime}$. We
denote the rate by $\frac{1}{\tau}$, where $\tau$ may be called the non-
radiative lifetime of the QD. The fluorescence life time can be calculated as
$\frac{\tau\tau^{\prime}}{\tau+\tau^{\prime}}$, where $\tau^{\prime}$ denotes
the radiative lifetime, which is to be treated in a future publication. In the
rest of this section, we describe the formalism for calculating $\tau$ as a
function of the system parameters.

###  II.1 Connection with the density-density response function

The initial state for the global system is written as
$|\uparrow,M\rangle=|\uparrow\rangle\otimes|M\rangle$, where the metal has
been assumed in state $|M\rangle$. Under the action of $H^{\prime}$, the QD
can jump to the lower state $|\downarrow\rangle$ by giving out energy to the
metal which synchronously transits to another state $|M^{\prime}\rangle$. The
transition amplitude reads

$\displaystyle\langle\downarrow,M^{\prime}|H^{\prime}|\uparrow,M\rangle=\int
d\mathbf{x}~{}\Phi^{*}(\mathbf{x})\hat{\rho}_{M^{\prime},M}(\mathbf{x}),$ (2)

where $\hat{\rho}_{M^{\prime},M}(\mathbf{x})=\langle
M^{\prime}|\hat{\rho}(\mathbf{x})|M\rangle$ and

$\Phi(\mathbf{x})=e\int
d\mathbf{y}\frac{\langle\uparrow|\hat{n}(\mathbf{y})|\downarrow\rangle}{|\mathbf{y}-\mathbf{x}|}$
(3)

is the potential generated by the QD. The total rate that the QD transfers its
energy to the metal can be obtained from Fermi’s golden rule, which gives

$\displaystyle\frac{1}{\tau}=\frac{2\pi}{\hbar}\sum_{M,M^{\prime}}p_{M}\mathinner{\\!\left\lvert\int
d\mathbf{x}~{}\Phi(\mathbf{x})\hat{\rho}_{M,M^{\prime}}(\mathbf{x})\right\rvert}^{2}\delta(\hbar\omega+\epsilon_{M}-\epsilon_{M^{\prime}}).$
(4)

Here we have assumed that the metal is initially in thermodynamic equilibrium
characterized by the Boltzmann factor $p_{M}=e^{-\varepsilon_{M}/k_{B}T}/Z$
with $Z$ being the partition function, $k_{B}$ the Boltzmann constant and $T$
the temperature, and $\delta(x)$ denotes the Dirac function.

We can relate $\frac{1}{\tau}$ as given in Eq. (4) to the charge density-
density response function, denoted by $\chi_{\omega}(\mathbf{x},\mathbf{y})$
of the metal. Physically, $\chi_{\omega}(\mathbf{x},\mathbf{y})$ gives the
charge density at $\mathbf{x}$ induced by an external oscillatory (at
frequency $\omega$) electrostatic potential localized at $\mathbf{y}$. Quantum
linear response theory gives Deng2020 ; lubensky

$\chi_{\omega}(\mathbf{x},\mathbf{y})=\sum_{M,M^{\prime}}p_{M}\left(\frac{\hat{\rho}_{M,M^{\prime}}(\mathbf{x})\hat{\rho}_{M^{\prime},M}(\mathbf{y})}{\hbar\omega+i0_{+}+\epsilon_{M}-\epsilon_{M^{\prime}}}-\frac{\hat{\rho}_{M,M^{\prime}}(\mathbf{y})\hat{\rho}_{M^{\prime},M}(\mathbf{x})}{\hbar\omega+i0_{+}+\epsilon_{M^{\prime}}-\epsilon_{M}}\right).$
(5)

Using Im$\left(\frac{1}{x+i0_{+}}\right)=-\pi\delta(x)$ and the fluctuation-
dissipation theorem lubensky , one may readily verify that

$\displaystyle\frac{1}{\tau}=\frac{2}{\hbar}\frac{1}{e^{-\beta\hbar\omega}-1}\mbox{Im}\left(\Gamma_{\omega}\right)\approx-\frac{2}{\hbar}~{}\mbox{Im}\left(\Gamma_{\omega}\right),$
(6)

where the approximation occurs at temperatures much lower than
$\hbar\omega/k_{B}$, which is assumed throughout the paper, and

$\Gamma_{\omega}=\int d\mathbf{x}\int
d\mathbf{y}~{}\Phi(\mathbf{x})\chi_{\omega}(\mathbf{x},\mathbf{y})\Phi^{*}(\mathbf{y}).$
(7)

Equations (6) and (7) are generally valid, not just for the QD coupled to a
metal but to any material, as long as the appropriate $\chi_{\omega}$ is
supplied. Equation (7) describes the net effect of the charges induced at
point $\mathbf{y}$ by the QD propagating to $\mathbf{x}$ whereby they interact
with the QD again.

Figure 2: Map of
$\log_{10}\left[\mbox{Im}\left(\frac{1}{\epsilon_{+}(k,\omega)}\right)+\mbox{Im}\left(\frac{1}{\epsilon_{-}(k,\omega)}\right)\right]$,
which reveals the dispersion of SPWs admitted on a slab of thickness $L$. SPW
branches are indicated by the numbers. The non-locality parameter $v$ is
expressed in unit of $\omega_{p}d$. Very small damping
$\gamma/\omega_{p}=0.01$ is used. Unlike LDM, HDM allows higher-frequency
modes to appear on a slab.

###  II.2 The dipole approximation

Considering that the QD is small, we may invoke the dipole approximation for
$\Phi(\mathbf{x})$, i.e.

$\Phi(\mathbf{x})\approx\mathbf{D}\cdot\partial_{\mathbf{x}_{0}}\frac{1}{\mathinner{\\!\left\lvert\mathbf{x}_{0}-\mathbf{x}\right\rvert}},$
(8)

where $\mathbf{x}_{0}=(\mathbf{r}_{0},z_{0})$ with $\mathbf{r}_{0}=(0,0)$ and
$z_{0}=-a<0$ locates the center of the QD, and
$\mathbf{D}=(D_{x},D_{y},D_{z})$ denotes the dipole transition element of the
QD, given by

$\mathbf{D}=e\int
d\mathbf{y}~{}\mathbf{y}\langle\uparrow|\hat{n}(\mathbf{y})|\downarrow\rangle,$
(9)

Note that the integral here is virtually carried over the volume of the QD.

We may make use of the translation symmetry along the slab surfaces by writing

$\chi_{\omega}(\mathbf{x},\mathbf{y})=\int\frac{d^{2}\mathbf{k}}{4\pi^{2}}~{}\chi_{\mathbf{k}\omega}(z,z^{\prime})e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}^{\prime})}.$
(10)

Here $\mathbf{x}=(\mathbf{r},z)$,
$\mathbf{y}=(\mathbf{r}^{\prime},z^{\prime})$ and $\mathbf{k}$ is a wave
vector. Also we introduce a double cosine Fourier transform, defined as

$\chi_{\mathbf{k}\omega}(z,z^{\prime})=\sum^{\infty}_{m,n=0}\frac{1}{L_{m}}\chi_{\mathbf{k}\omega}(n,m)\cos(q_{n}z)\cos(q_{m}z^{\prime}),$
(11)

Here $q_{n}=\frac{n\pi}{L}$ and $L_{m}=\frac{L}{2-\delta_{m,0}}$ with
$\delta_{m,m}$ being the Kronecker symbol. Physically,
$e^{i(\mathbf{k}\cdot\mathbf{r}-\omega t)}\chi_{\mathbf{k}\omega}(n,m)$ gives
the charge density induced by an external electrostatic potential of this
form: $e^{i(\mathbf{k}\cdot\mathbf{r}-\omega t)}\cos(q_{m}z)$. Via Eqs. (8),
(10) and (11), one arrives at

$\Gamma_{\omega}=\sum_{m,n}\sum_{\mu,\nu}D_{\mu}D^{*}_{\nu}\int
d^{2}\mathbf{k}\frac{\chi_{\mathbf{k}\omega}(n,m)}{L_{m}}\frac{\bar{k}_{\mu}e^{-ka}}{k^{2}+q^{2}_{n}}\frac{\bar{k}^{*}_{\nu}e^{-ka}}{k^{2}+q^{2}_{m}}\left(1-(-1)^{n}e^{-kL}\right)\left(1-(-1)^{m}e^{-kL}\right),$
(12)

where $k=|\mathbf{k}|$ and $\bar{\mathbf{k}}=(i\mathbf{k},k)$. Note that the
factor $e^{-ka}$ puts an effective cutoff on $k$: Fourier components with
$k>1/a$ cannot make significant contributions.

To make progress, we assume that the slab is invariant under reflection about
its mid-plane at $z=L/2$. It follows that quantities with different parities
cannot be mixed Deng2017b . In particular, $\chi_{\mathbf{k}\omega}(n,m)$ must
vanish if $m$ and $n$ have different parities. It suggests to split
$\chi_{\omega}$ into two sectors, the symmetric one
$\chi^{+}_{\mathbf{k}\omega}(l,l^{\prime})=\chi_{\mathbf{k}\omega}(2l,2l^{\prime})$
and the anti-symmetric one
$\chi^{-}_{\mathbf{k}\omega}(l,l^{\prime})=\chi_{\mathbf{k}\omega}(2l+1,2l^{\prime}+1)$,
where $l,l^{\prime}=0,1,...$. Additionally, we introduce the following
quantities,

$P^{s}_{k\omega}(l)=\sum^{\infty}_{l^{\prime}=0}\frac{\chi^{s}_{\mathbf{k}\omega}(l^{\prime},l)}{k^{2}+(q^{s}_{l^{\prime}})^{2}},~{}P_{s,\omega}(k)=\sum_{l}\frac{4\pi}{L^{s}_{l}}\frac{P^{s}_{k\omega}(l)}{k^{2}+(q^{s}_{l})^{2}},$
(13)

where $s=\pm$ and $q^{+}_{l}=q_{2l}$ and $q^{-}_{l}=q_{2l+1}$. We have
considered the fact that $\chi^{s}_{\mathbf{k}\omega}(l)$ does not depend on
the direction of $\mathbf{k}$ due to rotational symmetry about the $z$-axis.
These symmetry considerations allow us to finally get

$\frac{1}{\tau}=\frac{1}{2}\frac{e^{2}}{\hbar
d}\left(1+\cos^{2}\theta\right)R(\omega,a,L),$ (14)

where $\cos\theta=D_{z}/\mathinner{\\!\left\lvert\mathbf{D}\right\rvert}$ and
$d=\mathinner{\\!\left\lvert\mathbf{D}\right\rvert}/e$ is a length comparable
to the dimension of the QD, and

$R(\omega,a,L)=-\sum_{s}\int^{\infty}_{0}dk~{}k^{3}e^{-2ka}\left(1-se^{-kL}\right)^{2}~{}\mbox{Im}\left(P_{s,\omega}(k)\right)$
(15)

is a dimensionless quantity. In obtaining the expression for $R$, we have
performed an integration over the angle of $\mathbf{k}$ and, in each quantity
in the integrand the unit of length has been set with $d=1$ to simplify the
notation.

A few remarks are in order. Equation (15) shows that the QD probes the
quantity $P_{s,\omega}(k)$ rather than $\chi_{\mathbf{k}\omega}$. This accords
with the general observation Deng2020  that entities (e.g. the QD here)
residing outside a metal can only probe $P^{s}_{k\omega}$. Secondly, the
frequency, $\frac{e^{2}}{\hbar d}$, should be comparable to the QD intrinsic
frequency $\omega$, as the latter is roughly set by the size of the QD, which
is of the same order as $d$. Finally, Eq. (14) shows that the orientation of
$\mathbf{D}$ also affects the energy transfer rate, but only by a factor of
$2$ at the most.

Figure 3: Anti-symmetric modes manifested with
Im$\left({\frac{1}{\epsilon_{-}(k,\omega)}}\right)$ at fixed $kd=3$ for
various slab thickness $L$ and non-locality parameter $v$. Numbers in (a)
indicate SPW branches. $\gamma/\omega_{p}=0.01$. Legends for (a) and (b) are
the same as in (c).

###  II.3 Hydrodynamic model (HDM)

HDM treats the electrons in a metal as a fluid, governed by the linearized
Navier-Stokes equation barton1979 ; pendry2013 ; Deng2019 ; Deng2020

$mn_{0}\left(\partial_{t}+\gamma\right)\mathbf{V}=n_{0}e\mathbf{E}-mv^{2}\nabla
n.$ (16)

Here $m$ is the mass of an electron, $n_{0}$ the mean electron density, $n$
denotes the deviation from $n_{0}$, $\gamma$ is a phenomenological parameter
accounting for dissipation (viscous effects), $\mathbf{E}(\mathbf{x},t)$ the
total electric field felt by the electrons and $\mathbf{V}(\mathbf{x},t)$
their velocity field. In addition, $v$ is a parameter that controls non-local
effects. The local Drude model (LDM) is recovered as a limit of the HDM with
$v=0$.

In the Appendix, we show that for HDM

$\displaystyle
P^{s}_{k\omega}(l)=\frac{\omega^{2}_{p}}{4\pi}\frac{1}{\tilde{\omega}^{2}-\omega^{2}_{p}-v^{2}(q^{s}_{l})^{2}}\frac{1}{\epsilon_{s}(k,\omega)},$
(17)

where $\omega_{p}=\sqrt{4\pi n_{0}e^{2}/m}$ denotes the characteristic
frequency of the metal, and

$\epsilon_{s}(k,\omega)=\frac{\bar{\omega}^{2}-\Omega^{2}_{s}}{\bar{\omega}^{2}-\omega^{2}_{p}},\quad\Omega_{s}=\omega_{s}\sqrt{1+\frac{k}{Q}\frac{1+se^{-QL}}{1-se^{-QL}}}.$
(18)

Here $\omega_{s}=\frac{\omega_{p}}{\sqrt{2}}\left(1-se^{-kL}\right)$ and
$Q=v^{-1}\sqrt{\omega^{2}_{p}-\tilde{\omega}^{2}}$ with
$\tilde{\omega}^{2}=\bar{\omega}^{2}-k^{2}v^{2}$ and
$\bar{\omega}=\omega+i\gamma$. Using the identity that

$\sum^{\infty}_{l=0}\frac{1}{L^{s}_{l}}\frac{1}{q^{2}+(q^{s}_{l})^{2}}=\frac{1}{2q}\frac{1+se^{-qL}}{1-se^{-qL}},$
(19)

which is an even function of $q$, one easily finds

$P_{s,\omega}(k)=\left(1-\frac{1}{\epsilon_{s}(k,\omega)}\right)\frac{1}{k\left(1-se^{-kL}\right)}.$
(20)

As such, we get

$R(\omega,a,L)=\sum_{s}\int^{\infty}_{0}dk~{}k^{2}e^{-2ka}\left(1-se^{-kL}\right)~{}\mbox{Im}\left(\frac{1}{\epsilon_{s}(k,\omega)}\right).$
(21)

This equation, together with Eqs. (14) and (18), completely determines the
energy transfer rate $1/\tau$.

Equation (21) implies that the elementary excitations that can be detected by
the QD are all ciphered in the quantity $\epsilon_{s}(k,\omega)$, the zeros of
which in the complex $\omega$-plane locate the frequency and damping rate of
SPWs propagating with wavenumber $k$ along the slab surfaces, satisfying
$\bar{\omega}=\Omega_{s}(\bar{\omega})$.

The results for LDM are particularly simple: With $v=0$, $\Omega_{s}$ reduces
to $\omega_{s}$ and then

$\epsilon_{s}=\frac{\bar{\omega}^{2}-\omega^{2}_{s}}{\bar{\omega}^{2}-\omega^{2}_{p}},\quad
P^{s}_{k\omega}(l)=\frac{\omega^{2}_{p}}{4\pi}\frac{1}{\bar{\omega}^{2}-\omega^{2}_{s}}.$
(22)

The zeros of $\epsilon_{s}$ occur at $\omega=\omega_{s}-i\gamma$. Thus,
$\omega_{s}$ and $\gamma$ give the reputed frequency and the damping rate,
respectively, of SPWs admitted on a slab by LDM. For $k\gg L^{-1}$,
$\omega_{s}(k)$ reduces to $\omega_{p}/\sqrt{2}$ and becomes dispersion-less.
Splitting the integral in Eq. $(\ref{HDMR})$ into two parts, one over
$k\in[0,1/L]$ and the other $k\in[L^{-1},\infty]$, one may see that the latter
part gives a contribution to $R$ going as
$\frac{1}{2}\mbox{Im}\left(\frac{\omega^{2}_{p}}{\omega^{2}_{p}-2\bar{\omega}^{2}}\right)\left(\frac{d}{a}\right)^{3}$,
which blows up to infinity for $a$ tending to zero regardless of whether
$\omega$ is in resonance or not. Divergence of this kind is well known
horsley2014  and removed by nonlocal effects in HDM.

##  III Numerical results and discussions

###  III.1 Properties of surface plasma waves

The QD is localized in space and hence can give out energy to excite SPWs of
any wavenumber. As is clear from Eq. (21), the efficiency for SPWs of
wavenumber $k$ to get excited by the QD is largely set by the spectral weight

$\mbox{Im}\left(\frac{1}{\epsilon_{s}(k,\omega)}\right)=\frac{\epsilon^{\prime\prime}_{s}(k,\omega)}{\left[\epsilon^{\prime}_{s}(k,\omega)\right]^{2}+\left[\epsilon^{\prime\prime}_{s}(k,\omega)\right]^{2}},$
(23)

where we have put
$\epsilon_{s}=\epsilon^{\prime}_{s}-i\epsilon^{\prime\prime}_{s}$ with
$\epsilon^{\prime}_{s}$ and $-\epsilon^{\prime\prime}_{s}$ being the real and
imaginary parts of $\epsilon_{s}$, respectively. The imaginary part
$\epsilon^{\prime\prime}_{s}$ originates from dissipation of SPWs into heat,
which is ultimately traced back to the damping parameter $\gamma$ in HDM. As
$\gamma$ is typically small, $\epsilon^{\prime\prime}_{s}$ is also generally
small. Equation (23) then displays a resonance wherever
$\epsilon^{\prime}_{s}(k,\omega)=0$, which determines the SPW dispersion
relation.

To reveal the SPW dispersion, we have produced a map of dB
$=\log_{10}\left[\mbox{Im}\left(\frac{1}{\epsilon_{+}(k,\omega)}\right)+\mbox{Im}\left(\frac{1}{\epsilon_{-}(k,\omega)}\right)\right]$
in the $k-\omega$ plane, as exhibited in Fig. 2, for various values of the
non-locality parameter $v$. The SPW resonances coalescence into bright lines
of high dB, representing different branches.

It is usually thought that there are merely two branches of SPWs living on a
slab, stemming from the symmetric ($s=+$) and the anti-symmetric ($s=-$)
sectors, respectively. This is obviously true with LDM, corresponding to
$v=0$, for which the SPW has dispersion given by
$\omega_{s}=\frac{\omega_{p}}{\sqrt{2}}\sqrt{1-se^{-kL}}$ as already mentioned
in the preceding section [cf. Eq.(22)]. For small $kL$,
$\omega_{+}\propto\sqrt{kL}$ whereas $\omega_{-}\approx\omega_{p}$. For large
$kL$, both $\omega_{+}$ and $\omega_{-}$ tend to $\omega_{p}$. These features
are well borne out in the bottom row of panels (i) - (l) in Fig. 2. For
convenience, we shall label these ordinary branches by $l=+1$ and $l=-1$,
respectively, see the figure.

The SPW properties by HDM are strikingly different. In addition to the normal
branches seen in LDM ($l=\pm 1$), there are higher-frequency branches, which
we label as $l=\pm 2,\pm 3,...$ in the order of increasing frequency (see Fig.
2). The branches with $l>0$ all derive from the symmetric sector
$\epsilon_{+}$ while those with $l<0$ from the anti-symmetric sector
$\epsilon_{-}$. They have less intensity than the primary branches $l=\pm 1$
but can still be clearly seen in the top and middle rows of panels (a) - (h)
in Fig. 2. In Fig. 3, we plot as an illustration
$\mbox{Im}\left(\frac{1}{\epsilon_{-}(k,\omega)}\right)$ at fixed $kd=3$
versus $\omega$ to show that the branches $l<-1$ are absent from LDM. Similar
plot can be generated for the symmetric branches $l>1$, which are also absent
from LDM. Thus, the modes $|l|>1$ arise from hydrodynamic effects.

The symmetric and anti-symmetric branches are seen alternating with one
another, and they are rather evenly spaced in frequency. The spacing decreases
with increasing slab thickness $L$. Upon decreasing $L$, their frequencies all
(including that of mode $l=+1$) shift upward dramatically. Actually, for very
small $L$, even the lowest anti-symmetric branch $l=-1$ can be pushed way
higher than $\omega_{-}$ (see the first column of panels, (a), (e) and (i) in
Fig. 2). On the other hand, for large $L$, these higher-frequency branches get
closer and less pronounced and could be easily erased by damping. We note in
passing that the symmetric branches have less intensity (dB) at long
wavelengths (i.e. small $k$) than the anti-symmetric ones, a feature that
allows us to tell them apart by simple inspection.

Figure 4: Higher-frequency SPW modes appear wherever the factor $\exp(-QL)$
changes its phase from $\pi$ to $-\pi$, which gives rise to the anti-symmetric
modes $l<0$, or wherever its phase angle is zero, which gives rise to the
symmetric modes $l>1$. Numbers indicate branches. Parameters:
$v=0.1d\omega_{p}$, $L/d=1$, $\gamma=0.01\omega_{p}$ and $kd=3$.

Thanks to the hydrodynamic effects, the dispersion of the branch $l=-1$ is not
monotonic in HDM [Fig. 2 (d)]: the frequency decreases first and then
increases with increasing $k$. The branch $l=+1$ can also undergo a change in
its dispersion if $v$ is big, which is almost linear for small $L$ [Fig. 2
(a)] but becomes concave as $L$ increases [Fig. 2 (b) - (d)].

Properties of the higher-frequency branches can be understood from the
behaviors of the factor $\exp(-QL)$ in Eq. (18). As demonstrated in Fig. 4,
the branches $l>1$ occur wherever the phase angle of $\exp(-QL)$ approaches
zero, while the modes $l<-1$ occur wherever the phase angle crosses $\pi$. As
a result, the symmetric and anti-symmetric branches alternate in frequency.
The phase difference by $\pi$ comes from the minus sign picked up when going
from $s=+1$ to $s=-1$, as is clear from Eq. (18). For $\omega>\omega_{p}$,
$kv\ll\omega$ and $\gamma/\omega_{p}\ll 1$, we find
$Q\approx\frac{i\omega_{p}}{v}\sqrt{(\omega/\omega_{p})^{2}-1}$ having
neglected a small imaginary part $\sim\gamma/v$. Under such circumstance
$\exp(-QL)$ is purely oscillatory with $\omega$. For $\omega/\omega_{p}\gg 1$,
the oscillation period approximates a constant $\Delta\omega=2\pi v/L$, which
explains why these branches are nearly evenly spaced. Meanwhile, it shows that
the spacing decreases with decreasing $v/L$, as observed in Figs. 2 and 3.

Figure 5: Same as Fig. 2 (c) but for wider range of $k$, showing that higher-
frequency branches are dispersive at large $k$. $L=0.5d$ and $v$ is in unit of
$d\omega_{p}$.

If $kv>\omega_{p}$, $Q$ will depend on $k$ via $\omega_{p}$ replaced by
$\sqrt{\omega^{2}_{p}+k^{2}v^{2}}$. Then the spacing also depends on $k$ and
the branches become dispersive with frequencies varying with $k$. This is
demonstrated in Fig. 5.

Figure 6: Im$\left(\epsilon^{-1}_{s_{l}}(k,\omega)\right)$ at fixed $\omega$.
For strongly dispersive branch $l=+1$, it features a sharp peak centered at
$k_{\omega}$ satisfying $\omega=\omega_{l}(k_{\omega},\omega)$, where
$\omega_{l}$ is the dispersion relation of the branch $l$. For barely
dispersive branch $l=-1$, the peak is broad. The red curve (for $l=-1$) has
been multiplied by $20$ for visualization. Inset:
Im$\left(\epsilon^{-1}_{s_{l}}(k,\omega)\right)$ for branch $l=-1$ of LDM,
which is dispersion-less at large $k$. Parameters: $v=0.1d\omega_{p},L=0.1d$
and $\gamma=0.05\omega_{p}$. Figure 7: $\log_{10}R$ mapped unto the $a-L$
plane for various values of $\omega$. Panels (k) - (o) are cuts at fixed $L$
as indicated, showing that $R$ diverges in LDM but remains finite in HDM for
$a$ approaching zero. $\gamma=0.05\omega_{p}$ and $v$ is in unit of
$d\omega_{p}$.

For very large $L\gamma/v\gg 1$, $\exp(-QL)$ is negligibly small due to the
small but finite imaginary part of $Q$. In such case, the higher-frequency
branches shall disappear.

Figure 8: $\log_{10}R$ mapped unto the $\omega-a$ plane for various values of
$L$. $\gamma=0.05\omega_{p}$ and $v$ is in unit of $d\omega_{p}$.  Figure 9:
Resonances with higher-frequency SPW modes revealed in the oscillatory
behaviors of $R$. $R_{\text{max}}$ denotes the maximum of $R$ for the
corresponding curve. Parameters: $\gamma=0.01\omega_{p}$, $L=0.5d$ and $v$ is
in unit of $d\omega_{p}$.

###  III.2 Energy transfer rate $1/\tau$

After illustrating the structure of the SPWs supported on slab surfaces, we
now look into how energy is transmitted from the QD to the slab, generating
SPWs in the latter. The quantity of interest is $R$ [c.f. Eq. (21)], which
directly determines the energy transfer rate $\frac{1}{\tau}$ according to Eq.
(14).

Figure 10: $\log R$ mapped unto the $\omega-L$ plane for various values of
$a$. $\gamma=0.05\omega_{p}$ and $v$ is in unit of $d\omega_{p}$.

For HDM, as is clear from Fig. 2, the spectral weight
$\sum_{s}\mbox{Im}\left(\frac{1}{\epsilon_{s}(k,\omega)}\right)$ is
significant only for $(k,\omega)$ lying on one of the SPW branches (near or on
resonance). Suppose $\omega$ cuts the $l$-th branch at $k_{\omega}$, i.e.
$\omega=\omega_{l}(k_{\omega})$, where $\omega_{l}(k)$ denotes the dispersion
relation for this branch. Two situations might happen:

  1. (A)

If the branch is highly dispersive around $k_{\omega}$ (e.g. $l=+1$), then
only the modes with $k$ near $k_{\omega}$ can significantly contribute:
$\mbox{Im}\left(\frac{1}{\epsilon_{s}(k,\omega)}\right)$ displays a sharp peak
centered at $k_{\omega}$ with a small breadth $\Delta k$ that depends on the
group velocity $d\omega_{l}(k)/dk$ and the damping parameter $\gamma$ [Fig. 6
black curve]. From Eq. (21), we estimate that
$R(\omega,a,L)\approx\left(k_{\omega}de^{-k_{\omega}a}\right)^{2}\left(1-s_{l}e^{-k_{\omega}L}\right)W_{l}(\omega).$
Here $s_{l}$ denotes the sign of $l$ and $W_{l}(\omega)\approx d\Delta
k\cdot\mbox{Im}\left(\frac{1}{\epsilon_{s_{l}}(k_{\omega},\omega)}\right)$ is
the total spectral weight of the peak.

  2. (B)

If the branch is barely dispersive (e.g. low-$k$ sections of branches
$|l|>1$), a wide section of length $k_{c}$ could contribute [Fig. 6 red
curve]. We may then estimate that
$R(\omega,a,L)\approx\alpha\mbox{Im}\left[\epsilon^{-1}_{-}(k,\omega)\right]\int^{k_{c}}_{0}dk~{}k^{2}e^{-2ka}\left(1-s_{l}e^{-kL}\right)$,
where $\alpha$ is a parameter that depends on the shape of the section in
question. The integral can be explicitly evaluated but not needed here. As $L$
increases beyond $2a$, $R$ rises to a constant value
$\sim(d/a)^{3}\left(1-e^{-2k_{c}a}(1+2k_{c}a+2k^{2}_{c}a^{2})\right)\approx\frac{4k^{3}_{c}d^{3}}{3}e^{-2k_{c}a}$,
where the approximation occurs for $k_{c}a$ not too big.

In both situations, $R$ rises exponentially to a constant as $a$ goes to zero,
i.e. $\log R\propto a$.

In LDM, SPWs are strictly non-dispersive for $k\gg L^{-1}$ [inset of Fig. 6].
This is similar to situation (B) but with $k_{c}\rightarrow\infty$. As such,
$R$ diverges like $a^{-3}$ for $a$ approaching zero, as mentioned in Sec.
II.3. Note that this is so even for $\omega$ off resonance with any branches.

In Fig. 7, we display a map of $\log_{10}R$ in the $a-L$ plane at fixed
$\omega$ of various values ranging from $0.5\omega_{p}$ to $1.75\omega_{p}$.
It is clear – in particular from the panels (k) - (o) – that the above
scenarios, i.e. the blow-up of $R$ in LDM and its hinderance in HDM, are
vindicated. In comparison with LDM, $R$ displays a much richer structure in
HDM. For LDM, regarding its dependence on $L$, $R$ simply approaches a
constant $~{}\left(d/a\right)^{3}$ as $L$ increases [panels (f) - (j)]. For
HDM, however, the $L$ dependence of $R$ is not monotonic for any $\omega$. For
small $\omega$ that only resonates with the branch $l=+1$, $R$ displays a
single peak at some value $L$ [panels (a) and (b)]. This can be understood
from the SPW dispersion relations shown in Fig. 2. For very small $L$, the
dispersion $\omega_{+1}(k)$ resonates with $\omega$ at a very large
$k_{\omega}$, hence a very small $R$ due to the cutoff factor
$e^{-2k_{\omega}a}$. Increasing $L$ brings $k_{\omega}$ to less values and
larger $R$. However, as $L$ increases a crossover in the shape of
$\omega_{+1}(k)$ happens: it evolves from linear to concave, resulting in
smaller $\Delta k$ and a decrease in $R$ from its peak value located where the
crossover takes place. On the other hand, for $\omega$ above $\omega_{-1}(k)$,
$R$ displays an oscillatory structure featuring a number of peaks [panels (c)
- (e)], which each indicates a resonance with one of the higher-frequency
branches. Actually, as $L$ increases, these branches move down and closer and
then pass through $\omega$ one by one [see Fig. 2], producing the oscillations
seen in $R$.

The higher-frequency modes are even more manifest in the $\omega-a$ map of $R$
at fixed $L$, shown in Fig. 8 and Fig. 9. As seen in panels (a) - (d) here,
$R$ displays very clear oscillations signifying resonances with these modes at
not so small $a$. At very small $a$, though, the contributions from these
modes are blurred. The happens again due to the cutoff factor $e^{-2ka}$: for
very small $a$, this factor allows modes with very large $k$ to contribute,
which have much closer frequencies, and damping can easily blur their
individual contributions. For large $a$, only modes with small $k$ can
contribute significantly, and these modes are well spaced at least for not too
large $L$. If $L$ is very large, the higher-frequency modes are also closely
packed and damping erases the oscillatory structure in $R$, see panel (e). In
contrast, in LDM $R$ is rather pedestrian: it displays a two-branch structure
at small $L$ and $a$, which is blurred otherwise for analogous reasons, see
panels (f) - (j).

Finally, in Fig. 10 we display a map of $\log R$ in the $\omega-L$ plane at
fixed $a$. For LDM, we see the usual SPW branches in $R$ at small $L$, which
merge into a single one at large $L$ [panels (f) - (j)]. For HDM, the
oscillatory structure appears on the side where $\omega>\omega_{+1}$,
signifying resonances with the higher-frequency branches [panels (a) - (e)].
For large $a$, the maps for LDM and HDM look rather similar, except for the
higher-frequency modes – which are clearly visible for $L<\sim 5d$ – and the
following difference: for LDM $R$ strictly vanishes as soon as $\omega$ goes
above $\omega_{p}$ [panel (j)], whereas for HDM it extends beyond [panel (e)].
This is so because hydrodynamic effects push up the frequencies of SPWs, as
discussed in the preceding subsection.

From Fig. 10 we may also observe that, at smaller $a$ the oscillations on the
high-$\omega$ side have shorter periods [panel (b)], having signals from both
the symmetric and anti-symmetric higher-frequency modes; However, as $a$
increases, the signals from the symmetric modes fade out [panel (c) versus
(b)], leaving only the anti-symmetric ones still discernible at fairly large
$a$ [panel (e)]. At very small $a$ [panel (a)], the oscillations are not clear
on the $\log R$ map, but seeable on the $R$ map, which is not shown here.

##  IV Summary

To summarize, we have systematically examined the nonlocal effects and size
effects on the energy transfer from a dipole to a metal slab. The non-local
effects are implemented using a hydrodynamic model. We have shown that the
energy transfer rate, $\tau^{-1}$, can be related to the charge density-
density response function of the metal slab. For HDM and LDM (as a limit of
HDM), the response functions have been explicitly evaluated. The most
significant finding is the existence of higher-frequency SPW modes supported
on thin metal slabs. We have seen that these modes qualitatively modify the
structure of $\tau^{-1}$ in the parameter space $(\omega,a,L)$. Our results
demonstrate that it is imperative to build nonlocal effects into the modeling
necessitated in metal-induced energy transfer imaging and other areas such as
relaying information between QDs by means of metals or semiconductors
kimble2008 .

To experimentally detect the higher-frequency SPWs, one must ensure that the
slab thickness $L$ does not exceed $\sim v/\gamma$; otherwise, the
hydrodynamic effects would be erased by damping. Typically, $v\sim
0.1d\omega_{p}$ and hence $L$ should not be larger than $10d$ for
$\gamma=0.01\omega_{p}$. For most QDs, $d\sim 10~{}$nm, and thus $L\sim
100~{}$nm at the most for such parameters. Meanwhile, it is essential to have
the QD frequency well tunable above the characteristic frequency of the metal,
i.e. $\omega>\omega_{p}$ in order to reach resonance with the higher-frequency
modes. Most QDs have $\omega$ amounting to a few electron volts (eVs) in
energy. The slab material should then have $\omega_{p}$ around $1~{}$eV or
less. Often experimented metals, such as silver and gold, have $\omega_{p}$
much higher and are not suitable. In addition, in those metals electrodynamic
responses cannot be simply modeled by the HDM due to effects from inner atomic
shells Deng2019 ; liebsch1993 ; johnson1976 . One way out is to use doped
semiconductor slabs rather than metals. For semiconductors $\omega_{p}$ can be
tuned by doping, temperature and external field sha1995 , and they may be the
most promising candidate for experimentally studying the phenomena predicted
in the present work.

## Acknowledgement

We acknowledge the support of the Supercomputing Wales project, which is part-
funded by the European Regional Development Fund via Welsh Government.

*

##  Appendix A Slab density-density response function

Here we derive the electrodynamic response functions for HDM. Assuming a time
dependence $e^{-i\omega t}$, the dynamic equation (16) can be recast as

$-i\bar{\omega}mn_{0}\mathbf{V}=n_{0}e\mathbf{E}-mv^{2}\nabla
n,\quad\bar{\omega}=\omega+i\gamma.$ (24)

In terms of the current density $\mathbf{J}=n_{0}e\mathbf{V}$ and charge
density $\rho=ne$ due to the electrons,

$-i\bar{\omega}m\mathbf{J}=n_{0}e^{2}\mathbf{E}-mv^{2}\nabla\rho.$ (25)

Note that $\nabla\cdot\mathbf{E}=4\pi(\rho+\rho_{ext})$, where $\rho_{ext}$ is
the external charge density. Taking the divergence,

$-i\bar{\omega}\nabla\cdot\mathbf{J}=\omega^{2}_{p}(\rho+\rho_{ext})-v^{2}\nabla^{2}\rho.$
(26)

On the other hand, we have the equation of continuity, which reads
$i\bar{\omega}\rho=\nabla\cdot\mathbf{j}$ with
$\mathbf{j}=\mathbf{J}\left(\theta(z)-\theta(z-L)\right)$, where $\theta(z)$
denotes the Heaviside step function: $\theta(z\leq 0)=0$ and $\theta(z>0)=1$.
Combining this with Eq. (26) gives

$(v^{2}\nabla^{2}+\bar{\omega}^{2}-\omega^{2}_{p})\rho=\omega^{2}_{p}\rho_{ext}-i\bar{\omega}\left(J_{z}(\mathbf{x}_{0})\theta^{\prime}(z)-J_{z}(\mathbf{x}_{L})\theta^{\prime}(z-L)\right),$
(27)

where $\mathbf{x}_{0}=(\mathbf{r},0)$ and $\mathbf{x}_{L}=(\mathbf{r},L)$
denote a point on the surfaces at $z=0$ and $z=L$, respectively. From Eq.
(25),

$J_{z}(\mathbf{x})=\frac{in_{0}e^{2}}{m\bar{\omega}}E_{z}(\mathbf{x})+\frac{v^{2}}{i\bar{\omega}}\partial_{z}\rho(\mathbf{x}).$
(28)

In the literature barton1979 ; pendry2013 , it is often imposed that
$J_{z}(\mathbf{x}_{0})=J_{z}(\mathbf{x}_{L})=0$. Here we do not make such
imposition.

Making use of the translation symmetry along the surfaces and assuming an
overall dependence $e^{i\mathbf{k}\cdot\mathbf{r}}$ on the planar coordinates
$\mathbf{r}$, that is $\rho(\mathbf{r})=\rho(z)e^{i\mathbf{k}\cdot\mathbf{r}}$
and similarly for other field quantities, Eqs. (4) and (28) combine to yield

$\displaystyle(v^{2}\partial^{2}_{z}+\tilde{\omega}^{2}-\omega^{2}_{p})\rho(z)=\left(\frac{\omega^{2}_{p}}{4\pi}E_{z}(0)-v^{2}\rho^{\prime}(0)\right)\theta^{\prime}(z)$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad-\left(\frac{\omega^{2}_{p}}{4\pi}E_{z}(L)-v^{2}\rho^{\prime}(L)\right)\theta^{\prime}(z-L)+\omega^{2}_{p}\rho_{ext}.$

Doing a cosine Fourier transform to this equation leads to

$(\tilde{\omega}^{2}-\omega^{2}_{p}-q^{2}_{n}v^{2})\rho_{n}=\frac{\omega^{2}_{p}}{4\pi
L_{n}}\left[E_{z}(0)-(-1)^{n}E_{z}(L)\right]+\omega^{2}_{p}\rho_{ext,n},$ (30)

where $\rho_{n}$ are the cosine Fourier coefficients,

$\rho(z)=\sum_{n}\rho_{n}\cos(q_{n}z),\quad\rho_{n}=L^{-1}_{n}\int^{L}_{0}dz~{}\rho(z)\cos(q_{n}z).$
(31)

Note that the terms involving $\rho^{\prime}(0)$ or $\rho^{\prime}(0)$ in Eq.
(A) have cancelled out and hence disappeared from Eq. (30), which implies that
the results do not depend on whether $J_{z}$ vanishes at the surfaces or not.

One may show that

$E_{z}(0)=-\Phi^{\prime}_{ext}(0)-2\pi\xi_{0},\quad
E_{z}(L)=-\Phi^{\prime}_{ext}(L)+2\pi\xi_{L},$ (32)

where $\Phi_{ext}(z)$ is the external potential generated by the probe charge
$\rho_{ext}$, and

$\displaystyle\xi_{0}$ $\displaystyle=$
$\displaystyle\int^{L}_{0}dze^{-kz}\rho(z)=\sum_{n}\frac{k\rho_{n}}{k^{2}+q^{2}_{n}}\left(1-e^{-kL}(-1)^{n}\right),$
(33) $\displaystyle\xi_{L}$ $\displaystyle=$
$\displaystyle\int^{L}_{0}dze^{k(z-L)}\rho(z)=\sum_{n}\frac{k\rho_{n}(-1)^{n}}{k^{2}+q^{2}_{n}}\left(1-e^{-kL}(-1)^{n}\right).$

Considering that $(k^{2}-\partial^{2}_{z})\Phi_{ext}(z)=4\pi\rho_{ext}$, one
finds

$\frac{k^{2}+q^{2}_{n}}{4\pi}\Phi_{ext,n}=\rho_{ext,n}-\frac{\Phi^{\prime}_{ext}(0)-(-1)^{n}\Phi^{\prime}_{ext}(L)}{4\pi
L_{n}}.$ (35)

Using these results, Eq. (30) becomes

$\left(\tilde{\omega}^{2}-\omega^{2}_{p}-q^{2}_{n}v^{2}\right)\rho_{n}=\frac{\omega^{2}_{p}(k^{2}+q^{2}_{n})}{4\pi}\Phi_{ext,n}-\frac{\omega^{2}_{p}}{2L_{n}}(\xi_{0}+(-1)^{n}\xi_{L}).$
(36)

Note that

$\xi_{0}+(-1)^{n}\xi_{L}=\sum_{m}\frac{k\rho_{m}}{k^{2}+q^{2}_{m}}\left(1-e^{-kL}(-1)^{m}\right)\left(1+(-1)^{m+n}\right),$
(37)

which vanishes unless $n$ and $m$ have the same parity. This is a consequence
of the reflection symmetry possessed by the slab about its mid-plane. Thus, we
split $\rho(z)=\rho^{+}(z)+\rho^{-}(z)$, where
$\rho^{\pm}(z)=\sum^{\infty}_{l=0}\rho^{\pm}_{l}\cos(q^{\pm}_{l}z)$ with
$\rho^{+}_{l}=\rho_{2l}$, $\rho^{-}_{l}=\rho_{2l+1}$ and $q^{\pm}_{l}$ defined
in the main text, and introduce

$\xi_{\pm}=\frac{1}{2}(\xi_{0}\pm\xi_{L})=\sum_{l}\frac{k\rho^{\pm}_{l}}{k^{2}+(q^{\pm}_{l})^{2}}\left(1\mp
e^{-kL}\right).$ (38)

Under such provisions, Eq. (36) can be rewritten as

$\left(\tilde{\omega}^{2}-\omega^{2}_{p}-(q^{\pm}_{l})^{2}v^{2}\right)\rho^{\pm}_{l}=\frac{\omega^{2}_{p}(k^{2}+(q^{\pm}_{l})^{2})}{4\pi}\Phi^{\pm}_{ext,l}-\frac{\omega^{2}_{p}}{L^{\pm}_{l}}\xi_{\pm}.$
(39)

The above two equations can now be solved to produce

$\rho^{\pm}_{l}=\sum_{l^{\prime}}\chi^{\pm}_{k\omega}(l,l^{\prime})\Phi^{\pm}_{ext,l^{\prime}}$
(40)

with the density-density response function given by

$\displaystyle\chi^{\pm}_{k\omega}(l,l^{\prime})=\frac{\omega^{2}_{p}}{\tilde{\omega}^{2}-\omega^{2}_{p}-v^{2}(q^{\pm}_{l})^{2}}$
$\displaystyle\quad\times\left[\frac{k^{2}+(q^{\pm}_{l})^{2}}{4\pi}\delta_{l,l^{\prime}}-\frac{1}{\epsilon_{\pm}}\frac{\omega^{2}_{p}}{\tilde{\omega}^{2}-\omega^{2}_{p}-v^{2}(q^{\pm}_{l^{\prime}})^{2}}\frac{k(1\mp
e^{-kL})}{4\pi L^{\pm}_{l}}\right].$

For $L\rightarrow\infty$, this result reduces to what was derived in Ref.
[Deng2020, ] for a semi-infinite metal.

In Eq. A, we have reinstated the subscript $k\omega$, and the quantity

$\epsilon_{s}=1+\sum_{l}\frac{2\omega^{2}_{s}}{\tilde{\omega}^{2}-\omega^{2}_{p}-v^{2}(q^{s}_{l})^{2}}\frac{k/L^{s}_{l}}{k^{2}+(q^{s}_{l})^{2}}$
(42)

can be rewritten as

$\epsilon_{s}=1+\frac{2k\omega^{2}_{s}}{\bar{\omega}^{2}-\omega^{2}_{p}}\sum^{\infty}_{l=0}\frac{1}{L^{s}_{l}}\left(\frac{1}{k^{2}+(q^{s}_{l})^{2}}-\frac{1}{Q^{2}+(q^{s}_{l})^{2}}\right)$
(43)

with $Q^{2}=(\omega^{2}_{p}-\tilde{\omega}^{2})/v^{2}$. Using the identity Eq.
(19), we immediately get

$\epsilon_{s}(k,\omega)=1+\frac{\omega^{2}_{-s}}{\bar{\omega}^{2}-\omega^{2}_{p}}\left(1-\frac{k}{Q}\frac{1+se^{-QL}}{1-se^{-QL}}\frac{1-se^{-kL}}{1+se^{-kL}}\right),$
(44)

which after algebra gives Eq. (18). Using Eq. (A) in Eq. (13), one is led to
Eq. (17) for $P^{s}_{k\omega}(l)$ in the main text.

Setting $v=0$, one gets

$\displaystyle\chi^{\pm}_{k\omega}(l,l^{\prime})$ $\displaystyle=$
$\displaystyle\frac{\omega^{2}_{p}}{\bar{\omega}^{2}-\omega^{2}_{p}}\frac{k^{2}+(q^{\pm}_{l})^{2}}{4\pi}\delta_{l,l^{\prime}}$
(45) $\displaystyle\quad\quad\quad-\frac{k}{2\pi
L^{\pm}_{l}}\frac{\omega^{2}_{p}}{\bar{\omega}^{2}-\omega^{2}_{p}}\frac{\omega^{2}_{\pm
0}}{\bar{\omega}^{2}-\omega^{2}_{\pm 0}}$

as the response function for LDM.

## References

  * (1) H. Morawitz, Phys. Rev. 187, 1792 (1969); Phys. Rev. A 7, 1148 (1973).
  * (2) H. Kuhn, J. Chem. Phys. 53, 101 (1970).
  * (3) G. Barton,, Proc. R. Soc. A 320, 251 (1970).
  * (4) R. R. Chance, A. Prock and R. Silbey, J. Chem. Phys. 60, 2744 (1974).
  * (5) H. Morawitz and M. R. Phipott, Phys. Rev. B 10, 4863 (1974).
  * (6) R. R. Chance, J. Chem. Phys. 63, 1589 (1975).
  * (7) M. R. Philpott, J. Chem. Phys. 62, 1812 (1975).
  * (8) M. Moskovits, J. Chem. Phys. 69, 1 (1978).
  * (9) A. C. Pineda and D. Ronis, J. Chem. Phys. 83, 5330 (1985).
  * (10) V. N. Nustovit and T. V. Shahbazyan, J. Chem. Phys. 136, 204701 (2012).
  * (11) H. R. Stuart and D. G. Hall, Phys. Rev. Lett. 80, 5663 (1998).
  * (12) H. J. Kimble, Nature 453, 1023 (2008).
  * (13) A. G. Curto, G. Volpe, T. H. Taminiau, M. P. Kreuzer, R. Quidant and N. F. van Hulst, Science 329, 930 (2010).
  * (14) S. Dey, Y. Zhou, Y. Sun, J. Jenkins, D. Kriz, S Suib, O. Chen, S. Zou and J. Zhao, Nanoscale 10, 1038 (2018).
  * (15) H. Gro$\beta$, J. M. Hamm, T. Tufarelli, O. Hess and B. Hecht, Sci. Adv. 4, 4906 (2018).
  * (16) S. Chandra, H. Ahmed and S. McCormack, Opt. Express 28, 20191 (2020).
  * (17) H. Lu, Z. Huang, M. S. Martinez, J. C. Johnson, J. M. Luther and M. C. Beard, Energy Environ. Sci. 13, 1347 (2020).
  * (18) A. I. Chizhik, J. Rother, I. Gregor, A. Janshoff and J. Enderlein, Nat. Photonics 8, 124 (2014).
  * (19) N. Karedla, A. I. Chizhik, I. Gregor, A. M. Chizhik, O. Schulz and J. Enderlein, ChemPhysChem 15, 705 (2014).
  * (20) A. M. Chizhik, D. Ruhlandt, J. Pfaff, N. Karedla, A. I. Chizhik, I. Gregor, R. H. Kehlenbach and J. Enderlein, ACS Nano 11, 11839 (2017).
  * (21) T. Baronsky, D. Ruhlandt, B. R. Brükner, J. Schäfer, N. Karedla, S. Isbaner, D. Hähnel, I. Gregor, J. Enderlein, A. Janshoff, A. I. Chizhik, Nano Lett. 17, 3320 (2017).
  * (22) S. Isbaner, N. Karedla, I. Kaminska, D. Rnhlandt, M. Raab, J. Bohlen, A. Chizhik, I. Gregor, P. Tinnefeld, J. Enderlein and R. Tsukanov, Nano Lett. 18, 2616 (2018).
  * (23) A. Ghosh, A. Sharma, A. I. Chizhik, S. Isbaner, D. Ruhlandt, R. Tsukanov, I. Gregor, N. Karedla and J. Enderlein, Nat. Photonics 13, 860 (2019).
  * (24) Q. Cao, Z. Lin, X. Li and Z. Zhang, Phys. Rev. B 103, 035416 (2021).
  * (25) A. Vagov, I. A. Larkin, M. D. Croitoru and V. M. Axt, Phys. Rev. B 93, 195414 (2016).
  * (26) H.-Y. Deng, K. Wakabayashi and C.-H. Lam, Phys. Rev. B 95, 045428 (2017); H.-Y. Deng, J. Phys. Condens. Matter 29, 455002 (2017).
  * (27) X. Wu, Y. Sun and M. Pelton, Phys. Chem. Chem. Phys. 11, 5867 (2009).
  * (28) G. Barton, Rep. Prog. Phys. 42, 963 (1979).
  * (29) C. Ciracć, J. B. Pendry and D. R. Smith, ChemPhysChem 14, 1109 (2013).
  * (30) H.-Y. Deng, New J. Phys. 21, 043055 (2019); arXiv 2005.10760 (2020).
  * (31) H.-Y. Deng, Ann. Phys. 418, 168204 (2020).
  * (32) P. Chaikin and T. C. Lubensky, Principles of condensed matter physics (Cambridge University Press, Cambridge, 1995).
  * (33) H.-Y. Deng, Phys. Rev. B 95, 125442 (2017).
  * (34) S. A. R. Horsley and T. G. Philbin, New J. Phys. 16, 013030 (2014).
  * (35) A. Liebsch, Phys. Rev. Lett. 71, 145 (1993).
  * (36) P. B. Johnson and R. Christy, Phys. Rev. B 6, 4370 (1976).
  * (37) W. Sha, A. L. Smirl and W. F. Tseng, Phys. Rev. Lett. 74, 4273 (1995).


# Deploying Vaccine Distribution Sites for Improved Accessibility and Equity
to Support Pandemic Response

George Z. Li2,∗ Ann Li1 Madhav Marathe1,3 Aravind Srinivasan2 Leonidas
Tsepenekas2 Anil Vullikanti1,3

###### Abstract

In response to COVID-19, many countries have mandated social distancing and
banned large group gatherings in order to slow down the spread of SARS-CoV-2.
These social interventions along with vaccines remain the best way forward to
reduce the spread of SARS CoV-2. In order to increase vaccine accessibility,
states such as Virginia have deployed mobile vaccination centers to distribute
vaccines across the state. When choosing where to place these sites, there are
two important factors to take into account: accessibility and equity. We
formulate a combinatorial problem that captures these factors and then develop
efficient algorithms with theoretical guarantees on both of these aspects.
Furthermore, we study the inherent hardness of the problem, and demonstrate
strong impossibility results. Finally, we run computational experiments on
real-world data to show the efficacy of our methods.

††1Biocomplexity Institute and Initiative, University of Virginia††2Department
of Computer Science, University of Maryland††3Department of Computer Science,
University of Virginia††∗Correspondence to George Li at<EMAIL_ADDRESS>
##  1 Introduction

The COVID-19 pandemic continues to cause immense social, health, and economic
impact globally. As of writing this paper, the U.S. alone has seen over
850,000 deaths and over 65 million confirmed cases; see (vdh, ) for the latest
numbers. Vaccines have proven to be very effective in reducing the health
burden of the pandemic and continue to be the best strategy to control disease
spread and potentially end the pandemic in its current form. Despite the
effectiveness, administering COVID-19 vaccines to all eligible individuals in
the population continues to be a challenge. As of February 2022, only 64% of
the eligible population is fully vaccinated in the United States (nyt, b).
Furthermore, there is a significant disparity in vaccination rates between
demographics—the rate among Whites was 1.2 times that of African Americans and
1.1 times that of Hispanic people. The reasons why some people have not been
vaccinated include distrust and skepticism regarding COVID-19, accessibility
issues, and concerns about the cost (nyt, a). Lottery schemes, mandates,
vaccine clinics, and other strategies have been implemented to increase the
vaccination rate with varying levels of success. Since cost and accessibility
remain a challenge for a fraction of the population, especially minorities and
people in poorer neighborhoods, mobile vaccine clinics have been an important
part of the public health response strategy of government agencies. In this
paper, we study the problem of deploying mobile vaccine administration sites
with the goal of improving the accessibility of vaccines to individuals.

Deploying vaccination clinics is a form of a facility location problem (Celik
Turkoglu and Erol Genevois, 2020; Drezner, 1995), referred to as the
$k$-supplier problem, in which a limited set of $k$ facilities needs to be
placed so that every person (i.e., a client) is “close” to a facility; a
common metric to measure closeness is the maximum distance between a client
and their closest facility, though many other notions have been studied.
Facility location problems are well understood, and efficient approximation
algorithms and practical heuristics exist. However, deploying vaccine clinics
leads to a novel facility location problem (referred to as the
MobileVaccClinic problem) since people (clients) are mobile rather than
stationary. Suppose each person $p$ visits a set $S_{p}$ of locations during
the day; then it suffices to deploy a clinic close to at least one location in
$S_{p}$. Our contributions are the following:

  * •

We formalize the MobileVaccClinic problem for modeling the deployment of
mobile vaccine clinics in a way that takes into account human mobility
patterns (by considering the distance to a facility from any of the locations
visited by a person), fairness (by requiring that at least a fraction of
people in each demographic group have a nearby clinic), outliers (by allowing
partial coverage), and capacity constraints (by restricting the number of
people assigned to each clinic). We show that this problem is much harder than
the standard $k$-supplier problem and getting any bounded polynomial-time
approximation to the minimum distance is not possible, thus motivating
bicriteria algorithms.

  * •

We design two approximation algorithms. The first is a fixed-parameter
tractable algorithm that gives a $3$-approximation, where the parameterization
is on the number $u$ of locations where people travel. Note that even this is
non-trivial, since the possible locations $\mathcal{S}$ where we can place
facilities is still variable, so there are still $|\mathcal{S}|\choose k$
possible solutions. The second algorithm, based on covering problems, is a
$(1,\log{n}+1)$-bicriteria one, where $n$ is the number of people. This means
that if we violate the budget on the number of vaccine centers by a
$\log{n}+1$ multiplicative factor, we can find a solution that is optimal.
Finally, we extend our algorithms to have fairness guarantees in both the
original and outliers formulation of the problem.

  * •

We evaluate our algorithms for a realistic population of a county in Virginia.
We find that our algorithms generally give a significant improvement over
natural baselines. In particular, we see many shortcomings of only considering
a client’s home (rather than their entire travelling route), emphasizing the
importance of our problem formulation. Additionally, our algorithms allow us
to compute a tradeoff between the maximum distance to a clinic and the number
of clinics; this naturally enables us to give a recommendation to the
government on the most cost-effective budget policy. Finally, the solutions
computed by our algorithm have a useful “kernel” property—as the budget is
increased, the locations which were picked for a lower budget are still part
of the solution. This implies that an incrementally constructed solution
(which is how such facilities would be deployed in practice since the budget
is not known ahead of time) will still be good.

We remark that though our framework is motivated by the current COVID-19
pandemic, it can be generally applied to both epidemiological and non-
epidemiological settings. Examples within healthcare include placing testing
and treatment units (as deployed during the Ebola crisis) and delivering
healthcare in rural settings for resource-limited countries. Beyond
healthcare, the placement of mobile distribution centers arises in disaster-
management settings. For instance, shelters need to be set up for individuals
evacuating during a hurricane or forest fire, who might need food and other
basic survival kits. During such large events, mobile sites are also used to
place security posts and information kiosks.

##  2 Preliminaries

Recall that we wish to place vaccination centers such that vaccines are more
accessible to the population. This question is often formulated as an
appropriate variant of the facility location problem, which is well-studied in
the operations research literature (see Related Work). In our paper, we
introduce a new variant that follows a recent line of work on integrating the
mobility patterns of the population into disease models (Chang et al., 2021;
Wang et al., 2020). As is standard, we will use the distance from a
vaccination center as the metric for defining accessibility. The key change,
however, is that clients will be represented by a set of locations that they
visit (within a time period) instead of just one point. Though this will make
the problem much harder to solve efficiently, it will more strongly correlate
with the likelihood of a person going to a vaccine center.

Figure 1: An example of MobileVaccClinic. The different colors represent
different people and the circles represent the locations they visit (with the
bottom three being their homes). In this case, the blue location in the middle
is an optimal location to place a vaccination center. If we instead only
considered homes in the problem formulation, we would place the vaccination
center in the green circle marked with a star, which would require people to
deviate from their normal travels much more when getting a vaccine.

Problem Statement: We are given a set of locations $\mathcal{C}$ in a metric
space characterized by the distance function
$d:\mathcal{C}\times\mathcal{C}\mapsto\mathbb{R}_{\geq 0}$. We additionally
have a set of $n$ individuals/clients $\mathcal{P}$. Each individual
$p\in\mathcal{P}$ is associated with a set $S_{p}\subseteq\mathcal{C}$, which
we can interpret as the set of locations $p$ visits throughout the day.
Finally, the input also includes a positive integer $k$ constraining the
number of facilities we can place, and a set $\mathcal{S}\subseteq\mathcal{C}$
containing the locations where we are allowed to place facilities. The goal of
MobileVaccClinic is to choose a set $F\subseteq\mathcal{S}$ with $|F|\leq k$
to place facilities, such that for every $p\in\mathcal{P}$ we have
$d(S_{p},F)\leq R$, for the minimum $R$ possible. Here, we use the standard
notation where $d(S,F)=\min_{j\in S,j^{\prime}\in F}d(j,j^{\prime})$.
Intuitively, this objective tries to minimize the maximum distance between the
set of facilities placed and the locations visited by any client. We also
consider three natural extensions:

  * •

Outliers: in order to achieve herd immunity, we only need to vaccinate a large
portion of the population (rather than every single person). In order to model
this, we can take as input a parameter $q$, and seek to provide for only
$\lfloor qn\rfloor$ of the clients, thus ignoring the remaining ones.
Formally, the new objective is to minimize $R$ such that
$|\\{p\in\mathcal{P}:d(S_{p},F)\leq R\\}|\geq\lfloor qn\rfloor$.

  * •

Fairness: many studies have shown that COVID-19 disproportionately affects
some demographic groups (Tai et al., 2020). To counteract this, we seek to
guarantee that different demographic groups have similar accessibilities to
vaccines. As an example, when we solve the outliers formulation, we can
guarantee that we are covering the same proportion of each demographic group
when deciding the facility placements.

  * •

Capacity: it is natural to assume that the number of vaccines that can be
stored in each mobile facility is limited, say at most $L$. Therefore, in this
setting, we need to guarantee that every chosen facility will have at most $L$
people assigned to it.

##  3 Related Work

Due to its applications in a large number of domains, facility location and
broader location theory is a very well-studied area; see, e.g., the surveys by
(Drezner, 1995; Celik Turkoglu and Erol Genevois, 2020; Afshari and Peng,
2014). The general goal in this family of problems is to deploy facilities to
provide the best possible service to a set of clients. A huge number of
objectives have been considered, along with a plethora of variations such as
fairness variants and online or stochastic versions. The MobileVaccClinic
problem we study here is a generalization of the well-known $k$-center
problem, where the goal is to open at most $k$ centers while minimizing the
maximum distance of a point to its closest center. For this simple clustering
setting, there exist efficient $2$-approximation algorithms (Hochbaum and
Shmoys, 1985; Gonzalez, 1985). Furthermore, it is shown that unless P=NP this
is the best achievable approximation ratio (Hochbaum and Shmoys, 1986).

Location theory problems have also been considered in the area of healthcare,
e.g., (Nedjati and Valipour, 2012; Meskarian et al., 2017; de Vries et al.,
2020; Afshari and Peng, 2014). A lot of this work has been focused on placing
mobile clinics or temporary facilities to ensure good service, especially in
resource-poor countries. As mentioned in (Afshari and Peng, 2014), the
healthcare domain poses new challenges for location theory, such as
uncertainty, reliability, operation efficiency, patient safety, and cost-
effectiveness. Prior work has generally not considered the mobility of clients
at a detailed scale, which provides more flexibility in deploying facilities.
Our formulation of MobileVaccClinic explicitly models human mobility, thus
providing a realistic framework for public health agencies in their response
efforts.

##  4 Hardness Result

For our hardness result, we use the following problem studied in Anegg et al.
(2020), named $\gamma$-Colorful $k$-Center or $\gamma$C$k$C for short. This
problem is a generalization of the outliers version of $k$-center: in addition
to the classical constraints, colors (representing demographic groups) are
assigned to each client and the problem requires that a sufficient number of
points of each color is covered. The formal definition is given below:

######  Definition 4.1.

Let $\gamma\in\mathbb{Z}_{\geq 1}$ be the number of colors,
$k\in\mathbb{Z}_{\geq 1}$ be the budget, $\mathcal{C}$ be a set of points in a
metric space, and
$d:\mathcal{C}\times\mathcal{C}\xrightarrow{}{\mathbb{R}_{\geq 0}}$ be the
distance function on $\mathcal{C}$. For each $\ell\in[\gamma]$, let
$\mathcal{C}_{\ell}\subseteq\mathcal{C}$ be the points with color $\ell$ and
let $m_{\ell}\in\mathbb{Z}_{\geq 1}$ be the number of points with color $\ell$
which need to be covered. $\gamma$C$k$C asks for the minimum radius $R$
together with a set $F\subseteq\mathcal{C}$ with $|F|\leq k$, such that at
least $m_{\ell}$ points of $\mathcal{C}_{\ell}$ are covered within distance
$R$ by $F$. Formally, if $B(F,R)=\\{j\in\mathcal{C}:d(j,F)\leq R\\}$ then we
want $|B(F,R)\cap\mathcal{C}_{\ell}|\geq m_{\ell}$ for every
$\ell\in[\gamma]$.

In Anegg et al. (2020) the authors prove the following hardness result, which
we use to prove a hardness result for our problem later on.

######  Lemma 4.2.

When $\gamma$ is not a constant, there exist instances of $\gamma$C$k$C with
$m_{\ell}=1$ for all $\ell\in[\gamma]$, such that if $R^{*}$ is the optimal
value of the instance, the following hold:

  * •

For any $\rho>0$, it is NP-hard to find $F\subseteq\mathcal{C}$ with $|F|\leq
k$ and $|B(F,\rho R^{*})\cap\mathcal{C}_{\ell}|\geq m_{\ell}$ for all $\ell$.
In words, it is NP-hard to devise any approximation algorithm for
$\gamma$C$k$C.

  * •

For any $\rho>0$ and $\epsilon\in(0,1)$, it is NP-hard to find
$F\subseteq\mathcal{C}$ with $|F|\leq(1-\epsilon)\ln\gamma\cdot k$ and
$|B(F,\rho R^{*})\cap\mathcal{C}_{\ell}|\geq m_{\ell}$ for every
$\ell\in[\gamma]$. In words, it is NP-hard to devise any bicriteria
approximation for $\gamma$C$k$C, whose chosen centers will be at most
$(1-\epsilon)\ln\gamma\cdot k$.

  * •

These problematic instances consist of points on a line.

######  Remark 4.3.

Regarding the second statement in Lemma 4.2, the authors of Anegg et al.
(2020) show a bicriteria hardness result in terms of $\log|\mathcal{C}|$ and
not $\ln\gamma$. However, a closer look into their proof reveals that the
claim mentioned above follows trivially. We choose to present this form of the
claim because it better fits our narrative later on in the paper.

######  Theorem 4.4.

There exists a bicriteria preserving reduction of the problematic instances of
$\gamma$C$k$C described in Lemma 4.2 to instances of MobileVaccClinic with
$|\mathcal{P}|=\gamma$. Specifically, any $(\rho,\alpha)$-bicriteria
approximation for MobileVaccClinic translates to a $(\rho,\alpha)$-bicriteria
approximation for the problematic instances of $\gamma$C$k$C.

###### Proof.

Let
$(\mathcal{C},\mathcal{C}_{1},\ldots,\mathcal{C}_{\gamma},k,m_{1},\ldots,m_{\gamma})$
be a problematic instance of $\gamma$C$k$C as described in Lemma 4.2, and
recall that this instance has $m_{\ell}=1$ for all $\ell\in[\gamma]$. We will
now construct an instance of MobileVaccClinic as follows. The metric space for
MobileVaccClinic will be the same as in the $\gamma$C$k$C problem. That is, we
assume we have points $\mathcal{C}$ with a distance function $d$ on them. For
every $\ell\in[\gamma]$ construct a client $p_{\ell}$, and set
$S_{p_{\ell}}=\mathcal{C}_{\ell}$. The set of locations $\mathcal{S}$ for
MobileVaccClinic where we can place facilities will be the set of locations
$\mathcal{C}$ of $\gamma$C$k$C, and the value $k$ will stay the same for the
two problems.

Consider now the optimal solution $F^{*}$ of the $\gamma$C$k$C instance and
its corresponding value $R^{*}$. We claim that $F^{*}$ is a feasible solution
for the constructed MobileVaccClinic instance, and its value for that is
exactly the same. This is easy to see because $|F^{*}|\leq k$,
$|B(F^{*},R^{*})\cap\mathcal{C}_{\ell}|\geq 1$ for every $\ell\in[\gamma]$,
and $\mathcal{C}_{\ell}$ are exactly the locations visited by client
$p_{\ell}$. Hence if $R_{OPT}$ is the value of the optimal solution to the the
constructed MobileVaccClinic instance, we have $R_{OPT}\leq R^{*}$.

Take now any $(\rho,\alpha)$-bicriteria solution $F$ for MobileVaccClinic. At
first we trivially have $|F|\leq\alpha k$. Moreover, for every $\ell$ we can
express $d(F,S_{p_{\ell}})\leq\rho R_{OPT}$ (the condition guaranteed by the
$(\rho,\alpha)$-bicriteria solution $F$ for MobileVaccClinic) as $|B(F,\rho
R_{OPT})\cap\mathcal{C}_{\ell}|\geq 1$. Finally, because $R_{OPT}\leq R^{*}$,
we have $B(F,\rho R_{OPT})\subseteq B(F,\rho R^{*})$. Hence, $|B(F,\rho
R^{*})\cap\mathcal{C}_{\ell}|\geq|B(F,\rho R_{OPT})\cap\mathcal{C}_{\ell}|\geq
1$ for every $\ell\in[\gamma]$. The latter concludes the bicriteria preserving
reduction. ∎

######  Corollary 4.5.

Even when the metric space is the Euclidean line, we have the following for
MobileVaccClinic (unless P$=$NP):

  1. 1.

No approximation algorithm exists.

  2. 2.

Any bicriteria algorithm must use at least $k\ln n$ facilities.

##  5 Algorithms

In this section, we introduce efficient methods which give (approximately)
optimal facility placements, despite the hardness results. We also show how to
extend each of our algorithms to ignore outliers, incorporate fairness
constraints, and restrict the capacity of each facility.

###  5.1 Fixed-Parameter Tractability

Let $U=\bigcup_{p\in\mathcal{P}}S_{p}$ denote the set of all the locations
visited by the set of clients and $u=|U|$ be the number of locations in this
set. Due to potential privacy concerns, we can assume that the client
locations we have access to only include large public areas in the county such
as malls, shopping centers, etc. Hence, it is reasonable to conclude that $u$
is a fixed parameter, which we assume ranges from $15-30$. Given this fixed
parameter, we develop an efficient algorithm for our problem.

The main observation here is the following: consider an instance of
MobileVaccClinic and let $F^{*}$ be its optimal solution, whose maximum radius
we denote by $R^{*}$. For each $p\in\mathcal{P}$, we know that
$d(F^{*},S_{p})\leq R^{*}$, and hence there must exist a location $i_{p}\in
S_{p}$ with $d(i_{p},F^{*})\leq R^{*}$. See now that
$\\{i_{p}~{}|~{}p\in\mathcal{P}\\}\subseteq U$, and therefore
$|\\{i_{p}~{}|~{}p\in\mathcal{P}\\}|\leq u$. The latter implies that we can
guess, via an exhaustive search, the set $\\{i_{p}~{}|~{}p\in\mathcal{P}\\}$
in time at most $2^{u}$ (recall that since $u$ is considered a fixed
parameter, $2^{u}$ is thought of as a small constant). Let $A$ be the correct
guess for that set; we can think of $A$ as the set of locations through which
the optimal solution covers every client within distance $R^{*}$. Given $A$,
we see that the problem of computing $F^{*}$ reduces in a straightforward
manner to the well-known $k$-supplier problem (Hochbaum and Shmoys, 1985).

In $k$-supplier we have a set of points $\mathcal{X}$ and a set of locations
$\mathcal{Y}$ in a metric space with distance function $d$. The goal is to
choose $C\subseteq\mathcal{Y}$ with $|C|\leq k$, such that the maximum
distance of any point in $\mathcal{X}$ to its closest location of $S$ is
minimized. Hence, after correctly guessing $A$, we create an instance of
$k$-supplier where the points are the ones in $A$, and the set of locations
$\mathcal{Y}$ is $\mathcal{S}$. The previous discussion shows that $F^{*}$ is
a solution for this $k$-supplier instance, and its maximum radius will again
be $R^{*}$. Moreover, any $\rho$-approximate solution to the $k$-supplier
instance will trivially be a $\rho$-approximate solution for MobileVaccClinic.
Using the $3$-approximation algorithm from Hochbaum and Shmoys (1985) proves
the following theorem.

######  Theorem 5.1.

Algorithm 1 yields a $3$-approximation algorithm for MobileVaccClinic and runs
in time $2^{u}\operatorname*{poly}(n,|\mathcal{C}|)$.

Algorithm 1 FPT

1: for $A\in 2^{U}:|A\cap S_{p}|\neq 0,\forall p\in\mathcal{P}$ do

2: Obtain locations $F_{A}$ by running the $k$-supplier algorithm on the
appropriate instance discussed above.

3: Calculate the objective value for $F_{A}$.

4: end for

5: Pick the $F_{A}$ with the smallest objective value.

Moving forward, we see that the same approach of guessing the correct set of
client locations $A$ will also apply in different settings. In fact, the only
thing that may differ is the need for an alternative $k$-supplier algorithm
that can incorporate the specific constraints of each unique setting; we
survey some of these settings below.

Outliers: To modify our algorithm so that it only considers some fraction $q$
of the population, we only need to change the objective value evaluated in
line 3 of Algorithm 1. To improve efficiency, we can also only consider
guesses $A$ that contain locations from at least $\left\lfloor
qn\right\rfloor$ clients since the correct guess $A$ contains locations from
at least $\left\lfloor qn\right\rfloor$ clients. If we then feed $A$ to the
$k$-supplier algorithm in the exact same manner before, we will get a
$3$-approximation.

######  Corollary 5.2.

After changing the objective evaluated in line 3 to the partial objective,
Algorithm 1 gives a $3$-approximation for MobileVaccClinic with outliers.

Fairness: Although our algorithm provides an upper bound guarantee for the
maximum distance to a facility, the facility placement may significantly
differ between individuals, with some having a facility right next to them,
while others need to travel the whole $3R^{*}$ guarantee. Luckily, the vaccine
centers can vary from week to week or even day to day. Thus, we can use a
randomized algorithm such as the one given in Harris et al. (2020), to
guarantee that the re-provisioning of facilities over the course of many tries
will provide an improved per-point guarantee on expectation. Hence, we treat
the clients stochastically fairly.

######  Corollary 5.3.

When using the algorithm from Harris et al. (2020) instead of a simple
$k$-supplier algorithm, Algorithm 1 is able to output a distribution
$\tilde{\Omega}$ such that $\forall p\in\mathcal{P}$, we have
$\mathbb{E}_{F\sim\tilde{\Omega}}[d(S_{p},F)]\leq(1+2/e)R^{*}$ and
$\Pr[d(S_{p},F)\leq 3R^{*}]=1$.

Capacity: In this case, we assume that each facility we use has a capacity
$L$, i.e., at most $L$ clients can be assigned to it in any solution. Once
again, the FPT process we described earlier suffices to solve the problem.
Specifically, we can think of the set $A$ as the locations through which the
clients of $\mathcal{P}$ receive their service. Hence, as we did for the
regular case, we can create an instance of $k$-supplier where the points
requiring service are those of $A$, but this time each location of
$\mathcal{Y}$ for the $k$-supplier instance will have a capacity $L$. In other
words, this will be an instance of capacitated $k$-supplier. Furthermore, the
optimal solution of capacitated MobileVaccClinic will be a solution of the
same value for the capacitated $k$-supplier instance. Finally, it is also
trivial to see that any $\rho$-approximate solution for capacitated
$k$-supplier instance, will yield a $\rho$-approximate solution to capacitated
MobileVaccClinic.

######  Corollary 5.4.

When using the algorithm from An et al. (2013) instead of a simple
$k$-supplier algorithm, Algorithm 1 is an 11-approximation for capacitated
MobileVaccClinic.

###  5.2 Covering Algorithm

In Corollary 4.5, we show that any bicriteria algorithm needs to open at least
$k\ln n$ facilities in order to give a bounded approximation guarantee. Here,
we show that this is essentially tight: we give an algorithm that outputs a
set of locations of size at most $k(\ln n+1)$, while guaranteeing that the
objective value is at most that of an optimal solution.

Consider the related problem, which we call ClientCover, in which instead of
optimizing the radius $R$ given a budget $k$, we are given a target radius $R$
and want to choose a set $F\subseteq\mathcal{S}$ which minimizes $|F|$ and
guarantees that $d(S_{p},F)\leq R$ for each $p\in\mathcal{P}$. Notice that
this is just a standard Set Cover problem, where the sets are
$\\{p\in\mathcal{P}:d(S_{p},j)\leq R\\}$ for each $j\in\mathcal{S}$ and the
universe consists of the clients $\mathcal{P}$. Using a known greedy algorithm
for Set Cover (Williamson and Shmoys, 2011), we have an $H_{n}$-approximation
algorithm for ClientCover, where $H_{n}\leq\ln{n}+1$ is the $n$-th harmonic
number.

For generality, we will show how any $\alpha$-approximation algorithm for Set
Cover yields an $(1,\alpha)$-bicriteria algorithm for MobileVaccClinic via a
reduction to ClientCover. First, note that the optimal radius $R^{*}$ for an
instance of MobileVaccClinic is always the distance between some
$j\in\mathcal{C}$ and some $i\in\mathcal{S}$. Hence, there are at most
polynomially many options for it, specifically
$|\mathcal{C}|\cdot|\mathcal{S}|$. For each such option $R$, we create the
corresponding instance of ClientCover and run the set cover algorithm on it.
The final guarantees follow from the iteration when $R=R^{*}$. Observe at this
point that we can speed up the whole process by performing a binary search in
order to find $R^{*}$, and thus avoid the previously described exhaustive
search.

Algorithm 2 ClientCover Search

1: Binary search on the sorted list
$\\{d(i,j):j\in\mathcal{C},i\in\mathcal{S}\\}$, and let the current guess be
$R$:

2: Use $R$ to create the proper instance of ClientCover.

3: Obtain $\alpha$-approximate solution $F_{R}$ for that instance.

4: If $|F_{R}|>\alpha\cdot k$, increase $R$; else, decrease $R$.

5: Output $F_{R}$ for the minimum $R$ such that $F_{R}\leq\alpha\cdot k$.

######  Theorem 5.5.

Given an $\alpha$-approximation algorithm for set cover, Algorithm 2 gives an
$(1,\alpha)$-bicriteria algorithm for MobileVaccClinic.

###### Proof.

Let $R^{*}$ be the objective value for the optimal solution
$F^{*}\subseteq\mathcal{S}$ of MobileVaccClinic, where $|F^{*}|\leq k$. We
wish to show that ClientCover Search finds a radius $R$ in the list such that
$|F_{R}|\leq\alpha\cdot k$ and $R\leq R^{*}$. Consider an iteration of the
binary search where the radius guess is $R$. Suppose $R\geq R^{*}$; then there
must exist a solution of ClientCover of size at most $k$. The
$\alpha$-approximation algorithm will therefore output a set $F_{R}$ with
$F_{R}\leq\alpha\cdot k$ and $R$ will decrease. If $R<R^{*}$, then we either
find a solution with $F_{R}\leq\alpha\cdot k$, or we increase $R$ and move
closer to $R^{*}$. Finally, since $R^{*}$ is in the list
$\\{d(i,j):j\in\mathcal{C},i\in\mathcal{S}\\}$, the binary search necessarily
finds some $R\leq R^{*}$ with $F_{R}\leq\alpha\cdot k$. ∎

As in the case of our FPT algorithm, we can easily extend Algorithm 2 in order
to accommodate different settings. The only difference here lies at step 3,
where instead of a classic Set Cover algorithm we can run a different
algorithm.

Outliers: In order to modify our algorithm to only consider some fraction
$q\in(0,1)$ of the population, we can use some $\alpha$-approximation
algorithm for the Partial Set Cover problem, where the goal is to cover at
least a $q$-fraction of the universe elements. Hence, we naturally consider a
variant of ClientCover, which we call Partial ClientCover, that requires only
$\left\lfloor qn\right\rfloor$ points to be covered by balls of radius
$R^{*}$. Trivially, Partial ClientCover is a special case of Partial Set
Cover. Then the approach we described previously remains the same: we can
guess the optimal radius $R^{*}$ and obtain an $\alpha$-approximate solution
$F_{R^{*}}$ for the corresponding Partial ClientCover instance. This solution
will be optimal for MobileVaccClinic with outliers, while placing at most
$\alpha k$ facilities. In particular, we have the following corollary of
Theorem 5.5.

######  Corollary 5.6.

When using the greedy algorithm for Partial Set Cover (Williamson and Shmoys,
2011), Algorithm 2 gives a $(1,H_{\lfloor qn\rfloor})$-bicriteria algorithm
for MobileVaccClinic with outliers.

Fairness: When solving MobileVaccClinic with outliers, the algorithm may view
some demographic groups as outliers more often than others. To mitigate such
possibilities, we can use an algorithm for the Partition Set Cover problem
(Inamdar and Varadarajan, 2018) to guarantee that a large proportion of each
demographic group gets coverage. For example, we can guarantee that the
algorithm considers a proportional number of people from each (demographic)
group when choosing the vaccine center locations. The following approximation
guarantee will then follow directly from Inamdar and Varadarajan (2018) and
the outliers reduction before. We remark that even though the word “partition”
is in the name of the problem, the results of Inamdar and Varadarajan (2018)
extend to the case of overlapping demographic classes.

######  Corollary 5.7.

Let $C_{t}\subseteq\mathcal{P}$ for $t\in[r]$ be (not necessarily disjoint)
demographic classes and let $0\leq p_{t}\leq|C_{t}|$ be the coverage
requirements for each class. Using the algorithm of Inamdar and Varadarajan
(2018) at step 3, Algorithm 2 gives a $(1,O(\log n)+\log{r})$-bicriteria
algorithm while satisfying the coverage constraints.

Capacity: As before, we assume that each facility we use has capacity $L$. We
see that our general framework is still applicable: we can modify our
algorithm to satisfy these capacity constraints by replacing the Set Cover
algorithm with a Capacitated Set Cover algorithm when solving the ClientCover
problem. In fact, Wolsey (1981) shows that a greedy algorithm (similar to the
one for Set Cover) still gives a $\log{n}+1$ approximation algorithm in this
more general case.

######  Corollary 5.8.

When using the greedy algorithm given in Wolsey (1981) for Capacitated Set
Cover, Algorithm 2 gives a $(1,\log{n}+1)$-bicriteria algorithm while
satisfying the capacity constraints.

Budget: In the previous algorithms which solve ClientCover as a subroutine, we
violate the budget constraint $k$ by a non-trivial multiplicative factor,
which is a practical consideration that needs to be addressed. Luckily, it has
been shown that the greedy algorithm and other heuristics for Set Cover have
very small approximation ratios in practice (Lan et al., 2007). In fact, many
real-life instances of Set Cover are solved optimally or near optimally by the
greedy algorithm (Grossman and Wool, 1997). Given this empirical result (which
we also validate for our instances of the Set Cover problem), we get
$\alpha=1$ in our experiments when running ClientCover Search. In particular,
if we solve the ClientCover problem using a commercial mixed-integer linear
program (MILP) solver (Gurobi Optimization, 2021; Perron and Furnon, ), we can
solve the original problem to optimality. We emphasize that this is a non-
trivial contribution: directly formulating MobileVaccClinic as an MILP
requires $\Theta(n^{3})$ constraints, and we cannot even initialize the solver
using or-tools. In contrast, the Set Cover MILP only has $\Theta(n)$
constraints, which can be solved efficiently using commercial solvers. Hence,
when using an MILP solver to solve ClientCover, Algorithm 2 yields a practical
solution for solving MobileVaccClinic optimally.

Table 1: Network Information |  | Activity | Residential | Maximum | Measured
---|---|---|---|---|---
| Clients | Locations | Locations | Activity | Diameter (km)
Charlottesville City | 33156 | 5660 | 10038 | 9952 | 8.12
Albemarle County | 74253 | 9619 | 32981 | 24506 | 61.62

##  6 Experiments

###  6.1 Experimental Setup

Data: We run our experiments on the mobility data from Charlottesville City
and Albemarle County in Virginia. For these counties, we use synthetic data
constructed from the 2019 U.S. population pipeline (see Chen et al. (2021);
Machi et al. (2021) for details). This dataset was constructed by tracking the
week-long activity of county residents. Each resident is represented by a
sequence of activities, where each activity is described by duration, type,
and location in the county. The locations are given in geodetic coordinates
and are categorized as either a residential or activity (non-residential)
location. From this dataset, we can extract the locations visited by
individuals residing in the county and set all activity locations as potential
facility placements. A summary of the dataset is given in Table 1.

Baselines: We compare our algorithms with two heuristics: HomeCenters and
MostActive. In MostActive, we open vaccination centers at the $k$ most visited
locations. We set MostActive as the baseline because it is related to the
current heuristic used by the Virginia Department of Health. In HomeCenters,
we run k-supplier to place facilities at locations that minimize the maximum
distance from client homes. We compare with this baseline in order to show the
importance of considering mobility when placing the vaccination centers.

Objective: Recall that our objective is to minimize the maximum distance any
client needs to deviate from their path to reach some facility. Since our
location data is given in the geographic coordinate system, we approximate the
Earth as a sphere and use geodesic distance as our metric. In Section 6.2, we
notice that there is a sharp drop in the objective value if we only consider
$99\%$ of the population. As a result, we also evaluate the objective value of
our algorithms when $5\%$ of the people are considered outliers.

FPT details: When using FPT in our experiments, we pick $u=15$ locations that
cover the largest portion of the population (as given by the greedy algorithm
for the Maximum Coverage problem). We then run FPT using only knowledge of
these $u$ locations. The locations chosen are all popular public activity
locations, so we have a limited amount of privacy violation. As a result, the
performance of FPT is weaker on the full objective, but remains strong on
partial coverage (the outliers formulation). It is important to note that even
though we limit the knowledge of client-visited locations, FPT can still
choose to place facilities at any activity location in the dataset. For more
details on our implementation of FPT and the experiments, see our
GitHub111https://github.com/Ann924/MobileFacility.

###  6.2 Client Coverage Performance Analysis

First, we directly compare the performances between our algorithms and the
baselines. Because our objective value is defined by the maximum distance any
client must travel to reach their closest facility, it does not yield insight
into the distribution of travel distances. For this reason, we also assess how
large the radius around our placed facilities must be to cover proportion $p$
of the clients, for $p\in[0.8,1.0]$.

As seen in Figure 2, facility placements from HomeCenters and ClientCover
result in better full objectives while facility placements from MostActive and
FPT result in better partial objectives. This has a simple explanation: the
former two algorithms are forced to consider outliers since they directly
optimize the objective while the latter two optimize over only a portion of
the population. As a result, it makes sense to compare HomeCenters with
ClientCover and FPT with MostActive. We note that if we instead used the
outliers version of HomeCenters and ClientCover, we may not see this
disparity.

Figure 2: Charlottesville k = 10, Albemarle k = 20

###  6.3 Tradeoff between Radius and Budget

In addition to evaluating the performance of our algorithms at the current
budget, it is important to evaluate the sensitivity of our algorithms to an
increase in budget. That is, we want to know how much the objective value
would decrease if the county allocated more resources to deploy a greater
number of mobile facilities. This knowledge can influence policy decisions:
when an increase budget yields a sharp decrease in objective (rather than a
small decrease), the government has more incentive to fund another vaccination
center.

As seen in Figure 3, there is generally a sharp decrease in the objective
value when the budget is less than 6 for Charlottesville and 9 for Albemarle.
As the budget increases past those thresholds, the marginal returns become so
diminished that increasing the budget hardly changes the objective value. This
is especially prominent in the full objective performance of FPT and
MostActive. Hence, it is natural to recommend budgets of 6 and 9 to the
Charlottesville and Albemarle government, respectively.

Additionally, we wish to bring attention to the overall poor performance of
HomeCenters in these experiments. It is consistently outperformed by
ClientCover for the full objective and is outperformed by every algorithm when
evaluating the partial objective. Furthermore, HomeCenters does not exhibit
strong budget sensitivity when assessing the 95% objective. On Albemarle, its
performance plateaus around an objective value of 1.5km, which is more than 3
times larger than the objective of the algorithms that consider mobility
patterns.

A seemingly weird result from the experiment is the tradeoff curve for
HomeCenters. Though there is a general downward trend in the objective value
as the budget increases, there are cases in each county where increasing the
budget results in an increase in the objective value. This contradictory
phenomenon is caused by the limited correlation between the distance to homes
and our objective; as a result, noise/luck has a considerable effect. The
noisiness of HomeCenters emphasizes the importance of our work of modeling
mobile populations.

Figure 3: Tradeoff between maximum distance needed to travel and the number of
vaccine centers placed. In the 95% objective plots, ClientCover and
HomeCenters are run with the outliers formulation.

###  6.4 A Kernel Property

Through our experiments, we notice a nice (empirical) property of the vaccine
center locations selected by some of our algorithms. Imagine a case where we
(the government) have the funds to place five mobile vaccine centers and we
use our algorithms to pick the five locations to place them. Then, after two
weeks, the government decides that the disease is causing too much economic
devastation and, in turn, funds three more mobile vaccine centers. When we ask
our algorithms to place the eight vaccine centers (approximately) optimally,
it turns out that the eight chosen locations will often contain the original
five chosen locations as a subset. The original five locations are then called
a kernel.

Table 2: Charlottesville Kernel | Home |  | Client
---|---|---|---
| Center | FPT | Cover
3 $\rightarrow$ 4 | 1 | 1 | 2
4 $\rightarrow$ 5 | 1 | 1 | 4
5 $\rightarrow$ 6 | 2 | 1 | 5
6 $\rightarrow$ 7 | 3 | 2 | 6
7 $\rightarrow$ 8 | 0 | 1 | 6
8 $\rightarrow$ 9 | 0 | 1 | 5
9 $\rightarrow$ 10 | 3 | 1 | 7

Table 3: Albemarle Kernel | Home |  | Client
---|---|---|---
| Center | FPT | Cover
6 $\rightarrow$ 7 | 3 | 1 | 5
7 $\rightarrow$ 8 | 1 | 2 | 6
8 $\rightarrow$ 9 | 1 | 2 | 7
9 $\rightarrow$ 10 | 1 | 1 | 8
10 $\rightarrow$ 11 | 4 | 3 | 5
11 $\rightarrow$ 12 | 1 | 1 | 9
12 $\rightarrow$ 13 | 1 | 2 | 4
13 $\rightarrow$ 14 | 0 | 1 | 7
14 $\rightarrow$ 15 | 5 | 2 | 7
15 $\rightarrow$ 16 | 1 | 1 | 8
16 $\rightarrow$ 17 | 1 | 0 | 12
17 $\rightarrow$ 18 | 7 | 0 | 8
18 $\rightarrow$ 19 | 2 | 1 | 7
19 $\rightarrow$ 20 | 0 | 0 | 8

The ideal kernel property occurs when every set of chosen facilities of size
$k$ contains the set of chosen facilities for budget $k-1$. In order to
determine the presence of a kernel for each of our algorithms, we calculate
the number of facilities chosen with budget $k-1$ that are not also chosen
with budget $k$. These values populate Tables 3 and 3, where the leftmost
column denotes the budgets compared. By definition, MostActive has the kernel
property since it is a greedy algorithm. Our FPT algorithm also
(approximately) satisfies the kernel property while maintaining a stronger
performance than the baseline. The remaining two algorithms do not exhibit the
property: both HomeCenters and ClientCover pick (almost) completely different
locations upon increasing the budget. Because they require less relocation,
MostActive and FPT have advantageous properties when the budget is adaptive
and vaccine distribution is time-consuming.

We recognize that this is not necessarily applicable to our experimental
setting, COVID-19, since transportation of vaccines is (relatively) easy in
Virginia. However, for the Ebola outbreak in 2014, the kernel property was
recognized as an important property to have since vaccine distribution was a
much more costly process. Furthermore, we note that our algorithms are not
explicitly designed to have this property; it is only empirically verified.

###  6.5 Information Constraints

In our previous experiments with ClientCover, we assumed that we had full
knowledge of the locations each person visited throughout a day. Next, we wish
to understand how fine-grained this data needs to be in order for ClientCover
to outperform our other algorithms; this also addresses privacy concerns
raised when using the exact mobility data of individuals. In order to model
loss of detailed movement patterns, we cluster the locations within a given
radius $r$ together and apply ClientCover on the resulting cluster centers.
The details of the clustering algorithm can be found in our code, but the
general idea is to define each location to be a potential cluster center and
then use the greedy set cover algorithm to pick a minimum set of clusters
centers that cover all original locations with radius $r$. Using this general
method for both Charlottesville and Albemarle, we vary the radius $r$ between
100 and 600 meters to see how much privacy ClientCover can give while still
maintaining a superior performance over the baselines.

Figure 4: ClientCover performance under loss of information (k=20)

As we see in Figure 4, clustering with a radius of 0.1–0.48 km on
Charlottesville leads to a gradual increase in the objective value. At a
clustering radius of 0.48 km, there is a sharp increase in the objective value
from the facility placements. Even so, the resulting objective value is
significantly smaller than 2.39 km, as obtained by MostActive, and 2.52 km, as
obtained for FPT. Furthermore, even by clustering the data with a radius that
is 45% of the original ClientCover objective, we can still perform better than
the HomeCenters baseline. A similar trend occurs for Albemarle where the
change in the solution value is relatively small when clustering from
0.10–0.35 km but grows rapidly after 0.35 km. We conclude that even when
giving some privacy to individuals, ClientCover still performs much better
than FPT, MostActive, and HomeCenters.

##  7 Conclusion

Here, we introduce a generalization of the classical $k$-supplier problem
where we consider the mobility of populations when placing facilities. We show
that designing an approximation algorithm for this variant is NP-Hard, so we
turn to fixed-parameter tractability and bicriteria approximation algorithms
to get around our hardness result. Finally, we experimentally show the
efficacy of our algorithms in comparison to natural baselines. Since we have
demonstrated the importance of modeling mobile populations, a natural next
step is to extend other variants of the facility location problem to this
setting as well.

Acknowledgements: We express our sincere thanks to the AAMAS referees for
suggesting the experiments in Section 6.5 and the extension of capacity
constraints. We also thank members of the Biocomplexity COVID-19 Response Team
and the Network Systems Science and Advanced Computing (NSSAC) Division for
their thoughtful comments and suggestions related to epidemic modeling and
response support. George Li, Aravind Srinivasan, and Leonidas Tsepenekas were
supported in part by NSF award number CCF-1918749. Ann Li, Madhav Marathe, and
Anil Vullikanti were supported by DTRA (Contract HDTRA1-19-D-0007), University
of Virginia Strategic Investment Fund award number SIF160, National Institutes
of Health (NIH) Grants 1R01GM109718, 2R01GM109718, OAC-1916805 (CINES),
CCF-1918656 (Expeditions), CNS-2028004 (RAPID), OAC-2027541 (RAPID),
IIS-1908530, IIS-1955797, and IIS-2027848. The U.S. Government is authorized
to reproduce and distribute reprints for Governmental purposes notwithstanding
any copyright annotation thereon.

## References

  * Afshari and Peng [2014] Hamid Afshari and Qingjin Peng.  Challenges and solutions for location of healthcare facilities.  Industrial Engineering and Management, 2014:1–12, 2014.
  * An et al. [2013] Hyung-Chan An, Aditya Bhaskara, and Ola Svensson.  Centrality of trees for capacitated k-center, 2013.
  * Anegg et al. [2020] Georg Anegg, Haris Angelidakis, Adam Kurpisz, and Rico Zenklusen.  A technique for obtaining true approximations for $k$-center with covering constraints, 2020.
  * Celik Turkoglu and Erol Genevois [2020] Derya Celik Turkoglu and Mujde Erol Genevois.  A comparative survey of service facility location problems.  Annals of Operations Research, 292, 09 2020.
  * Chang et al. [2021] Serina Chang, Emma Pierson, Pang Wei Koh, Jaline Gerardin, Beth Redbird, David Grusky, and Jure Leskovec.  Mobility network models of covid-19 explain inequities and inform reopening.  Nature, 589, 2021.
  * Chen et al. [2021] Jiangzhuo Chen, Stefan Hoops, Achla Marathe, Henning Mortveit, Bryan Lewis, Srinivasan Venkatramanan, Arash Haddadan, Parantapa Bhattacharya, Abhijin Adiga, Anil Vullikanti, Aravind Srinivasan, Mandy L Wilson, Gal Ehrlich, Maier Fenster, Stephen Eubank, Christopher Barrett, and Madhav Marathe.  Prioritizing allocation of covid-19 vaccines based on social contacts increases vaccination effectiveness, 2021.
  * de Vries et al. [2020] Harwin de Vries, Joris van de Klundert, and Albert P.M. Wagelmans.  The roadside healthcare facility location problem a managerial network design challenge.  Production and Operations Management, 29(5):1165–1187, 2020.
  * Drezner [1995] Zvi Drezner.  Facility location: A survey of applications and methods, 1995.
  * Gonzalez [1985] Teofilo F Gonzalez.  Clustering to minimize the maximum intercluster distance.  Theoretical computer science, 38:293–306, 1985.
  * Grossman and Wool [1997] Tal Grossman and Avishai Wool.  Computational experience with approximation algorithms for the set covering problem.  European Journal of Operational Research, 101(1):81–92, 1997.
  * Gurobi Optimization [2021] LLC Gurobi Optimization.  Gurobi optimizer reference manual, 2021.
  * Harris et al. [2020] David G. Harris, Shi Li, Thomas Pensyl, Aravind Srinivasan, and Khoa Trinh.  Approximation algorithms for stochastic clustering, 2020.
  * Hochbaum and Shmoys [1985] Dorit S. Hochbaum and David B. Shmoys.  A best possible heuristic for the k-center problem.  Math. Oper. Res., 10(2), May 1985.
  * Hochbaum and Shmoys [1986] Dorit S. Hochbaum and David B. Shmoys.  A unified approach to approximation algorithms for bottleneck problems.  J. ACM, 33(3), May 1986.
  * Inamdar and Varadarajan [2018] Tanmay Inamdar and Kasturi Varadarajan.  On the partition set cover problem, 2018.
  * Lan et al. [2007] Guanghui Lan, Gail W. DePuy, and Gary E. Whitehouse.  An effective and simple heuristic for the set covering problem.  European Journal of Operational Research, 176, 2007.
  * Machi et al. [2021] Dustin Machi, Parantapa Bhattacharya, Stefan Hoops, Jiangzhuo Chen, Henning Mortveit, Srinivasan Venkatramanan, Bryan Lewis Mandy Wilson, Arindam Fadikar, Tom Maiden, Christopher L. Barrett, and Madhav V. Marathe.  Scalable epidemiological workflows to support covid-19 planning and response, 2021.
  * Meskarian et al. [2017] R. Meskarian, M. L. Penn, S. Williams, and T. A Monks.  Facility location model for analysis of current and future demand for sexual health services.  PLoS One, 2017.
  * Nedjati and Valipour [2012] Arman Nedjati and Mohammad Valipour.  Solving health care facility location problems with new heuristic algorithm method.  International journal of modeling and optimization, 12 2012.
  * nyt [a] “meet the four kinds of people holding us back from full vaccination”, new york times.  https://www.nytimes.com/interactive/2021/05/18/opinion/covid-19-vaccine-hesitancy.html.
  * nyt [b] See How Vaccinations Are Going in Your County and State.  https://www.nytimes.com/interactive/2020/us/covid-19-vaccine-doses.html.  Accessed: 2022-01-23.
  * [22] Laurent Perron and Vincent Furnon.  Or-tools.
  * Tai et al. [2020] Don Bambino Geno Tai, Aditya Shah, Chyke A. Doubeni, Irene G. Sia, and Mark L. Wieland.  The disproportionate impact of covid-19 on racial and ethnic minorities in the united states.  Infectious Diseases Society of America, 2020.
  * [24] COVID-19 Surveilance Dashboard.  https://nssac.bii.virginia.edu/covid-19/dashboard/.  Accessed: 2022-01-23.
  * Wang et al. [2020] Lijing Wang, Xue Ben, Anniruddha Adiga, Adam Sadilek, Ashish Tendulkar, Srinivasan Ventatramanan, Anil Vullikanti, Milind Tambe, and Madhav Marathe.  Using mobility data to understand and forecast covid 19 dynamics.  medRxiv, 2020.
  * Williamson and Shmoys [2011] David P. Williamson and David B. Shmoys.  The Design of Approximation Algorithms.  Cambridge University Press, USA, 1st edition, 2011.
  * Wolsey [1981] Laurence A. Wolsey.  An analysis of the greedy algorithm for the submodular set covering problem.  Combinatorica, 2, 1981.