# Imperfect Axion Precludes the Domain Wall Problem

Yue Zhang  Department of Physics, Carleton University, Ottawa, Ontario K1S
5B6, Canada

###### Abstract

The QCD axion needs not be an exact pseudoscalar for solving the strong CP
problem. Its imperfectness can play a profound role cosmologically. We propose
effective operators, where the Peccei-Quinn field linearly couples to Standard
Model particles, provide a dynamical solution to the domain wall problem that
prevails in post-inflationary axion models with discrete symmetry. Such
interactions generate a thermal potential that drives the axion field to a
universal value throughout the universe at high temperatures thus preventing
the birth of domain walls when the QCD potential switches on. We discuss
generic conditions for this mechanism to work and several concrete examples.
Combining with existing electric dipole moment and fifth force constraints, a
lower bound on the axion mass is obtained around $10^{-5}$ eV. Our findings
make a strong case for complementary axion searches with both quality
preserving and violating interactions.

The QCD axion Peccei and Quinn (1977a, b); Weinberg (1978); Wilczek (1978)
provides the leading solution to the strong CP problem in the Standard Model
and is subject to a rich program of terrestrial and astrophysical probes
Sikivie (2021); Irastorza and Redondo (2018). Embedding the spontaneous $U(1)$
Peccei-Quinn symmetry breaking theory in the context of expanding universe
brings about both opportunities and problems Marsh (2016); Di Luzio  _et al._
(2020). For ultraviolet complete models with more than one flavor of Peccei-
Quinn changed quarks, a discrete symmetry often remains after the breaking of
Peccei-Quinn symmetry Dine  _et al._ (1981); Zhitnitsky (1980). If it occurs
after inflation, axion domain walls separating different vacuum field values
populate the universe Sikivie (1982). Energy density of domain walls redshifts
slower than those of radiation and matter. If long-lived they would be
disastrous to the otherwise successful standard Big Bang cosmology.

A common solution to the domain problem is to introduce explicit symmetry
breaking effect through the scalar potential of the Peccei-Quinn field Sikivie
(1982); Barr and Seckel (1992); Sikivie (2008); Barr and Kim (2014); Caputo
and Reig (2019). The explicit breaking effect is usually expected to be highly
suppressed for the solution to strong CP problem to remain intact. Taking the
operator $\mathcal{L}=e^{i\delta}\phi^{n}/M_{\rm pl}^{n-4}$ as example, it
contributes an axion potential $V\sim(f_{a}^{n}/M_{\rm
pl}^{n-4})\cos(a/f_{a}+\delta)$ after the spontaneous symmetry breaking. The
phase $\delta$ has no reason to be close to $\bar{\Theta}_{\text{\sc qcd}}$
from the regular QCD axion potential. Neutron electric dipole moment thus sets
a limit $n\gtrsim 200/\ln(\sqrt{2}M_{\rm pl}/f_{a})$, e.g., $n\gtrsim 12$ for
$f_{a}=10^{11}\,$GeV Dine (2022); Banerjee  _et al._ (2022). Given the strong
Planck-scale suppression, such an explicit breaking term has little
cosmological impact at temperatures well above the QCD scale. After domain
walls appear, it can push the walls toward each other which eventually
annihilate away. Products from the domain walls are littered around the
universe. Excessive axion dark matter production and stellar cooling restrict
the axion mass within a rather narrow window Kawasaki  _et al._ (2015);
Ringwald and Saikawa (2016); Beyer and Sarkar (2022); Chang  _et al._ (1999).

In this letter, we aim for a more clean solution to the axion domain wall
problem by precluding their existence in the first place. A key observation is
that domain walls are not born immediately after the spontaneous Peccei-Quinn
symmetry breaking but have to wait until temperature of the universe cools
down to the QCD scale. This separation of scales allows one to imagine the
following scenario. At temperatures well above the QCD scale, explicit Peccei-
Quinn breaking effect generates a transient potential that allows the axion
field to settle down at a universal initial value everywhere. At lower
temperatures, this potential weakens and eventually gives way to the QCD axion
potential. As a result, the whole universe rolls toward today’s vacuum (that
solves the strong CP problem) starting from the universal initial value set
early on. No domain wall is produced in this case even if the low-energy axion
potential still possesses an approximate discrete symmetry. For such a
mechanism to work, the necessary conditions are: 1) axion potential generated
by explicit Peccei-Quinn breaking term is temperature dependent and more
important at higher temperatures; 2) the high temperature potential must not
have any discrete symmetry; 3) there should be enough time for the settling
down of axion field in the first step to complete. Clearly, the $\phi^{n}$
term mentioned above cannot fulfill this role for more than one reason. An
earlier attempt Barr and Kim (2014) resorts to a high-scale confining force to
generate the axion potential above the QCD scale. However, additional model
building acrobatics is needed for switching off the potential at lower
temperatures otherwise the strong CP problem would strike back.

We point out a unique class of explicit Peccei-Quinn breaking interactions
that can provide a dynamical, less harmful at low energy and still predictive
realization of the above idea,

$\mathcal{L}=\phi\mathcal{O}_{\rm SM}+{\rm h.c.}\ ,$ (1)

where $\mathcal{O}_{\rm SM}$ is a gauge-invariant operator made of Standard
Model fields. An obvious choice for $\mathcal{O}_{\rm SM}$ would be a term in
the Standard Model Lagrangian which is automatically neutral under the Peccei-
Quinn symmetry. At high temperatures the Standard Model particle plasma
contributes a potential of the form $\delta V\sim T^{4}\cos(a/f_{a})$, whereas
the zero-temperature counterpart is more suppressed, with $T$ replaced by one
of the known particle mass scales. Crucially, the high-temperature potential
has no discrete symmetry as long as the above interaction is linear in the
Peccei-Quinn field $\phi$.

For concreteness, we present several working examples for realizing the above
idea, where $\mathcal{O}_{\rm SM}$ is one of the regular Yukawa interaction
operators for fermion mass generation. First consider the electron Yukawa
case,

$\mathcal{L}=\frac{e^{i\delta}}{\Lambda}\frac{\sqrt{2}m_{e}}{v}\phi\bar{L}_{e}He_{R}+{\rm
h.c.}\ ,$ (2)

where $v=246\,$GeV, $L_{e}$ is the $SU(2)_{L}$ lepton doublet, and $H$ is the
Higgs doublet. For generality, we include a phase factor $\delta$ that differs
from $\bar{\Theta}_{\text{\sc qcd}}$ by $\gg 10^{-10}$. Hereafter, we focus on
the physics after spontaneous Peccei-Quinn symmetry breaking with the radial
mode of $\phi$ integrated out, such that $\phi=(f_{a}/\sqrt{2})e^{a/f_{a}}$.

Figure 1: Feynman diagrams that generate an axion field potential at high
temperatures before (left) and after (right) electroweak symmetry breaking.
The diamond vertex is given by the operator in Eq. (2) or (10) whereas the dot
vertex is the regular fermion Yukawa coupling. In the right diagram, the
diamond contains the Higgs condensate.

A thermal potential of axion can be generated by the explicit breaking term
Eq. (2), as shown by the diagrams in Figure 1. Because the Standard Model
particles are in the thermal plasma, the potential can be derived using the
imaginary time formalism Quiros (1999). At $\Lambda^{-1}$ order,

