# Extending preferred axion models via heavy-quark induced early matter
domination

Andrew Cheek  Jacek K. Osiński  Leszek Roszkowski

###### Abstract

We examine the cosmological consequences of the heavy quarks in KSVZ-type
axion models. We find that their presence often causes an early matter
domination phase, altering the evolution of the Universe. This extends the
axion mass into the region where standard cosmology leads to overproduction,
and allows for a greater number of axion models with non-renormalizable terms
to be viable. Quantitatively, we find that decays proceeding through effective
terms of up to dimension 9 ($d=9$) remain consistent with cosmological
constraints, in contrast with the result $d\leq 5$ previously found in the
literature. As a consequence, the heavy quarks can be much heavier and the
axion mass window with the correct relic density for dark matter is extended
by orders of magnitude, down to $m_{a}\approx 6\times 10^{-9}\,{\rm eV}$. This
is achieved without resorting to fine-tuning of the initial misalignment
angle, bolstering the motivation for many future axion haloscope experiments.
Additionally, we explore how these models can be probed through measurements
of the number of relativistic degrees of freedom at recombination.

##  1 Introduction

The dark matter problem in modern physics is one with a great number of
candidate solutions. Among these candidates, the QCD axion has emerged as a
compelling contender, offering a remarkable explanation to two of the biggest
challenges to fundamental physics: the dark matter, and the strong CP
problems. The solution to the latter can be achieved by the Peccei-Quinn (PQ)
mechanism [1, 2], where a spontaneously broken, global Abelian $U(1)_{\rm PQ}$
symmetry is added to the Standard Model. After QCD confinement, the Nambu-
Goldstone boson of the broken $U(1)_{\rm PQ}$, the axion, obtains a mass

$m_{a}=5.7\times 10^{-5}\,{\rm eV}\left(\frac{10^{11}\,{\rm
GeV}}{f_{a}}\right).$ (1.1)

The decay constant $f_{a}$ is related to the vacuum expectation value of the
scalar field associated with PQ symmetry, $\langle\phi\rangle=f_{a}/\sqrt{2}$.
Flavour-violating processes such as rare kaon decays require that $f_{a}$ is
much larger than the electroweak scale, thus making viable QCD axions very
light and hard to detect. For this reason, invisible axion models are
appealing as dark matter candidates. The attractiveness of axions as a dual
solution to the strong CP and the dark matter problems is intensified when one
realises that the cosmological evolution of the axion field leads to a cold
matter relic abundance in the early Universe [3, 4, 5].

Under Occam’s razor we are compelled to consider this scenario as preferred
until further evidence disfavors it. Additionally, the solution is highly
predictive, motivating a huge experimental effort to test the interaction of
axions with Standard Model (SM) particles, most notably the axion-photon
coupling, $g_{a\gamma}$. Astrophysical objects provide a powerful tool for
constraining the axion mass window [6, 7], but for the mass region most
relevant for axion dark matter, around $10\,{\rm\mu eV}$, laboratory-based
experiments are most promising [8]. In fact, haloscope experiments [9], such
as ADMX [10, 11] and HAYSTAC [12] were able to reach the QCD axion band for a
portion of the masses most relevant for dark matter. This is only the
beginning and there are several operating and proposed haloscopes that are
working to explore a much broader mass region, see Refs. [13, 14, 15, 16] for
reviews.

Within this picture however, it is important to take some deeper theoretical
considerations into account, such as the UV origin of the colour anomaly,
leading to additional particle content. Inspecting the effects of the
additional particles led to the proposal of the preferred axion models [17,
18], where some models were elevated by virtue of their consistency with the
current standard cosmology and behaviour up to the Planck scale, $M_{\rm Pl}$.

The cosmological arguments made for the preferred axion models largely rest on
models in which one introduces new heavy quarks, $Q$, charged under $U(1)_{\rm
PQ}$. In the post-inflationary scenario $Q$s will be produced through thermal
freeze-out and to align with observations, should be unstable. This leads to
the requirement that Lagrangian terms enabling $Q$ decay must proceed at a
renormalizable level ($d\leq 4$) or at the non-renormalizable level at most
through dimension-5 ($d\leq 5$) operators [17, 18].

In this paper we challenge this conclusion, showing that without any
additional particle content, models where $Q$ decays to the SM at the $d>5$
are indeed allowed by existing cosmological constraints. This is because the
new heavy quarks, although sufficiently unstable, in many cases will live long
enough to dominate the cosmic energy budget altering the evolution of the
Universe and, as a result, the dynamics of the misalignment mechanism [19, 20,
21, 22, 23, 24]. This has the consequence of increasing the viable masses for
the heavy quarks, as well as reopening the possibility for axion models that
were previously disregarded. Importantly for axion dark matter, this provides
a self consistent way to extend the axion mass to much lower values without
invoking additional fields or fine-tuning in a pre-inflationary symmetry
breaking scenario.

Additionally, we investigate the impact of the decay of these heavy quarks
into axions and how measurements of the number of relativistic degrees of
freedom, $N_{\rm eff}$, can help constrain their branching ratios. This
provides another avenue in which we alter the view on preferred axion models,
since some lead to large branching ratios into axions and others not.
Interestingly, we find that constraints are of a similar order, independent of
whether the $Q$ particles dominate early on.

This article is organised as follows, in section 2 we briefly review the
relevant aspects of the preferred axion model framework. In section 3 we argue
with approximations why and where we should expect early matter domination
(EMD), then we support our claims with a fully fledged numerical treatment,
showing that $Q$ decays of $d\geq 6$ are completely consistent with
observation. In section 4 we explore the consequences of heavy quark decay for
the additional amount of relativistic degrees of freedom, $\Delta N_{\rm
eff}$. We do this for all decays, including the standard $d\leq 5$ cases.
Finally, in section 5 we conclude the article, emphasising future directions
of study and implications on the search for axion dark matter.

##  2 Preferred axion models

The associated Noether current $J_{\mu}^{\rm PQ}$ for the Abelian $U(1)_{\rm
PQ}$ is conserved up to anomalies

$\partial^{\mu}J_{\mu}^{\rm
PQ}=\frac{g_{s}^{2}N}{16\pi^{2}}G\tilde{G}+\frac{e^{2}E}{16\pi^{2}}F\tilde{F},$
(2.1)

where the first term on the right is the colour anomaly, required for solving
the strong CP problem, the second term is the electromagnetic anomaly, $G$ and
$F$ are the field strength tensors of the gluon and electromagnetic fields ,
respectively, and $\tilde{G}$ and $\tilde{F}$ are their duals. The
coefficients $g_{s}$ and $e$ are the Standard Model (SM) coupling constants,
whereas $N$ and $E$ are the model-dependent anomaly coefficients.

There are two classes of phenomenological axion models that generate a colour
anomaly [25, 8]. The Kim–Shifman–Vainshtein–Zakharov (KSVZ) models, where SM
fields are uncharged under $U(1)_{\rm PQ}$, introduce the colour anomaly via a
heavy fermion, $Q$, which has colour charge [26, 27]. Alternatively, the
Dine–Fischler–Srednicki–Zhitnisky (DFSZ) models imbue SM quarks with
$U(1)_{\rm PQ}$ charges and introduce an additional weak doublet [28, 29]. The
preferred axion model framework [17, 18] is strongly motivated by the
cosmological implications of the heavy quarks and post-inflationary $U(1)_{\rm
PQ}$ breaking such that the reheating temperature is greater than the breaking
scale, $T_{R}>f_{a}$. Therefore, in this article we focus on the KSVZ setup,
which supplements the SM Lagrangian, $\mathcal{L}_{\rm SM}$ by,

$\mathcal{L}_{a}=\mathcal{L}_{\rm SM}+\mathcal{L}_{\rm
PQ}-V_{H\Phi}+\mathcal{L}_{\rm Qq}\,,$ (2.2)

where

$\mathcal{L}_{\rm
PQ}=|\partial_{\mu}\Phi|^{2}+\overline{Q}i\not{D}Q-(y_{Q}\overline{Q}_{L}Q_{R}\Phi+{\rm
H.c})\,.$ (2.3)

The potential $V_{H\Phi}$ contains the new scalar couplings and
$\mathcal{L}_{\rm Qq}$ contains terms that couple $Q$ to SM quarks,
$q=q_{L},d_{R},u_{R}$. We see in eq. (2.3) that the Yukawa term will give rise
to a mass term after PQ breaking, $m_{Q}=y_{Q}f_{a}$. For the purposes of this
work, we take $y_{Q}=1$. Since by construction, $Q_{L}$ and $Q_{R}$ are non-
trivially charged under $SU(3)_{c}$, they will strongly interact, annihilating
into quarks and hadrons. In the early Universe, at temperatures above the QCD
phase transition, $T_{\rm QCD}\sim 150\,{\rm MeV}$, annihilation can be
perturbatively calculated for free $Q$s and SM quarks,

$\langle\sigma
v\rangle_{\bar{Q}Q}=\frac{\pi\alpha_{s}^{2}}{16m_{Q}^{2}}\left(c_{f}n_{f}+c_{g}\right),$
(2.4)

where $n_{f}$ is the number of quark flavours that $Q$ can annihilate into,
and $\left(c_{f},c_{g}\right)=\left(2/9,\,220/27\right)$ for triplets and
$\left(c_{f},c_{g}\right)=\left(2/3,\,27/4\right)$ for octets. With
$T_{R}>f_{a}$, this interaction is strong enough to bring even very heavy $Q$s
into thermal equilibrium with the SM. As with any massive particle, when the
SM bath cools to $T\sim m_{Q}/25$, Boltzmann suppression starts to dominate
and the particle freezes out, leaving a relic abundance. In the case where $Q$
is stable or quasi-stable substantial cosmological issues are present, such as
overabundance or disruption of Big Bang Nucleosynthesis (BBN) [30, 31, 32,
33].

Neither $\mathcal{L}_{\rm PQ}$ nor the scalar potential $V_{H\Phi}$ contain
terms that allow for $Q$ decay, so a stable heavy quark is worth examining, as
indeed has been done in the literature [34, 35]. However, the strong
interaction rate of eq. 2.4 often overproduces $Q$ in the Universe such that
$\Omega_{Q}>\Omega_{\rm DM}$, requiring that $m_{Q}\lesssim 10\,{\rm TeV}$
[17, 18]. This situation is potentially alleviated by noting that at
temperatures below $T_{\rm QCD}$ hadrons start to form and non-perturbative
effects dominate the behaviour of strongly-charged objects. If the heavy
quarks are still prevalent in the Universe, exotic hadrons containing $Q$s may
form, which will likely enhance the annihilation cross-section thanks to
finite size effects. This could potentially restart $Q$-hadron annihilations,
altering the heavy number density, see Refs. [36, 37, 38, 39, 40, 41, 42, 35].
Despite this, the existence of heavy exotic hadrons is still experimentally
disfavored because of a number of phenomenological consequences such as; the
production of anomalously heavy isotopes [43], stellar and neutron star
stability [44, 45], as well as producing anomalously large heat flow in the
Earth’s interior [46].

Therefore, the heavy quarks should be unstable [17, 18]. The strongest
constraints on unstable particles, which at some point where in thermal
equilibrium, typically come from BBN. The limit of which is often taken in
terms of the lifetime of the unstable particle, $\tau\lesssim 0.01\,{\rm s}$
[33]. This limit relies on the assumption of radiation domination throughout
the early Universe [47]. A more robust limit is to ensure that the SM plasma
temperature at time of $Q$ decay is warmer than when BBN starts, i.e. $T_{\rm
decay}^{Q}\gtrsim T_{\rm BBN}\approx 3\,{\rm MeV}$ [48, 49, 50, 51, 52, 53].
With this in mind, we look to terms with heavy $Q$s and SM $q$’s at both the
renormalizable and non-renormalizable levels,

$\mathcal{L}_{Qq}=\mathcal{L}_{Qq}^{d\leq
4}+\mathcal{L}_{Qq}^{d>4}=\mathcal{L}_{Qq}^{d\leq
4}+\frac{1}{\Lambda^{(d-4)}}\mathcal{O}^{d>4}+{\rm h.c.}\,.$ (2.5)

Note that there are only a few charge configurations that lead to
$\mathcal{L}_{Qq}^{d\leq 4}\neq 0$ and the fact that both the $U(1)_{\rm PQ}$
and $U(1)_{Q}$ are expected to be broken at a scale up to the Planck scale,
non-renormalizable terms are also considered. In Refs. [17, 18] the
suppression scale was taken ‘conservatively’ to be $\Lambda=M_{\rm
Pl}=1.22\times 10^{19}\,{\rm GeV}$, but we stress that this is unknown.

The next restriction on axion models comes from the energy density of the
axion itself, which is produced via the misalignment mechanism, resulting in
an abundance [54, 55, 23]

$\displaystyle\Omega_{a}h^{2}\approx\begin{dcases}0.006\left(\frac{\theta_{\rm
i}}{1}\right)^{2}\left(\frac{5.6~{}\mu{\rm eV}}{m_{a}}\right)^{3/2}&\text{ for
}m_{a}\lesssim 3\,H(T_{\rm QCD})\,,\\\ 0.17\left(\frac{\theta_{\rm
i}}{1}\right)^{2}\left(\frac{5.6~{}\mu{\rm eV}}{m_{a}}\right)^{7/6}&\text{ for
}m_{a}\gtrsim 3\,H(T_{\rm QCD}),\end{dcases}$ (2.6)

