# Finding Closure Terms Directly from Coarse Data for 2D Turbulent Flow

Xianyang Chen  Jiacai Lu  Grétar Tryggvason  Department of Mechanical
Engineering
Johns Hopkins University, MD, USA

###### Abstract

Machine learning is used to develop closure terms for coarse grained model of
two-dimensional turbulent flow directly from the coarse grained data by adding
a source term to the Navier-Stokes equations to ensure that the coarse-grained
flow evolves in the correct way. The source term is related to the average
flow using a Neural Network with a relatively simple structure and smoothed
slightly to prevent instabilities in a posteriori test. The time dependent
coarse grained flow field is generated by filtering fully resolved results and
the predicted coarse field evolution agrees well with the filtered results,
both for the flow used to learn the closure terms and for flows not used for
the learning.

##  I Introduction

Direct Numerical Simulations (DNS) where every spatial and temporal scales in
turbulent flows are fully resolved usually require fine grids and long
computational times and are therefore too computationally intensive for most
flows of industrial interest. Thus, there is a need for reduced order
representations of the flow such as Large Eddy Simulation (LES), which are
able to predict some of the statistics by evolving only the large scales,
using subgrid models to account for the unresolved motion. Generating
appropriate subgrid models is still an active area of research and often
models perform well for one specific flow but not for others. The classical
Smagorinsky model, for example, is too dissipative in wall-bounded flows and
prevents laminar flows from transitioning to turbulence ([1]). DNS that fully
resolve all spatial and temporal scales for moderate Reynolds number flows
provides enormous amount of data which has been widely used for the
development of subgrid models, usually by adjusting coefficients in
relationships between the subgrid stresses and the coarse flow, proposed by
scaling and other considerations. The recent availability of methods to
“automatically” extract correlations from very large data-sets are offering
new ways to relate subgrid stresses to the resolved coarse flow. While current
focus is mostly on single phase flow, this should be particularly useful for
more complex flows, where the details of the subgrid scales are less well
known.

In LES the Navier-Stokes equations are filtered, which results in an
additional term $\mathcal{F}=\nabla\cdot{\boldsymbol{\tau}}_{SGS}$, where the
subgrid stresses are related to the filtered fully resolved flow by
${\boldsymbol{\tau}}_{SGS}=\widetilde{{\bf u}{\bf u}}-\tilde{\bf u}\tilde{\bf
u}$. Here, $\tilde{\bf u}$ is the filtered velocity and $\widetilde{{\bf
u}{\bf u}}$ is the filtered product. Comparisons of
${\boldsymbol{\tau}}_{SGS}$ computed from a model and computed from a filtered
fully resolved field are generally referred to as a priori tests. Here,
instead of using the filtered fully resolved solution, we compute
$\mathcal{F}$ so that the time dependent coarse flow evolves correctly. The
coarse flow is found by filtering a fully resolved one, but otherwise it is
not used (except for comparisons) and in principle the coarse flow could come
from other sources, such as experimental data, or from a field coarsened in
such a way that the connections with the filtered governing equations is not
obvious.

A number of researchers have recently explored the use of machine learning to
develop reduced order models in fluid mechanics ([2],[3], [4], [5], [6], [7],
[8], [9]). See [10] for a review. Most of the contributions have focused on
closure models for Reynolds-Averaged Navier-Stokes (RANS) equations, but some
efforts have also been made for LES. [11] used a Neural Network (NN) to
predict the turbulent viscosity coefficient for the mixed model in a turbulent
channel flow. [12] modeled the subgrid stress tensor directly using NN for
turbulent channel flows, without making any assumption about the model form,
but the results did not show obvious advantages over the Smagorinsky model.
[13] reported instabilities in a posteriori test with closure models using
Residual Neural Network (RNN), even though a priori test showed a high
correlation coefficient. Predicting the turbulent viscosity instead of the
subgrid stresses did, however, eliminate the instabilities. [14] included
second derivatives of the velocity field as inputs for NN to predicted subgrid
stresses, showing a better result than the classical Smagorinsky and the
dynamic Smagorinsky model for the energy dissipation. [15] used NN to model
the source term, rather than the subgrid stresses in homogeneous isotropic
turbulence and while they found a high correlation coefficients (0.99), they
had to add artificial diffusion to damp out instabilities. They also extended
their modeling to turbulence in compressible flows in [16]. We note that while
[15] work with the source term, they compute it from the subgrid stresses,
rather than finding it directly from the coarse field, as we do. For 2D
turbulence, [17] presented a priori test using NN models for deconvolution of
the flow field, as well as for 3D homogeneous isotropic turbulence, and
turbulence in stratified flow. In [18], they extended their work by conducting
a posteriori test for 2D turbulence, but had to truncate the predicted subgrid
source terms to supress instabilities.

##  II Method

###  II.1 The coarse field

We conduct fully resolved simulations of unsteady two dimensional flow and
filter the results to produce a coarse grained velocity and pressure fields.
The simulations are done using a standard projection method on a regular
staggered grid. A second order Runge-Kutta method is used for the time
integration, a second order Quick scheme is used for the advection terms and
second order centered approximations are used for all other spatial
derivatives. The coarse grained field, denoted by $\mathbf{\tilde{u}}$, is
generated by spatial filtering. Here we use a Gaussian:

$G_{\Delta}=(\frac{6}{\pi\Delta^{2}})^{3/2}exp(\frac{-6\|\mathbf{x}-\mathbf{x^{\prime}}\|^{2}}{\Delta^{2}}),$
(1)

where $\Delta$ is the cutoff length that separates the length scales.

To evolve the coarse field, we add a source term $\mathcal{F}$ to the Navier-
Stokes equations:

$\frac{\partial\mathbf{\tilde{u}}}{\partial
t}+\nabla\cdot\mathbf{\tilde{u}}\mathbf{\tilde{u}}=-\frac{1}{\rho}\nabla\tilde{p}+\nu\nabla^{2}\mathbf{\tilde{u}}+\mathcal{F}.$
(2)

We then estimate the time derivative by finding the difference in the velocity
at two-time levels and compute the source terms by:

$\mathcal{F}=\frac{\partial\mathbf{\tilde{u}}}{\partial
t}+\nabla\cdot\mathbf{\tilde{u}}\mathbf{\tilde{u}}+\frac{1}{\rho}\nabla\tilde{p}-\nu\nabla^{2}\mathbf{\tilde{u}},$
(3)

where $\tilde{p}$ is the filtered pressure. We note that this is different
from what is usually done for LES modeling using DNS results, where the
subgrid stresses are computed from the fully resolved results. While we have
the fully resolved flow field and could thus compute the subgrid stresses
directly in the usual way, our motivation is to test an approach applicable to
cases where we only have the coarse field. The extra term, $\mathcal{F}$, is
intended to ensure that equation (2) evolves the coarse field correctly.

###  II.2 Neural Network (NN) Architecture

The key assumption in LES is that the subgrid stresses can be related to the
coarse or large scale flow and we assume that the same is true for the source
term (even if it is not computed directly as the divergence of the stresses).
The source term cannot depend directly on the velocities since we can always
add a constant translation, so it must depend on the derivatives of the
velocity. For three-dimensional flows, [14] examined the relevance of a number
of variables as input for neural network for the subgrid stresses and
determined that the first and second derivative of the velocity were most
important. [15] determined the source term $\nabla\cdot\tau_{SGS}$ (which they
refer to as the SGS force) using the first derivatives, but included several
surrounding spatial points, instead of just the single point under
consideration. They found that working with the full source term, rather than
the stresses, resulted in a much better agreement with the filtered data.
Here, we assume that the source term depends on the local value of the first
and second derivatives: $\frac{\partial u}{\partial x}$, $\frac{\partial
u}{\partial y}$, $\frac{\partial v}{\partial x}$,
$\frac{\partial^{2}u}{\partial x^{2}}$, $\frac{\partial^{2}u}{\partial
x\partial y}$, $\frac{\partial^{2}u}{\partial y^{2}}$ and
$\frac{\partial^{2}v}{\partial x^{2}}$. Only 7 components are needed because
of the incompressibility condition. Because the flow is homogeneous and
without any directional preference, we can use the same expression for the two
components $\mathcal{F}_{x}$ and $\mathcal{F}_{y}$, since by rotating the
coordinate system by $90^{\circ}$ in counterclockwise direction, the new “x”
axis collapses with original “y” axis. Thus we can use the same mapping
function for the $y$ component, by slightly “rotating” the inputs. Thus

$\displaystyle\mathcal{F}_{x}=\mathcal{M}\Biggl{(}\frac{\partial u}{\partial
x},\frac{\partial u}{\partial y},\frac{\partial v}{\partial
x},\frac{\partial^{2}u}{\partial x^{2}},\frac{\partial^{2}u}{\partial
x\partial y},\frac{\partial^{2}u}{\partial
y^{2}},\frac{\partial^{2}v}{\partial x^{2}}\Biggr{)},$ (4)
$\displaystyle\mathcal{F}_{y}=\mathcal{M}\Biggl{(}\frac{\partial v}{\partial
y},-\frac{\partial v}{\partial x},-\frac{\partial u}{\partial
y},\frac{\partial^{2}v}{\partial y^{2}},-\frac{\partial^{2}v}{\partial
x\partial y},\frac{\partial^{2}v}{\partial
x^{2}},-\frac{\partial^{2}u}{\partial y^{2}}\Biggr{)},$ (5)

where $\mathcal{M}$ denotes the function found by the training. In this way we
preserve the rotational invariance property, are able to use the data for the
two components together in one big dataset, and train only once.

Figure 1: The structure of the Neural Network.

(a) $t^{*}=4.55$

(b) $t^{*}=22.73$

Figure 2: DNS vorticity at nondimensional times $t^{*}=4.55$ and $22.73$
(scaled by $\tau_{L}$).

The Neural Network is a multilayered structure, with a linear vector
transformation followed by a nonlinear activation function at each layer,
achieving complicated regression from inputs to outputs. Here, we use a Neural
Network with a relatively simple structure, containing one input layer with 7
neurons, 3 hidden layers with 20 neurons in each layer, and an output layer
with 1 neuron. Neurons in adjacent layers are fully connected to each other.
Figure 1 illustrates the structure of the Network. The transformation from
layer $l$ to layer $l+1$ can be expressed as

$\textbf{\em Y}^{l+1}=\sigma(\textbf{\em W}^{l}\textbf{\em X}^{l}+\textbf{\em
b}^{l}),$ (6)

where $\textbf{\em X}\in\mathbb{R}^{m\times 1}$ is the data in layer $l$,
$\textbf{\em Y}\in\mathbb{R}^{n\times 1}$ is the data in layer $l+1$,
$\textbf{\em W}\in\mathbb{R}^{n\times m}$ and $\textbf{\em
b}\in\mathbb{R}^{n\times 1}$ represents the weighting matrix and the bias term
that needs to be optimized. $\sigma(x)$ stands for the RELU activation
function, $\sigma(x)=max(x,0)$, used to introduce nonlinearity into the
mapping. For effective learning, the input data has been preprocessed in such
a way that $\textbf{\em X}^{0}=\frac{1}{\alpha}(\textbf{\em X}-\textbf{\em
X}_{mean})$, where $\textbf{\em X}_{mean}$ and $\alpha$ are the mean and
standard deviation of the raw data. The NN structure is built with Keras in
Python 3.7 and the mapping is optimized by stochastic gradient descent method
by monitoring the mean square error as the loss function. Based on our
experiment, 500 epochs are generally enough for finding the optimum and the
loss function does not decrease significantly with a further increase in the
number of epochs.

##  III Results and Discussion

Figure 3: Filtered vorticity field with filter width $\Delta=0.4$ at
$t^{*}=4.55$. Figure 4: Energy spectrum for the DNS field and the filtered DNS
field at $t^{*}=4.55$.

