# Current bounds on baryogenesis from complex Yukawa couplings of light
fermions

Shahaf Aharony Shapira<EMAIL_ADDRESS>Department of Particle
Physics and Astrophysics,
Weizmann Institute of Science, Rehovot, Israel 7610001

###### Abstract

We calculate the contribution to the baryon asymmetry of the Universe (BAU)
from a $\mathcal{CP}$-violating source of the light quarks (charm, strange,
down, up) and the electron, resulting from a dimension-six effective field
theory term. We derive relevant bounds from the electric dipole moments of the
electron and neutron to estimate the maximal contribution from each single
flavor modification. Current bounds show that the charm quark can generate at
most $\mathcal{O}(1\%)$ of the BAU, while the lighter quarks and the electron
contribute at much lower levels.

BAU

    baryon asymmetry of the Universe
SM

    Standard Model
NP

    New Physics
EDM

    Electric Dipole Moment
SMEFT

    Standard Model effective field theory
CPV

     $\mathcal{CP}$-violating
VEV

    vacuum expectation value
LHC

    Large Hadron Collider
eEDM

    electron-Electric Dipole Moment
nEDM

    Neutron-Electric Dipole Moment
dim-$6$

    dimension-six
dim-$4$

    dimension-four
C.L.

    confidence level

## Introduction

The BAU is defined and measured Tanabashi:2018oca  to be

$\displaystyle Y_{B}$
$\displaystyle\equiv\frac{n_{B}-n_{\bar{B}}}{s}\approx\left(8.6\pm
0.1\right)\cdot 10^{-11}\equiv Y_{B}^{\text{obs}}\,,$ (1)

where $n_{(\bar{B})B}$ is the (anti-)baryon number density and $s$ is the
entropy density of the Universe. A non-vanishing value can be either the
result of initial conditions, or dynamically generated during the early
Universe. The former requires fine tuning and is inconsistent with inflation.
The latter, which is the more acceptable mechanism to address the asymmetry,
is called Baryogenesis.

There are three necessary conditions, known as the Sakharov conditions
Sakharov:1967dj , that are required from any theory in order to explain such
an imbalance: Baryon number violation, $\mathcal{C}$-symmetry and
$\mathcal{CP}$-symmetry violation and interactions out of thermodynamic
equilibrium. Although the Standard Model (SM) meets all three criteria, the
rate at which it contributes is far too small to account for the observed
baryon asymmetry due to two factors: the smooth crossover of the electroweak
phase transition and the suppression from the Kobayashi-Maskawa (KM) mechanism
of $\mathcal{CP}$ violation. Thus, if the baryon asymmetry was generated via
electroweak baryogenesis, the electroweak phase transition had to be strongly
first order and a new sources of $\mathcal{CP}$ violation must exist at, or at
least not far above, the electroweak scale.

New Physics (NP) beyond the SM is highly motivated by several open questions
in Physics (e.g. dark matter, neutrino masses). However, despite the efforts
made to discover new particles, none were found up to the TeV scale. It is
then plausible, that NP is above the electroweak scale, and thus could be
integrated out. This allows us to use Standard Model effective field theory
(SMEFT) tools to explore higher order terms, without being model-dependent.

We add a $\mathcal{CP}$-violating (CPV) phase using a dimension-six (dim-$6$)
coupling of three Higgs fields to the SM charged fermions. The BAU is then
proportional to the CPV source, which could be constrained by the Electric
Dipole Moment (EDM) of both the electron and neutron and by the Higgs boson
decay. This was previously done for the third generation particles
Fuchs:2020uoc  and the muon Fuchs:2019ore . Of that list, it was shown that
the $\tau$ is the only sole-contributor that can provide the entire observed
value of the BAU. We generalize this procedure to evaluate the contribution of
the light SM quarks (charm, strange, down, up) and the electron and discuss
its results here.

This Letter is organized as follows. First, we describe the SMEFT framework,
including the complex dim-$6$ term and the CPV source it generates. We then
outline key points in the process of electroweak baryogenesis, which are
formulated by the two-step approach via the transport equations followed by
the sphaleron process. Next, we present our numerical results for the
contribution of a single flavor to the BAU. The contribution is later bounded
using the experimental measurements of the electron and neutron EDMs and
various Higgs boson decay processes. Finally, we discuss our results and
conclusions.

##  SMEFT framework

We examine the implications of adding the following effective dim-$6$ terms to
the SM,

$\displaystyle\mathcal{L}_{\text{eff Yuk}}$
$\displaystyle=-\left(y_{f}+\frac{|H|^{2}}{\Lambda^{2}}\left(X^{R}_{f}+i\gamma^{5}X^{I}_{f}\right)\right)\overline{\psi_{Lf}}\psi_{Rf}H+\text{h.c.}\,,$
(2)

where $y_{f}$ is the dimension-four (dim-$4$) Yukawa coupling, $H$ is the SM
Higgs field $H\sim\left(1,\ 2\right)_{+1/2}$, $\Lambda$ is the NP scale, $X$
is the dim-$6$ Wilson coefficient and $\psi$ is a SM fermion. In our notation,
the lower index, $f$, denotes the flavor whereas the upper index distinguishes
between the real ($R$) and imaginary ($I$) coefficients. We find it useful to
define

$\displaystyle T^{R,I}_{f}$
$\displaystyle\equiv\frac{v^{2}}{2\Lambda^{2}}\frac{X^{R,I}_{f}}{y_{f}}\,,$
(3)

where $v=246$ GeV is the vacuum expectation value (VEV) of the Higgs
background field. Accordingly, the mass ($m$) and effective Yukawa coupling
($\lambda$) of each flavor can be defined and explicitly written as

$\displaystyle\mathcal{L}_{\text{eff Yuk}}$ $\displaystyle\supset-
m_{f}\bar{\psi}_{Lf}\psi_{Ri}-\lambda_{f}\bar{\psi}_{Lf}\psi_{Rf}h+\text{h.c.}\,,$
(4) $\displaystyle m_{f}$
$\displaystyle=\frac{vy_{f}}{\sqrt{2}}\left(1+T^{R}_{f}+iT^{I}_{f}\right)\,,$
$\displaystyle\lambda_{f}$
$\displaystyle=\frac{y_{f}}{\sqrt{2}}\left(1+3T^{R}_{f}+i3T^{I}_{f}\right)\,.$
(5)

In the mass basis, where $m_{f}\in\mathbb{R}$, they are given by

$\displaystyle m_{f}^{\prime}$
$\displaystyle=\frac{vy_{f}}{\sqrt{2}}\sqrt{\left(1+T^{R}_{f}\right)^{2}+{T^{I}_{f}}^{2}}\,,$
$\displaystyle\lambda_{f}^{\prime}$
$\displaystyle=\frac{y_{f}}{\sqrt{2}}\frac{\left(T^{R}_{f}+1\right)\left(3T^{R}_{f}+1\right)+T^{I}_{f}\left(3T^{I}_{f}+2i\right)}{\sqrt{\left(1+T^{R}_{f}\right)^{2}+{T^{I}_{f}}^{2}}}\,.$
(6)

Throughout this work we will use the complex parameter $\kappa_{f}$,
representing the deviation from the SM in the mass basis,

$\displaystyle\kappa_{f}$
$\displaystyle\equiv\frac{\lambda_{f}^{\prime}v}{m_{f}^{\prime}}=\frac{\left(T^{R}_{f}+1\right)\left(3T^{R}_{f}+1\right)+T^{I}_{f}\left(3T^{I}_{f}+2i\right)}{\left(1+T^{R}_{f}\right)^{2}+{T^{I}_{f}}^{2}}\,,$
$\displaystyle\kappa_{f}^{I}$
$\displaystyle\equiv\imaginary(\kappa_{f})=\frac{2T^{I}_{f}}{\left(1+T^{R}_{f}\right)^{2}+{T^{I}_{f}}^{2}}\,.$
(7)

This notation is convenient since the baryon asymmetry is proportional to the
CPV source ($S_{f}$) and therefore linear in $\kappa_{f}^{I}$ Fuchs:2020uoc ,

$\displaystyle Y_{B}\propto
S_{f}\propto\imaginary(m_{f}^{*}m_{f}^{\prime})\propto\kappa_{f}^{I}\,.$ (8)

As we show below, $\kappa_{f}^{I}$ can be constrained by CPV observables, such
as the EDMs, and by Higgs-related measurements in colliders.

## Electroweak baryogenesis

Adding a complex effective Yukawa coupling to one of the fermions of the SM
will introduce a CPV source (the explicit form is given in Eqn. (• ‣ A)). The
source depends on the background Higgs boson field $\left(\phi_{b}\right)$,
which acquires a VEV. Assuming the electroweak phase transition is strongly
first order, e.g. assuming $\phi_{b}$ is described by the kink solution (Eqn.
(37)), results with a CPV source which peaks mostly inside the non-vanishing
VEV bubble. A schematic plot is given in Fig. 1. Namely, the CPV source
generates a chiral asymmetry, mainly inside the broken phase. That asymmetry
can be transformed into a an abundance of baryons via the weak sphaleron,
which is a non-perturbative effect of the SM. The rate of the weak sphaleron
process is given by $\Gamma_{ws}\sim e^{-\langle\phi_{b}\rangle/T}T$, where
$T$ is the temperature. Although the rate is exponentially suppressed inside
the bubble, it is fast outside the bubble, during the early Universe.
Therefore, for baryogenesis to occur at the electroweak phase transition, the
chiral asymmetry should have propagated outside the bubble and into the
symmetric phase. Finally, as the bubble continues to expand, it eventually
captures the resulting baryon asymmetry.

Figure 1: Schematic description of the kink solution and its resulting CPV
source along the distance from the bubble wall ($z=0$). Blue: the kink
solution for the background Higgs boson field, as given in Eqn. (37), orange:
the resulting source, normalized.

## Two-step approach

The dim-$6$ term described in Eqn. (2) will affect the dynamics of the number
densities, which are defined as the difference between the number densities of
particles and anti-particles. This effect is the first step of the approach
and is described via a set of transport equations. We generalize the set given
in deVries:2018tgs  to include all the fermions of the SM in addition to the
Higgs boson,

$\displaystyle\partial_{\mu}U_{i}^{\mu}\equiv\partial U_{i}$
$\displaystyle=-\Gamma_{M}^{U_{i}}\mu_{M}^{U_{i}}-\Gamma_{Y}^{U_{i}}\mu_{Y}^{U_{i}}+\Gamma_{ss}\mu_{ss}+S_{U_{i}}\,,$
$\displaystyle\partial D_{i}$
$\displaystyle=-\Gamma_{M}^{D_{i}}\mu_{M}^{D_{i}}-\Gamma_{Y}^{D_{i}}\mu_{Y}^{D_{i}}+\Gamma_{ss}\mu_{ss}+S_{D_{i}}\,,$
$\displaystyle\partial Q_{i}$ $\displaystyle=-\partial U_{i}-\partial
D_{i}\,,$ $\displaystyle\partial E_{i}$
$\displaystyle=-\Gamma_{M}^{E_{i}}\mu_{M}^{E_{i}}-\Gamma_{Y}^{E_{i}}\mu_{Y}^{E_{i}}+S_{E_{i}}\,,$
$\displaystyle\partial L_{i}$ $\displaystyle=-\partial E_{i}\,,$
$\displaystyle\partial h$
$\displaystyle=\sum_{i=1}^{3}\Gamma_{Y}^{U_{i}}\mu_{Y}^{U_{i}}-\sum_{i=1}^{3}\Gamma_{Y}^{D_{i}}\mu_{Y}^{D_{i}}-\sum_{i=1}^{3}\Gamma_{Y}^{E_{i}}\mu_{Y}^{E_{i}}\,,$
(9)

where $U_{R},D_{R},Q_{L}$ are the SM quark fields, $E_{R},L_{L}$ are the SM
charged lepton fields (the chirality is implicit hereafter), $\Gamma_{M,Y}$
are the relaxation and Yukawa rates, respectively, and $S_{f}$ are the CPV
sources. The chemical potentials are given by

$\displaystyle\mu_{Y}^{U_{i}}$
$\displaystyle=\frac{{U_{i}}}{k_{U_{i}}}-\frac{{Q_{i}}}{k_{Q_{i}}}-\frac{h}{k_{h}}\,,$
$\displaystyle\mu_{M}^{U_{i}}$
$\displaystyle=\frac{{U_{i}}}{k_{U_{i}}}-\frac{{Q_{i}}}{k_{Q_{i}}}\,,$
$\displaystyle\mu_{Y}^{D_{i}}$
$\displaystyle=\frac{{D_{i}}}{k_{D_{i}}}-\frac{Q_{i}}{k_{Q_{i}}}+\frac{h}{k_{h}}\,,$
$\displaystyle\mu_{M}^{D_{i}}$
$\displaystyle=\frac{{D_{i}}}{k_{D_{i}}}-\frac{{Q_{i}}}{k_{Q_{i}}}\,,$
$\displaystyle\mu_{Y}^{E_{i}}$
$\displaystyle=\frac{{E_{i}}}{k_{E_{i}}}-\frac{{L_{i}}}{k_{L_{i}}}+\frac{h}{k_{h}}\,,$
$\displaystyle\mu_{M}^{E_{i}}$
$\displaystyle=\frac{{E_{i}}}{k_{E_{i}}}-\frac{{L_{i}}}{k_{L_{i}}}\,,$
$\displaystyle\mu_{ss}$
$\displaystyle=\sum_{i=1}^{3}\left(\frac{2Q_{i}}{k_{Q_{i}}}-\frac{U_{i}}{k_{U_{i}}}-\frac{D_{i}}{k_{D_{i}}}\right)\,,$
(10)