with the standard cosmological evolution. Here $m_{a}$ is the zero-temperature
axion mass, $\theta_{i}$ is the initial angle of the axion potential and
$T_{\rm QCD}$ is the temperature of the QCD phase transition.

In the post-inflationary PQ symmetry breaking picture, different patches of
the Universe randomly take $\theta_{i}$ values on the unit circle, which is
then averaged to obtain $\langle\theta_{i}\rangle=\pi/\sqrt{3}$. By taking the
thermal relic value from PLANCK, $\Omega_{\rm DM}=0.12$ [56], we can see that
overabundance occurs for $m_{a}\lesssim 145\,{\rm\mu eV}$. Since $m_{a}$ is
determined by the QCD anomaly (1.1), we can read off an approximate limit for
the breaking scale, $f_{a}\lesssim 3\times 10^{11}\,{\rm GeV}$. By taking the
maximum heavy quark mass $m_{Q}\sim f_{a}$, Ref. [18] found that BBN
constraints and $\Omega_{a}<\Omega_{\rm DM}$ requires that $Q$ decays must
proceed through operators of dimension $d\leq 5$, thus providing a powerful
limit on KSVZ type axion models, and a central criterion for preferred axion
models.

This restriction on $d$ has the theoretically awkward consequence that the
$U(1)_{Q}$ symmetry present in $\mathcal{L}_{\rm PQ}$ must be of a low-
quality, i.e. terms that break $U(1)_{Q}$ must occur at low dimension even
with Planck-scale suppression. This is in contrast with the $U(1)_{\rm PQ}$
symmetry, which is required to be of high-quality, where terms that break PQ-
symmetry can only appear at $d\geq 11$ in order not to spoil the solution to
the strong CP problem [57, 58, 59]. Although this situation is not untenable
and can be explained by additional discrete symmetries [60, 17], it would be
theoretically appealing to increase the quality of $U(1)_{Q}$.

In this work we present our findings allowing the initial angle to vary.
$\theta_{i}\in\left[0.5,\pi/\sqrt{3}\right]$, we do this to be more agnostic
on the value of $\langle\theta_{i}\rangle$ and capture some of the effects
that could be relevant for the pre-inflationary breaking scenario without
resorting to a finely tuned angle of $\theta_{i}$, which is ultimately limited
by isocurvature fluctuations [61].

Another consequence of the post-inflationary PQ breaking is the expected
formation of topological defects such as cosmic strings and domain walls [62,
63, 64], both of which could contribute to the density of relic axions,
however, the exact contribution is still largely in dispute due to the
technical challenges when it comes to calculating the effect, see Refs. [65,
66, 67, 68, 69, 70]. Therefore we do not consider these contributions and
stick simply to the misalignment mechanism. Furthermore, domain walls can be
cosmologically dangerous due to their scaling behaviour [71]. This issue is
avoided by only considering models where the domain walls do not form,
specifically this is when the domain wall number is one, $N_{\rm DW}=1$ [72,
73]. However, since there are multiple ways around this problem [71, 74, 75],
it isn’t considered a strong requirement. Despite this, it is worth noting
that including this requirement with the finding of Refs. [17, 18], i.e. that
$Q$ decays must proceed through terms of dimension $d\leq 5$ one is left with
two hadronic axion models [76]

${\rm KSVZ-I}\,:\,\left(3,1,-1/3\right),\,\,\,\,\,{\rm or}\,\,\,\,\,{\rm KSVZ-
II}\,:\,\left(3,1,+2/3\right),$ (2.7)

where the brackets represent the SM charges, $\left({\rm SU}(3)_{c},\,{\rm
SU}(2),\,{\rm U}(1)_{Y}\right)$. From here one is free to choose the PQ
charges $\left(\chi_{L},\chi_{R}\right)$ of the left-handed and right-handed
heavy fermions, but this still is rather restrictive.

In the next section we will discuss how the unstable heavy fermions can affect
the cosmological evolution, the misalignment mechanism, and thus the window of
preferred axion models.

##  3 Early heavy Q domination and axion relic abundance

As described above, the axion models of interest here are those with very
heavy $Q$ fermions that have a charge under ${\rm SU}(3)_{c}$. The fast
interaction rate ensures that above the temperature $\sim m_{Q}$, the $Q$s
will be in thermal equilibrium. Once the temperature of the plasma cools below
$m_{Q}$ the heavy particle freezes out in the standard way described by the
Boltzmann equation for the number density of $Q$, $n_{Q}$,

$\frac{\differential n_{Q}}{\differential
t}+3Hn_{Q}=-\langle\sigma\,v\rangle\left[n_{Q}^{2}-(n_{Q}^{\rm
eq})^{2}\right],$ (3.1)

where for now we ignore the decay terms of $Q$. The superscript “${\rm eq}$”
denotes the equilibrium distribution, which gets Boltzmann suppressed at
temperatures below the particle mass, $T<m$. To get an idea about the
temperatures the thermal bath may have when a meta-stable $Q$ starts to
dominate the energy of the Universe, we can take the approximate textbook
result

$Y_{Q}^{\rm f}=\frac{n_{Q}}{s_{R}^{\rm SM}}\approx\frac{45}{2\pi^{2}g_{\star
S}(T_{\rm f})}\frac{x_{\rm f}^{3}H(T_{\rm f})}{m_{Q}^{3}\langle\sigma
v\rangle},$ (3.2)

where $s_{R}^{\rm SM}$ is the SM plasma entropy density and $T_{\rm
f}=m_{Q}/x_{\rm f}$ is the freeze-out temperature and numerically $x_{\rm
f}\approx 25$. The effective number of entropy degrees of freedom is given by
$g_{\star S}$. Assuming a radiation-dominated Universe prior to and at the
time of freeze-out one can obtain an estimate for the ratio of energy
densities for $Q$ and SM radiation

$\frac{\rho_{Q}}{\rho_{\rm R}^{\rm
SM}}=\frac{m_{Q}\,n_{Q}\left(m_{Q},T\right)}{\rho_{\rm R}^{\rm
SM}\left(T\right)}\sim 10^{10}\left(\frac{m_{Q}}{10^{12}\,{\rm
GeV}}\right)^{2}\left(\frac{1\,{\rm MeV}}{T}\right),$ (3.3)

where we have scaled $m_{Q}$ and $T$ values to the approximate mass scale for
the misalignment mechanism to produce axion dark matter and the temperature at
BBN ($T\sim 3\,{\rm MeV}$). This implies that the $Q$ particles of mass
$10^{12}\,{\rm GeV}$ will start to dominate the Universe at $T_{\rm eq.}\sim
10^{7}\,{\rm GeV}$, well above the temperature at the onset of BBN, therefore,
if the $Q$ decays are slow enough, many axion models in fact feature a period
of EMD.

To calculate the total decay width, $\Gamma_{Q}$, we approximate the behaviour
of $\mathcal{L}_{Qq}$, following Ref. [18] by assuming a constant matrix
element and massless final states to obtain an analytic phase space integral
[77],

$\Gamma_{d,n_{f}}=\frac{m_{Q}}{4\left(4\pi\right)^{2n_{f}-3}\left(n_{f}-1\right)!\left(n_{f}-2\right)!}\left(\frac{m_{Q}^{2}}{\Lambda^{2}}\right)^{d-4},$
(3.4)

where $n_{f}$ is the number of final state particles, for $d=5,6,7,8$ we have
at least $n_{f}=2,3,4,5$. This leads to the corresponding decay widths,

$\displaystyle\Gamma_{d=4}$ $\displaystyle=\frac{m_{Q}}{8\pi},\,\,\,\,$
$\displaystyle\Gamma_{d=5}$
$\displaystyle=\frac{m_{Q}^{3}}{16\pi\Lambda^{2}},$
$\displaystyle\Gamma_{d=6}$
$\displaystyle=\frac{m_{Q}^{5}}{512\pi^{3}\Lambda^{4}},\,\,\,\,$
$\displaystyle\Gamma_{d=7}$
$\displaystyle=\frac{m_{Q}^{7}}{49152\pi^{5}\Lambda^{6}},$
$\displaystyle\Gamma_{d=8}$
$\displaystyle=\frac{m_{Q}^{9}}{576\left(4\pi\right)^{7}\Lambda^{8}},\,\,\,\,$
$\displaystyle\Gamma_{d=9}$
$\displaystyle=\frac{m_{Q}^{11}}{11520\left(4\pi\right)^{9}\Lambda^{10}}\,,\,$
(3.5)

where, unlike Ref. [18] we do not set $\Lambda\rightarrow M_{\rm Pl}$,
although this is a reasonable upper limit to take because it is expected that
the PQ symmetry will be broken at or below such a scale.

To estimate the plasma temperature at time of $Q$ decay, $T^{Q}_{\rm decay}$,
we set $\Gamma_{Q}\equiv H_{\rm rad}\left(T=T^{Q}_{\rm decay}\right)$ where
$H_{\rm rad}$ is the Hubble parameter in the radiation-dominated Universe.
Generally this leads to the equation

$T^{Q}_{\rm decay}=\left(\frac{90M_{\rm
Pl}^{2}}{8\pi^{3}}\frac{(\Gamma^{Q}_{d})^{2}}{g_{\star}\left(T^{Q}_{\rm
decay}\right)}\right)^{1/4}.$ (3.6)

where numerically one has to iterate with the $T$ dependence in $g_{\star}$ to
get $T^{Q}_{\rm decay}$. For $d=4$ the decay rate is very fast and typically
$T^{Q}_{\rm decay}>T_{\rm f}$. However for dimension-5 decays, with maximum
suppression $\Lambda=M_{\rm Pl}$, we have

$T^{Q,\,d=5}_{\rm decay}\approx 3.2\times 10^{7}\,{\rm
GeV}\,\left(\frac{m_{Q}}{10^{12}\,{\rm
GeV}}\right)^{3/2}\,\left(\frac{1}{g_{\star}\left(T^{Q,\,d=5}_{\rm
decay}\right)}\right)^{1/4}.$ (3.7)

By inspecting when $T^{Q,\,d=5}_{\rm decay}\sim T_{\rm eq.}$, this equation
suggests that $d=5$ decays could produce EMD for masses $m_{Q}\gtrsim
10^{12}\,{\rm GeV}$. However, by comparing $T_{\rm eq}$ and $T_{\rm
decay}^{Q,\,d=5}$, we find that the period for $Q$-domination is very short.
On the other hand, models with dimension $d\geq 6$ decays have a much greater
time separation between the end of freeze-out, the temperature of quark
domination and decay. This means the period of EMD is much longer and thus has
a greater influence.

Figure 1: Heavy quark lifetimes shown in terms of the temperature of the SM
plasma at time of decay $T_{\rm decay}^{Q}$. Different colored lines show the
lifetimes via $d=5,\,6,\,7$ effective operators. LHC constraints are shaded
red and BBN constraints from $T\lesssim 3\,{\rm MeV}$ are hatched gray. The
orange shaded region shows where our estimated early $Q$ domination should
occur according to eq (3.3). The green hatched region is where the
misalignment mechanism overproduces axions, $\Omega_{a}>\Omega_{\rm DM}$,
assuming standard cosmology, if we take the $\theta_{i}=0.5$ value. The
temperature at which the misalignment mechanism starts to produce particles,
$T_{\rm osc}$, is shown with the black dotted line, also assuming standard
cosmology. Therefore we expect variation to the green hatched region and black
dotted line when they overlay with the orange region. The top axis is the
axion mass given by eq. (1.1) with $f_{a}=m_{Q}$.

We sketch this out in Fig. 1 where we show the $Q$ lifetimes for different
dimensions in terms of temperature at decay. Indeed our estimate for
$\rho_{Q}\approx\rho_{\rm SM}$ intersects with $T^{Q,\,d=5}_{\rm decay}$ at
around $m_{Q}\sim 10^{12}\,{\rm GeV}$ but the difference between $T_{\rm eq}$
and $T_{\rm decay}$ is not that large. We show where our above estimate
predicts early $Q$ domination by the orange shaded region and we see that for
$d=6$ & $7$ this should occur for a substantial period of time. We also show
the region of parameter space where misalignment would overproduce cold
axions, assuming radiation domination and initial angle $\theta_{i}=0.5$
111Ref. [18] chose $\theta_{i}=\pi/\sqrt{3}$ to define a more aggressive
constraint.. Importantly this figure suggests that the assumption of standard
early cosmology is not well founded for $d\geq 5$ for $m_{Q}$ above $\sim
10^{11}\,{\rm GeV}$.

It has been shown that a period of EMD can augment the dynamics of the
misalignment mechanism. The important physical quantity to evaluate when the
axion field starts to oscillate around its true minimum is the oscillation
temperature, $T_{\rm osc}$, defined by $3H\left(T_{\rm
osc}\right)=\tilde{m}_{a}\left(T_{\rm osc}\right)$. The temperature-dependent
axion mass $\tilde{m}_{a}(T)$ is given by

$\tilde{m}_{a}^{2}(T)=\frac{\chi(T)}{f_{a}^{2}}\,,$ (3.8)

where $\chi(T)$ is the QCD susceptibility, which has been calculated using
lattice QCD [78]. In the WKB approximation, $\theta\ll 1$ and slowly varying
$H$ and $\tilde{m}_{a}$, and a radiation-dominated early Universe, the
resulting oscillation temperature is [23]