$V(a,T)=\left\\{\begin{array}[]{cc}\displaystyle{-\frac{5m_{e}^{2}f_{a}T^{4}}{288\sqrt{2}v^{2}\Lambda}}\cos\left(\frac{a}{f_{a}}+\delta\right),&T\gg
v\ ;\vspace{2mm}\\\
\displaystyle{-\frac{m_{e}^{2}f_{a}T^{2}}{3\sqrt{2}\Lambda}}\cos\left(\frac{a}{f_{a}}+\delta\right),&\Lambda_{\rm
QCD}\ll T<v\ ,\\\ \end{array}\right.$ (3)

where all loop particles are treated massless. This serves as the only axion
potential at temperatures well above the QCD scale. The axion field’s
evolution follows the equation

$\ddot{a}+3H\dot{a}+\partial_{a}V(a,T)=0\ .$ (4)

The high-temperature potential drives the axion field value toward the minimum
with $a/f_{a}+\delta=0$, whereas the Hubble parameter exerts damping effect.
We focus on the zero mode evolution in each early patch of the universe that
is causally disconnected from others and could start with all possible initial
axion field values. Effects due to inhomogeneities are discussed near the end
of the paper. To solve the domain wall problem using the mechanism described
above, the axion field must roll sufficiently close to the bottom of the above
potential regardless of its initial condition. The temperature (time) window
for this to occur is from the spontaneous Peccei-Quinn symmetry breaking scale
to the QCD scale.

Figure 2: Parameter space of explicit Peccei-Quinn symmetry breaking effect
from the operators Eq. (2) (left) and (10) (right, strange quark case). The
axion domain wall problem is avoided below the orange shaded region using the
mechanism suggested in this work. The solid (dashed) orange curve corresponds
to a low energy axion potential with $Z_{N}$ discrete symmetry with
$N=2\,(6)$. Lower bounds on the cutoff scale $\Lambda$ come from the tests of
equivalence principle (green) and inverse square law of gravity (blue), red
giant cooling (purple), and the neutron electric dipole moment measurements
(gray) which exclude the shaded regions with various colors. The electric
dipole moment constraint becomes increasingly stronger for heavier fermions.
For the strange quark case, we find the lightest possible axion for solving
the axion domain wall problem and staying consistent with existing
experimental constraints, with mass $\sim 10^{-5}$ eV.

We introduce $\zeta=\ln R$ where $R$ is the scale factor of the expanding
universe. This allows Eq. (4) to be written as

$\frac{d^{2}}{d\zeta^{2}}\frac{a}{f_{a}}+\frac{d}{d\zeta}\frac{a}{f_{a}}+\frac{cM_{\rm
Pl}^{2}}{\Lambda f_{a}}\sin\left(\frac{a}{f_{a}}+\delta\right)=0\ ,$ (5)

where $M_{\rm Pl}$ is the Planck scale. The dimensionless parameter $c$ is a
constant $1.8\times 10^{-16}$ for $T\gg v$, whereas for $\Lambda_{\rm QCD}\ll
T<v$, $c$ is temperature dependent and given by $\simeq 2.1\times
10^{-10}\,{\rm GeV}^{2}e^{2\zeta}/f_{a}^{2}$. We assume that the Peccei-Quinn
symmetry is spontaneously broken at $T=f_{a}$ and set the corresponding
$\zeta$ to zero. Qualitatively, for the axion field to quickly roll to the
minimum of the above potential it requires

$\frac{cM_{\rm Pl}^{2}}{\Lambda f_{a}}\gtrsim 1\ ,$ (6)

leading to an upper bound on $\Lambda$ for given $f_{a}$. In practice, we take
a number of random initial values for $a/f_{a}$ ranging between 0 and $2\pi$
and numerically evolve above equation from $\zeta=0$ to
$\ln(f_{a}/\Lambda_{\text{\sc QCD}})$. The bound on $\Lambda$ is then derived
by requiring for all initial conditions,

$\left|\frac{a}{f_{a}}+\delta-\bar{\Theta}_{\text{\sc
qcd}}\right|\lesssim\frac{\pi}{N}\ ,$ (7)

at $\zeta=\ln(f_{a}/\Lambda_{\rm QCD})$. On the right-hand side, the factor
$N$ comes about assuming the low energy axion potential from QCD has a $Z_{N}$
discrete symmetry ($N>1$). The value $\pi/N$ represents the typical distance
of phase $\delta$ from the closest maximum of the QCD axion potential. It
gives a good estimate except that the high temperature minimum is located in
the vicinity of one of the zero temperature maxima. We disregard such fine-
tuned arrangements of parameters because they go against the spirit of
naturalness for solving the strong CP problem. Numerically, our result stays
robust against order one variations of $\pi/N$. Eq. (7) guarantees that all
early patches of the universe end up in the same vacuum thus avoiding the
production of any walls. The resulting upper bound on $\Lambda$ as a function
of $f_{a}$ is shown by the orange curves in the left panel of Figure 2, for
$N=2$ and 6 cases.

The explicit Peccei-Quinn symmetry breaking term Eq. (2) can affect a number
of low-energy probes. At zero temperature, it shifts the electron mass and in
turn sources a Coleman-Weinberg potential for the axion,

$V(a,T=0)=-\frac{|M_{e}(a)|^{4}}{64\pi^{2}}\left(\ln\frac{|M_{e}(a)|^{2}}{\mu^{2}}-\frac{3}{2}\right)\
,$ (8)

which adds to the regular QCD axion potential and can lower the axion quality.
Here
$M_{e}(a)=m_{e}\left[1+{f_{a}}e^{i(a/f_{a}+\delta)}/(\sqrt{2}\Lambda)\right]$
and $\mu$ is the renormalization scale. An order one difference between the
$\delta$ and $\bar{\Theta}_{\text{\sc qcd}}$ phases will result in a nonzero
contribution to the neutron electric dipole moment and requires
$\Lambda>5.3\times 10^{10}\,{\rm GeV}\left({10^{-5}\,\rm eV}/{m_{a}}\right)$
Zhang (2023).

In the presence of the explicit Peccei-Quinn breaking, the axion is no longer
a pure pseudoscalar boson. It induces a scalar coupling to electron, allowing
the axion to mediate a “fifth force” between macroscopic objects. Parameters
that control the violation of equivalence principle and the inverse square law
of gravity are Zhang (2023)

$\tilde{\alpha}=\frac{m_{e}^{2}\sin^{2}(\delta-\bar{\Theta}_{\text{\sc
qcd}})}{8\pi
Gu^{2}\Lambda^{2}},\quad\alpha=\frac{Z_{1}Z_{2}}{A_{1}A_{2}}\frac{m_{e}^{2}\sin^{2}(\delta-\bar{\Theta}_{\text{\sc
qcd}})}{8\pi Gu^{2}\Lambda^{2}}\ ,$ (9)

where $G$ is the gravitational constant, $u=931.5\,$MeV, $Z_{1,2}$ and
$A_{1,2}$ denote the atomic charge and weight of the gravitating objects.
Existing torsion balance experimental constraints on $\tilde{\alpha}$ and
$\alpha$ are presented in Adelberger  _et al._ (2009). In addition, heavier
axions with scalar coupling to electrons are constrained by stellar cooling
Hardy and Lasenby (2017). In Figure 2, experimental lower limits on $\Lambda$
are depicted, along with the upper bound for addressing the domain wall
problem. In this case, fifth force constraint dominates and we obtain a lower
bound on the axion mass $m_{a}\gtrsim 10^{-4}\,{\rm eV}$.

Next, we move on to the coupling of Peccei-Quinn field directly to quarks as
the source of explicit breaking,

$\mathcal{L}=\frac{e^{i\delta}}{\Lambda}\frac{\sqrt{2}m_{q}}{v}\phi\bar{Q}Hq_{R}+{\rm
h.c.}\ ,$ (10)

For up-type quark, $H$ is replaced by $i\sigma_{2}H^{*}$. At high temperatures
above the QCD scale, thermal loops similar to Figure 1 can generate an axion
potential that allows it to roll toward the universal value $a/f_{a}+\delta=0$
throughout the universe. It is useful to note that parametrically the thermal
potential is enhanced compared to Eq. (3) by a factor of
$(m_{q}/m_{e})^{2}N_{c}$, where $N_{c}=3$ is the color factor. The larger
coupling makes it relatively easier for the axion to settle down to the above
minimum to avoid domain walls.

For the low-energy axion potential, naively one may expected it to be
generated at tree level once the quark condensate turns on. However, note the
tree-level CP-violating effect of (10) is to modify the phase of quark masses,
which are part of $\bar{\Theta}_{\text{\sc QCD}}$. After minimizing the axion
potential, they make no net contribution to the neutron electric dipole moment
Zhang (2023). Instead, the leading axion potential is again of the Coleman-
Weinberg form, i.e., Eq. (8) but with the mass of electron replaced by the
corresponding quark mass along with a color factor. It is worth noting that
such a potential is proportional to fourth power of the fermion mass. The high
power of mass dependence implies that the electric dipole moment constraint
will stand out further for heavier quarks.

Through the $\phi$-quark operator, axion also mediates a fifth force by
coupling to the nuclei. In this case, the parameters $\tilde{a}$ and $a$ are
given by

$\tilde{\alpha}=\alpha\simeq\frac{f_{q}^{2}m_{N}^{2}\sin^{2}(\delta-\bar{\Theta}_{\text{\sc
qcd}})}{8\pi Gu^{2}\Lambda^{2}}\ ,$ (11)

where $m_{N}\simeq 938\,$MeV is the nucleon mass, $N$ represents proton or
neutron, and $f_{q}$ is defined through the hadronic matrix element, $\langle
N|m_{q}\bar{q}q|N\rangle=f_{q}m_{N}\bar{u}_{N}u_{N}$. We work in the isospin
conserving limit and use values of $f_{q}$ given in Cline  _et al._ (2013).
Compared to Eq. (9), $\tilde{a}$ and $a$ are enhanced by the square of nucleon
to electron mass ratio, but as a constraint on $\Lambda$ the mass effect is
linear.

Figure 2 (right panel) shows the results of Eq. (10) where the Peccei-Quinn
field couples to a the strange quark Yukawa term. The color schemes are the
same as before, with lower limit on $\Lambda$ set by various experimental
results and upper bound for solving the axion domain wall problem. Stellar
cooling only constrains $\Lambda$ up to $10^{10}$ eV in these cases Hardy and
Lasenby (2017) and lies well below the range of $\Lambda$ shown in the figure.
Compared to the electron case (left panel), we find that the viable window for
$\Lambda$ is pushed to higher values because the axion coupling is enhanced by
the larger fermion Yukawa coupling. The neutron electric dipole moment
constraint becomes significantly more important than equivalence principle and
inverse square law tests of gravity, due to the $m_{s}^{4}$ versus $m_{s}^{2}$
dependences discussed above. In the strange quark coupling case, we derive the
widest axion mass window with a lower bound

$m_{a}\gtrsim 10^{-5}\,{\rm eV}\ .$ (12)

We explore the couplings of Peccei-Quinn field to other quark and lepton
flavors as well, by turning on one coupling in each case. We find that fifth
force constraints are more important for light fermions ($u,d,e$) whereas EDM
dominates for heavier fermions. The strange quark case is roughly where the
transition occurs. We explore other dimension-5 operators where the Peccei-
Quinn field couples to the Higgs and gauge boson fields. We also consider $\mu
e^{i\delta}\phi H^{\dagger}H$ which is a renormalizable operator. For the
axion solution to the strong CP problem to work, the $\mu$ parameter needs to
be very small. To prevent domain walls $\mu$ is also bounded from below.
Combining the two requirements, we identify a window for $\mu$ as a function
of axion mass. Similar to $\Lambda$ for dimension-5 operators, the window
closes at a lower bound on the axion mass. All the results are listed in table
1. More calculation details are provided in the supplemental material. Among
all operators considered, the lowest possible axion mass is found for Peccei-
Quinn field coupling to the strange or muon Yukawa operator, or the gluonic
operator.

$\mathcal{O}_{\rm SM}$ | Yukawa
---|---
$u$ | $d$ | $s$ | $c$ | $b$ | $t$
$m_{a}^{\text{min}}$ (eV) | $10^{-4.3}$ | $10^{-4.4}$ | $10^{-5.2}$ | $10^{-4.0}$ | $10^{-3.5}$ | $10^{-2.0}$
$\Lambda$ (GeV) | $10^{17.5}$ | $10^{18.1}$ | $10^{20.0}$ | $10^{23.4}$ | $10^{25.0}$ | $10^{29.9}$
$\mathcal{O}_{\rm SM}$ | Yukawa | other operators
$e$ | $\mu$ | $\tau$ | $G_{\mu\nu}G^{\mu\nu}$ | $(H^{\dagger}H)^{2}$ | $H^{\dagger}H$
$m_{a}^{\text{min}}$ (eV) | $10^{-4.3}$ | $10^{-4.9}$ | $10^{-3.7}$ | $10^{-5.3}$ | $10^{-2.8}$ | $10^{-2.8}$
$\Lambda\,(\mu)$ (GeV) | $10^{15.6}$ | $10^{19.9}$ | $10^{23.6}$ | $10^{25.3}$ | $10^{33.1}$ | $10^{-28.6}$
Table 1: The lowest possible axion mass and corresponding cutoff scale
$\Lambda$ from precluding axion domain walls using the mechanism suggested in
this work. Results are presented for effective operators where the Peccei-
Quinn field linearly couples to various Standard Model gauge invariant
operators $\mathcal{O}_{\rm SM}$. In the last column, the $\phi H^{\dagger}H$
operator has dimension 3 and $\Lambda$ is replaced by a dimensionful coupling
$\mu$.

By precluding the birth of axion domain walls, explicit Peccei-Quinn breaking
effects introduced through the directly coupling with Standard Model particles
(Eq. (10) and (2)) can open up a much wider window of axion mass compared to
the Peccei-Quinn field self-interaction. In the latter case, domain walls
still occur and viable axion masses are constrained within a very narrow range
around tens of meV by the dark matter relic abundance and supernova cooling
Kawasaki  _et al._ (2015); Ringwald and Saikawa (2016); Beyer and Sarkar
(2022).

The viable axion mass range found above is consistent with the parameter space
where axion comprises all cold dark matter in the universe via the
misalignment mechanism Preskill  _et al._ (1983); Abbott and Sikivie (1983);
Dine and Fischler (1983). The initial axion field misalignment value is now
set as a by-product of the suggested domain wall problem solution,
$a/f_{a}=\delta-\bar{\Theta}_{\text{\sc qcd}}$. Moreover, the lower part of
the axion mass range can be probed by the major axion haloscope and helioscope
experiments Du  _et al._ (2018); Braine  _et al._ (2020); Rybka (2014);
Brubaker  _et al._ (2017); Zhong  _et al._ (2018); Alesini  _et al._ (2019);
Brun  _et al._ (2019); Lee  _et al._ (2020); Armengaud  _et al._ (2019). In
complementarity, the parameter space along the $\Lambda$ direction governs the
axion quality and will continue to be tested by the future electric dipole
moment Abel  _et al._ (2020); Ahmed  _et al._ (2019a, b), fifth force
experiments Adelberger  _et al._ (2009), and monopole-dipole interactions
Geraci  _et al._ (2018); Crescini  _et al._ (2017).

To summarize, whether the QCD axion is an exact pseudoscalar has been asked as
a phenomenological question in various contexts Moody and Wilczek (1984);
Kamionkowski and March-Russell (1992); Holman  _et al._ (1992); Raffelt
(2012); Mantry  _et al._ (2014); Bertolini  _et al._ (2021); Okawa  _et al._
(2022); Dekens  _et al._ (2022). This work explores the cosmological
significance of a new class of explicit PQ-symmetry breaking effective
operators of the form $\phi\mathcal{O}_{\rm SM}$ which extends previous
solutions to the axion domain wall problem. Such an operator generates a
thermal potential for the axion at temperatures well above the QCD scale, thus
allowing the axion to roll and settle down at a universal field value in all
patches of the early universe. After the universe cools further, the QCD axion
potential takes over but no domain walls is created even if there remains a
discrete symmetry, thanks to the above early focusing effect. The proposed
mechanism makes prediction in the axion quality and can be tested by a number
of ways, including the neutron electric dipole moment and fifth force
experiments. Existing limits set a lower bound on axion mass around $10^{-5}$
eV for this mechanism to work. Our solution to domain problem opens up a broad
window for cosmological viable axion models and strongly motivates a
complementary search of both the axion quality preserving and violating
interactions with known particles. The above mass bound does not apply if the
axion model is free from discrete symmetries Georgi and Wise (1982); Lazarides
and Shafi (1982); Barr  _et al._ (1982); Alonso-Álvarez  _et al._ (2023), or
the Peccei-Quinn symmetry is already broken before or during inflation Seckel
and Turner (1985); Linde (1985); Turner and Wilczek (1991); Dvali and
Senjanovic (1995); Redi and Tesi (2022).

We end by commenting on possible topological defects at very high
temperatures, around the spontaneous Peccei-Quinn symmetry breaking scale.
Axion string network could arise because the axion field initial values are
random among various early patches of the universe. The explicit Peccei-Quinn
breaking effect discussed in this work cannot only generate the high-
temperature potential for setting a common axion field value but also render
the axion strings unstable. Because the early potential features no discrete
global symmetry, any domain walls if ever formed along with the strings are
unzipped and destroyed right away Vilenkin and Everett (1982). All these would
happen at temperatures shortly after the spontaneous symmetry breaking. The
universe afterwards is filled with a common axion field value prior to the
QCD-scale temperatures, i.e., the possible existence of these topological
effects at early times does not affect the success of the suggested mechanism.
Their potential roles in the production of axion relic density and the extra
clustering of axion dark matter in structure formation are interesting open
questions but beyond the scope of this work.

###### Acknowledgements.

I acknowledge useful discussions with Gonzalo Alonso-Álvarez, Luca Di Luzio,
and Miha Nemevšek at various stages of this work. This work is supported by
the Arthur B. McDonald Canadian Astroparticle Physics Research Institute and
the Natural Sciences and Engineering Research Council of Canada.

## References

  * Peccei and Quinn (1977a) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn (1977b) R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977b).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Sikivie (2021) P. Sikivie, Rev. Mod. Phys. 93, 015004 (2021), arXiv:2003.02206 [hep-ph] .
  * Irastorza and Redondo (2018) I. G. Irastorza and J. Redondo, Prog. Part. Nucl. Phys. 102, 89 (2018), arXiv:1801.08127 [hep-ph] .
  * Marsh (2016) D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Di Luzio  _et al._ (2020) L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Dine  _et al._ (1981) M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
  * Zhitnitsky (1980) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 529 (1980).
  * Sikivie (1982) P. Sikivie, Phys. Rev. Lett. 48, 1156 (1982).
  * Barr and Seckel (1992) S. M. Barr and D. Seckel, Phys. Rev. D 46, 539 (1992).
  * Sikivie (2008) P. Sikivie, Lect. Notes Phys. 741, 19 (2008), arXiv:astro-ph/0610440 .
  * Barr and Kim (2014) S. M. Barr and J. E. Kim, Phys. Rev. Lett. 113, 241301 (2014), arXiv:1407.4311 [hep-ph] .
  * Caputo and Reig (2019) A. Caputo and M. Reig, Phys. Rev. D 100, 063530 (2019), arXiv:1905.13116 [hep-ph] .
  * Dine (2022) M. Dine, (2022), arXiv:2207.01068 [hep-ph] .
  * Banerjee  _et al._ (2022) A. Banerjee, J. Eby, and G. Perez, (2022), arXiv:2210.05690 [hep-ph] .
  * Kawasaki  _et al._ (2015) M. Kawasaki, K. Saikawa, and T. Sekiguchi, Phys. Rev. D 91, 065014 (2015), arXiv:1412.0789 [hep-ph] .
  * Ringwald and Saikawa (2016) A. Ringwald and K. Saikawa, Phys. Rev. D 93, 085031 (2016), [Addendum: Phys.Rev.D 94, 049908 (2016)], arXiv:1512.06436 [hep-ph] .
  * Beyer and Sarkar (2022) K. A. Beyer and S. Sarkar, (2022), arXiv:2211.14635 [hep-ph] .
  * Chang  _et al._ (1999) S. Chang, C. Hagmann, and P. Sikivie, Phys. Rev. D 59, 023505 (1999), arXiv:hep-ph/9807374 .
  * Quiros (1999) M. Quiros, in  _ICTP Summer School in High-Energy Physics and Cosmology_ (1999) pp. 187–259, arXiv:hep-ph/9901312 .
  * Zhang (2023) Y. Zhang, Phys. Rev. D 107, 055025 (2023), arXiv:2209.09429 [hep-ph] .
  * Adelberger  _et al._ (2009) E. G. Adelberger, J. H. Gundlach, B. R. Heckel, S. Hoedl, and S. Schlamminger, Prog. Part. Nucl. Phys. 62, 102 (2009).
  * Hardy and Lasenby (2017) E. Hardy and R. Lasenby, JHEP 02, 033 (2017), arXiv:1611.05852 [hep-ph] .
  * Cline  _et al._ (2013) J. M. Cline, K. Kainulainen, P. Scott, and C. Weniger, Phys. Rev. D 88, 055025 (2013), [Erratum: Phys.Rev.D 92, 039906 (2015)], arXiv:1306.4710 [hep-ph] .
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie (1983) L. F. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler (1983) M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Du  _et al._ (2018) N. Du _et al._ (ADMX), Phys. Rev. Lett. 120, 151301 (2018), arXiv:1804.05750 [hep-ex] .
  * Braine  _et al._ (2020) T. Braine _et al._ (ADMX), Phys. Rev. Lett. 124, 101303 (2020), arXiv:1910.08638 [hep-ex] .
  * Rybka (2014) G. Rybka (ADMX), Phys. Dark Univ. 4, 14 (2014).
  * Brubaker  _et al._ (2017) B. M. Brubaker _et al._ , Phys. Rev. Lett. 118, 061302 (2017), arXiv:1610.02580 [astro-ph.CO] .
  * Zhong  _et al._ (2018) L. Zhong _et al._ (HAYSTAC), Phys. Rev. D 97, 092001 (2018), arXiv:1803.03690 [hep-ex] .
  * Alesini  _et al._ (2019) D. Alesini _et al._ , Phys. Rev. D 99, 101101 (2019), arXiv:1903.06547 [physics.ins-det] .
  * Brun  _et al._ (2019) P. Brun _et al._ (MADMAX), Eur. Phys. J. C 79, 186 (2019), arXiv:1901.07401 [physics.ins-det] .
  * Lee  _et al._ (2020) S. Lee, S. Ahn, J. Choi, B. R. Ko, and Y. K. Semertzidis, Phys. Rev. Lett. 124, 101802 (2020), arXiv:2001.05102 [hep-ex] .
  * Armengaud  _et al._ (2019) E. Armengaud _et al._ (IAXO), JCAP 06, 047 (2019), arXiv:1904.09155 [hep-ph] .
  * Abel  _et al._ (2020) C. Abel _et al._ , Phys. Rev. Lett. 124, 081803 (2020), arXiv:2001.11966 [hep-ex] .
  * Ahmed  _et al._ (2019a) S. Ahmed _et al._ (TUCAN), Phys. Rev. C 99, 025503 (2019a), arXiv:1809.04071 [physics.ins-det] .
  * Ahmed  _et al._ (2019b) M. W. Ahmed _et al._ (nEDM), JINST 14, P11017 (2019b), arXiv:1908.09937 [physics.ins-det] .
  * Geraci  _et al._ (2018) A. A. Geraci _et al._ (ARIADNE), Springer Proc. Phys. 211, 151 (2018), arXiv:1710.05413 [astro-ph.IM] .
  * Crescini  _et al._ (2017) N. Crescini, C. Braggio, G. Carugno, P. Falferi, A. Ortolan, and G. Ruoso, Nucl. Instrum. Meth. A 842, 109 (2017), arXiv:1606.04751 [physics.ins-det] .
  * Moody and Wilczek (1984) J. E. Moody and F. Wilczek, Phys. Rev. D 30, 130 (1984).
  * Kamionkowski and March-Russell (1992) M. Kamionkowski and J. March-Russell, Phys. Lett. B 282, 137 (1992), arXiv:hep-th/9202003 .
  * Holman  _et al._ (1992) R. Holman, S. D. H. Hsu, T. W. Kephart, E. W. Kolb, R. Watkins, and L. M. Widrow, Phys. Lett. B 282, 132 (1992), arXiv:hep-ph/9203206 .
  * Raffelt (2012) G. Raffelt, Phys. Rev. D 86, 015001 (2012), arXiv:1205.1776 [hep-ph] .
  * Mantry  _et al._ (2014) S. Mantry, M. Pitschmann, and M. J. Ramsey-Musolf, Phys. Rev. D 90, 054016 (2014), arXiv:1401.7339 [hep-ph] .
  * Bertolini  _et al._ (2021) S. Bertolini, L. Di Luzio, and F. Nesti, Phys. Rev. Lett. 126, 081801 (2021), arXiv:2006.12508 [hep-ph] .
  * Okawa  _et al._ (2022) S. Okawa, M. Pospelov, and A. Ritz, Phys. Rev. D 105, 075003 (2022), arXiv:2111.08040 [hep-ph] .
  * Dekens  _et al._ (2022) W. Dekens, J. de Vries, and S. Shain, JHEP 07, 014 (2022), arXiv:2203.11230 [hep-ph] .
  * Georgi and Wise (1982) H. Georgi and M. B. Wise, Phys. Lett. B 116, 123 (1982).
  * Lazarides and Shafi (1982) G. Lazarides and Q. Shafi, Phys. Lett. B 115, 21 (1982).
  * Barr  _et al._ (1982) S. M. Barr, X. C. Gao, and D. Reiss, Phys. Rev. D 26, 2176 (1982).
  * Alonso-Álvarez  _et al._ (2023) G. Alonso-Álvarez, J. M. Cline, and T. Xiao, (2023), arXiv:2305.00018 [hep-ph] .
  * Seckel and Turner (1985) D. Seckel and M. S. Turner, Phys. Rev. D 32, 3178 (1985).
  * Linde (1985) A. D. Linde, Phys. Lett. B 158, 375 (1985).
  * Turner and Wilczek (1991) M. S. Turner and F. Wilczek, Phys. Rev. Lett. 66, 5 (1991).
  * Dvali and Senjanovic (1995) G. R. Dvali and G. Senjanovic, Phys. Rev. Lett. 74, 5178 (1995), arXiv:hep-ph/9501387 .
  * Redi and Tesi (2022) M. Redi and A. Tesi, (2022), arXiv:2211.06421 [hep-ph] .
  * Vilenkin and Everett (1982) A. Vilenkin and A. E. Everett, Phys. Rev. Lett. 48, 1867 (1982).