where $k$ counts the finite temperature degrees of freedom. All of the
benchmark parameters were taken from Ref. Fuchs:2020pun , and presented in
Appendix A. The solution of the transport equations is obtained by reduction
of order Fuchs:2020pun , i.e. $N$\- second order differential equations are
written as a set of $2N$\- first order differential equations, and numerical
diagonalization. By the second step of the approach the left-handed particles
participate in the sphaleron process which generates a baryon asymmetry (Eqn.
(58)).
The two-step procedure is summarized in more detail in Appendix B.

## Numerical result

Our numerical calculation for the BAU yields

$\displaystyle Y_{B}\left(\kappa_{f}^{I}\right)=-Y_{B}^{\text{obs}}\cdot$
$\displaystyle\left(-28\kappa_{t}^{I}+11\kappa_{\tau}^{I}+0.2\kappa_{b}^{I}+0.1\kappa_{\mu}^{I}\right.$
$\displaystyle\ \left.+0.03\kappa_{c}^{I}+2\cdot 10^{-4}\kappa_{s}^{I}+3\cdot
10^{-6}\kappa_{e}^{I}\right.$ $\displaystyle\ \left.+4\cdot
10^{-7}\kappa_{d}^{I}+9\cdot 10^{-8}\kappa_{u}^{I}\right)\,.$ (11)

This result agrees with previous analysis for the third generation particles
Fuchs:2020uoc  and the muon Fuchs:2019ore , but also includes the rest of the
charged SM fermions.
We point out that although quarks have (on average) stronger dim-$4$ Yukawa
couplings, which positively impact the CPV source, they also have more
washout. It is the result of lower diffusion and higher interaction rates, as
well as an additional interaction via the strong sphaleron. This is an
important difference between charged leptons and quarks.
Another interesting feature that holds only for the light fermions is that the
ratio between contributions of different flavors of either leptons or quarks
to the BAU is proportional to the SM Yukawa couplings squared, up to $5\%$.
Thus, for $f,f^{\prime}$ light flavors of either charged leptons or quarks

$\displaystyle\frac{Y_{B}^{f}}{Y_{B}^{f^{\prime}}}$
$\displaystyle\approx\left(\frac{y_{f}}{y_{f^{\prime}}}\right)^{2}\frac{\kappa_{f}^{I}}{\kappa_{f^{\prime}}^{I}}\,.$
(12)

This relation does not hold for the third generation particles, because their
interaction rates are orders of magnitude higher.
Finally, the desired $\kappa_{f}^{I}$ that saturates the contribution to the
BAU to its observed value could correspond to an unfavorable solution, when
demanding the theory to be perturbative. Let us denote the $\kappa_{f}^{I}$
which satisfies $Y_{B}=Y_{B}^{\text{obs}}$, according to Eqn. (Numerical
result), as $\kappa_{f}^{I*}$. By rearranging the definition of
$\kappa_{f}^{I}$ (Eqn. ( framework)), the solution is described by:

$\displaystyle\left(T_{f}^{R}+1\right)^{2}+\left(T_{f}^{I}-\frac{1}{\kappa_{f}^{I*}}\right)^{2}=\left(\frac{1}{\kappa_{f}^{I*}}\right)^{2}\,,$
(13)

See Fig. 2 for the geometrical interpretation.

Figure 2: The geometrical description of Eqn. (13) is a circle in
$\left(T_{f}^{R},T_{f}^{I}\right)$-space centered around $\left(-1,\
\frac{1}{\kappa_{f}^{I*}}\right)$ with radius
$\frac{1}{\left|\kappa_{f}^{I*}\right|}$.

For large values of $\kappa_{f}^{I*}$, this setting nearly cancels the dim-$4$
contribution, such that the mass of the fermions is effectively generated by
the complex part of the dim-$6$ term, i.e. it requires fine-tuning. Moreover,
for some cases, depending on the particle and the position on the circle, it
requires $y_{f}>4\pi$ which in non-perturbative (See Fig. 3), rendering the
analysis moot. We specifically point out the up and down quarks, for which
there is no perturbative theory that can account for the observed BAU.

Figure 3: Dim-$4$ Yukawa coupling $y_{f}$ as a function of $\theta$, the angle
on the circle described in Eqn. (13). The region above the dashed line
$y_{f}=4\pi$ is non-perturbative and is therefore considered unfavorable. Note
that $\theta=\frac{\pi}{2}$, which corresponds to $\left(-1,0\right)$ in
$\left(T_{f}^{R},T_{f}^{I}\right)$-space, is clearly unphysical and should be
excluded for all flavors. The up and down quarks are non-perturbative
throughout the entire range, which precludes them as sole-source of the
observed BAU.

## Bounds

The implications of a non-zero $\kappa_{f}^{I}$ are threefold; In addition to
the generation of baryon asymmetry, which was discussed above,
$\kappa_{f}^{I}$ contributes to the EDMs of the electron and the neutron, and
Higgs related measurements. In this Section we use these experimental results
to constrain the maximal value of $\kappa_{f}^{I}$. This will allow us to
infer the maximal contribution, from a single flavor modification, to the BAU.

##  electron-Electric Dipole Moment (eEDM) Bound

The upper bound obtained by the ACME collaboration on the eEDM Andreev:2018ayy
is

$\displaystyle\left|d_{e}^{\text{max}}\right|$ $\displaystyle=1.1\cdot
10^{-29}\ e\ \text{cm at $90\%$\leavevmode\nobreak\ C.L.}\,.$ (14)

The contribution of the SM fermions, other than the top and the electron, to
the eEDM is given by Brod:2018pli ; Panico:2018hal

$\displaystyle\frac{d_{e}}{e}$ $\displaystyle\simeq
4N_{c}Q_{f}^{2}\frac{\alpha}{\left(4\pi\right)^{3}}\frac{m_{e}m_{f}^{2}}{v^{2}m_{h}^{2}}\left(\ln^{2}\left(\frac{m_{f}^{2}}{m_{h}^{2}}\right)+\frac{\pi^{2}}{3}\right)\kappa_{f}^{I}\,.$
(15)

The contribution of the top quark to the eEDM is given by Brod:2013cka

$\displaystyle\frac{d_{e}}{e}$ $\displaystyle\simeq 9.4\cdot
10^{-27}\kappa_{t}^{I}\ \text{cm}\,.$ (16)

The contribution of the electron, via two-loop diagrams, to the eEDM is given
by Altmannshofer:2015qra

$\displaystyle\frac{d_{e}}{e}$ $\displaystyle\simeq 5.1\cdot
10^{-27}\kappa_{e}^{I}\ \text{cm}\,.$ (17)

##  Neutron-Electric Dipole Moment (nEDM) bound

The upper bound obtained by the nEDM collaboration on the nEDM Abel:2020gbr
is

$\displaystyle\left|d_{n}^{\text{max}}\right|$ $\displaystyle=1.8\cdot
10^{-26}\ e\ \text{cm at $90\%$\leavevmode\nobreak\ C.L.}\,.$ (18)

This bound constrains only the $\kappa_{f}^{I}$ of the quarks. The
implications of the nEDM bound on $\kappa^{I}_{u,d,s}$ were calculated in
Brod:2018lbf . Here we show the bound, updated to the latest measurement:

$\displaystyle\left|\kappa_{u}^{I}\right|$ $\displaystyle\lesssim 0.6\,,$
$\displaystyle\left|\kappa_{d}^{I}\right|$ $\displaystyle\lesssim 0.14\,,$
$\displaystyle\left|\kappa_{s}^{I}\right|$ $\displaystyle\lesssim 4.5\,.$ (19)

We also update the bounds on $\kappa_{b,c}^{I}$ given in Brod:2018pli ,
assuming $\real(\kappa_{b,c})=0$,

$\displaystyle\left|\kappa_{b}^{I}\right|$ $\displaystyle\lesssim\ 7.4\,,$
$\displaystyle\left|\kappa_{c}^{I}\right|$ $\displaystyle\lesssim 13.6\,.$
(20)

## Higgs boson decay

Assuming the production rate $\left(\sigma\right)$ of the Higgs boson is
unchanged due to the dim-$6$ operator, the signal strength $\mu_{h\to
f^{+}f^{-}}$ can be written as

$\displaystyle\mu_{i,h\to f^{+}f^{-}}$
$\displaystyle\equiv\frac{\sigma_{i}\left(pp\to h\right)\mathcal{B}\left(h\to
f^{+}f^{-}\right)}{\sigma_{i}^{\text{SM}}\left(pp\to
h\right)\mathcal{B}^{\text{SM}}\left(h\to f^{+}f^{-}\right)}$
$\displaystyle\approx\frac{\mathcal{B}\left(h\to
f^{+}f^{-}\right)}{\mathcal{B}^{\text{SM}}\left(h\to f^{+}f^{-}\right)}\,,$
(21)

where $\mathcal{B}$ is the branching ratio. It is then straightforward to
bound $\kappa_{f}$ using the signal strength. For the charged leptons, we
translate the upper bound on $\mu_{h\to\ell\ell}$ to the maximal value of
$|\kappa_{\ell}|$ via

$\displaystyle\mu_{i,h\to\ell^{+}\ell^{-}}$
$\displaystyle=\frac{\left|\kappa_{\ell}\right|^{2}}{1+\left(\left|\kappa_{\ell}\right|^{2}-1\right)\mathcal{B}\left(h\to\ell^{+}\ell^{-}\right)^{\text{SM}}}\,.$
(22)

For the light quarks, we bound $|\kappa_{q}|$ via its effect on the total
decay width of the Higgs boson. For that we used $\mu_{h\to b\bar{b}}$ which
was measured to be close to unity Sirunyan:2018kst . When NP interacts with
$q=u,c,d,s$, the signal strength of $h\to b\bar{b}$ is modified as

$\displaystyle\mu_{h\to b\bar{b}}$
$\displaystyle=\frac{1}{1+\left(\left|\kappa_{q}\right|^{2}-1\right)\mathcal{B}\left(h\to
q\bar{q}\right)^{\text{SM}}}\,.$ (23)

We can then use the experimental bound on $\mu_{h\to b\bar{b}}$ to constrain
the maximal branching fraction as follows:

$\displaystyle\mathcal{B}\left(h\to q\bar{q}\right)^{\text{max}}$
$\displaystyle\lesssim\frac{1}{\mu^{\text{min}}_{h\to b\bar{b}}}-1\approx
0.41\,,$ (24)

where we used $\mu_{h\to b\bar{b}}^{\text{min}}\approx 0.71$ at $95\%$
confidence level (C.L.). Assuming the only modification to the SM is the
dim-$6$ term $X$, we obtain

$\displaystyle\lambda_{q}\leq\sqrt{\frac{\mathcal{B}\left(h\to
q\bar{q}\right)^{\text{max}}}{\mathcal{B}\left(h\to
b\bar{b}\right)^{\text{min}}}}\
\frac{y_{b}}{\sqrt{2}}=\sqrt{\frac{\mathcal{B}\left(h\to
q\bar{q}\right)^{\text{max}}}{\mu_{h\to
b\bar{b}}^{\text{min}}\mathcal{B}\left(h\to b\bar{b}\right)^{\text{SM}}}}\
\frac{y_{b}}{\sqrt{2}}\,.$ (25)

Using the SM prediction $\mathcal{B}\left(h\to
b\bar{b}\right)^{\text{SM}}=0.58$, we get, up to $<1\%$

$\displaystyle\lambda_{q}\lesssim\frac{y_{b}}{\sqrt{2}}\,.$ (26)

Then, using eqn. ( framework),

$\displaystyle\left|\kappa_{q}\right|\lesssim\frac{m_{b}}{m_{q}}\,.$ (27)