$\displaystyle T_{\rm
osc}\begin{dcases}\left(\frac{1}{\pi}\sqrt{\frac{80\pi}{g_{\star}\left(T_{\rm
osc}\right)}}m_{a}M_{\rm Pl}\right)^{1/2}&\text{ for }T_{\rm osc}\leq T_{\rm
QCD}\,,\\\ \left(\frac{1}{\pi}\sqrt{\frac{80\pi}{g_{\star}\left(T_{\rm
osc}\right)}}m_{a}M_{\rm Pl}T_{\rm QCD}^{4}\right)^{1/6}&\text{ for
}m_{a}\gtrsim 3\,H(T_{\rm QCD}).\end{dcases}$ (3.9)

This estimate is depicted in Fig. 1 by a black dotted line. We can get an idea
of whether early $Q$ domination is likely to affect the dynamics of
misalignment by seeing if $T_{\rm osc}>T_{\rm decay}^{Q}$ inside the orange
region. One can see that for $d=5$ decay occurs at temperatures far too high
and thus the misalignment mechanism is unlikely to be altered, but for $d\geq
6$ we can expect some effects. Interestingly this will occur in the region
which is prohibited by axion overproduction in standard cosmology. The effects
of EMD have been shown to dilute the axion abundance if the matter decays
after $T_{\rm osc}$ [19, 20, 21]. Therefore, this picture suggests that the
axion models where $Q$ decays through $d\geq 6$ effective operators may in
fact not be inconsistent and should not be disregarded. We evaluate this more
precisely below.

###  3.1 Full numerical treatment

The estimates above only imply that $d\geq 6$ axion models should be
reassessed by the preferred axion models framework. To properly determine
whether this is the case we perform a full numerical evaluation of the coupled
Friedmann-Boltzmann equations along with the axion equation of motion. These
track the evolution of particle species that were produced thermally or
through $Q$ decay, or via the misalignment mechanism. The Friedmann equation
is

$\frac{3H^{2}M_{\rm Pl}^{2}}{8\pi}=\rho_{\mathrm{R}}^{\rm
SM}+\rho_{a}+\rho_{Q}\,,$ (3.10)

where $\rho$ is the energy density, denoted for each component via the
super/sub-scripts. The set of Boltzmann equations that describe the dynamics
of the SM plasma entropy density, $s_{R}^{\rm SM}$, the $Q$ number density,
$n_{Q}$, and the axion energy density, $\rho_{a}$, through particle
creation/annihilation processes are

$\displaystyle\frac{\differential s_{\mathrm{R}}^{\rm SM}}{\differential t}$
$\displaystyle=$ $\displaystyle-3Hs_{\mathrm{R}}^{\rm SM}+\frac{{\rm BR}_{\rm
SM}\Gamma_{Q}}{T}\rho_{Q}\,,$ (3.11a)
$\displaystyle\frac{\differential\rho_{a}}{\differential t}$ $\displaystyle=$
$\displaystyle-4H\rho_{a}+{\rm BR}_{a}\Gamma_{Q}\rho_{Q}+\langle E_{\rm
scat}^{a}\rangle\gamma_{a}\left(1-\frac{n_{a}}{n_{a}^{\rm eq}}\right)\,,$
(3.11b) $\displaystyle\frac{\differential n_{Q}}{\differential t}$
$\displaystyle=$
$\displaystyle-3Hn_{Q}-\Gamma_{Q}n_{Q}-\langle\sigma\,v\rangle\left[n_{Q}^{2}-(n_{Q}^{\rm
eq})^{2}\right]\,,$ (3.11c)

where the thermalized annihilation cross-section, $\langle\sigma v\rangle$,
has been introduced in eq. 2.4. The third term on the right hand side of eq.
(3.11b) determines whether the axions are in thermal equilibrium with the SM
plasma, which proceeds via the thermal axion production rate $\gamma_{a}$
[79]. The averaged axion scattering energy, $\langle E_{\rm scat}^{a}\rangle$,
describes the average energy released in each collision between the thermal
bath and the axion, $\sim 3T$ [80]. We find that when heavy $Q$ domination
occurs this term is irrelevant, but in section 4 we explore $d\leq 5$ cases
where it becomes important, and we discuss it in more detail there. We denote
‘simplified’ branching ratios for $Q$ decay into SM and axions by ${\rm
BR_{SM}}$ and ${\rm BR}_{a}$, respectively. These are simplified because we
have performed our analysis in a model-independent way and focus primarily on
the energy transfer between different species. For example, with the KSVZ-I
model which has $\left(\chi_{L},\chi_{R}\right)=(-1,-2)$, the lowest-order $Q$
decay term is

$\mathcal{L}_{Qq}=\Phi^{\dagger}\overline{Q}_{L}d_{R}+\frac{1}{\Lambda^{(d-4)}}\mathcal{O}^{d>4}$
(3.12)

therefore the conventional particle physics branching ratio for $Q\rightarrow
a+d$ would be $1$, but because the energy will be distributed evenly to both
the SM and axion energy densities, we will have ${\rm BR}_{\rm SM}={\rm
BR}_{a}=1/2$. We have performed the calculation this way to preserve model
independence and in section 4 we will elucidate implications for some specific
models. In this section however we assume ${\rm BR}_{\rm SM}=1.0$ such that
all axions are produced either thermally from bath particles or the
misalignment mechanism.

To model the misalignment mechanism one solves the equation of motion for the
axion field, $\tilde{a}$, which is expressed in terms of the axion angle,
$\theta=\tilde{a}/f_{a}$

$\left(\frac{\differential^{2}}{\differential
t^{2}}+3H(t)\frac{\differential}{\differential
t}\right)\theta(t)+\tilde{m}_{a}^{2}(t)\sin(\theta(t))=0\,,$ (3.13)

where $\tilde{m}_{a}(t)$ is the temperature- (and therefore time-) dependent
mass of the axion defined in eq. (3.8).

In practice, we solve the background cosmology described in equations (3.10)
and (3.11), obtaining the resultant $H$ up until the time of BBN. We then
solve eq. (3.13) using the dedicated misalignment solver MiMeS [81]. Indeed as
our estimates suggested we observe significant regions of parameter space
where heavy quarks dominate the early Universe and alter the cold relic axion
abundance. In Fig. 2 we show the axion abundance for different dimensions in
$\mathcal{L}_{Qq}$, $d=5,\,6,\,7,\,8$ in green, blue, red, and purple,
respectively. The bands extend over the domain of initial misalignment angles
which are not fine-tuned, $\theta_{i}\in\left[0.5,\pi/\sqrt{3}\right]$,
showing how $\Omega_{a}$ varies over $m_{Q}=f_{a}$. Remember that the
abundance $\Omega_{a}\propto\theta_{i}^{2}$ in the standard scenario so
$\theta_{i}=\pi/\sqrt{3}$ gives the lower $m_{Q}$ in the bands shown. For this
figure, we have set the cut-off scale for decays to $M_{\rm Pl}$ and convey
where BBN constraints become relevant (i.e. $T_{\rm decay}^{Q}\lesssim 3\,{\rm
MeV}$) by black hatching.

Figure 2: Misalignment abundance of axions $\Omega_{a}$ where heavy $Q$s can
decay via operators of different dimensions, $d\leq 5\,({\rm
green}),\,d=6\,({\rm blue}),\,d=7\,({\rm red}),\,{\rm and}\,d=8\,({\rm
purple})\,$. The horizontal dark grey band shows the value for $\Omega_{\rm
DM}$ as reported by the PLANCK collaboration [82].

As anticipated by the approximations in Fig. 1, our full calculations
exhibited in Fig. 2 confirm that the misalignment mechanism is not altered by
$d\leq 5$ decays, and what we observe is the result which corresponds to the
standard cosmology of radiation domination. However, $d=6,7,8$ do diverge from
the standard misalignment mechanism, in the process varying the $f_{a}$ (and
$m_{a}$) that corresponds to the correct relic abundance for axion dark
matter. We see that for all cases, when $m_{Q}$ is high enough they will tend
to the standard cosmology limit. This is consistent with Fig 1, where the
heavy quarks have decayed before the axion field starts to oscillate, $T_{\rm
decay}^{Q}>T_{\rm osc}$. Of course $T_{\rm osc}$ is itself also affected by
matter domination, but the rough picture still corresponds well. In both
figures we see that $m_{Q}\gtrsim 3\times 10^{12}\,{\rm GeV}\,{\rm
and}\,\gtrsim 7\times 10^{14}\,{\rm GeV}$ should correspond to the standard
misalignment result for $d=6$ and $7$, respectively.

By taking $\Lambda=M_{\rm Pl}$, and assuming standard cosmology, we see that
$d\geq 6$ would come into tension with BBN constraints especially if you take
the value $\theta_{i}=\pi/\sqrt{3}$. This is what is the primary motivator for
limiting the preferred axion models to $d\leq 5$. However, when including the
effects that frozen-out $Q$s will have on the early Universe the result is
that axions are only overproduced when $m_{Q}\gtrsim 2\times 10^{12}\,{\rm
GeV}$. For $d=7$ overproduction occurs when $m_{Q}\gtrsim 2\times
10^{14}\,{\rm GeV}$.

As mentioned above we include the standard $\theta_{i}$ range to explore how
reasonable variations on $\theta_{i}$ affect the results. We can see that for
axion dark matter, dimension-6 and 7 models have a much narrower mass window
than in the standard case. We even see that we can push the dimension of
effective operators to 8, but BBN limits occur almost at the same point as
$\Omega_{a}\approx\Omega_{\rm DM}$, making it impossible for
$\theta_{i}=\pi/\sqrt{3}$ to obtain the correct relic abundance. However, the
suppression scale $\Lambda$ is ultimately unknown. It is expected to break at
most at the Planck scale, but potentially below it.

In Fig. 3 we treat $\Lambda$ as a free parameter. Now the shaded regions show
what values of $m_{Q}$ and $\Lambda$ achieve the correct relic $\Omega_{\rm
DM}$ with the domain $\theta_{i}\in\left[0.5,\pi/\sqrt{3}\right]$. We see that
by varying $\Lambda$, decays at larger values of $d$ can occur without
altering the misalignment mechanism. This is shown by the large vertical
regions of each case where $m_{a}\sim\left[2\times 10^{-5}-2\times
10^{-4}\right]\,{\rm eV}$ corresponding to the standard-cosmology window. More
interesting is when the coloured shaded regions diverge from this pattern,
when EMD affects the misalignment mechanism. We see that the bands are much
more narrow, and occur for larger values of $\Lambda$ and $m_{Q}$. This is
because a larger $\Lambda$ suppresses the decay rate, making the
$Q$-domination period longer. At the same time, increasing $m_{Q}$ increases
$\Gamma_{Q}$ to a higher power than $\Lambda$. The driving effect for higher
$m_{Q}$ is that their freeze-out is sooner and the corresponding yield is
larger (see eq. (3.3)), further contributing to an increased period of EMD.

Figure 3: The values of $m_{Q}$ and $\Lambda$ that satisfy the correct relic
abundance through the misalignment mechanism assuming an initial misalignment
angle which is not fine-tuned, $\theta_{i}\in\left[0.5,\pi/\sqrt{3}\right]$ .
Colours blue, red, purple, and orange are for $d=6,7,8\,$ and $9$,
respectively.

Once again the left and right boundaries of the shaded regions correspond to
$\theta_{i}=\pi/\sqrt{3}$ and $0.5$, respectively. Above and to the left of
each region corresponds to underproduction of axions with respect to
$\Omega_{\rm DM}$ and therefore remains a viable area of parameter space for
axion models where dark matter is not fully axionic. This is limited from
above by the BBN constraints that are shown for different dimensions. Below
and to the right of the shaded colored regions are areas that would require a
fine-tuned $\theta_{i}<0.5$ to not overproduce the axions, helping to put
limits on the preferred axion window. We see also that by varying $\Lambda$ we
can include $d=9$ and the dark matter abundance can be realised at values of
$m_{Q}\sim 10^{15}\,{\rm GeV}$. As exhibited by the top axis in our figure,
the variation in $m_{Q}$ implies that the axion mass can vary by orders of
magnitude, bringing the benchmark region into the realms of low-mass axion
experiments such as ADMX [83], ABRACADABRA [84], FLASH [85], DMRadio [86], and
SRF [87]. This is the first time that the preferred axion models are shown to
actually predict masses down to $m_{a}\sim 10^{-8}\,{\rm eV}$ without invoking
fine-tuning in $\theta_{i}$ or additional fields that alter the cosmology.

To determine the relation between $\Lambda$ and $m_{Q}$ that these early
$Q$-dominated branches exhibit we refer to the expressions derived in Appendix
B of Ref. [23]. Our results are consistent with both the non-adiabatic and
adiabatic EMD expressions when $T_{\rm osc}>T_{\rm QCD}$: $\Lambda\propto
m_{Q}^{\frac{d-11/4}{d-4}}$ and $\Lambda\propto m_{Q}^{\frac{d-65/22}{d-4}}$,
respectively. This leads us to an important point that was alluded to in the
Introduction, namely the question of exotic hadron formation and how their
formation could affect the physics of the early Universe. Typically the
formation of these objects is considered in a radiation-dominated Universe and
would be worth exploring in the context of $Q$ domination, but this is beyond
the scope of this work. To estimate where exotic hadron formation may affect
the misalignment mechanism we have determined where $T_{\rm osc}=T_{\rm QCD}$,
shown as black dotted lines in Fig. 3. There are two lines shown, one each for
$d=8\,{\rm and}\,9$, because for the others, the line is in the region already
constrained by BBN. For $d=9$ we show a zoomed in region to show that there is
indeed a small region where $T_{\rm osc}<T_{\rm QCD}$, so that misalignment
may occur during this period where exotic hadrons have formed. In this period
of heavy $Q$ domination, the likely hadrons form into $\overline{Q}Q$ heavy
mesons, thus continuing the period of EMD and therefore not altering the
results shown in this section.

