# Anomaly Ratio Distributions of Hadronic Axion Models with Multiple Heavy
Quarks

Vaisakh Plakkot<EMAIL_ADDRESS>Sebastian Hoof
<EMAIL_ADDRESS>Institut für Astrophysik, Georg-August-Universität
Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany

(September 2021)

###### Abstract

We consider hadronic axion models that extend the Standard Model by one
complex scalar field and one or more new heavy quarks, i.e.
$N_{\mathcal{Q}}\geq 1$. We review previously suggested selection criteria as
well as categorize and catalog all possible models for $N_{\mathcal{Q}}\leq
9$. In particular, allowing for $N_{\mathcal{Q}}>1$ can introduce models that
spoil the axion solution of the strong CP problem. Demanding that Landau poles
do not appear below some energy scale limits the number of preferred models to
a finite number. For our choice of criteria, we find that $N_{\mathcal{Q}}\leq
28$ and only 820 different anomaly ratios $E/N$ exist (443 when considering
additive representations, 12 when all new quarks transform under the same
representation). We analyze the ensuing $E/N$ distributions, which can be used
to construct informative priors on the axion-photon coupling. The hadronic
axion model band may be defined as the central region of one of these
distributions, and we show how the band for equally probable, preferred models
compares to present and future experimental constraints.

##  I Introduction

QCD axions [1, 2], initially proposed as a solution to the strong CP problem
[3, 4], are excellent cold dark matter (DM) candidates [5, 6, 7, 8, 9].
Numerous experimental searches are currently underway to find such particles
[10]. One major challenge of axion detection is that the axion mass is set by
an unknown parameter, the axion decay constant $f_{a}$, which can range across
many orders of magnitude. Moreover, the axion’s interactions with the Standard
Model (SM) are usually model-dependent, and a UV axion model has to be
constructed in order to determine the exact relationship of $f_{a}$ and the
axion couplings.

One class of such UV models are hadronic (also called KSVZ-type) axion models
[11, 12], which extend the SM by a new complex scalar field and
$N_{\mathcal{Q}}\geq 1$ heavy, exotic quarks. For a given value of
$N_{\mathcal{Q}}$ there exist multiple, discrete models, which trace out lines
in the axion mass and axion-photon coupling parameter space. The locations of
these lines are determined by the anomaly ratio $E/N$ and a model-independent
contribution from axion-meson mixing.

To map and restrict the resulting landscape of axion models, it has been
suggested that phenomenological selection criteria can be used to single out
_preferred_ models [13, 14]. This allows us to restrict the parameter space
and helps experiments to assess their sensitivity requirements. However, so
far only the case of $N_{\mathcal{Q}}=1$ has been fully cataloged, which is
why we want to study models with $N_{\mathcal{Q}}>1$ as far as this is
feasible. First, we summarize the construction of KSVZ-type axion models and
phenomenological selection criteria in Secs. II and III. Subsequently, a
catalog of all possible models with $N_{\mathcal{Q}}\leq 9$ is presented and
the resulting $E/N$ distributions are discussed. We catalog all _preferred_
models, for which we find that the maximum possible number of $\mathcal{Q}$s
is $N_{\mathcal{Q}}=28$. In Sec. V we outline how the catalog of models can be
used to construct informative prior distributions on $E/N$. These can be used
to define the KSVZ axion model band and we show how it compares to current and
future experimental constraints. Finally, we summarize our work and end with
some closing remarks. Model catalogs and further supplementary material are
available on Zenodo [15].

##  II Hadronic axion models

Let us denote a representation of a particle as
$(\mathcal{C},\mathcal{I},\mathcal{Y})$, where $\mathcal{C}$ and $\mathcal{I}$
are the $\mathrm{SU}(3)_{\mathcal{C}}$ color and
$\mathrm{SU}(2)_{\mathcal{I}}$ isospin representations, respectively, while
$\mathcal{Y}$ denotes the particle’s $\mathrm{U}(1)_{\mathcal{Y}}$
hypercharge.

For example, the traditional KSVZ axion model contains a heavy chiral quark
$\mathcal{Q}=\mathcal{Q}_{L}+\mathcal{Q}_{R}\sim(3,1,0)$, charged under the
$\mathrm{U}(1)_{\text{PQ}}$ Peccei-Quinn (PQ) symmetry with charge
$\mathcal{X}=\mathcal{X}_{L}-\mathcal{X}_{R}=\pm 1$, and the complex scalar
field $\Phi\sim(1,1,0)$ with PQ charge normalized to $\mathcal{X}_{\Phi}=1$.
All SM fields are uncharged under the PQ symmetry in the KSVZ model, and the
relevant part of the Lagrangian is

$\displaystyle\mathcal{L}\supset\ $ $\displaystyle
i\,\overline{\mathcal{Q}}\,\gamma^{\mu}D_{\mu}\mathcal{Q}-(y_{\mathcal{Q}}\overline{\mathcal{Q}}_{L}\mathcal{Q}_{R}\Phi+\text{h.c.})$
$\displaystyle-\lambda_{\Phi}\left(|\Phi|^{2}-\frac{v_{a}^{2}}{2}\right)^{2},$
(1)

where $y_{\mathcal{Q}}$ is the Yukawa coupling constant and the last term is a
potential for the complex scalar field with order parameter $v_{a}$. The
Lagrangian is invariant under a chiral $\mathrm{U}(1)_{\text{PQ}}$
transformation $\Phi\mapsto\mathrm{e}^{i\alpha}\Phi$,
$\mathcal{Q}_{L/R}\mapsto\mathrm{e}^{\pm i\alpha/2}\mathcal{Q}_{L/R}$. The
field $\Phi$ attains a non-zero value at the minimum of the potential,
resulting in a spontaneously broken PQ symmetry. Expanding $\Phi$ around its
vacuum expectation value gives the axion as the corresponding angular degree
of freedom, with value in the interval $[0,2\pi v_{a})$. The mass of
$\mathcal{Q}$ is then $m_{\mathcal{Q}}=y_{\mathcal{Q}}v_{a}/\sqrt{2}$.

Performing a chiral $\mathrm{U}(1)$ transformation such that
$\mathcal{Q}_{L/R}\mapsto\mathrm{e}^{\pm ia/(2v_{a})}\mathcal{Q}_{L/R}$, the
mass term for $\mathcal{Q}$ can be made independent of the axion field phase.
This transformation adds an anomalous $G\widetilde{G}$ term to Eq. (1) as well
as an $F\widetilde{F}$ term, where $G$ and $F$ are the gluon and photon field
strength tensors, respectively, and the tilde denotes their duals. With the
electromagnetic (EM) and color anomaly contributions due to the
$\mathrm{U}(1)_{\text{PQ}}$ charged quarks labeled $E$ and $N$ respectively,
the coupling terms become

$\displaystyle\mathcal{L}$
$\displaystyle\supset\frac{N\alpha_{\text{s}}}{4\pi}\frac{a}{v_{a}}G\widetilde{G}+\frac{E\alpha_{\text{em}}}{4\pi}\frac{a}{v_{a}}F\widetilde{F}$
$\displaystyle=\frac{\alpha_{\text{s}}}{8\pi
f_{a}}aG\widetilde{G}+\frac{\alpha_{\text{em}}}{8\pi
f_{a}}\frac{E}{N}aF\widetilde{F}\,,$ (2)

where $f_{a}=v_{a}/(2N)$. The axion-photon coupling is thus parameterized by
the anomaly ratio $E/N$ alone.

More precisely, the mass and coupling to photons for QCD axion models are
given by [16, 17]

$\displaystyle m_{a}$ $\displaystyle=\frac{\chi_{0}^{2}}{f_{a}}=$5.69\pm
0.05\text{\,}\mathrm{\SIUnitSymbolMicro
eV}$\left(\frac{${10}^{12}\text{\,}\mathrm{GeV}$}{f_{a}}\right)\,,$ (3)
$\displaystyle g_{a\gamma\gamma}$
$\displaystyle=\frac{\alpha_{\text{em}}}{2\pi
f_{a}}\,C_{a\gamma\gamma}=\frac{\alpha_{\text{em}}}{2\pi
f_{a}}\left[\frac{E}{N}-C_{a\gamma\gamma}^{(0)}\right]$
$\displaystyle=\frac{\alpha_{\text{em}}}{2\pi
f_{a}}\left[\frac{E}{N}-($1.92\pm 0.04$)\right]\,.$ (4)

For some representation $r$ under which the heavy quark $\mathcal{Q}$ in the
KSVZ axion model transforms, the EM and color anomalies can be calculated as

$\displaystyle E$
$\displaystyle=\mathcal{X}\,d(\mathcal{C})\,\mathrm{tr}(q^{2})$
$\displaystyle=\mathcal{X}\,d(\mathcal{C})\,d(\mathcal{I})\left(\frac{d(\mathcal{I})^{2}-1}{12}+\mathcal{Y}^{2}\right)\,,$
(5a) $\displaystyle N$
$\displaystyle=\mathcal{X}\,d(\mathcal{I})\,T(\mathcal{C})\,,$ (5b)

where $d(\cdot)$ denotes the dimension of a representation,
$q=\mathcal{I}^{(3)}-\mathcal{Y}$ is the EM charge of $\mathcal{Q}$, and
$T(\mathcal{C})$ is the $\mathrm{SU}(3)_{\mathcal{C}}$ Dynkin index (see Ref.
[18]).

In KSVZ-type models, only $\mathcal{Q}$ is charged under the PQ symmetry
(apart from $\Phi$) and e.g. for $\mathcal{Q}\sim(3,1,0)$ we have
$N=\mathcal{X}/2$ and $E=3\mathcal{X}\,\mathrm{tr}(q^{2})$, using that
$T(3)=1/2$. In general, one finds for a single $\mathcal{Q}$ that

$\frac{E}{N}=6\,\mathrm{tr}(q^{2})=6q^{2}\,,$ (6)

where the last equality holds only when $\mathcal{Q}$ is a singlet under
$\mathrm{SU}(2)_{\mathcal{I}}$. This e.g. leads to the well-known result that
the original KSVZ model has $E/N=0$.

When considering models with multiple $\mathcal{Q}_{i}$, which have
representations $r_{i}$ and anomaly coefficients $E_{i}$ and $N_{i}$ given by
Eqs. (5a) and (5b), respectively, the overall anomaly ratio is simply

$\frac{E}{N}=\frac{\sum_{i}E_{i}}{\sum_{i}N_{i}}\,,$ (7)

where the index $i$ runs over the different quarks, labeled $i=1,\dots,n$.

