# Fluorescence quenching in plasmonic dimers due to electron tunneling

Henrikh M. Baghramyan1 Cristian Ciracì1<EMAIL_ADDRESS>1Center for
Biomolecular Nanotechnologies, Istituto Italiano di Tecnologia, Via Barsanti
14, 73010 Arnesano (LE), Italy

(August 27, 2024)

###### Abstract

Plasmonic nanoparticles provide an ideal environment for the enhancement of
fluorescent emission. On the one hand, they locally amplify the
electromagnetic fields, increasing the emitter excitation rate, and on the
other hand, they provide a high local density of states that accelerates
spontaneous emission. However, when the emitter is placed in close proximity
to a single metal nanoparticle, the number of nonradiative states increases
dramatically, causing the fluorescence to quench. It has been predicted
theoretically that, through a judicious placing of the emitter, fluorescence
in plasmonic nanocavities can be increased at monotonically. In this article,
we show that such monotonic increase is due to the use of local response
approximation in the description of the plasmonic response of metal
nanoparticles. We demonstrate that taking into account the electron tunneling
and the nonlocality of the surrounding system via the quantum hydrodynamic
theory results eventually in a quenching of fluorescence enhancement also when
the emitter is placed in a nanocavity, as opposed to local response and
Thomas-Fermi hydrodynamic theory results. This outcome marks the importance of
considering the quantum effects, in particular, the electron tunneling to
correctly describe the emission effects in plasmonic systems at nanoscale.

Quantum hydrodynamic theory, Plasmonics, Surface plasmons, Fluorescence,
quantum emitters

##  I Introduction

A quantum emitter in an excited states can undergo a decay by emitting a
photon into one of the modes provided by the surrounding electromagnetic
environment [1, 2]. Due to plasmonic effects, metal nanoparticles (NPs) can
provide a high local density of electromagnetic states in their vicinity. The
emission of a quantum emitter than can be highly amplified at the proximity of
plasmonic NPs [3]. Similarly, the fields of an external electromagnetic wave
can get localized and enhanced at the metal NP surface well beyond the
diffraction limit. A quantum emitter placed in the vicinity of a metal NP then
can also be more easily excited [4]. These two processes – the excitation and
the emission – contribute to the fluorescence of the emitter [5]. Plasmonic-
enhanced fluorescence has several applications [6], spanning through
biomedicine [7, 8, 9], sensors [10, 11], quantum information [12] and
bionanotechnolgies [13]. In particular, a hybrid plasmonic nanoantenna has
been suggested in [14] for the efficient enhancement of emission
directionality and keeping high quantum efficiency. Also, more than 100 times
fluorescence enhancements have been observed in [15], which is critical for
readout by inexpensive point-of-care detectors.

At very short distances from a single metal NP, however, fluorescence quenches
since the emitter couples mostly to nonradiative modes [16, 17, 18, 19]. To
suppress plasmonic quenching, one could place the emitter in a nanocavity [20,
21, 22, 23, 24, 25]. In the cavity formed by an NP coupled to a metal film,
for example, quenching disappears because higher-order modes gain a radiative
nature while the excitation rate increases [22, 23]. Within the local response
approximation (LRA), it is predicted that, as the NP-film distance is reduced,
the fluorescence rate monotonically increases, and quenching never occurs [22,
23]. This, however, could be the result of a simplistic description rather
than a physical consequence. In fact, local theories do not account for the
nonlocal nature of the dielectric function of the metal at short distances
[26, 27, 28]. Indeed, the hydrodynamic theory (HT) for the plasmonic response
predicts bound to the maximum field enhancement achievable in nanocavities
[26]. This is obtained by considering the electron pressure term when
evaluating the permittivity of the metal. The pressure term originates from
the Thomas-Fermi (TF) kinetic energy (KE) [29, 30]. However, the TF-HT is
usually used with hard-wall boundary conditions at the metal surface, which
does not allow electrons to spill outside of the metal [31, 32]. In addition,
TF KE is only a local functional of the density, i.e. it depends on the value
of the electron density at the point considered.

In Ref. [33], the authors consider the generalized nonlocal optical response
(GNOR) by accounting for both, nonlocal response and surface-enhanced Landau
damping, and, in fact, predict a reduction of the fluorescence enhancement as
compared to LRA calculations. However, even using the GNOR description there
is no evidence of quenching. A more accurate approach would require taking
into account the KE dependence on higher-order derivatives of the density. The
quantum hydrodynamic theory (QHT) adds such correction to the TF contribution,
namely the von Weizäcker (vW) KE term [34], which is equivalent to the quantum
potential [35, 36] and depends on the gradient of the electron density. In the
framework of TF and vW KE functionals, the spatial dependence of the electron
density in the metal is naturally taken into account and spill-out, and
electron tunneling effects can be easily considered in the calculation of
plasmon energy [37, 38] while still retaining retardation effects [38].
Nonetheless, the QHT in the TF-vW approximation may result in unstable
resonances due to the propagation nature of the asymptotic solutions at high
energies [38, 39]. This behaviour can be amended by considering a Laplacian
level correction in the KE [40, 39].

In this article, we apply the QHT to evaluate the fluorescence enhancement of
the emitter in the gap of the Na jellium sphere dimer in the weak coupling
regime. We show that for small enough gaps ( $\approx 0.5$ nm) the QHT
predicts fluorescence quenching, as opposed to the LRA [22, 23] calculations
where the fluorescence enhancement monotonically increases as the gap closes.
Our formalism allows distinguishing the two quantum corrections that are often
discussed in the literature: namely, spatial dispersion (or nonlocality) and
electronic spill-out. In this paper, we show that the former does not produce
any significant quenching [41, 42, 33, 43], while the latter does. Therefore,
both mechanisms manifest in fluorescence in a very different way. Finally, we
analyze the impact of Laplacian-level KE correction to the fluorescence
spectra.

##  II Theory

To study the fluorescence process, the excitation of the emitter and its
subsequent decay should be considered. In the weak-coupling regime, the
fluorescence enhancement, $\eta_{em}$, is the product of the excitation
enhancement, $\eta_{ex}$, at the location of the emitter, $\mathbf{r}_{e}$,
and quantum yield, $q$, of the emitter that accounts for the probability of
the photons to couple to radiative states. The excitation and the decay
problems are solved independently, as schematically shown in Fig. 1.

Figure 1: Schematic image of the excitation and decay phases. First, the
emitter is excited by a plane wave polarized along the $z$-axis impinging on
the NP dimer (I); then, the emitter’s dipole field interacts with the
electromagnetic states provided by the dimer (II).

The first step then is to calculate the excitation enhancement of the incident
radiation by the plasmonic environment. We consider an emitter placed at the
center of a Na jellium sphere dimer at equal distances from the spheres.
Because the field in the dimer gap is mostly oriented along the $z$-direction,
we only consider the case of an emitter vertically oriented. Moreover, we
assume a single frequency approximation and set the excitation and emission
frequencies to be equal[18]. At frequency $\omega$ and at the emitter location
$\mathbf{r}_{e}$ $\eta_{ex}$ can be calculated by the following formula [4]:

$\eta_{ex}=\frac{\gamma_{ex}}{\gamma_{ex}^{0}}=\frac{\left|\left(\mathbf{E}_{\rm{I}}\left(\mathbf{r}_{e},\omega\right)+\mathbf{E}_{0}\left(\mathbf{r}_{e},\omega\right)\right)\cdot\hat{\mathbf{n}}_{e}\right|^{2}}{\left|\mathbf{E}_{0}\left(\mathbf{r}_{e},\omega\right)\cdot\hat{\mathbf{n}}_{e}\right|^{2}}$
(1)

with $\gamma_{\mathrm{ex}}^{0}$ and $\gamma_{\mathrm{ex}}^{0}$ being the
excitation rate in vacuum and in our system respectively,
$\hat{\mathbf{n}}_{e}=\hat{\mathbf{z}}$ the dipole moment unit vector,
$\mathbf{E}_{0}$ is the incident field (a plane wave polarized along the
$z$-axis), and $\mathbf{E_{\rm{I}}}$ is the field scattered by the dimer.

The scattered field $\mathbf{E}_{\rm{I}}$ at the location $\mathbf{r}_{e}$ is
obtained by solving the Maxwell equations coupled to the quantum hydrodynamic
equation for the free electron dynamics in the Na NPs [37, 38]:

$\displaystyle\nabla\times\nabla\times\mathbf{E_{\rm{I}}}-\frac{\omega^{2}}{c^{2}}\mathbf{E_{\rm{I}}}=\omega^{2}\mu_{0}\textrm{{P}},$
(2a) $\displaystyle\frac{en_{0}}{m_{e}}\nabla\left(\frac{\delta
G\left[n\right]}{\delta
n}\right)_{1}+\left(\omega^{2}+i\gamma\omega\right)\textrm{{P}}=-\varepsilon_{0}\omega_{p}^{2}\left(\mathbf{E_{\rm{I}}}+\mathbf{E}_{0}\right),$
(2b)

with $\varepsilon_{0}$ and $\mu_{0}$ are respectively the vacuum permittivity
and permeability, P is polarization vector, $c$ is the speed of light in free
space, $e$ is the absolute value of the electron charge and $m_{e}$ is its
mass, $\gamma$ is the damping rate, and
$\omega_{p}\left(\textbf{r}\right)=\sqrt{e^{2}n_{0}\left(\textbf{r}\right)/\left(m_{e}\epsilon_{0}\right)}$
is the plasma frequency where $n_{0}\left(\textbf{r}\right)$ is the ground-
state electron density. For simplicity, in this work, we neglect nonlocal
damping [44, 45, 46, 33] and set the damping rate as
$\gamma=0.066\leavevmode\nobreak\ \textrm{eV}+v_{\textrm{F}}/R$ where $v_{F}$
is the Fermi velocity for Na and $R$ is the radius of the spheres. The
functional $G\left[n\right]$ accounts for the quantum nature of the electron
gas total internal energy. It is the sum of noninteracting KE functional,
$T_{s}\left[n\right]$, and exchange-correlation (XC) functional,
$E_{\textrm{XC}}^{\textrm{LDA}}\left[n\right]$, that is
$G\left[n\right]=T_{s}\left[n\right]+E_{\textrm{XC}}^{\textrm{LDA}}\left[n\right]$.
$T_{s}\left[n\right]$ can be approximated considering the contributions of the
TF, vW, Pauli-Gaussian second-order (PGS) and Laplacian-level functionals, as:

$T_{s}\left[n\right]=\int\tau\left(n,w,q\right)\textit{d}^{3}\textrm{{r}},$
(3)

where $w=\nabla n\cdot\nabla n$ and $q=\nabla^{2}n$, and

$\tau=\tau^{\mathrm{vW}}+\tau^{\mathrm{PGS}}+\tau^{\mathrm{TF}}\left[\beta
q_{r}^{2}+2\beta q_{0}^{2}\ln\left(1+q_{r}/q_{0}\right)\right]$ (4)

In Exp. (4), each $\tau$ denotes the KE density associated with each
functional [39].