##  4 Constraints from $\Delta N_{\rm eff}$

We have shown in the previous section that the heavy quarks required in KSVZ-
type axion models may produce a period of EMD which can alter the dynamics of
the misalignment mechanism. At the end of the early matter stage, the dominant
energy density is required to transfer predominantly to SM radiation in order
to fit with BBN and the formation of the Cosmic Microwave Background (CMB).
For the latter, any additional radiation alters small-scale anisotropies and
is parameterized by the effective number of relativistic degrees of freedom,
$N_{\rm eff}$. This quantity is highly sensitive to the thermal history of
particle physics content. At high temperatures, all SM particles, including
neutrinos, are in thermal equilibrium. Through expansion and cooling, most
particles thermally decouple when temperatures fall below their mass, $T<m$.
Only neutrinos decouple when $T\gg m_{\nu}$, and later entropy injections to
the SM bath separate the temperature of the photon, $T_{\gamma}$, and neutrino
$T_{\nu}$ baths.

As noted above, many axion models predict $Q$ decay into axions and ${\rm
BR}_{a}$ in the $\rho_{a}$ Boltzmann equation (eq. (3.11b)) will no longer be
zero. Given the huge mass difference between $Q$ and the axion, the resultant
axion will be relativistic, unlike the axions produced through misalignment.
Therefore, the energy density from eq. (3.11b) will contribute to the energy
density of radiation today as dark radiation. An additional amount of
relativistic degrees of freedom, $\Delta N_{\rm eff}$, is defined as,

$\rho_{\mathrm{R}}\equiv\rho_{\gamma}\left[1+\frac{7}{8}\left(\frac{T_{\nu}}{T_{\gamma}}\right)(N_{\mathrm{eff}}^{\rm
SM}+\Delta N_{\mathrm{eff}})\right]\,,$ (4.1)

where $N_{\rm eff}^{\rm SM}=3.045$ denotes the effective number of
relativistic neutrinos in the SM [88]. Using the above equation together with
$\rho_{\rm Rad}=\rho_{\rm R}^{\rm SM}+\rho_{a}$ allows us to solve for $\Delta
N_{\rm eff}$

$\Delta N_{\rm
eff}\equiv\left\\{\frac{8}{7}\left(\frac{4}{11}\right)^{-\frac{4}{3}}+N_{\rm
eff}^{\rm SM}\right\\}\frac{\rho_{a}}{\rho_{\rm R}^{\rm SM}}\,.$ (4.2)

Both $\rho_{\rm R}^{\rm SM}$ and $\rho_{a}$ are evaluated at the time of the
CMB, around matter-radiation equality $T_{\rm eq}=0.75$ eV. We solve
numerically the Friedmann-Boltzmann equations shown in Eq. (3.11) only until
BBN ($T_{\rm BBN}$). The thermally decoupled axions will evolve independently
of the SM and the energy densities of the two sectors can be easily
extrapolated using entropy conservation

$\displaystyle\Delta N_{\rm
eff}=\left\\{\frac{8}{7}\left(\frac{4}{11}\right)^{-\frac{4}{3}}+N_{\rm
eff}^{\rm SM}\right\\}\frac{\rho_{a}(T_{\rm BBN})}{\rho_{\rm R}^{\rm
SM}(T_{\rm BBN})}\left(\frac{g_{*}(T_{\rm BBN})}{g_{*}(T_{\rm
eq})}\right)\left(\frac{g_{*S}(T_{\rm eq})}{g_{*S}(T_{\rm
BBN})}\right)^{\frac{4}{3}}\,.$ (4.3)

Additionally eq. (3.11b) contains the scattering terms between axions, the SM
and heavy quarks. If these interactions are frequent enough, they will
thermalise. Below the $m_{Q}$ threshold the thermal production rate is
dominated by the dimension-5 anomalous term coupling the axion to the gluons
[89, 90]. Above $m_{Q}$, the scattering of the heavy quarks dominate even at
temperatures much greater than this threshold as shown in Ref. [79]. We have
used the KSVZ axion production rate $\gamma_{a}$ from gluon scattering given
in Ref. [79] as well as the cross-sections in order to determine $\gamma_{a}$
for any $m_{Q}$ value. We do this primarily for ease and because KSVZ-type
axion models are slightly more appealing because it is possible to achieve
$N_{\rm DW}=1$ with them. When we make the association $m_{Q}\sim f_{a}$ we
can find an approximate relation for when axion decoupling occurs,
$\gamma_{a}\left(T_{\rm dec}^{a}\right)\approx 3n_{a}^{\rm eq}H\left(T_{\rm
dec}^{a}\right)$,

$T_{\rm dec}^{a}\approx 3.2\times 10^{-3}\,(m_{Q})^{1.13}\,\,\,\,\,\,{\rm
for}\,\,\,\,m_{Q}\sim f_{a}\gtrsim 10^{10}\,{\rm GeV},$ (4.4)

where below $m_{Q}\sim 10^{10}\,{\rm GeV}$ the scaling is different, and is
irrelevant for this work. Furthermore, considering the fact that, through EMD,
the value for $f_{a}$ that results in the correct relic abundance is
increased, we should keep in mind that the maximum value of $f_{a}$ that can
support thermalization at some point in the history is approximately given by
[91]

$\displaystyle f_{a}\approx 3\times 10^{10}\,{\rm
GeV}\sqrt{\frac{m_{Q}}{10^{5}\,{\rm GeV}}}.$ (4.5)

With $m_{Q}\sim f_{a}$, axions will never equilibrate with the SM only when
$f_{a}=m_{Q}\gtrsim 10^{16}\,{\rm GeV}$, which is substantially higher than
the value considered in this paper. Therefore we find it important to
incorporate the axion production rate from the thermal plasma in our analysis.

###  4.1 $\Delta N_{\rm eff}$ with standard cosmology ($d\leq 5$ models)

By comparing the SM temperature at time of $Q$ decay $T_{\rm decay}^{Q}$,
defined by $\Gamma_{Q}\approx H(T_{\rm decay}^{Q})$, one can estimate whether
the heavy quarks will decay prior to or after the axions have decoupled.

For $d=4$, the relation is approximately

$\frac{T_{\rm decay}^{Q,\,d=4}}{T_{\rm
dec}^{a}}\approx\left(\frac{1}{8\pi}\right)^{1/2}\left(\frac{1}{g_{\star}}\right)^{1/4}\left(\frac{90M_{\rm
Pl}^{2}}{8\pi^{3}}\right)^{1/4}\left(\frac{1}{3.2\times
10^{-3}\,(m_{Q})^{1.13}}\right)\approx 16\sqrt{\frac{M_{\rm Pl}}{m_{Q}}}$
(4.6)

and since $M_{\rm Pl}>m_{Q}$, we have $T_{\rm decay}^{Q,\,d=4}>T_{\rm
dec}^{a}$, meaning that we expect the standard $\Delta N_{\rm eff}\approx
0.027$, which itself is actually in reach of the CMB-HD projections.
Furthermore, with the approximation $f_{a}\approx m_{Q}$, achieving the
appropriate relic density for axion dark matter ($f_{a}\sim 10^{12}\,{\rm
GeV}$) suggests $T_{\rm decay}^{Q,\,d=4}\sim 10^{4}\,T_{\rm dec}^{a}$.

Since the thermal axion production rate $\gamma_{a}$ is determined by both
$m_{Q}$ and $f_{a}$ independently, it is worth exploring the possibility where
$m_{Q}\ll f_{a}$. This way one can envisage a situation where axions never
couple through interactions mediated by the QCD anomaly term, such that
$T_{\rm dec}^{a}>T_{\rm f}^{Q}$. However, in the dimension-4 case the decay
rate is large enough to keep the axions in thermal equilibrium through decay
and inverse decay of $Q$s as described by the following Boltzmann equation
[92, 93, 94]

$\frac{\differential n_{a}}{\differential
t}+3Hn_{a}=\left(\frac{m^{2}T}{2\pi^{2}}K_{1}\left(\frac{m_{Q}}{T}\right)\Gamma_{Q\rightarrow
a+q}+\gamma_{a}\right)\left(1-\frac{n_{a}}{n_{a}^{\rm eq.}}\right)\,,$ (4.7)

where $K_{1}(m_{Q}/T)$ is the modified Bessel function of the second kind, and
the $Q$ particles are produced in equilibrium at temperatures above and around
freeze-out, $T_{\rm f}^{Q}\sim m_{Q}/25$. This has the effect of bringing the
axions into equilibrium regardless of the Peccei-Quinn breaking scale.
Therefore for $d=4$ preferred axion models the prediction of $\Delta N_{\rm
eff}=0.027$ is fairly robust. Of course in works such as Ref. [95, 79, 91], it
has been shown that in axion models where one does not explain dark matter one
can enhance the number of relativistic degrees of freedom, but typically this
is for values of $f_{a}$ substantially smaller than the values discussed in
this work.

For dimension-5 models we can once again compare $T_{\rm decay}^{Q}$ and
$T_{\rm dec}^{a}$,

$\frac{T_{\rm decay}^{Q,\,d=5}}{T_{\rm
dec}^{a}}\approx\left(\frac{m_{Q}^{3}}{16\pi\Lambda^{2}}\right)^{1/2}\left(\frac{1}{g_{\star}}\right)^{1/4}\left(\frac{90M_{\rm
Pl}^{2}}{8\pi^{3}}\right)^{1/4}\left(\frac{1}{3.2\times
10^{-3}\,(m_{Q})^{1.13}}\right)\approx 10\frac{\sqrt{m_{Q}M_{\rm
Pl}}}{\Lambda}$ (4.8)

which suggests that, when taking the value $f_{a}\sim m_{Q}\sim 10^{12}\,{\rm
GeV}$ that reproduces $\Omega_{\rm DM}$, $Q$ decay occurs after axions
thermally decouple when $\Lambda\gtrsim 10^{16}\,{\rm GeV}$. Instead, if we
set $\Lambda=M_{\rm Pl}$, we are met once again with a $m_{Q}/M_{\rm Pl}$
ratio, only this time it indicates that decays will occur after axion
decoupling as long as $m_{Q}\lesssim M_{\rm Pl}/100$. Therefore, the decays of
heavy $Q$s will directly alter the density of axions with respect to SM
radiation, despite never dominating the Universe.

Figure 4: Top row shows the evolution of energy density for SM radiation
(blue), heavy quarks, $Q$, (orange) and axions that contribute to dark
radiation (green). We show three cases for $Q$ decays that proceed via
dimension-5 non-renormalizable terms where the cutoff is $\Lambda=M_{\rm Pl}$
(solid and dashed lines for ${\rm BR}_{\rm SM}=1/2$ and ${\rm BR}_{\rm SM}=1$,
respectively) or $\Lambda=10^{16}\,{\rm GeV}$ (dotted lines). We show how the
branching ratio of $Q\rightarrow{\rm axions}$ affects the evolution. In the
bottom row we show the ratio $\rho_{a}/\rho_{\rm SM}$ which contributes to
$\Delta N_{\rm eff}$ Eq. (4.3). We take $m_{Q}\sim f_{a}$.

In Fig. 4 we show the evolution of the Universe for three $d=5$ cases,
differentiated by the branching ratio of decay to the SM and by the decay cut-
off scale, $\Lambda$, with all cases having $f_{a}=m_{Q}=10^{12}\,{\rm GeV}$.
The top row shows the SM, axion and heavy quark energy densities. On the left
column the x-axis is the inverse temperature, on the right, it is the
logarithm of the scale factor. Observe that with the maximum suppression
$M_{\rm Pl}$ for decay, the $Q$s get quite close to dominating, as anticipated
by our estimates in section 3. Therefore the $Q$ decays will produce some
meaningful impact on the SM radiation and dark radiation. We can see how the
decay of $Q$ changes the SM radiation energy density by observing that on the
top right panel the SM energy density in the case where ${\rm BR}_{\rm SM}=1$
is slightly higher than the ${\rm BR}_{\rm SM}=1/2$ case, because more of the
$Q$ energy is being pumped into the SM bath. On the left this effect can be
seen by noticing that the ${\rm BR}_{\rm SM}=1$ case has $\rho_{Q}$ depleting
at a faster rate, but in reality it is just the bath $T$ decreasing more
slowly. In the bottom row of Fig. 4 the energy density ratios of
$\rho_{a}/\rho_{\rm SM}$ are shown. The effect of $Q$ preferentially decaying
into SM particles can be seen by observing that the ${\rm BR}_{\rm SM}=1$ case
provides a substantially depleted ratio compared to the other two, which when
plugged into eq. (4.3) gives $\Delta N_{\rm eff}=0.0097$, under half of the
thermal, $d=4$ expectation. This situation is particularly relevant for many
$d=5$ models because many charge configurations produce, to the lowest order,
effective terms of the form

$\mathcal{L}_{Qf}^{d=5}=\overline{Q}_{\small L/R}\sigma_{\mu\nu}f_{\small
R/L}G^{\mu\nu}$ (4.9)