To generate a 2D turbulent flow in a doubly periodic domain, we start with
four perturbed shear layers that quickly break up by Kelvin-Helmholtz
instability into approximately isotropic flow as seen in figure 2, where the
vorticity is shown at two times. The domain size is $2\pi\times 2\pi$,
resolved by a $2048\times 2048$ grid. Due to the inverse energy cascade
characteristic of 2D turbulence, small vortices grab each other and form
larger ones, as shown in the evolution of the vorticity field. In this case,
turbulence energy dissipation rate is
$\epsilon=2\nu(\overline{S^{\prime}_{ij}}\overline{S^{\prime}_{ij}})=0.179$
and the Kolmogorov length scale is approximately
$\eta_{k}=(\nu^{3}/\epsilon)^{1/4}=0.009$. The Reynolds number defined by the
integral length scale is $Re_{L}=u_{RMS}L/\nu=2335$, where $u_{RMS}$ is the
RMS velocity $u_{RMS}=\sqrt{E_{k}}=3.28$ and the integral length scale can be
derived from that $L=({\pi/}{2u_{RMS}^{2}})\int\frac{E(k)}{k}\,dk=0.713$. The
large eddy turnover time scale is $\tau_{L}=L/u_{RMS}=0.22$, and the
simulation is run for over 45 large eddy turnover times.

Figure 5: The predicted value of the component of the source term versus the
true source terms. The correlation coefficients are shown on the top of the
frame. The $45^{\circ}$ red line shows a perfect fit.

###  III.1 NN learning

To generate a coarse flow field, we apply a Gaussian filter to the DNS flow
field so that the original $2048\times 2048$ field is coarsened to $128\times
128$ grid points, or $\Delta_{c}=16\Delta_{f}$, where $\Delta_{c}$ and
$\Delta_{f}$ denote the coarse grid spacing and the fine grid spacing,
respectively. The filter width in equation (1) is chosen to be $\Delta=0.4$,
which is about half the integral length scale $L$ and much larger than the
grid spacing of the coarse mesh ($\Delta_{c}=0.049$). Note that the filter
applied here is more aggressive than that in [17], who used
$\Delta_{c}=8\Delta_{f}$. The filter size is also much larger than the fine
grid spacing ($\Delta=130\Delta_{f}$), and more aggressive than in [15] where
$\Delta=32\Delta_{f}$. Figure 3 shows the filtered vorticity field on the
$128\times 128$ grid. We obtain the 1D energy spectrum
$E(k)=\frac{1}{2}|\hat{u^{\prime}}(k)|^{2}$, where $\hat{u^{\prime}}(k)$
represents the Fourier Transform of the velocity fluctuation field, by taking
the Fourier transform of the fluctuation along the x direction and average in
the y direction. Doing it in a reversed way (taking the Fourier transform
along the y direction and averaging in the x direction) gives similar results.
The energy spectrum for the fully resolved DNS solution and the filtered flow
are shown in figure 4. The spectrum shows the classical -3 power law in the
inertial range for 2D turbulence ([19]).

Here, we sample the data every $\Delta t^{*}=2.27$ from $t^{*}=2.27$ to
$t^{*}=27.27$ for the Neural Network training, giving a total of 12 filtered
DNS fields, each containing $128\times 128$ data points (total of $393,216$
samples by combining data for both components). Before the data was used to
train the NN, it was normalized by subtracting the mean and dividing by the
standard deviation, enforcing a mean and standard deviation of 0 and 1,
respectively. The data is split into a training set (90%) and a validation set
(10%) randomly in order to perform a cross-validation. Figure 5 shows the
predicted value of one component of the source term versus the true source
term computed from the data by equation (3). The correlation coefficients are
shown on the top of the frame. It achieves a value higher than 0.98 after 500
epochs. Increasing the number of epochs, as well as using different initial
guesses, has minimal impact on the results. The high correlation coefficient
in the cross-validation test is about the same as [15] found and higher than
what [14] found (0.7).

Figure 6: The vorticity contours at $t^{*}=11.36$. (a) Filtered DNS vorticity
field; (b) LES vorticity field with NN models; (c) Vorticity field on coarse
grid with no model.

The source term coming from the neural network is not a completely smooth
function and the noise is sometimes amplified when used in the modified
Navier-Stokes equations. To eliminate that we smooth the source term slightly,
using a box filter to smooth it over a 5 by 5 stencil on the course grid.
Using a Gaussian filter with the same cut-off length as used to coarsen the
fully resolved flow also gives similar results. We note that instabilities
were also seen by [18], [15] and [13]. In order to avoid the numerical
instabilities for 2D flows, [18] forced the subgrid stress to be only
dissipative at large scales by truncating the source term, assuming that the
small scales only cause dissipation of the kinetic energy. [20] pointed out
that forward- and backscatter are present in approximately equal amounts in
3D. This was also shown by [18], who found that even in 2D flows about half of
the grid points transfer energy from fine scales to large scales. Here, we
avoid excluding the backscatter mechanism. [15] prevented instabilities in 3D
by including artificial dissipation into the momentum equation to damp out the
high frequency fluctuations.

###  III.2 A posteriori tests

Studies of subgrid stresses in LES typically distinguish between a priori
comparison of predicted subgrid stresses with subgrid stresses computed by
filtering fully resolves results, and a posteriori comparison of how well the
flow field computed by LES matches the filtered fully resolved results. A
somewhat disappointing finding is that a successful a priori test does not
guarantee a good a posteriori comparison. In our case, where the closure terms
are computed directly from the coarse grained result, a priori test is
irrelevant. We are not trying to reproduce the subgrid stresses or the source
term computed from the filtered DNS solution, but to ensure that the coarsened
field evolves correctly. The a posteriori test on the other hand, where we
compare the evolution of the coarse grained field to the evolution of the
filtered fully resolved field, is of course the ultimate test.

####  III.2.1 The training data

The purpose of a subgrid model is to ensure that the coarse field reproduces
the evolution of the fully resolved field, filtered at each instance in time.
Since the flow is highly unsteady and slightly different initial conditions
diverge at long times, we do not expect the coarse solution to follow the
filtered field exactly, except initially. Thus, we can only do a detailed
comparison between the spatial fields at early times and must be content with
comparing statistical quantities at later times.

We use the filtered DNS field at $t^{*}=4.55$ (figure 3) as the initial
condition for the a posteriori test. The LES uses a $128\times 128$ grid, with
the NN predicted source terms $\mathcal{F}^{model}$ added on the right hand
side of the momentum equation, followed by a slight spatial smoothing of the
modeled source term field. We also ran an under-resolved DNS simulation on a
$128\times 128$ grid for comparison, also with the filtered fully resolved
solution at $t^{*}=4.55$ as initial conditions. Figure 6 shows the
instantaneous vorticity field as found by the filtered DNS, the LES and the
under-resolved DNS at $t^{*}=11.36$. The agreement between the filtered DNS
and the LES results is reasonably good and the coarse grid solution (under-
resolved DNS) is further away from the filtered solution. This suggests that
the NN model captures the evolution of the coarse filtered field, at least for
the first $6.81$ large eddy turnover times.

Figure 7: The 1D energy spectrum of the DNS field, the filtered DNS field, the
LES field and the coarse grid solution. Figure 8: The second order
nondimentional longitudinal structure function $D^{*}_{LL}$.

Since it is unrealistic to expect the coarse field obtained by filtering the
fully resolved results and the coarse field evolved by the model equations to
match completely at long times, we assess the performance of the model by
examining if statistical quantities are preserved. Here, we calculate the 1D
energy spectrum and the $2^{nd}$ order structure function for the filtered DNS
field, LES field with NN model, and the coarse grid solution with no models.
The 1D longitudinal structure function is defined as
$D_{LL}(r)=\overline{(u^{\prime}(\mathbf{x}+\mathbf{r})-u^{\prime}(\mathbf{x}))^{2}}$,
where $u^{\prime}$ denotes the horizontal velocity fluctuation. We average the
field at 11 different times and nondimensionalize the structure function by
$u^{2}_{RMS}$. Averaging with more samples or at other times does not
significantly change the structure function or the 1D spectrum. As shown in
figure 7, the energy spectrum obtained by the NN modeled LES follows the $-3$
power law well in the inertial range. The coarse grid solution is however, a
little bit away from the $-3$ power law at large scales. Figure 8 shows the 1D
$2^{nd}$ order nondimensional longitudinal structure function versus
$r/\Delta_{c}$. The NN modeled LES structure function agrees well with the
filtered fully resolved one, and is much closer than the coarse grid
solutions. We also plot the enstropy, defined as
$\Pi=\oint\omega\cdot\omega\,dS$, versus time for the three different cases in
figure 9. It can be seen that there is an obvious improvement in matching the
enstropy evolution by using the NN modeled LES compared to the under-resolved
solution. We note that it takes 25 days on a single Intel core to evolve the
DNS solution to $t^{*}=45.46$ but only less than half an hour for the NN
modeled LES.

Figure 9: The evolution of the enstropy for the filtered DNS, LES and an
under-resolved solution.

####  III.2.2 Non-training data

In the previous section we show that the Neural Network closure model works
well in terms of matching both the instantaneous flow field at early times as
well as the turbulence statistics at long times for the data used to train the
model. Here we apply the model to a different case generated by 4 vertical
shear layers separating flows with discontinuous velocities, still on a
$2048\times 2048$ grid for a domain size of $2\pi\times 2\pi$. The Reynolds
number, defined by the integral length scale, is about the same as for the
previous case, or $Re_{L}=2685$. The Kolmogorov length and the RMS velocity
fluctuations are $\eta_{k}=0.0082$ and $u_{RMS}=3.29$. The large eddy turnover
time is $\tau_{L}=0.24$. We choose the same filter size $\Delta=0.4$ and
coarsen the DNS fields onto $128\times 128$ grid points, since the NN model
has only been trained for one filter size.

We note that we keep the majority of the inputs to the NN within the range of
the training data, facilitating a nonlinear interpolation of the data. We
start the a posteriori test from the filtered DNS field at $t^{*}=4.17$, and
plot the instantaneous vorticity field at $t^{*}=10.42$ in figure 10, along
with the under-resolved solution on a $128\times 128$ grid. The LES vorticity
field shows a better agreement with the true field after nearly 7 large eddy
turnover times than the under-resolved one. Figure 11 shows that the LES
matches the energy spectrum well and is better than under-resolved solution in
terms of the structure function, shown in figure 12, and the evolution of the
enstropy in figure 13.

Figure 10: Contours at $t^{*}=10.42$ of (a) filtered DNS vorticity field (b)
LES vorticity field with NN model (c) under-resolved vorticity field on coarse
grid.

Figure 11: The 1D energy spectrum of the DNS field, the filtered DNS field,
the LES field and the coarse grid solution. Figure 12: The second order
nondimentional longitudinal structure function $D^{*}_{LL}$.

Figure 13: The evolution of the enstropy for the filtered DNS, the LES and the
under-resolved solution.

##  IV Conclusions

The main goal of the present paper is to explore how we can generate closure
terms for a reduced order, or coarse grained, model directly from the coarse
data, by finding what we need to add to the standard momentum equations so
that the coarse flow evolves in the same way as the filtered fully resolved
flow. The not unexpected conclusion is that it works, at least for the problem
examined here. Here we start with a fully resolved flow field and could, of
course, compute the standard subgrid stresses by filtering the fully resolved
velocities and the product of the fully resolved velocities separately and
subtract the product of the filtered velocities from the filtered product, as
is usually done (see [18]). However, if the fully resolved solution is not
available, or if the coarsening is carried out in such a highly nonlinear way
that is not possible to write down the analytical form of the source terms,
then the present approach should still work. When the subgrid stresses are
available they show up in the averaged momentum equations as fluxes and their
divergence ensures that no new momentum is added. This is not preserved when
we work with the source term directly, although we expect conservation to be
approximately satisfied. While the stresses can, in principle, be recovered by
integration at any given time, for flows in fully periodic domains those are
only specified up to a time dependent constant, potentially causing
difficulties with using stresses from different times for the learning. This
problem should not exist in wall-bounded flows, but we have not explored this
further at the present time. We also note that we took the subgrid stresses to
depend on any coarse quantity we could think of. It is likely that this number
can be reduced, but doing so is not needed for the purpose of this study.