# On the use of galaxies as clocks and the universal expansion

Anders Ahlström Kjerrgren1 and Edvard Mörtsell2
1KTH Royal Institute of Technology, SE 106 91 Stockholm, Sweden
2Oskar Klein Centre, Department of Physics, Stockholm University, SE 106 91
Stockholm, Sweden E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We set out to rederive the 8 Hubble parameter values obtained from estimated
relative galaxy ages by Simon et al. [Physical Review D, 71, 123001 (2005)].
We find that to obtain the level of precision claimed in $H(z)$,
unrealistically small galaxy age uncertainties have to be assumed. Also, some
parameter values will be correlated. In our analysis we find that the
uncertainties in the Hubble parameter values are significantly larger when 8
independent $H(z)$ are obtained using Monte Carlo sampling. Smaller
uncertainties can be obtained using Gaussian processes, but at the cost of
strongly correlated results. We do not obtain any useful constraints on the
Hubble parameter from the galaxy data employed.

###### keywords:

(cosmology:) cosmological parameters – (cosmology:) distance scale – methods:
data analysis

††pubyear: 2022††pagerange: On the use of galaxies as clocks and the universal
expansion–A.1

##  1 Introduction

In 2005 the paper titled “Constraints on the redshift dependence of the dark
energy potential” was published in Physical Review D (Simon et al., 2005). In
section IV.A of this paper 8 independent values of the Hubble parameter are
obtained at a large range of redshifts to a relatively high degree of
precision. These values are being used to test other aspects of cosmology, see
for instance Gómez-Valent et al. (2020); Odintsov et al. (2020); Chen et al.
(2017); Melia & Yennapureddy (2018); Gómez-Valent & Amendola (2018). Some
papers use the results, but cite compilations which include the results from
Simon et al. (2005) instead; examples include Renzi & Silvestri (2020) who
cited Stern et al. (2010), and Bonilla et al. (2021) who cited Moresco et al.
(2016).