where $\sigma_{\mu\nu}=\frac{i}{2}\left[\gamma_{\mu},\gamma_{\nu}\right]$ and
$G^{\mu\nu}$ is the strong field strength tensor, and $f$ here is a SM fermion
including leptons. These models predict ${\rm BR_{SM}}=1$ and therefore a
$\Delta N_{\rm eff}<0.027$.

The ${\rm BR}_{\rm SM}=1-{\rm BR}_{a}=1/2$ case produces an enhanced value of
$\rho_{a}/\rho_{\rm SM}$ which corresponds to $\Delta N_{\rm eff}=1.37$.
Current limits from PLANCK set $\Delta N_{\rm eff}\leq 0.244$ using the more
stringent result (TT,TE,EE+low E) [56] which is much lower than the value
predicted in our example, the corresponding limit on decays into axions is
then ${\rm BR}_{a}\leq{0.12}$. With more sensitive probes expected such as the
CMB-HD survey, substantial progress will be made. From Ref. [82] we take the
projection which brings the sensitivity to $\Delta N_{\rm eff}=0.014$. For our
example with $m_{Q}=10^{12}\,{\rm GeV}$ and $\Lambda=M_{\rm Pl}$, these future
limits will probe ${\rm BR}_{a}$ down to $\sim 0.01$, an order of magnitude
improvement. As a reminder, we defined the ‘simplified’ branching ratios,
${\rm BR}_{{\rm SM}/a}$, in a model-independent way, such that a value of
${\rm BR}_{a}=1/2$ corresponds to axion models with charge assignments such
that there was only one decay $Q\rightarrow a+q$, i.e., half of the energy
carried from a decay of $Q$ goes into the axion energy density. Because $Q$ is
a fermion and the axion is a scalar ${\rm BR}_{a}$ should maximally be $1/2$.

In Fig. 5 we show how ${\rm BR}_{a}$ is constrained by current limits on
$\Delta N_{\rm eff}$ provided by PLANCK [56] for different values of $m_{Q}$
and by assuming $\Lambda=M_{\rm Pl}$. We select $m_{Q}$ values in the region
where $\Omega_{a}>\Omega_{\rm DM}$ to show how this limits the parameter space
if we relax the restrictions on initial angle and allow for $\theta_{i}\ll 1$.
It also shows the general behaviour of the limit. If a different $\Lambda$ was
taken we expect the constraint to translate in the x-axis in the direction of
lower $m_{Q}$ values.

Figure 5: Existing limits and future projections on ${\rm BR}_{a}$ for
different values of the dimension $d$ through which $Q$ decays proceed. Both
limits and projections are taken from measurements of the CMB, current limits
are shaded while future projections are denoted by dotted lines. We only show
the specific limits for $d=5$ (green) and $d=6$ (blue) assuming
$\Lambda=M_{\rm Pl}$. We show the general result when early matter domination
(EMD) is achieved with brown lines because this is unchanged for $d\geq 6$
models. We indicate where axions would be overproduced unless one introduces
some fine-tuned $\theta_{i}$ for each dimension with a dashed vertical line
and an arrow.

###  4.2 $\Delta N_{\rm eff}$ with early matter domination (EMD) ($d\geq 6$
Models)

As shown in section 3 preferred axion models that have $\mathcal{L}_{Qq}$
terms of the dimension $6$ and above can lead to a significant cosmological
phase of EMD. As long as this happens, we can analytically approximate the
implications $Q$ decay into axions will have for $\Delta N_{\rm eff}$.
Furthermore, the extra powers of $\Lambda$ which suppress the decays of heavy
quarks ensure that axions are not in thermal equilibrium by the time of decay.
Assuming instantaneous decays, we can write the SM and axion energy density
Boltzmann equations,

$\displaystyle\frac{\differential\rho_{\mathrm{R}}^{\rm SM}}{\differential
t}+4H\rho_{\mathrm{R}}^{\rm SM}$ $\displaystyle=$ $\displaystyle{\rm BR}_{\rm
SM}\Gamma_{Q}\rho_{Q}\,,$ (4.10a)
$\displaystyle\frac{\differential\rho_{a}}{\differential t}+4H\rho_{a}$
$\displaystyle=$ $\displaystyle{\rm BR}_{a}\Gamma_{Q}\rho_{Q}\,.$ (4.10b)

as in eq. 3.11 but we have dropped the axion annihilation term. Then by
rewriting in terms of comoving energy density $a^{4}\rho$, we are left with

$\frac{\differential\left(a^{4}\rho_{\mathrm{R}}^{\rm
SM}\right)}{\differential t}=\frac{{\rm BR}_{\rm SM}}{{\rm
BR}_{a}}\frac{\differential\left(a^{4}\rho_{a}\right)}{\differential t},$
(4.11)

using the fact that both $\Gamma_{Q}$ and the branching ratios are not time-
dependent. We integrate both sides of the equation, and use ${\rm BR}_{\rm
SM}=1-{\rm BR}_{a}$ to get

$\frac{\rho_{\mathrm{R}}^{\rm SM}}{\rho_{a}}=\frac{1-{\rm BR}_{a}}{{\rm
BR}_{a}}\,.$ (4.12)

With this relation we can use eq. (4.3) to calculate $\Delta N_{\rm eff}$ at
recombination. Note that in this approximation $\rho_{\rm R}^{\rm
SM}/\rho_{a}$ is set at time of $Q$ decay, $T_{\rm decay}^{Q}$, not $T_{\rm
BBN}$ as we do with our numerical calculation. The predicted $\Delta N_{\rm
eff}$ is not very sensitive to the precise value of $T_{\rm decay}^{Q}$
because the main variation comes from the temperature dependence of
$g_{\star}$ and $g_{\star S}$. There are three values of significance, $T_{\rm
decay}^{Q}\sim 100\,{\rm GeV}$, for temperatures above which there is a
minimal variation in the $g_{\star}$ values, $T_{\rm decay}^{Q}\sim T_{\rm
QCD}$, below which bound state formation may alter the number density of $Q$,
and $T_{\rm decay}^{Q}\sim T_{\rm BBN}$, below which constraints from BBN make
the scenario unviable. Current limits from PLANCK set $\Delta N_{\rm eff}\leq
0.244$, which corresponds to an upper limit on the axion branching ratio ${\rm
BR}_{a}\leq 0.077,\,0.057,\,0.038$ for $T_{\rm decay}^{Q}\ =100\,{\rm
GeV},\,T_{\rm QCD},\,T_{\rm BBN}$, respectively. For the CMB-HD survey, the
projection of a sensitivity down to $\Delta N_{\rm eff}=0.014$ [82] translates
to ${\rm BR}_{a}\leq 0.0048,\,0.0036,\,0.0023,$ for the same respective
temperatures. In Fig. 5 we show this approximate result for $Q$ domination
with $T_{\rm decay}^{Q}=100\,{\rm GeV}$ with a brown shaded region (currently
excluded) and brown dotted line (projection).

Additionally, in Fig. 5 we plot the result for $d=6$ and $\Lambda=M_{\rm Pl}$
obtained by fully solving the system numerically in blue. We do this to
compare our approximate result with the full calculation and see a remarkable
level of agreement. The constraints from $d=6$ do start to drop below our
approximation at around $\sim 5\times 10^{12}\,{\rm GeV}$, but this is because
$T_{\rm decay}^{Q}<100\,{\rm GeV}$ so there are fewer SM species in thermal
equilibrium.

For $d>6$ the results look similar to $d=6$, just the BBN constraint varies
according to $d$ and $\Lambda$, so we only show the approximate result for EMD
for these cases. As stressed above, $\Lambda$ is unknown and could be small
enough for all dimensions considered to allow for $m_{Q}\sim 10^{11}\,{\rm
GeV}$, see Fig. 3. Additionally, according to eq. 3.3 and Fig. 1 early $Q$
domination for this mass begins around $\sim 10^{5}\,{\rm GeV}$. So as long as
$m_{Q}\gtrsim 10^{6}\,{\rm GeV}$ such that freeze-out occurs above the
temperature of radiation and heavy quark equality, this limit may be
applicable. Although with our assumption that $f_{a}\sim m_{Q}$, experimental
constraints coming from supernovae [96, 97, 98] and neutron star cooling [99,
100] limit us to be above $m_{Q}\gtrsim 10^{9}\,{\rm GeV}$.

We also use our numerical solver to check what the predicted $\Delta N_{\rm
eff}$ is for the case when ${\rm BR_{SM}}=1$. Here because the dominating
quarks reheat the SM plasma and not the axions, the ratio $\rho_{a}/\rho_{\rm
R}^{\rm SM}$ is going to be extremely suppressed. We find the result can be
very low indeed, $\Delta N_{\rm eff}\approx 10^{-12}$, making such axions
virtually undetectable. Naively one assumes that if the axion is ever in
equilibrium is has at least $\Delta N_{\rm eff}\sim 0.027$, here we have an
example where the axion was thermally coupled but diluted substantially by the
decay of the heavy $Q$s.

The vertical dashed lines in Fig. 5 indicate where $\theta_{i}=0.5$ produces
too many axions via misalignment $\Omega_{a}>\Omega_{\rm DM}$ for $d=6,7,8,9$
in colors blue, red, purple and orange, respectively. These values correspond
to the points where the coloured shaded regions in Fig. 3 meet their relevant
BBN constraint.

##  5 Conclusions

In this article we explore the role that heavy quarks required in many axion
models may have on the evolution of the early Universe. In particular we find
that these heavy quarks can dominate in the early Universe and dilute the
relic cold axion abundance produced through misalignment. This has great
consequence for our understanding of preferred axion models, a subgroup of
axion models that meet desirable criteria. Assuming standard early cosmology,
the preferred axion models were believed to be limited to have heavy quark
decay at most by dimension-5 effective operators. We have demonstrated that
this is not the case and decays that proceed through effective operators up to
dimension 9 still survive all known constraints.

Our results here motivate the return to axion models where decays occur with
effective operators of dimension $d\geq 6$, since these should no longer be
considered disfavored. In this article we have not attempted to collate these
models in order to remain as general as possible. When we do study higher-
dimensional axion models, we may find more models that meet with more
restrictive criteria than considered here, such as no domain walls $N_{\rm
DW}=1$.

Figure 6: The extended region of preferred axion dark matter masses denoted by
the colored lines, the previously understood window is shown in green,
referring to decays via $d\leq 5$ effective operators and standard cosmology
(SC). Decays that proceed through higher-dimensional effective operators and
therefore produce a period of EMD are shown in blue ($d=6$), red ($d=7$),
purple ($d=8$) and orange ($d=9$). We have remained agnostic about the
suppression scale $\Lambda$. The maximum $m_{a}$ is the same for each case,
and given by the rightmost green marker. The existing limits and future
projections are taken from Ref. [101].

The effects of EMD on misalignment have been studied before, but the
remarkable feature of this work is that we have not invoked an arbitrary
matter field that dominates and decays before BBN, nothing additional needs to
be added beyond what is already present in the axion model. Instead the heavy
quarks required by many axion models actually produce this matter domination
and it is therefore inaccurate to make conclusions about which models to
discard without taking this into account. Importantly, this result has
substantial consequences for axion dark matter and its expected mass range. In
Fig. 6 we show how the preferred axion dark matter region is extended by this
new insight, and it extends by three orders of magnitude on the low $m_{a}$
end. An exciting consequence of this result is that a large area of this
parameter space is still viable but will be probed in the future by haloscope
experiments. These experiments probe the axion-photon coupling $g_{a\gamma}$,
thus we plot our results along the KSVZ line in the $|g_{a\gamma}|,$ vs.
$m_{a}$ plane. Our results provide a consistent and well-motivated class of
axion dark matter models to test. Furthermore, these mass values are achieved
without fine-tuning of the initial misalignment angle and are entirely
consistent with post-inflationary $U(1)_{\rm PQ}$ breaking.

The axion limits and projections of Fig. 6 have been plotted using Ref. [101]
and are indicated by dark and light gray shaded regions, respectively. The
limits come from microwave cavity experiments ADMX [10, 11] and HAYSTAC [12],
whereas the future projections comprise of results from a broad group of
operational or proposed experiments [13, 14, 15, 16], some of which are based
on the microwave cavity technique, i.e. CULTASK [102, 103], KLASH [85], ORGAN
[104], RADES [105], and QUAX [106]. Experiments using a different technology
but probing a similar mass range are MADMAX [107], BABYIAXO [108] and FLASH
[85]. For lower axion masses in the $10^{-9}$ to $10^{-7}$ eV range,
experiments such as SRF [87, 109], DMRadio [86], and ABRACADABRA [84] show
particular promise.

Since the matter-dominating field in this scenario may actually decay into
axions, the phase of early $Q$ domination will alter the expectation for the
axion contribution to the number of relativistic degrees of freedom. This
actually gives us another experimental window into preferred axion models with
$d\geq 5$. Probing the relative decays into axions and SM particles may allow
further discrimination of models. In particular, all $d=4$ models are expected
to reproduce the equilibrium result $\Delta N_{\rm eff}\approx 0.027$, but
many $d=5$ predict a ${\rm BR_{SM}}=1$ and predict a value less than half of
the equilibrium result, $\Delta N_{\rm eff}\approx 0.0097$. As mentioned
above, further study on $d\geq 6$ is required but models may indicate a
specific value of ${\rm BR_{SM}}$ which could already be constrained by
experiment or be in reach of future probes. In the event of discovering axion
dark matter in the non-standard mass window, this information would be pivotal
for determining whether $U(1)_{\rm PQ}$ breaking occurred before or after
inflation.