Table 1:  The BAU calculated following the full set of transport equations. $Y_{B}^{f}$ is the BAU resulting from $S_{f}$. Collider constraints are at $\sim 95\%$ C.L. Fuchs:2020uoc ; Aad:2020xfq ; ICHEP:2020 ; Sirunyan:2018kst ; Aad:2019ojw . EDM constraints are at $90\%$ C.L. Andreev:2018ayy ; Abel:2020gbr  for all, except for the bottom and charm, for which the nEDM constraints are at $68\%$ C.L. Brod:2013cka . Source |  BAU, $Y_{B}^{f}$ | Higgs decay | eEDM | nEDM | maximal % $Y_{B}^{\text{obs}}$
---|---|---|---|---|---
$S_{\tau}$ | $-9.9\cdot 10^{-10}\cdot\kappa^{I}_{\tau}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{\tau}\right|\lesssim 1.1}$ | $\left|\kappa_{\tau}^{I}\right|\lesssim 0.3$ | - | $337\%$
$S_{\mu}$ | $-1.0\cdot 10^{-11}\cdot\kappa^{I}_{\mu}$ | $\left|\kappa_{\mu}\right|\lesssim 1.4$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{\mu}^{I}\right|\lesssim 31}$ | - | $17\%$
$S_{b}$ | $-2.1\cdot 10^{-11}\cdot\kappa^{I}_{b}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{b}\right|\lesssim 1.7}$ | $\left|\kappa_{b}^{I}\right|\lesssim 0.2$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{b}^{I}\right|\lesssim\ 7.4}$ | $5.8\%$
$S_{t}$ | $2.4\cdot 10^{-9}\cdot\kappa^{I}_{t}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{t}\right|\lesssim 1.2}$ | $\left|\kappa_{t}^{I}\right|\lesssim 1.2\cdot 10^{-3}$ | - | $3.3\%$
$S_{c}$ | $-2.7\cdot 10^{-12}\cdot\kappa^{I}_{c}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{c}\right|\lesssim 3.3}$ | $\left|\kappa_{c}^{I}\right|\lesssim 0.4$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{c}^{I}\right|\lesssim\ 13.6}$ | $1.1\%$
$S_{s}$ | $-1.6\cdot 10^{-14}\cdot\kappa^{I}_{s}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{s}\right|\lesssim 43}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{s}^{I}\right|\lesssim 109}$ | $\left|\kappa_{s}^{I}\right|\lesssim 4.5$ | $0.08\%$
$S_{d}$ | $-3.8\cdot 10^{-17}\cdot\kappa^{I}_{d}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{d}\right|\lesssim 890}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{d}^{I}\right|\lesssim 2.3\cdot 10^{4}}$ | $\left|\kappa_{d}^{I}\right|\lesssim 0.14$ | $6\cdot 10^{-6}\%$
$S_{u}$ | $-8.1\cdot 10^{-18}\cdot\kappa^{I}_{u}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{u}\right|\lesssim 1900}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{u}^{I}\right|\lesssim 2.2\cdot 10^{4}}$ | $\left|\kappa_{u}^{I}\right|\lesssim 0.6$ | $6\cdot 10^{-6}\%$
$S_{e}$ | $-2.5\cdot 10^{-16}\cdot\kappa^{I}_{e}$ | $\mathbin{\color[rgb]{.5,.5,.5}\left|\kappa_{e}\right|\lesssim 266}$ | $\left|\kappa_{e}^{I}\right|\lesssim 2.2\cdot 10^{-3}$ | - | $6\cdot 10^{-7}\%$

## Results

We present our results in Table 1. The first prominent result is that no light
charged fermion could give the dominant contribution to the BAU. Of all the SM
charged fermions only the $\tau$ could produce $100\%$ of the observed BAU
Fuchs:2020uoc . The next in importance can be the $\mu$ Fuchs:2019ore , which
brings us to consider the relatively negligible effect of the quarks. As
mentioned, although the CPV term is proportional to the dim-$4$ Yukawa
coupling squared, which tends to be larger for quarks, they also have more
washout than leptons which leads to lower contribution to the BAU. In
addition, the bounds are comparable, and thus their maximal percentage is
relatively small.
That being said, one could consider two flavor modification, in which case the
electron has a special feature. Assuming the interplay of different species
does not dramatically change the numerical calculation, the electron could
cancel the contribution of another particle to the eEDM, while leaving the
contribution to the BAU essentially unchanged. When combined, the interference
allows particles that are constrained mostly by the eEDM, such as the top
Fuchs:2020uoc , to account for the BAU, e.g. for $\kappa_{t}^{I}\approx 0.04$
the top generates $Y_{B}^{\text{obs}}$, while $\kappa_{e}^{I}\approx-0.06$
cancels the top’s contribution to the eEDM (See Fig. 4).

Figure 4: The interplay between the electron and the top could relax the
stringent bound on the top, which is from the eEDM measurement. This allows
the top to saturate the BAU. Since the next leading bound on both particles is
orders of magnitude weaker, such cancellation would be difficult to detect.

Because of possible cancellation, it is much harder to exclude such an elusive
hypothetical scenario.

## Conclusions

We consider the CPV source resulting from a dim-$6$ SMEFT term which couples
three Higgs fields to the SM fermions. We generalize the procedure described
in Ref. Fuchs:2020pun  to evaluate the complete set of single flavor
modifications from all of the SM charged fermions.
We deduce that although a stronger dim-$4$ Yukawa coupling enhances the CPV
source, quarks have more washout than leptons and are therefore less favorable
candidates to produce the BAU. Moreover, to saturate $Y_{B}^{\text{obs}}$,
some of the particles require non-perturbative dim-$4$ Yukawa couplings, e.g.
the up and down quarks, and are therefore unequivocally ruled out as sole-
contributors via this mechanism.
Constrained by the upper bounds of the electron and the neutron EDMs and the
Higgs boson decay measurements we evaluate the maximal contribution from each
single flavor modification (See Table 1). We conclude that the $\tau$ is the
only candidate able to produce $100\%\ Y_{B}^{\text{obs}}$ Fuchs:2020uoc .
Other than the $\mu$, which could provide up to $17\%\ Y_{B}^{\text{obs}}$
Fuchs:2019ore , the rest of the charged fermions produce negligible
contributions (less than $6\%\ Y_{B}^{\text{obs}}$).
An interplay of different flavors could relax current bounds to allow the
observed baryon asymmetry be accounted for by a CPV Yukawa coupling.
Specifically, the interplay of the top and the electron could allow the top to
saturate $Y_{B}^{\text{obs}}$. In addition, off-diagonal couplings should also
be considered in a more complete description.

## Acknowledgments

We are grateful to Yossi Nir and Yehonatan Viernik for helpful discussions.

##  Appendix A Benchmark parameters

The input used for this work is the following:

  * •

Coupling constant at nucleation temperature:

$\displaystyle g_{s}$ $\displaystyle=1.23\,,$ $\displaystyle g$
$\displaystyle=0.65\,,$ $\displaystyle g^{\prime}$ $\displaystyle=0.36\,.$
(28)

  * •

Bubble wall velocity and width:

$\displaystyle v_{w}$ $\displaystyle=0.05\,,$ $\displaystyle L_{w}$
$\displaystyle=0.11\ \text{GeV}^{-1}\,.$ (29)

  * •

VEV during nucleation and at $0$ temperature:

$\displaystyle v_{N}$ $\displaystyle=152\text{ GeV}\,,$ $\displaystyle v_{0}$
$\displaystyle=246\text{ GeV}\,.$ (30)

  * •

The SM fermion masses were taken from Tanabashi:2018oca .

  * •

The diffusion coefficients are calculated in Joyce:1994zn :

$\displaystyle D_{\ell_{L}}$ $\displaystyle=D_{h}=\frac{100}{T}\,,$
$\displaystyle D_{\ell_{R}}$ $\displaystyle=\frac{380}{T}\,,$ $\displaystyle
D_{q_{L}}$ $\displaystyle=D_{q_{R}}=\frac{6}{T}\,.$ (31)

  * •

The Mass and Yukawa rates are given in Table 2.

Table 2: Relaxation rates (for the broken phase), and Yukawa rate (for both phases), calculated from Fuchs:2020pun . Particle |  $\Gamma_{M}^{B}$ (GeV) |  $\Gamma_{Y}$ (GeV)
---|---|---
$\tau$ | $4.9\cdot 10^{-3}$ | $5.6\cdot 10^{-4}$
$\mu$ | $1.7\cdot 10^{-5}$ | $2.0\cdot 10^{-6}$
$e$ | $3.9\cdot 10^{-10}$ | $4.4\cdot 10^{-11}$
$t$ | $102$ | $2.6$
$c$ | $4.7\cdot 10^{-3}$ | $1.6\cdot 10^{-4}$
$u$ | $1.4\cdot 10^{-8}$ | $4.7\cdot 10^{-10}$
$b$ | $5.3\cdot 10^{-2}$ | $1.7\cdot 10^{-3}$
$s$ | $2.7\cdot 10^{-5}$ | $9.0\cdot 10^{-7}$
$d$ | $6.5\cdot 10^{-8}$ | $2.1\cdot 10^{-9}$

These interaction rates are calculated for $T_{I}=T_{R}=0$.

  * •

The weak sphaleron rate

$\displaystyle\Gamma_{ws}$ $\displaystyle=6\underbrace{\kappa}_{\sim
20}\underbrace{\alpha_{w}}_{\alpha_{w}=\frac{g^{2}}{4\pi}}^{5}T=120T\left(\frac{g^{2}}{4\pi}\right)^{5}$
$\displaystyle\overset{\text{nucleation}}{\longrightarrow}120T_{N}\left(\frac{0.65^{2}}{4\pi}\right)^{5}\approx
4.5\cdot 10^{-4}\ \text{GeV}\,.$ (32)

  * •

The strong sphaleron rate

$\displaystyle\Gamma_{ss}$ $\displaystyle=14\alpha_{s}^{4}T\approx 0.26\
\text{GeV}\,.$ (33)

  * •

Thermal width:

$\displaystyle\Gamma_{\text{leptons}}$ $\displaystyle\approx 2\cdot 10^{-3}\
T\,,$ $\displaystyle\Gamma_{\text{quarks}}$ $\displaystyle\approx 0.16\ T\,.$
(34)

The temperature during nucleation is $T_{N}=88$ GeV, and for the SM we use
$\mathcal{R}=\frac{15}{4}$.

Thermal functions:

  * •

$J$ is given by

$\displaystyle J_{e}(T)$
$\displaystyle=\int\limits_{0}^{\infty}\frac{k^{2}dk}{\omega_{L}^{e}\omega_{R}^{e}}\imaginary\left[\frac{n_{e}\left(\mathcal{E}_{L}^{e}\right)-n_{e}\left(\mathcal{E}_{R}^{e*}\right)}{\left(\mathcal{E}_{L}^{e}-\mathcal{E}_{R}^{e*}\right)^{2}}\left(\mathcal{E}_{L}^{e}\mathcal{E}_{R}^{e*}-k^{2}\right)\right.$
$\displaystyle\left.+\frac{n_{e}\left(\mathcal{E}_{L}^{e}\right)+n_{e}\left(\mathcal{E}_{R}^{e}\right)}{\left(\mathcal{E}_{L}^{e}+\mathcal{E}_{R}^{e}\right)^{2}}\left(\mathcal{E}_{L}^{e}\mathcal{E}_{R}^{e}+k^{2}\right)\right]\,.$
(35)

  * •

The source is given by

$\displaystyle S_{f}(z;T)$
$\displaystyle=\frac{v_{w}N_{c}^{f}}{\pi^{2}}\imaginary\left(m_{f}^{\prime}m_{f}^{*}\right)J_{f}(T)$
$\displaystyle=\frac{v_{w}N_{c}^{f}{y^{f}_{\text{SM}}}^{2}}{2\pi^{2}v_{0}^{2}}J_{f}(T)\phi_{b}^{3}(z)\phi^{\prime}_{b}(z)\times\kappa_{I}^{f}\,.$
(36)

For the background Higgs boson field we use the kink solution:

$\displaystyle\phi_{b}$
$\displaystyle=\frac{v_{N}}{2}\left(1+\tanh(\frac{z}{L_{w}})\right)\,.$ (37)

  * •

The frequencies, energies and Fermi- Dirac distributions are

$\displaystyle\omega_{L,R}^{f_{i}}(k)$ $\displaystyle=\sqrt{k^{2}+\real(\delta
m_{f_{i\ L,R}}^{2}(T))}\,,$ $\displaystyle\mathcal{E}_{L}^{f_{i}}$
$\displaystyle=\omega_{L,R}^{f_{i}}(k)-i\Gamma_{f_{i}}\,,$ $\displaystyle
n_{f}\left(\mathcal{E}\right)$
$\displaystyle=\frac{1}{e^{\frac{\mathcal{E}}{T}}+1}$ (38)

  * •

The thermal masses are given by