## Appendix: A priori test

Although a priori tests are generally irrelevant to our approach, where the
closure terms are found directly from the coarse data, we expect the source
term computed using the coarse field by equation (3), to be the same as the
one computed from filtering the fully resolved results. Thus, we show in
figure A.1 a comparison of the source term predicted by the Neural Network and
the source term computed from the fully resolved flow at one time for the
training data. Note that the second one is not the source term used to train
the Neural Network—which is computed from the evolution of the filtered field
by equation (3). The agreement is reasonably good for both components, with a
mean square error (MSE) of 1.103 and 0.467 for the x and the y components
respectively. The MSE is defined by
$MSE=\frac{1}{N}\sum\|\mathcal{F}^{true}-\mathcal{F}^{modeled}\|^{2}$.

Figure A.1: A comparison of the instantaneous true $\mathcal{F}$ field (left)
and the predicted one (right) for both x (top) and y components (bottom) at
$t^{*}=36.36$.

###### Acknowledgements.

This research was supported in part by the National Science Foundation Grant
CBET-1953082. We also acknowledge the use of computational resources at the
Maryland Advanced Research Computing Center (MARCC).

## References

  * Piomelli and Zang [1991] U. Piomelli and T. A. Zang, Large-eddy simulation of transitional channel flow, Computer physics communications 65, 224 (1991).
  * Ling  _et al._ [2016] J. Ling, A. Kurzawski, and J. Templeton, Reynolds averaged turbulence modelling using deep neural networks with embedded invariance, Journal of Fluid Mechanics 807, 155 (2016).
  * Xiao  _et al._ [2016] H. Xiao, J.-L. Wu, J.-X. Wang, R. Sun, and C. Roy, Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier–Stokes simulations: A data-driven, physics-informed Bayesian approach, Journal of Computational Physics 324, 115 (2016).
  * Wang  _et al._ [2017] J.-X. Wang, J.-L. Wu, and H. Xiao, Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data, Physical Review Fluids 2, 034603 (2017).
  * Wu  _et al._ [2019] J. Wu, H. Xiao, R. Sun, and Q. Wang, Reynolds-averaged Navier–Stokes equations with explicit data-driven Reynolds stress closure can be ill-conditioned, Journal of Fluid Mechanics 869, 553 (2019).
  * Tracey  _et al._ [2015] B. D. Tracey, K. Duraisamy, and J. J. Alonso, A machine learning strategy to assist turbulence model development, in  _53rd AIAA aerospace sciences meeting_ (2015) p. 1287.
  * Weatheritt and Sandberg [2017] J. Weatheritt and R. D. Sandberg, Hybrid Reynolds-averaged/Large-Eddy Simulation methodology from symbolic regression: formulation and application, AIAA Journal 55, 3734 (2017).
  * Ma  _et al._ [2015] M. Ma, J. Lu, and G. Tryggvason, Using statistical learning to close two-fluid multiphase flow equations for a simple bubbly system, Physics of Fluids 27, 092101 (2015).
  * Ma  _et al._ [2016] M. Ma, J. Lu, and G. Tryggvason, Using statistical learning to close two-fluid multiphase flow equations for bubbly flows in vertical channels, International Journal of Multiphase Flow 85, 336 (2016).
  * Duraisamy  _et al._ [2019] K. Duraisamy, G. Iaccarino, and H. Xiao, Turbulence modeling in the age of data, Annual Review of Fluid Mechanics 51, 357 (2019).
  * Sarghini  _et al._ [2003] F. Sarghini, G. De Felice, and S. Santini, Neural networks based subgrid scale modeling in Large Eddy Simulations, Computers & fluids 32, 97 (2003).
  * Gamahara and Hattori [2017] M. Gamahara and Y. Hattori, Searching for turbulence models by artificial neural network, Physical Review Fluids 2, 054604 (2017).
  * Beck  _et al._ [2019] A. Beck, D. Flad, and C.-D. Munz, Deep neural networks for data-driven LES closure models, Journal of Computational Physics 398, 108910 (2019).
  * Wang  _et al._ [2018] Z. Wang, K. Luo, D. Li, J. Tan, and J. Fan, Investigations of data-driven closure for subgrid-scale stress in Large-Eddy Simulation, Physics of Fluids 30, 125101 (2018).
  * Xie  _et al._ [2020a] C. Xie, J. Wang, and E. Weinan, Modeling subgrid-scale forces by spatial artificial neural networks in Large Eddy Simulation of turbulence, Physical Review Fluids 5, 054606 (2020a).
  * Xie  _et al._ [2020b] C. Xie, J. Wang, H. Li, M. Wan, and S. Chen, Spatially multi-scale artificial neural network model for Large Eddy Simulation of compressible isotropic turbulence, AIP Advances 10, 015044 (2020b).
  * Maulik and San [2017] R. Maulik and O. San, A neural network approach for the blind deconvolution of turbulent flows, Journal of Fluid Mechanics 831, 151 (2017).
  * Maulik  _et al._ [2019] R. Maulik, O. San, A. Rasheed, and P. Vedula, Subgrid modelling for two-dimensional turbulence using neural networks, Journal of Fluid Mechanics 858, 122 (2019).
  * Kraichnan [1967] R. H. Kraichnan, Inertial ranges in two-dimensional turbulence, The Physics of Fluids 10, 1417 (1967).
  * Piomelli  _et al._ [1991] U. Piomelli, W. H. Cabot, P. Moin, and S. Lee, Subgrid-scale backscatter in turbulent and transitional flows, Physics of Fluids A: Fluid Dynamics 3, 1766 (1991).


# Towards Quantum Telescopes: Demonstration of a Two-Photon Interferometer for
Quantum-Assisted Astronomy

Jesse Crawford a  Denis Dolzhenko a  Michael Keach a  Aaron Mueninghoff b
Raphael A. Abrahao a  Julian Martinez-Rincon a  Paul Stankus a  Stephen
Vintskevichc  Andrei Nomerotski a  aBrookhaven National Laboratory, Upton NY
11973, USA
bStony Brook University, Stony Brook, NY 11794, USA
cIndependent researcher, Ras Al Khaimah, United Arab Emirates  *e-mails for
correspondence<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

Classical optical interferometery requires maintaining live, phase-stable
links between telescope stations. This requirement greatly adds to the cost of
extending to long baseline separations, and limits on baselines will in turn
limit the achievable angular resolution. Here we describe a novel type of two-
photon interferometer for astrometry, which uses photons from two separate sky
sources and does not require an optical link between stations. Such techniques
may make large increases in interferometric baselines practical, even by
orders of magnitude, with corresponding improvement in astrometric precision
benefiting numerous fields in astrophysics. We tested a benchtop analogue
version of the two-source interferometer and unambiguously observe correlated
behavior in detections of photon pairs from two thermal light sources, in
agreement with theoretical predictions. This work opens new possibilities in
future astronomical measurements.

##  1 Introduction and Basic Concepts

Classical optical Michelson interferometers collect photons from a sky source
into two or more sub-apertures, and these are then transported through optical
links and brought into interference at a common point. For a pair of sub-
apertures separated by a baseline $B$ the measured interference pattern is
sensitive to features of the source with angular size
$\delta\theta\sim\lambda/B$, where $\lambda$ is the photon wavelength. The
optical link between the stations must meet demanding requirements, with
lengths remaining stable to within a fraction of a wavelength. This makes the
Michelson interferometry expensive and difficult to extend to long baselines,
and this in turn limits the achievable resolution [1, 2].

The use of quantum optics to improve the precision of astronomical
measurements is a long-desired goal of both the optical and astronomical
communities. In particular, the seminal Gottesman-Jennewein-Croke (GJC)
proposal [3] attracted attention as a way to build a quantum-enhanced
telescope, i.e. a very-long-baseline interferometer enabled by the use of
quantum optical effects. However, the GJC proposal is dependent on the use of
quantum repeaters [4, 5, 6, 7], a technology which still requires a
substantial amount of development, limiting practical implementations of the
technique. Alternatively, the proposed Stankus-Nomerotski-Slosar-Vintskevich
(SNSV) approach [8] relies on enhancements of the Hanbury Brown and Twiss
(HBT) effect. It simplifies requirements for an optical quantum-assisted
astronomical measurement, providing a practical pathway to achieve the goal of
more precise astrometric measurements for astronomical objects, and is based
on the quantum interference phenomena, while avoiding employment of quantum
repeaters. Other theoretical techniques for improved resolution on imaging the
starlight were also introduced and some proof-of-principle demonstrations were
successfully performed [9, 10, 11, 12, 13].

Here we describe a proof-of-principle experiment based on the SNSV scheme [8,
14, 15], shown schematically in the left part of Figure 1. In brief, photons
are collected from two sky sources in two observing stations. Light from the
sources is directed into beamsplitters for interference. Correlations of
photon detections between the four beamsplitter outputs are sensitive to the
relative phase difference of incoming photons from the two sources. From this
the opening angle between the sources can be determined, allowing for longer
baselines when compared to the classical interferometry, since it entirely
removes the need for a physical optical path between the stations. Therefore,
such an approach could improve astrometric precision.

Figure 1: $(a)$ Concept of the proposed Stankus-Nomerotski-Slosar-Vintskevich
(SNSV) two-photon interferometer, which interferes and detects photons, shown
as flat waves, from two astronomical sources. $(b)$ Equivalent scheme for the
tabletop implementation. The argon lamps 1 and 2 indicate the two ports
associated with input spatial and polarization modes. Each output detection
port $D_{1},D_{2},D_{3},$ and $D_{4}$ corresponds to detectors in two
observation stations $L$ and $R$ in (a).

The measured coincidence rates of the outputs can be described by a second-
order intensity correlation function $\Gamma_{jk}$ which has the following
simplified form:

$\Gamma_{jk}=A+B\cos(\delta_{j}-\delta_{k})$ (1)

where A and B are the coefficients dependent on the photon polarization and
coherence properties, and $\delta_{j,k}$ are the phases for given modes.

We provide a detailed theoretical derivation of the intensity correlation
function in quasi-monochromatic approximation in the Supplement focusing on
the polarization degrees of freedom, and on the relationship between the mode
indistinguishability and interference effects. There, we predict that the
output channels from Figure 1 should form pairs such that channels 1&3 should
be anti-correlated with channel pairs 1&4 and 2&3 but correlated with 2&4,
while the channel pair 1&4 would be correlated with the channel pair 2&3 but
anti-correlated with 1&3 and 2&4\. We also predict that coincidences of the
channel pairs 1&2 and 3&4 would have a stationary rate without oscillations, a
non-trivial correlation pattern.

The basis of the proposed technique relies on two phenomena: the Hanbury
Brown-Twiss (HBT) [16] and Hong-Ou-Mandel (HOM) [17] effects. We expect the
two photons to bunch within their coherence time as in the HBT effect,
resulting in the HBT peaks in the time difference distributions. In addition
to this, we predict that the population of the HBT peak would depend on the
relative phase difference between the two photons. The HOM effect would play a
role when correlations between two outputs of the same beamsplitter are
considered, leading to two indistinguishable photons to coalesce at a
beamsplitter output. Note that the presented scheme has similarities with the
Noh-Fougères-Mandel experiments [18, 19, 20, 21] but, in contrast, it is
focused on the second-order correlation analysis, which describes coincidences
between pairs of detectors. The optical effects discussed above is a
manifestation of the interference of two indistinguishable photons. To be
indistinguishable the photons must have similar frequencies and arrive at the
beamsplitters at the same time. These arguments establish requirements on the
timing and spectral resolution for the interferometer instrumentation [14, 22,
23, 24].