Note that, when labeling a tuple of $\mathcal{Q}$s in a model, there exists a
“relabeling symmetry.” For example, assume that two $\mathcal{Q}$s with the
same $\mathrm{U}(1)_{\text{PQ}}$ charge respectively transform under
representations $r_{1}$ and $r_{2}$, denoted by $r_{1}\oplus r_{2}$. Then
there is an equivalency relation such that $r_{1}\oplus r_{2}\sim r_{2}\oplus
r_{1}$, in the sense that they trivially give the same anomaly ratio $E/N$.
Similarly, we can also consider combinations of representations with
“$\ominus$”, the symbol we use to denote $\mathcal{Q}$s with opposite
$\mathrm{U}(1)_{\text{PQ}}$ charges such that $r_{i}\ominus
r_{j}\Rightarrow\mathcal{X}_{i}=-\mathcal{X}_{j}$. Here we have e.g.
$r_{1}\oplus r_{2}\ominus r_{2}\sim r_{1}\ominus r_{2}\oplus r_{2}\sim
r_{2}\ominus\left(r_{1}\oplus r_{2}\right)$, as all three models trivially
give the same overall anomaly ratio.

The relabeling symmetry allows us to simplify the presentation of the catalog,
and we refer to a list of models where this symmetry has been accounted for as
“non-equivalent.” It may also play a role in the statistical interpretation of
the catalog: if not all $\mathcal{Q}$s are indistinguishable, the multiplicity
arising from the equivalency relation must be taken into account. We comment
on this further in Section V.1.

##  III Phenomenological selection criteria

Let us now review the various selection criteria for _preferred_ axion models,
most of which have already been proposed and discussed extensively in Refs.
[13, 14]. Here, we focus on the applicability in the pre- and post-
inflationary PQ symmetry breaking scenarios and observe that
$N_{\mathcal{Q}}~{}>~{}1$ allows for the existence of a new criterion related
to the axion’s ability to solve the strong CP problem.

###  III.1 Dark matter constraints

A natural requirement is to demand that axions do not produce more DM than the
observed amount, $\Omega_{\text{c}}h^{2}\lesssim 0.12$ [19]. For QCD axions
this results in an upper bound on $f_{a}$ and previous studies of _preferred_
axion models used $f_{a}<$5\text{\times}{10}^{11}\text{\,}\mathrm{GeV}$$ [13,
14], assuming a post-inflationary cosmology with realignment axion production.
Let us extend this discussion and make a few comments regarding the different
cosmological scenarios and their impact on the $f_{a}$ bound.

First, in the pre-inflationary PQ symmetry breaking scenario, the initial
misalignment angle of the axion field, denoted by $\theta_{\text{i}}$, is a
random variable. Since any topological defects are inflated away, realignment
production is the only relevant contribution and the limit on $f_{a}$ depends
on its “naturalness.” While this is not a uniquely defined concept, using the
usual assumption of uniformly distributed angles,
$\theta_{\text{i}}\sim\mathcal{U(-\pi,\pi)}$ the code developed in Ref. [20]
finds $f_{a}<$4\text{\times}{10}^{12}\text{\,}\mathrm{GeV}$$ for the 95%
credible region of posterior density.111Note that we used a prior of
$\log_{10}(f_{a}/$\mathrm{G}\mathrm{e}\mathrm{V}$)\sim\mathcal{U}(6,16)$,
which introduces some prior dependence, and also included QCD nuisance
parameters [20]. This limit on $f_{a}$ effectively relies on the naturalness
being encoded automatically in prior on $\theta_{\text{i}}$.

Second, when topological defects can be neglected in the post-inflationary
symmetry breaking, the relic axion density is determined by an average of
misalignment angles over many causally-disconnected patches. This corresponds
to the benchmark scenario of Refs. [13, 14]. Again using the code developed in
Ref. [20], we obtain $f_{a}<$2\text{\times}{10}^{11}\text{\,}\mathrm{GeV}$$
(at the 95% CL).

The third and last case is the post-inflationary scenario including a
significant contribution from topological defects i.e. cosmic strings and
domain walls (DWs). In fact, recent studies indicate that the production of
axions via topological defects dominates the vacuum realignment production
[21, 22]. For models with domain wall number $N_{\text{\tiny DW}}\equiv 2N=1$
(cf. Section III.5), the authors find that
$f_{a}\lesssim${10}^{10}\text{\,}\mathrm{GeV}$$, while models with
$N_{\text{\tiny DW}}>1$ reduce the value of $f_{a}$ by a factor
$\mathcal{O}(N_{\text{\tiny DW}})$ [22]. For the _preferred_ models considered
in this work, $N_{\text{\tiny DW}}\leq 28$ such that the bound might be
loosened to about
$f_{a}\lesssim$3\text{\times}{10}^{8}\text{\,}\mathrm{GeV}$$. It should be
noted that these results rely on extrapolating the outcome of numerical
simulations more than 60 orders of magnitude, and they hence are potentially
subject to large systematic uncertainties.

In summary, the upper limit on $f_{a}$, and hence the results presented in
what follows, very much depend on the cosmological scenario at hand. To
simplify the discussion, to avoid the potentially large uncertainties
mentioned above, and to better compare with previous work of Ref. [14], we
also adopt $f_{a}<$5\text{\times}{10}^{11}\text{\,}\mathrm{GeV}$$.

However, we stress again that a different choice of $f_{a}$ will affect the
number of _preferred_ models, as $f_{a}$ is one of the factors that determines
the value of $m_{\mathcal{Q}}$. This is because
$m_{\mathcal{Q}}=y_{\mathcal{Q}}\,v_{a}/\sqrt{2}=y_{\mathcal{Q}}\,N_{\text{\tiny
DW}}f_{a}/\sqrt{2}$, such that $f_{a}$ provides an upper bound on
$m_{\mathcal{Q}}$. Moreover, a universal bound on the $m_{\mathcal{Q}}$ (up to
the Yukawa couplings) requires that all $\mathcal{Q}$s are coupled to the
$\Phi$ field in the same way to get a single $v_{a}$ parameter. So long as the
coupling $y_{\mathcal{Q}}\sim\mathcal{O}(1)$ or lower, the upper bound on
$f_{a}=v_{a}/N_{\text{\tiny DW}}$ is indeed an upper limit to
$m_{\mathcal{Q}}$. Larger values of the coupling require fine-tuning of
parameters, and are hence deemed undesirable from a theoretical viewpoint. In
what follows, we choose
$m_{\mathcal{Q}}=$5\text{\times}{10}^{11}\text{\,}\mathrm{GeV}$$ as a
conservative value for all $\mathcal{Q}$ masses (see Sec. III.4 for more
details on the influence on Landau pole constraints).

Finally, note that the $\mathcal{Q}$s themselves contribute to the matter
content in the Universe, and we need to consider the possibility that their
abundance exceeds $\Omega_{\text{c}}h^{2}$. Since this issue can be avoided if
the lifetime of the quarks is short enough, we discuss this in the next
section.

###  III.2 Lifetimes

Other than the possibility that the $\mathcal{Q}$s’ abundances exceed
$\Omega_{\text{c}}h^{2}$, there also exist additional experimental and
observational constraints, which have already been discussed before [13, 14].

To avoid the DM constraints, we require the $\mathcal{Q}$s to decay into SM
particles with a reasonably low lifetime. Heavy quarks with
$m_{\mathcal{Q}}\gg$1\text{\,}\mathrm{TeV}$$ and lifetimes
$$0.01\text{\,}\mathrm{s}$<\tau_{\mathcal{Q}}<${10}^{12}\text{\,}\mathrm{s}$$
are severely constrained, as they would also affect Big Bang Nucleosynthesis
and observations of the Cosmic Microwave Background [23, 24]. Fermi-LAT
excludes
$${10}^{13}\text{\,}\mathrm{s}$<\tau_{\mathcal{Q}}<${10}^{26}\text{\,}\mathrm{s}$$,
thus excluding lifetimes greater than even the age of the Universe
($\sim${10}^{17}\text{\,}\mathrm{s}$$) [25]. As a result, for heavy quarks
($m_{\mathcal{Q}}\gg$1\text{\,}\mathrm{TeV}$$), only representations with
$\tau_{\mathcal{Q}}<${10}^{-2}\text{\,}\mathrm{s}$$ are considered to be a
part of the _preferred_ window. Lighter relics would be excluded from
experimental bounds e.g. at the LHC [26].

Such a constraint on the $\mathcal{Q}$ lifetime, when applied to the heavy
quark decay rate, translates to restrictions on the dimensionality of the
possible $\mathcal{Q}$ to SM fermion decay operators. With
$m_{\mathcal{Q}}\lesssim$5\text{\times}{10}^{11}\text{\,}\mathrm{GeV}$$, the
lifetime constraints in turn constrain operators to have dimensions
$d\leq~{}5$ [13, 14]. This implies a total of 20 possible representations for
$\mathcal{Q}$, all charged under $\mathrm{SU}(3)_{\mathcal{C}}$ and
$\mathrm{U}(1)_{\mathcal{Y}}$. The lifetime constraint has no further
consequence on cases with $N_{\mathcal{Q}}>1$ under the assumption that the
different $\mathcal{Q}_{i}$ do not interact among themselves or decay into
particles other than SM fermions.

As noted before [14], the lifetime constraints are typically not required in
the pre-inflationary PQ symmetry breaking scenario. This is because the
$\mathcal{Q}$s can get diluted by inflation, which prevents them from becoming
cosmologically dangerous relics after they freeze out. Without these
constraints, many more models with even higher-dimensional operators can
exist, and restricting ourselves to at most five-dimensional operators
therefore only becomes an assumption in this case.

###  III.3 Failure to solve the strong CP problem

This criterion is specific to models with $N_{\mathcal{Q}}>1$ that allow the
$\mathcal{Q}$s to have opposite $\mathrm{U}(1)_{\text{PQ}}$ charges. It is
clear from Eq. 5b that the addition of multiple heavy quarks can lead to a
smaller overall $N$ than the individual $N_{i}$, but only when one or more of
the quarks have a (relative) negative $\mathrm{U}(1)_{\text{PQ}}$ charge. In
some cases a total cancellation of the $N_{i}$ terms occurs ($N=0$). While
these models give rise to massless axion-like particles with a coupling to
photons governed by $E$, they do not solve the strong CP problem: as can be
seen from Eq. 2, $N=0$ means that there is no $G\widetilde{G}$ contribution in
the Lagrangian. Considering that the primary objective of QCD axion models is
to solve the strong CP problem, we propose that only models with $N\neq 0$
should be considered _preferred_.

###  III.4 Landau poles

The single most powerful criterion amongst the ones proposed by Refs. [13, 14]
in the context of this work comes from the observation that representations
with large $\mathcal{C}$, $\mathcal{I}$, or $\mathcal{Y}$ can induce Landau
poles (LPs) at energies well below the Planck mass. At an LP, the value of a
coupling mathematically tends to infinity, signaling a breakdown of the
theory. Since quantum gravity effects are only expected to appear at energies
near the Planck mass, a breakdown of the theory before that point can be
regarded as problematic or undesirable.

It has thus been proposed that _preferred_ models have LPs at energy scales
$\Lambda_{\text{LP}}\gtrsim${10}^{18}\text{\,}\mathrm{GeV}$$. From the 20
representations mentioned previously, only 15 fulfil this criterion [13, 14];
we refer to these as “LP-allowed” models and label them $r_{1}$ to $r_{15}$
(as per Table II in Ref. [13]).