$\displaystyle\real(\delta m_{L_{Li}}^{2}(T))$
$\displaystyle=\left(\frac{3}{32}g^{2}+\frac{1}{32}{g^{\prime}}^{2}+\frac{1}{16}y_{e_{i}}^{2}\right)T^{2}\equiv
a_{L_{Li}}^{2}T^{2}\,,$ $\displaystyle\real(\delta m_{e_{iR}}^{2}(T))$
$\displaystyle=\left(\frac{1}{8}{g^{\prime}}^{2}+\frac{1}{8}y_{e_{i}}^{2}\right)T^{2}\equiv
a_{e_{Ri}}^{2}T^{2}\,,$ $\displaystyle\real(\delta m_{q_{Li}}^{2}(T))$
$\displaystyle=\left(\frac{1}{6}g_{s}^{2}+\frac{3}{32}g^{2}+\frac{1}{288}{g^{\prime}}^{2}+\frac{1}{16}y_{u_{i}}^{2}\right.$
$\displaystyle\left.+\frac{1}{16}y_{d_{i}}^{2}\right)T^{2}\equiv
a_{q_{Li}}^{2}T^{2}\,,$ $\displaystyle\real(\delta m_{u_{iR}}^{2}(T))$
$\displaystyle=\left(\frac{1}{6}g_{s}^{2}+\frac{1}{18}{g^{\prime}}^{2}+\frac{1}{8}y_{u_{i}}^{2}\right)T^{2}\equiv
a_{u_{Ri}}^{2}T^{2}\,,$ $\displaystyle\real(\delta m_{d_{iR}}^{2}(T))$
$\displaystyle=\left(\frac{1}{6}g_{s}^{2}+\frac{1}{72}{g^{\prime}}^{2}+\frac{1}{8}y_{d_{i}}^{2}\right)T^{2}\equiv
a_{d_{Ri}}^{2}T^{2}\,,$ $\displaystyle\real(\delta m_{h}^{2}(T))$
$\displaystyle=\left(\frac{3}{16}g^{2}+\frac{1}{16}{g^{\prime}}^{2}+\frac{1}{12}\sum_{i=e,\mu,\tau}y_{e_{i}}^{2}\right.$
$\displaystyle\left.+\frac{1}{4}\sum_{i=u,c,t}y_{u_{i}}^{2}+\frac{1}{4}\sum_{i=d,s,b}y_{d_{i}}^{2}\right)T^{2}\equiv
a_{h}^{2}T^{2}\,.$ (39)

  * •

The finite temperature degrees of freedom are given by

$\displaystyle k_{f_{i}}(a_{f_{i}})$
$\displaystyle=k_{0}^{f_{i}}\frac{6}{\pi^{2}}\int\limits_{a_{f_{i}}}^{\infty}dx\frac{xe^{x}}{\left(e^{x}\pm
1\right)^{2}}\sqrt{x^{2}-a_{f_{i}}^{2}}\,.$ (40)

where $k_{0}^{f}$ is the number of degrees of freedom, and $+$ ($-$) is for
fermions (the Higgs boson).

##  Appendix B Solving the set of transport equations: the two-step approach

In this appendix we provide further information regarding the analytic
techniques we used to calculate the produced baryon asymmetry, similarly to
Ref. Fuchs:2020pun . We used the two-step approach: First we solve the set of
transport equations, as given in Eqn. (Two-step approach). Then, the baryon
asymmetry is obtained by summing over the left-handed number densities and
solving an equation which describes the weak sphaleron process.

## Reduction of dimensions

The left hand side of Eqn. (Two-step approach) can be written as a one
dimensional - second order differential equation with respect to the bubble
wall dimension denoted $z$:

$\displaystyle\partial f\equiv\partial_{\mu}f^{\mu}$
$\displaystyle=\frac{\partial f^{0}}{\partial
t}-\vec{\nabla}\cdot\vec{f}=\underbrace{\frac{\partial z}{\partial t}}_{\equiv
v_{w}}\underbrace{\frac{\partial f^{0}}{\partial
z}}_{f^{\prime}}+\vec{\nabla}\cdot\left(-D_{f}\vec{\nabla}f^{0}\right)$
$\displaystyle=v_{w}f^{\prime}-D_{f}\underbrace{\laplacian{f^{0}}}_{\equiv
f^{\prime\prime}}=v_{w}f^{\prime}-D_{f}f^{\prime\prime}\,,$ (41)

where we used Fick’s first law and the diffusion approximation.

## Reduction of order

We can solve this set of $N=16$\- second order differential equations by
reduction of order:

$\displaystyle g_{f_{i}}$ $\displaystyle\equiv f^{\prime}_{i}\,,$
$\displaystyle\vec{f}$
$\displaystyle=\begin{pmatrix}\vec{U}_{R}&\vec{D}_{R}&\vec{Q}_{L}&\vec{E}_{R}&\vec{L}_{L}&h\end{pmatrix}\,,$
$\displaystyle\vec{\chi}$
$\displaystyle=\begin{pmatrix}\vec{f}&\vec{g}_{f}\end{pmatrix}^{T}\,.$
$\displaystyle\vec{\chi}^{\prime}=\begin{pmatrix}0&\mathbbm{1}_{NxN}\\\
\hat{\Gamma}&\hat{V}\end{pmatrix}\vec{\chi}+\vec{S}\equiv\hat{K}\vec{\chi}+\vec{S}\,.$
(42)

We are left with $2N$-first order differential equations. Then, $\hat{K}$ can
be diagonalized numerically using MATLAB.

## Numerical diagonalization

The solution for the symmetric phase is given by

$\displaystyle{\vec{\chi}}^{S}$
$\displaystyle=\sum_{i=1}^{2N}C_{i}^{S}e^{\lambda_{i}^{S}z}\vec{u}_{i}^{S}\equiv\hat{\Phi}^{S}(z)\vec{C}^{S}\,,$
$\displaystyle\hat{\Phi}^{X}_{i,j}\left(z\right)$
$\displaystyle=e^{\lambda_{jj}^{X}z}\left(\vec{u}_{j}^{X}\right)_{i}\,.$ (43)

where $C_{j}$’s are constants, $\lambda_{j}$ are the eigenvalues, and
$\vec{u}_{j}$’s are the eigenvectors. We define

$\displaystyle\hat{\lambda}$
$\displaystyle=\text{diag}\left(\lambda_{i}\right)$ $\displaystyle i$
$\displaystyle=[1:2N]\,.$ (44) $\displaystyle\hat{\phi}$
$\displaystyle=\begin{pmatrix}|&|&&|\\\
\vec{u}_{1}&\vec{u}_{2}&\dots&\vec{u}_{2N}\\\ |&|&&|\end{pmatrix}$ (45)

Accordingly,

$\displaystyle\hat{\Phi}_{i,j}\left(z\right)$
$\displaystyle=\hat{\phi}_{ij}e^{\hat{\lambda}_{jj}z}=\begin{pmatrix}|&|&&|\\\
e^{\lambda_{1}z}\vec{u}_{1}&e^{\lambda_{2}z}\vec{u}_{2}&\dots&e^{\lambda_{2N}z}\vec{u}_{2N}\\\
|&|&&|\end{pmatrix}$ (46)

The full solution in the broken phase is obtained by variation of parameters
to be

$\displaystyle{\vec{\chi}}^{B}$
$\displaystyle=\hat{\Phi}^{B}(z)\vec{C}^{B}+\hat{\Phi}^{B}(z)\int\limits_{0}^{z}\left(\hat{\Phi}^{B}(x)\right)^{-1}\vec{S}(x)\
dx\,.$ (47)

## Boundary conditions

  * •

The integration constants of the divergent modes in the symmetric phase (which
correspond to $\lambda^{S}_{j}\leq 0$) are set to zero,

$\displaystyle C^{S}_{0-}$ $\displaystyle=0\,.$ (48)

  * •

The positive eigenvalues in the broken phase, $C^{B+}_{j}$ (which correspond
to $\lambda^{B}_{j}>0$), are chosen such that they cancel the divergent part
of the full solution at infinity:

$\displaystyle\vec{C}^{B}_{+}$
$\displaystyle=-\int\limits_{0}^{\infty}\left(\hat{\Phi}^{B}_{+}(x)\right)^{-1}\vec{S}(x)\
dx\,.$ (49)

  * •

We demand continuity at $z=0$.

    * –

In the symmetric phase we have

$\displaystyle\vec{\chi}^{S}_{i}\left(z\to 0^{-}\right)$
$\displaystyle=\hat{\phi}^{S}_{ij}\vec{C}^{S}_{j}=\hat{\phi}^{S}_{i+}\vec{C}^{S}_{+}\,.$
(50)

    * –

In the broken phase we have

$\displaystyle\vec{\chi}^{B}_{i}\left(z\to 0^{+}\right)$
$\displaystyle=\hat{\phi}^{B}_{ij}\vec{C}^{B}_{j}=\hat{\phi}^{B}_{i(0-)}\vec{C}^{B}_{0-}+\underbrace{\hat{\phi}^{B}_{i+}\vec{C}^{B}_{+}}_{\equiv
b_{i}}$ $\displaystyle=\hat{\phi}^{B}_{i(0-)}\vec{C}^{B}_{0-}+b_{i}\,.$ (51)

Continuity at $z=0$ is then

$\displaystyle\hat{\phi}^{S}_{i+}\vec{C}^{S}_{+}$
$\displaystyle\overset{!}{=}\hat{\phi}^{B}_{i(0-)}\vec{C}^{B}_{0-}+b_{i}\longrightarrow\hat{\phi}^{S}_{i+}\vec{C}^{S}_{+}-\hat{\phi}^{B}_{i(0-)}\vec{C}^{B}_{0-}=b_{i}\,.$
(52)

We obtain a linear set of equations,

$\displaystyle\hat{\phi}_{SB}$
$\displaystyle\equiv\left(\begin{array}[]{c|c}|&|\\\
\hat{\phi}^{S}_{i+}&\hat{\phi}^{B}_{i(0-)}\\\ |&|\end{array}\right)\,,$
$\displaystyle\vec{C}_{SB}$
$\displaystyle\equiv\begin{pmatrix}\vec{C}^{S}_{+}\\\
-\vec{C}^{B}_{0-}\end{pmatrix}\,,$ (56) $\displaystyle\hat{\phi}_{SB}\
\vec{C}_{SB}\overset{!}{=}b_{i}\,.$ (57)

Solving it sets the rest of the coefficients.

## The solution

The BAU is then given by

$\displaystyle Y_{B}$
$\displaystyle=\frac{3\Gamma_{\text{ws}}}{2D_{q}\alpha_{+}s}\int\limits_{0}^{-\infty}e^{-\alpha_{-}x}n_{L}(x)\
dx\,,$ (58)

where $n_{L}$ is the density of left handed particles in the symmetric phase
(where the weak sphaleron process is efficient),

$\displaystyle n_{L}(z)$
$\displaystyle=\sum_{i=1}^{3}\left(Q_{Li}(z)+L_{Li}(z)\right)\,,$ (59)

and

$\displaystyle\alpha_{\pm}$
$\displaystyle=\frac{1}{2D_{q}}\left(v_{w}\pm\sqrt{4D_{q}\Gamma_{\text{ws}}\mathcal{R}+v_{w}^{2}}\right)$
(60)

## References

  * (1) M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98, no.3, 030001 (2018)
  * (2) A. D. Sakharov, Pisma Zh. Eksp. Teor. Fiz. 5, 32 (1967) [JETP Lett. 5, 24 (1967)] [Sov. Phys. Usp. 34, no. 5, 392 (1991)] [Usp. Fiz. Nauk 161, no. 5, 61 (1991)].
  * (3) E. Fuchs, M. Losada, Y. Nir and Y. Viernik, JHEP 05, 056 (2020) [arXiv:2003.00099 [hep-ph]].
  * (4) E. Fuchs, M. Losada, Y. Nir and Y. Viernik, Phys. Rev. Lett. 124, no.18, 181801 (2020) [arXiv:1911.08495 [hep-ph]].
  * (5) J. De Vries, M. Postma and J. van de Vis, JHEP 04, 024 (2019) [arXiv:1811.11104 [hep-ph]].
  * (6) E. Fuchs, M. Losada, Y. Nir and Y. Viernik, [arXiv:2007.06940 [hep-ph]].
  * (7) V. Andreev et al. [ACME], Nature 562, no.7727, 355-360 (2018)
  * (8) G. Panico, A. Pomarol and M. Riembau, JHEP 04, 090 (2019) [arXiv:1810.09413 [hep-ph]].
  * (9) J. Brod and E. Stamou, [arXiv:1810.12303 [hep-ph]].
  * (10) J. Brod, U. Haisch and J. Zupan, JHEP 11, 180 (2013) doi:10.1007/JHEP11(2013)180 [arXiv:1310.1385 [hep-ph]].
  * (11) W. Altmannshofer, J. Brod and M. Schmaltz, JHEP 05, 125 (2015) [arXiv:1503.04830 [hep-ph]].
  * (12) C. Abel et al. [nEDM], Phys. Rev. Lett. 124, no.8, 081803 (2020) [arXiv:2001.11966 [hep-ex]].
  * (13) J. Brod and D. Skodras, JHEP 01, 233 (2019) [arXiv:1811.05480 [hep-ph]].
  * (14) A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 121, no.12, 121801 (2018) [arXiv:1808.08242 [hep-ex]].
  * (15) G. Aad et al. [ATLAS], [arXiv:2007.07830 [hep-ex]].
  * (16) R. Carlin [CMS], 40th International Conference on High Energy Physics (ICHEP2020), 3 Aug 2020, virtual conference.
  * (17) G. Aad et al. [ATLAS], Phys. Lett. B 801, 135148 (2020) [arXiv:1909.10235 [hep-ex]].
  * (18) M. Joyce, T. Prokopec and N. Turok, Phys. Rev. D 53, 2930-2957 (1996) [arXiv:hep-ph/9410281 [hep-ph]].

  *[BAU]: baryon asymmetry of the Universe
  *[SM]: Standard Model
  *[NP]: New Physics
  *[SMEFT]: Standard Model effective field theory
  *[CPV]: CP-violating
  *[dim-$6$]: dimension-six
  *[EDM]: Electric Dipole Moment
  *[dim-$4$]: dimension-four
  *[VEV]: vacuum expectation value
  *[]:
  *[eEDM]: electron-Electric Dipole Moment
  *[nEDM]: Neutron-Electric Dipole Moment
  *[C.L.]: confidence level