## Acknowledgments

AC would like to thank Dorian Amaral, Ken Mimasu, Matthew Kirk, Luca Di Luzio,
and Bradley Kavanagh for useful discussions. JO would like to thank Paola
Arias and Moira Venegas for discussions regarding thermally produced axions.
This work is supported by the grant “AstroCeNT: Particle Astrophysics Science
and Technology Centre" carried out within the International Research Agendas
programme of the Foundation for Polish Science financed by the European Union
under the European Regional Development Fund.

## References

  * [1] R.D. Peccei and H.R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * [2] R.D. Peccei and H.R. Quinn, _Some Aspects of Instantons_ , _Nuovo Cim. A_ 41 (1977) 309.
  * [3] J. Preskill, M.B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 127.
  * [4] M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137.
  * [5] L.F. Abbott and P. Sikivie, _A Cosmological Bound on the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 133.
  * [6] G.G. Raffelt, _Astrophysical methods to constrain axions and other novel particle phenomena_ , _Phys. Rept._ 198 (1990) 1.
  * [7] G.G. Raffelt, _Astrophysical axion bounds_ , _Lect. Notes Phys._ 741 (2008) 51 [hep-ph/0611350].
  * [8] L. Di Luzio, M. Giannotti, E. Nardi and L. Visinelli, _The landscape of QCD axion models_ , _Phys. Rept._ 870 (2020) 1 [2003.01100].
  * [9] P. Sikivie, _Experimental Tests of the Invisible Axion_ , _Phys. Rev. Lett._ 51 (1983) 1415.
  * [10] ADMX collaboration, _A SQUID-based microwave cavity search for dark-matter axions_ , _Phys. Rev. Lett._ 104 (2010) 041301 [0910.5914].
  * [11] ADMX collaboration, _Extended Search for the Invisible Axion with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 124 (2020) 101303 [1910.08638].
  * [12] HAYSTAC collaboration, _Results from phase 1 of the HAYSTAC microwave cavity axion experiment_ , _Phys. Rev. D_ 97 (2018) 092001 [1803.03690].
  * [13] I.G. Irastorza and J. Redondo, _New experimental approaches in the search for axion-like particles_ , _Prog. Part. Nucl. Phys._ 102 (2018) 89 [1801.08127].
  * [14] P. Sikivie, _Invisible Axion Search Methods_ , _Rev. Mod. Phys._ 93 (2021) 015004 [2003.02206].
  * [15] J. Billard et al., _Direct detection of dark matter—APPEC committee report_ , _Rept. Prog. Phys._ 85 (2022) 056201 [2104.07634].
  * [16] C.B. Adams et al., _Axion Dark Matter_ , in _Snowmass 2021_ , 3, 2022 [2203.14923].
  * [17] L. Di Luzio, F. Mescia and E. Nardi, _Redefining the Axion Window_ , _Phys. Rev. Lett._ 118 (2017) 031801 [1610.07593].
  * [18] L. Di Luzio, F. Mescia and E. Nardi, _Window for preferred axion models_ , _Phys. Rev. D_ 96 (2017) 075003 [1705.05370].
  * [19] P.J. Steinhardt and M.S. Turner, _Saving the Invisible Axion_ , _Phys. Lett. B_ 129 (1983) 51.
  * [20] G. Lazarides, R.K. Schaefer, D. Seckel and Q. Shafi, _Dilution of Cosmological Axions by Entropy Production_ , _Nucl. Phys. B_ 346 (1990) 193.
  * [21] M. Kawasaki, T. Moroi and T. Yanagida, _Can decaying particles raise the upper bound on the Peccei-Quinn scale?_ , _Phys. Lett. B_ 383 (1996) 313 [hep-ph/9510461].
  * [22] N. Ramberg and L. Visinelli, _Probing the Early Universe with Axion Physics and Gravitational Waves_ , _Phys. Rev. D_ 99 (2019) 123513 [1904.05707].
  * [23] P. Arias, N. Bernal, D. Karamitros, C. Maldonado, L. Roszkowski and M. Venegas, _New opportunities for axion dark matter searches in nonstandard cosmological models_ , _JCAP_ 11 (2021) 003 [2107.13588].
  * [24] P. Arias, N. Bernal, J.K. Osiński and L. Roszkowski, _Dark matter axions in the early universe with a period of increasing temperature_ , _JCAP_ 05 (2023) 028 [2207.07677].
  * [25] J.E. Kim and G. Carosi, _Axions and the Strong CP Problem_ , _Rev. Mod. Phys._ 82 (2010) 557 [0807.3125].
  * [26] J.E. Kim, _Weak Interaction Singlet and Strong CP Invariance_ , _Phys. Rev. Lett._ 43 (1979) 103.
  * [27] M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, _Can Confinement Ensure Natural CP Invariance of Strong Interactions?_ , _Nucl. Phys. B_ 166 (1980) 493.
  * [28] A.R. Zhitnitsky, _On Possible Suppression of the Axion Hadron Interactions. (In Russian)_ , _Sov. J. Nucl. Phys._ 31 (1980) 260.
  * [29] M. Dine, W. Fischler and M. Srednicki, _A Simple Solution to the Strong CP Problem with a Harmless Axion_ , _Phys. Lett. B_ 104 (1981) 199.
  * [30] 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].
  * [31] K. Jedamzik, _Big bang nucleosynthesis constraints on hadronically and electromagnetically decaying relic neutral particles_ , _Phys. Rev. D_ 74 (2006) 103509 [hep-ph/0604251].
  * [32] K. Jedamzik, _Bounds on long-lived charged massive particles from Big Bang nucleosynthesis_ , _JCAP_ 03 (2008) 008 [0710.5153].
  * [33] M. Kawasaki, K. Kohri, T. Moroi and Y. Takaesu, _Revisiting Big-Bang Nucleosynthesis Constraints on Long-Lived Decaying Particles_ , _Phys. Rev. D_ 97 (2018) 023502 [1709.01211].
  * [34] V. De Luca, A. Mitridate, M. Redi, J. Smirnov and A. Strumia, _Colored Dark Matter_ , _Phys. Rev. D_ 97 (2018) 115024 [1801.01135].
  * [35] C. Gross, A. Mitridate, M. Redi, J. Smirnov and A. Strumia, _Cosmological Abundance of Colored Relics_ , _Phys. Rev. D_ 99 (2019) 016024 [1811.08418].
  * [36] C.B. Dover, T.K. Gaisser and G. Steigman, _COSMOLOGICAL CONSTRAINTS ON NEW STABLE HADRONS_ , _Phys. Rev. Lett._ 42 (1979) 1117.
  * [37] K. Griest and M. Kamionkowski, _Unitarity Limits on the Mass and Radius of Dark Matter Particles_ , _Phys. Rev. Lett._ 64 (1990) 615.
  * [38] E. Nardi and E. Roulet, _Are exotic stable quarks cosmologically allowed?_ , _Phys. Lett. B_ 245 (1990) 105.
  * [39] A. Arvanitaki, C. Davis, P.W. Graham, A. Pierce and J.G. Wacker, _Limits on split supersymmetry from gluino cosmology_ , _Phys. Rev. D_ 72 (2005) 075011 [hep-ph/0504210].
  * [40] J. Kang, M.A. Luty and S. Nasri, _The Relic abundance of long-lived heavy colored particles_ , _JHEP_ 09 (2008) 086 [hep-ph/0611322].
  * [41] C. Jacoby and S. Nussinov, _The Relic Abundance of Massive Colored Particles after a Late Hadronic Annihilation Stage_ , 0712.2681.
  * [42] M. Kusakabe and T. Takesako, _Resonant annihilation of long-lived massive colored particles through hadronic collisions_ , _Phys. Rev. D_ 85 (2012) 015005 [1112.0860].
  * [43] M.L. Perl, P.C. Kim, V. Halyo, E.R. Lee, I.T. Lee, D. Loomba et al., _The Search for stable, massive, elementary particles_ , _Int. J. Mod. Phys. A_ 16 (2001) 2137 [hep-ex/0102033].
  * [44] M.P. Hertzberg and A. Masoumi, _Astrophysical Constraints on Singlet Scalars at LHC_ , _JCAP_ 04 (2017) 028 [1607.06445].
  * [45] A. Gould, B.T. Draine, R.W. Romani and S. Nussinov, _Neutron Stars: Graveyard of Charged Dark Matter_ , _Phys. Lett. B_ 238 (1990) 337.
  * [46] G.D. Mack, J.F. Beacom and G. Bertone, _Towards Closing the Window on Strongly Interacting Dark Matter: Far-Reaching Constraints from Earth’s Heat Flow_ , _Phys. Rev. D_ 76 (2007) 043523 [0705.4298].
  * [47] R. Allahverdi et al., _The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe_ , 2006.16182.
  * [48] M. Kawasaki, K. Kohri and N. Sugiyama, _MeV scale reheating temperature and thermalization of neutrino background_ , _Phys. Rev. D_ 62 (2000) 023506 [astro-ph/0002127].
  * [49] S. Hannestad, _What is the lowest possible reheating temperature?_ , _Phys. Rev. D_ 70 (2004) 043506 [astro-ph/0403291].
  * [50] K. Ichikawa, M. Kawasaki and F. Takahashi, _The Oscillation effects on thermalization of the neutrinos in the Universe with low reheating temperature_ , _Phys. Rev. D_ 72 (2005) 043522 [astro-ph/0505395].
  * [51] K. Ichikawa, M. Kawasaki and F. Takahashi, _Constraint on the Effective Number of Neutrino Species from the WMAP and SDSS LRG Power Spectra_ , _JCAP_ 05 (2007) 007 [astro-ph/0611784].
  * [52] P.F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti, _Bounds on very low reheating scenarios after Planck_ , _Phys. Rev. D_ 92 (2015) 123534 [1511.00672].
  * [53] T. Hasegawa, N. Hiroshima, K. Kohri, R.S.L. Hansen, T. Tram and S. Hannestad, _MeV-scale reheating temperature and thermalization of oscillating neutrinos by radiative and hadronic decays of massive particles_ , _JCAP_ 12 (2019) 012 [1908.10189].
  * [54] P. Sikivie, _Axion Cosmology_ , _Lect. Notes Phys._ 741 (2008) 19 [astro-ph/0610440].
  * [55] D.J.E. Marsh, _Axion Cosmology_ , _Phys. Rept._ 643 (2016) 1 [1510.07633].
  * [56] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [57] M. Kamionkowski and J. March-Russell, _Planck scale physics and the Peccei-Quinn mechanism_ , _Phys. Lett. B_ 282 (1992) 137 [hep-th/9202003].
  * [58] R. Holman, S.D.H. Hsu, T.W. Kephart, E.W. Kolb, R. Watkins and L.M. Widrow, _Solutions to the strong CP problem in a world with gravity_ , _Phys. Lett. B_ 282 (1992) 132 [hep-ph/9203206].
  * [59] S.M. Barr and D. Seckel, _Planck scale corrections to axion models_ , _Phys. Rev. D_ 46 (1992) 539.
  * [60] A. Ringwald and K. Saikawa, _Axion dark matter in the post-inflationary Peccei-Quinn symmetry breaking scenario_ , _Phys. Rev. D_ 93 (2016) 085031 [1512.06436].
  * [61] M.P. Hertzberg, M. Tegmark and F. Wilczek, _Axion Cosmology and the Energy Scale of Inflation_ , _Phys. Rev. D_ 78 (2008) 083507 [0807.1726].
  * [62] T.W.B. Kibble, _Topology of Cosmic Domains and Strings_ , _J. Phys. A_ 9 (1976) 1387.
  * [63] T.W.B. Kibble, G. Lazarides and Q. Shafi, _Walls Bounded by Strings_ , _Phys. Rev. D_ 26 (1982) 435.
  * [64] A. Vilenkin and E.P.S. Shellard, _Cosmic Strings and Other Topological Defects_ , Cambridge University Press (7, 2000).
  * [65] C. Hagmann, S. Chang and P. Sikivie, _Axion radiation from strings_ , _Phys. Rev. D_ 63 (2001) 125018 [hep-ph/0012361].
  * [66] O. Wantz and E.P.S. Shellard, _Axion Cosmology Revisited_ , _Phys. Rev. D_ 82 (2010) 123508 [0910.1066].
  * [67] T. Hiramatsu, M. Kawasaki, T. Sekiguchi, M. Yamaguchi and J. Yokoyama, _Improved estimation of radiated axions from cosmological axionic strings_ , _Phys. Rev. D_ 83 (2011) 123531 [1012.5502].
  * [68] M. Kawasaki, K. Saikawa and T. Sekiguchi, _Axion dark matter from topological defects_ , _Phys. Rev. D_ 91 (2015) 065014 [1412.0789].
  * [69] M. Gorghetto and G. Villadoro, _Topological Susceptibility and QCD Axion Mass: QED and NNLO corrections_ , _JHEP_ 03 (2019) 033 [1812.01008].
  * [70] M. Buschmann, J.W. Foster, A. Hook, A. Peterson, D.E. Willcox, W. Zhang et al., _Dark matter from axion strings with adaptive mesh refinement_ , _Nature Commun._ 13 (2022) 1049 [2108.05368].
  * [71] P. Sikivie, _Of Axions, Domain Walls and the Early Universe_ , _Phys. Rev. Lett._ 48 (1982) 1156.
  * [72] A. Vilenkin and A.E. Everett, _Cosmic Strings and Domain Walls in Models with Goldstone and PseudoGoldstone Bosons_ , _Phys. Rev. Lett._ 48 (1982) 1867.
  * [73] S.M. Barr, K. Choi and J.E. Kim, _Some aspects of axion cosmology in unified and superstring models_ , _Nucl. Phys. B_ 283 (1987) 591.
  * [74] M.S. Turner, _Cosmic and Local Mass Density of Invisible Axions_ , _Phys. Rev. D_ 33 (1986) 889.
  * [75] K. Choi, H.B. Kim and J.E. Kim, _Axion cosmology with a stronger QCD in the early universe_ , _Nucl. Phys. B_ 490 (1997) 349 [hep-ph/9606372].
  * [76] G. Alonso-Álvarez, J.M. Cline and T. Xiao, _The Flavor of QCD Axion Dark Matter_ , 2305.00018.
  * [77] L. Di Luzio, R. Gröber, J.F. Kamenik and M. Nardecchia, _Accidental matter at the LHC_ , _JHEP_ 07 (2015) 074 [1504.00359].
  * [78] S. Borsanyi et al., _Calculation of the axion mass based on high-temperature lattice quantum chromodynamics_ , _Nature_ 539 (2016) 69 [1606.07494].
  * [79] F. D’Eramo, F. Hajkarim and S. Yun, _Thermal QCD Axions across Thresholds_ , _JHEP_ 10 (2021) 224 [2108.05371].
  * [80] G.F. Giudice, E.W. Kolb and A. Riotto, _Largest temperature of the radiation era and its cosmological implications_ , _Phys. Rev. D_ 64 (2001) 023508 [hep-ph/0005123].
  * [81] D. Karamitros, _MiMeS: Misalignment mechanism solver_ , _Comput. Phys. Commun._ 275 (2022) 108311 [2110.12253].
  * [82] CMB-HD collaboration, _Snowmass2021 CMB-HD White Paper_ , 2203.05728.
  * [83] I. Stern, _ADMX Status_ , _PoS_ ICHEP2016 (2016) 198 [1612.08296].
  * [84] C.P. Salemi et al., _Search for Low-Mass Axion Dark Matter with ABRACADABRA-10 cm_ , _Phys. Rev. Lett._ 127 (2021) 081801 [2102.06722].
  * [85] D. Alesini, D. Babusci, D. Di Gioacchino, C. Gatti, G. Lamanna and C. Ligi, _The KLASH Proposal_ , 1707.06010.
  * [86] DMRadio collaboration, _Projected sensitivity of DMRadio-m3: A search for the QCD axion below 1 $\mu$eV_, _Phys. Rev. D_ 106 (2022) 103008 [2204.13781].
  * [87] A. Berlin, R.T. D’Agnolo, S.A.R. Ellis and K. Zhou, _Heterodyne broadband detection of axion dark matter_ , _Phys. Rev. D_ 104 (2021) L111701 [2007.15656].
  * [88] P.F. de Salas and S. Pastor, _Relic neutrino decoupling with flavour oscillations revisited_ , _JCAP_ 07 (2016) 051 [1606.06986].
  * [89] P. Graf and F.D. Steffen, _Thermal axion production in the primordial quark-gluon plasma_ , _Phys. Rev. D_ 83 (2011) 075011 [1008.4528].
  * [90] A. Salvio, A. Strumia and W. Xue, _Thermal axion production_ , _JCAP_ 01 (2014) 011 [1310.6982].
  * [91] P. Arias, N. Bernal, J.K. Osiński, L. Roszkowski and M. Venegas, _Revisiting signatures of thermal axions in nonstandard cosmologies_ , 2308.01352.
  * [92] L.J. Hall, K. Jedamzik, J. March-Russell and S.M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080 [0911.1120].
  * [93] F. D’Eramo, N. Fernandez and S. Profumo, _Dark Matter Freeze-in Production in Fast-Expanding Universes_ , _JCAP_ 02 (2018) 046 [1712.07453].
  * [94] Y. Du, F. Huang, H.-L. Li, Y.-Z. Li and J.-H. Yu, _Revisiting dark matter freeze-in and freeze-out through phase-space distribution_ , _JCAP_ 04 (2022) 012 [2111.01267].
  * [95] F. D’Eramo, F. Hajkarim and S. Yun, _Thermal Axion Production at Low Temperatures: A Smooth Treatment of the QCD Phase Transition_ , _Phys. Rev. Lett._ 128 (2022) 152001 [2108.04259].
  * [96] J.W. Brockway, E.D. Carlson and G.G. Raffelt, _SN1987A gamma-ray limits on the conversion of pseudoscalars_ , _Phys. Lett. B_ 383 (1996) 439 [astro-ph/9605197].
  * [97] J.A. Grifols, E. Masso and R. Toldra, _Gamma-rays from SN1987A due to pseudoscalar conversion_ , _Phys. Rev. Lett._ 77 (1996) 2372 [astro-ph/9606028].
  * [98] A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi and A. Ringwald, _Revisiting the SN1987A gamma-ray limit on ultralight axion-like particles_ , _JCAP_ 02 (2015) 006 [1410.3747].
  * [99] K. Hamaguchi, N. Nagata, K. Yanagi and J. Zheng, _Limit on the Axion Decay Constant from the Cooling Neutron Star in Cassiopeia A_ , _Phys. Rev. D_ 98 (2018) 103015 [1806.07151].
  * [100] L.B. Leinson, _Impact of axions on the Cassiopea A neutron star cooling_ , _JCAP_ 09 (2021) 001 [2105.14745].
  * [101] C. O’HARE, _cajohare/axionlimits: Axionlimits_ , July, 2020.  10.5281/zenodo.3932430.
  * [102] S. Lee, S. Ahn, J. Choi, B.R. Ko and Y.K. Semertzidis, _Axion Dark Matter Search around 6.7 $\mu$eV_, _Phys. Rev. Lett._ 124 (2020) 101802 [2001.05102].
  * [103] CAPP collaboration, _First Results from an Axion Haloscope at CAPP around 10.7 $\mu$eV_, _Phys. Rev. Lett._ 126 (2021) 191802 [2012.10764].
  * [104] B.T. McAllister, G. Flower, J. Kruger, E.N. Ivanov, M. Goryachev, J. Bourhill et al., _The ORGAN Experiment: An axion haloscope above 15 GHz_ , _Phys. Dark Univ._ 18 (2017) 67 [1706.00209].
  * [105] A.A. Melcón et al., _Axion Searches with Microwave Filters: the RADES project_ , _JCAP_ 05 (2018) 040 [1803.01243].
  * [106] D. Alesini et al., _Search for invisible axion dark matter of mass m ${}_{a}=43~{}\mu$eV with the QUAX–$a\gamma$ experiment_, _Phys. Rev. D_ 103 (2021) 102004 [2012.09498].
  * [107] MADMAX collaboration, _A new experimental approach to probe QCD axion dark matter in the mass range above 40 $\mu$eV_, _Eur. Phys. J. C_ 79 (2019) 186 [1901.07401].
  * [108] S. Ahyoune et al., _A proposal for a low-frequency axion search in the 1-2 $\mu$eV range and below with the BabyIAXO magnet_, 2306.17243.
  * [109] Z. Tang et al., _SRF Cavity Searches for Dark Photon Dark Matter: First Scan Results_ , 2305.09711.