The running of the couplings are computed at two-loop level with the
renormalization group equation [27, 28]

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}\alpha_{i}^{-1}$
$\displaystyle=-a_{i}-\frac{b_{ij}}{4\pi}\alpha_{j}\,\,,$ (8)

where

$\displaystyle a_{i}$
$\displaystyle=-\frac{11}{3}\,C_{2}(G_{i})+\frac{4}{3}\sum_{F}\kappa\,T(F_{i})+\frac{1}{3}\sum_{S}\eta\,T(S_{i})\,,$
(9a) $\displaystyle b_{ij}$
$\displaystyle=\left[-\frac{34}{3}\,\big{(}C_{2}(G_{i})\big{)}^{2}+\sum_{F}\left(4C_{2}(F_{i})+\frac{20}{3}C_{2}(G_{i})\right)\kappa\,T(F_{i})+\sum_{S}\left(4C_{2}(S_{i})+\frac{2}{3}C_{2}(G_{i})\right)\eta\,T(S_{i})\right]\delta_{ij}$
$\displaystyle+4\left(1-\delta_{ij}\right)\left[\sum_{F}\kappa\,C_{2}(F_{j})\,T(F_{i})+\sum_{S}\eta\,C_{2}(S_{j})\,T(S_{i})\right]\,,$
(9b)

with $i,j\in\\{1,2,3\\}$ for the three gauge groups,
$\alpha_{i}=g_{i}^{2}/4\pi$, $\mathrm{t}=\frac{1}{2\pi}\ln(\mu/m_{Z})$ for
energy scale $\mu$ and $Z$ boson mass $m_{Z}$, while $a_{i}$ and $b_{i}$ are
the one- and two-loop beta functions. $C_{2}$ and $T$ are the quadratic
Casimir and Dynkin indices of the corresponding gauge group, respectively, and
$F$ and $S$ denote fermionic and scalar fields. $G_{i}$ denotes the adjoint
representation of the gauge group, and $\kappa=\frac{1}{2},1$ for Weyl and
Dirac fermions, while $\eta=1$ for complex scalars.222The case of
$\eta=\frac{1}{2}$ for real scalars is not relevant for the present study.
Also note that the expression for $b_{ij}$ in Ref. [28] is slightly erroneous
since the second term applies only to the non-diagonal elements of $b_{ij}$,
as found when comparing with the SM beta functions in Ref. [27]. Adding
multiple $\mathcal{Q}$s to the theory increases the coefficients of beta
functions through the fermionic terms. As a consequence, the couplings diverge
faster i.e. induce LPs at lower energy scales, as has been anticipated before
[14].

Since the addition of more particles with a given representation into a gauge
theory only worsens the running of the corresponding gauge coupling, it is
possible to find the number of copies of a particle that can be included in
the theory before it induces an LP below ${10}^{18}\text{\,}\mathrm{GeV}$.
This drastically reduces the number of LP-allowed combinations possible.
Integrating all $\mathcal{Q}_{i}$ in at
$m_{\mathcal{Q}}=$5\text{\times}{10}^{11}\text{\,}\mathrm{GeV}$$, we find that
there are 59,066 non-equivalent combinations of $\mathcal{Q}_{i}$ from the
representations $r_{1},r_{2},\dots,r_{15}$ that do not induce LPs below
${10}^{18}\text{\,}\mathrm{GeV}$.

As the $\mathcal{Q}_{i}$ contribute to the beta functions above the energy
scale $m_{\mathcal{Q}}$, the running of the gauge coupling begins to deviate
from the SM only at this scale. Different values of $m_{\mathcal{Q}}$ are
bound to produce different results for the LPs; the lower $m_{\mathcal{Q}}$
is, the earlier an LP appears. As an example, consider
$m_{\mathcal{Q}}=${10}^{10}\text{\,}\mathrm{GeV}$$: for $N_{\mathcal{Q}}=3$,
we find that 888 models are _preferred_ (they have
$\Lambda_{\text{LP}}>${10}^{18}\text{\,}\mathrm{GeV}$$ and $N\neq 0$ as per
the discussion in III.3), compared to 1,442 models when
$m_{\mathcal{Q}}=$5\text{\times}{10}^{11}\text{\,}\mathrm{GeV}$$. Furthermore,
we use the same mass for all $\mathcal{Q}_{i}$ in the models, which may not be
the case in reality (due to different $y_{\mathcal{Q}_{i}}$ or e.g. in multi-
axion models). However, setting the masses to the highest possible value in
the _preferred_ window allows us to keep the number of disfavored models to a
minimum. Without further information on the values of $f_{a}$ and individual
$m_{\mathcal{Q}_{i}}$, excluding fewer models may be advantageous in the sense
of presenting a more inclusive $E/N$ catalog.

###  III.5 Other interesting model properties

Let us summarize a few other possible model properties, already discussed in
Refs. [13, 14].

Since the axion is the angular degree of freedom of the PQ scalar field, it
has a periodic potential and several degenerate vacua, given by the domain
wall number $N_{\text{\tiny DW}}=2N$. During PQ symmetry breaking, the axion
field can settle into any of these degenerate minima in different Hubble
patches, giving rise to domain walls. The energy density contained in such
topological defects can far exceed the energy density of the Universe [29] in
the post-inflationary PQ breaking scenario. However, in models with
$N_{\text{\tiny DW}}=1$, the string-domain wall configuration would be
unstable [30], which presents a possible solution and makes $N_{\text{\tiny
DW}}=1$ a desirable property of such models.

However, the DW problem can be avoided by allowing for a soft breaking of the
PQ symmetry [29]. Moreover, in a pre-inflationary PQ symmetry breaking
scenario, the patches and the topological defects are inflated away [31]. In
line with Refs. [13, 14], we therefore do not impose this criterion.

Among the 15 LP-allowed representations, only two have $N_{\text{\tiny
DW}}=1$. When all $\mathcal{Q}_{i}$ have the same $\mathrm{U}(1)_{\text{PQ}}$
charges, such a restriction would forbid any models with multiple heavy
quarks. With this in mind, a constraint on $N_{\text{\tiny DW}}$ is not used
to exclude $N_{\mathcal{Q}}>1$ models. In cases where the $\mathcal{Q}_{i}$
are permitted to have opposite $\mathrm{U}(1)_{\text{PQ}}$ charges, more
complicated models with $N_{\text{\tiny DW}}=1$ can be built by choosing the
$\mathcal{Q}_{i}$ such that $\sum_{i}N_{i}=1/2$. Even then, the number of such
models is few in comparison to the whole set of LP-allowed models.

Another intriguing property is the unification of the gauge couplings due to
the presence of the $\mathcal{Q}$s. The authors of Refs. [13, 14] note that
one of the 15 LP-allowed representations induces a significant improvement in
unification. While we do not investigate this further, we expect to find more
models that improve unification for higher $N_{\mathcal{Q}}$, which might be
an interesting topic for a future study.

##  IV Model catalog and anomaly ratio distributions

Table 1: Selected statistics for the complete set of models with $N_{\mathcal{Q}}\leq 9$. We include information about the $E/N$ ratios that give rise to the largest axion-photon coupling i.e. $\widehat{E/N}\equiv\mathrm{argmax}_{E/N}(|E/N-1.92|)$, photophobic models ($|E/N-1.92|<0.04$), and _preferred_ (LP-allowed and $N\neq 0$) models.  $N_{\mathcal{Q}}$ | Total #models | $\widehat{E/N}$ | LP-allowed | $N\neq 0$ | #_preferred_ | $\widehat{E/N}$ | photophobic
---|---|---|---|---|---|---|---
|  |  | fraction of total [%] |  | among _preferred_
1 | 20 | $\phantom{-00}44/3$ | 75.00 | 100.00 | 15 | $\phantom{-0}44/3$ | 0.00%
2 | 420 |  0$-184/3$ | 49.52 | 91.67 | 189 | $\phantom{-}122/3$ | 1.59%
3 | 5,740 | $\phantom{-0}368/3$ | 25.98 | 97.40 | 1,442 | $\phantom{-}170/3$ | 1.11%
4 | 61,810 |  0$-538/3$ | 11.60 | 97.37 | 6,905 | $-136/3$ | 1.29%
5 | 543,004 | $\phantom{-0}698/3$ | 4.42 | 98.13 | 23,198 | $-148/3$ | 1.27%
6 | 4,073,300 |  0$-928/3$ | 1.50 | 98.32 | 58,958 | $-160/3$ | 1.28%
7 | 26,762,340 | $-1108/3$ | 0.47 | 98.55 | 120,240 | $\phantom{-}164/3$ | 1.33%
8 | 157,233,175 | $\phantom{-}1292/3$ | 0.14 | 98.68 | 207,910 | $-166/3$ | 1.34%
9 | 838,553,320 | $-1312/3$ | 0.04 | 98.79 | 312,360 | $-142/3$ | 1.37%
Figure 1: Number of non-equivalent models with different properties as a
function of $N_{\mathcal{Q}}$. We show the number of all possible, additive,
LP-allowed, $N_{\text{\tiny DW}}=1$, $N=0$, photophobic ($|E/N-1.92|<0.04$)
models, as well as the number of different and unique (no other non-equivalent
model has the same $E/N$ value such that the underlying model is uniquely
identifiable) $E/N$ values.

Let us discuss a few key findings and properties of the model catalog created
in this work, which we summarize in Table 1 and Fig. 1.

To structure the discussion, we single out two subsets of the total model
space: one where all $\mathcal{Q}_{i}$ transform under the same representation
and one where the representations are arbitrary but the
$\mathrm{U}(1)_{\text{PQ}}$ charges of the quarks have the same sign (we call
these “additive models”).

###  IV.1 Subset I. Identical representations

First, consider the case where only representations of the form
$\bigoplus_{j=1}^{N_{\mathcal{Q}}}r_{i}$ with fixed $i\in[1,20]$ are allowed.
The number of possible models for a given $N_{\mathcal{Q}}$ is then simply
$N_{r}=20$, such that the total number of models up to and including some
$N_{\mathcal{Q}}$ is $N_{\text{tot}}=N_{r}\,N_{\mathcal{Q}}$.

Given that all quarks in such models have the same representation and
$\mathrm{U}(1)_{\text{PQ}}$ charge, only twelve discrete values of $E/N$ are
allowed when the LP criterion is taken into account [13]. However, the
relative distribution is determined by the effect of each representation on
the gauge group beta functions. We find that $\bigoplus_{j=1}^{28}r_{1}$ is
the only LP-allowed model for $N_{\mathcal{Q}}=28$ and that there are in total
79 _preferred_ models in this subset.

###  IV.2 Subset II. Allowing different additive representations

Next, consider the case where we can have arbitrary additive representations,
written in such a way that they respect the relabeling symmetry:
$\bigoplus_{i=1}^{20}\bigoplus_{j}^{n_{i}}r_{i}$, where
$\sum_{i}n_{i}=N_{\mathcal{Q}}$ with $n_{i}\geq 0$. The number of models in
this subset is