# Adaptive Epidemic Forecasting and Community Risk Evaluation of COVID-19

Vishrawas Gopalakrishnan, Sayali Navalekar, Pan Ding, Ryan Hooley, Jacob
Miller, Raman Srinivasan, Ajay Deshpande, Xuan Liu, Simone Bianco, James H.
Kaufman  IBM, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
(2018)

###### Abstract.

Pandemic control measures like lock-down, restrictions on restaurants and
gatherings, social-distancing have shown to be effective in curtailing the
spread of COVID-19. However, their sustained enforcement has negative economic
effects. To craft strategies and policies that reduce the hardship on the
people and the economy while being effective against the pandemic, authorities
need to understand the disease dynamics at the right geo-spatial granularity.
Considering factors like the hospitals’ ability to handle the fluctuating
demands, evaluating various reopening scenarios, and accurate forecasting of
cases are vital to decision making. Towards this end, we present a flexible
end-to-end solution that seamlessly integrates public health data with
tertiary client data to accurately estimate the risk of reopening a community.
At its core lies a state-of-the-art prediction model that auto-captures
changing trends in transmission and mobility. Benchmarking against various
published baselines confirm the superiority of our forecasting algorithm.
Combined with the ability to extend to multiple client-specific requirements
and perform deductive reasoning through counter-factual analysis, this
solution provides actionable insights to multiple client domains ranging from
government to educational institutions, hospitals, and commercial
establishments.

epidemiology, COVID-19, compartmental models

††copyright: acmcopyright††journalyear: 2018††doi:
10.1145/1122445.1122456††conference: Woodstock ’18: ACM Symposium on Neural
Gaze Detection; June 03–05, 2018; Woodstock, NY††booktitle: Woodstock ’18: ACM
Symposium on Neural Gaze Detection, June 03–05, 2018, Woodstock, NY††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Mathematics of computing Time
series analysis††ccs: Computing methodologies Supervised learning††ccs:
Computing methodologies Modeling and simulation††ccs: Applied computing
Bioinformatics

##  1\. Introduction

Due to the COVID-19 pandemic, it is estimated that half of the world’s 3.3
billion global workforces were/are at risk of losing their livelihoods (ILO
and WHO, 2020) and in the U.S itself, there has been a -11% change in the
number of low-wage category jobs, and overall a -6% decrease across
categories. Similar disruptions are reported across different work sectors
with AHA estimating losses of $202.6 billion for America’s hospitals and
health systems for the period of March-1 2020 to June-30 2020 (Meredith et
al., 2020) and over 1,300 colleges and universities across all 50 states
canceled in-person classes affecting the quality of education as well as
leading to financial losses (Smalley, 2020). To alleviate these issues,
studies have recommended that governments and industries should, as soon as
possible, try to resume operations even in a limited capacity and then adapt
their operations based on the evolving situation (Shen et al., 2020).

Figure 1. Overview Diagram of the Solution

Towards this goal, our solution - part of IBM’s Watson Works, provides timely
information on the pandemic’s current and future state in a region. As a
solution framework, with an objective to cater to different specific client
domains, it provides actionable insights that support planning safe reopening
of offices, crafting school schedules, meeting hospitalization/ICU demands,
planning for clinical vaccine trials, etc. Additionally, it has also helped
local officials and state government test hypotheses and crafting their
intervention strategies for the COVID-19 pandemic. Several efforts like (Li et
al., 2020; Team, 2020; Brooks et al., [n.d.]; Zou et al., 2020; COVID et al.,
2020) provide forecasts of COVID-19. However, the superiority of our method
lies in its ability to provide a solution to:

  * •

Model at a hyper-local level.

  * •

Translate the current and predicted case load numbers to community reopening
risk metrics in accordance with national/regional laws.

  * •

Easily adapt to business specific tasks and objectives like hospital/ICU
demand projections or identifying appropriate location for vaccine trials.

  * •

Scale globally through distributed processing.

  * •

Interact and simulate what-if scenarios.

  * •

Incorporate tertiary signals apart from cases and deaths data to better
estimate the future trends.

As shown in Section 4, our system consistently performs better in both short
and long-term projections, enabling users to take suitable corrective steps.
The “base” prediction model is easily enhanced to handle user-specific tasks
by plugging in custom analytics. For instance, this solution helped a major
hospital in south Florida to accurately estimate hospital bed and ICU demands;
thus, allowing them to plan/reschedule elective procedures and surgeries,
which otherwise could have got cancelled. It also helped a pharmaceutical
company to identify countries for their vaccine’s clinical trials.

Apart from customizing some of the known methodologies for data orchestration,
the deployed solution also makes a significant contribution to core
epidemiology by providing a template for models to incorporate important
tertiary signals without increasing the parameter search space. Further, it
also demonstrates a novel way of allowing parameters (e.g., transmission rate)
to assume multiple values over the course of the pandemic to better reflect
changes in the pandemic, including interventions, policy enactments, and their
time-varying efficacy. We also report learnings and insights regarding
variation in the relative importance of different input signals over time, the
impact of geographic level (county vs. state) on the prediction, etc.

##  2\. System and Deployment Overview

Figure 1 provides a holistic view of the data flow of our solution. The
solution begins by first connecting to ‘The Weather Channel’ (TWC) API and
consuming daily published COVID-19 cases, deaths, and other pertinent
statistics. TWC scrapes data from individual federal/state/county government
websites and stores data, including individual cases, test positivity rate,
hospitalization data, etc. After performing audits on data quality, TWC pushes
data into their database for downstream consumption.

The next step in the pipeline is Data Curation wherein raw COVID-19 cases and
death numbers are pre-processed and merged with tertiary data, including
mobility and other client-specific data such as individual hospital and ICU
demand data. The latter is optional as the Client Specific Analysis module can
be disabled. As shown in Figure 1, pre-processing steps include adjusting
spatial granularity, handling reporting irregularities, data smoothening, and
_inflection point detection_.

The pre-processed data is then passed to the Prediction Module which provides
case and death projections with upper and lower interval bounds. The current
case and death data, along with the projections from the model, are provided
to the Client Independent Analytics module to calculate various
epidemiological metrics and derive a score between 1-6 (lower is better)
reflecting the community transmission risk in a region. At the same time, the
model prediction output along with supplemental client data drives Client
Specific Analytics insights like hospitalization/ICU demand or simulating the
effect of an event (holiday travel, sports event) or intervention (government
restrictions) on the pandemic’s transmission rate.

The entire process is wrapped in an Apache NiFi data flow pipeline that
invokes relevant modules and manages the data flow. The pipeline is triggered
once a day after the latest COVID-19 data is populated by TWC. The ‘Prediction
Module’ is triggered once every three days as the projections rarely deviate a
lot in such a short period. However, the ‘Client Independent Analytics’, which
includes computing community risk, is triggered every day to provide up-to-
date statistics like doubling period, 14-day case trends, 1-week, and 2-week
ahead community risk levels, etc. The ‘Client Specific Analytics’ is run as
per client requirements.

The following sections focus on a pipeline configured for U.S.A data, although
the solution itself is applicable globally with relevant changes to account
for the availability of tertiary data.

##  3\. Methodology

This section details various algorithms and steps used in the different
modules shown in Figure 1 by first introducing the motivation/source of the
problem, the solution intuition, and lastly, by providing specific operational
examples.

###  3.1. Data Pre-Processing

This step entails two important tasks: (i) Denoising the COVID-19 data and
(ii) Identifying approximate timestamps of the inflection points in disease
dynamics.

####  3.1.1. Denoising Public Health Data:

##### Geo-Spatial Noise:

Any pandemic is defined by its transmission rate - the average rate at which
individuals pass on the infection. This, in turn, is determined by the
population contact or _mixing rate_ in the region. Assuming a similar contact
rate across regions will lead to an incorrect susceptible population. For
instance, counties of New York City-Newark-Jersey City, NY-NJ-PA Metropolitan
Statistical Area have a higher mixing rate than rural counties of Alaska or
North Dakota, which are sparsely populated and self-contained. Thus, the first
step is to normalize regions/counties by aggregating them into a self-
contained geo-unit of analysis using a metric that reflects the population
mixing characteristics. Additionally, it is often found that in an area where
there is a major hospital center, reported cases around that area are skewed
as patients from nearby counties are reported at the major hospital. For
example, COVID cases that occurred in Bronx County, N.Y. were reported as
originating in Westchester County, N.Y.

(a) Alaska

(b) New York

Figure 2. County Grouping by Residence-Work Commute

Towards this end, our solution uses Louvain’s method of graph clustering
(Blondel et al., 2008) to group adjacent counties together. The algorithm uses
the census bureau’s residence to workplace commute data
2011-2015111www.census.gov/data/tables/2015/demo/metro-micro/commuting-
flows-2015.html to create an adjacency matrix where each node is a county and
edges represent the strength of commute. Note that we symmetrize the adjacency
matrix to reflect the return journey as well. Using this approach, we were
able to group 3,132 counties into super-groups of 685 county clusters. There
are about 4.59 counties per cluster, and 19% of the total clusters are
singletons. The intersection of these county clusters with statistical regions
like Metropolitan Statistical Areas is high. Further, the clustering method
allows us to cover the entire USA as opposed to the limited commercially
active zones covered by the census bureau’s statistical area maps.
Additionally, it provides the flexibility to impose state border conditions on
the adjacency matrix; thus, supporting the roll-up of statistics to a state-
level when required. Figure 2 shows the county cluster maps of Alaska and New
York.

##### Reporting Noise:

Often, public health data sources avoid making retrospective changes to
correct reporting errors. As an example, Figure 3 shows that the cumulative
incidence data for Bristol County, RI, (FIPS code 44001) is not monotonically
increasing. The reason for such aberrations can be numerous, ranging from
sources incorrectly assigning new cases to a hospital (instead of the
patient’s home address) to changes in the way cases are counted. To resolve
these issues, we employ Isotonic regression (Chakravarti, 1989) defined by:

(1) $\min\sum_{i=1}^{n}w_{i}(x_{i}-a_{i})^{2}{\displaystyle{\text{ subject to
}}x_{i}\leq x_{j}{\text{ for all }}(i,j)\in E.}$

Figure 3. Isotonic Regression on Cumulative Incidence Cases

Isotonic regression involves finding a weighted least-squares fit
$x\in\mathbb{R}^{n}$ to a vector $a\in\mathbb{R}^{n}$ with weights vector
$w\in\mathbb{R}^{n}$ subject to a set of non-contradictory constraints of the
kind $x_{i}\leq x_{j}$. Figure 3 also shows the fit of Isotonic regression on
the noisy data.

While the Isotonic regression ensures the cumulative number of cases is
monotonically increasing, it does not safeguard against sudden spikes that
arise periodically due to a backlog in lab tests reported (e.g., on the
weekend) or correction in previously under-reported numbers. Since we model
both daily and cumulative incidences, these fluctuations would affect the
model performance. Towards that end, we use an adaptive-degree polynomial
filter (ADPF) (Barak, 1995) to smoothen the daily cases. It is shown that that
ADPF performs nearly as well as the optimally chosen fixed-degree Savitzky-
Golay filter and outperforms sub-optimally chosen Savitzky-Golay filters
(Barak, 1995). Note that we only calculate cumulative deaths as the daily
death numbers are small, and stochasticity makes daily optimization
problematic for deaths.

####  3.1.2. Auto Change Detection in Disease Dynamics:

As mentioned earlier, transmission rate defines the dynamics of the pandemic,
and this varies over time. The change in the transmission rate (and other
disease parameters) primarily stems from Non-Pharmaceutical Interventions
(NPIs) and changes in human behavior. NPIs include social distancing
practices, stay-at-home orders, night-time movement restrictions, restricted
public gatherings, etc. Unlike pharmaceutical interventions like vaccines that
reduce the susceptible population by moving individuals from $S$ to a
‘removed’ compartment (e.g., $R$), NPIs aim to lower the contact rate between
infectious and susceptible individuals.