The purpose of this work is to investigate the reliability of the resulting
values of the Hubble parameter. The practical implementation of the method
used in Simon et al. (2005) is not trivial. First, it requires determining the
ages of galaxies accurately and precisely to obtain the dependence on the
redshift, $t(z)$. The second step is to infer the derivative
$\operatorname{d\\!}t/\operatorname{d\\!}z$ of this relation, a task more
difficult than inferring the function values in the sense that it, in general,
leads to larger uncertainties. Assuming that the first step is done
flawlessly, we show that the results in Simon et al. (2005) assumes
unrealistically small galaxy age uncertainties due to the difficulty in making
precise inferences of derivatives. 8 values of the Hubble parameter are
obtained, and therefore 8 derivatives are computed, with uncertainties between
10 and 20$\%$ using 32 galaxies, indicating that age uncertainties much
smaller than the reported 12% have been assumed.111This work is based on the
master’s thesis Ahlström Kjerrgren (2021), written by Anders Ahlström
Kjerrgren and supervised by Edvard Mörtsell.

##  2 The Hubble parameter from galaxy data

It was proposed by Jimenez & Loeb (2002) that if passively evolving galaxies
at different redshifts are considered the Hubble parameter can be obtained by
computing the age difference of the galaxies. The equation

$H(z)=-\frac{1}{1+z}\derivative{z}{t}$ (1)