$\displaystyle N(N_{\mathcal{Q}})$
$\displaystyle=\binom{N_{\mathcal{Q}}+N_{r}-1}{N_{\mathcal{Q}}}\,,$ (10)
$\displaystyle N_{\text{tot}}$
$\displaystyle=\sum_{n}\binom{n+N_{r}-1}{n}=\binom{N_{\mathcal{Q}}+N_{r}}{N_{\mathcal{Q}}}\,.$
(11)

We find that, after applying the selection criteria, there are 59,066
_preferred_ models for $N_{\mathcal{Q}}\leq 28$. In particular, for
$N_{\mathcal{Q}}=28$, there are only nine LP-allowed models, none of which can
be extended by another quark while preserving the criterion. The highest
freedom in this subset is found for $N_{\mathcal{Q}}=10$, where 5,481 models
fall in the _preferred_ region.

Among these models, the smallest and largest anomaly ratios are 1/6 and 44/3
respectively, both of which come from $N_{\mathcal{Q}}=1$ models. The median
of the distribution of this set of models is $\mathrm{med}(E/N)\approx 1.87$,
indicating that $|C_{a\gamma\gamma}|\sim 0$ is a real possibility for a larger
fraction of the model space. Indeed, there are several models that have an
$E/N$ ratio close to the nominal value of the model-independent parameter
$C_{a\gamma\gamma}^{(0)}$. We define models as “photophobic” if their $E/N$
ratio is within one standard deviation of the nominal
$C_{a\gamma\gamma}^{(0)}$ value:

$\left|E/N-1.92\right|<0.04\,.$ (12)

We find that 3,255 models ($\approx 5.5\%$) among the 59,066 non-equivalent
models are photophobic. Considering all _preferred_ additive models up to
$N_{\mathcal{Q}}\leq 28$, there are 443 different $E/N$ values. Out of these,
28 are unique in the sense that they are uniquely identifiable since their
anomaly ratio $E/N$ is different from any other non-equivalent model.

###  IV.3 Complete set

Figure 2: Example histogram of the anomaly ratio $E/N$ for non-equivalent
$N_{\mathcal{Q}}=5$ models. Blue bars correspond to the “additive” subset and
red bars to the complete set of models i.e. also allowing for opposite
$\mathrm{U}(1)_{\text{PQ}}$ charges.

Finally, let us comment on the complete set of possible models where we may
also subtract representations, denoted by “$\ominus$.” Allowing
$\mathrm{U}(1)_{\text{PQ}}$ charges to have one of the two possible values for
each $\mathcal{Q}_{i}$, we open the window to a much wider range of possible
$E/N$ values. In particular, the anomaly ratio, and thus the axion-photon
coupling, can become negative (see Fig. 2) and, as mentioned before, the
solution to the strong CP problem can be spoilt in models with $N=0$.

For $n_{\oplus}+n_{\ominus}=N_{\mathcal{Q}}$, where $n_{\oplus}$ and
$n_{\ominus}$ are the number of $\mathcal{Q}$s with “positive” and “negative”
$\mathrm{U}(1)_{\text{PQ}}$ charges,333We remind the reader that “positive”
and “negative” are only relative concepts, in the sense that we consider two
models equivalent if the only difference between them is that the
$\mathrm{U}(1)_{\text{PQ}}$ charges of _all_ quarks get flipped going from one
to the other. respectively, the number of models with $n_{\oplus}>n_{\ominus}$
is simply

$\displaystyle
N(n_{\oplus},n_{\ominus})=\binom{n_{\oplus}+N_{r}-1}{n_{\oplus}}\,\binom{n_{\ominus}+N_{r}-1}{n_{\ominus}}\,.$
(13)

In the case where $n\equiv n_{\oplus}=n_{\ominus}$, accounting for the fact
that the anomaly ratio depends on the relative $\mathrm{U}(1)_{\text{PQ}}$
charges of the $\mathcal{Q}$s such that we have an equivalence of the type
$(r_{i}~{}\oplus~{}r_{j})~{}\ominus~{}(r_{k}~{}\oplus~{}r_{l})~{}\sim~{}(r_{k}~{}\oplus~{}r_{l})~{}\ominus~{}(r_{i}~{}\oplus~{}r_{j})$,
we also need to take care not to double-count models exhibiting this symmetry,
giving

$\displaystyle
N(n,n)=\frac{1}{2}\,\binom{n+N_{r}-1}{n}\left[\binom{n+N_{r}-1}{n}+1\right]\,.$
(14)

With this, we find that the number of models grows very fast as
$N_{\mathcal{Q}}$ increases. This also makes it computationally difficult to
compute and store all of the different combinations – let alone check the
criteria for _preferred_ models. We therefore restrict the complete analysis
in this case to $N_{\mathcal{Q}}\leq 9$.

The anomaly ratio distribution in the complete set exhibits a peak near zero,
and we expect the trend to continue even for larger $N_{\mathcal{Q}}$.
However, in general care should be taken when interpreting the “trends”
visible in Fig. 1. For example, the number of LP-allowed models will
eventually go down again as we move towards $N_{\mathcal{Q}}=28$, despite the
quickly growing total number of possible models. One may speculate that the
number of uniquely identifiable $E/N$ ratios could exhibit a similar behavior
as the number of LP-allowed models, while the number of different $E/N$ might
eventually saturate.

Allowing for opposite $\mathrm{U}(1)_{\text{PQ}}$ charges gives rise to models
with large axion-photon coupling; the largest and smallest values of $E/N$
found, $170/3$ and $-166/3$ respectively, give larger $|C_{a\gamma\gamma}|$
than what is possible in the previously discussed subsets. Note that the
$N_{\mathcal{Q}}=8$ model for $E/N=-166/3$ ($r_{2}\oplus r_{2}\oplus
r_{5}\oplus r_{6}\oplus r_{7}\ominus r_{1}\ominus r_{9}\ominus r_{9}$) was not
reported in Refs. [13, 14] as giving the highest possible
$|C_{a\gamma\gamma}|$; instead the authors indicated that $E/N=170/3$ led to
the largest absolute value of the coupling. We find that among the complete
set of 5,753,012 _preferred_ models, there are 81,502 photophobic models and
820 different anomaly ratios, with 79 out of those also being from uniquely
identifiable models.

##  V Impact on axion searches

In this section, we discuss possible statistical interpretations of the
hadronic axion model catalog and show the impact of these on the mass-coupling
parameter space.

###  V.1 On constructing E/N prior distributions

The catalog of KSVZ models – even after applying the selection criteria – is
but a list of _possible_ models. It does not inherently contain information
about how _probable_ each model is. The model with $E/N=-166/3$ gives the
largest $|C_{a\gamma\gamma}|\approx 57$, which will place an upper bound on
the axion-photon coupling and delimit the upper end of the KSVZ axion band. On
the other end, complete decoupling with photons ($C_{a\gamma\gamma}\approx 0$)
is also possible within the theoretical errors. Since any of the models might
be realized in Nature, perhaps due to a deeper underlying reason that is not
obvious at present, one might be satisfied with this picture.

However, the boundaries of the band are extreme cases and do not take into
account where the bulk of possible models can be found. For example, defining
a desired target sensitivity for an experiment becomes non-trivial in the face
of $C_{a\gamma\gamma}$ potentially being extremely close to zero. We propose
instead that covering a certain fraction of all possible models or
constructing a prior volume might be more meaningful ways to define such a
target.

To directly interpret an $E/N$ histogram as a distribution implicitly makes
the assumption that each model is equally likely to be realized in Nature.
While this interpretation might be considered “fair,” one could argue that
models with many $\mathcal{Q}$s are more “contrived” and consequently
introduce a weighting factor that penalizes models with $N_{\mathcal{Q}}\gg
1$. This could be achieved with e.g. exponential suppression via a weighting
factor $\propto\mathrm{e}^{-N_{\mathcal{Q}}}$, or $\propto
2^{-N_{\mathcal{Q}}}$. Another option could be to choose models that are
minimal extensions ($N_{\mathcal{Q}}=1$) or similar to the family structure of
the SM ($N_{\mathcal{Q}}=3$ or e.g. a weighting $\propto
3^{N_{\mathcal{Q}}}/N_{\mathcal{Q}}!$).

Such consideration are aligned with the Bayesian interpretation of statistics,
and will probably meet criticism for this reason. However, as pointed out in
Ref. [20], at least in the pre-inflationary PQ symmetry breaking scenario,
which is fundamentally probabilistic in nature, the Bayesian approach is well
motivated. Furthermore, Ref. [20] also proposed that the discrete nature of
KSVZ models should be reflected in the prior choice of $E/N$. Such a
physically-motivated prior should further reflect the combinatorics of KSVZ
model building by including the multiplicity of $E/N$ ratios. As mentioned at
the end of Section II, this multiplicity also depends on whether or not the
$\mathcal{Q}_{i}$ are distinguishable by e.g. having different masses.

Figure 3: Anomaly ratio distributions for all _preferred_ additive KSVZ
models, using different weightings. For equal weighting, we show the
underlying histogram (blue shading) and a smooth Gaussian kernel density
estimate of the distribution (blue line), while for others we only show the
latter for simplicity.

With this in mind, we show different statistical interpretations of the
anomaly ratio in Fig. 3. For visualization purposes, we show kernel density
estimates of the distributions for different weighting factors mentioned
above, while reminding the reader that the underlying histograms and
distributions are actually discrete and not continuous.

From Fig. 3 it becomes clear that the different weightings can change the
width of the distribution, introducing a prior dependence in an analysis.
However, the modes of the distributions remain around $E/N\sim 2$, which means
that a partial cancellation of the axion-photon coupling $C_{a\gamma\gamma}$
is typically possible, as already observed in Fig. 2.

###  V.2 Experimental constraints on _preferred_ KSVZ axion models

Figure 4: The KSVZ axion band as defined by the 68% and 95% central regions of
$|C_{a\gamma\gamma}|=|E/N-C_{a\gamma\gamma}^{(0)}|$, drawing $E/N$ from a
distribution of all _preferred_ KSVZ axion models (each representation assumed
to be equally probable). The grey line marks the highest possible absolute
value of the coupling ($E/N=-166/3$), while the black line indicates the
classical KSVZ model ($E/N=0$). For context, we show various present (shaded
regions) and future (dashed lines) haloscope (blue) and helioscope (red)
limits and forecasts [32] as well as bounds from hot dark matter [33], energy
loss in SN1987A [34], and recent string simulations [22].

Of course, a possible partial cancellation of the axion-photon coupling has
consequences on the various astrophysical, cosmological, and laboratory
searches (see e.g. Ref. [10]) for axions. The most powerful analyses combine
the results of different experiments to place joint limits on the properties
of different types of axions (e.g. Refs. [35, 36, 20]).