While there are efforts to capture the impact and effectiveness of NPIs in
curtailing COVID-19 cases (Suryanarayanan et al., 2020; Sharov, 2020; Schwartz
et al., 2020), there are challenges in using them in a live system. Notably,
there are region-specific factors like long weekends or festivities, which
affect the disease parameters that are not usually part of the database. There
is also an overhead of maintaining an up-to-date list of all interventions.
Furthermore, transmission rate changes may happen due to human behavior as
well. Lastly, there is often a time delay between the implementation (or
relaxation) of an NPI and its effect on the observed cases and deaths. This
time delay may be unknown until changes in the cases and deaths are observed.
Our decision to auto-detect changes in the transmission rate stems from this
last observation. Namely, the entire history of the pandemic in a region is
split into multiple timepieces, where each piece represents the observed
variable of parameters in that time range. Note that in this work, we focus
only on transmission rate changes due to NPIs as there is not yet enough data
on vaccination.

Figure 4. Comparison of auto-detected intervention and fitted transmission
rate ($\beta$) with local NPIs and holidays

To identify changes in model parameters, we use Kneedle algorithm (Satopaa et
al., 2011) on the COVID-19 case data to detect inflection points. The
algorithm is run on smoothened daily case data. If, despite smoothening, the
algorithm detects multiple inflection points in close vicinity (within 7
days), they are merged. We also take into consideration the ‘significance’ of
each inflection point, wherein we not only measure the change in the slope of
the curve around the inflection point but also ensure that if a ‘knee’ is
detected, the subsequent values over two weeks are decreasing (vice-versa for
an ‘elbow’). Figure 4 illustrates the auto-detected intervention for county-
level data using our algorithm and superimposes the local government imposed
NPI.

These inflection points provide an approximate suggestion to our model
(Section 3.2.3) regarding times where it might have to (re)adjust parameters
to get an optimal fit. Note that this approach is just one of the two parts to
dynamically adapt time sensitive parameters like transmission rate. This part
handles scenarios where the effect of an intervention (or relaxation) is
already visible. The scenario where the effect is anticipated is handled by
incorporating mobility as a tertiary source (refer Section 3.2.2).

###  3.2. Base Model

####  3.2.1. Model Basics:

We use a compartment structure - a generic framework of models extensively
used in epidemiology (Kermack and McKendrick, 1927; Hu et al., 2013;
Balabdaoui and Mohr, 2020; Tsay et al., 2020), to model COVID-19 cases and
deaths. A compartment model includes the transition/flow of people across
different stages of a disease outbreak. At the minimum, a compartment flow
model contains two states - $S$ representing a group of susceptible people in
a population and $I$ representing infectious people. Depending on the type of
micro-parasitic infection, epidemiologists add compartments based on the
organism’s life cycle, patient state, and transmission mode. It is worth
noting that the ecosystem can span beyond humans; for instance, in the case of
Avian Influenza there is an interplay between humans and birds (and
mosquitoes) with different sets of parameters. In the case of COVID-19, we
would need at minimum two additional compartments to reflect people who
recovered ($R$) and died ($D$). It is to be noted that the compartment $R$
should be thought to be as ‘Removed’ compartment; i.e., if the $R$ compartment
is initialized with say 20% of the population in the region, it would then
amount to saying that 20% of the people would never get affected - indirectly
accounting for a concept called ‘herd immunity’. So depending on the
biological property of the virus, the ‘herd immunity’ ratio is set - for
COVID-19, we initialize $R$ with 20% of the total population (Fox et al.,
2020).

Figure 5. SAYCIWRDS Compartmental Model Flow

In our model, we add two additional ‘infectious’ compartments to account for
Asymptomats, including unreported infectious people, ($A$) and Pre-symptomats
($C$). Like the normal flu, it is reported (Johansson et al., 2021) that many
people infected by COVID-19 do not exhibit any symptoms. However, they
continue to shed viruses and contribute to the force of infection. Conversely,
some infected people with mild symptoms may not visit a doctor to be counted
in the COVID-19 infected list. The $A$ compartment represents these patients.

There is an incubation time for an individual to begin exhibiting the
symptoms. Even during the time when they are ‘pre-symptomatic’, such people
may shed the Sars-CoV-2 virus and contribute to the force of infection. It is
necessary to account for this ‘pre-symptomatic’ period in our model. It is
important to note that someone being marked COVID-19 positive corresponds to
the transition from $C\rightarrow I$. The time it takes for a person to move
from $C\rightarrow I$ corresponds to symptom appearance rate ($\alpha$).
Normally that represents just the incubation time, but in this case, it also
accounts for latency in testing and reporting.

Figure 5 represents our compartment model. The $S$ compartment holds the
susceptible population after accounting for the percentage of ‘removed’
population used to initialize the $R$ compartment. From the $S$ compartment,
an infected person can be either in $A$ or in $C$, which is determined by two
things: (a) Likelihood of the person getting infected - represented by the
transmission rate ($\beta$), and (b) Likelihood that patient ends up being
tested if not showing any symptoms - represented by the case reporting rate
($\xi$). As shown in Figure 5, we also introduced a compartment $W$, which
does not contribute to the force of infection. The $W$ compartment accounts
for patients who become worse and require extra attention and possibly
hospitalization. We assume the people in $W$ are isolated in hospital and
hence do not contribute to the force of infection. As this base model uses
only reported incidence and death data for fitting, $W$ acts as a delay
compartment to account for the time delay ($\omega$) between infection and
death. People can recover from $A$, $I$, $W$; albeit with different rates
$\gamma_{A},\gamma_{I},\gamma_{W}$. People enter the death compartment $D$ at
a rate of $\mu_{D}$. COVID-19 being similar to influenza, also has an
‘immunity loss rate’ ($\rho$), denoting the rate at which a person previously
infected may lose their immunity and become susceptible again. We set it to
1/(10 months) (Kissler et al., 2020).

####  3.2.2. Incorporating Mobility:

Figure 5 also shows an isolation compartment ($Y$). The $S$ and $Y$
compartments have a cyclic connection. The idea is to stage $Y$ as a transient
state to simulate people who are not mixing with the population. No
restriction on population mixing is an assumption made by standard compartment
models. However, with the stay-at-home order and decreased
mobility/interaction, this assumption is violated. Continuing that assumption
would lead to faster depletion of the susceptible population and
overestimation of the asymptomats (note that we fit only with respect to cases
and deaths, and so the $A$ compartment is indirectly determined by the
testing/case reporting rate $\xi$). To account for the population ‘hidden’
from virus exposure and scenarios where effects of relaxation of an
intervention are yet to reflect on observed cases (recalling discussions from
Section 3.1.2), we set the transition rate from $S\leftrightarrow Y$ to be
equal to rate of change in the daily mobility. This value is obtained from
tertiary sources like Apple222https://covid19.apple.com/mobility,
Google333https://www.google.com/covid19/mobility/, or Descartes
Lab444https://github.com/descarteslabs/DL-COVID-19. If the normalized
difference between two days is positive, it means the mobility is on the rise
and the flow from $Y\rightarrow S$ is triggered, else from $S\rightarrow Y$.
Thus, at any point in time, the transition between $S$ and $Y$ is always
unidirectional.

The advantage of this model is that we are able to include mobility data
without any additional tunable parameter(s). Note that in theory, the
isolation compartment $Y$ can be attached to $A$ and $C$ with exits to $R$ and
$I$, respectively. However, this will lead to over-complication of the model
and an additional parameter $Y_{A}\rightarrow R$.

####  3.2.3. Model Formulation:

Here, we provide the mathematical description of the compartmental flow
illustrated in Figure 5. These differential equations represent changes in the
state of population (each compartment) at any time instance. As the time scale
is continuous, we need to add interpolation functions on any of our tertiary
sources that directly control the flow - e.g., changes in the mobility that
controls $S\leftrightarrow Y$. Let $f_{mob}$ be a function fitted on the
mobility, then the equation from $S\leftrightarrow Y$ can be described as:

(2) $\displaystyle T_{Y\rightarrow S}=max(0,f_{mob}(t))$ (3) $\displaystyle
T_{S\rightarrow Y}=max(0,-1*f_{mob}(t))$

The fitted function $f_{mob}$ can also be used to extrapolate mobility data to
get likely mobility values in the future. Now, the compartmental flow rates
can be defined as:

(4) $\displaystyle\frac{dS}{dt}$ $\displaystyle=$
$\displaystyle-\beta*\frac{S*(I+A+C)}{N}+\rho*R+min(Y,T_{Y\rightarrow
S}*(S+Y))$ $\displaystyle-min(S,T_{S\rightarrow Y}*(S+Y))$ (5)
$\displaystyle\frac{dY}{dt}$ $\displaystyle=$ $\displaystyle
min(S,T_{S\rightarrow Y}*(S+Y))-min(Y,T_{Y\rightarrow S}*(S+Y))$ (6)
$\displaystyle\frac{dA}{dt}$ $\displaystyle=$
$\displaystyle(1-\xi)*(\beta*S*(I+A+C)/N)-\gamma_{A}*A$ (7)
$\displaystyle\frac{dC}{dt}$ $\displaystyle=$
$\displaystyle(\xi)*(\beta*S*(I+A+C)/N)-\alpha*C$ (8)
$\displaystyle\frac{dI}{dt}$ $\displaystyle=$
$\displaystyle\alpha*C-(\gamma_{I}+\omega)*I$ (9) $\displaystyle\frac{dW}{dt}$
$\displaystyle=$ $\displaystyle\omega*I-(\mu_{d}+\gamma_{W})*W$ (10)
$\displaystyle\frac{dR}{dt}$ $\displaystyle=$
$\displaystyle\gamma_{I}*I+\gamma_{A}*A+\gamma_{W}*W-\rho*R$ (11)
$\displaystyle\frac{dD}{dt}$ $\displaystyle=$ $\displaystyle\mu_{D}*W$

Note that the $\beta$ in the above equations are time-varying parameters as
explained in Section 3.1.2 \- i.e., there will be one $\beta$ for each of the
inflection-detected time intervals, and selecting which beta to modify is
based on the time stamp $t$ of the ordinary differential equation (ODE)
solver. Recall that the approach in Section 3.1.2 only provides an initial
estimate for the time of the inflection. Since the ODE solver iterates over
time steps, it can find the ‘optimal’ $t$ from the rough estimate to fit the
piece-wise curve. Through this, we are able to mimic the natural evolution of
the disease, and unlike a regression algorithm, this does not violate any
epidemiological principle. We also set $\xi$ ($S\rightarrow A,C$), $\omega$
($I\rightarrow W$), $\mu_{d}$ ($W\rightarrow D$) to be a time-sensitive
parameter similar to $\beta$.

####  3.2.4. Model Solver:

We use an LSODA algorithm (Petzold, 1983; Hindmarsh, 1983; Hindmarsh and
Petzold, 2005) to solve the ODEs. LSODA automatically selects between non-
stiff (Adams) and stiff (BDF) methods. It uses the non-stiff method initially
and dynamically monitors data to decide which method to use. The parameter
estimation is done by Levenberg–Marquardt algorithm (Levenberg, 1944;
Marquardt, 1963); selected mainly for its simplicity and speed. The error
function to minimize includes daily case NRMSE, cumulative case NRMSE, and
cumulative death NRMSE. Because the entire system of equations receives
continuous external ‘shocks’ via the mobility values, the solver takes
additional time to converge. However, as discussed in Section 4.1.2, including
this external data provides valuable information and improves performance.

###  3.3. Client Independent Analytics

Our solution provides a set of analytics regardless of which client use-case
it is serving. These typically include epidemiological metrics like $R_{0}$,
doubling time, 14-day rolling average, etc. Amongst these, $R_{0}$ is one of
great importance. It provides a likelihood of how much an infected individual
is likely to infect others. It is directly connected to the transmission rate
($\beta$) and recovery time ($\gamma$). However, it is also model specific.
Solving Equations 6 to Equation 11, under steady-state, we derive $R_{0}$:

(12) $\displaystyle R_{0}$
$\displaystyle=\frac{\beta}{\gamma_{A}}\left\\{1+\left(\frac{\gamma_{A}}{\gamma_{I}+\omega}\right)-\left[1-\frac{\gamma_{A}}{\alpha}\right]\xi\right\\}$

where the $\beta$, $\gamma_{A}$, and $\gamma_{I}$ in the above equation refer
to the first fitted value of $\beta$, $\gamma_{A}$, and $\gamma_{I}$. Recall
that we fit multiple $\beta$ mimicking the different transmission rates over
the course of the pandemic in a region. If we calculate $R$ at each of those
different times, then using the fraction of susceptibles during each of those
times, we get $R-$effective ($R_{eff}$). $R_{t}$ is used to denote the latest
$R_{eff}$.

Together with doubling time (Equation 13), these two statistics provide an
accurate sense of the pandemic situation on the ground.

(13) ${\text{doubling time}}={\frac{\ln(2)}{\text{growth\\_rate}}}$