$\displaystyle{\color[rgb]{0,0,0}T_{\rm TF}[n]=\int\tau^{\rm
TF}\left(n\right)\textit{d}^{3}\textrm{{r}}=\int\left(E_{h}a_{0}^{2}\right)\frac{3}{10}\left(3\pi^{2}\right)^{2/3}n^{5/3}\textit{d}^{3}\textrm{{r}},}$
(5a) $\displaystyle{\color[rgb]{0,0,0}T_{\rm vW}[n,w]=\int\tau^{\rm
vW}\left(n,w\right)\textit{d}^{3}\textrm{{r}}=\int\left(E_{h}a_{0}^{2}\right)\frac{w}{8n}\textit{d}^{3}\textrm{{r}},}$
(5b) $\displaystyle{\color[rgb]{0,0,0}T_{\rm
PG\alpha}[n,w]=\int\tau^{\textrm{PG}\alpha}\left(n,w\right)\textit{d}^{3}\textrm{{r}}=\int\tau^{\rm
TF}\left(n\right)e^{-\alpha Cn^{-8/3}w}\textit{d}^{3}\textrm{{r}},}$ (5c)

where $E_{h}$ is the Hartee energy, $a_{0}$ is the Bohr radius, $\alpha=40/27$
for PGS, $\beta$ and $q_{0}$ are the free parameters of the Laplacian-level
term and can be set as $\beta=0.25$ and $q_{0}=700$ as in [39], and
$q_{r}=3\nabla^{2}n/(40\tau^{TF})$ is the reduced Laplacian. The expressions
of the potentials in $\left(\delta G\left[n\right]/\delta
n\right)_{1}=\left(\delta E_{\mathrm{XC}}^{\mathrm{LDA}}/\delta
n\right)_{1}+\left(\delta T_{s}/\delta n\right)_{1}$ can be found in the
Section 2 in the SM.

For the equilibrium density, we use the following model expression:

$n_{0}\left(\rho,z\right)=\frac{f_{0}}{1+e^{\kappa\left[r_{\uparrow}\left(\rho,z\right)-R\right]}}+\frac{f_{0}}{1+e^{\kappa\left[r_{\downarrow}\left(\rho,z\right)-R\right]}},$
(6)

where
$r_{\uparrow}\left(\rho,z\right)=\sqrt{\rho^{2}+\left(z-g/2-R\right)^{2}}$ and
$r_{\downarrow}\left(\rho,z\right)=\sqrt{\rho^{2}+\left(z+g/2+R\right)^{2}}$
are for the upper and lower spheres respectively, and $g$ is the vertical
separation measured from the surfaces of the spheres (see Fig. 1). Since we
assume a jellium description, the surface of the ion background is located at
a distance of $r_{s}$, the Wigner-Seitz radius, from the last atom such that a
$g$ gap would correspond to a distance of $2r_{s}+g$ between the closest
atoms. In Exp. (6) $\kappa=1.05/a_{0}$ is fitted to the asymptotic decay of
Kohn-Sham densities [38], with $a_{0}$ being the Bohr radius and $f_{0}$ is
the normalization constant that is found from the condition:

$\int_{\Omega}n_{0}dV=2\times\frac{4}{3}\pi R^{3}n^{+}$ (7)

where $n^{+}=\left(4/3\pi r_{s}^{3}\right)^{-1}/a_{0}^{3}$ denotes the density
of positive charges with $r_{s}=4$ a.u. Wigner-Seitz radius for Na. The
expressions for the KE functionals and linear potentials associated with them
can be found in [39, 38].

After calculating the plasmonic response of the dimer to the incident
radiation, and hence the field that excites the emitter, we proceed to the
second step that is the evaluation of the decay rates.

The rate of spontaneous emission is a sum of the radiative and nonradiative
decays, i. e. $\gamma_{\rm sp}=\gamma_{\rm r}+\gamma_{\rm nr}$, and can be
found from the relation:

$\gamma_{\mathrm{sp}}=\frac{2\pi^{2}c}{\hbar\varepsilon_{0}}\rho_{0}\left[{\mathbf{p}}\cdot\operatorname{Im}\\{\mathbf{G}(\mathbf{r}_{e},\mathbf{r}_{e})\\}\cdot{\mathbf{p}}\right],$
(8)

where $\hbar$ is the reduced Planck constant,
$\rho_{0}=\frac{\omega^{2}}{\pi^{2}c^{3}}$ is the LDOS in vacuum, $\mathbf{p}$
is the transition dipole moment and $\bf G$ is the dyadic Green’s function of
the system [43].

The expression in Eq. (8) can be evaluated by using the relation:

$\mathbf{E_{\rm{II}}}(\mathbf{r},\omega)=\mu_{0}\omega^{2}\mathbf{G}\left(\mathbf{r},\mathbf{r}_{e}\right)\cdot\mathbf{p}_{c},$
(9)

where $\mathbf{p}_{c}=2\mathbf{p}$ is the classical dipole moment set to
$\left|\mathbf{p}_{c}\right|=1\rm{D}$ [47]. The field $\mathbf{E_{\rm{II}}}$
is found by solving a system of equations similar to Eqs. $\eqref{fq01}$,
where $\mathbf{E}_{0}$ is substituted by the field associated with a
$z$-oriented dipole radiating in free-space [43].

The radiative decay rate, $\gamma_{\mathrm{r}}$ can be found by evaluating the
energy flow radiating out from the system:

$\gamma_{r}=\frac{1}{2}\frac{\gamma^{0}}{W^{0}}\int_{\partial\Omega}\operatorname{Re}\left\\{\mathbf{E}_{\rm{II}}\times\mathbf{H}^{*}_{\rm{II}}\right\\}\mathrm{d}A$
(10)

where
$\gamma^{0}=\omega^{3}|\mathbf{p}|^{2}/\left(3\hbar\pi\varepsilon_{0}c^{3}\right)$
is the radiative decay rate and
$W^{0}=\omega^{4}|\mathbf{p}_{c}|^{2}/\left(12\pi\varepsilon_{0}c^{3}\right)$
the radiated power in vacuum. Having found $\gamma_{\rm r}$, the nonradiative
part is obtained as
$\gamma_{\mathrm{nr}}=\gamma_{\mathrm{sp}}-\gamma_{\mathrm{r}}$, and the
quantum yield $q$ can be found via:

$q=\frac{\gamma_{\mathrm{r}}}{\gamma_{\mathrm{sp}}}.$ (11)

Finally, the fluorescence enhancement of the emitter in the plasmonic
environment of the dimer is given by:

$\eta_{em}=\eta_{ex}q.$ (12)

In this article, we will consider the following cases with increasing
complexity in the description of the electron dynamics in the Na spheres:

  1. I:

LRA, where the metal is described through the conventional Drude model,
corresponding to $G\left[n\right]=0$ in Eq. (2b);

  2. II:

TF-HT, in which we consider the TF pressure term in the hydrodynamic equation,
i.e. $G\left[n\right]=T_{\rm TF}[n]$;

  3. III:

QHT, where we account for the vW KE and XC corrections,
$G\left[n\right]=T_{\rm TF}[n]+T_{\rm vW}[n,w]+E_{\rm XC}[n]$, and consider a
spatial dependent equilibrium charge density;

  4. IV:

QHT-PGSL, in which the KE functional of the QHT is improved by adding the PGS
and Laplacian terms, i. e.,
$\tau=\tau^{\mathrm{vW}}+\tau^{\mathrm{PGS}}+\tau^{\mathrm{TF}}\beta
q_{r}^{2}$;

  5. V:

QHT-PGSLN, where we add the Laplacian-dependent logarithm correction to the
QHT-PGSL, i. e.,
$\tau=\tau^{\mathrm{vW}}+\tau^{\mathrm{PGS}}+\tau^{\mathrm{TF}}\left[\beta
q_{r}^{2}+2\beta q_{0}^{2}\ln\left(1+q_{r}/q_{0}\right)\right]$.

##  III Methods

The system of Eqs. (2), where ${\bf E}_{0}$ is either a plane wave or a dipole
filed, is solved using a finite element method through a customized
implementation into a commercial software. Moreover, to reduce the
computational cost, we take advantage of the symmetry of the system, by
expanding the fields in a series of cylindrical harmonics, each of which can
propagate independently [48, 49]. Following this method, it is possible to
reduce a three-dimensional calculation to few two-dimensional ones. In fact,
because the systems considered here are very small it is possible to limit the
maximum azimuthal number $|m_{max}|\leq 1$. Note that since the field
associated with a vertically oriented dipole is axially symmetric, the
emission properties can be exactly evaluated by only considering the azimuthal
number $m=0$ [46]. More detail on the FEM implementation is given in the
Appendix Section A.

##  IV Discussion of results

Let us start by comparing the enhancement of the incident field
$\mathbf{E}_{0}$ inside the dimer gap when different models are considered.

Figure 2: Norm of the total electric field (a-c) and the absolute value of the
induced change density (d) as obtained from the LRA, TF and QHT at
$\lambda=400$ nm for $g=0.1$ nm and $R=10$ nm. The induced charge density is
normalized to the maximum value in each scenario.

In Fig. 2, we show the distribution of the norm of the total field in the gap
region for the LRA, TF-HT, and QHT approaches. Note that the spatial range is
different for each model. For the LRA (Fig. 2 (a)), the field is extremely
localized and greatly enhanced compared to other methods. This happens because
the response of the system is purely local and the induced charges can
accumulate on an infinitesimally small layer at the metal surface confining
the field in a much tighter region, as can be seen from Fig. 2 (a). So, by
shrinking the gap, we get increasing values for the field and the induced
density, as can be seen in Fig. 2 (d). The situation is different for the TF-
HT as depicted in Fig. 2 (b). In this case, the induced charge spreads into
the bulk region since the optical response is no longer completely local. Note
that the field is localized into a volume of about 10 times wider than the
previous case. For the TF-HT, we used $n_{0}=n^{+}$ and imposed hard-wall
boundary conditions, i.e., $\mathbf{P}\cdot\hat{\mathbf{n}}=0$ [38], at the
metal surface, which forces the induced charge density
$n_{1}\left(\textrm{{r}}\right)$ to spread only inside the metal as shown in
Fig. 2 (d).

The most physical scenario is observed for the QHT in Fig. 2 (c). At optical
frequencies, the electron spill-out from the metal surface [50, 51] results in
the electron tunneling, which is the QHT case. This process is different to
the one discussed in [52] where the tunneling happens through the emitter as
the latter serves as a conductive linker . Both the equilibrium density
$n_{0}$ and the induced charge density $n_{1}$ are not constrained by the
surface of the metal: we use Eq. (6) for $n_{0}$ where the $s_{0}$ parameter
controls the spill-out, and hence the tunneling of the density out of the
metal and the hard-wall boundary condition on $n_{1}$ is lifted. As a result,
the electric field goes into the metal and moves out of the gap center
($r=0,z=0$) by having maximum values in the regions to the left and right of
the gap center. Indeed, the behavior of the field resembles the movement of
the liquid as it could be naturally foreseen considering the hydrodynamic
approach to the problem (see Eq. (2b)). The spill-out and tunneling affect the
current density $\mathbf{J}=\partial\mathbf{P}/\partial t=-i\omega\mathbf{P}$,
at the surface of the spheres for each scenario as it can be seen from Fig. S4
in SM.

Figure 3: The dependence of the excitation enhancement $\eta_{ex}$ and quantum
yield $q$ on the gap size ($0.25\leq g\leq 100$ nm) as calculated from the
LRA, TF-HT and QHT for different radii at $\lambda=400$ nm.

Next, let us analyze how the terms in Eq. (12), i.e., the excitation
enhancement, $\eta_{ex}$, and quantum yield, $q$, change by varying the size
of the gap for the different models, as shown in Fig. 3 for a fixed
wavelength. For large gaps, the electric field is localized only at the
surface of the NPs (especially for the smaller systems) and, hence,
$\eta_{ex}$ is relatively small and not much affected by the model used (the
differences are due mostly to the shift of the plasmon resonance).