To investigate this further, consider e.g. a prior on $E/N$ where all
_preferred_ (LP-allowed models with $d\leq 5$ operators and $N\neq 0$), non-
equivalent KSVZ models are considered equally probable.444Recall that the
$d\leq 5$ condition is due to the lifetime constraints (see Section III.2) in
the post-inflationary scenario, while it is only an assumption for the pre-
inflationary case (potentially reasonable for as being a minimal extension of
the SM). We can then generate samples for
$C_{a\gamma\gamma}=E/N-C_{a\gamma\gamma}^{(0)}$, where $E/N$ is drawn from its
discrete distribution and
$C_{a\gamma\gamma}^{(0)}\sim\mathcal{N}(1.92,\,0.04)$ i.e. follows a normal
distribution with mean 1.92 and standard deviation 0.04.

We find that the central 68% region of the ensuing distribution corresponds to
$|C_{a\gamma\gamma}|\in[0.39,5.22]$, while the 95% region is
$|C_{a\gamma\gamma}|\in[0.06,17.30]$. The corresponding model bands in the
mass-coupling plane are shown in Fig. 4, and the bulk of these models can be
constrained by present and future experiments. In fact, while complete
cancellation of $C_{a\gamma\gamma}$ is possible within the theoretical
uncertainty for some $E/N$ values, we find that the bulk of models is at worst
somewhat suppressed. This is very encouraging for experimental searches.

Had we only considered additive models, the 95% region would be
$|C_{a\gamma\gamma}|\in[0.02,1.67]$, such that the upper end of the band would
be lower than the traditional KSVZ model with $E/N=0$. This can be readily
understood from the $E/N$ distributions in Fig. 3, whose mode is typically
close to $C_{a\gamma\gamma}^{(0)}$ such that the value of
$|g_{a\gamma\gamma}|$ is lower than what would be expected for
$|C_{a\gamma\gamma}|\sim\mathcal{O}(1)$. In this case, the planned future
experiments would _not_ be able to probe large parts of the band, indicating
that the choice of prior – even if physically-motivated – can induce a
noticeable impact on the results.

##  VI Summary and conclusions

We provide a catalog of all hadronic, or KSVZ, axion models with
$N_{\mathcal{Q}}\leq 9$, featuring 1,027,233,129 non-equivalent models in
total. When we apply the selection criteria for _preferred_ models, we find a
limit of $N_{\mathcal{Q}}\leq 28$ and that only 5,753,012 non-equivalent
models with 820 different $E/N$ values exist (59,066 non-equivalent models
with 443 different $E/N$ values for additive representations). While relaxing
existing or adding new criteria can increase or reduce these numbers, we
generically expect that the Landau pole (LP) criterion will be a powerful tool
to limit the number of possible models – even with modified constraints or in
other axion models. This is similar to the Standard Model, where the number of
families can also be restricted by demanding that LPs do not appear below some
energy scale. We further propose that only models with QCD anomaly $N\neq 0$
be considered _preferred_.

Our model catalog can be a useful, searchable database for researchers wishing
to study the KSVZ axion model space. It allows to e.g. make statements about
what fraction of possible models a given experiment is sensitive to. We made
catalogs, histograms, and example Python scripts available on the Zenodo
platform for this purpose [15].

Some models in the catalog might be considered “contrived” as they add many
new particles to the theory. Of course, in case of a discovery or if any other
appealing reason for a seemingly more complicated models is put forward, this
perception might change. In absence of such reasons, the $E/N$ values may be
interpreted as statistical distributions, which encode assumptions about the
probability of the different models. We generally outlined how prior
distributions can be constructed from the catalog and gave concrete examples
of such choices.

For the specific choice of equally probable _preferred_ models, we consider
the consequences for axion searches and the definition of the KSVZ axion band.
Here we suggest that the latter may be defined as the central 95% region of
all models, taking into account uncertainties from the model-independent
contribution to the axion-photon coupling. If only “additive models” are
considered, the bulk of the _preferred_ models can unfortunately not be probed
by current or future experiments since the anomaly ratio distributions in this
case tend to peak around $E/N\sim 2$. In general, using the discrete $E/N$
distributions improves on unphysical prior choices considered in the past
(e.g. Ref. [20]).

Even when ignoring the statistical perspective, it is useful for axion
searches to know that the _preferred_ models only admit 820 different $E/N$
values. In case of an axion detection, one may therefore test these discrete
models against each other to see which models are most compatible with the
detected signal. One could further test them against a generic axion-like
particle or other QCD axion models. In an ideal scenario, this might even
allow an experiment to infer the underlying high-energy structure of a model,
which highlights the known property of axion models to connect high-energy
physics to low-energy observables.

In summary, the powerful LP criterion restricts the number of KSVZ models to a
finite value. In that sense, the catalog presented here is a complete list of
all _preferred_ KSVZ models, which may be used as input for axion searches and
forecasts. Since KSVZ models could e.g. be extended by also considering
multiple complex scalar fields or feature more complex couplings to the SM,
and since there are other kinds of QCD axion models such as the DFSZ-type
models, this work presents another step forward in mapping the landscape of
all phenomenologically interesting axion models.

###### Acknowledgements.

We thank Maximilian Berbig, Joerg Jaeckel, and David ‘Doddy’ J. E. Marsh for
useful comments and discussions. This paper is based on results from VP’s
ongoing M.Sc. project. SH is supported by the Alexander von Humboldt
Foundation and the German Federal Ministry of Education and Research. We used
the Scientific Computing Cluster at GWDG, the joint data center of Max Planck
Society for the Advancement of Science (MPG) and the University of Göttingen.

## References

  * [1] S. Weinberg, A new light boson?, _Physical Review Letters_ 40 (1978) 223.
  * [2] F. Wilczek, Problem of strong P and T invariance in the presence of instantons, _Physical Review Letters_ 40 (1978) 279.
  * [3] R.D. Peccei and H.R. Quinn, CP conservation in the presence of pseudoparticles, _Physical Review Letters_ 38 (1977) 1440.
  * [4] R.D. Peccei and H.R. Quinn, Constraints imposed by CP conservation in the presence of pseudoparticles, _Phys. Rev. D_ 16 (1977) 1791.
  * [5] J. Preskill, M.B. Wise and F. Wilczek, Cosmology of the invisible axion, _Physics Letters B_ 120 (1983) 127.
  * [6] L.F. Abbott and P. Sikivie, A cosmological bound on the invisible axion, _Physics Letters B_ 120 (1983) 133.
  * [7] M. Dine and W. Fischler, The not-so-harmless axion, _Physics Letters B_ 120 (1983) 137.
  * [8] M.S. Turner, Coherent scalar-field oscillations in an expanding universe, _Phys. Rev. D_ 28 (1983) 1243.
  * [9] M.S. Turner, Cosmic and local mass density of “invisible” axions, _Physical Review D_ 33 (1986) 889.
  * [10] I.G. Irastorza and J. Redondo, New experimental approaches in the search for axion-like particles, _Progress in Particle and Nuclear Physics_ 102 (2018) 89 [1801.08127].
  * [11] J.E. Kim, Weak-interaction singlet and strong CP invariance, _Phys. Rev. Lett._ 43 (1979) 103.
  * [12] M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, Can confinement ensure natural CP invariance of strong interactions?, _Nuclear Physics B_ 166 (1980) 493.
  * [13] L. Di Luzio, F. Mescia and E. Nardi, Redefining the Axion Window, _Physical Review Letters_ 118 (2017) 031801 [1610.07593].
  * [14] L. Di Luzio, F. Mescia and E. Nardi, Window for preferred axion models, _Phys. Rev. D_ 96 (2017) 075003 [1705.05370].
  * [15] V. Plakkot and S. Hoof, “Model catalogues and histograms of KSVZ axion models with multiple heavy quarks.” Published on Zenodo, 2021.  DOI: 10.5281/zenodo.5091707.
  * [16] G.G. di Cortona, E. Hardy, J.P. Vega and G. Villadoro, The QCD axion, precisely, _JHEP_ 1 (2016) 34 [1511.02867].
  * [17] M. Gorghetto and G. Villadoro, Topological susceptibility and QCD axion mass: QED and NNLO corrections, _JHEP_ 2019 (2019) 33 [1812.01008].
  * [18] R. Slansky, Group theory for unified model building, _Phys. Rep._ 79 (1981) 1.
  * [19] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi et al., Planck 2018 results. VI. Cosmological parameters, _A &A_ 641 (2020) A6 [1807.06209].
  * [20] S. Hoof, F. Kahlhoefer, P. Scott, C. Weniger and M. White, Axion global fits with Peccei-Quinn symmetry breaking before inflation using GAMBIT, _JHEP_ 3 (2019) 191 [1810.07192].
  * [21] M. Gorghetto, E. Hardy and G. Villadoro, Axions from strings: the attractive solution, _JHEP_ 7 (2018) 151 [1806.04677].
  * [22] M. Gorghetto, E. Hardy and G. Villadoro, More Axions from Strings, _SciPost Physics_ 10 (2021) 050 [2007.04990].
  * [23] M. Kawasaki, K. Kohri and T. Moroi, Big-Bang nucleosynthesis and hadronic decay of long-lived massive particles, _Phys. Rev. D_ 71 (2005) 083502 [astro-ph/0408426].
  * [24] J. Chluba, Distinguishing different scenarios of early energy release with spectral distortions of the cosmic microwave background, _Mon. Not. Roy. Astron. Soc._ 436 (2013) 2232 [1304.6121].
  * [25] Fermi-LAT collaboration, Fermi LAT Search for Dark Matter in Gamma-ray Lines and the Inclusive Photon Spectrum, _Phys. Rev. D_ 86 (2012) 022002 [1205.2739].
  * [26] S. Jäger, S. Kvedaraitė, G. Perez and I. Savoray, Bounds and prospects for stable multiply charged particles at the LHC, _JHEP_ 04 (2019) 041 [1812.03182].
  * [27] M.E. Machacek and M.T. Vaughn, Two Loop Renormalization Group Equations in a General Quantum Field Theory. 1. Wave Function Renormalization, _Nucl. Phys. B_ 222 (1983) 83.
  * [28] L. Di Luzio, R. Gröber, J.F. Kamenik and M. Nardecchia, Accidental matter at the LHC, _JHEP_ 07 (2015) 074 [1504.00359].
  * [29] P. Sikivie, Axions, Domain Walls, and the Early Universe, _Phys. Rev. Lett._ 48 (1982) 1156.
  * [30] S.M. Barr, K. Choi and J.E. Kim, Axion Cosmology in Superstring Models, _Nucl. Phys. B_ 283 (1987) 591.
  * [31] J.E. Kim, Light pseudoscalars, particle physics and cosmology., _Phys. Rep._ 150 (1987) 1.
  * [32] C. O’Hare, “cajohare/AxionLimits: AxionLimits.” Published on Zenodo, 2020.  DOI: 10.5281/zenodo.3932430.
  * [33] W. Giaré, E.D. Valentino, A. Melchiorri and O. Mena, New cosmological bounds on hot relics: Axions & Neutrinos, _MNRAS_ (2021)  [2011.14704].
  * [34] P. Carenza, T. Fischer, M. Giannotti, G. Guo, G. Martínez-Pinedo and A. Mirizzi, Improved axion emissivity from a supernova via nucleon-nucleon bremsstrahlung, _JCAP_ 2019 (2019) 016 [1906.11844].
  * [35] M. Giannotti, I.G. Irastorza, J. Redondo, A. Ringwald and K. Saikawa, Stellar recipes for axion hunters, _JCAP_ 10 (2017) 010 [1708.02111].
  * [36] L. Visinelli and S. Vagnozzi, Cosmological window onto the string axiverse and the supersymmetry breaking scale, _Phys. Rev. D_ 99 (2019) 063517 [1809.06382].