The SNSV technique has parallels with the GJC approach [3], which employs a
path-entangled source of single photons distributed between two observing
stations, and its variations [25, 26, 27, 28, 29, 30]. In the SNSV proposal,
photons provided by the second thermal light source would replace photons from
a single photon source, typically a spontaneous parametric down conversion
(SPDC) source [31]. While an on-sky light source is inherently thermal, we
note that this has important practical advantages, compared to the above-
mentioned quantum source, such as the complete absence of the direct optical
link between the stations and uniformity of instrumentation for observations.

The two-photon interferometer presented here allows for precision improvement,
in principle, by orders of magnitude, which could benefit numerous fields in
cosmology and astrophysics. There are many scientific opportunities that would
benefit from substantial improvements in astrometric precision such as testing
theories of gravity by direct imaging of black hole accretion discs, precision
parallax for the cosmic distance ladder, mapping microlensing events, peculiar
motions and dark matter; see Ref. [8] for a more comprehensive discussion. As
a numerical example, it was determined in Ref. [8] that a nominal precision on
the order of 10 $\mu$as on the opening angle between two bright stars in a
single night’s observation could be reached.

In the following, we start with discussing the benchtop experimental setup of
the interferometer in Section 2. We then focus on the methods and results in
Section 3. Section 4 provides a discussion of the results. Lastly, we present
the conclusions and the future outlook in Section 5. Detailed theoretical
considerations and extended experimental setup description are presented in
the Supplement.

##  2 Experimental Setup

The experimental setup utilized in the measurements is shown in the right part
of Figure 1. Two argon lamps with isolated 794.82 nm spectral lines were used
as sources of thermal photons. The photons were directed to four 50:50 non-
polarizing beamsplitters to arrange interference as in the original scheme in
the left part of Figure 1. The photon phases remained stable to the
environmental disturbances for extended periods of time and could be adjusted
deterministically with specialized phase shifters to study the phase
dependence. The shifters were implemented as small angle glass wedges, which
can be moved laterally in fine steps. The four outputs of the interferometer
were instrumented with fast single-photon detectors. Two types of detectors
were used, the single photon avalanche detectors (SPAD) [32, 33, 34] and
superconducting nanowire single-photon detectors (SNSPD) [35, 36, 37], both
with temporal resolution of the order of 100 ps. The SPAD/SNSPD digital output
signals were then read out by a TDC module. A detailed description of the
setup is provided in the Supplement.

We employed the following experimental procedure. Firstly, one of the phase
shifters was moved in small steps by a distance of about 0.45 mm over the
duration of 15 minutes, which corresponded to a shift by five wavelengths,
followed by a pause of two minutes. Then, the second phase shifter was moved
in the same manner followed by another two-minute pause. The total duration of
the undisturbed measurement was about 35 minutes. The photon time-stamps of
four output channels were continuously logged on disk for post-processing. A
variety of runs were performed with two different detector types and varying
polarizer configurations. These configurations included experiments with
unpolarized photons, when the polarizers were removed from the beam paths, and
experiments with different respective polarizations for the two lamps.

##  3 Methods and Results

The main goal of the analysis was to determine the dependence of two-photon
correlations on the relative phase of the photons. Algorithms were developed
to condition the raw data by removing the afterpulses, then to find
coincidences of photon pairs in different channels to identify the HBT peaks,
and to determine the dependence of the peak population on the photon phase.

###  3.1 HBT Peaks

The HBT peaks appear in the time difference distributions of channel pairs due
to the two-photon interference yielding the photon bunching [16]. The
characteristic shape of the distribution is determined by the convolution of
the corresponding photon coherence time due to the spectral width of the argon
line, and timing resolution [14]. We study this effect by analyzing the
distribution of photon detection time differences by combining various
detector channels. The analysis algorithm searches for pairs of single photons
detected within a $\pm$ 20 ns window from each other.

In previous work [14], we also determined that the timing resolution is the
predominant contribution to the HBT peak temporal width, so we model the HBT
peaks with a normalized $g^{(2)}$ autocorrelation function of this form:

$g_{ij}^{(2)}(\Delta t,t_{0},p)=1+V_{HBT}e^{-\frac{\left(\Delta
t-t_{0}\right)^{2}}{2\sigma^{2}}},\ p\in\\{VV,VH\\},\
ij\in\\{12,13,14,23,24,34\\},$ (2)

where $V_{HBT}$ is the visibility of HBT effect, $\Delta t$ is the time
difference with an offset $t_{0}$, $\sigma$ is the standard deviation, $p$ is
the polarization, and indices $i$ and $j$ are labeling the detector pairs; see
Eqs. (S10) - (S14) in the Supplement for the detailed derivation.

As expected, prominent peaks are seen for most of the six combinations of
channels, as shown in Figure 2. The left and right parts of Figure 2 show the
HBT peaks with vertical-vertical, or VV, and vertical-horizontal, or VH, as
the orientation of the polarizers at the two inputs. Time offsets in the peak
positions are due to the varying paths in different channels including small
differences in the length of optical fibers. Note that for the VH case the
interference is happening for two photons from the same source and the VH
photon pairs from the different sources do not contribute. We discuss
properties of the HBT peak distributions in Section 4.

Figure 2: Normalized coincidence count rates of two-photon detections as a
function of the time difference between them, $\Delta$t, for different output
channel combinations and different input polarizations, as labeled. The peak
in each case indicates the enhanced correlation between two photons,
calibrating where simultaneous pairs will appear. Left: results with both
polarizers aligned vertically, called VV configuration. Right: results with
one polarizer aligned vertically and the other polarizer aligned horizontally,
called VH configuration. The parameters $t_{0}$ and $\sigma$ are the time
offset and the Gaussian fit width, respectively, according to Equation 2. The
peaks are signatures of the Hanbury Brown-Twiss (HBT) effect, appearing due to
the photon bunching.

###  3.2 Coincidence rates

As the first step, the HBT peaks were fit with a Gaussian function to
determine the peak width $\sigma$ and the central value of the entire 35 min
dataset. Then the coincidence rate for channel pairs was determined using two
techniques. In the first, simple approach a window of $\pm 1.5\sigma$ around
the HBT peak central value was selected, where $\sigma$ was taken from the
corresponding fit of the peak. Then the number of entries in this window in a
predefined time bin, typically 10 - 30 sec, was determined and plotted as a
function of time as shown in Figure 3 for all six channel combinations for the
SNSPD data set with polarizers. The left and right parts of Figure 3 show,
respectively, the number of coincidences in the 20 sec time bins with VV and
VH polarizer orientation at the interferometer inputs. It is expected that
these time trending plots will have oscillatory behavior in the VV case due to
the advance of the photon phase caused by the phase shifters as predicted in
Equation 1. We indeed observe this behavior but defer a detailed discussion of
the main features of these measurements to Section 4.

---
Figure 3:  Number of two-photon detections within a $\pm$ 1.5$\sigma$ window
about $t_{0}$ of the HBT peak in 20 sec time bins plotted versus time for the
SNSPD data set with polarizers. The time accounts for the real time during
which the phase shifters were slowly moving a preset distance corresponding to
five wavelengths, then paused for two minutes, then moved again for five
wavelengths. The graphs show results for six combinations of channel pairing,
1&2, 1&3, 1&4, 2&3, 2&4 and 3&4\. Left: Results with VV configuration of
polarizers. Right: Results with VH configuration of polarizers. See the text
for discussion.

In a more sophisticated approach, the HBT peak in each predefined, sequential
time bin was fit with a Gaussian profile, only allowing to vary the peak
amplitude and background detection level. The peak center $t_{0}$ and standard
deviation $\sigma$ of these fits were held constant at the values determined
from the initial fit of the overall HBT peak. Then the area under each of
these Gauss peaks within $\pm 1.5\sigma$ of the peak center was determined,
effectively subtracting the background detection rate, and plotted as a
function of time. We show an example of such fit in Figure 4. Mathematically
this result can be described as the convolution of the second-order
correlation function with a filtering function, see Eqs.(S13) and (S14) in
Supplement. All fits in this method were made using LMFIT [38], which was also
used to explicitly calculate uncertainties on the parameter best-fit values.
This technique should give better statistical accuracy since it uses more
statistics for the flat background and can also account for a slow drift of
the flat background.

To determine the visibility and relative phase, the oscillatory behavior of
trending plots due to the phase evolution was fit with a cosine function:

$\Gamma_{ip_{i}jp_{j}}(t)=\langle A\rangle_{ip_{i}jp_{j}}+\langle
B\rangle_{ip_{i}jp_{j}}\cos{\left(\frac{2\pi}{T}t-\Delta\delta_{ij}\right)}.$
(3)

where $\langle B\rangle$ denotes the signal amplitude, $T$ is the
characteristic period of slow phase adjustments by the phase shifters, $t$ is
the data time stamp, $\Delta\delta_{ij}\equiv\delta_{i}-\delta_{j}$ is the
constant relative phase, and $\langle A\rangle$ is the background level.
Labels $p_{i}\in\\{H,V\\}$ and $p_{j}\in\\{H,V\\}$ are labels of polarization
modes tailored to spatial modes $i$ and $j$ respectively. We define the
corresponding visibility as:

${\mathcal{V}_{OSC}}=\frac{\langle B\rangle_{ip_{i}jp_{j}}}{\langle
A\rangle_{ip_{i}jp_{j}}}.$ (4)

Figure 4: Two-photon coincidences count rates for the oscillations for
channels 1&3 in the SNSPD data set with with VV polarization fit to Equation
3. The coincidences rates were determined fitting a Gaussian peak in each 10
sec time bin. Data points are presented together with one standard deviation
error bars.

##  4 Analysis and Discussion

Oscillations in the trending plots were fit using a cosine function as
explained above, from which the main parameters were extracted. These
parameters included the oscillation period, visibility and phase. Out of
these, the latter two, visibility and relative phase, are critical parameters,
unambiguously predicted by the theory.

###  4.1 Oscillations in Coincidence Rate

We show the coincidence rates for all channel combinations in Figure 3. As
previously shown in Equation 1, the channel pairs on the opposite
interferometer arms (pairs 1&3, 1&4, 2&3, and 2&4) are expected to oscillate
in or out of phase, while channel pairs on the same interferometer arm (pairs
1&2, 3&4) are not expected to oscillate. Since for the HOM effect only
indistinguishable photons are supposed to interfere, we do not expect any
oscillations in the VH configuration of polarizations.

Only the first five oscillation periods were used for the cosine fits in these
results extracting the relative phase, which was calculated with respect to
the channel pair 1&4 phase, and the visibility, calculated using Eq. 4. This
analysis was performed for several datasets: VV polarized and unpolarized
datasets with SNSPD detectors, and unpolarized dataset with SPAD detectors. We
summarize the measured visibility for the cross-station channel combinations
and for the datasets in the top part of Figure 5 and the relative phases, in
radians, in the bottom part of Figure 5. The measured parameters all behave as
expected. As can be seen the visibility was higher for the polarized dataset,
as expected, by approximately a factor of two compared to the unpolarized
dataset. In the unpolarized case, the visibility of oscillations for the SPAD
detectors were consistently smaller compared to the visibility of the peaks
for the SNSPD detectors. This can be explained by the inferior timing
resolution of SPAD detectors.

Figure 5: Top: Visibility of oscillations for the two-photon coincidences
count rates. Bottom: Relative phases, in radians, calculated with respect to
the channel pair 1&4 phase. Colors correspond to different datasets. Pol and
noPol are the VV poparized and unpolarized datasets, respectively. Data points
are presented together with one standard deviation error bars.