where $growth\\_rate$ is defined as

(14) ${\text{growth\\_rate}}=\frac{\ln\left(N(t)/N(0)\right)}{t}$

####  3.3.1. Community Risk Evaluation

The community risk evaluation uses a rule-based approach to convert the
epidemiological metrics into risk scores for a region between 1-6, with 1
being the safest. The features that determine the score include recent
incidence trends, the case prevalence in the population, and local government
reopening policy. Of these, the case prevalence and government reopening
policies determine the manual time-dependant adjustable thresholds in our
algorithm. Algorithm 1 describes the algorithm to determine the community
risk. In our algorithm, there are 3 thresholds - a soft ($\kappa$), a hard
($\lambda$), and a flat threshold ($\tau$), that help merge the three features
into one value. For the U.S.A, based on White House reopening guidelines, we
identified $\kappa$ to be 10 cases/100K, $\lambda$ to be 5 cases/100K, and
$\tau$ to be 2 cases/100K below which week-over-week change is considered
flat. These values were adjusted to account for local reopening guidelines
like the traffic light system used in Mexico and 5 alert levels used in the
U.K.

Algorithm 1 Community Risk Determination Algorithm

1:Input: 21 days historical data ($A_{1}$, $A_{2}$, $A_{3}$, with $A_{3}$
being the latest weekly average) and predictions
($A^{{}^{\prime}}_{1}$,$A^{{}^{\prime}}_{2}$,$A^{{}^{\prime}}_{3}$), $\kappa$,
$\lambda$, $\tau$

2:Output: Current and 3 weeks projected community risk scores

3:is_strict_decr = $A_{3}<A_{2}\land A_{2}<A_{1}$

4:is_strict_incr = $A_{3}>A_{2}\land A_{2}>A_{1}$

5:if $A_{3}$, $A_{2}$, $A_{1}$ ¡ $\kappa$ then

6: if $A_{3}$, $A_{2}$¡$\lambda$ then

7: is_flat = $abs(A_{3}-A_{2}),abs(A_{2}-A_{1}),abs(A_{3}-A_{1})<=\tau$

8: is_flat_decr = $abs(A_{2}-A_{1})<=\tau\land(A_{3}<A_{2})$

9: if $is\\_flat\lor is\\_strict\\_decr\lor is\\_flat\\_decr$ then

10: curr_risk_score = 1.0

11: else curr_risk_score = 2.0

12: else if $is\\_strict\\_decr$ then return curr_risk_score = 2.0

13: else curr_risk_score = 3.0

14:else

15: if $is\\_strict\\_decr$ then curr_risk_score = 4.0

16: else if  $is\\_strict\\_incr$ then curr_risk_score = 6.0

17: else curr_risk_score = 5.0

18:Repeat Steps 3-17 for
$A^{{}^{\prime}}_{1}$,$A^{{}^{\prime}}_{2}$,$A^{{}^{\prime}}_{3}$ to get
projected risk scores

###  3.4. Client Specific Analytics

As mentioned above, our base model could be easily extended for client-
specific use cases. In one of the client engagements, we augmented the base
model with hospitalization data to predict both hospitalization and ICU
demand. Figure 6 shows the model extension. The input to this extension is the
predicted incidence from our base model. Note that we do not use the predicted
deaths from the base model as-is, because in this case, we need to infer the
mortality rate from the ICU training data. The hospitalization, ICU, and death
rates are learnt using the same infrastructure used in solving the base model.
The equations can be derived similarly to our base model described in Section
3.2.3.

As shown in Section 4.2.3, extending to simulate what-if scenarios are also
straightforward, as one can adjust the learnt parameters and run the model to
solve the ODEs in a ‘scoring’ mode.

Figure 6. Hospitalization Model Add-on | AUGUST - 3rd | OCTOBER - 26th
---|---|---
Rank | Week-1 | Week-1+2 | Week-1+2+3 | Week-1 | Week-1+2 | Week-1+2+3
1 | IBM (0.24) | IBM (0.23) | LANL (0.26) | IBM (0.19) | IBM (0.162) | IBM (0.223)
2 | LANL (0.42) | LANL (0.27) | IBM (0.27) | ESG (0.216) | Oliver Wyman (0.174) | Oliver Wyman (0.228)
3 | LNQ (0.52) | LNQ (0.31) | LNQ (0.31) | COVID-19 Sim (0.217) | Karlen (0.1851) | LNQ (0.244)
4 | DDS-NBDS (0.54) | DDS-NBDS (0.34) | DDS-NBDS (0.34) | UCLA (0.226) | LNQ (0.1853) | Karlen (0.248)
5 | ERDC (0.56) | ERDC (0.36) | USC (0.35) | UGA-CEID (0.23) | Google-HSPH (0.22) | DDS (0.27)
Table 1. Incidence MAPE

##  4\. Results

Having described the relevant algorithms and modules, we now present a
detailed evaluation analysis of the results of the various modules in the
system. Our experiments can be broken down into two broad categories wherein
we first measure the prediction performance of our base/core system, and then
through various user story scenarios, explain the performance, usefulness, and
dexterity of the various analytic modules.

##### Base Module Evaluation Scheme:

We use Mean Absolute Percentage Errors (MAPE) to quantitatively measure the
prediction performance. We compare our results with the baselines presented in
CDC555www.cdc.gov/coronavirus/2019-ncov/covid-data/mathematical-modeling.html.
To avoid bias in measuring performance, we measure the accuracy at both the
falling phase (end of first wave ~Aug-3rd) and the rising curve (start of the
second wave ~Oct-26th). Note that CDC did not report case prediction before
Aug-1st. We also perform a time series ablation study to measure the
importance of mobility data and then wrap up the base module analysis by
presenting geo-spatial and running time analysis.

##### Analytic Module Evaluation Scheme:

We consider real use-cases to characterize the performance of the ‘Client
Independent Analytics’ and ‘Client Specific Analytics’ modules. In one of the
client engagements pertaining to the state of Rhode-Island, we highlight how
the compartment model helps track the disease progression in the region and
compare the validity of the estimated community risk scores with the situation
on the ground. Through our support for a major hospital in Tampa, FL, we
evaluate the hospital and ICU demand prediction module. As part of our
engagement, we also simulate what-if scenarios in response to the Thanksgiving
and holiday season, enumerating the effectiveness of various counter-measures.

###  4.1. Base Model Evaluation

####  4.1.1. Prediction Evaluation:

Table 1 shows the cumulative weekly incidence MAPEs calculated from the end of
training dates (August-3rd and October-26th). As we can see, the ‘IBM’ system
consistently comes in top-3, and for the week-1 and week-1+2, it is at the
top. The numbers in the bracket indicate the MAPE, and we can see that for the
week-1, the difference is substantial. Note that for October, the rest of the
top-5 baselines do not make it into top-5 for the remaining weeks. It is also
to be noted that some of the methodologies like Oliver Wyman started
publishing their results in September and hence do not appear in the August
comparison.

Table 2 shows the cumulative deaths MAPE. In the interest of space, we report
only the October-26th value. We observe a similar trend like Table 1 where our
model consistently appears in top-2.

We also noted that some reported methodologies predict only cases or deaths
(e.g., YYG). Based on Table 1 and 2, IBM and Oliver Wyman are comparable
models with IBM doing particularly well in the 1st week. There are two
important points to remember while interpreting these tables. First, death
numbers are substantially smaller than incidence, so the sensitivity of MAPE
to incorrect predictions is higher. Second, as the MAPE is calculated based on
weekly values, a small deviation in the reported ground truth numbers at the
end of an evaluation week may adversely impact the MAPE calculation itself.

Rank |  Week-1 |  Week-1+2 |  Week-1+2+3
---|---|---|---
1 |  IBM (0.0146) |  Columbia (0.0146) |  IBM (0.026)
2 |  Oliver Wyman (0.0167) |  IBM (0.0189) |  MSRA (0.36)
3 |  CovidComplete (0.017) |  MOBS (0.0245) |  Oliver Wyman (0.037)
4 |  MOBS (0.01745) |  CovidComplete (0.0247) |  Karlen (0.038)
5 |  USC (0.01779) |  Oliver Wyman (0.025) |  UMass-MB (0.0384)
Table 2. Death MAPE for October-26th

####  4.1.2. Impact of Mobility:

Table 3 compares the MAPE of the model with and without mobility for multiple
runs. While there is an improvement across the board, we note a 11%
improvement in the pandemic’s early days. Figure 7 shows how the mobility data
is able to adjust the trajectory more accurately. While some studies have
already suggested that mobility is a leading indicator for COVID-19 cases
(Miller et al., 2020), in this work, we further show how to incorporate that
information in practice and its usefulness in not only improving the accuracy
but also understanding the population mixing and the susceptible pool. We also
observed that the contribution of mobility data towards prediction has waned
over time. We suspect this may be due to increased awareness, mask adoption,
and increasing herd immunity (i.e., as the population develops antibodies from
infection or vaccination). Thus, it is our learning that using mobility data,
if available, is a must to get accurate predictions during the early stages of
a pandemic. It provides a more precise estimate of susceptibles and the
magnitude of the potential second wave.

| Week-1 | Week-2 | Week-3 | Week-4
---|---|---|---|---
Aug-3 | 9.25 | 10.38 | 12.78 | 15.29
Sept-14 | 2.63 | 1.82 | 1.54 | 1.12
Oct-26 | 2 | 1.59 | 0.71 | 0.36
Table 3. Percentage Improvement In Incidence MAPE by Incorporating Mobility
Data Figure 7. Predictions with and without Mobility

####  4.1.3. Impact of Geo-Location:

The choice of geo-unit for analysis is important. As shown in Table 4,
modeling at a state level as opposed to aggregating the projections from
county-cluster leads to sub-par MAPE. This is because the local intricacies of
population mixing are lost when the entire state is considered a well-mixed
unit. Thus, crafting an accurate policy at a higher geo-level requires hyper-
local modeling. A similar conclusion also holds true when using the state
projections to estimate county clusters using area population to drill down.
In general, we recommend modeling at county-clusters granularity during the
pandemic’s initial stages to avoid data sparsity issues. Once the counties
within the cluster enter the exponential growth phase (vs. sporadic or
intermittent case reports), modeling at the county level will prove
beneficial, albeit with a trade-off on the running time.

Average[(State MAPE)/ (County Cluster Aggregated to State MAPE)] |  Average[(State Drilled Down to County Cluster MAPE)/(County Cluster MAPE)]
---|---
Incidence | Deaths | Incidence | Death
1.44 | 1.105 | 1.818 | 9.27
Table 4. Comparison of Effect of Geo-Unit on Incidence and Death MAPEs

####  4.1.4. Running Time Efficiency:

The running time of the pipeline is determined by the choice of parameter
optimization algorithm and the external “shocks” introduced through tertiary
sources like mobility data. These shocks cause a sudden change in the
equilibrium of the system. Furthermore, with the increase in the number of
timestamps (and consequently the number of time-sensitive parameters), the
effect is non-linear. Table 5 shows the running time of finding the best
parameters for 20 initializer combinations.

Time (mins) | 05/01 | 08/03 | 09/14 | 10/26
---|---|---|---|---
No Mobility | 1.317 | 13.33 | 19.82 | 30.03
Mobility | 40.61 | 171.54 | 273.03 | 404.33
Table 5. Running Time Comparison With and Without Mobility Data for 20
Initializer Combinations

The time-sensitive parameters also make our ODEs stiff at certain places and
our experiments with different differential equation solvers, including Runge-
Kutta, BDF, and LSODA, showed that LSODA gave the best performance.

We also evaluated different parameter estimation algorithms. Table 6 shows the
efficiency of the different algorithms with respect to Levenberg–Marquardt
algorithm, determined to be the best choice. Note that we ensured convergence
of approximately the same training error when measuring the running time
efficiency.

| Levenberg – Marquardt | SLSQP | Nelder-Mead | BFGS
---|---|---|---|---
Times Slower | 1X | 3.5X | 5X | 5X
Table 6. Running Time Comparison of Different Parameter Estimation Algorithms

###  4.2. Analytics Module Evaluation

####  4.2.1. Community Risk Prediction - Case Study:

The tri-state region of Rhode Island, Connecticut, and Massachusetts share a
common border and worked in tandem to enact restrictions in March-April to
control the pandemic. However, as shown in Figure 8, their end date of
restrictions are not similar. In the majority of the cases, RI had relaxed
their restriction 2-3 weeks earlier.

Based on the caseloads and predictions at the time, our solution estimated
that in mid-June early-July, CT and MA were safe to reopen with risk values of
around 1-2. RI was still marked as unsafe with a score of 5. These scores were
validated when the cases started increasing for RI in August while CT and MA
were still doing okay.