# A simple iterative map forecast of the COVID-19 pandemic

André E. Botha  Department of Physics, University of South Africa, Private Bag
X6, Florida 1710, South Africa<EMAIL_ADDRESS>Wynand Dednam
Departamento de Fisica Aplicada, Universidad de Alicante, San Vicente del
Raspeig, 03690 Alicante, Spain

###### Abstract

We develop a simple 3-dimensional iterative map model to forecast the global
spread of the coronavirus disease. Our model contains at most two fitting
parameters, which we determine from the data supplied by the world health
organisation for the total number of cases and new cases each day. We find
that our model provides a surprisingly good fit to the currently-available
data, which exhibits a cross-over from exponential to power-law growth, as
lock-down measures begin to take effect. Before these measures, our model
predicts exponential growth from day 30 to 69, starting from the date on which
the world health organisation provided the first ‘Situation report’ (21
January 2020 – day 1). Based on this initial data the disease may be expected
to infect approximately 23% of the global population, i.e. about 1.76 billion
people, taking approximately 83 million lives. Under this scenario, the global
number of new cases is predicted to peak on day 133 (about the middle of May
2020), with an estimated 60 million new cases per day. If current lock-down
measures can be maintained, our model predicts power law growth from day 69
onward. Such growth is comparatively slow and would have to continue for
several decades before a sufficient number of people (at least 23% of the
global population) have developed immunity to the disease through being
infected. Lock-down measures appear to be very effective in postponing the
unimaginably large peak in the daily number of new cases that would occur in
the absence of any interventions. However, should these measure be relaxed,
the spread of the disease will most likely revert back to its original
exponential growth pattern. As such, the duration and severity of the lock-
down measures should be carefully timed against their potentially devastating
impact on the world economy.

## Introduction

On 11 March 2020, the world health organisation (WHO) characterised the 2019
outbreak of coronavirus disease (COVID-19) as a pandemic, referring to its
prevalence throughout the whole world[1]. The outbreak started as a pneumonia
of an unknown cause, which was first detected in the city of Wuhan, China. It
was reported as such to the WHO on the 31st December 2019, and has since
reached epidemic proportions within China, where it has infected more than 80
000 citizens, to date. During the first six weeks of 2020 the disease spread
to more than 140 other countries, creating wide-spread political and economic
turmoil, due to unprecedented levels of spread and severity.

The rapid spread of COVID-19 is fuelled by the fact that the majority of
infected people do not experience severe symptoms, thus making it more likely
for them to remain mobile, and hence to infect others[2]. At the same time the
disease can be lethal to some members of the population, having a globally
averaged fatality ratio of about $6\%$, at present (16 April 2020).
Furthermore the new virus seems capable of surviving for unusually long
periods on plastic and metal surfaces, which are both, frequently encountered
in everyday life. Detailed analysis of the virus also revealed that its outer
surface consists of club-like “spikes” that are about four times more
effective at establishing the infection than in the closely related
coronavirus that causes severe acute respiratory syndrome (SARS) in
2002-2003[3]. It is most likely this particular combination of traits that has
made the COVID-19 outbreak one of the largest in recorded history.

While there are a number of models available for the global spread of
infectious diseases[4], some even containing very sophisticated traffic
layers[5], relatively few researchers are making use of simpler models that
can provide the big picture without being difficult to interpret
unambiguously. In the latter category of relatively simple models we could
find only a discrete epidemic model for SARS [6], and more recently, a
comparison of the logistic growth and susceptible-infected-recovered (SIR)
models for COVID-19[7].

In our present work we develop a simple discrete 3-dimensional iterative map
model, which shares some similarities with the classic SIR model. We show that
our model can fit both the initial exponential growth in the number of cases,
before lock-down measures began to take effect (about day 69), as well as
subsequent power-law growth. In our view the current interventions are
necessary to prevent the unimaginably fast (exponential) spread of the disease
until some solution (perhaps in the form of a vaccine) can be developed.

## Results

### Simple model

As a simple exponential growth model for the global data we initially develop
a 3-dimensional iterative map model given by

$\displaystyle x_{i+1}$ $\displaystyle=$ $\displaystyle x_{i}+\alpha
y_{i}\left(z_{0}-x_{i}\right)/z_{0}$ $\displaystyle y_{i+1}$ $\displaystyle=$
$\displaystyle\alpha y_{i}\left(z_{0}-x_{i}\right)/z_{0}$ (1) $\displaystyle
z_{i+1}$ $\displaystyle=$ $\displaystyle z_{i}-cy_{i},$

where $x_{i}$ is the total number of confirmed cases, $y_{i}$ is the number of
new cases and $z_{i}$ is the global population, on any given day $i$. We
denote the only fitting parameter by $\alpha$, while $c$ is a fixed parameter
equal to the fraction of people who have died from the disease. According to
the latest available global data from the WHO (see the last row of Table 1 in
Methods ), $c=0.06343$.

We briefly describe the physical content of Eqs. (1). The first equation
simply updates the total number of cases by setting it equal to the previous
total number of cases, plus the number of new cases. Here the factor of
$1/z_{0}$ has been introduced for convenience, to ensure that the
proportionality constant $\alpha$ remains close to unity. In the second
equation see that the number of new cases is assumed to be proportional to the
previous number of new cases multiplied by the previous number of susceptible
people. For simplicity we assume here the the number of susceptible people are
all those who have not yet had the disease. The third equation is merely to
keep track of the global population, by subtracting the estimated number of
people who have died each day, based on the fraction $c$.

We have found that our simple model provides a very good fit to the global
data between the 29th and 70th days of the pandemic. By using
Levenberg–Marquardt (least squares) optimisation[8] (see Methods), we find
$\alpha=1.1393$, for the initial condition $x_{30}=75194.3$ $y_{30}=349.477$
and $z_{30}=z_{0}=7.7000\times 10^{9}$.

Figure 1 shows a comparison of the data with the model, as well as a forecast
made up to the 200th day.

Figure 1: A comparison of the data with Eqs. (1). (a) Close up view of the
extent to which the model fits the global data for the total number of
confirmed cases $x_{i}$ and new cases $y_{i}$, so far. (b) A forecast of how
the disease is likely to spread globally, up to the 200th day, i.e. 7 August
2020. As the inset shows, the peak number of new infections is predicted to
occur on day 133 (19 May 2020), when about 60 million new cases may be
expected.

As we see in Figure 1(a), the model provides a good fit to the data between
the 29th and 70th days, having a correlation coefficient (R-squared value) of
0.99432 or higher (see Table 1). The forecast made in Figure 1(b)
(corresponding to the 69th ’Last day’ in Table 1), predicts that approximately
a quarter of the worlds population, i.e. $\approx 1.76/7.7=0.23$, would have
had COVID-19 by the 200th day. The peak of the pandemic is expected to occur
on day 133, when about 60 million daily new cases can be expected. We also
predict that by the beginning of August 2020, hardly any new cases should
occur; however, the total number of lives lost by then could be as high as 83
million, i.e. a number which is almost the same as the global population
increase for the year 2019[9].

Last day | $\alpha$ | $x_{200}\times 10^{9}$ | $\max\\{y\\}\times 10^{6}$ | Day of max$\\{y\\}$ | Deaths$\times 10^{6}$ | $c$ | R-squared
---|---|---|---|---|---|---|---
59 | 1.12640 | 1.628 | 50.487 | 140 | 68.087 | 0.04183 | 0.99880
60 | 1.12874 | 1.653 | 52.218 | 138 | 69.506 | 0.04204 | 0.99833
61 | 1.13131 | 1.682 | 54.179 | 137 | 70.702 | 0.04203 | 0.99684
62 | 1.13256 | 1.693 | 55.054 | 136 | 74.092 | 0.04376 | 0.99615
63 | 1.13483 | 1.719 | 56.817 | 135 | 74.895 | 0.04358 | 0.99478
64 | 1.13652 | 1.737 | 58.143 | 134 | 75.637 | 0.04354 | 0.99396
65 | 1.13751 | 1.746 | 58.804 | 134 | 77.747 | 0.04452 | 0.99347
66 | 1.13861 | 1.758 | 59.711 | 133 | 79.136 | 0.04503 | 0.99361
67 | 1.13849 | 1.755 | 59.571 | 133 | 80.423 | 0.04583 | 0.99318
68 | 1.13923 | 1.762 | 60.100 | 133 | 81.663 | 0.04634 | 0.99402
69 | 1.13930 | 1.761 | 60.097 | 133 | 83.118 | 0.04719 | 0.99432
70 | 1.13799 | 1.746 | 59.065 | 133 | 83.399 | 0.04776 | 0.99280
71 | 1.13589 | 1.722 | 57.375 | 134 | 83.498 | 0.04848 | 0.98965
72 | 1.13505 | 1.712 | 56.673 | 135 | 84.377 | 0.04929 | 0.98965
73 | 1.13362 | 1.694 | 55.490 | 135 | 86.020 | 0.05078 | 0.98828
74 | 1.13197 | 1.674 | 54.213 | 136 | 86.655 | 0.05175 | 0.98633
75 | 1.13021 | 1.651 | 52.778 | 137 | 89.481 | 0.05419 | 0.98396
76 | 1.12831 | 1.629 | 51.307 | 138 | 90.189 | 0.05538 | 0.98096
77 | 1.12565 | 1.599 | 49.314 | 140 | 89.228 | 0.05582 | 0.97413
78 | 1.12189 | 1.556 | 46.592 | 142 | 88.258 | 0.05674 | 0.95830
79 | 1.11901 | 1.521 | 44.480 | 144 | 89.030 | 0.05855 | 0.94451
Table 1:  Variation in the predictions made by the model using only data up to
and including the ’Last day’, as indicated in the first column. The predicted
total number of ’Deaths’ indicated in column 6 is equal to $z_{30}-z_{200}$,
i.e. the difference between the initial population and the remaining
population on day 200.

In Table 1 we see that the fitting parameter $\alpha$, and hence the
predictions made by the model, do change somewhat as more of the available
data is used in the fitting procedure. To see the variation in $\alpha$ more
clearly we have plotted the first and second columns of Table 1 in Figure 2.
It shows that, as more data is used, there seems to be a general upward trend
in $\alpha$, until the 69th ’Last day’. Beyond this our ’Simple model’ no
longer provides a good fit to the data, as is evident from the rapidly
declining R-squared values that are given in the far right column of the
Table.