relates the relevant quantities. The derivative
$\operatorname{d\\!}z/\operatorname{d\\!}t$ is taken with respect to the
cosmic time, however if the galaxies considered are born at the same cosmic
time the derivative eliminates the constant offset and galaxy ages can thus be
used as a proxy for the cosmic time. The derivative is then further
approximated as the ratio $\Delta z/\Delta t$.

The method requires the identification of the same type of galaxies which at
lower redshifts are older versions of those at higher redshifts (Weinberg et
al., 2013). Identifying such galaxies yields a relation between galaxy age and
redshift. This relation is visible in the data used by Simon et al. (2005),
shown in the left panel of Fig. 1.

By binning the data and using equation (1) 8 determinations of the Hubble
parameter are obtained by Simon et al. (2005), shown in the right panel of
Fig. 1, with 10 or 20% 1-$\sigma$ uncertainties.

Figure 1: Data and results in Simon et al. (2005). Error bars correspond to
1-$\sigma$ uncertainties. Left: The galaxy age and redshift data. All age
uncertainties are 12%, as given to Samushia et al. (2010) by the authors in a
private communication. Right: The Hubble parameter values.

##  3 Attempt at reproducing the results

We first attempt to reproduce the results as close as possible. For clarity we
focus on a few select $H(z)$: the ones at the lowest, second lowest, and
highest redshift. See appendix A for a full attempt at reproducing the
results.