We also checked the behavior of the interferometer for varying polarizations
of the two beams. Figure 6 shows the visibilities of two-photon detections as
a function of polarization angle, for channel pairs 1&3, 1&4, 2&3, and 2&4\.
Here, one input polarization angle was varied in 22.5∘ increments from 0∘ to
180∘ starting from either VV or VH configuration, while the other polarizer
was not moved. The dependence on the relative angle for the two polarizations
was changing in anti-phase for these two cases, as expected. We compare the
experimental data to the theoretical predictions, as in Eq. 5, derived in
Supplement, and find a good agreement.

$\displaystyle{V_{VV}(\theta)}=\frac{r_{VV}\sin^{2}(\theta)}{1+r_{VV}\sin^{2}{\theta}+r_{VV}^{2}\sin^{4}{\theta}+r_{HV}^{2}\cos^{4}{\theta}+r_{HV}\cos^{2}{\theta}+r_{HV}r_{VV}\sin^{2}{\theta}\cos^{2}{\theta}}$
$\displaystyle{V_{VH}(\theta)}=\frac{r_{HH}\cos^{2}(\theta)}{1+r_{HH}\cos^{2}{\theta}+r_{HH}^{2}\cos^{4}{\theta}+r_{VH}^{2}\sin^{4}{\theta}+r_{VH}\sin^{2}{\theta}+r_{HH}r_{VH}\sin^{2}{\theta}\cos^{2}{\theta}},$
(5)

where $r_{ij}$ are the normalized rates defined in Supplement.

Figure 6: Visibility of two-photon detections as a function of polarization
angle, for channel pairs 1&3, 1&4, 2&3, and 2&4\. Here, one input polarization
angle was varied in 22.5∘ increments from 0∘ to 180∘ starting from either VV
or VH configuration, while the other polarizer was not moved. Different colors
refer to different configurations as indicated. Data points are presented
together with one standard deviation error bars. The measured experimental
data is fitted by Eq. 5. To fit the results we added to the fit function a
small constant offset, which described the accidental coincidences counts.

###  4.2  HBT peak visibility, cancellation of HOM and HBT effects

As already discussed in Section 3 the HBT peaks were fit with a Gaussian
function. The resulting timing resolution ($\sigma$) was found to be equal to
280 ps and 140 ps, respectively for SPAD and SNSPD detectors. This is in
agreement with expected timing resolution of those detectors assuming a time
difference measurement for two independent photons. The HBT peak visibility
was introduced in Eq. 2. A graph summarizing the visibilities for all
combinations of channel pairs in the main datasets is shown in Figure 7.
Similarly to the oscillation visibility, the HBT visibility was significantly
higher for the polarized case, decreasing for the unpolarized case and for the
SPAD detectors, which have worse timing resolution.

Figure 7:  Visibility of the HBT peaks presented in Figure 2 for various
datasets. Colors correspond to different datasets. Pol and noPol are the VV
poparized and unpolarized datasets, respectively. Data points are presented
together with one standard deviation error bars.

We also note that the HBT peak visibility is almost zero for the 1&2 and 3&4
pair combinations in the VV configuration. These particular channels evaluate
coincidences of two photons exiting the two opposite sides of the same
beamsplitter so their coincidence rate will have a dip due to the Hong-Ou-
Mandel (HOM) effect. This will cancel out the photon bunching peak due to the
HBT effect. However, the results in Figures 2 and 7 show it’s not an exact
cancellation, which is likely due to the non-ideal equalization of the
interferometer arms, also a valid explanation for the spread of visibility for
the middle four channel combinations in Figures 5 and 7.

##  5 Conclusions and Outlook

In this work we implemented a proof-of-principle demonstration of the SNSV
proposal for quantum-assisted astronomy [8] and described bench-top
experiments with a two-photon interferometer. The results, in particular, the
observed phase dependence, confirm the predicted functionality of the proposed
instrument, suggesting that it is a viable experimental approach that can
improve the astrometric precision.

The next steps in the exploration of the approach would be observations of on-
sky light sources. We note that the phase shifting performed in the above
experiments and corresponding oscillations in the coincidence rates are
directly analogous to the Earth rotation fringe shifting [39]. Due to Earth’s
rotation, the effective baseline between the two stations change, which
induces changes in the interference patterns, in particular, in the frequency
of corresponding fringe oscillations. The frequency value is predicted in [8]
to be proportional to the opening angle between the stars. The uncertainty of
the fringe rate has more favorable scaling with observation period than simple
photon statistic methods, so it can be considered a promising observable for
the quantum astrometry approach.

The star spectra are typically broadband, so another obvious extension of the
technique is spectroscopic binning, which would allow us to estimate the
observables in numerous spectral bins. Each bin would act as an independent
experiment, so the sensitivity of the interferometer would improve as
$\sqrt{N}$, where $N$ is the number of spectral bins. The sensitivity also
improves with a larger number of stations and better timing resolution, as
discussed in Ref. [8].

Though we operate here with the thermal states of light, which have the
classical Gaussian quasiprobability distribution, we consider a quantum
description of these phenomena, as in this work, to be very instructive. It
emphasizes the role of mode (or path) indistinguishability in the quantum
interference phenomena. Moreover, this description can be extended further by
employing the Quantum Continuous Variable (CV) formalism [40]. The CV
formalism will play an important role in expanding the presented technique to
multiple observing stations assisted by auxiliary ground-based quantum states
such as squeezed states and non-trivial quantum channels [41, 42, 43, 44, 45,
46], providing novel opportunities in astronomy for extraction of valuable
information. At the same time, the complexity of analytical descriptions and
data processing would grow tremendously with an increased number of stations
and auxiliary quantum states involved in the measurements, including entangled
ones [3, 25, 26], so an important future research goal is to provide a
theoretical description of such expansions based on both entities, the
Gaussian states and quantum channels, powered by machine-learning methods [47,
48].

In summary, we built and characterized a tabletop prototype of a two-photon
interferometer, which could improve the astrometric precision by orders of
magnitude, by means of enabling extra long baselines between observing
stations. The implemented interferometer allowed us to test the important
features of the SNSV proposal, in particular to demonstrate that the relative
phase of two photons from two independent thermal sources has a direct effect
on their bunching due to the HBT effect. The approach demonstrated here is
technically feasible with existing technologies of single photon detection
[14, 39] and allows us to move towards measurements with on-sky sources. This
work represents a major step towards quantum-assisted astronomy.

##  6 Acknowledgments

This work was supported by the U.S. Department of Energy QuantISED award and
BNL LDRD grants 19-30 and 22-22. A.M. acknowledges support under the Lourie
Fellowship from the Stony Brook University Department of Physics and
Astronomy. We are grateful to Jonathan Schiff and Rom Simovitch for the
software development and support.

##  7 Disclosures

The authors declare no conflicts of interest.