In Figure 2 we also plot (blue solid line) the mean value
$\bar{\alpha}=1.1389$ over ’Last day’ 66-69. The value of $\alpha$ is
essentially constant over these four days.

Figure 2: Variation of the fitting parameter $\alpha$ as more and more of the
available data is used in the fitting procedure. We see that the value of
$\alpha$ reaches its maximum on ’Last day’ 69. The average value of $\alpha$
over ’Last day’ 66-69 is given by $\bar{\alpha}=1.1389$.

### Transition from exponential to power-law growth

China is currently the only country in which the rapid further spread of the
virus appears to be under control. An investigation of the data for China
observed that the total number of infections as well as the number of
recoveries and the number of deaths followed a power law growth, rather than
growing exponentially[10, 11]. A wider study for 25 different countries found
that other countries also showed scale free power law growth behaviour with
different country-specific scaling exponents[12]. Subsequent in-depth analysis
of more available data of country-specific growth behaviours revealed,
however, that different countries exhibit generally different growth
behaviours[13]. It was found that while indeed some countries such as for
example the Netherlands or Norway exhibit a power law behaviour over nearly
the whole range of the available data, other countries such as for example the
US or France display long periods of very few infections followed by a massive
exponential surge but subsequently quickly reduce the growth rate again: a
behaviour that is clearly not describable by a power law. In general, however,
the spread of the disease within individual countries initially appears to be
exponential, followed by power law growth once national lock-down measures
take effect. We may therefore expect a similar trend in the global data, which
is more or less what we observe from day 70 onward. Note that day 70
corresponds to the 30th of March 2020, which is two to three weeks after many
of the hardest hit European countries implemented their national lock-downs:
Italy (10 March), Spain (14 March), France (17 March), etc.

One can rationalise the transition from exponential to power law growth as
follows. Any person living his/her life in the city has on any given day a
large number of contacts with people and infectious objects. Of these
interactions, relatively few can be classified as relations that are of
particular meaning to that individual, such as contact with a spouse,
children, friends, family relations, coworkers, etc. Thus most contacts are
completely random interactions, such as through public transport, contacts
with clients, public parks, shops, entertainment facilities, etc. For the
spread of the virus it is irrelevant whether the susceptible person’s
interaction with another individual is meaningful or random, as long as the
virus can be transmitted.

Since random interactions are generally never with the same people or
infectious objects, they generally give rise to an exponential spread of the
disease. The introduction of governmental measures, such as a ’lock-down’,
severely reduce random interactions whilst leaving more open the possibility
of relationship interactions. Human societies (friendships, relations,
collaborations, social networks etc.) have been found to be organised in small
world or scale-free networks[14, 15] exhibiting power law behaviours. Given
that the further spreading of the epidemic is now effectively being
restricted, on a global scale, to only those relation networks, it follows
that the global spread pattern should also change from exponential to a power
law growth, as has been shown in countries such as China, the Netherlands,
Norway and others.

In light of this expectation we have expanded our original simple model to
reflect the change of the growth behaviour after world wide governmentally
instituted lock-down measures have started to take effect, i.e.on about day
69, according to the global data.

### Improved model

As described in Methods we have derived an improved model which can take into
account both exponential and power law growth. This model is given by

$\displaystyle x^{\prime}_{i+1}$ $\displaystyle=$ $\displaystyle
x_{i}^{\prime}+y_{i}^{\prime}\exp\left(b\left[\frac{a}{y_{i}^{\prime}+a}\right]^{\frac{1}{b}}\right)\left(z_{i}^{\prime}-x_{i}^{\prime}\right)/z_{0}$
$\displaystyle y^{\prime}_{i+1}$ $\displaystyle=$ $\displaystyle
y^{\prime}_{i}\exp\left(b\left[\frac{a}{y^{\prime}_{i}+a}\right]^{\frac{1}{b}}\right)\left(z^{\prime}_{i}-x^{\prime}_{i}\right)/z_{0}$
(2) $\displaystyle z^{\prime}_{i+1}$ $\displaystyle=$ $\displaystyle
z^{\prime}_{i}-cy^{\prime}_{i},$

where there are now two fitting parameters, $a$ and $b$.

Equation (2) is more general than (1). In the limit,
$a>>\max\\{y^{\prime}_{i}\\}$, it reduces to the same form as equation (1),
since the exponential factor essentially becomes independent of
$y^{\prime}_{i}$, i.e. $a/(y^{\prime}_{i}+a)\xrightarrow[\infty]{a}1$ and, by
comparison to equation (1), $\alpha=\exp\left(b\right)$.

The main advantage of equation (2) is that it can take into account the
observed cross-over from exponential growth during the initial phase of the
pandemic, to a power law growth, now that measures to curb the spread are
starting to take effect on a global scale.

Since our ’Improved model’ can take into account the effects of the
interventions put in place to curb the global spread of the virus, we were
expecting it to provide a more optimistic forecast than that of our ’Simple
model’. However, although our ’Improved model’ does predict a vast reduction
in the daily number of new cases, it also predicts a prolonged spread of the
disease, as long as the interventions are strictly enforced. This is shown by
the green curve in Figure 3.

Figure 3: A comparison of the forecasts given by the ’Simple’ (red and blue
curves) and ’Improved’ (green and black curves) models. With the ’Improved’
model (Eqn. (2) with the fitting parameters $a=85710$ and $b=0.23756$) we find
that the growth in the number of cases increases according to a power law with
exponent less than 1. However, should all lock-down measures be lifted on day
100 (i.e. the by the end of April 2020), one can expect the growth to revert
back to its original exponential behaviour (Eqn. (1) with $\alpha=1.1393$).

The problem however is that this type of power-law growth is in fact so
gradual, that it would have to be maintained for several decades before a
significant number of the population become immune to the disease, i.e. as a
result of having contracted the virus. At the same time, should the stringency
of the interventions be relaxed at any time, the spread will revert to its
original exponential growth behaviour, implying that the spread of the disease
could quickly get out of hand, again and again. In Figure 3 we assume, for
example, that all restrictions are lifted at the end of April (day 100),
giving rise to the exponential increase shown by the red curve. The blue curve
shows the daily number of new cases undergoes hardly any reduction in its
maximum peak height, though its peak is shifted forward by about a month, in
comparison to the forecast given in Figure 1.

## Discussion

One can of course try to answer more specific questions with more
sophisticated models, like the discrete model we mentioned for SARS [6];
however, here we have been more interested in developing a very simple model
that brushes over the details and only captures the essential, large scale
behaviour.

As we have seen, there are essentially two possible behaviours: either the
virus can proceed naturally, leading to exponential growth in the number of
cases; or else its progression can be slowed down to a more manageable power-
law growth, i.e. for as long as the present, stringent, lock-down measures can
be imposed. But for how long will it be necessary (or possible) to impose the
present measures? On the near horizon there does not seem to be any easy
answer to this question, nor is there any clear solution to the pandemic
itself. A vaccine is unlikely to be developed for several months, and there
seems to be only a remote possibility that the virus may mutate naturally into
something less lethal to humans.

The direct human cost of an unchecked (exponential) spread of the virus could
be truly devastating[16]. On the one hand, it could result in a catastrophic
loss of tens of millions of lives, as our model predicts, but on the other
hand it will end the pandemic once and for all. The current measures, leading
to more manageable power-law growth, come at a very high economic cost and are
therefore not sustainable for much longer. So far these measures have included
enforced quarantine, which has led to a severe slowdown in economic activity
and manufacturing production, principally due to declining consumption and
disrupted global supply chains[17]. (As an example of the severity of the
slowdown in production, several major car manufacturers are gradually halting
production in major manufacturing hubs throughout the developed world [18].)
This decline, coupled with the associated economic uncertainty, has had knock
on effects in the form of historically unprecedented stock market falls[19].
Although the stock market is more of an indicator of the future value of the
profits of listed corporations, their collapsed share prices could trigger
severe financial crises because of a spike in bankruptcies. (The debt of US
corporations is the highest it has ever been[20].) The inevitable loss of jobs
will also lead to an inability to pay bills and mortgages, increased levels of
crime, etc. In principle, such a major decline in economic conditions could
also result in an equally large-scale loss of life/livelihoods, albeit over a
more prolonged period of time.

## Methods

### The global WHO data

We have fitted our models to the global data extracted from the daily
‘Situation reports’ made available through the world health organisation’s
webpage [21]. For convenience we reproduce this data in Table 2.

Day | Cases | New | Deaths
---|---|---|---
30 | 75204 | 1872 | 2009
31 | 75748 | 548 | 2129
32 | 76769 | 1021 | 2247
33 | 77794 | 599 | 2359
34 | 78811 | 1017 | 2462
35 | 79331 | 715 | 2618
36 | 80239 | 908 | 2700
37 | 81109 | 871 | 2762
38 | 82294 | 1185 | 2804
39 | 83652 | 1358 | 2858
40 | 85403 | 1753 | 2924
41 | 87137 | 1739 | 2977
42 | 88948 | 1806 | 3043
43 | 90869 | 1922 | 3112
44 | 93091 | 2223 | 3198
45 | 95324 | 2232 | 3281
46 | 98192 | 2873 | 3380
47 | 101927 | 3735 | 3486
48 | 105586 | 3656 | 3584
49 | 109577 | 3993 | 3809

Day | Cases | New | Deaths
---|---|---|---
50 | 113702 | 4125 | 4012
51 | 118319 | 4620 | 4292
52 | 125260 | 6741 | 4613
53 | 132758 | 7499 | 4955
54 | 142539 | 9764 | 5392
55 | 153517 | 10982 | 5735
56 | 167511 | 13903 | 6606
57 | 179112 | 11526 | 7426
58 | 191127 | 15123 | 7807
59 | 209839 | 16556 | 8778
60 | 234073 | 24247 | 9840
61 | 266073 | 32000 | 11184
62 | 292142 | 26069 | 12784
63 | 332930 | 40788 | 14510
64 | 372757 | 39827 | 16231
65 | 414179 | 40712 | 18440
66 | 462684 | 49219 | 20834
67 | 509164 | 46484 | 23335
68 | 571678 | 62514 | 26494
69 | 634835 | 63159 | 29957

Day | Cases | New | Deaths
---|---|---|---
70 | 693224 | 58411 | 33106
71 | 750890 | 57610 | 36405
72 | 823626 | 72736 | 40598
73 | 896450 | 72839 | 45525
74 | 972303 | 75853 | 50321
75 | 1051635 | 79332 | 56985
76 | 1133758 | 82061 | 62784
77 | 1133758 | 82061 | 62784
78 | 1279722 | 68766 | 72614
79 | 1353361 | 73639 | 79235
80 | 1436198 | 82837 | 85522
81 | 1521252 | 85054 | 92798
82 | 1610909 | 89657 | 99690
83 | 1696588 | 85679 | 105952
84 | 1773084 | 76498 | 111652
85 | 1844863 | 71779 | 117021