With the four galaxies at the lowest redshifts one obtains the effective
redshift ${z\simeq 0.17}$ corresponding to the Hubble parameter value at the
lowest redshift. Determining $H(z)$ requires the determination of the
derivative $\operatorname{d\\!}z/\operatorname{d\\!}t$. By fitting a straight
line to the data using $\chi^{2}$ minimisation we identify the slope as the
reciprocal derivative $\operatorname{d\\!}t/\operatorname{d\\!}z$. Combining
this with equation (1) yields the same central value as Simon et al. (2005),
namely $H(z\simeq 0.17)\simeq 83$ km/s/Mpc. However, using 12% uncertainties
in the galaxy ages results in $>100\%$ uncertainty in the Hubble parameter. It
is also obtained when computing a weighted average of the first and second
galaxy, and the third and fourth galaxy, and then calculating $H(z)$ with
$\Delta z/\Delta t$ using these two weighted averages and regular uncertainty
propagation. To obtain 10% uncertainty in $H(z)$, as Simon et al. (2005), the
age uncertainties must be lowered to $\simeq 1.2\%$, i.e. a factor of ten.

The Hubble parameter at the second lowest redshift must be computed using the
third, fourth, and fifth galaxy (counting from lowest redshift), as they yield
the effective redshift $z\simeq 0.29$ and $H(z)\simeq 70$ km/s/Mpc. For this
value the uncertainty is in the order of 80%. To obtain 20% uncertainty, as
Simon et al. (2005), the age uncertainties must be lowered to $\simeq 3\%$.

The observant reader will notice that the third and fourth galaxies are used
for both the first and second Hubble parameter values, meaning that the
results are not independent. However, this correlation is not quantified by
Simon et al. (2005).

For the $H(z)$ value at the highest redshift we find no combination of
galaxies which yields the same effective redshift as the one presented by
Simon et al. (2005).

###  3.1 The galaxy age uncertainties

The above analysis demonstrates that a significant reduction in galaxy age
uncertainties is necessary to obtain the same precision in $H(z)$ as Simon et
al. (2005). It is claimed that some of the systematic errors present in the
absolute age determinations cancel when computing the relative ages (Verde,
2021). However, as we argue below, we find that the assumed size of the
cancellation underestimates the relative age uncertainties.

First, it is unclear to what degree error bars on the galaxy ages in Fig. 1
(their fig. 1) represent statistical or systematic uncertainties. The ages are
determined using the SPEED model (Jimenez et al., 2004) from which one obtains
a likelihood surface over metallicity and age. Marginalising over metallicity,
the age and its _statistical_ uncertainty is obtained. Furthermore, systematic
errors are discussed after the analysis in their fig. 1 is complete,
indicating that systematic effects are ignored, not cancelled.

Since galaxy age determinations are quite model dependent, we find it unlikely
that such a large portion of the uncertainties cancel. 20 of the 32 galaxies
in Simon et al. (2005) are taken from the Gemini Deep Deep Survey (GDDS)
(McCarthy et al., 2004) in which ages are estimated. These galaxy ages are
reanalysed by Simon et al. (2005) using the SPEED model. It is claimed that
the ages do not change by much, but Fig. 2 shows that this is not the case,
and that using a different model can greatly change the estimated ages. The
relative ages may also differ drastically, as seen in Fig. 3, which shows the
percentage change in $\Delta t$ when going from the GDDS data set to the data
set in Simon et al. (2005). We compute $\Delta t$ for galaxies adjacent in
redshift where matching redshifts exist, marked with dashed blue lines in Fig.
2, as both data sets presumably refer to the same galaxies at these redshifts.

Figure 2: Comparison between the data in Simon et al. (2005) and the data
presented in McCarthy et al. (2004) for the GDDS data set. Dashed blue lines
indicate where the data sets have points with the same redshift. The two
galaxies 53W069 and 53W091 are also included. Figure 3: The percentage changes
in $\Delta t$ when going from the GDDS data set in McCarthy et al. (2004) to
the data set in Simon et al. (2005). Here, $\Delta t$ is only computed for
galaxies adjacent in redshift where matching redshifts exist, marked with
dashed blue lines in Fig. 2. At these redshifts, both data sets presumably
refer to the same galaxies. Galaxy pairs are counted from lowest redshift,
i.e. pair 1 is computed with the galaxies at the lowest and second lowest
redshifts, pair 2 with the galaxies at the second and third lowest redshifts,
and so on.

The fact that obtaining the precise $H(z)$ in Simon et al. (2005) requires
unrealistically small age uncertainties is corroborated when considering that
other papers which also compute the Hubble parameter using galaxy ages have
one thing in common: many more galaxies are used for each $H(z)$, without
obtaining significantly smaller uncertainties. To name a few examples: 1
$H(z)$ is obtained with 7% precision using $\sim$14000 galaxies by Moresco et
al. (2011), 8 $H(z)$ are obtained with 5–17% precision using $\sim$11000
galaxies by Moresco et al. (2012), 2 $H(z)$ are obtained with 21 & 27%
precision using more than 30 galaxies by Moresco (2015), 5 (or 1) $H(z)$ are
obtained with 11–16% (or 6%) precision using $\sim$130000 galaxies by Moresco
et al. (2016), and 1 $H(z)$ is obtained with 75% precision using 16 galaxies
by Ratsimbazafy et al. (2017).

Furthermore, we simulate an idealised case where galaxy ages are placed on the
$\Lambda$CDM prediction for $t(z)$ and galaxy redshifts are placed equidistant
between $z=0$ and $z=0.25$.222Assuming that the data set contains galaxies
with $z\in[0,2]$ and one wants 8 $H(z)$ this is the bin width to use. We find
that $\sim 70$ ($\sim 280$) galaxies with 12% age uncertainties are necessary
in this bin to obtain $H(z)$ with 20% (10%) precision. For bins at higher
redshifts even more galaxies are needed.

##  4 Reanalysing the data

We have shown that the resulting Hubble parameter values in Simon et al.
(2005) are too precise given the employed data. In the following sections we
use the same data and reanalyse it with two different methods, Monte Carlo
sampling and Gaussian processes, to show specifically how one can obtain
$H(z)$ from $t(z)$ and what uncertainties should be expected using this data
set.

###  4.1 Monte Carlo sampling

We place the galaxy data in bins with width $\Delta z$. In each bin we sample
the ages from normal distributions and perform a least squares linear fit
where the slope is identified as $\operatorname{d\\!}t/\operatorname{d\\!}z$.
With equation (1) we obtain a distribution for $H(z)$ from which we identify
the Hubble parameter as the median and the uncertainty limits by the central
region containing 68.27% of the data.