## References

  * [1] E. Pedretti, J. D. Monnier, T. ten Brummelaar, and N. D. Thureau, “Imaging with the CHARA interferometer,” New Astronomy Reviews 53, 353–362 (2009).
  * [2] T. A. ten Brummelaar, H. A. McAlister, S. T. Ridgway, J. W. G. Bagnuolo, N. H. Turner, L. Sturmann, J. Sturmann, D. H. Berger, C. E. Ogden, R. Cadman _et al._ , “First results from the CHARA array. II. a description of the instrument,” The Astrophysical Journal 628, 453–465 (2005).
  * [3] D. Gottesman, T. Jennewein, and S. Croke, “Longer-baseline telescopes using quantum repeaters,” Phys. Rev. Lett. 109, 070503 (2012).
  * [4] N. Sangouard, C. Simon, H. de Riedmatten, and N. Gisin, “Quantum repeaters based on atomic ensembles and linear optics,” Rev. Mod. Phys. 83, 33–80 (2011).
  * [5] W. J. Munro, K. Azuma, K. Tamaki, and K. Nemoto, “Inside quantum repeaters,” IEEE Journal of Selected Topics in Quantum Electronics 21, 78–90 (2015).
  * [6] E. K. Carlson, “The key device needed for a quantum internet,” Physics 13, 104 (2020).
  * [7] D. Gottesman, “Quantum telescopes,” in _Optical and Infrared Interferometry and Imaging VII,_ vol. 11446 P. G. Tuthill, A. Mérand, and S. Sallum, eds., International Society for Optics and Photonics (SPIE, 2020), p. 1144615.
  * [8] P. Stankus, A. Nomerotski, A. Slosar, and S. Vintskevich, “Two-photon amplitude interferometry for precision astrometry,” The Open Journal of Astrophysics 5 (2022).
  * [9] M. Tsang, “Quantum nonlocality in weak-thermal-light interferometry,” Phys. Rev. Lett. 107, 270402 (2011).
  * [10] M. Tsang, “Resolving starlight: a quantum perspective,” Contemporary Physics 60, 279–298 (2019).
  * [11] M. E. Pearce, E. T. Campbell, and P. Kok, “Optimal quantum metrology of distant black bodies,” Quantum 1, 21 (2017).
  * [12] L. A. Howard, G. G. Gillett, M. E. Pearce, R. A. Abrahao, T. J. Weinhold, P. Kok, and A. G. White, “Optimal imaging of remote bodies using quantum detectors,” Phys. Rev. Lett. 123, 143604 (2019).
  * [13] Z. S. Tang, K. Durak, and A. Ling, “Fault-tolerant and finite-error localization for point emitters within the diffraction limit,” Opt. Express 24, 22004–22012 (2016).
  * [14] A. Nomerotski, P. Stankus, A. Složar, S. Vintskevich, S. Andrewski, G. A. Carini, D. Dolzhenko, D. England, E. V. Figueroa, S. Gera, J. Haupt, S. Herrmann, D. Katramatos, M. Keach, A. Parsells, O. Saira, J. Schiff, P. Svihra, T. Tsang, and Y. Zhang, “Quantum-assisted optical interferometers: instrument requirements,” in _Optical and Infrared Interferometry and Imaging VII,_ A. Mérand, S. Sallum, and P. G. Tuthill, eds. (2020), Proc. SPIE.
  * [15] M. Keach, S. Bellavia, Z. Chen, J. Crawford, D. Dolzhenko, E. Figueroa, A. Mueninghoff, A. Nomerotski, J. Schiff, R. Simovitch, A. Slosar, P. Stankus, and S. Vintskevich, “Increasing baselines and precision of optical interferometers using two-photon interference effects,” in _Optical and Infrared Interferometry and Imaging VIII,_ A. Mérand, S. Sallum, and J. Sanchez-Bermudez, eds. (SPIE, 2022).
  * [16] R. H. Brown and R. Twiss, “A test of a new type of stellar interferometer on sirius,” Nature 178, 1046–1048 (1956).
  * [17] C. K. Hong, Z. Y. Ou, and L. Mandel, “Measurement of subpicosecond time intervals between two photons by interference,” Phys. Rev. Lett. 59, 2044–2046 (1987).
  * [18] J. W. Noh, A. Fougères, and L. Mandel, “Measurement of the quantum phase by photon counting,” Phys. Rev. Lett. 67, 1426–1429 (1991).
  * [19] J. W. Noh, A. Fougères, and L. Mandel, “Operational approach to the phase of a quantum field,” Phys. Rev. A 45, 424–442 (1992).
  * [20] J. W. Noh, A. Fougères, and L. Mandel, “Further investigations of the operationally defined quantum phase,” Phys. Rev. A 46, 2840–2852 (1992).
  * [21] J. W. Noh, A. Fougères, and L. Mandel, “Measurements of the probability distribution of the operationally defined quantum phase difference,” Phys. Rev. Lett. 71, 2579–2582 (1993).
  * [22] A. Nomerotski, M. Keach, P. Stankus, P. Svihra, and S. Vintskevich, “Counting of hong-ou-mandel bunched optical photons using a fast pixel camera,” Sensors 20 (2020).
  * [23] K. M. Jordan, R. A. Abrahao, and J. S. Lundeen, “Quantum metrology timing limits of the hong-ou-mandel interferometer and of general two-photon measurements,” Phys. Rev. A 106, 063715 (2022).
  * [24] F. Bouchard, A. Sit, Y. Zhang, R. Fickler, F. M. Miatto, Y. Yao, F. Sciarrino, and E. Karimi, “Two-photon interference: the hong–ou–mandel effect,” Reports on Progress in Physics 84, 012402 (2020).
  * [25] E. T. Khabiboulline, J. Borregaard, K. De Greve, and M. D. Lukin, “Quantum-assisted telescope arrays,” Phys. Rev. A 100, 022316 (2019).
  * [26] E. T. Khabiboulline, J. Borregaard, K. De Greve, and M. D. Lukin, “Optical interferometry with quantum networks,” Phys. Rev. Lett. 123, 070504 (2019).
  * [27] M. Brown, V. Thiel, M. Allgaier, M. Raymer, B. Smith, P. Kwiat, and J. Monnier, “Interferometry-based astronomical imaging using nonlocal interference with single-photon states,” in _Frontiers in Optics $+$ Laser Science 2021,_ (OSA, 2021).
  * [28] M. Brown, V. Thiel, M. Allgaier, M. Raymer, B. Smith, P. Kwiat, and J. Monnier, “Long-baseline interferometry using single photon states as a non-local oscillator,” in _Quantum Computing, Communication, and Simulation II,_ P. R. Hemmer and A. L. Migdall, eds. (SPIE, 2022).
  * [29] D. Diaz, Y. Zhang, V. O. Lorenz, and P. G. Kwiat, “Emulating quantum-enhanced long-baseline interferometric telescopy,” in _Frontiers in Optics $+$ Laser Science 2021,_ (Optica Publishing Group, 2021), p. FTh6D.7.
  * [30] M. Brown, V. Thiel, M. Allgaier, M. Raymer, B. Smith, P. Kwiat, and J. Monnier, “Proof-of-principle laboratory demonstration of long-baseline interferometric imaging using distributed single-photons,” in _Quantum 2.0 Conference and Exhibition,_ (Optica Publishing Group, 2022), p. QM3C.1.
  * [31] A. Christ, A. Fedrizzi, H. Hübel, T. Jennewein, and C. Silberhorn, “Parametric down-conversion,” in _Experimental Methods in the Physical Sciences,_ vol. 45 (Elsevier, 2013), pp. 351–410.
  * [32] L. Gasparini, B. Bessire, M. Unternahrer, A. Stefanov, D. Boiko, M. Perenzoni, and D. Stoppa, “SUPERTWIN: towards 100kpixel CMOS quantum image sensors for quantum optics applications,” in _Quantum Sensing and Nano Electronics and Photonics XIV,_ M. Razeghi, ed. (2017), Proc. SPIE.
  * [33] M. Perenzoni, L. Pancheri, and D. Stoppa, “Compact spad-based pixel architectures for time-resolved image sensors,” Sensors 16 (2016).
  * [34] M. Lee and E. Charbon, “Progress in single-photon avalanche diode image sensors in standard CMOS: From two-dimensional monolithic to three-dimensional-stacked technology,” Japanese Journal of Applied Physics 57, 1002A3 (2018).
  * [35] A. Divochiy, F. Marsili, D. Bitauld, A. Gaggero, R. Leoni, F. Mattioli, A. Korneev, V. Seleznev, N. Kaurova, O. Minaeva, G. Gol'tsman, K. G. Lagoudakis, M. Benkhaoul, F. Lévy, and A. Fiore, “Superconducting nanowire photon-number-resolving detector at telecommunication wavelengths,” Nature Photonics 2, 302–306 (2008).
  * [36] D. Zhu, M. Colangelo, C. Chen, B. A. Korzh, F. N. C. Wong, M. D. Shaw, and K. K. Berggren, “Resolving photon numbers using a superconducting nanowire with impedance-matching taper,” Nano Letters (2020).
  * [37] B. Korzh, Q. Zhao, J. P. Allmaras, S. Frasca, T. M. Autry, E. A. Bersin, A. D. Beyer, R. M. Briggs, B. Bumble, M. Colangelo, G. M. Crouch, A. E. Dane, T. Gerrits, A. E. Lita, F. Marsili, G. Moody, C. Peña, E. Ramirez, J. D. Rezac, N. Sinclair, M. J. Stevens, A. E. Velasco, V. B. Verma, E. E. Wollman, S. Xie, D. Zhu, P. D. Hale, M. Spiropulu, K. L. Silverman, R. P. Mirin, S. W. Nam, A. G. Kozorezov, M. D. Shaw, and K. K. Berggren, “Demonstration of sub-3 ps temporal resolution with a superconducting nanowire single-photon detector,” Nature Photonics 14, 250–255 (2020).
  * [38] M. Newville, R. Otten, A. Nelson, A. Ingargiola, T. Stensitzki, D. Allan, A. Fox, F. Carter, Michał, R. Osborn, D. Pustakhod, lneuhaus, S. Weigand, Glenn, C. Deil, Mark, A. L. R. Hansen, G. Pasquevich, L. Foks, N. Zobrist, O. Frost, A. Beelen, Stuermer, azelcer, A. Hannum, A. Polloreno, J. H. Nielsen, S. Caldwell, A. Almarza, and A. Persaud, “lmfit/lmfit-py: 1.0.3,” Zenodo (2021).
  * [39] Z. Chen, A. Nomerotski, A. Složar, P. Stankus, and S. Vintskevitch, “Astrometry in two-photon interferometry using earth rotation fringe scan,” arXiv preprint arXiv:2205.09091 (2022).
  * [40] A. Serafini, _Quantum continuous variables: a primer of theoretical methods_ (CRC press, 2017).
  * [41] Y. Xia, W. Li, W. Clark, D. Hart, Q. Zhuang, and Z. Zhang, “Entangled sensor networks empowered by machine learning,” in _Optical Fiber Communication Conference (OFC) 2021,_ (Optica Publishing Group, 2021), p. Th3A.4.
  * [42] S. V. Vintskevich, D. A. Grigoriev, and S. N. Filippov, “Effect of an incoherent pump on two-mode entanglement in optical parametric generation,” Phys. Rev. A 100, 053811 (2019).
  * [43] A. J. Brady, C. Gao, R. Harnik, Z. Liu, Z. Zhang, and Q. Zhuang, “Entangled sensor-networks for dark-matter searches,” PRX Quantum 3, 030333 (2022).
  * [44] B.-H. Wu, S. Guha, and Q. Zhuang, “Entanglement-assisted multi-aperture pulse-compression radar for angle resolving detection,” arXiv preprint arXiv:2207.10881 (2022).
  * [45] Q. Zhuang and Z. Zhang, “Physical-layer supervised learning assisted by an entangled sensor network,” Phys. Rev. X 9, 041023 (2019).
  * [46] A. Cox, Q. Zhuang, C. Gagatsos, B. Bash, and S. Guha, “Transceiver designs to attain the entanglement assisted communications capacity,” arXiv preprint arXiv:2208.07979 (2022).
  * [47] M. Krenn, J. Landgraf, T. Foesel, and F. Marquardt, “Artificial intelligence and machine learning for quantum technologies,” arXiv:2208.03836 (2022).
  * [48] S. Vintskevich, N. Bao, A. Nomerotski, P. Stankus, and D. Grigoriev, “Classification of four-qubit entangled states via machine learning,” arXiv preprint arXiv:2205.11512 (2022).


# Modeling of Asymptotically Periodic Outbreaks:
a long-term SIRW2 description of Covid-19?

Alex Viguerie
Division of Mathematics
Gran Sasso Science Institute
Viale Francesco Crispi 7
L’Aquila, AQ 67100, Italy
<EMAIL_ADDRESS>
&Margherita Carletti
Department of Pure and Applied Sciences
University of Urbino C. Bo
Piazza della Repubblica 13
Urbino (Pu) 61029, Italy
<EMAIL_ADDRESS>
&Alessandro Veneziani
Department of Mathematics
Department of Computer Science
Emory University
400 Dowman Drive
Atlanta, GA 30322, USA
<EMAIL_ADDRESS>
&Guido Silvestri
Yerkes Division of Microbiology and Immunology
Yerkes National Primate Research Center
Yerkes Nonhuman Primate Genomics Core
Emory University
Atlanta, GA 30329, USA
<EMAIL_ADDRESS>

###### Abstract

As the outbreak of COVID-19 enters its third year, we have now enough data to
analyse the behavior of the pandemic with mathematical models over a long
period of time. The pandemic alternates periods of high and low infections, in
a way that sheds a light on the nature of mathematical model that can be used
for reliable predictions. The main hypothesis of the model presented here is
that the oscillatory behavior is a structural feature of the outbreak, even
without postulating a time-dependence of the coefficients. As such, it should
be reflected by the presence of limit cycles as asymptotic solutions. This
stems from the introduction of (i) a non-linear waning immunity based on the
concept of immunity booster (already used for other pathologies); (ii) a fine
description of the compartments with a discrimination between individuals
infected/vaccinated for the first time, and individuals already
infected/vaccinated, undergoing to new infections/doses. We provide a proof-
of-concept that our novel model is capable of reproducing long-term
oscillatory behavior of many infectious diseases, and, in particular, the
periodic nature of the waves of infection. Periodic solutions are inherent to
the model, and achieved without changing parameter values in time. This may
represent an important step in the long-term modeling of COVID-19 and similar
diseases, as the natural, unforced behavior of the solution shows the
qualitative characteristics observed during the COVID-19 pandemic.

##  1 Introduction

It has been two years now that the world experienced the COVID-19 pandemic. We
have enough data to identify some patterns in the evolution of the outbreak
over time. In particular, in many situations it is evident that the outbreak
features an oscillatory behavior, showing an alternating behavior of high and
low infected populations. For example, in Fig. 1, we plot the infected curve
in Italy (source: [1]). We also plot a Fourier spectral analysis of the data,
showing that the dominant frequencies are around 1 year-1 and three months-1.
Whether this is the effect of a seasonal infection-rate, government measures,
or a structural nature of the pandemic is currently an open question.

The calibration of mathematical models with a predictive purpose can be widely
educated by the answer to this question. If we assume that the periodicity is
a structural feature of the pandemic, a mathematical model should reflect this
in the nature of its asymptotic solutions. On the contrary, if the periodicity
is the result of seasonality, or government responses to the pandemic (such as
lockdowns and masking mandates), the model should incorporate this effect in
its parameter calibration.

In the majority of work modeling COVID-19, this latter assumption is employed.
Such works typically adopt classical compartmental models like SIR
(Susceptible-Infected-Recovered, see [2, 3] for detailed discussion), and its
derivations. These models may lead to excellent short-term predictions;
however, as they do not have periodic solutions as asymptotic time limit, they
will ultimately fail to yield a periodic solution when using time-independent
parameters - they do not have periodic asymptotic solutions [4]. In order to
recover the observed periodic behavior, then, such models must employ time-
dependent parameter fitting, generally attributing the periodicity to
underlying changes in population behavior, seasonality, government measures,
and other factors [5, 6, 7, 8, 9]. In contrast, in order to recover a truly
periodic solution without regular re-parameterization, the presence of limit-
cycles becomes a mandatory aspect of long-term modeling.

It is worth noting that the two different possible answers (nature of the
problem vs. time-dependence of the parameters) are not mutually exclusive.
Looking at the spectral analysis of the Italian data, we argue that it may be
a combination of different effects on long and short time scales.