Table 2:  Data used for the Levenberg–Marquardt (least squares) optimisation
of the parameter $\alpha$ in Eqs. 1. The columns contain the day since the
first situation report (21 January 2020), the total number of cases (column
2), the number of new cases (column 3), and the total number of deaths (column
4).

### Python script for fitting the parameters

For the reader’s convenience, the complete python script for the optimisation
of the single parameter $\alpha$ of the simple model is provided on the
following page. In this script, the function `leastsq()`, imported from the
module `scipy.optimize`[22], uses Levenberg–Marquardt optimization to minimize
the residual vector returned by the function `ef()`. The function `leastsq()`
is called from within `main()`, which reads in the data and sets up the
initial parameter and the other two quantities (the initial values `x[0]` and
`y[0]`) for optimisation. These three quantities are then passed to
`leastsq()`, via the vector `v0`.

For the data in Table 2, the output from the script should be:

    R-squared = 0.99432
    alpha = 1.1393
    x[29] = 75194.3
    y[29] = 349.477
    z[29] = 7700000000.0

The python scripts that were used to produce the Figures and Table 1 in this
paper are available from the corresponding author upon request.

    from scipy import linspace,sqrt,dot,concatenate,zeros,loadtxt
    from scipy.optimize import leastsq
    #
    def ef(v,x,y,z,c,xdata,ydata):
        ’’’
        # Residual (error function)
        ’’’
        x[0] = v[1]; y[0] = v[2]
        for i in range(len(xdata)):
            y[i+1] = v[0]*y[i]*(z[i]-x[i])/z[0]
            x[i+1] = x[i] + y[i+1]
            z[i+1] = z[i] - c*y[i]
        er = concatenate(((x[0:-1]-xdata)/xdata,(y[0:-1]-ydata)/ydata))
        return er
    #
    def main():
        f = loadtxt(’data.dat’)     # (WHO data as in Table 1)
        xdata = f[0:41,1]
        ydata = f[0:41,2]
        c = f[-1,3]/f[-1,1]         # fraction of deaths occuring
        n = len(xdata)+1            # number of data points
        m = int(f[0,0])             # day of first data point
        #
        x = zeros(n,’d’)            # for storing total cases
        y = zeros(n,’d’)            # for storing new cases
        z = zeros(n,’d’)            # for storing total population
        #
        x[0] = 76000                # initial guess -- to be optimized
        y[0] = 300                  # initial guess -- to be optimized
        z[0] = 7.7e9                # estimated total global population
        #
        v0 = zeros(3,’d’)           # quantities to be be optimised
        v0[0] = 1.14                # initial guess for alpha
        v0[1] = x[0]
        v0[2] = y[0]
        #
        v,cov,infodict,mesg,ier = leastsq(ef,v0,args=(x[0:n],y[0:n],z[0:n],c,xdata
                                   ,ydata),full_output=True,ftol=1e-13,maxfev=500)
        #
        ys = concatenate((x[0:n-1], y[0:n-1]))
        er = concatenate(((x[0:-1]-xdata),(y[0:-1]-ydata)))
        ssErr = dot(er,er)
        ssTot = ((ys-ys.mean())**2).sum()
        rsquared = 1.0-(ssErr/ssTot )
        x[m] = v[1]; y[m] = v[2]; z[m] = 7.7e9
        print(’R-squared = ’ + str(round(rsquared,5)))
        print(’alpha = ’ + str(round(v[0],5)))
        print(’x[’+str(m)+’] = ’ + str(round(x[m],1)))
        print(’y[’+str(m)+’] = ’ + str(round(y[m],3)))
        print(’z[’+str(m)+’] = ’ + str(round(z[m],1)))

    if __name__ == "__main__":
       main()

### Improved model equations

The form of the equations in our improved model was motivated by a Taylor
series expansion of the power law

$y\left(t\right)=at^{b}\text{, }$ (3)

where $a$ and $b$ are real positive constants. For convenience we first take
the natural logarithm on both sides of (3),

$\ln y\left(t\right)=\ln a+b\ln\left(t\right)\text{.}$ (4)

For $\Delta t\ll t$, we then expand about $t$, to get

$\displaystyle\ln y\left(t+\Delta t\right)$ $\displaystyle=$ $\displaystyle\ln
a+b\ln\left(t+\Delta t\right)$ (5) $\displaystyle\approx$ $\displaystyle\ln
a+b\left(\ln t+\frac{\Delta t}{t}-\frac{1}{2}\left(\frac{\Delta
t}{t}\right)^{2}+\frac{1}{3}\left(\frac{\Delta
t}{t}\right)^{3}-\ldots\right)\text{.}$

Subtracting (4) from (5) produces, to leading order in $\Delta t/t$,

$\ln\left(\frac{y\left(t+\Delta t\right)}{y\left(t\right)}\right)\approx
b\frac{\Delta t}{t}\text{.}$ (6)

From (4) we also see that

$\ln\left(t\right)=\frac{\ln y\left(t\right)-\ln
a}{b}=\ln\left(\frac{y\left(t\right)}{a}\right)^{\frac{1}{b}}\text{,}$ (7)

or

$\frac{1}{t}=\left(\frac{a}{y\left(t\right)}\right)^{\frac{1}{b}}\text{.}$ (8)

After substituting (8) into (6), and taking the exponential on both sides of
the resulting equation, we arrive at

$y\left(t+\Delta t\right)\approx y\left(t\right)\exp\left(b\Delta
t\left(\frac{a}{y\left(t\right)}\right)^{\frac{1}{b}}\right)$ (9)

In terms of an iterative map the power law growth given by (9) suggests the
form

$y_{i+1}=y_{i}\exp\left(b\left[\frac{a}{y_{i}+a}\right]^{\frac{1}{b}}\right)$
(10)

where we have put $\Delta t=1$ and added a shift in the denominator
($y_{i}\rightarrow y_{i}+a$) for convenience. This shift allows the iterative
map model to change smoothly from pure exponential growth ($y_{i}<<a$), to
something that approximates power law growth for values of $y_{i}$ that become
comparable to $a$. Note that, because of the latter approximation, the (now)
fitting parameter $b$ is no longer equal to the power indicated by $b$ in Eqn.
(3).

## References

  * [1] WHO Director-General’s opening remarks at the media briefing on COVID-19 - 11 march 2020.  https://www.who.int/dg/speeches/detail/who-director-general-s-opening-remarks-at-the-media-briefing-on-covid-19—11-march-2020 (Last accessed 03/2020).
  * [2] Covert coronavirus infections could be seeding new outbreaks.  Nature Magazine (20 March 2020) https://www.scientificamerican.com/article/covert-coronavirus-infections-could-be-seeding-new-outbreaks/.
  * [3] Shang, J., Ye, G., Shi, K. _et al._ Structural basis of receptor recognition by SARS-CoV-2.  _Nature_ (2020).
  * [4] Rabajante, J. F.  Insights from early mathematical models of 2019-ncov acute respiratory disease (COVID-19) dynamics.  https://arxiv.org/abs/2002.05296 (Last accessed 03/2020).
  * [5] GLEAMviz: The global epidemic and mobility model.  http://www.gleamviz.org/ (Last accessed 03/2020).
  * [6] Zhou, Y. & Ma, Z.  A discrete epidemic model for SARS transmission and control in China.  _Math. Comp. Modelling_ 40, 1491 (2009).
  * [7] Batista, M.  Estimation of the final size of the COVID-19 epidemic.  _MedRXiv pre-print_ DOI: https://doi.org/10.1101/2020.02.16.20023606 (2020).
  * [8] Marquardt, D. W.  An algorithm for least–squares estimation of nonlinear parameters.  _SIAM J. Appl. Math._ 11, 431 (1963).
  * [9] Worldometer.  World population (2020 and historical).  https://www.worldometers.info/world-population (Last accessed 04/2020).
  * [10] Ziff, A. L. & Ziff, R. M.  Fractal kinetics of covid-19 pandemic, DOI: 10.1101/2020.02.16.20023820 (2020).  %****␣covid19.bbl␣Line␣75␣****https://www.medrxiv.org/content/early/2020/03/03/2020.02.16.20023820.full.pdf.
  * [11] Li, M., Chen, J. & Deng, Y.  Scaling features in the spreading of covid-19.  https://arxiv.org/abs/2002.09199 (2020).  2002.09199.
  * [12] Singer, H. M.  Short-term predictions of country-specific covid-19 infection rates based on power law scaling exponents.  https://arxiv.org/abs/2003.11997 (2020).  2003.11997.
  * [13] Singer, H. M.  The covid-19 pandemic: growth patterns, power law scaling, and saturation.  https://arxiv.org/abs/arXiv:2004.03859 (2020).  2004.03859.
  * [14] Watts, D. J. & Strogatz, S. H.  Collective dynamics of ’small-world’ networks.  _Nature_ 393, 440–442, DOI: 10.1038/30918 (1998).
  * [15] Watts, D. J.  _Small worlds: the dynamics of networks between order and randomness_ , vol. 9 (Princeton University Press, 2004).
  * [16] Virus could have killed 40 million without global response.  (Nature News, 27 March 2020) https://www.nature.com/articles/d41586-020-00154-w (Last accessed 03/2020).
  * [17] A Covid-19 Supply Chain Shock Born in China Is Going Global - 20 March 2020.  https://www.bloomberg.com/news/articles/2020-03-20/a-covid-19-supply-chain-shock-born-in-china-is-going-global (Last accessed 03/2020).
  * [18] Coronavirus: Car production halts at Ford, VW and Nissan - 18 march 2020.  https://www.fleetnews.co.uk/news/manufacturer-news/2020/03/18/coronavirus-car-production-halts-at-ford-vw-and-nissan (Last accessed 03/2020).
  * [19] Coronavirus: FTSE 100, Dow, S&P 500 in worst day since 1987 - 12 March 2020.  https://www.bbc.com/news/business-51829852 (Last accessed 03/2020).
  * [20] A Modern Jubilee as a cure to the financial ills of the Coronavirus - 3 March 2020.  https://www.patreon.com/posts/modern-jubilee-34537282 (Last accessed 03/2020).
  * [21] Coronavirus disease (covid-2019) situation reports.  https://www.who.int/emergencies/diseases/novel-coronavirus-2019 (Last accessed 03/2020).
  * [22] Langtangen, H. P.  _Python scripting for computational science_ (Springer, Berlin, 2008), 3rd edn.

## Acknowledgements

A. E. B. would like to acknowledge M. R. Kolahchi, H. M. Singer and V. Hajnová
for helpful discussions about this work. Both authors wish to thank A. Thomas
for uncovering some of the related literature.

## Author contributions statement

A. E. B devised the research project, developed the models, performed the
numerical simulations and produced all the figures. Both authors analysed the
results and wrote the paper.

## Ethics declarations

The authors declare no competing interests.