# CC-GPX: Extracting High-Quality Annotated Geospatial Data from Common Crawl

Ilya Ilyankou<EMAIL_ADDRESS>0009-0008-7082-7122 ,  Meihui Wang
<EMAIL_ADDRESS>0000-0001-8420-2141 ,  James Haworth
<EMAIL_ADDRESS>0000-0001-9506-4266 UCL SpaceTimeLabLondonUK and  Stefano
Cavazzi<EMAIL_ADDRESS>0000-0003-3575-0365 Ordnance
SurveySouthamptonUK

(2024; 17 May 2024)

###### Abstract.

The Common Crawl (CC) corpus is the largest open web crawl dataset containing
9.5+ petabytes of data captured since 2008. The dataset is instrumental in
training large language models, and as such it has been studied for
(un)desirable content, and distilled for smaller, domain-specific datasets.
However, to our knowledge, no research has been dedicated to using CC as a
source of annotated geospatial data. In this paper, we introduce an efficient
pipeline to extract annotated user-generated tracks from GPX files found in
CC, and the resulting multimodal dataset with 1,416 pairings of human-written
descriptions and MultiLineString vector data from the 6 most recent CC
releases. The dataset can be used to study people’s outdoor activity patterns,
the way people talk about their outdoor experiences, and for developing
trajectory generation or track annotation models. Our reproducible code is
available on GitHub111https://github.com/ilyankou/cc-gpx.

Common Crawl, GPS, GPX, GIS, hiking, user-generated routes

††copyright: acmlicensed††journalyear: 2024††doi: XXXXXXX.XXXXXXX††isbn:
978-1-4503-XXXX-X/18/06††ccs: Information systems Geographic information
systems††ccs: Information systems Web mining††ccs: Computing methodologies
Natural language processing††ccs: Human-centered computing Geographic
visualization

##  1\. Introduction

The Common Crawl (CC) corpus222https://commoncrawl.org/ is the largest
publicly available web crawl dataset containing 9.5+ petabytes of data dating
back from 2008 (Baack, 2024). Due to its size, the dataset is widely used to
train large language models, including GPT-3 (Thompson, 2022).

In the past, researchers have studied the contents of Common Crawl for
undesirable content (Luccioni and Viviano, 2021), extracted documents to
create the largest Spanish-language crawling corpus (Gutiérrez-Fandiño et al.,
2022), extracted parallel texts for machine translation (Smith et al., 2013),
PDF documents (Turski et al., 2023) and Word documents with layout annotations
(Weber et al., 2023), and built the ‘Colossal Clean Crawled Corpus’ (C4)
dataset333https://paperswithcode.com/dataset/c4 among other things. However,
to our knowledge, no research was dedicated to analysing the geospatial
contents or extracting annotated geospatial data from CC.

While websites and apps such as Strava444https://www.strava.com/,
AllTrails555https://www.alltrails.com/, and
Outdooractive666https://www.outdooractive.com/en/ contain large collections of
curated and user-generated routes that are often accompanied by textual
descriptions, these typically come with substantial licensing restrictions.
CC, on the other hand, provides free access to all its data collected in
accordance with website policies (subject to fair use restrictions), and
exposes researchers to niche and varied websites that are harder to come
across any other way.

In this paper, we introduce an efficient pipeline to identify, download, and
clean GPX files that contain user-generated tracks, such as those produced by
recording hiking, running, or cycling activities, along with human-written
textual descriptions of those tracks. We run our process on the six most
recent CC releases (spanning just over a year) to generate a multi-lingual,
multi-activity dataset containing over 1,400 pairs of annotated
MultiLineString features. When extended to all CC data back to 2008, the
dataset is likely to contain over 10,000 samples.

The resulting dataset can be used to study the way people talk about outdoor
activities, to facilitate the generation of human-like descriptions of
activity tracks, or to study outdoor activity trajectories.

##  2\. Data Collection and Processing

Hikers, cyclists, and in particular runners often track their activities using
GPS-enabled mobile devices, such as Garmin smart watches (Karahanoğlu et al.,
2021), and sometimes add detailed annotations describing their experience.
Alternatively, people may plan their active journey beforehand using websites
such as cycle.travel, and add notes before or after the journey. In this
paper, we decided to focus on the popular XML-based
.gpx777http://www.topografix.com/GPX/1/1/ (GPX Exchange Format) files that are
most commonly used to share the activities, and ignore less popular formats
such as .fit or .tcx which are very rare in Common Crawl.

###  2.1. Downloading GPX files from CC

We use index tables from the 6 most recent releases, CC-MAIN-2023-* and CC-
MAIN-2024-10, to identify GPX files. Each index table contains a list of URLs
and MIME-types of all web pages and files scraped by the CCBot in that
release, and is broken into 300 parts. We use duckdb888https://duckdb.org/ to
query each part to locate files whose detected MIME type is ilike ‘%gpx%’ or
whose filename extension is ‘.gpx’ (typically both conditions are true at
once).

For each GPX file, we record its Web Archive (WARC) file location, together
with the offset and length in bytes. The offset and length indicate the
location of the GPX file within a WARC file. This allows us to use Python’s
requests999https://github.com/psf/requests library and the Range request
header parameter to download individual GPX files from CC without the need to
download large WARC files first:

    requests.get(
        f’https://data.commoncrawl.org/{warc_file}’
        headers={
            ’Range’: f’bytes={offset}-{offset+length-1}’
        })

This massively speeds up the data acquisition process; in fact, it takes
approximately 40 minutes to download all GPX files from a single CC release,
which typically contains around 90 terabytes of compressed data spread across
3 billion files. We identified 112,953 GPX files across the six most recent CC
releases. We were able to successfully download 111,102 files (98.4% of those
identified). Among downloaded files, 102,103 (or 91.9%) came from unique URLs,
indicating that very few GPX files appear in more than one CC release. After
removing duplicate GPX files that come from different URLs, we are left with
94,170 GPX files.

###  2.2. Identifying high-quality tracks