Figure 1: Infected population in Italy over the time range Feb 24 2020- Mar 9
2022. Top: linear scale, center: log-scale. We observe a clear oscillatory,
periodic dynamics. Bottom: Fourer spectral analysis showing the dominant
frequencies (1/370 days and 1/148 days).

In this paper, we focus on a model where the wavy pattern is primarily
intrinsic to the infection qualitatively, acknowledging that behavioral,
seasonal, and other external aspects (e.g., higher infectivity in Winter,
lower in Summer, different infectivity of different variants, lockdowns)
impact the pattern quantitatively.

Our starting point is the concept of waning immunity advocated for other
pathologies [10]. Waning Immunity has been already considered for COVID-19,
see, e.g., [11, 12, 13, 14] with different mathematical descriptions. The
basic idea is that the immunity induced by recovery or vaccination declines
over time. This was largely demonstrated by the cases of
reinfection/infection-after-vaccination documented for the COVID-19 outbreak.
However, the reduction of the immunity does not follow a simple pattern,
because there is an immunity booster effect: an individual with a declining
immunity exposed to an infected person may get a beneficial effect,
reinforcing (or boosting) her/his immune system. This leads to a model
sometimes called SIRW, where the compartment W (Waning) added to the classical
SIR denotes the “weakly immune” individuals that may go back to the
Susceptible compartment or get an immunity booster by contact with infected or
vaccination. It has been demonstrated by a geometric approach [10] that this
model does have limit cycles for particular combinations of the parameters.
One of the important aspects of the parameter calibration is the “multiscale-
in-time” nature of the model, where the natural growth of the population
(according to a Malthusian model) is, in fact, a slower process compared to
the infection/vaccination rate.

In this brief note, we propose a new extension of this model. We further
stratify all four subpopulations of the SIRW into “never previously exposed”
and “previously exposed" subpopulations. The rationale is that, while a
vaccinated or previously infected individual may lose their immunity from
infection over time, some baseline protection is nonetheless maintained; the
illness is less likely to be severe and is more likely to pass quickly. A
similar concept in the case of pertussis (with no immunity booster) was
advocated in [15]. As the model is formally a duplication of the four SIRW
compartments (plus the Deceased), we call this model SIRW2 (SIRWsquare).

The main purpose of this communication is to provide a proof-of-concept that
the SIRW2 model can actually describe the wavy behavior even in the absence of
a time-dependent calibration of the parameters. The presented work presents
preliminary results that suggest that this model can provide a reliable long-
term description of the COVID-19 outbreak.

The document is organized as follows. We first introduce the waning immunity
and explain the SIRW model (studied in [10]), in order to describe the immune-
boosting process and provide the motivation behind the SIRW2 model. We then
introduce the full SIRW2 model, providing the model equations, a description
of the dynamics they are intended to model, and a complete list of the
relevant parameters and their purpose. We then numerically demonstrate with
constant parameters that this model may in principle reproduce wave-like
dynamics. We conclude with a series of perspectives on future steps in terms
of both theoretical analysis and validation. We aim at testing the hypothesis
that the periodic dynamics of COVID-19 pandemic is intrinsic to the nature of
the infection, so to be properly described (and quantitatively predicted) by
the SIRW2 model.

##  2 The Mathematical Model

###  2.1 The SIRW model

The SIRW model shown in [10] considers the susceptible $s$, infected $i$,
recovered $r$, and waning $w$ compartments and reads as follows:

$\displaystyle d_{t}s$ $\displaystyle=-\beta si+\kappa w,$ (1) $\displaystyle
d_{t}i$ $\displaystyle=\beta si-\phi i,$ (2) $\displaystyle d_{t}r$
$\displaystyle=\phi i-\xi r+\chi iw,$ (3) $\displaystyle d_{t}w$
$\displaystyle=\xi r-\chi iw-\kappa w.$ (4)

A full list of the variables and parameters in (1)-(4) is provided in Table 1.

The model (1)-(4) describes the passage of individuals throughout the four
compartments. A susceptible individual ($s$), upon contact with an infected
individual $i$, may become infected (the rate of which is controlled by the
contact rate $\beta$ and proportional to the product $si$), moving into the
compartment $i$. Infected individuals recover with a rate $\phi$, moving into
the $r$ compartment, where they acquire immunity from infection. At a rate of
$\xi$, recovered individuals begin to lose their immunity, entering the waning
compartment $w$.

While the immunity is waning, two things can happen:

  1. (i)

The waning individual can come into contact with an infected individual again,
with contact rate $\chi$, boosting their immunity, returning them to the $r$
compartment, or;

  2. (ii)

The waning individual, with no additional contact with the virus, loses their
immunity and returns to the $s$ compartment at rate $\kappa$.

In [10], the model (1)-(4) was shown to admit periodic-orbit solutions in
time, suggesting a correspondence between the modeled immune-boosting process
and the possible oscillatory dynamics. Thus, we believe this model, and the
underlying immune-boosting mechanism, may provide a key towards a more
realistic model of COVID-19, and the immune process generally.

Parameter/variable | Name | Units
---|---|---
$s$ | Susceptible individuals | Persons
$i$ | Infected individuals | Persons
$r$ | Recovered individuals | Persons
$w$ | Waning-immunity individuals | Persons
$m$ | Deceased individuals | Persons
$\beta$ | Contact rate between $s$ and $i$ | Persons${}^{-1}\cdot$ Days-1
$\phi$ | Recovery rate for | Days-1
$\xi$ | Loss-of-immunity initiation rate | Days-1
$\chi$ | Immunity boosting-contact rate | Persons${}^{-1}\cdot$ Days-1
$\kappa$ | Loss-of-immunity rate | Days-1
Table 1: Relevant variables, parameters, parameter names, and units for the
systems (1)-(4).

###  2.2 The SIRW2 model

While the more detailed immunological process described in model (1)-(4)
provides a sensible starting point for the modeling of COVID-19, preliminary
numerical simulations show that the oscillatory behavior is monotonically
damped, a circumstance that is not corresponded by data. In fact, the model
(1)-(4) does not consider vaccinations and does not consider mortality from
the disease. The addition of several new terms for vaccination and mortality,
as well as an additional compartment for deceased is sufficient for their
inclusion.

The more important modification necessary, however, is that the loss of
immunity is not a complete and total elimination of all immunity. In reality,
individuals maintain some baseline of immunological protection, as explained
above, in the form of, for instance, Memory B- cells [16]. These protections
both protect against severe disease, and lead to a faster recovery process
upon reinfection, when compared to infection in an individual with no prior
immunological protection.

To this end, we introduce an extension of (1)-(4) incorporating these crucial
effects: vaccination, disease mortality, and distinction between populations
with and without prior immunological exposure. To model vaccination, we add a
vaccination term, which can both grant and boost immunity (i.e. enable the
movement $s\to r$ and $w\to r$). For the addition of mortality, we add another
compartment, denoted $m$ (for ‘mortality’), to track deceased individuals.
Finally, we stratify the populations based on whether or not they have prior
immunological protection, either from vaccination or prior infection. We also
incorporate a Malthus demographic independent of the disease into the model.
The extended model, which we denote by SIRW2 (SIRW-squared) model, reads:

$\displaystyle d_{t}s_{u}$
$\displaystyle=-\beta_{u}^{u}s_{u}i_{u}-\beta_{u}^{e}s_{u}i_{e}-\alpha
s_{u}+\psi n-\rho s_{u},$ (5) $\displaystyle d_{t}i_{u}$
$\displaystyle=\beta_{u}^{u}s_{u}i_{u}+\beta_{u}^{e}s_{u}i_{e}-\phi_{u}i_{u}-\gamma_{u}i_{u}-\rho
i_{u},$ (6) $\displaystyle d_{t}r_{u}$ $\displaystyle=\phi_{u}i_{u}+\alpha
s_{u}-\xi_{u}r_{u}+\chi_{u}^{e}w_{u}i_{e}+\chi_{u}^{u}w_{u}i_{u}+\alpha
w_{u}-\rho r_{u},$ (7) $\displaystyle d_{t}w_{u}$
$\displaystyle=\xi_{u}r_{u}-\chi_{u}^{e}w_{u}i_{e}-\chi_{u}^{u}w_{u}i_{u}-\kappa_{u}w_{u}-\alpha
w_{u}-\rho w_{u},$ (8) $\displaystyle d_{t}s_{e}$
$\displaystyle=-\beta_{e}^{u}s_{e}i_{u}-\beta_{e}^{e}s_{e}i_{e}+\kappa_{u}w_{u}+\kappa_{e}w_{e}-\alpha
s_{e}-\rho s_{e},$ (9) $\displaystyle d_{t}i_{e}$
$\displaystyle=\beta_{e}^{u}s_{e}i_{u}+\beta_{e}^{e}s_{e}i_{e}-\phi_{e}i_{e}-\gamma_{e}i_{e}-\rho
i_{e},$ (10) $\displaystyle d_{t}r_{e}$
$\displaystyle=\phi_{e}i_{e}-\xi_{e}r_{e}+\chi_{e}^{e}w_{e}i_{e}+\chi_{e}^{u}w_{e}i_{u}+\alpha
w_{e}-\rho r_{e},$ (11) $\displaystyle d_{t}w_{e}$
$\displaystyle=\xi_{e}r_{e}-\chi_{e}^{e}w_{e}i_{e}-\chi_{e}^{u}w_{e}i_{u}-\kappa_{e}w_{e}-\alpha
w_{e}-\rho w_{e},$ (12) $\displaystyle d_{t}m$
$\displaystyle=\gamma_{u}i_{u}+\gamma_{e}i_{e}.$ (13)

Here, $n=s_{u}+i_{u}+r_{u}+w_{u}+s_{e}+i_{e}+r_{e}+w_{e}$, is the size of the
living population. The subscript $u$ indicates ‘unexposed’, meaning this class
of individuals has not been previously infected with, or vaccinated against,
the disease. The subscript $e$ indicates ‘previously exposed’, meaning that
this class of individuals has been infected previously, or vaccinated against,
the disease.

Verbally, the model (5)-(13) describes the population of persons who have
never been infected with the disease $s_{u}$, who may become infected via
contact with either a first-time infected ($i_{u}$) or reinfected ($i_{e}$)
person. After moving to the first-infection compartment $i_{u}$, a person may
recover (with rate $\phi_{u}$) or die (rate $\gamma_{u}$). Assuming recovery,
the movement is into the compartment $r_{u}$, indicating recovery with
immunity. Alternatively, an individual may be vaccinated (with rate $\alpha$)
and move directly from $s_{u}$ to $r_{u}$ without any previous infection.

Once in the the recovered state $r_{u}$, an individual begins to lose
immunity, and, at a rate of $\xi_{u}$, the individual moves to the waning-
state $w_{u}$. In the waning-state, a few things can happen:

  1. (i)

An individual can come into contact with an infected individual (either
$i_{u}$ or $i_{e}$, with contact rate $\chi_{u},\,\chi_{e}$ respectively,
boosting the immunity and returning to the compartment $r_{u}$;

  2. (ii)

An individual can receive a vaccine booster (rate $\alpha$), also boosting
their immunity and returning to $r_{u}$;

  3. (iii)

In the absence of additional contact with the virus, in the form of exposure
to an infected individual or a vaccine booster, the immunity can be lost with
a rate $\kappa_{u}$.

In this last case, the switch is to the compartment $s_{e}$. From there, the
model progression works similarly to in the first-infection (sub-scripted with
$u$) case, in the movement between the $i_{e}$, $r_{e}$, and $w_{e}$ states.
However, the previously-exposed infection progression differs from the first-
time infection progression in that the relevant parameters are different. For
most diseases we postulate that:

  1. (iv)

The mortality rate is substantially lower: $\gamma_{e}$ << $\gamma_{u}$;

  2. (v)

The recovery rate is substantially higher: $\phi_{e}$ >> $\phi_{u}$.

This is in line with observation and general consensus. Our model evidence
also suggests that, perhaps:

  1. a)