In each bin a constant derivative is assumed, which introduces an error. This
error is quantified using the $\Lambda$CDM prediction for $H(z)$ as the
fractional error $H(z_{\mathrm{max}})/H(z_{\mathrm{min}})-1$, with
$z_{\mathrm{max}}$ ($z_{\mathrm{min}}$) being the maximum (minimum) galaxy
redshift in the bin. Cosmological parameters for the $\Lambda$CDM model are
taken from Zyla et al. (2020).

We present the results in Fig. 4. When 8 $H(z)$ are obtained the uncertainties
are larger than $\sim 90\%$. To obtain $\sim 20\%$ uncertainties much larger
bin widths must be used, such that only 2 $H(z)$ are obtained. However, then
the binning error is too large for the inference to be useful.

Figure 4: Relative uncertainty in $H(z)$ and the fractional binning errors vs
bin width, obtained through Monte Carlo sampling. Relative uncertainty defined
as half the width of the 68.27$\%$ area divided by the median of the samples.

To visualise the precision and accuracy of the $H(z)$ obtained we show the
actual values in Fig. 5 for a few select bin widths. For comparison we include
the $\Lambda$CDM prediction and vary $\Omega_{\Lambda}$ by $\pm 10\%$. The
resulting $H(z)$ are not precise enough to determine $\Omega_{\Lambda}$ with
this precision and are unlikely to be useful when constraining other
cosmological parameters. Furthermore, when 8 $H(z)$ are obtained we can not
confidently conclude that the universe is expanding, as the vertical error
bars extend to negative values.

Figure 5: $H(z)$ obtained using Monte Carlo sampling for a few bin widths,
with vertical error bars showing the central 68.27% limits. The dashed black
line is the $\Lambda$CDM prediction, and the grey area is obtained by changing
$\Omega_{\Lambda}$ by $\pm 10\%$ and using
$\Omega_{\mathrm{M}}=1-\Omega_{\Lambda}$. In no case is the result good enough
to distinguish between different values of $\Omega_{\Lambda}$ in this range.
The dashed vertical lines mark the bins, whereas the horizontal error bars
indicate the minimum and maximum redshift of galaxies in the bins.
_Top left:_ $\Delta z=0.2$. _Top right:_ $\Delta z=0.4$. _Bottom left:_
$\Delta z=0.6$. _Bottom right:_ $\Delta z=1.2$.

###  4.2 Gaussian processes

The method of Gaussian processes (GPs) generates a continuous inference where
every inferred point uses all of the data to some extent. The trade-off is
that the inferences at different points are correlated. This method has been
used in cosmology for instance to infer properties of dark energy (Holsclaw et
al., 2011; Seikel et al., 2012; Wang & Meng, 2017; Zhang & Li, 2018; Seikel &
Clarkson, 2013), the Hubble constant333Note that many papers use the resulting
Hubble parameter values from Simon et al. (2005) in their inference of
$H_{0}$. (Gómez-Valent & Amendola, 2018; Liao et al., 2019; Busti et al.,
2014; Bengaly et al., 2020), and other applications (Haridasu et al., 2018;
Shafieloo et al., 2012; Bengaly, 2020; Belgacem et al., 2020). Thus, we
consider this method to infer the Hubble parameter.

A Gaussian process is completely determined by its mean function and
covariance function, both of which have to be chosen. We follow convention and
choose the zero mean function. As for the covariance function we use five
different ones: the squared exponential (SE), and four versions from the so
called Matérn class functions.

For an extensive overview of Gaussian processes we refer to Rasmussen &
Williams (2006). However, in this work we follow the method outlined in
section 2 of Seikel et al. (2012), more specifically section 2.2 and 2.3 where
the derivative of a function is reconstructed and then used to construct a
function of said derivative by sampling from a predicted distribution of the
derivative. We reconstruct $\operatorname{d\\!}t/\operatorname{d\\!}z$ and use
equation (1) to obtain $H(z)$.

The results are shown in Fig. 6, where the results from Simon et al. (2005)
and the $\Lambda$CDM prediction are included for comparison. The uncertainties
are very large for larger redshifts and relatively large for the smallest
redshifts, but here the GPs divert from the $\Lambda$CDM prediction. For the
centre redshifts the GP uncertainties are of comparable magnitude to the ones
in Simon et al. (2005). However, considering Fig. 7 we see that these smaller
uncertainties come at the cost of strong correlations for inferences at nearby
redshifts and quite strong anti-correlations for separated inferences.

Figure 6: $H(z)$ obtained by sampling from the predicted distribution of
$\operatorname{d\\!}t/\operatorname{d\\!}z$. Connected lines are the median of
the distribution at each redshift $z$. The central 68.27% of the samplings are
within the dashed lines. Also included are the results from Simon et al.
(2005) and the $\Lambda$CDM prediction for $H(z)$. Inset shows a larger range
in $H(z)$ at the higher redshifts. Figure 7: Correlation matrices for the
derivative $\operatorname{d\\!}t/\operatorname{d\\!}z$, similar to those for
$H(z)$. SE refers to the squared exponential, and $\nu$ refers to the
corresponding Matérn class function.

##  5 Summary and conclusion

We have reanalysed the work by Simon et al. (2005) in which galaxy age and
redshift data is used to obtain the Hubble parameter by computing relative
galaxy ages. We have shown that the results are too precise given the
precision in the galaxy ages. To obtain the same precision, age uncertainties
have to decrease from 12% to 1–3%, which is unlikely even if some of the
systematic errors in the absolute age determinations cancel when computing the
relative ages. Fig. 3 indicates that systematic effects between different
relative age determinations are in the order of 90%. Simulating an idealised
scenario we found that many more galaxies are needed, where the number is more
in line with the number of galaxies used in other papers. Furthermore, some
galaxies in Simon et al. (2005) are used to compute multiple $H(z)$ meaning
that the results are not independent.

Comparing the data in Simon et al. (2005) with the cited data we found a quite
large difference. Nevertheless, we chose to ignore this discrepancy and
reanalysed the data presented in Simon et al. (2005). This was done with two
methods: Monte Carlo sampling, and Gaussian processes. The former method
showed that obtaining 8 $H(z)$ yields uncertainties larger than $\sim 90\%$.
By increasing the bin width 20% uncertainties could be obtained, but then only
2 $H(z)$ were obtained and the error introduced by the binning was too large
for the inference to be useful. The latter method yielded smaller
uncertainties for much of the redshift interval, but this came at the cost of
strong correlations and anti-correlations.

The analysis done in this work raises doubts about the accuracy of the
resulting Hubble parameter values and corresponding uncertainties by Simon et
al. (2005) and for this reason, we do not recommend their use for constraining
cosmological models.

## Data Availability

The data in this article is derived from a source in the public domain, where
the data being analysed is available in the public domain as figures made by
Simon et al. (2005) or in a table composed by Samushia et al. (2010).