Figure 4: The fluorescence enhancement for different radii, $g$ and $\lambda$
evaluated for the LRA, TF-HT and QHT. Figure 5: The dependence of the
fluorescence enhancement for $R=10$ nm gap size $g$ and the wavelength
$\lambda$ as computed from LRA, TF-HT and QHT.

As the gap size shrinks, $\eta_{ex}$ reaches very high values for the LRA and
TF-HT compared to the QHT as expected. In addition, for QHT, we observe a drop
at very small values of the gap since the field is mostly concentrated off the
location of the emitter (see Fig. 2 (c)). An opposite behavior is observed for
the quantum yield. The smaller the gaps, the closer the emitter is to the NPs.
At short distances, the evanescent field of the emitter can more easily
coupled to plasmonic nonradiative modes channeling more energy toward Ohmic
losses, and so reducing the quantum yield. Note that the nonmonotonic behavior
in $\eta_{ex}$ and $q$ is mainly due to the shift of the plasmon resonance as
the gap varies.

Having seen singularly the behavior of $\eta_{ex}$ and $q$, we can now analyze
the fluorescence enhancement shown in Fig. 4 for the LRA, TF-HT and QHT,
respectively. For gaps $g>1$ nm, our results are in line with previously
published results for the LRA [22] and TF-HT [41, 42, 43, 46]. Namely, the
maximum fluorescence enhancement is found in the correspondence of the dipolar
plasmon resonance (DPR). At smaller gaps ($g\simeq 0.5$ nm), high-order modes
become more and more relevant in terms of fluorescence enhancement. This
result however is not general and should be material-dependent. We would like
to mention that charge transfer modes can be observed in dimers, but they
appear at energies below 1 eV [53] and are not present in the spectral range
shown in Figs. 4 and 5. As the gap closes, in fact, the DPR is shifted to much
longer wavelengths for which Na ohmic losses increase, and both quantum yield
and field enhancement reduce (see Fig. S1 and Fig. S2 in the Supplemental
Material (SM). This is clearly visible for the LRA, although the TF-HT and QHT
show similar results but with more spread-out resonances due to the effect of
nonlocality. Note that for a fixed gap, increasing the NP radius has the same
effect as reducing the gap [41]. This is formally true for the LRA in the
quasi-static limit, where spectra are exactly the same except for an overall
factor (related to the size of the cavity) [54, 55].

At $g\lesssim 0.5$ nm and for $R=2.5$ nm, the maximum value of $\eta_{em}$ is
at the dipolar plasmon resonance (DPR) for the TF-HT and QHT, but for the LRA,
it is at the higher-order mode with a lower wavelength. If we remain in the
same interval $g\lesssim 0.5$ but increase the radius, the maximum value
gradually moves to the higher-order modes. These results mark the importance
of considering the contribution of higher-order modes for smaller gaps, as it
is also confirmed by TD-DFT calculations [53]. At $g\gtrsim 0.5$, the main
contribution to the fluorescence comes from the DPR for all cases. In
addition, the LRA predicts almost more than an order of magnitude higher
values as compared to the TF-HT and QHT. Also, the LRA results in more peaks
than the TF-HT and QHT for fixed radius. This could be attributed to the fact
that the $\eta_{ex}$ is higher for the LRA, and hence higher-order modes can
be distinguished from the DPR. As the radius increases, more higher-order
modes appear in all cases [56], and the DPR moves to the longer wavelengths.
On the other hand, by widening the gap, the spheres become more separated, and
the DPR moves to shorter wavelengths.

Let us now focus on the case with $R=10$ nm shown in Fig. 5. The key
difference of the QHT from the LRA and TF-HT is the quenching of fluorescence
at gap sizes below $\approx 0.5$ nm. This is somewhat an expected result from
a physics standpoint since fluorescence cannot be arbitrarily enhanced. The
QHT then predicts the maximum fluorescence enhancement to occur at $g\approx
0.5$ nm before the electron tunneling starts to have a detrimental impact.
Indeed, corrections introduced by the TF pressure are not enough to account
for the quantum nature of the electrons at the metal surface, as it is
confirmed by the monotonic increase of $\eta_{em}$ in Figs. 5(a). In other
words, just considering the nonlocality via electron pressure term still
results in an overestimated fluorescence enhancement. This has also been
confirmed by considering spill-out and nonlocality via Feibelman
$d$-parameters [28] and has resulted in the broadening of the Purcell
enhancement spectra near the plasmon frequency in emitter Na sphere and
emitter Na surface interactions. The LRA, on the other hand, is likely to
largely overestimate the attainable fluorescence enhancement at small gaps.
These result then indicates the necessity of considering both nonlocality and
electron tunneling effects to correctly predict the fluorescence enhancement
at sub-nanometer gaps.

Figure 6: The dependence of fluorescence enhancement on the wavelength for
different gap values as calculated from the QHT, PGSL and PGSLN. The radius of
spheres is $R=2.5$ nm.

Finally, we compare the fluorescence enhancement as calculated from the QHT,
QHT-PGSL, and QHT-PGSLN in Fig. 6 for a sphere dimer with $R=2.5$ nm for
different values of the gap. In our previous article [39], we showed that the
PGSL and PGSLN do not give the nonphysical oscillations in the absorption
spectrum of sphere dimers as compared to the QHT. Nevertheless, the
contribution of those oscillations to the fluorescence spectra is very small
compared to other modes. These oscillations should appear before $\lambda=350$
nm in Figs. 4 (g)-(i) and Fig. 5 (c) but clearly, we do not observe them. This
is the reason why we choose the QHT over Laplacian models to compare with the
LRA and TF-HT. It is interesting to observe that for all approaches the
fluorescence is quenched as opposed to the LRA and TF-HT results. This result
again emphasizes the importance of considering the spill-out of the electron
density at the metal surface and tunneling effects between closely spaced NPs,
as it is incorporated in all quantum approaches (QHT, QHT-PGSL and QHT-PGSLN).
As we can see, the biggest enhancement is achieved at $g=0.5$ nm for all
scenarios and the highest value is obtained for the QHT-PGSL in Fig. 6 (b).
The latter result is expected since the QHT-PGSL slightly overestimates the
field inside the gap [39]. On the other hand, the QHT-PGSLN does not suffer
from such overestimation and, although more complex to implement, it should be
considered as the most accurate. Interestingly, these results do not differ
much from the results obtained with the much simpler QHT model.

##  V Conclusions and Future Perspectives

To conclude, we demonstrated that the consideration of nonlocality and
tunneling in the framework of the QHT is crucial for the correct evaluation of
fluorescence enhancement in plasmonic sphere dimers. While the LRA predicts
increasing enhancement [22, 23] as the NPs distance is decreased, the use of
the QHT results in the quenching of the fluorescence of the emitter. In
addition, only accounting for the nonlocality in the framework of the TF-HT
still results in an increasing rate of fluorescence. Results similar to the
QHT are obtained if the second-order derivatives of the density are included
in the definition of the KE functional associated to the electron dynamic in
the metallic volumes in the form of PGSL and PGSLN functionals. The
fluorescence quenching happens at approximately $g=0.5$ nm for all scenarios.
These results should also be valid for similar systems, such as the NP on film
[22].

More generally, our method has allowed to introduce a fundamental mechanism
for fluorescence quenching, which marks the limits for the local and
hydrodynamic theories. Even though our work is purely theoretical,
experimental validation of such mechanism could be in principle achieved by
using analogous larger scale systems, such as a heavily-doped semiconductor
systems, or by using single atom emitters at cryogenic temperatures, by
exploiting the strong gradients of the electric field at the dimer gap in
order to create an optical trap. Another possibility, is that of exploiting
the weaker transverse dipole of traditional molecules, which would allow the
molecule to lay flat in the gap, thus reaching distances below $\sim 0.5$ nm.

Although we have considered only quantum effects in the dielectric response of
the NPs, while assuming a point-dipole approximation for the emitter, we
believe that equally important contribution could be introduced by taking into
account the spatial extension of the quantum emitter beyond the point-dipole
approximation [57] and the quantum nature of light [58].

## Acknowledgments

The authors thank Antonio I. Fernández-Domínguez for insightful discussions
and helpful comments.

##  Appendix A FEM implementation

We use FEM in COMSOL Multiphysics to solve the involved partial differential
equations. Following the usual FEM procedure (2a) and (2b) are multiplied by
the test functions $\textrm{{E}}_{\rm{I}}$ and P̃, and integration by parts is
performed. Also, since $\left(\frac{\delta G\left[n\right]}{\delta
n}\right)_{1}$ can have fourth-order of $n_{1}$ derivatives for PGSL and PGSLN
[see Expr. S5 in the SM], we introduce $\textrm{{F}}=\nabla n_{1}$ and
$\textrm{{O}}=\nabla\left(\nabla^{2}n_{1}\right)=\nabla^{2}\textrm{{F}}$
variables to lower the order to the first. In addition, since the geometry of
dimer is axisymmetric, the 2.5D technique is used [48, 38, 59]. Following it,
the vector fields are expressed in cylindrical harmonics
$v_{\rho,z,\phi}(\rho,z,\phi)=\sum_{m}v_{\rho,z,\phi}^{(m)}(\rho,z)e^{-im\phi}$,
where $\textrm{{v}}=\textrm{{E}};\textrm{{P}};\textrm{{F}};\textrm{{O}}$, and
$m\in\mathbb{Z}$. The test functions depend on $m$ via $e^{im\phi}$ form.
Thus, we will have $2m_{\textrm{max}}+1$ two-dimensional problems
($m_{\textrm{max}}$ is the maximum value for $m$), instead of the original
three-dimensional problem. The final system of equations for all scenarios is
the following:

$\displaystyle{\color[rgb]{0,0,0}\int\left\\{\left(\nabla\times\textrm{{E}}^{\left(m\right)}_{\rm{I}}\right)\cdot\left(\nabla\times\mathrm{\mathbf{\tilde{E}}}^{\left(m\right)}_{\rm{I}}\right)-\left(\frac{\omega^{2}}{c^{2}}\textrm{{E}}^{(m)}_{\rm{I}}+\mu_{0}\omega^{2}\textrm{{P}}^{(m)}\right)\cdot\mathrm{\mathbf{\tilde{E}}}^{\left(m\right)}_{\rm{I}}\right\\}\rho
d\rho dz=0,}$ (13a)
$\displaystyle{\color[rgb]{0,0,0}\int\left\\{-\frac{e}{m_{e}}\left(\frac{\delta
G\left[n\right]}{\delta
n}\right)_{1}^{\left(m\right)}\left(\nabla\cdot\mathrm{\mathbf{\tilde{P}}}^{\left(m\right)}\right)+\frac{1}{n_{0}}\left[\left(\omega^{2}+i\gamma\omega\right)\textrm{{P}}^{(m)}+\epsilon_{0}\omega_{p}^{2}\left(\textrm{{E}}^{(m)}_{\rm{I}}+\textrm{{E}}^{(m)}_{\textrm{0}}\right)\right]\cdot\mathrm{\mathbf{\tilde{P}}}^{\left(m\right)}\right\\}\rho
d\rho dz=0,}$ (13b)
$\displaystyle{\color[rgb]{0,0,0}\int\left\\{\left(\nabla\cdot\textrm{{P}}^{\left(m\right)}\right)\left(\nabla\cdot\mathrm{\mathbf{\tilde{F}}}^{\left(m\right)}\right)+e\textrm{{F}}^{\left(m\right)}\cdot\mathrm{\mathbf{\tilde{F}}}^{\left(m\right)}\right\\}\rho
d\rho dz=0,}$ (13c)
$\displaystyle{\color[rgb]{0,0,0}\int\left\\{\left(\nabla\cdot\textrm{{F}}^{\left(m\right)}\right)\left(\nabla\cdot\mathrm{\mathbf{\tilde{O}}}^{\left(m\right)}\right)+\textrm{{O}}^{\left(m\right)}\cdot\mathrm{\mathbf{\tilde{O}}}^{\left(m\right)}\right\\}\rho
d\rho dz=0.}$ (13d)

For the LRA, we solved only the wave equation:

${\color[rgb]{0,0,0}\int\left\\{\left(\nabla\times\textrm{{E}}^{\left(m\right)}_{\rm{I}}\right)\cdot\left(\nabla\times\mathrm{\mathbf{\tilde{E}}}^{\left(m\right)}_{\rm{I}}\right)-\frac{\omega^{2}}{c^{2}}\left(\textrm{{E}}^{(m)}_{\rm{I}}-\left(\varepsilon_{r}-1\right)\textrm{{E}}^{(m)}_{\textrm{0}}\right)\cdot\mathrm{\mathbf{\tilde{E}}}^{\left(m\right)}_{\rm{I}}\right\\}\rho
d\rho dz=0,}$ (14)

with Drude permittivity
$\varepsilon=1-\frac{\omega_{p}^{2}}{\omega^{2}+i\omega\gamma}$. In the case
of the TF-HT, Eqs. (13a) and (13b) need to be solved.

The similar procedure should be done for the evaluation of fields for decay
processes, with the exception that now $\mathbf{E}_{0}$ is the $z$-oriented
dipole field in free-space [43].

We used perfectly matched layers (PMLs) to emulate an infinite domain and
avoid unwanted reflections, and a zero flux boundary condition is imposed on
the electric field at the outer boundary of the PML. Also, the Dirichlet
boundary conditions $\textrm{{P}}=0,\textrm{{F}}=0$, and $\textrm{{O}}=0$ are
set on the polarization domain (PD). For the TF-HT, an additional hard-wall
boundary condition $\mathbf{P}\cdot\hat{\mathbf{n}}=0$ [38] is set on the NP
surface.

The meshing is done by using denser mesh on the surface of the dimer and on
the artificial region around the dipole (see Figs. 3b, and 3c in the SM). We
found that using a mapped mesh with rectangular elements at the location of
the dipole results in a better convergence for the QHT, PGSL and PGSLN (see
Fig. 3c in the SM). For other scenarios and other regions of the calculation
domain we used only triangular meshing. For the PML, the mapped mesh with
quadrilateral elements is used in all cases.

## References

  * Purcell  _et al._ [1946] E. M. Purcell, H. C. Torrey, and R. V. Pound, Resonance absorption by nuclear magnetic moments in a solid, Phys. Rev. 69, 37 (1946).
  * Goy  _et al._ [1983] P. Goy, J. M. Raimond, M. Gross, and S. Haroche, Observation of Cavity-Enhanced Single-Atom Spontaneous Emission, Phys. Rev. Lett. 50, 1903 (1983).
  * Tam  _et al._ [2007] F. Tam, G. P. Goodrich, B. R. Johnson, and N. J. Halas, Plasmonic Enhancement of Molecular Fluorescence, Nano Lett. 7, 496 (2007).
  * Novotny and Hecht [2012] L. Novotny and B. Hecht,  _Principles of Nano-Optics_, 2nd ed. (Cambridge University Press, 2012).
  * Lakowicz [2006a] J. R. Lakowicz,  _Principles of Fluorescence Spectroscopy_, 3rd ed. (Springer, New York, 2006).
  * Gaponenko [2015] S. Gaponenko, Fluorescence and Raman Scattering in Plasmonic Nano-structures: From Basic Science to Applications, in  _Nano-Structures for Optics and Photonics_, edited by B. Di Bartolo, J. Collins, and L. Silvestri (Springer Netherlands, Dordrecht, 2015) pp. 323–338.
  * Lakowicz [2006b] J. R. Lakowicz, Plasmonics in Biology and Plasmon-Controlled Fluorescence, Plasmonics 1, 5 (2006b).
  * Lakowicz  _et al._ [2006] J. R. Lakowicz, M. H. Chowdury, K. Ray, J. Zhang, Y. Fu, R. Badugu, C. R. Sabanayagam, K. Nowaczyk, H. Szmacinski, K. Aslan, and C. D. Geddes, Plasmon-controlled fluorescence: A new detection technology, in  _Biomedical Optics 2006_, edited by T. Vo-Dinh, J. R. Lakowicz, and Z. Gryczynski (San Jose, CA, 2006) p. 609909.
  * Puerto  _et al._ [2021] A. Puerto, C. López-Fernández, J. L. Bella, I. Elvira, G. Minguez-Vega, A. García-Cabañes, and M. Carrascosa, Fluorescence enhancement based on plasmonic nanoparticle structures on ferroelectric platforms for bioimaging applications, in  _Optical Methods for Inspection, Characterization, and Imaging of Biomaterials V_, Vol. 11786, edited by P. Ferraro, S. Grilli, M. Ritsch-Marte, and C. K. Hitzenberger, International Society for Optics and Photonics (SPIE, 2021) pp. 227 – 234\.
  * Li  _et al._ [2015] M. Li, S. K. Cushing, and N. Wu, Plasmon-enhanced optical sensors: A review, Analyst 140, 386 (2015).
  * Langer  _et al._ [2015] J. Langer, S. M. Novikov, and L. M. Liz-Marzán, Sensing using plasmonic nanostructures and nanoparticles, Nanotechnology 26, 322001 (2015).
  * Fernández-Domínguez  _et al._ [2018] A. I. Fernández-Domínguez, S. I. Bozhevolnyi, and N. A. Mortensen, Plasmon-Enhanced Generation of Nonclassical Light, ACS Photonics 5, 3447 (2018).
  * Darvill  _et al._ [2013] D. Darvill, A. Centeno, and F. Xie, Plasmonic fluorescence enhancement by metal nanostructures: Shaping the future of bionanotechnology, Phys. Chem. Chem. Phys. 15, 15709 (2013).
  * Yang  _et al._ [2020a] G. Yang, Q. Shen, Y. Niu, H. Wei, B. Bai, M. H. Mikkelsen, and H.-B. Sun, Unidirectional, Ultrafast, and Bright Spontaneous Emission Source Enabled By a Hybrid Plasmonic Nanoantenna, Laser & Photonics Reviews 14, 1900213 (2020a).
  * Cruz  _et al._ [2020] D. F. Cruz, C. M. Fontes, D. Semeniak, J. Huang, A. Hucknall, A. Chilkoti, and M. H. Mikkelsen, Ultrabright Fluorescence Readout of an Inkjet-Printed Immunoassay Using Plasmonic Nanogap Cavities, Nano Lett. 20, 4330 (2020).
  * Dulkeith  _et al._ [2002] E. Dulkeith, A. C. Morteani, T. Niedereichholz, T. A. Klar, J. Feldmann, S. A. Levi, F. C. J. M. van Veggel, D. N. Reinhoudt, M. Möller, and D. I. Gittins, Fluorescence Quenching of Dye Molecules near Gold Nanoparticles: Radiative and Nonradiative Effects, Phys. Rev. Lett. 89, 203002 (2002).
  * Kühn  _et al._ [2006] S. Kühn, U. Håkanson, L. Rogobete, and V. Sandoghdar, Enhancement of Single-Molecule Fluorescence Using a Gold Nanoparticle as an Optical Nanoantenna, Phys. Rev. Lett. 97, 017402 (2006).
  * Anger  _et al._ [2006] P. Anger, P. Bharadwaj, and L. Novotny, Enhancement and Quenching of Single-Molecule Fluorescence, Phys. Rev. Lett. 96, 113002 (2006).
  * Bharadwaj and Novotny [2007] P. Bharadwaj and L. Novotny, Spectral dependence of single molecule fluorescence enhancement, Opt. Express 15, 14266 (2007).
  * Faggiani  _et al._ [2015] R. Faggiani, J. Yang, and P. Lalanne, Quenching, Plasmonic, and Radiative Decays in Nanogap Emitting Devices, ACS Photonics 2, 1739 (2015).
  * Kongsuwan  _et al._ [2018a] N. Kongsuwan, A. Demetriadou, R. Chikkaraddy, J. J. Baumberg, and O. Hess, Fluorescence enhancement and strong-coupling in faceted plasmonic nanocavities, EPJ Appl. Metamat. 5, 6 (2018a).
  * Kongsuwan  _et al._ [2018b] N. Kongsuwan, A. Demetriadou, R. Chikkaraddy, F. Benz, V. A. Turek, U. F. Keyser, J. J. Baumberg, and O. Hess, Suppressed Quenching and Strong-Coupling of Purcell-Enhanced Single-Molecule Emission in Plasmonic Nanocavities, ACS Photonics 5, 186 (2018b).
  * Yang  _et al._ [2020b] B. Yang, G. Chen, A. Ghafoor, Y. Zhang, Y. Zhang, Y. Zhang, Y. Luo, J. Yang, V. Sandoghdar, J. Aizpurua, Z. Dong, and J. G. Hou, Sub-nanometre resolution in single-molecule photoluminescence imaging, Nat. Photonics 14, 693 (2020b).
  * Cuartero-González and Fernández-Domínguez [2020] A. Cuartero-González and A. I. Fernández-Domínguez, Dipolar and quadrupolar excitons coupled to a nanoparticle-on-mirror cavity, Physical Review B 101, 035403 (2020).
  * Cuartero-González and Fernández-Domínguez [2018] A. Cuartero-González and A. I. Fernández-Domínguez, Light-forbidden transitions in plasmon-emitter interactions beyond the weak coupling regime, ACS Photonics 5, 3415 (2018).
  * Ciracì  _et al._ [2012] C. Ciracì, R. T. Hill, J. J. Mock, Y. Urzhumov, A. I. Fernandez-Dominguez, S. A. Maier, J. B. Pendry, A. Chilkoti, and D. R. Smith, Probing the Ultimate Limits of Plasmonic Enhancement, Science 337, 1072 (2012).
  * Raza  _et al._ [2013] S. Raza, N. Stenger, S. Kadkhodazadeh, S. V. Fischer, N. Kostesha, A.-P. Jauho, A. Burrows, M. Wubs, and N. A. Mortensen, Blueshift of the Surface Plasmon Resonance in Silver Nanoparticles Studied with EELS, Nanophotonics 2, 131 (2013).
  * Gonçalves  _et al._ [2020] P. A. D. Gonçalves, T. Christensen, N. Rivera, A.-P. Jauho, N. A. Mortensen, and M. Soljačić, Plasmon–Emitter Interactions at the Nanoscale, Nat Commun 11, 366 (2020).
  * E. Fermi [1927] E. Fermi, Un Metodo Statistico per la Determinazione di alcune Prioprietà dell’Atomo, Rend. Accad. Naz. Lincei 6, 602 (1927).
  * Thomas [1927] L. H. Thomas, The Calculation of Atomic Fields, Math. Proc. Camb. Philos. Soc. 23, 542 (1927).
  * Raza  _et al._ [2011] S. Raza, G. Toscano, A.-P. Jauho, M. Wubs, and N. A. Mortensen, Unusual Resonances in Nanoplasmonic Structures due to Nonlocal Response, Phys. Rev. B 84, 121412(R) (2011).
  * Ciracì  _et al._ [2013] C. Ciracì, J. B. Pendry, and D. R. Smith, Hydrodynamic Model for Plasmonics: A Macroscopic Approach to a Microscopic Problem, ChemPhysChem 14, 1109 (2013).
  * Tserkezis  _et al._ [2017] C. Tserkezis, N. A. Mortensen, and M. Wubs, How nonlocal damping reduces plasmon-enhanced fluorescence in ultranarrow gaps, Physical Review B 96, 085413 (2017).
  * C. F. von Weizsäcker [1935] C. F. von Weizsäcker, Zur Theorie der Kernmassen, Z. Phys. 96, 431 (1935).
  * Bohm [1952a] D. Bohm, A Suggested Interpretation of the Quantum Theory in Terms of ”Hidden” Variables. II, Phys. Rev. 85, 180 (1952a).
  * Bohm [1952b] D. Bohm, A Suggested Interpretation of the Quantum Theory in Terms of ”Hidden” Variables. I, Phys. Rev. 85, 166 (1952b).
  * Toscano  _et al._ [2015] G. Toscano, J. Straubel, A. Kwiatkowski, C. Rockstuhl, F. Evers, H. Xu, N. Asger Mortensen, and M. Wubs, Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics, Nat Commun 6, 7132 (2015).
  * Ciracì and Della Sala [2016] C. Ciracì and F. Della Sala, Quantum Hydrodynamic Theory for Plasmonics: Impact of the Electron Density Tail, Phys. Rev. B 93, 205405 (2016).
  * Baghramyan  _et al._ [2021] H. M. Baghramyan, F. Della Sala, and C. Ciracì, Laplacian-Level Quantum Hydrodynamic Theory for Plasmonics, Phys. Rev. X 11, 011049 (2021).
  * Constantin  _et al._ [2018] L. A. Constantin, E. Fabiano, and F. Della Sala, Semilocal Pauli–Gaussian Kinetic Functionals for Orbital-Free Density Functional Theory Calculations of Solids, J. Phys. Chem. Lett. 9, 4385 (2018).
  * Christensen  _et al._ [2014] T. Christensen, W. Yan, S. Raza, A.-P. Jauho, N. A. Mortensen, and M. Wubs, Nonlocal Response of Metallic Nanospheres Probed by Light, Electrons, and Atoms, ACS Nano 8, 1745 (2014).
  * Tserkezis  _et al._ [2016] C. Tserkezis, N. Stefanou, M. Wubs, and N. A. Mortensen, Molecular fluorescence enhancement in plasmonic environments: exploring the role of nonlocal effects, Nanoscale 8, 17532 (2016).
  * Jurga  _et al._ [2017] R. Jurga, S. D’Agostino, F. Della Sala, and C. Ciracì, Plasmonic Nonlocal Response Effects on Dipole Decay Dynamics in the Weak- and Strong-Coupling Regimes, J. Phys. Chem. C 121, 22361 (2017).
  * Mortensen  _et al._ [2014] N. A. Mortensen, S. Raza, M. Wubs, T. Søndergaard, and S. I. Bozhevolnyi, A generalized non-local optical response theory for plasmonic nanostructures, Nat Commun 5, 3809 (2014).
  * Ciracì [2017] C. Ciracì, Current-Dependent Potential for Nonlocal Absorption in Quantum Hydrodynamic Theory, Phys. Rev. B 95, 245434 (2017).
  * Ciracì  _et al._ [2019] C. Ciracì, R. Jurga, M. Khalid, and F. Della Sala, Plasmonic Quantum Effects on Single-Emitter Strong Coupling, Nanophotonics 8, 1821 (2019).
  * Ciracì  _et al._ [2014] C. Ciracì, A. Rose, C. Argyropoulos, and D. R. Smith, Numerical studies of the modification of photodynamic processes by film-coupled plasmonic nanoparticles, J. Opt. Soc. Am. B 31, 2601 (2014).
  * Ciracì  _et al._ [2013] C. Ciracì, Y. A. Urzhumov, and D. R. Smith, Effects of Classical Nonlocality on the Optical Response of Three-Dimensional Plasmonic Nanodimers, J. Opt. Soc. Am. B 30, 2731 (2013).
  * Khalid and Ciracì [2019] M. Khalid and C. Ciracì, Numerical Analysis of Nonlocal Optical Response of Metallic Nanoshells, Photonics 6, 39 (2019).
  * Zhu  _et al._ [2016] W. Zhu, R. Esteban, A. G. Borisov, J. J. Baumberg, P. Nordlander, H. J. Lezec, J. Aizpurua, and K. B. Crozier, Quantum Mechanical Effects in Plasmonic Structures with Subnanometre Gaps, Nat Commun 7, 11495 (2016).
  * Esteban  _et al._ [2012] R. Esteban, A. G. Borisov, P. Nordlander, and J. Aizpurua, Bridging Quantum and Classical Plasmonics with a Quantum-Corrected Model, Nat Commun 3, 825 (2012).
  * Benz  _et al._ [2015] F. Benz, C. Tserkezis, L. O. Herrmann, B. de Nijs, A. Sanders, D. O. Sigle, L. Pukenas, S. D. Evans, J. Aizpurua, and J. J. Baumberg, Nanooptics of Molecular-Shunted Plasmonic Nanojunctions, Nano Lett. 15, 669 (2015).
  * Barbry  _et al._ [2015] M. Barbry, P. Koval, F. Marchesin, R. Esteban, A. G. Borisov, J. Aizpurua, and D. Sánchez-Portal, Atomistic Near-Field Nanoplasmonics: Reaching Atomic-Scale Resolution in Nanooptics, Nano Lett. 15, 3410 (2015).
  * Mayergoyz  _et al._ [2005] I. D. Mayergoyz, D. R. Fredkin, and Z. Zhang, Electrostatic (plasmon) resonances in nanoparticles, Phys. Rev. B 72, 155412 (2005).
  * Aubry  _et al._ [2010] A. Aubry, D. Y. Lei, A. I. Fernández-Domínguez, Y. Sonnefraud, S. A. Maier, and J. B. Pendry, Plasmonic Light-Harvesting Devices over the Whole Visible Spectrum, Nano Lett. 10, 2574 (2010).
  * García de Abajo [1999] F. J. García de Abajo, Multiple scattering of radiation in clusters of dielectrics, Phys. Rev. B 60, 6086 (1999).
  * Neuman  _et al._ [2018] T. Neuman, R. Esteban, D. Casanova, F. J. García-Vidal, and J. Aizpurua, Coupling of Molecular Emitters and Plasmonic Cavities beyond the Point-Dipole Approximation, Nano Letters 18, 2358 (2018).
  * Feist  _et al._ [2020] J. Feist, A. I. Fernández-Domínguez, and F. J. García-Vidal, Macroscopic QED for quantum nanophotonics: emitter-centered modes as a minimal basis for multiemitter problems, Nanophotonics 10, 477 (2020).
  * Ciracì  _et al._ [2013] C. Ciracì, Y. Urzhumov, and D. R. Smith, Far-Field Analysis of Axially Symmetric Three-Dimensional Directional Cloaks, Opt. Express 21, 9397 (2013).


#  The Impact of Vaccine Hesitancy on Epidemic Spreading

C. H. Leung, María E. Gibbs, and Philip E. Paré *C. H. Leung and P. E. Paré
are with the School of Electrical and Computer Engineering, Purdue University.
M. E. Gibbs is with the Department of Aerospace Engineering and Engineering
Mechanics, The University of Texas at Austin. This work was supported in part
by the National Science Foundation, grants NSF-CNS #2028738 and NSF-ECCS
#2032258.

###### Abstract

The COVID-19 pandemic has devastated the world in an unprecedented way,
causing enormous loss of life. Time and again, public health authorities have
urged people to become vaccinated to protect themselves and mitigate the
spread of the disease. However, vaccine hesitancy has stalled vaccination
levels in the United States. This study explores the effect of vaccine
hesitancy on the spread of disease by introducing an SIRS-Vκ model, with
compartments of susceptible (S), infected (I), recovered (R), and vaccinated
(V). We leverage the concept of carrying capacity to account for vaccine
hesitancy by defining a vaccine confidence level $\kappa$, which is the
maximum number of people that will become vaccinated during the course of a
disease. The inverse of vaccine confidence is vaccine hesitance,
$(\frac{1}{\kappa})$. We explore the equilibria of the SIRS-Vκ model and their
stability, and illustrate the impact of vaccine hesitance on epidemic spread
analytically and via simulations.

##  I Introduction

The origins of epidemiological modeling extend as far back as the 18th century
[1], with Bernoulli’s 1760 treatise on smallpox considered one of its
foundational documents [2]. In the 20th century, the field experienced
significant growth, including the introduction of approaches that employed
dynamic systems theory [1]. In particular, Kermack and McKendrick’s 1932
Contributions to the Mathematical Theory of Epidemics introduced compartmental
models, which are heavily used to this day [3, 4]. With this method, members
of the population are compartmentalized based on their current state of
health. The most common models are the SIS, SIR, and SEIR models, in which
people are grouped as: susceptible (S), exposed (E), infected (I), or
recovered (R) [4]. Expansions of the aforementioned models have also been
developed. For instance, in light of the COVID-19 pandemic, Giordano et al.
constructed a SIDARTHE-V model that subdivides the infected population
according to the detection and severity of their symptoms and includes
vaccination levels [5].

We include a vaccinated compartment in our model because vaccines are a
powerful tool in mitigating epidemic spread and saving lives [1]. The modern
practice of vaccination is rooted in Edward Jenner’s 1796 experiments on
smallpox [6, 7]. Since then, incredible breakthroughs have been made in
vaccinations, including the use of vaccines to eradicate smallpox [8]. Today,
vaccines are widely used around the world as early as infancy to prevent
infections [9]. Due to the power of vaccines, the drastic death rates caused
by the spread of the novel SARS-CoV-2 virus have spurred worldwide vaccination
efforts to combat the virus [10, 11]. In the U.S., three vaccines are
currently available: Pfizer-BioNTech, Moderna, and Johnson & Johnson’s Janssen
[12].

Nevertheless, vaccine hesitancy remains a real, growing concern for public
health [13]. An individual’s vaccine hesitancy is affected by context,
interpersonal experiences, and disease-specific metrics [14]. In the United
States, vaccine hesitancy is especially apparent regarding the reception of
the COVID-19 vaccine. As of September 11, 2021, 53% of Americans have been
fully vaccinated against COVID-19, with an additional 9.2% partially
vaccinated [11]. However, these levels are lower than those predicted by a
February 2021 survey by Pew Research Center. According to the survey, 69% of
Americans reported they would probably, definitely, or had already been at
least partially vaccinated [15]. In fact, the number of vaccines administered
daily in the U.S. in September 2021 is less than half of those administered
daily in mid-April 2021 [11].

There have been dozens of studies on COVID-19 vaccine hesitance around the
world [16, 17, 18]. There are also studies on how vaccination levels affect
the dynamics of disease spread. Pires and Crokidakis explore the effect of
vaccine opinions in the spread of a disease [19]. Similarly, [20] considers
the dynamics when opinions are continuous. However, to the best of the
authors’ knowledge, the correlation between a maximum vaccination capacity
$\kappa$ and the equilibria of system dynamics [21] has not yet been
established. This work analytically illustrates the effect of $\kappa$ on the
dynamics of COVID-19 spread. Specifically, $\kappa$ is shown to determine the
endemic and disease-free equilibria of the SIRS-Vκ system.

##  II SIRS-Vκ Model

In this paper, we take an SIRS model and add a vaccinated compartment with a
vaccine confidence $\kappa$ as the upper-bound of vaccination level. The
result is an SIRS-Vκ model: susceptible (S), infected (I), recovered (R),
vaccinated (V). Figure 1 gives a graphical representation of the compartmental
flow.

######  Definition II.1 (Vaccine Hesitance).

We define the vaccine hesitance as $\kappa^{-1}$, where $\kappa$ is the
vaccine confidence.

Each compartment in the grouped SIRS-Vκ model is a scalar between $0$ and $1$
modeling the fraction of the population in each epidemic compartment at time
$t$. The following are the equations describing the flow from one compartment
to another in the group SIRS-Vκ model:

$\displaystyle\dot{S}$ $\displaystyle=-\beta
SI-\rho\left(1-\frac{V}{\kappa}\right)S+\omega R$ (1a) $\displaystyle\dot{I}$
$\displaystyle=\beta SI-\gamma I$ (1b) $\displaystyle\dot{R}$
$\displaystyle=\gamma I-\omega R-\rho\left(1-\frac{V}{\kappa}\right)R$ (1c)
$\displaystyle\dot{V}$
$\displaystyle=\rho\left(1-\frac{V}{\kappa}\right)(S+R).$ (1d)

The parameters are defined as follows:

  * •

$\beta$ is the frequency-dependent transmission rate

  * •

$\gamma$ is the recovery rate

  * •

$\rho$ is the rate of vaccination roll-out

  * •

$\kappa$ is the vaccine confidence of the population

  * •

$\omega$ is the rate of waning natural immunity.

All of the above parameters are defined in the range of $(0,\infty)$.

Figure 1: SIRS-Vκ Compartments.

######  Assumption 1 (Frequency dependent and population preserving system).

We assume $S(t_{0}),I(t_{0}),R(t_{0})\in[0,1]$,
$V(t_{0})\in[0,\underline{\kappa}]$, where
$\underline{\kappa}:=\min\\{1,\kappa\\}$, and
$S(t_{0})+I(t_{0})+R(t_{0})+V(t_{0})=1$ $\forall i\in\\{1,...,n\\}$.
Furthermore, all parameters of the system are strictly positive.

######  Lemma II.1 (Permissible range of $V$).

Let Assumption 1 hold, then $0\leq V(t)\leq\underline{\kappa}\ \ \forall t\geq
t_{0},\forall i\in\\{1,...,n\\}$.

###### Proof.

Suppose Assumption 1 holds, we first note that
$\dot{V}(t_{0})=\rho(S(t_{0})+R(t_{0}))\geq 0$ when $V(t_{0})=0$. Furthermore,
$\lim_{V(t_{0})\downarrow 0}\dot{V}(t_{0})=\rho(S(t_{0})+R(t_{0}))\geq 0$.
Then, consider the upper bound of $V(t)$, $\dot{V}(t_{0})=0$ when
$V(t_{0})=\underline{\kappa}$. Moreover, for any
$V(t_{0})\leq\underline{\kappa}$,
$\lim_{V(t_{0})\uparrow\kappa}\dot{V}(t_{0})\geq 0$. Therefore,
$v(t)\in[0,\underline{\kappa}]$ $\forall t\geq t_{0}$ for all permissible
initial values in Assumption 1.

The following lemma shows that the bounds on the initial state apply to all
states for all $t\geq t_{0}$.

######  Lemma II.2 (Admissible range of states values).

Let Assumption 1 hold. Then $S(t),I(t),R(t),V(t)\in[0,1]$ and
$S(t)+I(t)+R(t)+V(t)=1$ $\forall t\geq 0,\forall i\in\\{1,...,n\\}$.

###### Proof.

We first observe that $\dot{S}(t)+\dot{I}(t)+\dot{R}(t)+\dot{V}(t)=0$ $\forall
t\geq 0$ by summing over (1a)-(1d), then by taking the integration:

$\displaystyle\int_{t_{0}}^{t}\dot{S}(\tau)+\dot{I}(\tau)+\dot{R}(\tau)+\dot{V}(\tau)d\tau$
$\displaystyle=0$ $\displaystyle S(t)+I(t)+R(t)+V(t)$ $\displaystyle=1\ \
\forall t\geq 0$

since the initial condition is given as
$S(t_{0})+I(t_{0})+R(t_{0})+V(t_{0})=1$.

To show that $S(t),I(t),R(t),V(t)\in[0,1]$ $\forall t\geq 0$, we rewrite the
state vector as $z(t)\in\mathbb{R}^{4}$, then we note that for any
$j\in\\{1,\ldots,4\\}$, $\dot{z}_{j}=-z_{j}f_{j}(z_{k\neq j})+g_{j}(z)$ where
$f_{j}(z_{k\neq j}),g_{j}(z)\geq 0$ $\forall z\geq 0$, also $z_{k\neq j}$ is
the state vector $z$ without its $j^{th}$ entry. This applies to
$\dot{z}_{j}=\dot{V}$ also, as we have demonstrated in Lemma II.1. Then, we
observe that $\lim_{z_{j}\downarrow 0}\dot{z}_{j}(t)=g_{j}(z)\geq 0\ \ \forall
z\geq 0$. Therefore, $z(t)\geq 0$ when $z(t_{0})\geq 0$. Then $z_{j}(t)\leq 1$
follows directly from $z_{j}(t)\geq 0$ and $S(t)+I(t)+R(t)+V(t)=1$.

We denote the set of admissible states of (1) by
$\mathcal{S}=\\{(S(t),I(t),R(t),V(t)):\ S(t)+I(t)+R(t)+V(t)=1\ \wedge\
S(t),I(t),R(t)\in[0,1]\ \wedge\ V(t)\in[0,\underline{\kappa}]\ \forall t\geq
t_{0}\\}$.

##  III Stability Analysis

This section will lay down the existence and uniqueness of the two equilibria
of the SIRS-Vκ model and their stability conditions. We denote an equilibrium
of the SIRS-Vκ model as $(S^{*},I^{*},R^{*},V^{*})$ where $S^{*}$, $I^{*}$,
$R^{*}$, and $V^{*}$ are the steady-state values of (1) as $t\to\infty$.

######  Definition III.1 (Disease-Free Equilibrium).

A disease-free equilibrium (DFE) is an equilibrium with the steady-state
infected level $I^{(d)}=0$ $\forall i\in\\{1,...,n\\}$, where
$\dot{S}^{(d)}=\dot{I}^{(d)}=\dot{R}^{(d)}=\dot{V}^{(d)}=0$.

Similarly, we define the endemic equilibrium as follows:

######  Definition III.2 (Endemic Equilibrium Point).

An endemic equilibrium point (EEP) is an equilibrium with the steady-state
infected level $I^{(e)}>0$, where
$\dot{S}^{(e)}=\dot{I}^{(e)}=\dot{R}^{(e)}=\dot{V}^{(e)}=0$.

######  Lemma III.1 (Strictly Positive EEP).

If there exists an endemic equilibrium $(S^{(e)},I^{(e)},R^{(e)},V^{(e)})$,
then it is strictly greater than $0$, i.e. $z^{(e)}>0$ $\forall
z\in\\{S,I,R,V\\}$. Furthermore, $V^{(e)}=\underline{\kappa}$.

###### Proof.

Assume, by way of contradiction, that $S^{(e)}=0$, then $\dot{S}^{(e)}=-\gamma
I\neq 0$, which contradicts the definition of EEP. Therefore, $S^{(e)}>0$. A
similar argument can be made by substituting $R^{(e)}=0$ into (1c) to get
$R^{(e)}>0$. Therefore, from (1d), we have

$0=\rho\left(1-\frac{V^{(e)}}{\kappa}\right)(S^{(e)}+R^{(e)})$ (2)

and the only way for the R.H.S. to equal zero is if $V^{(e)}=\kappa$. Notice
that (2) cannot be satisfied when $\kappa\geq 1$, and therefore, EEP only
exists when $\kappa<1$, concluding the proof.

Lastly, we reiterate the definition of the basic reproduction number through
the next generation matrix method [22, 23]:

######  Definition III.3 (Basic and Effective Reproduction Number).

Let $\dot{z_{i}}=f_{i}(z)-g_{i}(z)$ be the differential equation of the
$i^{th}$ infected compartment, where $f_{i}(z)$ is the function governing the
rate of appearance, $g_{i}(z)$ is the function governing the rate of
transferring into other compartments, and $i\in\\{1,\ldots,m\\}$ for $m$
infected compartments out of $n$ total compartments. Let
$F=\left[\frac{\partial f_{i}(z_{0})}{\partial z_{j}}\right]$ and
$G=\left[\frac{\partial g_{i}(z_{0})}{\partial z_{j}}\right]$ be the Jacobians
of $[f_{i}]$ and $[g_{i}]$, $\forall j\in\\{1,\ldots,m\\}$, then $FV^{-1}$ is
the next generation matrix, and its spectral radius $\sigma(FV^{-1})$ is the
basic reproduction number $\mathcal{R}_{0}$ of the system. $z_{0}$ is the
initial condition where full susceptibility of population is assumed
($S(t_{0})\approx 1$). The effective reproduction number is defined as
$\mathcal{R}_{t}:=S(t)\mathcal{R}_{0}$.

The $ij^{th}$ entry of $FV^{-1}$ is the expected number of secondary cases in
the $i^{th}$ compartment produced by an infected individual in the $j^{th}$
compartment [24]. In the case of a single infected compartment, it is simply
the average number of secondary infected cases introduced by one infected
individual.

######  Proposition 1.

The basic reproduction number $\mathcal{R}_{0}$ of (1) is
$\frac{\beta}{\gamma}$, and the effective reproduction number
$\mathcal{R}_{t}$ is upper bounded by $(1-V(t))\frac{\beta}{\gamma}$.

###### Proof.

From (1b), $f(I)=\beta SI$, $g(I)=\gamma I$. By computing the Jacobians, we
get $F=\frac{\partial\beta SI(t_{0})}{\partial I}=\beta S(t_{0})$ and
$G=\frac{\partial\gamma I}{\partial I}=\gamma$. Therefore,
$\mathcal{R}_{0}=\sigma(FG^{-1})=\frac{\beta}{\gamma}$ since $\mathcal{R}_{0}$
is defined for $S(t_{0})\approx 1$. Further,
$\mathcal{R}_{t}=S(t)\frac{\beta}{\gamma}=(1-V(t)-R(t))\frac{\beta}{\gamma}\leq(1-V(t))\frac{\beta}{\gamma}$.

Notice that we use the word uniqueness in the following text not to refer to
the uniqueness of the set of all equilibria of (1) but of a particular class
of equilibria, i.e. the DSF and EEP are sets of cardinality at most one. As we
will see in the following, the EEP could be an empty set under certain
conditions.

######  Proposition 2 (Existence and uniqueness of DFE).

There always exists a unique disease-free equilibrium
$(S^{(d)},I^{(d)},R^{(d)},V^{(d)})$ =
$(1-\underline{\kappa},0,0,\underline{\kappa})$.

###### Proof.

By the definition of a disease-free equilibrium and from (1b), we know the
$j^{th}$ DFE is characterized by $I^{(d_{j})}=0$, and $\dot{S}^{(d_{j})}$ $=$
$\dot{I}^{(d_{j})}$ $=$ $\dot{R}^{(d_{j})}$ $=$ $\dot{V}^{(d_{j})}$ $=$ $0$
$\forall j\in\mathbb{R}$. Equation (1c) can be evaluated as
$0=\left(\omega+\rho\left(1-\frac{V^{(d_{j})}}{\kappa}\right)\right)R^{(d_{j})}$
by substitution, which implies either:

$\displaystyle V^{(d_{1})}=\kappa\left(1+\frac{\omega}{\rho}\right)$ or (3)
$\displaystyle R^{(d_{2})}=0.$ (4)

In the first case (3), we can substitute $V^{(d_{1})}$ into (1d) and solve for
$S^{(d_{1})}+R^{(d_{1})}=0$. By Assumption 1 and Lemma II.2,
$S^{(d_{1})}=R^{(d_{1})}=0$ is the only solution to the equation. Together
with the presumptions $V^{(d_{1})}=\kappa\left(1+\frac{\omega}{\rho}\right)$
and $I^{(d_{1})}=0$, we realize that
$S^{(d_{1})}+I^{(d_{1})}+R^{(d_{1})}+V^{(d_{1})}\neq 1$, which is not in the
permissible domain of states by Assumption 1 and Lemma II.2. Therefore, DFE
$(d_{1})$ is not a feasible equilibrium.

Now, by (4), we can substitute $R^{(d_{2})}=0$ into (1a), then we have
$0=\rho\left(1-\frac{V^{(d_{2})}}{\kappa}\right)S^{(d_{2})}$, which leads to
two subcases:

$\displaystyle S^{(d_{21})}=0$ or (4a) $\displaystyle V^{(d_{22})}=\kappa.$
(4b)

In the first subcase (4a), since $I^{(d_{21})}=0$ by the definition of the DFE
and $R^{(d_{21})}=0$ by (4), we have
$(S^{(d_{21})},I^{(d_{21})},R^{(d_{21})},V^{(d_{21})})$ $=$ $(0,0,0,1)$, which
is true when $\kappa\geq 1$ by recalling Lemma II.1 and that
$\underline{\kappa}=\min\\{1,\kappa\\}$. The second subcase (4b), combined
with $I^{(d_{22})}=0$, (4), and Lemma II.2, gives
$(S^{(d_{22})},I^{(d_{22})},R^{(d_{22})},V^{(d_{22})})=(1-\kappa,0,0,\kappa)$
which is true when $\kappa\in(0,1]$. DFE $(d_{21})$ and $(d_{22})$ can be
combined as $(1-\underline{\kappa},0,0,\underline{\kappa})$ to be applicable
to the whole range of $\kappa\in(0,\infty)$. Since
$(1-\underline{\kappa},0,0,\underline{\kappa})$ is the only DFE satisfying
Definition III.1 and (1), it is unique under Assumption 1. Since
$(1-\underline{\kappa},0,0,\underline{\kappa})$ is a valid state under
Assumption 1 for all admissible values of $\kappa$, the DFE always exists.

######  Remark.

Note that $\kappa$ ceases to act as an upper bound on the vaccination state
while still affecting the effective rate of vaccination
$\rho\left(1-\frac{V}{\kappa}\right)$ when $\kappa$ is greater than $1$.
Moreover, $\rho\left(1-\frac{V}{\kappa}\right)\to\rho$ as $\kappa\to\infty$.

######  Proposition 3 (Uniqueness of EEP).

The endemic equilibrium point of (1):

$\left(S^{(e)},I^{(e)},R^{(e)},V^{(e)}\right)=\left(\frac{\gamma}{\beta},\frac{1-\kappa-\gamma/\beta}{1+\gamma/\omega},\frac{1-\kappa-\gamma/\beta}{1+\omega/\gamma},{\kappa}\right)$

is unique. Moreover,

$\frac{I^{(e)}}{R^{(e)}}=\frac{\omega}{\gamma}.$

###### Proof.

If an endemic equilibrium exists, by Definition III.2, $I^{(e)}>0$. Applying
this fact and setting (1b) to zero gives $S^{(e)}=\frac{\gamma}{\beta}$.
Therefore, the result from Lemma III.1, with $n=1$, implies $V^{(e)}=\kappa$.
By substituting $V^{(e)}=\kappa$ into (1c), we will have
$\frac{I^{(e)}}{R^{(e)}}=\frac{\omega}{\gamma}$. To solve for $I^{(e)}$ and
$R^{(e)}$, recall the global preservation of population result in Lemma II.2,
we have $I^{(e)}+R^{(e)}=1-\kappa-\frac{\gamma}{\beta}$. We can combine this
equation with $\frac{I^{(e)}}{R^{(e)}}=\frac{\omega}{\gamma}$, then solving
for the system of two equations with two unknowns gives us the desired result.
Therefore, we have shown the uniqueness of permissible EEP under Assumption 1.

Notice that since $I^{*}$ can either be $0$ or $(0,1]$, the DFE and the EEP
are the only two possible admissible equilibria of (1) under Assumption 1.

######  Lemma III.2 (Necessary and sufficient condition for existence of EEP).

The endemic equilibrium exists in (1) if and only if
$\kappa<1-\frac{\gamma}{\beta}$, which is equivalent to
$(1-\kappa)\frac{\beta}{\gamma}>1$.

###### Proof.

From Assumption 1 and Lemma II.2, the EEP in Proposition 3 exists if and only
if the range of parameters satisfy the four inequality:

$\begin{cases}0<\frac{\gamma}{\beta}<1\\\
0<\frac{1-\kappa-\gamma/\beta}{1+\gamma/\omega}<1\\\
0<\frac{1-\kappa-\gamma/\beta}{1+\omega/\gamma}<1\\\ 0<\kappa<1.\end{cases}$
(5)

The first and forth constraints are absorbed by Assumption 1 and the second
constraint. The second and third constraints can be reduced to
$0<1-\kappa-\gamma/\beta$, which is equivalent to $\kappa<1-\gamma/\beta$ and
$(1-\kappa)\frac{\beta}{\gamma}>1$.

Fig. 2 illustrates the set of states $\mathcal{S}$, which is a 3-dimensional
simplex with six vertices $(S,I,R,D,F,G)$ with coordinates $(1,0,0,0)$,
$(0,1,0,0)$, $(0,0,1,0)$, $(1-\kappa,1,0,\kappa)$, $(0,1-\kappa,0,\kappa)$,
and $(0,0,1-\kappa,\kappa)$, respectively. The larger simplex is a tetrahedron
because of Lemma II.2, which states that $S(t)+I(t)+R(t)+V(t)=1$ and
$(S(t),I(t),R(t),V(t))\in[0,1]^{4}$. To see this clearly, consider an SIR
compartmental model with three states, $(S,I,R)$. If $S+I+R=1$, then $(S,I,R)$
lie only on the surface of the unit circle defined by the 1-norm. Furthermore,
if $(S,I,R)\in[0,1]^{3}$, then $(S,I,R)$ lies only on the surface belonging to
the first quadrant of the 3-dimensional space, where $S$, $I$, $R$ are non-
negative. Notice that the surface is an equilateral triangle with its three
vertices $S=(1,0,0)$, $I=(0,1,0)$, and $R=(0,0,1)$. Generalizing this notion
to the 4-dimensional case gives us Fig. 2.

######  Theorem 1 (DFE GAS Conditions).

The DFE of (1) is globally asymptotically stable (GAS) if and only if
$(1-\kappa)\frac{\beta}{\gamma}\leq 1$.

###### Proof.

Our goal is to show by exhaustion of cases that all the states in
$\mathcal{S}$ are either the DFE or will approach the DFE as $t\to\infty$
under the condition $(1-\kappa)\frac{\beta}{\gamma}\leq 1$. Fig. 2 visualizes
the simplex of states, i.e. the range of admissible states $\mathcal{S}$, of
which we want to show GAS around the DFE.

Suppose $(1-\kappa)\frac{\beta}{\gamma}\leq 1$. If $S(t)+R(t)=0$, which is
equivalent to $S(t)=R(t)=0$, then either $(S(t)=R(t)=I(t)=0\wedge V(t)=1)$,
which is a special case of DFE when $\kappa\geq 1$, or $(S(t)=R(t)=0\wedge
I(t)>0)$. If $(S(t)=R(t)=0\wedge I(t)>0)$, then by (1c):

$\displaystyle\dot{R}(t)$ $\displaystyle=\gamma I(t)>0.$

Hence, every admissible state in (1) which satisfies $S(t)+R(t)=0$ will have
$\dot{R}(t)>0$, hence evolving to another state such that $S(t)+R(t)\neq 0$.
In other words, states on the $IV$ edge in Fig. 2 are not stable and will move
toward vertex $R$ as time unfolds. The only state on $IV$ which violates this
tendency of evolving towards $R$ is the vertex $V$, which is already the DFE
if it lies within $\mathcal{S}$. For every state that satisfies
$(S(t)+R(t)>0\wedge V(t)<\kappa)$,

$\displaystyle\dot{V}(t)$
$\displaystyle=\rho\left(1-\frac{V(t)}{\kappa}\right)(S(t)+R(t))$
$\displaystyle>0$

by substituting $S(t)+R(t)>0$ into (1d). Therefore, every state with
$V(t)<\kappa$ will satisfy $V(t)=\kappa$ as $t\to\infty$ or it is already the
DFE.

If $(I(t)=R(t)=0\wedge V(t)=\kappa)$, then $S(t)=1-\kappa$ by Assumption 1 and
Lemma II.2, which state that $S(t)+I(t)+R(t)+V(t)=1$ $\forall t\geq t_{0}$,
which is the DFE. Otherwise, $(V(t)=\kappa\wedge I(t)+R(t)>0)$. Thus, by
substituting $V(t)=\kappa$ into $S(t)+I(t)+R(t)+V(t)=1$, we have

$\displaystyle S(t)$ $\displaystyle=1-\kappa-I(t)-R(t)$
$\displaystyle\leq\frac{\gamma}{\beta}-I(t)-R(t)$ (6)
$\displaystyle<\frac{\gamma}{\beta},$ (7)

where (6) holds because the condition $(1-\kappa)\frac{\beta}{\gamma}\leq 1$
is equivalent to $\frac{\gamma}{\beta}\geq 1-\kappa$ and (7) holds because
$I(t)+R(t)>0$. Therefore, from (1b) and (7), we have

$\displaystyle\ \ \ \dot{I}(t)$
$\displaystyle=\beta\frac{\gamma}{\beta}I(t)-\gamma I(t)$
$\displaystyle=-2\gamma I(t)$ $\displaystyle\leq 0,$

where equality holds only if $I(t)>0$. Therefore, every state which satisfies
$V(t)=\kappa$ will satisfy $(V(t)=\kappa\wedge I(t)=0)$ as $t\to\infty$ or it
is already the DFE.

Lastly, for any state that satisfies $(V(t)=\kappa\wedge I(t)=0\wedge
R(t)>0)$, from (1c), we have

$\displaystyle\dot{R}(t)=-\omega R(t)<0.$

Thus, every state that satisfies $(V(t)=\kappa\wedge I(t)=0\wedge R(t)>0)$
will satisfy $(V(t)=\kappa\wedge I(t)=0\wedge R(t)=0)$ as $t\to\infty$ or it
is already the DFE. Therefore, the DFE is GAS if
$(1-\kappa)\frac{\beta}{\gamma}\leq 1$.

If $(1-\kappa)\frac{\beta}{\gamma}>1$, then there exists an EEP in
$\mathcal{S}$ by Lemma III.2, which concludes the other direction of the
proof.

Figure 2: Admissible range of state values of (1) when EEP exists. $D$ is the
DFE, and $E$ is the EEP. The plane $ABC$ represents all states with
$S(t)=\frac{\gamma}{\beta}$, and the plane $DFG$ represents all states with
$V(t)=\kappa$. The ratio between the segments
$\frac{\|EQ\|}{\|PE\|}=\frac{\omega}{\gamma}$.

Referring to the points in Fig. 2, $(1-\kappa)\frac{\beta}{\gamma}\leq 1$ is
achieved when the plane $ABD$ slides up or $DFG$ slides back, such that the
segment $PQ$ vanishes. One can show that all states above the $ABC$ plane have
$\dot{I}>1$ and below $\dot{I}<1$. An intuitive way to interpret Theorem 1 is
that when $ABC$ is strictly above $DFG$, all states are pushed away from $I$
and pulled toward $V$, until they reach the DFE.

######  Remark.

Notice that the necessary and sufficient GAS condition of the disease-free
equilibrium is not $\mathcal{R}_{0}<1$ for the SIRS-Vκ model, because the
basic reproduction number only guarantees sufficient local asymptotic
stability when it is less than one. Therefore, it might at times over/under-
estimate the threshold condition of the spread process.

##  IV Simulations

In this section, we investigate the impact of the vaccine confidence $\kappa$
on the behavior of the SIRS-Vκ model. Fig. 3 compares the trajectories of the
vaccination level $V(t)$ with different values of $\kappa$, ranging from $0.1$
to $\infty$. Note that, consistent with the analysis results, $\kappa$ slows
the rate of convergence and, if $\kappa<1$, can lower the limit of $V(t)$.
Fig. 3 and Fig. 4 were plotted with initial condition
$(S(t_{0}),I(t_{0}),R(t_{0}),V(t_{0}))=(0.54,0.41,0.05,0)$. Fig. 5 and Fig. 6
were plotted with initial condition
$(S(t_{0}),I(t_{0}),R(t_{0}),V(t_{0}))=(0.99,0.01,0,0)$, and the initial
condition of Fig. 7 is $(S(t_{0}),I(t_{0}),R(t_{0}),V(t_{0}))=(0.7,0.3,0,0)$.
All plots use the parameter set:
$(\beta=1.6,\gamma=.8,\rho=.12,\kappa=.8,\omega=.2)$ except for Fig. 4 and
Fig. 3 we have used $\omega=3$.

Figure 3: The impact of varying $\kappa$ on $V(t)$

Fig. 4 illustrates how varying $\kappa$ affects the infection $I(t)$. Since
$\frac{\beta}{\gamma}=2$ in this case, $\kappa=0.5$ barely satisfies the
condition for GAS in Theorem 1 to reach the DFE as time goes to infinity. We
also notice that $\kappa$ slows the rate of convergence of $I(t)$ when
comparing the cases $\kappa=(0.5,0.6,1.2,\infty)$.

Figure 4: The impact of varying $\kappa$ on $I(t)$

Fig. 5 compares $I(t)$ and $R(t)$ between the SIRS, SIRSV, and SIRS-Vκ models.
Notice the key difference between the three models is that SIRSV and SIRS-Vκ
has $\rho>0$ while $SIRS$ has $\rho=0$, and SIRSVκ has $\kappa=0.8<\infty$
while SIRSV assumes $\kappa\to\infty$. In this particular setup, both the
SIRSV and SIRSVκ models reach the DFE before $t=40$, while the SIRS model
settles at the EEP. Furthermore, we observe a higher infection peak with
slower peak time in the SIRS-V$\kappa$ model.

Figure 5: Comparation between SIRS, SIRSV, and SIRS-Vκ

models.

To investigate how the vaccine confidence affects the maximum peak infection
value and time, we plot the maximum peak infection value and time versus
$\kappa$ in Fig. 6. Note that the maximum peak infection value decreases
monotonically with $\kappa$. On the other hand, the peak infection time
reaches its maximum when $\kappa=0.2$ before decaying, which is not expected
and will need further investigation.

Figure 6: Maximum Infection value and peak time versus $\kappa$.

Fig. 7 plots the earliest time when $I(t)<0.001$ with respect to different
$(1-\kappa)\frac{\beta}{\gamma}$ through varying $\kappa$. We use the
condition $I(t)<0.001$ as a way to investigate the asymptotic behavior of the
system under different parameter settings. The program returns $-10$ if $I(t)$
fails to converge close enough to $0$ before $t=100$. The sharp drop at
$(1-\kappa)\frac{\beta}{\gamma}\approx 1$ aligns with our finding in Lemma
III.2 and Theorem 1 that the DFE ceases to be GAS when
$(1-\kappa)\frac{\beta}{\gamma}$ is larger than $1$. Note that the drop occurs
slightly before $1$ because we do not allow $t>100$, which is much smaller
than asymptotic time. Finally, all the simulations indicate that the EEP has a
large region of convergence as long as it exists in $\mathcal{S}$. More
investigation is required to firmly conclude these findings analytically.

Figure 7: Threshold behavior with respect to $(1-\kappa)\frac{\beta}{\gamma}$.

##  V Conclusion and Future Work

In this paper, we have proposed an SIRS-Vκ model, where $\kappa$ is the
vaccine confidence level. We have proved the existence of unique endemic and
disease-free states, both of which depend on $\kappa$. Furthermore, $\kappa$
acts as an important component in determining the necessary and sufficient
condition of the global asymptotic stability, namely
$(1-\kappa)\frac{\beta}{\gamma}\leq 1$, of the disease-free equilibrium. From
the perspective of control, manipulating the transmission rate $\beta$,
recovery rate $\gamma$, and the vaccine confidence $\kappa$ are equally viable
ways of mitigating epidemic spread. From the disease modeling perspective, we
have shown through analytical and numerical methods that ignoring vaccine
confidence in the model will introduce significant biases when determining the
system’s stability condition, maximum peak infection time, and its threshold
behaviors. The dependence of COVID-19 prevalence on vaccine hesitance should
encourage Americans to become vaccinated if they have not done so already.

For future work, we are interested in extending our findings from the scalar
case to the networked case to further inspect the impact of varying vaccine
confidence across different geographic/demographic populations. Moreover, we
also would like to provide analysis on the region of convergence of EEP and
DFE when $(1-\kappa)\frac{\beta}{\gamma}>1$, and the effect of vaccine
confidence on the maximum peaking time of $I(t)$.

## References

  * [1] M. J. Keeling and P. Rohani, _Modeling Infectious Diseases in Humans and Animals_. Princeton University Press, 2008\.
  * [2] K. Dietz and J. Heesterbeek, “Daniel Bernoulli’s epidemiological model revisited,” _Mathematical Biosciences_ , vol. 180, no. 1, pp. 1–21, 2002\.
  * [3] W. O. Kermack and A. G. McKendrick, “Contributions to the mathematical theory of epidemics. ii.—the problem of endemicity,” _Proceedings of the Royal Society of London. Series A, containing papers of a mathematical and physical character_ , vol. 138, no. 834, pp. 55–83, 1932.
  * [4] P. E. Paré, C. L. Beck, and T. Başar, “Modeling, estimation, and analysis of epidemics over networks: An overview,” _Annual Reviews in Control_ , vol. 50, pp. 345–360, 2020.
  * [5] G. Giordano, M. Colaneri, A. Di Filippo, F. Blanchini, P. Bolzern, G. De Nicolao, P. Sacchi, P. Colaneri, and R. Bruno, “Modeling vaccination rollouts, SARS-CoV-2 variants and the requirement for non-pharmaceutical interventions in Italy,” _Nature Medicine_ , vol. 27, no. 6, pp. 993–998, 2021.
  * [6] E. Jenner, _An inquiry into the causes and effects of the Variolae Vaccinae, a disease discovered in some of the western counties of England, particularly Gloucestershire, and known by the name of the cow pox_. printed, for the author, by Sampson Low, 1800\.
  * [7] M. R. Hilleman, “Vaccines in historic evolution and perspective: A narrative of vaccine discoveries,” _Vaccine_ , vol. 18, no. 15, pp. 1436–1447, 2000\.
  * [8] D. A. Henderson, “The eradication of smallpox–an overview of the past, present, and future,” _Vaccine_ , vol. 29, pp. D7–D9, 2011.
  * [9] “Immunization coverage.” [Online]. Available: https://www.who.int/news-room/fact-sheets/detail/immunization-coverage
  * [10] E. Dong, H. Du, and L. Gardner, “An interactive web-based dashboard to track COVID-19 in real time,” _The Lancet Infectious Diseases_ , vol. 20, no. 5, pp. 533–534, 2020.
  * [11] E. Mathieu, H. Ritchie, E. Ortiz-Ospina, M. Roser, J. Hasell, C. Appel, C. Giattino, and L. Rodés-Guirao, “A global database of COVID-19 vaccinations,” _Nature Human Behaviour_ , pp. 1–7, 2021.
  * [12] “Different COVID-19 vaccines.” [Online]. Available: https://www.cdc.gov/coronavirus/2019-ncov/vaccines/different-vaccines.html
  * [13] E. Dubé, C. Laberge, M. Guay, P. Bramadat, R. Roy, and J. A. Bettinger, “Vaccine hesitancy: An overview,” _Human Vaccines & Immunotherapeutics_, vol. 9, no. 8, pp. 1763–1773, 2013.
  * [14] N. E. MacDonald and the SAGE Working Group on Vaccine Hesitancy, “Vaccine hesitancy: Definition, scope and determinants,” _Vaccine_ , vol. 33, no. 34, pp. 4161–4164, 2015.
  * [15] C. Funk and A. Tyson, “Growing share of Americans say they plan to get a COVID-19 vaccine–or already have,” _Pew Research Center Science & Society.[(accessed on 2 April 2021)]_, 2021.
  * [16] P. Sowa, Ł. Kiszkiel, P. P. Laskowski, M. Alimowski, Ł. Szczerbiński, M. Paniczko, A. Moniuszko-Malinowska, and K. Kamiński, “COVID-19 vaccine hesitancy in Poland—multifactorial impact trajectories,” _Vaccines_ , vol. 9, no. 8, p. 876, 2021.
  * [17] B. L. C. A. d. Oliveira, M. A. G. Campos, R. C. d. S. Queiroz, B. F. d. Souza, A. M. d. Santos, A. A. M. d. Silva _et al._ , “Prevalence and factors associated with COVID-19 vaccine hesitancy in Maranhão, Brazil,” _Revista de Saúde Pública_ , vol. 55, 2021.
  * [18] J. Aw, J. J. B. Seng, S. S. Y. Seah, and L. L. Low, “COVID-19 vaccine hesitancy—a scoping review of literature in high-income countries,” _Vaccines_ , vol. 9, no. 8, p. 900, 2021.
  * [19] M. A. Pires and N. Crokidakis, “Dynamics of epidemic spreading with vaccination: Impact of social pressure and engagement,” _Physica A: Statistical Mechanics and its Applications_ , vol. 467, pp. 167–179, 2017.
  * [20] M. A. Pires, A. L. Oestereich, and N. Crokidakis, “Sudden transitions in coupled opinion and epidemic dynamics with vaccination,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2018, no. 5, p. 053407, 2018\.
  * [21] A. A. Dhondt, “Carrying capacity: a confusing concept.” _ACTA OECOL.(OECOL. GEN.)._ , vol. 9, no. 4, pp. 337–346, 1988.
  * [22] O. Diekmann, J. A. P. Heesterbeek, and J. A. Metz, “On the definition and the computation of the basic reproduction ratio $R_{0}$ in models for infectious diseases in heterogeneous populations,” _Journal of Mathematical Biology_ , vol. 28, no. 4, pp. 365–382, 1990.
  * [23] P. Van den Driessche and J. Watmough, “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission,” _Mathematical Biosciences_ , vol. 180, no. 1-2, pp. 29–48, 2002.
  * [24] P. Van den Driessche, “Reproduction numbers of infectious disease models,” _Infectious Disease Modelling_ , vol. 2, no. 3, pp. 288–303, 2017.