Immunity is lost more quickly after reinfection: $\xi_{e}$ > $\xi_{u}$;

  2. b)

however, it is also boosted more easily among those with prior immunity:
$\chi_{e}^{e}>\chi_{u}^{e}$, $\chi_{e}^{u}>\chi_{u}^{u}$.

A full list and short description of the different variables, parameters, and
their units is given in Table (2).

Parameter/variable | Name | Units
---|---|---
$s_{u,e}$ | Susceptible individuals (no prior immunity, prior immunity) | Persons
$i_{u,e}$ | Infected individuals (no prior immunity, prior immunity) | Persons
$r_{u,e}$ | Recovered individuals (no prior immunity, prior immunity) | Persons
$w_{u,e}$ | Waning-immunity individuals (no prior immunity, prior immunity) | Persons
$m$ | Deceased individuals | Persons
$\beta_{u}^{u,e}$, $\beta_{e}^{u,e}$ | Contact rate between $s_{u}$ and $i_{u},i_{e}$, $s_{e}$ and $i_{u},i_{e}$. | Persons${}^{-1}\cdot$ Days-1
$\phi_{u,e}$ | Recovery rate for $i_{u}$, $i_{e}$ | Days-1
$\xi_{u,e}$ | Loss-of-immunity initiation rate for $r_{u}$, $r_{e}$ | Days-1
$\chi_{u}^{u,e}$, $\chi_{e}^{u,e}$ | Immunity boosting rates for $r_{u}$, $r_{e}$ | Persons${}^{-1}\cdot$ Days-1
$\kappa_{u,e}$ | Loss-of-immunity rate for $w_{u}$, $w_{e}$ | Days-1
$\gamma_{u,e}$ | Mortality rate from disease, for $i_{u}$, $i_{e}$ | Days-1
$\alpha$ | Vaccination rate | Days-1
$\psi$ | Natality rate | Days-1
$\rho$ | General mortality rate | Days-1
Table 2: Relevant variables, parameters, parameter names, and units for the
system (5)-(13).

##  3 Numerical results

We show a brief demonstration of the model’s capability by running it for a
hypothetical scenario. We take parameter values as shown in Table 3. We note
that these values do not necessarily correspond to any particular outbreak
values; they are designed to showcase the type of solution that the model
(5)-(13) can admit, in particular on a long term. We simulate the model over a
two-year period. The model is implemented in MATLAB [17]. Since the the
underlying equations are stiff, we utilize the ode23s solver.

Parameter/variable | Units | Value
---|---|---
$s_{u_{0},e_{0}}$ | Persons | 9999, 0
$i_{u_{0},e_{0}}$ | Persons | 8, 0
$r_{u_{0},e_{0}}$ | Persons | 0, 0
$w_{u_{0},e_{0}}$ | Persons | 0, 0
$m$ | Persons | 0
$\beta_{u}^{u,e}$ | Persons${}^{-1}\cdot$ Days-1 | 1.54e-5, 4e-5
$\beta_{e}^{u,e}$ | Persons${}^{-1}\cdot$ Days-1 | 5.5e-5, 5.5e-5
$\phi_{u,e}$ | Days-1 | .13, .14
$\xi_{u,e}$ | Days-1 | .005, .005
$\chi_{u}^{u,e}$ | Persons${}^{-1}\cdot$ Days-1 | 7.7.e-5, 2.e-4
$\chi_{e}^{u,e}$ | Persons${}^{-1}\cdot$ Days-1 | 2.8e-4, 2.8e-4
$\kappa_{u,e}$ | Days-1 | .005, .005
$\gamma_{u,e}$ | Days-1 | .003, 4.4e-5
$\alpha$ | Days-1 | .002
$\psi$ | Days-1 | 5e-5
$\rho$ | Days-1 | 5e-5
Table 3: Parameter values for the shown simulation.

The results of the simulation are shown in Fig. 2. There are several important
characteristics to note: the first is that, as shown in the second and third
figures on the top row, the solution exhibits a periodic behavior, similar to
the observed dynamic. The second is that, despite the total number of
infections remaining high, in the third wave, we see a marked drop in
mortality; this is due to larger percentages of the population having gained
immunity, either through prior exposure to the disease, vaccination, or both.

In Fig. 3, we show the number of patients hospitalized with COVID-19 (left)
and total deceased from COVID-19 (right) in Italy [1]. We note that large
changes in testing protocols throughout the pandemic make a direct use of
case-counts unreliable, and hence we consider hospitalizations, as this data
is less subject to outside influence. We note that both the hospitalization
and death data display the same general qualitative trends as in the
simulation; the hospitalization data shows a clear periodic behavior, and the
death curve shows intermittent periods of stalling and growth, with the
overall trend flattening in time. We believe that, as in our model, this
reduction in severe illness, despite continuing high levels of the disease, is
due to the large amount of previously acquired immunity in the population. We
may expect that, even in the presence of continued periodic surges, the
mortality rate will remain low.

We stress that the agreement is, at this stage, only qualitative and very
preliminary. A proper parameter fitting is required in order for the model to
produce meaningful predictions. However, the results shown indicate that,
given such a parameter fitting, the proposed model is capable of naturally
producing the complex periodic dynamics observed during COVID-19, even without
reparameterization. As such dynamics are beyond the reach of most similar
models employed for COVID-19, we believe that such a result, despite its
preliminary nature, is significant. A fine tuning of the parameters should
include possible seasonality and lockdown dynamics.

Figure 2: Simulation results with the model (5)-(13) with the parameters found
in Table 2. With constant parameter values, the model is able to reproduce the
qualitative behavior of the COVID-19 pandemic. Figure 3: Number of
hospitalized persons (left) and deaths (right) in Italy. We see that, when
compared to Fig. 2, many of the key features are preserved: the qualitative
behavior of the infection and deceased curves show a good resemblance with the
observed behavior of the COVID-19 pandemic.

##  4 Conclusions and Perspectives

In this brief note, we have introduced an SIRW2 model, related to the models
shown in [10, 15]. The proposed model is intended to describe a more
sophisticated immune-system process than that found in most models: by the
incorporation of immune boosting, and by stratification in terms of persons
with and without previous exposure to the disease (either by vaccination or
prior infection). We then demonstrated, through proof-of-concept simulations,
that this modeling immunity may naturally recreate wave-like behaviors. This
is consistent with the long-term results of COVID-19 outbreak, without
requiring time-dependent parameter adjustment. Qualitatively, we obtain good
agreement with observed COVID-19 infection dynamics, providing strong evidence
that, after a full parameter-fitting process, the proposed model may be well-
suited for COVID-19, as well as other infectious diseases showing a similar
behavior. Much work has to be done for further developments. In particular,
the theory of the model must be further developed, including the analysis of a
suitable definition of $R_{0}$ and $R_{t}$, stability and equilibrium
solutions, and bifurcation analyses. Rigorous derivations describing the model
behavior (for instance, the oscillation period) in terms of the model
parameters is also important. Most importantly, an application of the model on
real-world data, incorporating a rigorous parameter fitting process, is
necessary as a quantitative validation. Finally, this modeling framework can
of course be developed even further via the incorporation of additional
factors, such as time-lags, hospitalizations, differentiation between
symptomatic and asymptomatic infections, and other important epidemiological
considerations.

## References

  * [1] Protezione Civile Italiana.  Covid-19.  https://github.com/pcm-dpc/COVID-19, 2020.
  * [2] William Ogilvy Kermack and Anderson G McKendrick.  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(772):700–721, 1927.
  * [3] Dimitri Breda, Odo Diekmann, WF De Graaf, A Pugliese, and R Vermiglio.  On the formulation of epidemic models (an appraisal of Kermack and McKendrick).  Journal of biological dynamics, 6(sup2):103–117, 2012.
  * [4] Herbert W Hethcote.  The mathematics of infectious diseases.  SIAM review, 42(4):599–653, 2000.
  * [5] Alex Viguerie, Guillermo Lorenzo, Ferdinando Auricchio, Davide Baroli, Thomas JR Hughes, Alessia Patton, Alessandro Reali, Thomas E Yankeelov, and Alessandro Veneziani.  Simulating the spread of COVID-19 via a spatially-resolved susceptible–exposed–infected–recovered–deceased (SEIRD) model with heterogeneous diffusion.  Applied Mathematics Letters, 111:106617, 2021.
  * [6] Malú Grave, Alex Viguerie, Gabriel F Barros, Alessandro Reali, and Alvaro LGA Coutinho.  Assessing the spatio-temporal spread of COVID-19 via compartmental models with diffusion in Italy, USA, and Brazil.  Archives of Computational Methods in Engineering, 28(6):4205–4223, 2021.
  * [7] Sheng Zhang, Joan Ponce, Zhen Zhang, Guang Lin, and George Karniadakis.  An integrated framework for building trustworthy data-driven epidemiological models: Application to the COVID-19 outbreak in New York City.  PLoS computational biology, 17(9):e1009334, 2021.
  * [8] Nicola Parolini, Luca Dede’, Paola F Antonietti, Giovanni Ardenghi, Andrea Manzoni, Edie Miglio, Andrea Pugliese, Marco Verani, and Alfio Quarteroni.  Suihter: A new mathematical model for COVID-19. application to the analysis of the second epidemic outbreak in Italy.  Proceedings of the Royal Society A, 477(2253):20210027, 2021.
  * [9] Mohamed Aziz Bhouri, Francisco Sahli Costabal, Hanwen Wang, Kevin Linka, Mathias Peirlinck, Ellen Kuhl, and Paris Perdikaris.  COVID-19 dynamics across the US: A deep learning study of human mobility and social behavior.  Computer Methods in Applied Mechanics and Engineering, 382:113891, 2021.
  * [10] Hildeberto Jardón-Kojakhmetov, Christian Kuehn, Andrea Pugliese, and Mattia Sensi.  A geometric analysis of the SIS, SIRS and SIRWS epidemiological models.  Nonlinear Analysis: Real World Applications, 58:103220, 2021.
  * [11] Nursanti Anggriani, Meksianis Z Ndii, Rika Amelia, Wahyu Suryaningrat, and Mochammad Andhika Aji Pratama.  A mathematical COVID-19 model considering asymptomatic and symptomatic classes with waning immunity.  Alexandria Engineering Journal, 61(1):113–124, 2022.
  * [12] Thomas Crellen, Li Pi, Emma L Davis, Timothy M Pollington, Tim CD Lucas, Diepreye Ayabina, Anna Borlase, Jaspreet Toor, Kiesha Prem, Graham F Medley, et al.  Dynamics of SARS-CoV-2 with waning immunity in the UK population.  Philosophical transactions of the royal society b, 376(1829):20200274, 2021.
  * [13] Daniel M Altmann and Rosemary J Boyton.  Waning immunity to SARS-CoV-2: implications for vaccine booster strategies.  The Lancet Respiratory Medicine, 9(12):1356–1358, 2021.
  * [14] Elie Dolgin et al.  COVID vaccine immunity is waning-how much does that matter.  Nature, 597(7878):606–607, 2021.
  * [15] Michiel van Boven, Hester E de Melker, Joop FP Schellekens, and Mirjam Kretzschmar.  Waning immunity and sub-clinical infection in an epidemic model: implications for pertussis in The Netherlands.  Mathematical Biosciences, 164(2):161–182, 2000.
  * [16] Brian J Laidlaw and Ali H Ellebedy.  The germinal centre B cell response to SARS-CoV-2.  Nature Reviews Immunology, 22(1):7–18, 2022.
  * [17] MATLAB.  version R2019b.  The MathWorks Inc., Natick, Massachusetts, 2019.