# Mapping Incidence and Prevalence Peak Data for SIR Forecasting Applications

Alexander C<EMAIL_ADDRESS>[ G. Casey
<EMAIL_ADDRESS>[ Lauren J.
<EMAIL_ADDRESS>[ Nishant
<EMAIL_ADDRESS>[ Lauren A.
<EMAIL_ADDRESS>[ Sara Del
<EMAIL_ADDRESS>[ Dave
<EMAIL_ADDRESS>[ Statistical Sciences
(CCS-6), Computer, Computational, and Statistical Sciences Division, Los
Alamos National Laboratory, Los Alamos, NM, 87545 presep=, ]e1,e2,e3,e7  Data
Analytics & Forecasting, A-1: Information Systems & Modeling, Los Alamos
National Laboratory, Los Alamos, NM, 87545 presep=, ]e5,e6  Information
Sciences (CCS-3), Theoretical Division, Los Alamos National Laboratory, Los
Alamos, NM, 87545 presep=, ]e4

###### Abstract

Infectious disease modeling and forecasting have played a key role in helping
assess and respond to epidemics and pandemics. Recent work has leveraged data
on disease peak infection and peak hospital incidence to fit compartmental
models for the purpose of forecasting and describing the dynamics of a disease
outbreak. Incorporating these data can greatly stabilize a compartmental model
fit on early observations, where slight perturbations in the data may lead to
model fits that project wildly unrealistic peak infection. We introduce a new
method for incorporating historic data on the value and time of peak incidence
of hospitalization into the fit for a Susceptible-Infectious-Recovered (SIR)
model by formulating the relationship between an SIR model’s starting
parameters and peak incidence as a system of two equations that can be solved
computationally. This approach is assessed for practicality in terms of
accuracy and speed of computation via simulation. To exhibit the modeling
potential, we update the Dirichlet-Beta State Space modeling framework to use
hospital incidence data, as this framework was previously formulated to
incorporate only data on total infections.

Compartmental Models,

Disease Forecasting,

Hospital Incidence,

Prevalence,

###### keywords:

, , , , , , and

##  1 Introduction

Compartmental models have seen broad usage at the onset of several disease
outbreaks in the last century as a means to project expected numbers of
infected and to drive healthcare response. Broadly speaking, a compartmental
model describes the dynamics of a disease spread by breaking a population into
set categories and modeling the process by which individuals transition
through these categories. Perhaps the most basic is the Kermack–McKendrick
model, often called the SIR or Susceptible-Infectious-Recovered model, which
models the movement of subjects from being Susceptible, to Infected (and
contagious), and then to Removed (Kermack, McKendrick and Walker, 1927). Since
the development of the SIR model, further research has extended the idea to
include additional compartments, such as the Exposed category in the SEIR
model – where a subject has been exposed to the disease but is not yet
infectious – and the ability to move back to the Susceptible category in the
SEIS model (see Walter and Contreras (1999); Hethcote (2000) for an overview
of each). Compartmental models have been applied to projection tasks such as
the 2014-15 Ebola epidemic (Chretien, Riley and George, 2015), the 2009 A/H1N1
influenza pandemic (Nsoesie et al., 2014), several HIV outbreaks (Anderson,
1988; Golub, Gorr and Gould, 1993; Nyabadza, Mukandavire and Hove-Musekwa,
2011), the recent COVID-19 pandemic (Zhao and Chen, 2020; Cooper, Mondal and
Antonopoulos, 2020; Zhang et al., 2022), and many other epidemiological
projection tasks in the last century (see, for instance, Guanghong et al.
(2004); LaDeau et al. (2011); Zhan et al. (2019) and references therein).

The elegance of compartmental models is in their succinct ability to describe
the state of a population in terms of how subjects transfer to and from the
different categories, and thus fitting these models requires estimation of
interpretable quantities such as rates of infection and recovery. For
instance, denote the proportion of the population in the susceptible,
infectious, and removed compartments by $S_{t},I_{t},$ and $R_{t},$
respectively, such that $S_{t}+I_{t}+R_{t}=1$ for all $t$. Then the SIR model
is determined by the equations

$\displaystyle\frac{dS_{t}}{dt}$ $\displaystyle=-\beta S_{t}I_{t},$ (1a)
$\displaystyle\frac{dI_{t}}{dt}$ $\displaystyle=\beta S_{t}I_{t}-\gamma
I_{t},$ (1b) $\displaystyle\frac{dR_{t}}{dt}$ $\displaystyle=\gamma I_{t},$
(1c)

where $\beta>0$ is the disease transmission rate and $\gamma>0$ is the rate of
recovery. If one knew these two rates, and the initial number of individuals
in each category – $S_{0},I_{0}$ and $R_{0}$ – the numbers
$S_{t},I_{t},\text{and}\;R_{t}$ could be numerically simulated for any time-
point $t$. While modern computational tools do make simulating these
quantities feasible, the need to simulate the entire system numerically to get
the exact counts in each compartment is challenging. That is, fitting these
models to data and can be computationally expensive for more intricate
compartmental structures. Researchers have determined analytic solutions for
the entire SIR model (explicit forms/approximations to the number of
susceptible, infected, and removed at a specific time) that do not require
numeric simulation (Harko, Lobo and Mak, 2014; Schlickeiser and Kröger, 2021;
Carvalho and Gonçalves, 2021). However, these solutions involve
reparameterizing the time axis, and explicit calculations back onto the
original time axis require numeric integration or approximation methods.
Similar approaches are also used to obtain exact and approximate solutions to
the more complicated compartmental models, such as the SEIR model (Wang, Wei
and Zhang, 2014; Piovella, 2020) and the SIRS model (Acedo, González-Parra and
Arenas, 2010).

Analytic maps from the initial starting parameters of an SIR curve to
quantities of interest (QoIs) – such as the value of the peak of the infected
curve and the limiting number of susceptible individuals – were developed as
early as the mid-1900s (Kendall, 1956; Bailey, 1957; Hethcote, 1976). The
focus on these quantities has generally been motivated by their usefulness to
making public health decisions. For instance, in papers such as Hethcote
(2000) and Weiss (2013), the maximum value of the infected compartment (not
the time of maximum infection) was studied for the purpose of informing public
health officials of the maximum number of infections after the initial
estimation of the disease transmission and recovery rates, since knowing this
maximum quantity informs how many hospital beds might be needed. The authors
in Castro et al. (2020) point out that the time of peak infection is also
informative for healthcare officials, and they develop an analytic form for
the peak time of infection using an approximation of the SIR curve. An exact
form for the peak infection time is available in Kendall (1956) and Deakin
(1975), albeit in terms of an integral without a closed form. Several modern
papers study fast approximations to the integral form for peak infection time
(Cadoni and Gaeta, 2020; Turkyilmazoglu, 2021), using the previously-mentioned
analytic solutions to the SIR model. These different approximations for peak
infection time are evaluated and compared in Kröger, Turkyilmazoglu and
Schlickeiser (2021).

Some recent papers study SIR curve QoIs for inferential tasks related to
modeling an epidemic (Miller, 2009; Lang et al., 2018). In Amaro (2023), the
Gumbel distribution is suggested as a good approximation of the infection
curve in an SIR model, and maps between the SIR curve parameters and peak
value/time are used to develop Method of Moments estimators of the Gumbel
parameters. The Gumbel distribution is then used to approximate the exact
solution of the SIR model. In Osthus et al. (2017), the relationship between
peak infection value and the SIR parameters is used to incorporate historical
data on epidemic peaks into the inferential task of fitting an SIR curve
during the early stages of an epidemic. The authors point out that an SIR
curve is sensitive to small perturbations in the transmission and recovery
rates, and that incorporating these data discourages models that well-
represent early data but drastically over-predict the peak infection value. In
McAndrew et al. (2024), both historic and surveyed QoIs are used to constrain
an SEIRH model (where the added “H” category refers to people hospitalized at
a given time) in much the same way as Osthus et al. (2017).

Figure 1: Two SIR models with incidence and starting values $S_{0}=0.9$ and
$I_{0}=R_{0}=0.05$. The plot on the left uses parameters
$(\beta,\gamma)=(1.137,0.446)$; the plot on the right uses $(2.592,1.058)$.

Both Amaro (2023) and Osthus et al. (2017) use information on peak infection
value and time to inform the modeling of a pandemic and epidemic,
respectively. Since most modern monitoring systems approximate the daily
number of new cases of a disease, called incidence, this is not the most
practically useful development. The number of infected individuals at any
given time, called prevalence, is typically an unobserved quantity (Noordzij
et al., 2010). For examples of all quantities in an SIR curve on a fixed time
axis, see Figure 1. Note that peak prevalence and peak incidence need not
occur at the same time, and that it is mathematically possible for the
incidence to be greater than prevalence (as visualized on the right side of
Figure 1). This being said, in most disease outbreaks of note, prevalence is
typically greater than incidence.

The model in McAndrew et al. (2024) incorporates incidence QoI data by
extending the parameter space to include peak incidence value (PIV) and peak
incidence time (PIT) and by defining a prior on these values using either
historic or survey data. While this model achieves a similar goal to the
ultimate goal of this paper, its main distinction is that that it requires an
importance sampling scheme to fit, which we outline here for completeness:
priors on the rate parameters and initial values for the compartmental model
are used to sample proposal values, the entire system is numerically simulated
to determine the (PIV, PIT), then the likelihood of this system is determined
via the prior probability of (PIV, PIT). Note that while this need not be true
for the SEIRH model considered in McAndrew et al. (2024) (since it is not an
SIR model), defining both a prior on the rate parameters and on PIV and PIT
for an SIR model is redundant given the initial values $S_{0}$, $I_{0}$, and
$R_{0}$, since a given set of starting parameters should immediately determine
these QoIs.

Incorporating historical peak values and times of an outbreak into a Bayesian
compartmental forecasting model via prior specifications is highly useful to
constrain forecasts111For a review of the differences between disease
forecasts and disease projections, see Massad et al. (2005)., especially early
in the outbreak (see Osthus et al. (2017), Osthus et al. (2019), and McAndrew
et al. (2024)). Relating these historical QoIs requires determining a
relationship between historical peak values and times as well as the
parameters and initial conditions of the compartmental model (e.g., the SIR
model). The challenge lies in this relationship. Analytically, this
relationship is known for observed peak prevalence value (PPV) and peak
prevalence time (PPT). However, almost all infectious disease data is of
incidence. Despite being principally unsound, the data application in Osthus
et al. (2017) treated incidence data as if it were prevalence in order to use
the known analytic relationships. In this paper, we make three contributions.
First, we develop the methods to map peak values and times of incidence to the
parameters of the SIR model, given the initial conditions. Second, we
demonstrate how to incorporate these new mappings into the model of Osthus et
al. (2017) and provide the code used to do so. Third, we compare the forecasts
of the modeling framework of Osthus et al. (2017) using both incidence data as
incidence, and using the misspecified prevalence data. In addition to
comparing these forecasts, we show that if SIR parameter inference is desired,
then mistaking incidence data for prevalence data will result in biased
estimation.

This paper is outlined as follows. In Section 2, we review existing methods
for mapping SIR parameters to peak/time of prevalence. In Section 3, we
develop identical maps for incidence, then provide methods for inverting these
maps for both incidence and prevalence. In Section 4, we develop and improve
upon the modeling framework in Osthus et al. (2017) to incorporate historic
incidence peak/time data when fitting an SIR curve. We apply this model to
influenza data in Section 5.

##  2 Analytic Solutions to the SIR Equations for Incorporating Prevalence
Data

The system (1a)-(1c) can be solved analytically for all time by
reparameterizing the time axis to be instead in terms of the number of
individuals removed from the system beyond the initial amount in the Removed
category $R_{0}$,

$\displaystyle S_{t}$ $\displaystyle=S_{0}e^{-\beta\tau_{t}},$ (2a)
$\displaystyle I_{t}$
$\displaystyle=S_{0}+I_{0}-S_{0}e^{-\beta\tau_{t}}-\gamma\tau_{t},$ (2b)
$\displaystyle R_{t}$ $\displaystyle=R_{0}+\gamma\tau_{t},$ (2c)

where $\tau_{t}$ is the inverse of the map

$t_{\tau}:=\int_{0}^{\tau}\left[\frac{d\tau^{\prime}}{S_{0}+I_{0}-S_{0}e^{\beta\tau^{\prime}}-\gamma\tau^{\prime}}\right]^{+},$
(3)

where $[x]^{+}:=\max\\{x,0\\}.$ The above form for the SIR dynamics is
particularly useful since it allows one to calculate the number of individuals
in each category without numerically simulating the entire system. The major
drawback of this form is that mapping values back to an interpretable time
requires one to approximate the integral in (3), since this integral is
nonelementary. There exist other methods that derive an analytic solution to
the SIR system (Harko, Lobo and Mak, 2014). We use this formulation since it
is well-known (having been originally developed in Kendall (1956)), and
because it simplifies the calculation of SIR curve QoIs, such as PPV and PPT.
For instance, the maximum of (2b) with respect to $\tau$ occurs at

$\tau^{*}=\frac{1}{\beta}\log\left(\frac{\beta S_{0}}{\gamma}\right),$ (4)

which can be derived via straightforward calculus. PPT can then be
approximated from (3), and PPV calculated from (2b).

The process above of calculating PPV and PPT, denoted
$(I_{t_{\tau^{*}}},t_{\tau^{*}}),$ from the initial starting parameters can be
reversed for fixed $S_{0},I_{0},R_{0}$. That is, given
$(I_{t_{\tau^{*}}},t_{\tau^{*}})$, we propose a method to calculate
$(\beta,\gamma)$. Define $\rho$ as

$\rho=\frac{\beta}{\gamma}.$ (5)

The quantity $S_{0}\rho$ is more widely known as the basic reproduction
number, which measures the average number of additional infections generated
by a single new infection (Cadoni and Gaeta, 2020). Plugging (4) into (2b)
gives

$I_{t}=S_{0}+I_{0}-\frac{1}{\rho}-\frac{\log(\rho S_{0})}{\rho}.$ (6)

From here, setting $I_{t}:=I_{t_{\tau^{*}}}$ and solving for $\rho$ is a
(piecewise) convex optimization problem, which is typically fast and accurate
computationally. Note that a change of variables for (3) gives

$\beta
t_{\tau}=\int_{0}^{\beta\tau}\left[\frac{d\hat{\tau}}{S_{0}+I_{0}-S_{0}e^{-\hat{\tau}}-\hat{\tau}/\rho}\right]^{+},$
(7)

and that $\beta\tau^{*}=\log(\rho S_{0})$. Thus, once the value of $\rho$ is
approximated, (7) can be evaluated to get $\beta t_{\tau^{*}}$. Since
$t_{\tau^{*}}$ is assumed known a priori, the values of $\beta$ and $\gamma$
can be obtained via arithmetic.

We introduce the above derivations not only because they will be used in the
next Section, but also to motivate the type of calculation we aim to develop
in this paper. While the above approach still requires computational methods
to map from PPV and PPT to the SIR model parameters, this approach is much
more direct than the brute force method of simulating several $(\beta,\gamma)$
combinations until a prevalence curve with a peak sufficiently near
$(I_{t_{\tau^{*}}},t_{\tau^{*}})$ is discovered (Prangle, 2016).

##  3 Mapping Incidence Data to SIR Parameters

As mentioned above, a primary focus of this work will be to develop maps
between PIV and PIT and the initial SIR curve parameters. In this direction,
one can derive an equivalent formulation of the SIR dynamics in (1a) - (1c) by
replacing (1b) with

$\displaystyle I_{t}$ $\displaystyle=I_{t-1}+\beta S_{t-1}I_{t-1}-\gamma
I_{t-1}.$ (8)

Putting (8) into words, prevalence is equal to the prevalence at the last time
step, plus those in the infected category that infect those in the susceptible
category, minus the number in the infected category that are removed from the
system. The form for prevalence in (8) is particularly useful for this
application since it contains a term that explicitly models incidence at time
$t$: $\beta S_{t-1}I_{t-1}$.

###  3.1 Mapping SIR Parameters to Peak Incidence Value and Time

Reparameterizing the time axis for the term $\beta S_{t}I_{t}$ using (2a) and
(2b) gives the following form

$\beta(S_{0}e^{-\beta\tau_{t}})(S_{0}+I_{0}-S_{0}e^{-\beta\tau_{t}}-\gamma\tau_{t}).$
(9)

The value for $\tau$ that maximizes (9) also satisfies the equation,

$-(S_{0}+I_{0})+\gamma\tau+2S_{0}e^{-\beta\tau}-\frac{1}{\rho}=0.$ (10)

Solving (10) is also a convex optimization problem on a single parameter, and
is thus feasible and accurate to do numerically. A solution to (10) gives the
$\tau$ for the timepoint directly before the time of max incidence,
$\tau_{t^{*}-1}$, where $t^{*}$ denotes the timepoint of max incidence. Of
course, this is all that is required to calculate the PIV, $\beta
S_{t-1}I_{t-1},$ by a direct application of (9). The calculation of PIT is
similarly straightforward. Using (2b), the prevalence at $\tau_{t^{*}-1}$ can
be directly calculated. The prevalence value for $\tau_{t^{*}-1}$ and (2c)
then gives $\tau_{t^{*}},$ and (3) can then be used to calculate $t^{*}$,

$t^{*}=\int_{0}^{\tau_{t^{*}}}\left[\frac{d\tau^{\prime}}{S_{0}+I_{0}-S_{0}e^{\beta\tau^{\prime}}-\gamma\tau^{\prime}}\right]^{+}.$

###  3.2 Mapping Incidence Peak Value and Time to SIR Parameters

The proposed method to map PIV and PIT back to the parameters of an SIR curve
solves the following system of equations implied by (9) and (10):

$\displaystyle\beta(S_{0}e^{-\beta\tau^{\star}})(S_{0}+I_{0}-S_{0}e^{-\beta\tau^{\star}}-\gamma\tau^{\star})$
$\displaystyle=\text{PIV}$ (11)
$\displaystyle-(S_{0}+I_{0})+\gamma\tau^{\star}+2S_{0}e^{-\beta\tau}-\frac{1}{\rho}$
$\displaystyle=0,$ (12)

where $\tau^{\star}$ solves the equation

$\int_{0}^{\tau^{\star}}\left[\frac{d\tau^{\prime}}{S_{0}+I_{0}-S_{0}e^{\beta\tau^{\prime}}-\gamma\tau^{\prime}}\right]^{+}=\text{PIT}.$
(13)

While solving a system of two equations with two unknowns (for $\beta$ and
$\gamma$) is generally feasible computationally, the major bottleneck for this
problem is the need to invert (13). A brute-force computational approach to
solving this system of equations would require one to both invert and solve
(13) for every value of $(\beta,\gamma)$ investigated. While experiments in
this direction have proven to be surprisingly fast and accurate, the
confounding computational approximations encourage a more analytic solution,
or alternative computational strategies. We consider three possible
alternatives to computationally estimating the integral (referred to the
“Compute Integral” method in the subsequent).

####  3.2.1 Taylor Approximation

Given a candidate $(\beta,\gamma)$ in any numerical solver, we approximate
$\tau^{\star}$ by taking the second degree Taylor expansion of the integral in
(13) and solving for $\tau^{\star}$ algebraically (Murray, 2002). This leads
to the following closed-form approximation,

$\tau^{\star}=\frac{\beta^{2}}{S_{0}}\left[\left(\rho
S_{0}-1\right)+\kappa\tanh\left(\frac{\gamma\kappa(\text{PIT})}{2}-\phi\right)\right]+R_{0},$
(14)

where

$\displaystyle\kappa$
$\displaystyle=\sqrt{\left(S_{0}\rho-1\right)^{2}+2S_{0}I_{0}\rho^{2}},$
$\displaystyle\phi$
$\displaystyle=\frac{1}{\kappa}\operatorname{arctanh}\left[S_{0}\rho-1\right].$

Using this approximation for $\tau^{\star}$, we numerically solve the system
of equations expressed by (11) and (12).

####  3.2.2 Single ODE Approximation

Instead of using a Taylor approximation to estimate $\tau^{\star}$ it is
possible to do so by numerically solving an Ordinary Differential Equation
(ODE). Using the definition of $\tau$, we combine equations (2a) and (2c) to
get

$S_{t}=S_{0}\exp\left(-\rho(R_{t}-R_{0})\right).$ (15)

From here, we use (1c) and the assumption that $I_{t}=1-S_{t}-R_{t}$ to get
the ODE,

$\frac{dR_{t}}{dt}=\gamma\left[1-S_{0}e^{-\rho(R_{t}-R_{0})}-R_{t}\right].$
(16)

The authors in Cadoni and Gaeta (2020) point out that while (16) is a
transcendental equation, it can be solved numerically, and has a single,
unique solution by the general existence and uniqueness theorem for the
solutions of ODEs. Using this approximation for $\tau^{\star}$, we numerically
solve the system of equations expressed by (11) and (12).

####  3.2.3 Full ODE Approximation

As a final computational method for mapping PIV and PIT to $\beta$ and
$\gamma$, we numerically solve the system of ODEs described in (1a) - (1c). In
an optimization algorithm, this would require the ODE to be solved for every
possible $(\beta,\gamma)$ pair queried. Much like the brute force
computational method, we expect this method to be accurate, but to come at a
high computational cost.

####  3.2.4 Comparision via Simulation

We compare all the above approaches via simulation, repeating the following
steps 1000 times:

  1. 1.

Sample a (PPV, PPT) from the bivariate normal $\mathcal{N}(\mu,\Sigma)$ where
$\mu=(0.0144,17.9)$,

$\Sigma=\begin{pmatrix}0.000036&0.0187\\\ -0.0187&16.09\end{pmatrix},$

truncated so that PPV $\in(\theta_{I_{0}},1)$ and PPT $\in(1,35)$. This
corresponds to the set of feasible values and sampling distribution described
in Osthus et al. (2017);

  2. 2.

Use the method from Section 2 to map this back to $(\beta,\gamma)$, then
numerically simulate the system in (1a)-(1c) to get the “true” values for PIT
and PIV;

  3. 3.

With PIT and PIV, use each method described above to find approximate values
$(\hat{\beta},\hat{\gamma})$;

  4. 4.

For all approximation methods, compare the estimated $\widehat{\text{PIT}}$
and $\widehat{\text{PIV}}$ (gotten by numerically simulating the system from
the appropriate $(\hat{\beta},\hat{\gamma})$) against the true PIT and PIV.

We outline the results from the above simulation in Table 1. While the fastest
method is the Taylor approximation, this method is also the least accurate.
This is as expected, since this approximation is only accurate for
sufficiently small values of $\rho(R_{t}-R_{0})$ (Murray, 2002). The Single
ODE approximation method is slightly more accurate than the Taylor
approximation, but at a higher computational cost. The Compute Integral and
the Full ODE approximation are comparable, with the Compute Integral method
being more accurate and the Full ODE method being faster.

While the Compute Integral method is the slowest of these approaches, it is
still somewhat fast (taking around half a second on average), and it is the
most accurate overall (since the Taylor approximation and Single ODE methods
are greatly off for PIV). Since the data sizes for PIT and PIV data are not
exorbitantly large, this computational burden would be acceptable in
converting a data set of $(\text{PIV},\text{PIT})$ values to
$(\hat{\beta},\hat{\gamma})$ values. For the applications in this paper, we
will use the Compute Integral method. Future work might examine alternative
approximation methods, especially for approximating the inverse of the
integral in (13).

Quantity | Compute Integral | Taylor Approx. | Single ODE | Full ODE
---|---|---|---|---
Avg. PIV Error | 4.07e-4 | 1.51e-4 | 1.36e-4 | 16.02e-4
Std. Dev. PIV Error | 1.836e-4 | 0.566e-4 | 0.660e-4 | 6.848e-4
Avg. PIT Error | 0 | 3.209 | 2.236 | 0
Std. Dev. PIT Error | 0 | 1.888 | 1.647 | 0
Avg. Runtime | 5.626e-1 | 0.035e-1 | 2.017e-1 | 0.605e-1
Std. Dev. Runtime | 2.060e-1 | 0.015e-1 | 0.510e-1 | 0.324e-1
Table 1: Errors and runtimes (in seconds) for the four approximation methods
described in Section 3.2, using the simulation study described in Section
3.2.4. The methods are the Compute Integral approach (top of 3.2), the Taylor
approximation (3.2.1), the single ODE approximation (3.2.2), and the full ODE
approximation (3.2.3).

##  4 A Bayesian State-Space SIR Model

In this section, we introduce the Dirichlet-Beta state-space model (DBSSM)
from Osthus et al. (2017) and update it to incorporate historic PIV and PIT
data. The original formulation of the DBSSM was to answer an observed issue
associated with using the SIR model for early-pandemic forecasting tasks.
Namely, that two SIR curves that reasonably fit early count data may lead to
drastically different PPV predictions. In a simulation example, Osthus et al.
(2017) show two such SIR curves that have peaks that differ by 30% of the
entire population, even though they have a nearly-identical fit to the early-
pandemic data observations (see Figure 3 in the cited paper). To address this
stability issue, the DBSSM incorporates historic PPV and PPT data into the
prior specifications to discourage SIR curve fits with peak values that are
greatly above reasonable expectations. As mentioned previously, this
incorporation of prevalence data is not the most practical approach, since
incidence data is generally the observed quantity. A further shortcoming in
the original DBSSM formulation is that it learns a map between PPT and the SIR
curve parameters, rather than using an analytic map. After introducing this
model, we will identify ways that the methodology developed in this paper will
improve these issues for the DBSSM.

Let $y_{t}$ be the observed proportion in a population that tested positive
for some disease at some timepoint $t$, and let
$\theta_{t}=(S_{t},I_{t},R_{t})^{\prime}$. Then the DBSSM is defined as,

$y_{t}|\theta_{t},\phi\sim\text{Beta}\left(\lambda\text{In}_{t},\lambda(1-\text{In}_{t})\right)$
(17) $\theta_{t}|\theta_{t-1},\phi\sim\text{Dirichlet}\left(\iota
f(\theta_{t-1},\beta,\gamma)\right),$ (18)

where $\phi=\\{S_{0},I_{0},R_{0},\beta,\gamma,\lambda,\iota\\}$,
$\lambda,\iota$ are variance control parameters, and $f$ is a map that
propagates the SIR system determined by $(\theta_{t-1},\beta,\gamma)$ forward
one step according to (1a)-(1c). Note that, by this set up, the set of
parameters
$\theta_{0:t^{\prime}}=\\{\theta_{0},\theta_{1},\dots,\theta_{t^{\prime}}\\}$
is a first-order Markov chain, and that for all $s\neq t$, the data
observations $y_{s}$ and $y_{t}$ are independent given $\theta_{t}$. The
variable $\text{In}_{t}$ denotes the incidence of the system at time $t$; in
the original formulation of this model, the prevalence – $I_{t}$ – was used
here. The incidence at time $t$ is directly calculated using $\theta_{t}$ and
(8).

The conditional expectations of the model described by (17) and (18) are
unbiased,

$E(y_{t}|\theta_{t},\phi)=\text{In}_{t}$ (19)
$E(\theta_{t}|\theta_{t-1},\phi)=f(\theta_{t-1},\beta,\gamma)$ (20)

while their respective variances reduce to zero as $\lambda,\iota\to\infty$.
Of course, the conditional mean in (20) is dependant upon the accuracy of $f$
in propogating the latent space $\theta_{t-1}$ forward one time step. The
authors in Osthus et al. (2017) used a fourth-order Runge-Kutta approximation
and observed reasonable accuracy.

We review the full Bayesian framework of the DBSSM and provide the prior
specification in Appendix A. The main innovation of this model is that the
parameter space is expanded to include the latent variable $z=(PPT,PPV)$, and
this latent variable is given a prior that incorporates historic PPT and PPV
data. In the following, we review this prior, $\pi(z|\theta_{0})$, and the
mechanism by which this prior informs the SIR parameters $\beta$ and $\gamma$.
These priors are then each updated according to the theory developed in this
paper.

###  4.1 Specification of $\pi(z|\theta_{0})$

The prior on $\pi(z|\theta_{0})$ in the DBSSM is a minimal-assumption
distribution on historic data on PPT and PPV. Note that this is the mechanism
by which the authors in Osthus et al. (2017) directly address the
aforementioned stability issues with fitting an SIR curve with early pandemic
data. They do so by fitting a normal distribution to historic influenza QoI
data, truncated to enforce that an epidemic will occur (the lower bound on
$PPV$ was set to $I_{0}$) and so that the peak happens within the influenza
forecasting season ($PPT$ was required to be between the $1^{st}$ and
$35^{th}$ weeks). While somewhat loose, this prior gives very small (or zero)
probability to values of $(PPV,PPT)$ that are drastically outside of
historically observed pandemics.

In the formulation of the DBSSM developed in this paper, the same prior used
on PPT and PPV is now used on PIT and PIV. To connect this constraint into the
model, we must next define how the assumption on this latent space affects the
learning of the SIR parameters $\beta$ and $\gamma$.

###  4.2 Specification of $\pi(\beta,\gamma|z,\theta_{0})$

With the addition of the latent variable $z$, the prior needed for the SIR
parameters is $\pi(\beta,\gamma|z,\theta_{0})$. In the original formulation of
the DBSSM, this prior is reparameterized according to $(\rho,\gamma)$, then
factorized. Thus, priors are instead given to $\pi(\rho|z,\theta_{0})$ and
$\pi(\gamma|\rho,z,\theta_{0})$. This additional formulation is done to
utilize the following analytic relationship from Weiss (2013),

$PPV=g_{1}(S_{0},I_{0},\rho)=I_{0}+S_{0}-S_{0}\rho\left[\log(S_{0})+1-\log(S_{0}\rho)\right].$
(21)

By inverting this relationship, samples from the latent quantity $z$
immediately determine the corresponding value of $\rho$. The appropriate prior
on this quantity would then be

$\pi(\rho|z,\theta_{0})~{}\propto~{}\delta(S_{0}\rho-
g_{1}^{-1}(PPV,S_{0},I_{0})),$ (22)

where $\delta$ is the Dirac delta function. For the prior on $\gamma$, the map
between $\gamma$ and PPT is estimated using a simulated data set of 5250 SIR
curves. This map,

$PPT=g_{2}(S_{0},I_{0},\rho,\gamma),$ (23)

was then used in lieu of any analytic form, and the prior on $\gamma$ was set
to

$\pi(\rho|z,\theta_{0})~{}\propto~{}\delta(S_{0}\rho-
g_{2}^{-1}(PPV,S_{0},I_{0},\rho)),$ (24)

We reiterate that there are two major shortcomings to the above prior
specifications. First, the above priors assume that there is access to
historical PPV and PPT data, which is typically not the case, as public health
data are generally on incidence, not prevalence. In the original paper,
incidence data were used instead of prevalence data without explicit
justification. Second, a map between PPT and the SIR parameters is estimated
even though an analytic map between these quantities exists (see (3) and (4)),
unnecessarily introducing a source of uncertainty.

The methods developed in this paper correct the limitations found in the
original formulation of the DBSSM, since they provide maps to replace
$g_{1},g_{2}$ above with

$(PIV,PIT)=h(S_{0},I_{0},\beta,\gamma),$ (25)

where $h^{-1}$ denotes the algorithm discussed in Section 3.2. Thus, the joint
prior used for $\beta,\gamma$ in this updated version of the DBSSM is

$\pi(\beta,\gamma|z,\theta_{0})~{}\propto~{}\delta\left(||(\beta,\gamma)-h^{-1}(PPV,S_{0},I_{0})||_{1}\right),$
(26)

where $||\cdot||_{1}$ is the 1-norm.

The naive treatment of incidence data as prevalence data (as was done in
Osthus et al. (2017)), need not necessarily lead to a loss of forecasting
accuracy in the final model. An incidence curve can well approximate, or even
be equivalent to, a prevalence curve. For instance, consider the SIR model
constrained so that $\gamma=1$, which corresponds to the Reed-Frost model
(Abbey, 1952). In this case, incidence is precisely equal to prevalence, and
thus either method for incorporating historic QoI data should yield an
equivalent model. This insight leads to the following Remark:

######  Remark 1.

Using incidence QoI data in place of prevalence QoI data naively leads to an
SIR curve where the $I$ compartment – which normally corresponds to prevalence
– now models the progression of incidence. Using historic incidence QoI data
in the way outlined in this manuscript uses the incidence curve to model
incidence as is desired. While either can be viable for the purposes of
forecasting, only the method described in this paper leads to rate parameter
estimates that can be interpreted as infection and recovery rates.

##  5 Application to Seasonal Influenza Data

We recreate the data application from Osthus et al. (2017), using the updated
model from Section 4. The aim of this application is not to improve the
forecasting in the original formulation of the DBSSM (see Remark 1). Rather,
we will demonstrate that the forecasting capabilities of this model remain the
same, while we also observe different estimates for the infection rate
$\beta$, the recovery rate $\gamma$ and the basic reproduction number $\rho$.

The source data modified and then used for this application are counts of
patients seen in the US with an influenza-like illness (ILI), where ILI is
defined as having a temperature of at least 100 degrees Fahrenheit, a cough
and/or a sore throat, and no known cause for those symptoms other than
influenza (CDC, 2024). These data are collected weekly, where more than 3400
outpatient healthcare providers report to the CDC the number of patients with
ILI they treated (CDC, 2023).

The number of patients reported as having ILI will naturally also include
cases of respiratory illnesses other than influenza. Following the approach of
Shaman et al. (2013), we use virologic surveillance data (where patients are
actually tested for influenza) to estimate the proportion of ILI patients with
influenza, then multiply ILI data by this proportion. This corrected data is
referred to as ILI+. For more details on this adjustment, see Shaman et al.
(2013). Note that the ILI+ data estimates the weekly incidence of influenza
cases – not prevalence.

Figure 2: The DBSSM fit to the 2010 US nationwide influenza outbreak starting
at week 13 (left column) and week 22 (right column). The top plots reproduces
the forecasting model of Osthus et al. (2017), which incorrectly treated
incidence data as though it were prevalence data. The bottom plots fit the
forecasting model described in Section 3, which correctly treats incidence
data as incidence data. Earlier on in the outbreak, the model developed in
this paper has a tighter prediction interval. The forecasts starting at week
22 are quite similar.

To fit both versions of the DBSSM, we use ten influenza seasons: the seasons
that started in the years 2002-2007, and the seasons the started during
2010-2013. The years 2008 and 2009 were omitted to be consistent with Osthus
et al. (2017); these two years correspond to a pandemic and the focus of that
work was to forecast seasonal influenza. Each season is defined as 35
consecutive weeks starting on roughly the first week of October (epidemiology
week 40, treated as $t=1$).

For a estimated proportion of individuals in a population infected with
influenza at timepoints $\mathbb{T}=\\{1,\dots,T\\}$, suppose only the ILI+
data up through $t^{\prime}\in\mathbb{T}$ are observed. Given this, we
simulate 62500 from the posterior
$\pi(\boldsymbol{\theta_{{1:t^{\prime}}}},\boldsymbol{\phi}|y_{1:t^{\prime}})$
for four separate chains, discarding the first 12500 as burn-in and thinning
out all but every tenth observation in the remaining samples. Given these
draws from the posterior distribution, the posterior predictive density,
$\pi(y_{(t^{\prime}+1):T}|y_{1:t^{\prime}})$, is used to estimate “future”
observations of ILI+ data.

We perform two separate fits of this posterior model, on the first thirteen
days ($t^{\prime}=13$) and on the first twenty two days ($t^{\prime}=22$), for
the considered ILI+ data for the 2010 influenza season in the United States.
These fits are performed both using the original formulation of the DBSSM,
which naively uses incidence data directly in place of prevalence data, and
using the new formulation developed in this paper, which uses the new maps
developed in Section 3 for a more principled treatment of incidence data.
These fits and forecasts are outlined in Figure 2. The dark shaded grey
regions prior to $t^{\prime}$ mark the 95 percentiles of the posterior
density, while the lighter grey shaded regions after $t^{\prime}$ make the 95%
prediction intervals. The forecast using incidence data up until
$t^{\prime}=13$ has a narrower prediction interval than the one using
prevalence data; each of the forecasts that use data up through
$t^{\prime}=22$ are laregly comparable.

Parameter |  $\beta$ |  $\gamma$ |  $S_{0}\beta/\gamma$
---|---|---|---
Prevalence Median |  2.15 |  1.60 |  1.21
Incidence Median |  3.21 |  2.66 |  1.09

Figure 3: Infection rates, recovery rates, and reproductions numbers drawn
from the Gibbs sampler used to fit the DBSSM using both specifications of
incidence data. Using incidence data specified as prevalence leads to
different estimates for these parameters. These values were calculated using
the data up through timepoint 22 ($t^{\prime}=22$).

In addition to the slight improvements on forecasting we observe in Figure 2,
this method also has strong implications for the interpretability of
$\beta,\gamma$ for the fitted model. Indeed, only the updated version of the
DBSSM developed in this paper leads to realizations of these parameters that
can accurately be interpreted as the infection rate ($\beta$), recovery rate
($\gamma$), and the basic reproduction number ($S_{0}\rho$).

##  6 Discussion

The main contribution of this paper is the development of methods to map the
time and value of peak incidence to the SIR curve parameters, and vice versa,
for the purpose of forecasting tasks and inference on disease rate parameters.
We do this by computationally solving a system of equations ((9) and (10)).
There are several impactful uses of these maps in the context of previous
literature. First, much like how the peak prevalence value (PPV) and time
(PPT) are useful for public health response to an epidemic (Weiss, 2013), the
analogous quantities for incidence are also useful, since they describe the
influx new patients entering the hospital system on a given day. Second, this
work improves upon existing work that uses historical prevalence data to model
epidemics by creating a map from PIT and PIV to the SIR parameters, since
incidence is typically the data that is available for ongoing epidemics
(Osthus et al., 2017; Amaro, 2023). In the case of the application in Osthus
et al. (2017), where incidence data were used in place of prevalence data
without justification, we have shown that this leads to biased SIR parameter
estimates (see Figure 3). Furthermore, our results indicate that forecasts
performed using the erroneous data specification leads to larger prediction
intervals earlier on in the outbreak, although this forecast is largely
comparable for the correct specification later on in the outbreak (see Figure
2 and Remark 1). Of course, it remains more correct in principle to use
incidence data appropriately when fitting compartment models for forecasting
with ongoing incidence data (Nsoesie, Mararthe and Brownstein, 2013; Chowell
et al., 2016; Abolmaali and Shirzaei, 2021). We have provided a modeling
framework that incorporates this data appropriately (see Section 4). Lastly,
while the methods discussed in McAndrew et al. (2024) do incorporate incidence
data to fit forecasting models, this paper uses a Bayesian framework and
importance sampling. Since the maps developed here are deterministic, they can
be used in both a Bayesian and a Frequentist framework.

As a direction for future work, it would be useful to investigate better
approximations for the solution to (13). The Taylor Approximation in Section
3.2.1 was by far the fastest computationally, but it came with the highest
error on PIT. Finding a fast and accurate approximation to this equation would
greatly increase the runtime for applications where the map between the SIR
parameters and PIT/PIV must be evaluated several hundreds of thousands of
times. However, for most applications (including the one in this paper), the
Compute Integral approximation is sufficiently fast.

As a second direction for future work, it would be interesting to investigate
the analogous maps for more complicated compartmental models, such as the SEIR
model, the SIRS model, and the SEIRH model.

## Open Research Section

All software used to perform the simulations and studies in this paper are
publicly available at https://github.com/lanl/precog.

## Acknowledgments

Research presented in this article was partially supported by the Laboratory
Directed Research and Development program of Los Alamos National Laboratory
under project number 20240066DR. Los Alamos National Laboratory is operated by
Triad National Security, LLC, for the National Nuclear Security Administration
of U.S. Department of Energy (Contract No. 89233218CNA000001).

This research was partially funded by NIH/NIGMS under grant R01GM130668-01
awarded to Sara Y. Del Valle.

## Conflict of Interest

The authors have no non-financial or other financial competing interests to
declare that are relevant to the content of this article other than the
aforementioned declared funding sources.

## References

  * Abbey (1952) [author] Abbey, HelenH. (1952). AN EXAMINATION OF THE REED-FROST THEORY OF EPIDEMICS. Human Biology 24 201–233.
  * Abolmaali and Shirzaei (2021) [author] Abolmaali, SainaS. and Shirzaei, SamiraS. (2021). A comparative study of SIR Model, Linear Regression, Logistic Function and ARIMA Model for forecasting COVID-19 cases. AIMS Public Health 8 598–613. 10.3934/publichealth.2021048
  * Acedo, González-Parra and Arenas (2010) [author] Acedo, L.L., González-Parra, GilbertoG. and Arenas, Abraham J.A. J. (2010). An exact global solution for the classical SIRS epidemic model. Nonlinear Analysis: Real World Applications 11 1819-1825. https://doi.org/10.1016/j.nonrwa.2009.04.007
  * Amaro (2023) [author] Amaro, J. E.J. E. (2023). Systematic description of COVID-19 pandemic using exact SIR solutions and Gumbel distributions. Nonlinear Dynamics 111 1947–1969. 10.1007/s11071-022-07907-4
  * Anderson (1988) [author] Anderson, Roy MR. M. (1988). The role of mathematical models in the study of HIV transmission and the epidemiology of AIDS. Journal of Acquired Immune Deficiency Syndromes 1 241–256.
  * Bailey (1957) [author] Bailey, Norman T. J.N. T. J. (1957). The mathematical theory of epidemics. Griffin.
  * Cadoni and Gaeta (2020) [author] Cadoni, MarianoM. and Gaeta, GiuseppeG. (2020). Size and timescale of epidemics in the SIR framework. Physica D: Nonlinear Phenomena 411 132626\. https://doi.org/10.1016/j.physd.2020.132626
  * Carvalho and Gonçalves (2021) [author] Carvalho, Alexsandro M.A. M. and Gonçalves, SebastiánS. (2021). An analytical solution for the Kermack–McKendrick model. Physica A: Statistical Mechanics and its Applications 566 125659\. https://doi.org/10.1016/j.physa.2020.125659
  * Castro et al. (2020) [author] Castro, MarioM., Ares, SaúlS., Cuesta, José A.J. A. and Manrubia, SusannaS. (2020). The turning point and end of an expanding epidemic cannot be precisely forecast. Proceedings of the National Academy of Sciences 117 26190-26196. 10.1073/pnas.2007868117
  * CDC (2023) [author] CDC, (2023). U.S. Influenza Surveillance: Purpose and Methods. data retrieved from Centers for Disease Control and Prevention, National Center for Immunization and Respiratory Diseases (NCIRD), https://www.cdc.gov/flu/weekly/overview.htm.
  * CDC (2024) [author] CDC, (2024). Glossary of Influenza (Flu) Terms. data retrieved from Centers for Disease Control and Prevention, National Center for Immunization and Respiratory Diseases (NCIRD), https://www.cdc.gov/flu/about.
  * Chowell et al. (2016) [author] Chowell, GerardoG., Sattenspiel, LisaL., Bansal, ShwetaS. and Viboud, CécileC. (2016). Mathematical models to characterize early epidemic growth: A review. Phys Life Rev 18 66–97. 10.1016/j.plrev.2016.07.005
  * Chretien, Riley and George (2015) [author] Chretien, Jean-PaulJ.-P., Riley, StevenS. and George, Dylan BD. B. (2015). Mathematical modeling of the West Africa Ebola epidemic. Elife 4 e09186.
  * Cooper, Mondal and Antonopoulos (2020) [author] Cooper, IanI., Mondal, ArghaA. and Antonopoulos, Chris GC. G. (2020). A SIR model assumption for the spread of COVID-19 in different communities. Chaos Solitons Fractals 139 110057\. 10.1016/j.chaos.2020.110057
  * Deakin (1975) [author] Deakin, Michael A. B.M. A. B. (1975). A standard form for the Kermack-McKendrick epidemic equations. Bulletin of Mathematical Biology 37 91-95. https://doi.org/10.1016/S0092-8240(75)80011-5
  * Golub, Gorr and Gould (1993) [author] Golub, AndrewA., Gorr, Wilpen LW. L. and Gould, Peter RP. R. (1993). Spatial diffusion of the HIV/AIDS epidemic: modeling implications and case study of AIDS incidence in Ohio. Geographical analysis 25 85–100.
  * Guanghong et al. (2004) [author] Guanghong, DingD., Chang, LiuL., Jianqiu, GongG., Ling, WangW., Ke, ChengC. and Di, ZhangZ. (2004). SARS epidemical forecast research in mathematical model. Chinese Science Bulletin 49 2332–2338.
  * Harko, Lobo and Mak (2014) [author] Harko, TiberiuT., Lobo, Francisco S. N.F. S. N. and Mak, M. K.M. K. (2014). Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates. Applied Mathematics and Computation 236 184-194. https://doi.org/10.1016/j.amc.2014.03.030
  * Hethcote (1976) [author] Hethcote, Herbert W.H. W. (1976). Qualitative analyses of communicable disease models. Mathematical Biosciences 28 335-356. https://doi.org/10.1016/0025-5564(76)90132-2
  * Hethcote (2000) [author] Hethcote, Herbert WH. W. (2000). The mathematics of infectious diseases. SIAM review 42 599–653.
  * Kendall (1956) [author] Kendall, David GD. G. (1956). Deterministic and stochastic epidemics in closed populations. In Proceedings of the third Berkeley symposium on mathematical statistics and probability 4 149–165. University of California Press Berkeley.
  * Kermack, McKendrick and Walker (1927) [author] Kermack, William OgilvyW. O., McKendrick, A. G.A. G. and Walker, Gilbert ThomasG. T. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 115 700-721. 10.1098/rspa.1927.0118
  * Kröger, Turkyilmazoglu and Schlickeiser (2021) [author] Kröger, MartinM., Turkyilmazoglu, MustafaM. and Schlickeiser, ReinhardR. (2021). Explicit formulae for the peak time of an epidemic from the SIR model. Which approximant to use? Physica D: Nonlinear Phenomena 425 132981\. https://doi.org/10.1016/j.physd.2021.132981
  * LaDeau et al. (2011) [author] LaDeau, Shannon LS. L., Glass, Gregory EG. E., Hobbs, N ThompsonN. T., Latimer, AndrewA. and Ostfeld, Richard SR. S. (2011). Data–model fusion to better understand emerging pathogens and improve infectious disease forecasting. Ecological Applications 21 1443–1460.
  * Lang et al. (2018) [author] Lang, John CJ. C., De Sterck, HansH., Kaiser, Jamieson LJ. L. and Miller, Joel CJ. C. (2018). Analytic models for SIR disease spread on random spatial networks. Journal of Complex Networks 6 948-970. 10.1093/comnet/cny004
  * Massad et al. (2005) [author] Massad, EduardoE., Burattini, Marcelo NM. N., Lopez, Luis FL. F. and Coutinho, Francisco A BF. A. B. (2005). Forecasting versus projection models in epidemiology: the case of the SARS epidemics. Med Hypotheses 65 17–22. 10.1016/j.mehy.2004.09.029
  * McAndrew et al. (2024) [author] McAndrew, ThomasT., Gibson, Graham C.G. C., Braun, DavidD., Srivastava, AbhishekA. and Brown, KateK. (2024). Chimeric Forecasting: An experiment to leverage human judgment to improve forecasts of infectious disease using simulated surveillance data. Epidemics 47 100756\. https://doi.org/10.1016/j.epidem.2024.100756
  * Miller (2009) [author] Miller, Joel CJ. C. (2009). Spread of infectious disease through clustered populations. J R Soc Interface 6 1121–1134. 10.1098/rsif.2008.0524
  * Murray (2002) [author] Murray, James DicksonJ. D. (2002). Mathematical Biology: An introduction. Springer.
  * Noordzij et al. (2010) [author] Noordzij, MarliesM., Dekker, Friedo WF. W., Zoccali, CarmineC. and Jager, Kitty JK. J. (2010). Measures of disease frequency: prevalence and incidence. Nephron Clin Pract 115 c17-20. 10.1159/000286345
  * Nsoesie, Mararthe and Brownstein (2013) [author] Nsoesie, ElaineE., Mararthe, MadhavM. and Brownstein, JohnJ. (2013). Forecasting peaks of seasonal influenza epidemics. PLoS Curr 5\. 10.1371/currents.outbreaks.bb1e879a23137022ea79a8c508b030bc
  * Nsoesie et al. (2014) [author] Nsoesie, Elaine OE. O., Brownstein, John SJ. S., Ramakrishnan, NarenN. and Marathe, Madhav VM. V. (2014). A systematic review of studies on forecasting the dynamics of influenza outbreaks. Influenza and other respiratory viruses 8 309–316.
  * Nyabadza, Mukandavire and Hove-Musekwa (2011) [author] Nyabadza, F.F., Mukandavire, Z.Z. and Hove-Musekwa, S. D.S. D. (2011). Modelling the HIV/AIDS epidemic trends in South Africa: Insights from a simple mathematical model. Nonlinear Analysis: Real World Applications 12 2091-2104. https://doi.org/10.1016/j.nonrwa.2010.12.024
  * Osthus et al. (2017) [author] Osthus, DaveD., Hickmann, Kyle SK. S., Caragea, Petruţa CP. C., Higdon, DaveD. and Del Valle, Sara YS. Y. (2017). Forecasting seasonal influenza with a state-space SIR model. Ann Appl Stat 11 202–224. 10.1214/16-AOAS1000
  * Osthus et al. (2019) [author] Osthus, DaveD., Gattiker, JamesJ., Priedhorsky, ReidR. and Valle, Sara Y. DelS. Y. D. (2019). Dynamic Bayesian Influenza Forecasting in the United States with Hierarchical Discrepancy (with Discussion). Bayesian Analysis 14 261 – 312. 10.1214/18-BA1117
  * Piovella (2020) [author] Piovella, NicolaN. (2020). Analytical solution of SEIR model describing the free spread of the COVID-19 pandemic. Chaos, Solitons & Fractals 140 110243\. https://doi.org/10.1016/j.chaos.2020.110243
  * Prangle (2016) [author] Prangle, DennisD. (2016). Lazy ABC. Statistics and Computing 26 171–185. 10.1007/s11222-014-9544-3
  * Schlickeiser and Kröger (2021) [author] Schlickeiser, RR. and Kröger, MM. (2021). Analytical solution of the SIR-model for the temporal evolution of epidemics: part B. Semi-time case. Journal of Physics A: Mathematical and Theoretical 54 175601\. 10.1088/1751-8121/abed66
  * Shaman et al. (2013) [author] Shaman, JeffreyJ., Karspeck, AliciaA., Yang, WanW., Tamerius, JamesJ. and Lipsitch, MarcM. (2013). Real-time influenza forecasts during the 2012–2013 season. Nature communications 4 2837\.
  * Turkyilmazoglu (2021) [author] Turkyilmazoglu, MustafaM. (2021). Explicit formulae for the peak time of an epidemic from the SIR model. Physica D: Nonlinear Phenomena 422 132902\. https://doi.org/10.1016/j.physd.2021.132902
  * Walter and Contreras (1999) [author] Walter, Gilbert G.G. G. and Contreras, MarthaM. (1999). Introduction to Compartmental Models 111–123. Birkhäuser Boston, Boston, MA. 10.1007/978-1-4612-1590-5_13
  * Wang, Wei and Zhang (2014) [author] Wang, XiaoyunX., Wei, LijuanL. and Zhang, JuanJ. (2014). Dynamical analysis and perturbation solution of an SEIR epidemic model. Applied Mathematics and Computation 232 479-486. https://doi.org/10.1016/j.amc.2014.01.090
  * Weiss (2013) [author] Weiss, Howard HowieH. H. (2013). The SIR model and the foundations of public health. Materials matematics 0001–17.
  * Zhan et al. (2019) [author] Zhan, ZhichengZ., Dong, WeihuaW., Lu, YongmeiY., Yang, PengP., Wang, QuanyiQ. and Jia, PengP. (2019). Real-Time Forecasting of Hand-Foot-and-Mouth Disease Outbreaks using the Integrating Compartment Model and Assimilation Filtering. Scientific Reports 9 2661\. 10.1038/s41598-019-38930-y
  * Zhang et al. (2022) [author] Zhang, PeijueP., Feng, KairuiK., Gong, YuqingY., Lee, JieonJ., Lomonaco, SaraS. and Zhao, LiangL. (2022). Usage of Compartmental Models in Predicting COVID-19 Outbreaks. The AAPS Journal 24 98\. 10.1208/s12248-022-00743-9
  * Zhao and Chen (2020) [author] Zhao, ShileiS. and Chen, HuaH. (2020). Modeling the epidemic dynamics and control of COVID-19 outbreak in China. Quantitative Biology 8 11-19. https://doi.org/10.1007/s40484-020-0199-0

##  Appendix A Specification of the Dirichlet-Beta State-Space Model

The DBSSM is fit using a fully Bayesian framework, where the interest lies in
estimating the joint posterior density of the latent space
$\theta_{1:t^{\prime}}$ and $\phi$ given the observed data $y_{1:t^{\prime}}$.
Using the conditional independence assumption described above, this density
can be written as follows:

$\pi(\theta_{1:t^{\prime}},\phi|y_{1:t^{\prime}})~{}\propto~{}\pi(\phi)\pi(y_{1:t^{\prime}},\theta_{1:t^{\prime}}|\phi)=\pi(\phi)\prod_{t=1}^{t^{\prime}}\mathcal{L}(y_{t}|\theta_{t},\phi)\pi(\theta_{t}|\theta_{t-1},\phi),$
(27)

where $\pi(\phi)$ is some prior on $\phi$,
$\mathcal{L}(y_{t}|\theta_{t},\phi)$ is the data likelihood determined by
(17), and the distribution $\pi(\theta_{t}|\theta_{t-1},\phi)$ is determined
by (18). To perform forecasts on observations $y_{t^{\prime}:T}$, where $T$ is
the final timepoint of the outbreak, one uses the posterior predictive
distribution, where the model and latent-space parameters are integrated out:

$\pi(y_{(t^{\prime}+1):T}|y_{1:t^{\prime}})=\int\int\pi(y_{(t^{\prime}+1):T},\theta_{1:T},\phi|y_{1:t^{\prime}})d\theta_{1:T},d\phi.$
(28)

To complete the specification of the DBSSM, one must determine what priors to
put on the model parameters in $\phi$. It is here that the authors in Osthus
et al. (2017) directly address the aforementioned stability issues with
fitting an SIR curve with early pandemic data. Define the latent variable
$z=(PPT,PPV)$, and let $\boldsymbol{\phi}$ be the expanded model parameter
vector that includes $z$. We factorize the new prior distribution on
$\boldsymbol{\phi}$ to get

$\pi(\boldsymbol{\phi})=\pi(\iota)\pi(\lambda|\iota)\pi(\theta_{0}|\lambda,\iota)\pi(z|\theta_{0},\lambda,\iota)\pi(\beta,\gamma|z,\theta_{0},\lambda,\iota).$
(29)

Several modeling assumptions on this conditional distribution give the
abbreviated form,

$\pi(\boldsymbol{\phi})=\pi(\iota)\pi(\lambda)\pi(\theta_{0})\pi(z|\theta_{0})\pi(\beta,\gamma|z,\theta_{0}).$
(30)

Specification of the individual distributions in (30) is what remains to fully
define the DBSSM. The priors $\pi(\iota),\pi(\lambda),$ and $\pi(\theta_{0})$
remain unchanged and can be found in the original paper. In Section 4, the
priors on $\pi(z|\theta_{0})$ and $\pi(\beta,\gamma|z,\theta_{0})$ are
described and updated, when necessary, with the theory developed in this
paper.


11institutetext: Center for Photonic Innovations and Institute for Laser
Science, University of Electro-Communications, Chofu-Shi, Tokyo 182-8585,
Japan. 22institutetext: NS Materials Inc., Tsukushino, Fukuoka 818-0042,
Japan.

# Efficient fiber in-line single photon source based on colloidal single
quantum dots on an optical nanofiber

K. Muhammed Shafi 11 Kali P. Nayak _Corresponding author:_
<EMAIL_ADDRESS>Akiharu Miyanaga 22 Kohzo Hakuta 11

(Received: date / Revised version: date)

###### Abstract

We demonstrate a fiber in-line single photon source based on a hybrid system
of colloidal single quantum dots deposited on an optical nanofiber and cooled
down to cryogenic temperature (3.7 K). We show that a charged state (trion) of
the single quantum dot exhibits a photo-stable emission of single photons with
high quantum efficiency, narrow linewidth (3 meV FWHM) and fast decay time
($10.0\pm 0.5$ ns). The single photons are efficiently coupled to the guided
modes of the nanofiber and eventually to a single mode optical fiber. The
brightness (efficiency) of the single photon source is estimated to be $16\pm
2\%$ with a maximum photon count rate of $1.6\pm 0.2$ MHz and a high single
photon purity ($g^{2}(0)=0.11\pm 0.02$). The device can be easily integrated
to the fiber networks paving the way for potential applications in quantum
networks.

††journal: Applied Physics B

##  1 Introduction

Single photons play a central role in the realization of quantum information
science kimble ; OBrien . Therefore, development of an ideal single photon
source (SPS) is a long-standing frontier of research and currently a
technology on high demand. There have been significant advances and
innovations in this direction, albeit an ideal and on-demand SPS is yet to be
realized Aharonovich ; Senellart .

A straightforward approach to realize an SPS is to collect emission from an
isolated quantum emitter. Neutral atoms or ions may be an orthodox choice for
a quantum emitter. However, isolating single atoms or ions still requires a
complex experimental apparatus. In this context, solid-state quantum emitters
like quantum dots (QDs) or atom-like defect in crystalline hosts Aharonovich ;
Senellart ; Michler  are one of the promising choices for practical SPS, based
on the excellent emission properties and easier techniques to isolate single
emitters. However, the presence of the complex mesoscopic solid-state host
induces various challenges, e.g. inhomogeneous spectral broadening, emission
intermittency and non-identical emitters, etc. In the last decade, there have
been significant efforts and developments to overcome such challenges
Aharonovich ; Senellart ; Park ; Lin .

On the other hand, a key challenge in this approach is to efficiently collect
the single photons emitted from the individual emitter into a single spatial
mode. Moreover, from the view point of applications in quantum networks, the
single photons must be efficiently coupled to a single mode optical fiber
(SMF). Advances in the nanophotonics and nanofabrication have shown promising
developments in this regard. Excellent ideas and technologies have emerged to
manipulate the fluorescence of single emitters using nanophotonic platforms
like nano-waveguides, photonic crystal waveguides or cavities, plasmonic
structures and other microresonators Aharonovich ; Senellart ; BramatiNanorod
; BramatiNanorod2 ; Sipahigil ; Daveau ; Hoang ; Waks . Efficient collection
of single emitter fluorescence into the nanophotonic platforms has been
successfully demonstrated. However, coupling single photons from nanophotonic
platforms to a SMF with high efficiency still remains a challenge Daveau ;
Waks .

In this context, tapered optical fiber with subwavelength diameter waist,
optical nanofiber (ONF), provides a unique fiber in-line platform for
collecting single emitter fluorescence famsan1 ; Kali1 . The distinct point of
the technique is that along with the strong photonic confinement, the guided
mode of the ONF can adiabatically evolve to the SMF mode with near unity
efficiency. This provides an automatic and alignment-free fiber-coupled
platform for an SPS. In the last decade, there has been significant
development to interface the ONF platform with single emitters in the form of
laser-cooled atoms Kali2 ; Kali3 , QDs Yalla2 ; FujiwaraQD , color centers in
nano-diamond ONFNDarno ; Fujiwara , 2D materials Schell2D  and molecules Skoff
.

In particular, using colloidal single QDs on the ONF, it was demonstrated that
the channeling efficiency can reach the theoretical limit of 22% famsan1 ;
Yalla1 , which is determined by the transverse confinement of the photonic
mode in the ONF. Furthermore, it was demonstrated that the channeling
efficiency can be enhanced to 65% using a composite photonic crystal ONF
cavity famsan2 ; Yalla3 . However, the experiments were performed in the room
temperature resulting in a broadband emission from the QD.

The key requirements for an efficient SPS are that the quantum emitter must be
bright, producing high single photon count rate with high quantum efficiency
(QE) and preferably in a narrow spectral window Aharonovich ; Senellart .
Various quantum photonics protocols require narrow band (in ideal case,
indistinguishable) single photons Aharonovich ; Senellart . In this context,
cryogenic operation of the ONF/QD hybrid system may suppress the phonon
mediated broadening of the QD emission, leading to a narrow band SPS.

We have reported the cryogenic cooling of the ONF/QD hybrid system to a
temperature of 3.7 K Shafi1 ; Shafi2 . A gradient thick shell type CdSe QDs
with an outermost ZnS shell, was used for the experiments. The spectroscopic
studies at the cryogenic temperature had revealed that the photo-luminescence
(PL) emission occurs not only from the neutral exciton but also from the
charged exciton known as trion Shafi1 ; Shafi2 ; Louyer . Moreover, from the
PL decay profile, it was inferred that the upper state of the neutral exciton,
consists of two levels with a dark state. This leads to an effectively slow
decay component for the neutral exciton Shafi1 ; Biadala . While, the trion
state behaves effectively as a two-level system and the decay rate is one
order faster than the effective slow decay component of the neutral exciton
Shafi1 ; Fernee1 ; Javaux . Therefore, the presence of the dark state of the
neutral exciton hinders the implementation of the system as an efficient SPS.

Recently, we have found that at cryogenic temperature, the CdSe QDs can be
photo-charged by irradiating with shorter wavelength lasers. In particular,
the QDs can be prepared in a permanently charged state using a 355 nm laser
which eradicates the issues related to the dark state of the neutral exciton.
A systematic investigation of the photo-charging behavior of the QDs measured
by varying the wavelength and the intensity of the irradiating laser will be
reported elsewhere.

Here, we demonstrate the capability of the ONF/QD hybrid system as an
efficient fiber in-line SPS. We show that at cryogenic temperature (3.7 K),
the charged single QD exhibits a photo-stable emission of single photons with
high QE, narrow linewidth (3 meV FWHM) and fast decay time ($10.0\pm 0.5$ ns).
The single photons are efficiently coupled to the guided modes of the ONF and
eventually to a SMF. The brightness (efficiency) of the SPS is estimated to be
16$\pm 2\%$ with a maximum photon count rate of 1.6$\pm 0.2$ MHz and a high
single photon purity ($g^{2}(0)=0.11\pm 0.02$).

Figure 1: Schematic diagram of the fiber in-line single photon source. The
ONF/QD hybrid system is cooled down in a cryostat. The single QD is irradiated
with a pulsed laser launched through the optical window of the cryostat. The
emission from the single QD is coupled to the guided mode of the ONF and
eventually to a single mode optical fiber (SMF). The temporal and spectral
characteristics of the fiber-guided photons are measured simultaneously at the
opposite ends of the fiber using a time correlated single photon counter
(TCSPC) and an optical multi-channel spectrum analyzer (OMA), respectively.
The abbreviations LPF, IF, BS and APD denote long-pass filter, interference
filter, beam splitter and avalanche photodiode (single photon counting
module), respectively.

##  2 Experimental setup

Figure 1 shows the schematic diagram of the experiment. The experimental setup
is based on a hybrid system of single QDs deposited on an ONF and cooled down
inside a cryostat. The QDs used for the experiment were colloidal gradient
thick shell type CdSe QDs with an outermost ZnS shell. The QDs were
synthesized using a method described in Ref. Shafi1 . At room temperature, the
QDs emit at a wavelength of 640 nm and the emission QE was measured to be
85$\pm$5% Shafi1 . ONFs with optical transmission of $>$99% were fabricated by
adiabatically tapering commercial SMFs (SM 600, Fibercore) using a heat and
pull technique. ONFs used for the experiment had a waist diameter of
310$\pm$10 nm and a uniform waist length of 2.5 mm. The ONF diameter was
chosen to maximize the channeling of the single QD emission into the ONF
guided modes with an efficiency as high as 22% famsan1 ; Yalla1 .

The QDs were deposited on the ONF using a computer-controlled subpico-liter-
needle-dispenser system combined with an inverted microscope Yalla2 ; Shafi1 .
The dispenser consists of a tapered glass-tube containing the QD solution and
a tungsten needle with a tip diameter of 5 $\mu$m. Once the needle tip passes
through the glass-tube, it carries a small amount of QD solution at its edge
and a tiny amount of this solution is deposited on the ONF surface by bringing
the tip into contact with the ONF. To achieve reproducible deposition of
single QD, the concentration of the QDs was adjusted by diluting the colloidal
solution of the QDs in toluene. From various trials, we have estimated an
overall success probability of depositing single QD for each trial to be about
60%.

For the ONF sample presented in this work, ten depositions were performed with
a spatial period of 150 $\mu$m. The transmission of the ONF was $>$98% after
the deposition procedure. After the depositions, the ONF/QD hybrid system was
installed into a custom-designed optical cryostat Shafi1 . The ONF
transmission was reduced to 90% during the installation process. The ONF/QD
hybrid system was cooled down to 3.7 K via He buffer gas cooling Shafi1 . The
transmission was maintained throughout the cooling process and at the
cryogenic temperatures.

The emission characteristics of the QD were measured from the PL photons
coupled into the guided mode of the ONF and eventually to a SMF. For measuring
the emission characteristics, the QD was excited with a picosecond pulsed-
laser with a pulse width of 20 ps FWHM and a wavelength of 532 nm. For photo-
charging, the QD was irradiated with a cw laser at a wavelength of 355 nm. The
excitation lasers were launched through the optical window of the cryostat and
irradiated perpendicular to the ONF. The excitation lasers were focused on the
ONF to a spot-size of 19 $\mu$m FWHM using a combination of lenses placed
outside the cryostat. The polarization of the excitation lasers was linear
with a polarization axis perpendicular to the ONF axis.

The spectral and temporal characteristics of the fiber-guided photons were
measured simultaneously at the opposite ends of the fiber. At one end of the
fiber, the PL spectra of the QD was measured using an optical multi-channel
analyser (OMA) with a resolution of 1.0 meV (0.3 nm). The scattering from the
excitation laser coupled into the fiber guided mode was filtered out using a
560 nm long-pass filter (O56, HOYA). At the other end of the fiber, the
temporal characteristics were measured using single photon counting modules
(APDs, Perkin Elmers) and a time correlated single photon counting system
(TCSPC, PicoHarp 300). The single photon characteristics of the emitted
photons were inferred from the photon correlations measured using a Hanbury-
Brown-Twiss (HBT) setup with two APDs and the TCSPC. As illustrated in Fig. 1,
the fiber output was first filtered using a band-pass interference filter (2
nm or 15 nm FWHM) to remove the background photons from the signal, then split
into two channels using a beam splitter. The photons in the two channels were
measured using two separate APDs and the arrival times of the photons were
recorded using a TCSPC. The photon correlations were extracted from the
coincidence counts between the two channels for varying delay time. The PL
decay profile was measured by observing the temporal correlations between the
excitation laser pulse and the PL photons. The repetition rate of the pulsed
laser was set to 10 MHz for the measurements of PL spectrum, photon count rate
and photon correlations. While, only for the measurement of PL decay profiles,
the repetition rate was set to 500 kHz. All data presented here were measured
at a temperature of 3.7 K.

##  3 Results

The blue trace in Fig. 2(a) shows the PL spectrum of the single QD measured at
3.7 K with an integration time of 120 s. The spectrum displays well-resolved
four peaks labeled as 1, 2, 3 and 4. The observed linewidth (FWHM) of the peak
1 is 3.2 meV and that of the peak 2 is 1.9 meV and they are separated by 16.5
meV. This leads to the assignment of the peaks 1 and 2, as the zero phonon
lines (ZPLs) for the exciton state of a neutral CdSe nanocrystal and for the
trion state of a charged CdSe nanocrystal, respectively Shafi1 ; Louyer ;
Biadala . The linewidth of the peak 1 is broader than that of the peak 2,
since the upper state of the exciton consists of a bright and a dark exciton
state separated by 1.5 meV which was not resolved. The peaks 3 and 4 are red-
shifted by 26.5 meV from the peaks 1 and 2, respectively and are attributed to
the longitudinal optical phonon (LO-phonon) replica of the exciton and the
trion ZPLs Shafi1 ; Shafi2 ; Louyer . Note that the emission from the exciton
and the trion state occurs almost randomly and the occurrence probability
reaches a stationary value for the integration time as long as 120 s Shafi1 .
We estimated the trion occurrence probability (TOP) to be 44$\pm 3$% from the
ratio between the photon counts emitted by the trion state and the total
photon counts from both the trion and the exciton states.

Figure 2: PL characteristics of a single QD in the ONF/QD hybrid system at 3.7
K. (a) PL spectrum, (b) PL decay profile and (c) photon count rate for a
single QD measured before photo-charging. The blue traces are the
experimentally measured data. The yellow curve in (b) shows the fit using a
sum of three exponential functions. (d) PL spectrum, (e) PL decay profile and
(f) photon count rate for a single QD measured after photo-charging. The red
traces are the experimentally measured data. The green curve in (e) shows the
fit using a single exponential function. The vertical axes in (b) and (e) are
in log-scale. All data were measured for an excitation fluence of 6.4
$\mu$J/cm2.

The blue trace in Fig. 2(b) shows the PL decay profile of the single QD at 3.7
K. The decay curve was acquired for an integration time of 200 s with a
resolution of 64 ps. The observed decay curve can be well fitted with a sum of
three exponential decay functions corresponding to the decay processes for the
trion state and the effective fast and slow decay components of the neutral
exciton state. The yellow curve in Fig. 2(b) shows the fit. From the fit, we
extract three decay times of 2.3$\pm$0.3 ns, 9.6$\pm$0.3 ns and 120$\pm$5 ns
corresponding to the effective fast decay component of the neutral exciton
state, the decay of the trion state and the effective slow decay component of
the neutral exciton state, respectively Shafi1 ; Biadala ; Labeau .

Figure 2(c) shows the photon count rate of the PL of the single QD observed at
3.7 K. The resolution of the time axis is 17 ms. It may be seen that the
photon count rate is stable without any intermittency. The observed count rate
is 13.0$\pm$0.5 kHz. It should be noted that at room temperature the PL
emission shows intermittency due to the lower QE of the trion Shafi1 ; Javaux
. The absence of intermittency at cryogenic temperature indicates that both
the exciton and trion state have the same QE Shafi1 ; Javaux . As detailed in
Ref. Shafi1 , the typical overall QE for such QDs is higher than 94%.

Figure 2(d) shows the PL spectrum of the single QD after the QD was irradiated
with 355 nm laser at 5 W/cm2 for about 10 min. The observed spectrum shows
only a single emission line corresponding to the trion energy of 1.99 eV (623
nm) along with a LO phonon replica. The absence of the neutral exciton peak
suggests that the emission occurs only from the trion state with 100% TOP.
Therefore, it is inferred that the irradiation with 355 nm resulted in
permanent photo-charging of the QD.

To further confirm the charged situation, we have measured the PL decay
profile of the photo-charged QD. The red trace in Fig. 2(e) displays the PL
decay profile of the trion emission. The observed PL decay profile shows a
single exponential decay which is consistent with the single upper state model
of the trion state Javaux . The decay profile is fitted by a single
exponential function (green curve) yielding a decay time of 10.0$\pm$0.5 ns.
The decay time is one order shorter than the effective slow decay component of
the neutral exciton state.

Figure 2(f) displays the typical photon count rate for the photo-charged QD.
It may be seen that there is no significant change in the photon count rate
compared to Fig. 2(c). This suggests that there is no change in the overall QE
of the QD Shafi1 ; Fernee1 ; Javaux .

The purity of the fiber in-line SPS was measured by carrying out photon
correlation measurements of the fiber-guided PL photons from the single QD. To
ensure complete relaxation of the excitation between the subsequent laser
pulses, the excitation repetition rate was set to be slower than the PL decay
rate. The red trace in Fig. 3(a) shows the photon correlations measured for an
excitation fluence of 6.4 $\mu$J/cm2 at 10 MHz repetition rate. The
coincidences between the outputs of the beam splitter are plotted against the
delay time ($\tau$) with a time-bin of 1 ns. It may be seen that periodic
peaks with an interval of 100 ns are observed, well corresponding to the
repetition rate. However, the central peak at $\tau=0$ is strongly suppressed
clearly indicating the single photon characteristics of the source. We
estimate the second-order correlation value at $\tau=0$ to be g2(0)= 0.06$\pm
0.01$. The value of g2(0) was estimated from the ratio between the area of the
peak at $\tau=0$ and the averaged area of side peaks on each side at $\pm 100$
ns. The value of g2(0) gives a measure of the purity of the fiber in-line SPS.

Figure 3: Photon statistics of the PL of the photo-charged single QD measured
at 3.7 K. (a) Photon correlation measurement of the fiber-guided PL photons
for an excitation fluence of 6.4 $\mu$J/cm2. (b) Photon count rates (black
squares with error bars) at one-end of the fiber measured for different
excitation fluences. The red curve is the fit using the saturation model given
in Eq. (1). (c) Photon correlation measurement for an excitation fluence of
108 $\mu$J/cm2 which is larger than the ${F}_{sat}$. (d) g2(0) for different
excitation fluences.

To estimate the brightness (efficiency) of the SPS, the excitation fluence was
increased until the photon count rate shows saturation behavior. Figure 3(b)
shows the observed photon count rates (black squares with error bars) at one
side of the fiber for different excitation fluences. It may be seen that for
higher excitation fluence the photon count rate shows a saturation behavior.
The observations are fitted using a saturation model for two-level system
Daveau ; Waks  given as

${I(F)}={I}_{max}[1-\exp(-\frac{F}{F_{sat}})]$ (1)

where ${I}_{max}$ is the saturated photon count rate, ${F}$ and ${F}_{sat}$
are the excitation fluence and the saturation fluence, respectively. The red
curve in Fig. 3(b) shows the fit. From the fit, we obtained a maximum count
rate of single photons at one side of the fiber to be ${I}_{max}$= $327\pm 22$
kHz and a saturation fluence of ${F}_{sat}=79\pm 6$ $\mu$J/cm2.

To confirm the purity of the single photons at the saturation condition, we
have measured the photon correlations at higher excitation fluence. Figure
3(c) shows the photon correlations measured at a fluence of 108 $\mu$J/cm2
which is well above the ${F}_{sat}$. It may be seen that the central peak at
$\tau=0$ is still strongly suppressed indicating that the single photon
characteristics of the SPS is maintained at saturation condition. From the
photon correlations, we estimated a g2(0) = $0.11\pm 0.02$ at this excitation
fluence. Figure 3(d) displays the values of g2(0) estimated for different
excitation fluences. It may be seen that g2(0) gradually increases from 0.06
to 0.11 at the saturation condition. This indicates that although the multi-
photon emission probability is not significant for the present QD, but at
higher excitation fluence it is not negligible.

We have also recorded the spectral and photon count rate characteristics of
the photo-charged single QD at different excitation conditions. Figure 4(a)
displays the linewidth (red squares with error bars) and the peak energy
(black squares with error bars) of the trion PL spectra measured at different
excitation fluences. The linewidths are obtained after correcting for the
instrumental response function with a width of 1 meV. It may be seen that
there is a systematic increase in the linewidth with increasing excitation
fluence. The linewidth increased from 1.1 meV (at 6.4 $\mu$J/cm2) to 3.2 meV
(at 108 $\mu$J/cm2). On the other hand, the peak energy shift is negligible
(within the resolution limit). Figure 4(b) shows the measured photon count
rate at a fluence of 108 $\mu$J/cm2. It may be seen that the fluctuations in
the photon count rate are larger compared to that shown in Fig. 2(f). The
average photon count rate is estimated to be $280\pm 30$ kHz.

Figure 4: PL characteristics of the photo-charged single QD at higher
excitation fluence. (a) The linewidth (red squares with error bars) and the
peak energy (black squares with error bars) of the trion PL spectra measured
at different excitation fluences. The linewidths are obtained after correcting
for the instrumental response function with an width of 1 meV. (b) Photon
count rate measured at an excitation fluence of 108 $\mu$J/cm2 which is larger
than the ${F}_{sat}$.

##  4 Discussion

The total detection efficiency (${\alpha}$) of the single photons from the ONF
to the APD was estimated to be 40% by taking into account various efficiencies
of the optical system which include the ONF transmission (90%), the filter
transmission efficiency (83%), the fiber-coupling efficiency (90.5%) to the
APD and the detection efficiency of the APD at 620 nm (60%). By correcting for
the transmission and detection efficiencies, we estimate a maximum photon
count rate (${\Gamma_{SP}}$ = $\frac{2{I}_{max}}{\alpha}$) of 1.6$\pm$0.2 MHz
for the fiber in-line SPS. Here the factor two accounts for the total
channeling into both ends of the fiber. The brightness (source efficiency,
$\beta$) of the fiber in-line SPS is estimated by relating ${\Gamma_{SP}}$ to
the repetition rate ($R$) of the excitation laser as $\beta$ =
$\frac{\Gamma_{SP}}{R}$ = 16$\pm$2%. It should be noted that $\beta$ can be
related to the QE ($\eta_{q}$) of the QD and the channeling efficiency
($\eta$) into the ONF guided modes given as $\beta=\eta_{q}\eta$. It is
reasonable to assume a maximum value for the channeling efficiency
($\eta=22\%$) famsan1 ; Yalla1  which would yield a higher $\beta$-value. The
reduction in the observed $\beta$ may indicate an overall reduction in the QE
of the QD at higher excitation fluence.

As shown in Fig. 4(b), at high excitation fluence, the photon count rate shows
larger fluctuations with a characteristic intermittency behavior. The
intermittency may indicate the presence of another emission channel with lower
QE. Therefore, the intermittency can effectively reduce the QE of the single
QD at high excitation fluence. The intermittency might be related to the
emission from Auger-assisted multi-exciton process at high excitation fluence
Javaux ; Canneson . This is also evident from the increase in the $g^{2}(0)$
value at high excitation fluence. The increase in the linewidth of the PL
spectra and onset of the intermittency may also indicate local heating at high
excitation fluence Shafi1 ; Javaux . However, the negligible shift in the peak
energy suggests that the local temperature was well maintained Shafi1  and the
effect of temperature change may not explain the broadening of the spectral
width. The observed increase in the spectral width may be attributed to the
enhancement of the spectral diffusion process due to the onset of multi-
exciton process Empedocles .

Another issue to address at higher excitation fluence is the background
induced by the PL of the ONF itself. The level of such background depends on
the spot size of the excitation laser and choice of appropriate band-pass
filter. In the present cryostat setup, the spot size is limited to 19 $\mu$m
and at high excitation fluence, the PL of the ONF induces background photons.
However using a 2 nm band-pass filter we can reasonably suppress the
background. This is evident from the photon coincidences plotted in Fig. 3(c).
We have found that under similar conditions, using a 15 nm band-pass filter
will induce background photons which adds a constant offset to the photon
coincidences. This offset is $<10$% of the coincidence peak height. In
principle using an objective lens inside the cryostat can reduce the spot-size
to $\sim$ 1 $\mu$m to further suppress the background level. It should be
noted that it is one of the advantages of using the ONF platform compared to
emitters on the bulk substrates where the excitation volume is much larger
resulting in significant background from the PL of the substrate itself.

Out of the ten deposited positions on the ONF, we have found single QDs at
five positions with similar emission characteristics. The data presented here
are representative data measured at one of the positions. For low excitation
fluence, we have observed no significant degradation of the emission from the
single QDs at least up to 10 days of operation at cryogenic temperatures.
However, for excitation fluence higher than the saturation fluence, we have
observed degradation of the QD emission within 2 h of continuous operation. To
realize a practical device for SPS, further engineering of the QDs will be
essential for robust operating conditions.

It should be noted that the QD emission into the ONF guided modes, should
depend on the excitation polarization. However, the information about the
dipole orientation and azimuthal position of the QD will play a significant
role, which is currently not fully explored. At different deposited positions
on the ONF, we have observed a variation in the visibility ranging from 10 to
60%, measured by changing the excitation polarization from perpendicular to
parallel with respect to the ONF axis. The polarization axis was chosen
perpendicular to the ONF axis to maximize the photon count rate into the ONF
guided modes.

To enhance the performance of the hybrid ONF/QD system as an SPS, an important
step would be to incorporate a cavity structure on the ONF and establish an
ONF cavity QED system at cryogenic temperatures. Although Purcell enhancement
of solid-state quantum emitters on ONF cavities has been successfully
demonstrated Yalla3 ; Schell , the experiments were performed at room
temperature. As a result, the emission spectral widths were 20-40 nm broad and
the cavity enhancement was limited only to a narrow region of the emission
spectrum. This limits the overall efficiency of the system as an SPS. In the
present hybrid system at cryogenic temperatures, the spectral diffusion width
saturates at 3 meV ($\sim$1.0 nm) FWHM, which is comparable to the typical
cavity linewidths reported for the ONF cavity QED experiments Yalla3 ; Schell
; SileFIB ; TakeuchiFIB . As a result, the entire emission band can be
strongly enhanced into the ONF cavity mode leading to higher source
efficiency. Due to the faster cavity enhanced decay rate, one may also expect
a higher photon count rate. Moreover using a one-sided cavity, all photons can
be extracted at one end of the fiber. Another key aspect of an SPS is the
indistinguishability of the single photons. This will require further
narrowing of the spectral width by employing near resonant excitation scheme
Fernee2  and much narrow linewidth ONF cavities as in Ref. Kali4 .

##  5 Conclusions

In conclusion, we have demonstrated a fiber in-line single photon source based
on a hybrid system of colloidal single quantum dots deposited on an optical
nanofiber and cooled down to cryogenic temperature (3.7 K). We have shown that
the charged single quantum dot exhibits a photo-stable emission of single
photons with high quantum efficiency, narrow linewidth (3 meV FWHM) and fast
decay time ($10.0\pm 0.5$ ns). The single photons were efficiently coupled to
the guided modes of the nanofiber and eventually to a single mode optical
fiber. The brightness (efficiency) of the single photon source was estimated
to be 16$\pm 2\%$ with a maximum photon count rate of 1.6$\pm 0.2$ MHz and a
high single photon purity ($g^{2}(0)=0.11\pm 0.02$). The device can be easily
integrated to the fiber networks paving the way for potential applications in
quantum networks.

##  6 Acknowledgement

This work was supported by the Japan Science and Technology Agency (JST)
through Strategic Innovation Program (Grant No. JPMJSV0918). We acknowledge
the contributions of Kazunori Iida and Emi Tsutsumi, in the preparation of the
thick shell quantum dot samples.

## References

  * (1) H. J. Kimble, The quantum internet, Nature 453, 1023 (2008).
  * (2) J. O’Brien, A. Furusawa, and J. Vuckovic, Photonic quantum technologies, Nat. Photon 3, 687 (2009).
  * (3) I. Aharonovich, D. Englund, and M. Toth, Solid-state single-photon emitters, Nat. Photonics 10, 631 (2016).
  * (4) P. Senellart, G. Solomon, and A. White, High-performance semiconductor quantum-dot single-photon sources, Nat. Nanotechnol. 12, 1026 (2017).
  * (5) P. Michler, A. Kiraz, C. Becher, W. Schoenfeld, P. Petroff, L. Zhang, E. Hu, and A. Imamoglu, A quantum dot single-photon turnstile device, Science 290, 2282 (2000).
  * (6) Y.-S. Park, J. Lim, and V. I. Klimov, Asymmetrically strained quantum dots with non-fluctuating single-dot emission spectra and subthermal room-temperature linewidths, Nat. Mater. 18, 249 (2019).
  * (7) X. Lin, X. Dai, C. Pu, Y. Deng, Y. Niu, L. Tong, W. Fang, Y. Jin, and X. Peng, Electrically-driven single-photon sources based on colloidal quantum dots with near-optimal antibunching at room temperature, Nat. Commun. 8, 1132 (2017).
  * (8) M. Manceau, S. Vezzoli, Q. Glorieux, F. Pisanello, E. Giacobino, L. Carbone, M. De Vittorio and A. Bramati, Effect of charging on CdSe/CdS dots-in-rods single-photon emission, Phys. Rev. B 90, 035311 (2014).
  * (9) S. Vezzoli, M. Manceau, G. Lemenager, Q. Glorieux, E. Giacobino, L. Carbone, M. De Vittorio and A. Bramati, Exciton fine straucture of CdSe/CdS nanocrystals determined by polarization microscopy at room temperature, ACS Nano 9, 7992 (2015).
  * (10) A. Sipahigil, R. E. Evans, D. D. Sukachev, M. J. Burek, J. Borregaard, M. K. Bhaskar, C. T. Nguyen, J. L. Pacheco, H. A. Atikian, C. Meuwly, R. M. Camacho, F. Jelezko, E. Bielejec, H. Park, M. Loncar, and M. D. Lukin, An integrated diamond nanophotonics platform for quantum-optical networks, Science 354, 847 (2016).
  * (11) R. S. Daveau, K. C. Balram, T. Pregnolato, J. Liu, E. H. Lee, J. D. Song, V. Verma, R. Mirin, S. W. Nam, L. Midolo, S. Stobbe, K. Srinivasan, and P. Lodahl, Efficient fiber-coupled single-photon source based on quantum dots in a photonic-crystal waveguide, Optica 4, 178 (2017).
  * (12) T. B. Hoang, G. M. Akselrod, and M. H. Mikkelsen, Ultrafast room-temperature single photon emission from quantum dots coupled to plasmonic nanocavities, Nano Lett. 16, 270 (2016).
  * (13) C.-M. Lee, M. A. Buyukkaya, S. Aghaeimeibodi, A. Karasahin, C. J. K. Richardson, and E. Waks, A fiber-integrated nanobeam single photon source emitting at telecom wavelengths, Appl. Phys. Lett. 114,, 171101 (2019).
  * (14) F. L. Kien, S. D. Gupta, V. I. Balykin, and K. Hakuta, Spontaneous emission of a cesium atom near a nanofiber: efficient coupling of light to guided modes, Phys. Rev. A 72, 032509 (2005).
  * (15) K. P. Nayak, M. Sadgrove, R. Yalla, F. L. Kien and K. Hakuta, Nanofiber quantum photonics, J. Opt. 20, 073001 (2018).
  * (16) K. P. Nayak, P. N. Melentiev, M. Morinaga, F. L. Kien, V. I. Balykin, and K. Hakuta, Optical nanofiber as an efficient tool for manipulating and probing atomic fluorescence, Opt. Express 15, 5431 (2007).
  * (17) K. P. Nayak, J. Wang, and J. Keloth, Real-time observation of single atoms trapped and interfaced to a nanofiber cavity, Phys. Rev. Lett. 123, 213602 (2019).
  * (18) R. Yalla, K. P. Nayak, and K. Hakuta, Fluorescence photon measurements from single quantum dots on an optical nanofiber, Opt. Express 20, 2932 (2012).
  * (19) M. Fujiwara, K. Toubaru, T. Noda, H.-Q. Zhao, and S. Takeuchi, Highly efficient coupling of photons from nanoemitters into single-mode optical fibers, Nano Lett. 11, 4362 (2011).
  * (20) L. Liebermeister, F. Petersen, A. V. Münchow, D. Burchardt, J. Hermelbracht, T. Tashima, A. W. Schell, O. Benson, T. Meinhardt, A. Krueger, A. Stiebeiner, A. Rauschenbeutel, H. Weinfurter, and M. Weber, Tapered fiber coupling of single photons emitted by a deterministically positioned single nitrogen vacancy center, Appl. Phys. Lett. 104, 031101 (2014).
  * (21) M. Fujiwara, H.-Q. Zhao, T. Noda, K. Ikeda, H. Sumiya, and S. Takeuchi, Ultrathin fiber-taper coupling with nitrogen vacancy centers in nanodiamonds at cryogenic temperatures, Opt. Lett. 40, 5702 (2015).
  * (22) A. W. Schell, H. Takashima, T. T. Tran, I. Aharonovich, and S. Takeuchi, Coupling quantum emitters in 2D materials with tapered fibers, ACS Photonics 4, 761–767 (2017).
  * (23) S. M. Skoff, D. Papencordt, H. Schauffert, B. C. Bayer, and A. Rauschenbeutel, Optical-nanofiber-based interface for single molecules, Phys. Rev. A 97, 043839 (2018).
  * (24) R. Yalla, F. L.Kien, M. Morinaga, and K. Hakuta, Efficient channeling of fluorescence photons from single quantum dots into guided modes of optical nanofiber, Phys. Rev. Lett. 109, 063602 (2012).
  * (25) F. L. Kien and K. Hakuta, Cavity-enhanced channeling of emission from an atom into a nanofiber, Phys. Rev. A 80, 053826 (2009).
  * (26) R. Yalla, M. Sadgrove, K. P. Nayak, and K. Hakuta, Cavity quantum electrodynamics on a nanofiber using a composite photonic crystal cavity, Phys. Rev. Lett. 113, 143601 (2014).
  * (27) K. M. Shafi, W. Luo, R. Yalla, K. Iida, E. Tsutsumi A. Miyanaga and K. Hakuta, Hybrid system of an optical nanofibre and a single quantum dot operated at cryogenic temperatures, Sci. Rep. 8, 13494 (2018).
  * (28) K. M. Shafi, W. Luo, R. Yalla and K. Hakuta, Emission characteristics for a single CdSe quantum dot on an optical nanofiber at cryogenic temperatures, J. Phys. Conf. Ser. 1220, 012025 (2019).
  * (29) Y. Louyer, L. Biadala, P. Tamarat, and B. Lounis, Spectroscopy of neutral and charged exciton states in single CdSe/ZnS nanocrystals, Appl. Phys. Lett. 96, 203111 (2010).
  * (30) L. Biadala, Y. Louyer, P. Tamarat and B. Lounis, Direct observation of the two lowest exciton zero-phonon lines in single CdSe/ZnS nanocrystals, Phys. Rev. Lett. 103, 037404 (2009).
  * (31) M. J. Fernee, C. Sinito, Y. Louyer, C. Potzner, T. Nguyen, P Mulvaney, P. Tamarat and B. Lounis, Magneto-optical properties of trions in non-blinking charged nanocrystals reveal an acoustic phonon bottleneck, Nat. Commun. 3, 1287 (2012).
  * (32) C. Javaux, B. Mahler, B. Dubertret, A. Shabaev, A. V. Rodina, Al. L. Efros, D. R. Yakovlev, F. Liu, M. Bayer, G. Camps, L. Biadala, S. Buil, X. Quelin and J-P. Hermier, Thermal activation of non-radiative Auger recombination in charged colloidal nanocrystals, Nature Nanotech. 8, 206 (2013).
  * (33) O. Labeau,P. Tamarat, and B. Lounis, Temperature dependence of the luminescence lifetime of single CdSe/ZnS quantum dots, Phys. Rev. Lett. 90, 257404 (2003).
  * (34) D. Canneson, L. Biadala, S. Buil, X. Quelin, C. Javaux, B. Dubertret, and J.-P. Hermier, Blinking suppression and biexcitonic emission in thick-shell CdSe/CdS nanocrystals at cryogenic temperature, Phys. Rev. B 89, 035303 (2014).
  * (35) S. A. Empedocles and M. G. Bawendi, Influence of spectral diffusion on the line shapes of single CdSe nanocrystallite quantum dots, J. Phys. Chem. B 103, 1826 (1999).
  * (36) A. W. Schell, H. Takashima, S. Kamioka, Y. Oe, M. Fujiwara, O. Benson and S. Takeuchi, Highly efficient coupling of nanolight emitters to a ultra-wide tunable nanofibre cavity, Sci. Rep. 5, 9619 (2016).
  * (37) W. Li, J. Du, V. G. Truong, and S. Nic Chormaic, Optical nanofiber-based cavity induced by periodic air-nanohole arrays, Appl. Phys. Lett. 110, 253102 (2017).
  * (38) H. Takashima, A. Fukuda, H. Maruya, T. Tashima, A. W. Schell, and S. Takeuchi, Fabrication of a nanofiber bragg cavity with high quality factor using a focused helium ion beam, Opt. Express 27, 6792 (2019).
  * (39) M. J. Fernee, C. Sinito, Y. Louyer, P. Tamarat and B. Lounis, The ultimate limit to the emission linewidth of single nanocrystals, Nanotechnology 24, 465703 (2013).
  * (40) K. P. Nayak, P. Zhang, and K. Hakuta, Optical nanofiber-based photonic crystal cavity, Opt. Lett. 39, 232 (2014).