## References

  * Ahlström Kjerrgren (2021) Ahlström Kjerrgren A., 2021, Master’s thesis, KTH Royal Institute of Technology
  * Belgacem et al. (2020) Belgacem E., Foffa S., Maggiore M., Yang T., 2020, Physical Review D, 101, 063505
  * Bengaly (2020) Bengaly C. A., 2020, Monthly Notices of the Royal Astronomical Society: Letters, 499, L6
  * Bengaly et al. (2020) Bengaly C. A., Clarkson C., Maartens R., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 053
  * Bonilla et al. (2021) Bonilla A., Kumar S., Nunes R. C., 2021, The European Physical Journal C, 81, 1
  * Busti et al. (2014) Busti V. C., Clarkson C., Seikel M., 2014, Monthly Notices of the Royal Astronomical Society: Letters, 441, L11
  * Chen et al. (2017) Chen Y., Kumar S., Ratra B., 2017, The Astrophysical Journal, 835, 86
  * Gómez-Valent & Amendola (2018) Gómez-Valent A., Amendola L., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 051
  * Gómez-Valent et al. (2020) Gómez-Valent A., Pettorino V., Amendola L., 2020, Physical Review D, 101, 123513
  * Haridasu et al. (2018) Haridasu B. S., Luković V. V., Moresco M., Vittorio N., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 015
  * Holsclaw et al. (2011) Holsclaw T., Alam U., Sansó B., Lee H., Heitmann K., Habib S., Higdon D., 2011, Physical Review D, 84, 083501
  * Jimenez & Loeb (2002) Jimenez R., Loeb A., 2002, The Astrophysical Journal, 573, 37
  * Jimenez et al. (2004) Jimenez R., MacDonald J., Dunlop J. S., Padoan P., Peacock J. A., 2004, Monthly Notices of the Royal Astronomical Society, 349, 240
  * Liao et al. (2019) Liao K., Shafieloo A., Keeley R. E., Linder E. V., 2019, The Astrophysical Journal Letters, 886, L23
  * McCarthy et al. (2004) McCarthy P. J., et al., 2004, The Astrophysical Journal Letters, 614, L9
  * Melia & Yennapureddy (2018) Melia F., Yennapureddy M. K., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 034
  * Moresco (2015) Moresco M., 2015, Monthly Notices of the Royal Astronomical Society: Letters, 450, L16
  * Moresco et al. (2011) Moresco M., Jimenez R., Cimatti A., Pozzetti L., 2011, Journal of Cosmology and Astroparticle Physics, 2011, 045
  * Moresco et al. (2012) Moresco M., et al., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 006
  * Moresco et al. (2016) Moresco M., et al., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 014
  * Odintsov et al. (2020) Odintsov S. D., Gómez D. S.-C., Sharov G. S., 2020, Physical Review D, 101, 044010
  * Rasmussen & Williams (2006) Rasmussen C. E., Williams C. K. I., 2006, Gaussian Processes for Machine Learning. The MIT Press
  * Ratsimbazafy et al. (2017) Ratsimbazafy A., Loubser S., Crawford S., Cress C., Bassett B., Nichol R., Väisänen P., 2017, Monthly Notices of the Royal Astronomical Society, 467, 3239
  * Renzi & Silvestri (2020) Renzi F., Silvestri A., 2020, arXiv preprint arXiv:2011.10559
  * Samushia et al. (2010) Samushia L., Dev A., Jain D., Ratra B., 2010, Physics Letters B, 693, 509
  * Seikel & Clarkson (2013) Seikel M., Clarkson C., 2013, arXiv preprint arXiv:1311.6678
  * Seikel et al. (2012) Seikel M., Clarkson C., Smith M., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 036
  * Shafieloo et al. (2012) Shafieloo A., Kim A. G., Linder E. V., 2012, Physical Review D, 85, 123530
  * Simon et al. (2005) Simon J., Verde L., Jimenez R., 2005, Physical Review D, 71, 123001
  * Stern et al. (2010) Stern D., Jimenez R., Verde L., Kamionkowski M., Stanford S. A., 2010, Journal of Cosmology and Astroparticle Physics, 2010, 008
  * Verde (2021) Verde L., 2021, Private Communication
  * Wang & Meng (2017) Wang D., Meng X.-H., 2017, Physical Review D, 95, 023508
  * Weinberg et al. (2013) Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., Riess A. G., Rozo E., 2013, Physics reports, 530, 87
  * Zhang & Li (2018) Zhang M.-J., Li H., 2018, The European Physical Journal C, 78, 460
  * Zyla et al. (2020) Zyla P. A., et al., 2020, Progress of Theoretical and Experimental Physics, 2020, 083C01

##  Appendix A Full reproduction attempt

We begin by attempting to reproduce the binning by following the instructions
in Simon et al. (2005). We separate the description into three steps, which
are visualised in Fig. 8 (step zero is the unprocessed redshift data for the
galaxies in their sample), where the colour coding is described in the figure
caption.

The first step is to simply group together galaxies that are within $\Delta
z=0.03$. After this step it is said that most groups contain more than one
galaxy, but we obtain nineteen groups with only eight groups containing two or
more galaxies. Nevertheless, in the next step galaxies which are $>2\sigma$
away from oldest galaxy in a group are supposed to be discarded. We find no
such galaxies and therefore do not include this step in Fig. 8.

The second step is to compute the relative ages for groups which are separated
by $\Delta z\in[0.1,0.15]$. This effectively merges the groups into the final
bins, and is what we visualise in Fig. 8.

For each bin obtained we compute the effective redshift as the mean of the
redshifts in the bin. This is step three in Fig. 8, where we also include the
resulting effective redshifts in Simon et al. (2005) for comparison. We obtain
a quite different set of redshifts and for this reason do not compute the
Hubble parameter using the obtained bins. The difference is already large
enough to conclude that we find it difficult to reproduce the results by
following the instructions.

Figure 8: Visualisation of our attempt of reproducing the binning in Simon et
al. (2005). _Step 0:_ Unprocessed redshift data for the 32 galaxies. _Step 1:_
Galaxies within $\Delta z=0.03$ have been grouped together. Different colours
correspond to different groups. _Step 2:_ Groups from the previous step which
are separated by at least 0.1 but no more than 0.15 have been placed in the
same bin. Data points with the same shape and colour correspond to the same
bin. Note that some data points are in two different bins. _Step 3:_ The
effective redshifts in Simon et al. (2005) (black and marked with stars), and
the effective redshifts computed for every bin from the previous step
(coloured as the bins in the previous step and marked with crosses).

###  A.1 Matching the effective redshifts

As we could not reproduce the binning by following the instructions in Simon
et al. (2005) we will in this section disregard it and instead attempt to find
the bins which yield the same effective redshifts as those presented in Simon
et al. (2005).

To this end we first construct all subsets of the 32 data points which have a
size ranging from two to ten data points, and for which the difference between
the minimum and maximum redshifts is smaller than $0.2$. This limit is chosen
since it is slightly more generous than than the limit of $\Delta z=0.15$
combined with the grouping criteria of $\Delta z=0.03$, discussed above.444The
reason we do not set it to $0.18$ is because we do not want to be too
restrictive, as there is something which we do not understand in binning
process in Simon et al. (2005), as demonstrated in the previous section. After
this filtering, the effective redshifts are computed to two decimals for all
remaining subsets, defined as the arithmetic mean of the redshifts in each
subset. As a last step we compare the obtained effective redshifts with the
ones in Simon et al. (2005) and choose the bins with matching effective
redshift.555In the case of multiple matches, we choose the bin which produces
the closest effective redshift, using all decimals, among the bins with the
most amount of data points. The resulting bins are shown in the middle panel
of Fig. 9. We are able to match all redshifts, except for the last where we
instead choose the bin producing the closest effective redshift.

Figure 9: _Top:_ Unprocessed redshift data. _Middle:_ The bins which produce
(almost) the same effective redshifts as in Simon et al. (2005). Data
corresponding to the same bin have the same colour. Note that some data points
are used twice, i.e are in two bins. _Bottom:_ The effective redshifts of the
above bins, with the same colour coding. Only the last redshift was not
managed to be matched.

Using these bins we compute the derivative
$\operatorname{d\\!}t/\operatorname{d\\!}z$ using a linear $\chi^{2}$
minimisation fit in each bin. The result is shown in Fig. 10, which
demonstrates uncertainties in the order of 100%. The bins used here should be
credible to be close to the ones used in Simon et al. (2005), but the obtained
uncertainties are significantly larger.

Figure 10: The derivative $\operatorname{d\\!}t/\operatorname{d\\!}z$ for the
bins which match the effective redshifts in Simon et al. (2005), obtained
using a linear $\chi^{2}$ minimisation fit.