We use gpxpy101010https://github.com/tkrajina/gpxpy Python library to parse
response strings into GPX objects that can be analysed. A typical GPX file
consists of tracks; tracks consist of segments; segments consist of points.
GPX files representing simple recorded point-to-point activities, such as
hiking or cycling trips, would normally contain a single track that consists
of one or more segments (for example, a runner may choose to record each lap
as a separate segment). Multi-track GPX files that we inspected were typically
used to record multi-day activities such as races around Europe. We decided to
keep only single-track GPX files that represent activities no longer than 100
km (62.1 mi) because we believe longer activities cannot be adequately
described in a few paragraphs, and are unlikely to be relevant to most people.
We also remove activities shorter than 0.5 km (0.31 mi).

We keep tracks that have at least 1 GPS point per 100 m (328 ft) on average
($\sim$87% of all tracks satisfy this constraint). The threshold was
determined empirically by analysing individual GPX files.

###  2.3. Identifying activity descriptions

We extract descriptions from each track’s <desc> tag. We use
BeautifulSoup111111https://www.crummy.com/software/BeautifulSoup/ Python
library to remove occasional HTML tags present in text, replace characters
such as newlines and tabs with single spaces, and use regular expressions to
remove text inside square and curly brackets (tags likely added by some apps).

We only keep the tracks that have an associated description between 50 and
2000 characters long after the manipulations above. The upper character limit
represents $\sim$$99^{th}$ percentile of the description length distribution.

To identify high-quality track and activity descriptions, we employ Llama-3 8B
Instruct, one of the most capable open-source language models as of April 2024
(AI, 2024). We use the prompt Does the text in triple quotes represent a high-
quality and insightful route or track description, or an activity description
such as hiking, cycling, or racing? Respond with ‘True’ or ‘False’. If you are
unsure, say ‘False’. Text: “‘{text}”’, and remove tracks for which False is
returned.

Positive examples (tracks included in the dataset):

  * •

‘A lovely 4 hour walk from Ockley railway station to Dorking, via the
beautiful view from Leith Hill Tower.’

  * •

‘Quieter roads and backstreets, quirky interest but still direct.’

  * •

‘This trail leading to the borough of Kukleny takes in significant industrial
buildings evoking the glory of the iron and leather industries, as well as
family houses, a funeral hall, and a Cubist vocational school.’

Negative examples (tracks excluded from the dataset):

  * •

‘Cheese-making is a centuries-old tradition in Gruyère and firmly rooted in
the heritage. Immerse yourself in a walk with a historic feel.’

  * •

‘File with points/tracks from Locus Map Classic/3.65.2’

  * •

‘Note check opening times of Manor Park Cremitorium’

###  2.4. Removing personal information

The GPX files available in CC would typically be uploaded and openly shared by
the users in the open web. However, it is our responsibility as researchers to
respect user privacy and remove any information that can be deemed personal.

We use regular expressions to mask emails and URLs (very occasionally found in
descriptions) with tokens ¡EMAIL¿ and ¡URL¿. We also identify phone numbers
using Google’s libphonenumber121212https://github.com/google/libphonenumber
Python library, and mask those with ¡TELEPHONE¿. We also remove timestamps
from points due to both privacy and data quality reasons (many tracks lack
timestamps, and, unlike missing elevation data, timestamps are impossible to
recreate).

We then use Llama-3 8B Instruct with the prompt Does the text in triple quotes
contain any personally identifiable information, such as someone’s address or
name? Respond with ‘True’ or ‘False’. If you are unsure, say ‘True’. Text:
“‘{text}”’ to further identify and remove tracks whose descriptions may
contain personal information.

###  2.5. Translating descriptions to English

The majority of GPX files we identified are from European websites, and their
descriptions are not in English. We use
pycld2131313https://github.com/aboSamoor/pycld2 Python library to identify the
original language. We remove tracks whose language is identified as Unknown as
those often contain majority numbers or incoherent text, as well as all tracks
whose descriptions are in languages that have 5 or fewer samples.

Initially, we experimented with several popular open-source LLMs, including
Mistral and the Llamas (2 & 3), to translate descriptions into English, a task
LLMs are shown to be good at (Zhu et al., 2023). Unfortunately, we discovered
that LLMs occasionally modify formatting and add unnecessary context (e.g.,
‘The text is in Ukrainian. Here’s the English translation: …’). Extensive
prompt-engineering was not able to remedy these issues. We ultimately chose to
use argos-translate141414https://github.com/argosopentech/argos-translate, a
popular offline Python translation library based on
OpenNMT151515https://opennmt.net/, an open source neural machine translation
system, to translate all non-English descriptions into English.

###  2.6. Adding missing elevation

Just over half of all relevant GPX files contained device-recorded point
elevation data. To produce a consistent dataset with all tracks represented by
3D points (lat, lon, elevation), we use the Shuttle Radar Topography
Mission161616https://www.earthdata.nasa.gov/sensors/srtm digital elevation
model, acquired via the SRTM.py171717https://github.com/tkrajina/srtm.py
Python library, to approximate point elevations in tracks where elevation data
is not originally available. We note the source of the elevation data, GPS
(device-recorded) or DEM (SRTM), in the elev_source field (see Table 1).

##  3\. Dataset

The final dataset consists of 1,416 tracks (or around 1.5% of the initially
deduplicated GPX files) and accompanying descriptions in both the original
language and English, along with other properties described in Table 1. The
examples of two tracks with descriptions are shown in Figures 1 and 2.

Figure 1. A 18.2 km (11.3 mi) route in the UK. The description reads:
‘Longbridge to Whalley Slowway following part of the Ribble Way. Difficult to
find a good crossing of the A59. The crossing chosen crosses the road from
footpath to footpath in a place with good visibility. The road
junctions/bridges were actually worse as would need to walk along a fast road
with no pavement rather than just cross once at right angles. This crossing
sets up good sections without roads. Good spacing of waypoints at Old Langho
and Ribchester.’

Example of a linear route from the dataset.

Figure 2. A 13.8 km (8.6 mi) circular route in Germany. The description in
German reads: ‘Der Weg ist sehr gut gekennzeichnet mit einem schwarzen
Hirschkäfer (Hootzemann) auf weißem Grund. Mein Start- und Zielpunkt war das
Schützenhaus Eiweiler in der Nähe der Großwald-Brauerei.’ The English
translation is: ‘The path is very well marked with a black deer beetle
(Hootzemann) on white ground. My starting and finishing point was the
Schützenhaus Eiweiler near the Großwald brewery.’

Example of a circular route from the dataset.

The GPX files come from 135 unique domain names. Track descriptions come in 11
languages, with the majority being French (802), German (330), and English
(110). The original descriptions are between 50 and 1999 characters long (for
translated, between 34 and 1866), with a median of 303 characters (for
translated, 285.5). The tracks originate in 25 countries, with most popular
being France (710), Austria (164), Germany (133), Switzerland (116), and Italy
(81). The shortest track is 607 m (1991 ft) while the longest is 99.3 km (61.7
mi). The average track length is 20.5 km (12.7 mi), while the median is much
shorter at 12.1 km (7.5 mi). Histograms of route length (in km) and
description length (in characters) are shown in Figure 3.

Figure 3. Select dataset properties.

Select dataset properties.

Table 1. Dataset property descriptions # | Property |  Description
---|---|---
1 | url |  URL of the GPX file
2 | warc_file |  CC WARC file with GPX file
3 | warc_offset |  GPX file position in WARC
4 | warc_len |  GPX file byte length
5 | country |  Country name as determined by the first point in the track intersecting boundaries from https://www.geoboundaries.org/
6 | desc |  Original track description
7 | desc_lang |  Track description language code, as determined by pycld2
8 | desc_en |  Track description translated into English
9 | elev_source |  GPS if elevation is recorded by device; DEM if determined later from Shuttle Radar Topography Mission
10 | elev_highest |  Track’s highest point, m
11 | elev_lowest |  Track’s lowest point, m
12 | uphill |  Cumulative elevation gain, m
13 | downhill |  Cumulative elevation loss, m
14 | length_2d |  Track length disregarding elevation, m
15 | length_3d |  Track length accounting for elevation, m
16 | is_circular |  True if start and end points are within 350 m from each other, False otherwise
17 | geometry |  MultiLineString Z geometry in GPS coordinates: (lat, lon, elevation)

##  4\. Applications

The tracks in the resulting dataset represent real-life outdoor activities
that were either recorded (i.e., completed) or planned using GIS software.
Thus, these routes can be used in training of GPS trajectory generation models
that are focused on recreational activities.

The description-MultiLineString pairs can be used to study how people describe
their outdoor activities, such as what kinds of landmarks, landscapes, and
experiences they choose to mention or ignore, and get an insight into turn-by-
turn decisions of users. The immediate application can be to fine-tine
language models to generate human-like descriptions of existing or
synthetically composed routes. This dataset can become especially powerful
when combined with Point of Interest (POI) and/or road segment datasets, such
as those from OpenStreetMap181818https://welcome.openstreetmap.org/ or
Overture191919https://overturemaps.org/.

##  5\. Conclusion

In this paper, we demonstrated that Common Crawl can be used as a novel source
of annotated geospatial data. As a case study, we built an efficient pipeline
to extract annotated outdoor activity tracks from GPX files that can be
executed on a medium-powered laptop in a matter of hours. The resulting
multimodal dataset contains textual descriptions of outdoor activities paired
with MultiLineString vector geospatial data and can be used to study people’s
outdoor activity habits, and the relationship between what people experience
and what they say they experience. The dataset can be useful as a fine-tuning
set for various GeoAI models, from trajectory to human-like description
generators. We recognise that our case study captures just a fraction of all
annotated geospatial data available in CC, and we intend to extend the scope
in future work.

###### Acknowledgements.

This work was supported by the Ordnance Survey, and the Engineering and
Physical Sciences Research Council [grant no. EP/Y528651/1].

## References

  * (1)
  * AI (2024) Meta AI. 2024.  Introducing Meta Llama 3: The most capable openly available LLM to date.  https://ai.meta.com/blog/meta-llama-3/
  * Baack (2024) Stefan Baack. 2024.  Training Data for the Price of a Sandwich: Common Crawl’s Impact on Generative AI.  (Feb. 2024).
  * Gutiérrez-Fandiño et al. (2022) Asier Gutiérrez-Fandiño, David Pérez-Fernández, Jordi Armengol-Estapé, David Griol, and Zoraida Callejas. 2022.  esCorpius: A Massive Spanish Crawling Corpus.  http://arxiv.org/abs/2206.15147 arXiv:2206.15147 [cs].
  * Karahanoğlu et al. (2021) Armağan Karahanoğlu, Rúben Gouveia, Jasper Reenalda, and Geke Ludden. 2021.  How Are Sports-Trackers Used by Runners? Running-Related Data, Personal Goals, and Self-Tracking in Running.  _Sensors_ 21, 11 (Jan. 2021), 3687.  https://doi.org/10.3390/s21113687 Number: 11 Publisher: Multidisciplinary Digital Publishing Institute.
  * Luccioni and Viviano (2021) Alexandra Sasha Luccioni and Joseph D. Viviano. 2021.  What’s in the Box? A Preliminary Analysis of Undesirable Content in the Common Crawl Corpus.  http://arxiv.org/abs/2105.02732 arXiv:2105.02732 [cs].
  * Smith et al. (2013) Jason R. Smith, Herve Saint-Amand, Magdalena Plamada, Philipp Koehn, Chris Callison-Burch, and Adam Lopez. 2013.  Dirt cheap web-scale parallel text from the Common Crawl. In _Smith, Jason R; Saint-Amand, Herve; Plamada, Magdalena; Koehn, Philipp; Callison-Burch, Chris; Lopez, Adam (2013). Dirt cheap web-scale parallel text from the Common Crawl. In: 51st Annual Meeting of the Association for Computational Linguistics, Sofia, Bulgaria, August 2013. Association for Computational Linguistics, 1374-1383._ Association for Computational Linguistics, Sofia, Bulgaria, 1374–1383.  https://doi.org/10.5167/uzh-80038
  * Thompson (2022) Alan D Thompson. 2022.  What’s in my AI? A Comprehensive Analysis of Datasets Used to Train GPT-1, GPT-2, GPT-3, GPT-NeoX-20B, Megatron-11B, MT-NLG, and Gopher.  (2022).
  * Turski et al. (2023) Michał Turski, Tomasz Stanisławek, Karol Kaczmarek, Paweł Dyda, and Filip Graliński. 2023.  CCpdf: Building a High Quality Corpus for Visually Rich Documents from Web Crawl Data. In _Document Analysis and Recognition - ICDAR 2023_ , Gernot A. Fink, Rajiv Jain, Koichi Kise, and Richard Zanibbi (Eds.). Springer Nature Switzerland, Cham, 348–365.  https://doi.org/10.1007/978-3-031-41682-8_22
  * Weber et al. (2023) Maurice Weber, Carlo Siebenschuh, Rory M Butler, Anton Alexandrov, Valdemar R Thanner, Georgios Tsolakis, Haris Jabbar, Ian Foster, Bo Li, and Rick Stevens. 2023.  WordScape: a Pipeline to extract multilingual, visually rich Documents with Layout Annotations from Web Crawl Data.  (2023).
  * Zhu et al. (2023) Wenhao Zhu, Hongyi Liu, Qingxiu Dong, Jingjing Xu, Shujian Huang, Lingpeng Kong, Jiajun Chen, and Lei Li. 2023.  Multilingual Machine Translation with Large Language Models: Empirical Results and Analysis.  http://arxiv.org/abs/2304.04675 arXiv:2304.04675 [cs].