To validate these numbers, we examine the model’s compartmental values around
the time of June reopening. Figure 8 shows the population normalized fitted
compartment values of $S$ (susceptible), $Y$ (isolation), and $I+W$ (total
infectious) from March-1st to October-20th. It also superimposes the end date
of the restrictions. Comparing the prevalence (normalized infection load)
across states at the time of reopening, we can observe that RI’s prevalence
was high, and it stayed relatively high throughout the summer. Positivity
rates also remained high (note that this also depends on state policy on who
gets tested). Looking at the normalized fraction in $S$ and $Y$, we see that
all states are comparable. Thus, the higher numbers in RI (from August
onwards) have to do with the high sustained prevalence, and this was in part
because of early relaxation of restrictions (15 days to 1 month in advance).
This not only validates that the system correctly estimated the risk scores
for RI but also shows that CT would most likely have a second wave earlier
than MA again. As expected, the risk score of RI further worsened to 6, while
CT and MA eventually slipped to risk level 6 from 3 on October-27th and
October-29th, respectively.

Figure 8. Population Normalized $S$, $I$ and $I+W$ Compartment Values for RI,
CT & MA with reopening dates

####  4.2.2. Predicting Hospitalization and ICU Demand:

Figure 9 shows the fitting and the prediction of hospitalization and ICU beds
for Hillsborough County, FL. The figure is for an engagement between April to
November 2020, where we provided projections at a regular cadence (indicated
in the legend). The model accuracy predicted a peak incidence on 7/20/2020 and
peak ICU on 7/26/2020. The client used these projections to chart a hospital
capacity forecast and take steps to handle the surge. The blue zone above and
below the total ICU prediction represents the MAPE (5.9% on average). Based on
the projections, the hospital decided not to cancel other elective procedures
and surgeries in the time period shown.

Figure 9. Predicting Hospitalization & ICU Demands

####  4.2.3. Counterfactual Analysis - Hillsborough Use Case:

Counterfactual analysis is important in simulating what-if scenarios. Clients,
especially local authorities like hospitals or government, would like to
understand the effects of future events (including seasonal holidays) on the
state of the pandemic. At the same time, it also allows them to plan and take
remedial steps. For instance, Figure 10 shows the base scenario wherein we
show the prediction based on mobility and training as of October-$26^{th}$,
2020 for Hillsborough County, FL. With the Thanksgiving holidays around the
end of November, local government officials may be interested in knowing ways
to reduce the spread of COVID-19. As mentioned earlier, the transmission of a
virus is primarily based on population mixing, and one potent way to curb the
population mixing is to enforce lock-downs. However, lock-downs have a
detrimental effect on the economy, and so one needs to thoroughly understand
the implication of different degrees of enforcement on the spread so that one
can take the most appropriate step without overreacting.

Figure 10 shows the change in daily COVID-19 cases with varying mobility. Note
that we do not change any other parameters obtained from the base scenario.
The most potent intervention for this state of the epidemic appears to be a
10% reduction in the mobility values observed on October-$26^{th}$. However,
that would require a significant level of curtailment. By comparing the
historical case trends, the 7% decrement scenario matches structurally to the
events observed in July. A 5% reduction achieves the effect of stay-at-home;
however, it leads to a slower decline in the subsequent epidemic wave (gray
curve). The figure also reveals a minimum threshold level for mobility
curtailment of ¿~2%. In the case of 2% reduction, although the mobility value
has decreased, there is still an insufficient movement of people into the $Y$
compartment to significantly reduce future incidence over the base scenario.
Only for higher values does mobility reduction impact future cases. This form
of what-if analysis allows the authorities to assess multiple possible
interventions and provides insight into their relative efficacies based on
measurements of the effects of past NPI’s.

Figure 10. Simulating what-if Scenarios with Different Mobility Restrictions -
Hillsborough County, FL

##  5\. Deployment and Conclusion

The solution uses Python for modeling and JAVA for data orchestration. We
built a Docker container and deployed it in IBM Cloud Kubernetes Cluster so
that each geo-unit can be processed by an image independently. Apache NiFi
orchestrates the entire data flow. There are two NiFi pipelines, one for the
core epidemiological module, which is run once in three days. The second is
for the ‘Client Independent Analytics’, a.k.a community risk pipeline, which
is run nightly. A manifest configuration file containing the list of geo-ids
and any list of hyper-parameters like model compartment initializer values is
used to spin one new pod on the cluster for the data flow of a geo-unit. This
step is throttled to 18 geo-units at a time due to cluster performance
limitations. Note that each pod itself is highly parallelized to solve
multiple compartment initializer combinations concurrently. On completion of a
pod, the NiFi pipeline stores the output of the epidemiology run in an IBM COS
bucket. The Community Risk NiFi pipeline, which is independently triggered
nightly, collects the results from the COS bucket and then triggers the
process to download the latest case data and merge it with the epidemiological
predictions. After computing the statistics and risk scores, it finally pushes
the output to DB2 On Cloud through a Data Access Layer (DAL) for clients to
consume. The ‘Client Specific Analytics’ has a separate private pipeline that
uses these outputs with any additional client data. Figure 11 shows the final
UI presentation to the end-user.

In the near future, we intend to extend this solution with additional signals
like vaccines, demography, regional comorbidity factors, etc. We strongly
believe that the solution and findings of this paper would render
epidemiologists and decision-makers better equipped under similar conditions
in the future.

Figure 11. User Interface of the Solution

## References

  * (1)
  * Balabdaoui and Mohr (2020) Fadoua Balabdaoui and Dirk Mohr. 2020.  Age-stratified discrete compartment model of the COVID-19 epidemic with application to Switzerland.  _Scientific reports_ 10, 1 (2020), 1–12.
  * Barak (1995) Phillip Barak. 1995\.  Smoothing and differentiation by an adaptive-degree polynomial filter.  _Analytical Chemistry_ 67, 17 (1995), 2758–2762.
  * Blondel et al. (2008) Vincent D Blondel, Jean-Loup Guillaume, Renaud Lambiotte, and Etienne Lefebvre. 2008.  Fast unfolding of communities in large networks.  _Journal of Statistical Mechanics: Theory and Experiment_ 2008, 10 (2008), P10008.
  * Brooks et al. ([n.d.]) Logan C Brooks, Evan L Ray, Jacob Bien, Johannes Bracher, Aaron Rumack, Ryan J Tibshirani, and Nicholas G Reich. [n.d.].  Comparing ensemble approaches for short-term probabilistic COVID-19 forecasts in the US.  ([n. d.]).
  * Chakravarti (1989) Nilotpal Chakravarti. 1989\.  Isotonic median regression: a linear programming approach.  _Mathematics of operations research_ 14, 2 (1989), 303–308.
  * COVID et al. (2020) IHME COVID, Christopher JL Murray, et al. 2020\.  Forecasting COVID-19 impact on hospital bed-days, ICU-days, ventilator-days and deaths by US state in the next 4 months.  _MedRxiv_ (2020).
  * Fox et al. (2020) Spencer J Fox, Pratyush Potu, Ravi Srinivasan, Michael Lachmann, and Lauren Ancel Meyers. 2020\.  The COVID-19 herd immunity threshold is not low: A re-analysis of European data from spring of 2020.  _medRxiv_ (2020).
  * Hindmarsh and Petzold (2005) AC Hindmarsh and LR Petzold. 2005.  LSODA, ordinary differential equation solver for stiff or non-stiff system.  (2005).
  * Hindmarsh (1983) Alan C Hindmarsh. 1983\.  ODEPACK, a systematized collection of ODE solvers.  _Scientific computing_ (1983), 55–64.
  * Hu et al. (2013) Kun Hu, Christian Thoens, Simone Bianco, Stefan Edlund, Matthew Davis, Judith Douglas, and James Kaufman. 2013.  Modeling the dynamics of dengue fever. In _International Conference on Social Computing, Behavioral-Cultural Modeling, and Prediction_. Springer, 486–494.
  * ILO and WHO (2020) IFAD ILO, FAO and WHO. 2020\.  _Impact of COVID-19 on people’s livelihoods, their health and our food systems_.  Retrieved January 17, 2021 from https://www.who.int/news/item/13-10-2020-impact-of-covid-19-on-people%27s-livelihoods-their-health-and-our-food-systems
  * Johansson et al. (2021) Michael A Johansson, Talia M Quandelacy, Sarah Kada, Pragati Venkata Prasad, Molly Steele, John T Brooks, Rachel B Slayton, Matthew Biggerstaff, and Jay C Butler. 2021\.  SARS-CoV-2 Transmission From People Without COVID-19 Symptoms.  _JAMA network open_ 4, 1 (2021), e2035057–e2035057.
  * Kermack and McKendrick (1927) William Ogilvy Kermack and Anderson G McKendrick. 1927.  A contribution to the mathematical theory of epidemics.  _Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character_ 115, 772 (1927), 700–721.
  * Kissler et al. (2020) Stephen M Kissler, Christine Tedijanto, Edward Goldstein, Yonatan H Grad, and Marc Lipsitch. 2020\.  Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period.  _Science_ 368, 6493 (2020), 860–868.
  * Levenberg (1944) Kenneth Levenberg. 1944\.  A method for the solution of certain non-linear problems in least squares.  _Quarterly of applied mathematics_ 2, 2 (1944), 164–168.
  * Li et al. (2020) Michael Li, Hamza T Bouardi, Omar Lami, Nikolaos Trichakis, Thomas Trikalinos, Mohammad Fazel Zarandi, and Dimitris Bertsimas. 2020\.  _Overview of DELPHI Model V3 - COVIDAnalytics_.  MIT Research Lab Technical Report. Massachusetts Institute of Technology, Cambridge, MA.
  * Marquardt (1963) Donald W Marquardt. 1963\.  An algorithm for least-squares estimation of nonlinear parameters.  _Journal of the society for Industrial and Applied Mathematics_ 11, 2 (1963), 431–441.
  * Meredith et al. (2020) J Wayne Meredith, Kevin P High, and Julie Ann Freischlag. 2020\.  Preserving elective surgeries in the COVID-19 pandemic and the future.  _Jama_ 324, 17 (2020), 1725–1726.
  * Miller et al. (2020) Andrew C Miller, Nicholas J Foti, Joseph A Lewnard, Nicholas P Jewell, Carlos Guestrin, and Emily B Fox. 2020.  Mobility trends provide a leading indicator of changes in SARS-CoV-2 transmission.  _medRxiv_ (2020).
  * Petzold (1983) Linda Petzold. 1983\.  Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations.  _SIAM journal on scientific and statistical computing_ 4, 1 (1983), 136–148.
  * Satopaa et al. (2011) Ville Satopaa, Jeannie Albrecht, David Irwin, and Barath Raghavan. 2011.  Finding a” kneedle” in a haystack: Detecting knee points in system behavior. In _2011 31st international conference on distributed computing systems workshops_. IEEE, 166–171.
  * Schwartz et al. (2020) Ira B Schwartz, James H Kaufman, Kun Hu, and Simone Bianco. 2020\.  Predicting the impact of asymptomatic transmission, non-pharmaceutical intervention and testing on the spread of COVID19 COVID19.  _medRxiv_ (2020).
  * Sharov (2020) Konstantin S Sharov. 2020\.  Creating and applying SIR modified compartmental model for calculation of COVID-19 lockdown efficiency.  _Chaos, Solitons & Fractals_ 141 (2020), 110295.
  * Shen et al. (2020) Huayu Shen, Mengyao Fu, Hongyu Pan, Zhongfu Yu, and Yongquan Chen. 2020.  The impact of the COVID-19 pandemic on firm performance.  _Emerging Markets Finance and Trade_ 56, 10 (2020), 2213–2230.
  * Smalley (2020) Andrew Smalley. 2020\.  _Higher Education Responses to Coronavirus (COVID-19)_.  Retrieved January 17, 2021 from https://www.ncsl.org/research/education/higher-education-responses-to-coronavirus-covid-19.aspx
  * Suryanarayanan et al. (2020) Parthasarathy Suryanarayanan, Ching-Huei Tsou, Ananya Poddar, Diwakar Mahajan, Bharath Dandala, Piyush Madan, Anshul Agrawal, Charles Wachira, Osebe Mogaka Samuel, Osnat Bar-Shira, et al. 2020\.  WNTRAC: Artificial Intelligence Assisted Tracking of Non-pharmaceutical Interventions Implemented Worldwide for COVID-19.  _arXiv preprint arXiv:2009.07057_ (2020).
  * Team (2020) Oliver Wyman Team. 2020\.  _COVID-19 Pandemic Navigator Core Model: Overview of methodology and use cases_.  Working Paper. Oliver Wyman.
  * Tsay et al. (2020) Calvin Tsay, Fernando Lejarza, Mark A. Stadtherr, and Michael Baldea. 2020.  Modeling, state estimation, and optimal control for the US COVID-19 outbreak.  _Scientific Reports_ 10, 1 (01 Jul 2020), 10711\.
  * Zou et al. (2020) Difan Zou, Lingxiao Wang, Pan Xu, Jinghui Chen, Weitong Zhang, and Quanquan Gu. 2020\.  Epidemic Model Guided Machine Learning for COVID-19 Forecasts in the United States.  _medRxiv_ (2020).