11institutetext: Center for Cosmology and Particle Physics, Department of
Physics, New York University, New York, NY 10003, USA22institutetext:
Department of Physics, University of California, Berkeley, California 94720,
USA33institutetext: Theoretical Physics Group, Lawrence Berkeley National
Laboratory, Berkeley, California 94720, USA44institutetext: Department of
Physics, University of Chicago, Chicago, IL 60637, USA55institutetext: Enrico
Fermi Institute, University of Chicago, Chicago, IL 60637, USA66institutetext:
Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL
60637, USA77institutetext: Kavli Institute for the Physics and Mathematics of
the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The
University of Tokyo, Kashiwa, Chiba 277-8583, Japan88institutetext:
Theoretical Physics Department, CERN, Geneva, Switzerland

#  A Heavy QCD Axion and the Mirror World

David I. Dunsky 2,3 Lawrence J. Hall 4,5,6,7,8 Keisuke Harigaya

###### Abstract

We study the mirror world with dark matter arising from the thermal freeze-out
of the lightest, stable mirror particle – the mirror electron. The dark matter
abundance is achieved for mirror electrons of mass 225 GeV, fixing the mirror
electroweak scale near $10^{8}$ GeV. This highly predictive scenario is
realized by an axion that acts as a portal between the two sectors through its
coupling to the QCD and mirror QCD sectors. The axion is more massive than the
standard QCD axion due to additional contributions from mirror strong
dynamics. Still, the strong CP problem is solved by this ‘heavy’ axion due to
the alignment of the QCD and mirror QCD potentials. Mirror entropy is
transferred into the Standard Model sector via the axion portal, which
alleviates overproduction of dark radiation from mirror glueball decays. This
mirror scenario has a variety of signals: (1) primordial gravitational waves
from the first-order mirror QCD phase transition occurring at a temperature
near 35 GeV, (2) effects on large-scale structure from dark matter self-
interactions from mirror QED, (3) dark radiation affecting the cosmic
microwave background, and (4) the rare kaon decay,
$K^{+}\rightarrow(\pi^{+}+\rm{axion})$. The first two signals do not depend on
any fundamental free parameters of the theory while the latter two depend on a
single free parameter, the axion decay constant.

##  1 Introduction

A unified theory of nature may fragment into a low-energy effective theory
having multiple sectors, each with its own gauge symmetry. If the fermions
fall into sets carrying gauge charges of only a single sector, then the
lightest fermion(s) in each sector are stable. In the Standard Model (SM)
sector these are the lightest neutrino, the electron, and the proton. The
cosmological dark matter may be the lightest fermion of some other sector, the
dark sector. While this framework is plausible, there are many such schemes
each with a variety of free parameters, so that it is hard to construct
realistic theories that can be tested.

A unique and highly predictive scheme is the mirror world, where a $Z_{2}$
symmetry leads to an exact copy of the sector describing the directly observed
particles and interactions. Indeed, before it was known that our sector was
described by gauge symmetry, this mirror world was introduced as a way to
preserve parity in nature Lee:1956qn ; Kobzarev:1966qya  as the $Z_{2}$
symmetry can also send $\bar{r}\rightarrow-\bar{r}$ and flip fermion
chirality. Over successive decades, interest in the mirror world, and the
possibility of mirror dark matter, has increased Okun:2006eb ; Foot:2014mia .
Furthermore, with the discovery of the possibility of kinetic mixing between
photons of different $U(1)$ gauge groups Holdom:1985ag , it was proposed that
kinetic mixing could probe dark sectors Goldberg:1986nk . However, mirror dark
matter, with the mirror spectrum identical to the SM spectrum, faces several
challenges. Measurements of the cosmic microwave background require the mirror
sector to be at a lower temperature than the SM sector, both from mirror
neutrino contributions to dark radiation and from mirror baryon acoustic
oscillations Cyr-Racine:2012tfp ; Bansal:2022qbi . In addition, the $Z_{2}$
symmetry must be broken by interactions beyond the SM, or by initial
conditions, so that mirror baryons have a larger cosmological abundance than
SM baryons. Furthermore, the mirror matter halos should not dissipate and
collapse into disks like SM matter, and must avoid the constraints on self-
interactions from observations of the Bullet cluster.

There is an alternative route to the mirror world. Beyond kinetic mixing,
there is one other interaction of dimension 4 that couples the two sectors,
$|H|^{2}|H^{\prime}|^{2}$, where $H$ and $H^{\prime}$ are the SM and mirror
Higgs doublets. This term is not excluded by any symmetry, and it changes
everything! If it is absent, as discussed above, the two Higgs vevs are equal,
$v^{\prime}=v$, so that the mirror spectrum is identical to the SM spectrum.
Two cases emerge if this additional Higgs portal operator is significant. In
the first, the vacuum is determined at tree-level and has two phases
Foot:2000tp . The symmetric phase has the two vevs equal, and is excluded by
Higgs physics at the LHC. In the asymmetric phase one vev vanishes, so that
mirror matter is extremely light and it is unclear how it yields dark matter.
A realistic possibility for mirror dark matter emerges from adding a soft
$Z_{2}$ breaking Higgs interaction, giving a hierarchy of vevs with
$v^{\prime}>v$ Barbieri:2005ri . The second case, more minimal and
constraining, has the vacuum determined by radiative corrections, via the
Higgs Parity mechanism Hall:2018let . This gives one vev much larger than the
other, $v^{\prime}\gg v$, with the ratio determined by the measured values of
the Higgs boson mass, the top quark mass, and the QCD coupling constant. With
current data the central value of $v^{\prime}$ is $10^{12}$ GeV, with a
$3\sigma$ lower bound of $10^{9}$ GeV.

In a previous paper Dunsky:2019upk , we explored this mirror world from Higgs
Parity and found successful cosmologies with thermally produced mirror
electron, $e^{\prime}$, dark matter. With a high reheat temperature after
inflation, the freeze-out abundance of $e^{\prime}$ leads to the observed dark
matter if $v^{\prime}=10^{8}$ GeV, well below the predicted range of
$v^{\prime}$. Furthermore, after the mirror QCD phase transition the lightest
mirror glueball decays to mirror photons leading to too much dark radiation.
Remarkably, a realistic cosmology results because, with $v^{\prime}$ of order
$(10^{9}-10^{10})$ GeV, the $e^{\prime}$ and dark radiation can be
sufficiently diluted by entropy created by the decays of mirror neutrinos.
This cosmology requires the top quark mass to be about $2\sigma$ high.
Alternatively, if only the SM sector is produced after inflation with a low
reheat temperature, $e^{\prime}$ dark matter can arise from freeze-in, via
either the Higgs or kinetic mixing portal, for a wide range of $v^{\prime}$.

In this paper, we study the above mirror world with the addition of a Peccei-
Quinn (PQ) symmetry Peccei:1977hh ; Peccei:1977ur  giving a KSVZ axion that is
even under the $Z_{2}$ symmetry (for the case of a Weinberg-Wilczek axion see
Berezhiani:2000gh ). This leads to a successful cosmology with $e^{\prime}$
dark matter from freeze-out without the need for dilution. The new colored
states carrying the PQ symmetry, the PQ quarks, have couplings with the Higgs
boson that modify the Higgs Parity mechanism to allow $v^{\prime}=10^{8}$ GeV.
Furthermore, the mirror glueball decays to axions that deposit their energy in
the SM sector, avoiding over-production of dark radiation.

While parity may solve the strong CP problem, it does not do so in the mirror
world; it simply relates the non-zero strong CP parameters of the two sectors.
A $Z_{2}$-even axion solves the strong CP problem and is much heavier than the
conventional axion Rubakov:1997vp ; Fukuda:2015ana ; Hook:2019qoh , since the
dominant contribution to its mass arises from mirror QCD which gets strong
near 30 GeV. If such a heavy axion decays after neutrino decoupling it heats
up the photon bath relative to the neutrinos, suppressing the neutrino
contribution to the radiation energy. To avoid this the axion scale must be
low, $f_{a}<10^{5}$ GeV, and the axion mass high, of order 10 MeV. Such a
low-$f_{a}$, high-mass axion greatly ameliorates the quality problem of the PQ
symmetry, which can be more easily understood as an accidental consequence of
other symmetries.

The mirror world has also been motivated as a setting for the Twin Higgs
mechanism Chacko:2005pe , which improves the naturalness of the weak scale. In
this case $v^{\prime}$ is around the TeV scale, far below that predicted by
the exact mirror world, so that soft breaking of $Z_{2}$ in the Higgs
potential is required Barbieri:2005ri . Such soft breaking could arise from
spontaneous breaking in some other sector. We also consider this possibility
for obtaining $v^{\prime}=10^{8}$ GeV, needed for $e^{\prime}$ dark matter
from freeze-out, without the need of any Yukawa coupling between the PQ quarks
and the Higgs boson. Refs. Barbieri:2016zxn ; Fukuda:2017ywn ;
Barbieri:2017opf  consider various dark matter candidates for the theories
with soft and hard $Z_{2}$ breaking.

In addition to the SM parameters, the mirror world studied in this paper has
three parameters relevant for cosmological and particle physics signals,
$v^{\prime},\epsilon$ and $f_{a}$. Without dilution from mirror neutrino
decay, the $e^{\prime}$ dark matter abundance fixes $v^{\prime}=10^{8}$ GeV.
Thus, the effects of self-interactions on the dark matter halo ellipticity
from mirror electromagnetism involve no free parameters. Similarly, gravity
waves from the mirror QCD phase transition involve no free parameters of the
underlying particle theory, since the phase transition temperature can be
computed. Furthermore, the gravity wave signal is not diluted by entropy
generated from neutrino decay, as in Dunsky:2019upk . On the other hand, the
predicted rate for $K\rightarrow\pi a$ depends on $f_{a}$. The amount of dark
radiation is sensitive to axion physics and also depends on $f_{a}$; in
addition it is sensitive to the electromagnetic anomaly of the PQ symmetry.
Direct detection of $e^{\prime}$ dark matter can occur via kinetic mixing and
depends on $\epsilon$. If the SM is embedded in a unified theory at scale
$v_{G}$, then $\epsilon$ depends sensitively on $v_{G}$, so that the direct
detection signal is correlated with the proton decay rate Dunsky:2019upk .

We present the theory in section 2, including portal operators and the PQ
sector, and discuss the spontaneous breaking of the $Z_{2}$ via Higgs vevs. In
section 3 we elaborate on key aspects of particle phenomenology, including the
spectrum of mirror fermions, properties of the axion and the lightest mirror
glueball, and the generation of neutrino and mirror neutrino masses. The
cosmological history is presented in section 4: thermal decoupling between the
two sectors is studied, as well as the parameter space that yields the
observed dark matter via mirror electron freeze-out. In section 5, we present
signals of this scheme, arising from dark radiation, rare $K$ decays, dark
matter self-interactions, and gravity waves. In section 6 we explore the
possibility of increasing $v^{\prime}$ so that the freeze-out abundance of
mirror electrons is too high, and subsequently reduced by entropy production
from mirror neutrino decay. In the appendix, we present a model where the
axion quality problem is solved such that the PQ symmetry is an accidental
consequence of anomaly-free discrete symmetries.

##  2 The Mirror Theory with Exact Parity and an Axion

Motivated by mirror dark matter and the axion quality problem, we study a
theory containing the Standard Model (SM), its $Z_{2}$ symmetric mirror (SM′),
and a single axion field coupled to both sectors. In this section, we discuss
the Lagrangrian of the two sectors and determine the spectrum of particles.

###  2.1 Lagrangian

The $Z_{2}$ symmetry maps the SM gauge group and particles into their $Z_{2}$
mirrors

$\displaystyle SU(3)\times SU(2)\times U(1)\;\;\leftrightarrow\;\;$
$\displaystyle SU(3)^{\prime}\times SU(2)^{\prime}\times U(1)^{\prime}$
$\displaystyle q,\bar{u},\bar{d},\ell,\bar{e}\;\;\leftrightarrow\;\;$
$\displaystyle
q^{\prime},\bar{u}^{\prime},\bar{d}^{\prime},\ell^{\prime},\bar{e}$
$\displaystyle H\;\;\leftrightarrow\;\;$ $\displaystyle H^{\prime}$
$\displaystyle F^{\mu\nu},W^{\mu\nu},G^{\mu\nu}\;\;\leftrightarrow\;\;$
$\displaystyle F^{\prime\mu\nu},W^{\prime\mu\nu},G^{\prime\mu\nu}$
$\displaystyle\Psi,\bar{\Psi}\;\;\leftrightarrow\;\;$
$\displaystyle\Psi^{\prime},\bar{\Psi}^{\prime}$ $\displaystyle
P\;\;\leftrightarrow\;\;$ $\displaystyle P\,,$ (1)

where primes indicate mirror fields, and matter is described by 2-component,
left-handed, Weyl fields. $\Psi$ and its mirror, $\Psi^{\prime}$, are heavy
quarks charged under the PQ symmetry while
$P=\frac{1}{\sqrt{2}}(s+f_{a})e^{ia/f_{a}}$ is the PQ breaking field that
contains the axion, $a$, as its angular mode and the saxion, $s$, as its
radial mode.

The Lagrangian of the theory is given by

$\displaystyle\mathcal{L}=\mathcal{L}_{\rm SM}+\mathcal{L}_{\rm
SM^{\prime}}+\mathcal{L}_{\rm Portal}+\mathcal{L}_{\rm PQ},$ (2)

where $\mathcal{L}_{\rm SM}$ is the SM Lagrangrian up to dimension-5
containing Yukawa interactions, the canonical field kinetic energies, and the
Weinberg neutrino operator $\ell\ell HH$. Similarly, $\mathcal{L}_{\rm
SM^{\prime}}$ is the associated mirror Lagrangian which is related to
$\mathcal{L}_{\rm SM}$ by the mapping given in (2.1). The portal,
$\mathcal{L}_{\rm Portal}$, contains the $Z_{2}$ symmetric operators that mix
the SM and mirror sectors

$\displaystyle\mathcal{L_{\rm
Portal}}=\frac{\epsilon}{2}B_{\mu\nu}B^{\prime\mu\nu}-\lambda^{\prime}|H|^{2}|H^{\prime}|^{2}+\left(\frac{\xi_{ij}}{M_{D}}\ell^{\prime
i}H\;\ell^{j}H^{\prime}+\rm{h.c.}\right),$ (3)

which we call the kinetic mixing portal, Higgs portal, and neutrino portal,
respectively.

Last, the Lagrangian for the PQ breaking field and the PQ quarks is

$\displaystyle\mathcal{L}_{\rm
PQ}=\xi(\bar{\Psi}\Psi+\bar{\Psi}^{\prime}\Psi^{\prime})P$
$\displaystyle+\text{h.c.}+\lambda(|P|^{2}-f_{a}^{2}/2)^{2}+|{D}^{\mu}P|^{2}+i\Psi^{\dagger}\not\\!{\partial}\Psi+\text{h.c.}$
(4)

In Appendix A, we present a model where an exact, extra $Z_{2n+1}$ symmetry
ensures the quality of the PQ symmetry.

In the remainder of this work, we consider the effective theory below the
scale $\langle P\rangle=f_{a}/\sqrt{2}$ which, after integrating out $\Psi$,
generates the axion interactions

$\displaystyle\mathcal{L}_{a}=\frac{1}{32\pi^{2}}\frac{a}{f_{a}}\left(g_{3}^{2}\,G_{\mu\nu}\tilde{G}^{\mu\nu}+g_{3}^{\prime
2}\,G^{\prime}_{\mu\nu}\tilde{G^{\prime}}^{\mu\nu}\right)+\frac{1}{32\pi^{2}}\frac{E}{N}\frac{a}{f_{a}}\left(e^{2}\,F_{\mu\nu}\tilde{F}^{\mu\nu}+e^{\prime
2}\,F^{\prime}_{\mu\nu}\tilde{F^{\prime}}^{\mu\nu}\right)\,.$ (5)

Note that above the mirror symmetry breaking scale $\langle
H^{\prime}\rangle=v^{\prime}\gg v$, the $Z_{2}$ symmetry enforces
$g_{3}=g_{3}^{\prime}$ and $e=e^{\prime}$. Below this scale, the
renormalization group running of the two sectors differs so that
$g_{3}^{\prime}$ and $e^{\prime}$ can diverge from their SM counterparts.

###  2.2 Breaking of Mirror and SM Electroweak Symmetries

If the $Z_{2}$ symmetry of (2.1) is exact, the potential for $H$ and
$H^{\prime}$ is

$\displaystyle
V_{H,H^{\prime}}\;=\;V_{\rm{SM}}(H)+V_{\rm{SM^{\prime}}}(H^{\prime})+\lambda^{\prime}|H|^{2}|H^{\prime}|^{2}.$
(6)

The vacuum, determined at loop level by the Higgs Parity mechanism
Hall:2018let , has one vev much larger than the other. Defining $H^{\prime}$
to have the larger vev, $v^{\prime}$, the low-energy effective potential of
$H$ at tree-level is

$\displaystyle V_{\rm L.E.}=-2\lambda v^{\prime
2}\left(1-\frac{\lambda^{\prime}}{2\lambda}\right)|H|^{2}+\lambda\left(1-\frac{{\lambda^{\prime}}^{2}}{4\lambda^{2}}\right)|H|^{4}\,.$
(7)

Identifying the quadratic term with the SM Higgs mass squared,
$|\lambda^{\prime}-2\lambda|\ll 1$ is required to give $v\ll v^{\prime}$. Eq.
(7), then implies the quartic term of the Higgs potential is nearly zero at
the scale $v^{\prime}$ as it is proportional to $\lambda^{\prime}-2\lambda$
Hall:2018let .

In the SM, the renormalization of the Higgs quartic is dominated by quantum
corrections from the top quark which causes $\lambda$ to decrease from low to
high energies and eventually become zero. From a perspective of running the
Higgs quartic from low to high energies, we can thus identify $v^{\prime}$
with the energy scale $\mu$ at which the Higgs quartic vanishes
$\lambda(\mu=v^{\prime})\simeq 0$. In the SM and in a pure mirror model, where
the only particle content is the SM and its $Z_{2}$ mirror, the scale at which
$\lambda=0$ occurs is at $\mu\approx 10^{12}$ GeV for $m_{\rm top}$ and
$\alpha_{s}$ at their central experimental values Buttazzo:2013uya ;
Workman:2022ynf ; Dunsky:2019upk . Furthermore, $\mu$ is above
$10^{9}(10^{10})$ GeV at $3\sigma(2\sigma)$.

However, as discussed in Sec. 4.2, $v^{\prime}>10^{8}$ GeV is problematic from
a cosmological perspective because freeze-out of the stable mirror electron
overproduces dark matter. One possibility is to dilute the mirror electrons by
large entropy generation, for example from mirror neutrino decay
Dunsky:2019upk . Another possibility is to include particle content beyond the
SM with interactions that makes $\lambda$ run faster so that it vanishes at
the scale $\mu\approx 10^{8}$ GeV. In the mirror model, with an axion having
$f_{a}<v^{\prime}$ as considered in this paper, the PQ quarks $\Psi$ of Eq.
(4) are present in the effective theory below $v^{\prime}$. Furthermore, they
have Yukawa interactions with $H$ if their gauge charges are the same as one
of the SM quark species, $q,\bar{u}$ or $\bar{d}$. For example, choosing the
PQ quark to be an SU(2) singlet with hypercharge $1/3$, and calling it $D$,
allows the Yukawa interaction

$\displaystyle\mathcal{L}_{D}=y_{D}(q\bar{D}H+q^{\prime}\bar{D^{\prime}}H^{\prime})\,.$
(8)

Below $v^{\prime}$, the first operator of Eq. (8) generates a quantum
correction to the beta function of $\lambda\propto-y_{D}^{4}$. Reducing the
scale at which $\lambda=0$ from $10^{12}$ GeV to $10^{8}$ GeV can easily be
accomplished for $y_{D}\sim\mathcal{O}(1)$. Note that $\bar{D}$ is defined as
the linear combination of anti-down type quarks that couples to $P$ via the
coupling $\xi$ of (4). Dark radiation bounds as discussed in Sec. 5.1
constrain $f_{a}$ so that the values of $f_{a}$ of interest to us are of order
30 TeV. Consequently, with $\xi\sim\mathcal{O}(1)$ the mass of $D$ is also of
order 30 TeV.

On the other hand, the mass of $D^{\prime}$ is complicated by the large vev of
$H^{\prime}$ which strongly mixes $D^{\prime}$ with the down-type component of
$q^{\prime}$, $d^{\prime}$, which may be the mirror down, strange, or bottom
quark. The associated mass matrix is

$\displaystyle
M_{d^{\prime}D^{\prime}}=\begin{pmatrix}d^{\prime}&&D^{\prime}\end{pmatrix}\begin{pmatrix}y_{d}v^{\prime}&&y_{D}v^{\prime}\\\
0&&\xi f_{a}\end{pmatrix}\begin{pmatrix}\bar{d^{\prime}}\\\
\bar{D^{\prime}}\end{pmatrix}+\text{h.c.}$ (9)

Eq. (9) possesses a heavy eigenvalue of mass $m_{D,\rm heavy}\sim
y_{D}v^{\prime}\sim v^{\prime}$ and a light eigenvalue of mass $m_{D,\rm
light}\sim\xi f_{a}y_{d}/y_{D}\sim y_{d}f_{a}$. As long as $y_{d}$ is the
bottom quark Yukawa coupling, then $m_{D,\rm light}$ is greater than the
mirror QCD scale (see Eq. (11)) and the $SU(3)^{\prime}$ sector remains an
effective zero-flavor Yang-Mills theory. In this scenario, the
$SU(3)^{\prime}$ phase transition is first order and a strong gravitational
wave signal from bubble nucleation is generated as discussed in Sec. 5.4.

Another option for obtaining $v^{\prime}=10^{8}$ GeV, needed for $e^{\prime}$
dark matter from freeze-out, is to add a soft $Z_{2}$-breaking term to the
Higgs potential of (7): $\delta^{2}(|H|^{2}-|H^{\prime}|^{2})$. In this case,
the PQ quarks need not couple to the Higgs; indeed, they may be neutral under
the electroweak gauge group. This is important as the dark radiation signal,
discussed in Sec. 5.1, is very sensitive to the electromagnetic anomaly of the
PQ symmetry.

##  3 Particle Spectra and Lifetimes

Figure 1: Left panel: spectrum of mirror quarks and leptons for
$v^{\prime}=8\times 10^{7}$ GeV, which yields the observed dark matter
abundance in $e^{\prime}$ from freeze-out. The mass of the lightest mirror
glueball, $S^{\prime}$, is also shown, as is the relevant range of axion
masses. Right panel: important temperatures in the thermal history of the
universe. The mirror electroweak scale, $v^{\prime}$, is also shown.

###  3.1 Mirror Fermions

Below the mirror symmetry breaking scale, $v^{\prime}$, charged mirror
fermions acquire a mass $m^{\prime}=yv^{\prime}$ where $y$, their Yukawa
coupling with $H^{\prime}$, is identical to the Yukawa coupling of their SM
counterparts, up to renormalization differences below $v^{\prime}$.
Consequently, charged mirror fermions, such as the mirror electron,
$e^{\prime}$, are heavier than their SM counterparts by approximately the
ratio of electroweak vacuum expectation values, $v^{\prime}/v\gg 1$. The
spectrum of mirror particles, including the renormalization of the Yukawa
couplings below $v^{\prime}$, is shown in Fig. 1.

Note that all mirror fermions are unstable except for $e^{\prime}$ and
$u^{\prime}$ which cannot decay since they are the lightest particles charged
under the unbroken $U(1)_{\rm EM}^{\prime}$ and $SU(3)^{\prime}$ symmetries,
respectively. Since $e^{\prime}$ is stable and has suppressed interactions
with SM particles, its relic cosmological abundance may provide the dark
matter. As discussed in Sec. 4.2, thermal freeze-out of $e^{\prime}$ leads to
the observed dark matter abundance if $v^{\prime}\simeq 8\times 10^{7}$ GeV,
setting our normalization for $v^{\prime}$ in all future equations. The mass
of $e^{\prime}$ is

$\displaystyle m_{e}^{\prime}=y_{e}v^{\prime}\simeq 225\;{\rm
GeV}\left(\frac{v^{\prime}}{8\times 10^{7}\,{\rm GeV}}\right).$ (10)

Note that while $u^{\prime}$ is also stable and a potential dark matter
candidate, its relic abundance is always subdominant to that of $e^{\prime}$.
This is because below the mirror QCD phase transition, mirror hadrons composed
of $s^{\prime}$, $u^{\prime}$ and $d^{\prime}$ efficiently annihilate into
$\gamma^{\prime}$ or decay into the hadron $u^{\prime}u^{\prime}u^{\prime}$
whose abundance is small compared to $e^{\prime}$ Dunsky:2019upk .

###  3.2 Axion

At energy scales above $v^{\prime}$ the $Z_{2}$ symmetry ensures that the QCD
and QCD′ gauge couplings are equal, $g_{3}=g_{3}^{\prime}$. However, below
$v^{\prime}$ the $Z_{2}$ symmetry is broken, and the large mass hierarchy
between the Standard and mirror fermions significantly changes the
renormalization flow of $g_{3}^{\prime}$ away from $g_{3}$. As the heavier
mirror fermions decouple in the IR, $g_{3}^{\prime}$ becomes larger than
$g_{3}$, leading to a higher confinement scale in the $SU(3)^{\prime}$ sector.
The value of $g_{3}$ and the renormalization scale $\mu=v^{\prime}=8\times
10^{7}$ GeV sets the initial value of $g_{3}^{\prime}$ at the same scale. We
then run the 2-loop beta function for $g_{3}^{\prime}$ Caswell:1974gg ;
Jones:1974mm ; Machacek:1983tz  to low energies, taking into account the
mirror quark thresholds. We define $\Lambda_{\rm QCD^{\prime}}$ as the scale
where $g_{3}^{\prime}(\mu)$ that is defined in the $\overline{MS}$ scheme and
computed at two-loop level diverges. Around $v^{\prime}\sim 10^{8}$ GeV, it is
given by

$\displaystyle\Lambda_{\rm QCD^{\prime}}\simeq 27{\,\rm
GeV}\left(\frac{v^{\prime}}{8\times 10^{7}\,\rm GeV}\right)^{4/11}\,.$ (11)

Consequently, the single axion in the theory, which couples to both Standard
Model and mirror sectors, acquires a large mass from mirror QCD instantons,

$\displaystyle m_{a}=0.6\frac{\Lambda_{\rm QCD^{\prime}}^{2}}{f_{a}}\simeq
9\,{\rm MeV}\left(\frac{v^{\prime}}{8\times 10^{7}\,\rm
GeV}\right)^{8/11}\left(\frac{f_{a}}{5\times 10^{4}\,\rm GeV}\right)^{-1}\,,$
(12)

where we have used the lattice calculations for the topological susceptibility
of $N_{f}=0$ flavor $SU(3)$ and the relation between $\Lambda_{\rm
QCD^{\prime}}$ and the Sommer scale Durr:2006ky ;
FlavourLatticeAveragingGroupFLAG:2021npn  to relate the mass of the axion with
$\Lambda_{\rm QCD^{\prime}}$. The axion mass is thus much larger than the
standard QCD axion mass $m_{a,0}\sim m_{\pi}f_{\pi}/f_{a}$. Moreover, the
mapping given in (2.1) ensures the phases of the potential generated by
$SU(3)$ and $SU(3)^{\prime}$ dynamics are aligned so that the strong CP
problem is still solved by the axion, despite its non-standard mass. The
heavier axion mass relaxes the PQ quality problem that plagues the standard
QCD axion Kamionkowski:1992mf ; Barr:1992qq ; Ghigna:1992iv ; Holman:1992us .

For $m_{a}\ll m_{\pi^{0}}\simeq 135$ MeV, the lifetime of the axion is set by
the $a\rightarrow 2\gamma$ and $a\rightarrow 2\gamma^{\prime}$ decay channels.
The total decay rate is

$\displaystyle\Gamma_{a}=\frac{g_{\gamma}^{2}+g_{\gamma^{\prime}}^{2}}{64\pi}\frac{m_{a}^{3}}{f_{a}^{2}}\simeq\left(1.5\times
10^{-1}\,\rm s\right)^{-1}\left(\frac{m_{a}}{9\,\rm
MeV}\right)^{3}\left(\frac{f_{a}}{5\times 10^{4}\,\rm GeV}\right)^{-2}\,,$
(13)

where the axion-photon and dark photon couplings are
$g_{\gamma}\simeq(E/N-5/3-\mathcal{F_{\theta}}(m_{a}))\alpha/2\pi$ and
$g_{\gamma^{\prime}}\simeq(E/N)\alpha^{\prime}/2\pi$, respectively. Here,
$E/N$ is the ratio of the charge to color anomaly of the PQ quarks, $5/3$
comes from the axial rotation that removes the axion-gluon coupling,
$\mathcal{F}_{\theta}$ from axion-meson mixing Aloni:2018vki , and
$\alpha\,(\alpha^{\prime})$ the (mirror) fine-structure constant. $E/N=2/3$
when the PQ quarks have the same charges as down-quarks and $8/3$ if they have
the same charges as up-quarks or in a complete $SU(5)$ multiplet. For
$m_{a}\ll m_{\pi}$, as occurs for the axion in most of our parameter space,
$(5/3+\mathcal{F}_{\theta})\simeq 2.03$. For Eq. (13) and the remainder of
this paper, we take $\alpha^{\prime}=\alpha$ due to the small renormalization
differences between the SM and SM′ electromagnetic sectors.

###  3.3 Mirror Glueballs

Mirror glueballs materialize after the QCD′ phase transition. Self scattering
quickly leaves the glueball gas in a state dominantly composed of the lightest
mirror glueball, $S^{\prime}$, with a mass Durr:2006ky

$\displaystyle m_{S^{\prime}}\simeq 6.8\;\Lambda_{\rm QCD^{\prime}}\simeq
185\,{\rm GeV}\left(\frac{v^{\prime}}{8\times 10^{7}\,{\rm
GeV}}\right)^{4/11}.$ (14)

The mirror glueball can decay into two Higgses, two mirror photons, or two
axions; it can also self-scatter with itself to reduce its abundance. Each
rate is strongly dependent on $v^{\prime}$, but for $v^{\prime}=10^{8}$ GeV,
the rates increasingly dominate in the aforementioned order. The decay rate
into electroweak gauge bosons is given by Dunsky:2019upk

$\displaystyle\Gamma_{S^{\prime}\rightarrow
W^{+}W^{-}/ZZ}\simeq\frac{1}{8\pi}\left(\frac{2.7}{16\pi^{2}}\right)^{2}\frac{m_{S^{\prime}}^{5}}{v^{\prime
4}}\simeq\left(11\,{\rm s}\right)^{-1}\left(\frac{v^{\prime}}{8\times
10^{7}\,\rm GeV}\right)^{-\frac{24}{11}}\,.$ (15)

Here we ignored the phase-space suppression and the actual decay rate is even
smaller. The decay rate into $2\gamma^{\prime}$ is given by

$\displaystyle\Gamma_{S^{\prime}\rightarrow\gamma^{\prime}\gamma^{\prime}}\simeq\frac{1}{16\pi}\left(\frac{2.7\alpha}{270\pi}\right)^{2}\frac{m_{S^{\prime}}^{9}}{m_{u^{\prime}}^{8}}\simeq\left(1.1\times
10^{-12}\,{\rm s}\right)^{-1}\left(\frac{v^{\prime}}{8\times 10^{7}\,\rm
GeV}\right)^{-\frac{52}{11}}\,.$ (16)

The decay rate into $2a$ is given by

$\displaystyle\Gamma_{S^{\prime}\rightarrow
aa}\approx\frac{1}{8\pi}\frac{m_{S^{\prime}}^{5}}{f_{a}^{4}}\simeq\left(4.8\times
10^{-16}\,{\rm s}\right)^{-1}\left(\frac{f_{a}}{5\times 10^{4}\,\rm
GeV}\right)^{-4}\left(\frac{v^{\prime}}{8\times 10^{7}\,\rm
GeV}\right)^{\frac{20}{11}}\,.$ (17)

Last, the annihilation of two lightest glueballs (CP-even states) into an
axion and the lightest CP-odd glueball is of similar importance in reducing
the glueball abundance, and has a rate

$\displaystyle\Gamma_{S^{\prime}S^{\prime}\rightarrow S^{\prime}a}$
$\displaystyle\sim\frac{\left(m_{S^{\prime}}T\right)^{\frac{3}{2}}}{f_{a}^{2}}e^{-\frac{m_{S}^{\prime}}{T}}$
$\displaystyle\approx\left(7.1\times 10^{-22}\,{\rm
s}\right)^{-1}\left(\frac{f_{a}}{5\times 10^{4}\,\rm
GeV}\right)^{-2}\left(\frac{v^{\prime}}{8\times 10^{7}\,\rm
GeV}\right)^{\frac{12}{11}}\left(\frac{T}{m_{S}^{\prime}}\right)^{\frac{3}{2}}e^{-\frac{m_{S}^{\prime}}{T}}.$
(18)

###  3.4 Mirror Neutrinos and Neutrinos

A successful theory of $e^{\prime}$ dark matter from freeze-out has important
implications for neutrino masses because late decays of mirror neutrinos can
dilute the $e^{\prime}$ and can affect nucleosynthsis. Parity implies that the
SM and SM′ Weinberg operators give correlated masses for the SM and mirror
neutrinos. In addition, the neutrino portal operator of (3) gives a neutrino
Yukawa coupling. Hence, below $v^{\prime}$ the relevant EFT for neutrino
masses is

$\displaystyle{\cal
L}_{\nu,\nu^{\prime}}=\frac{1}{2}m_{\nu^{\prime}_{i}}\left(\nu^{\prime}_{i}\nu^{\prime}_{i}+\ell_{i}\ell_{i}\,\frac{HH}{v^{\prime
2}}\right)+\nu^{\prime}_{i}\,y_{ij}\,\ell_{j}\,H+{\rm h.c.},\qquad
y_{ij}=\frac{v^{\prime}}{M_{D}}\xi_{ij}.$ (19)

The light neutrino mass matrix is therefore

$\displaystyle m_{\nu_{ij}}=\frac{v^{2}}{v^{\prime
2}}\,m_{\nu^{\prime}_{i}}\;\delta_{ij}-y^{T}_{ik}\,\frac{1}{m_{\nu^{\prime}_{k}}}\,y_{kj}\;v^{2}.$
(20)

We call the first term the “direct” contribution and the second the “seesaw”
contribution.

It is useful to study the direct and seesaw contributions to the light
neutrino mass matrix from each $\nu^{\prime}_{i}$. Each $\nu^{\prime}_{i}$
couples to a single combination of $\ell_{j}$, which we call
$\tilde{\ell}_{i}$, so that the Yukawa coupling of $\nu^{\prime}_{i}$ can be
written as

$\displaystyle{\cal
L}_{y_{i}}=y_{i}\;\nu^{\prime}_{i}\,\tilde{\ell_{i}}\,H+{\rm
h.c.},\qquad\tilde{\ell_{i}}=y_{ij}\,\ell_{j}/y_{i},\qquad
y_{i}^{2}\equiv\sum_{j}|y_{ij}|^{2}.$ (21)

Thus, each $\nu^{\prime}_{i}$ gives mass contributions to two different
states, a direct one for $\nu_{i}$ and a seesaw one for $\tilde{\nu}_{i}$. If
large leptonic mixing angles arise from the neutrino sector, these two states
are expected to be very different, although generically they are not
orthogonal. Consequently, the Lagrangian for the light neutrino masses can be
written as a sum of three such terms, one from each $\nu^{\prime}_{i}$

$\displaystyle{\cal L}_{\nu}=\frac{1}{2}\sum_{i}\left(m_{\nu_{i}}^{\rm
dir}\;\nu_{i}\nu_{i}-m_{\nu_{i}}^{\rm
ss}\;\tilde{\nu}_{i}\tilde{\nu}_{i}\right)+{\rm
h.c.}=\frac{1}{2}\sum_{i}\left(\frac{v^{2}}{v^{\prime
2}}\;m_{\nu^{\prime}_{i}}\;\nu_{i}\nu_{i}-\;\frac{y_{i}^{2}v^{2}}{m_{\nu^{\prime}_{i}}}\;\tilde{\nu}_{i}\tilde{\nu}_{i}\right)+{\rm
h.c.}.$ (22)

The resulting neutrino mass matrix therefore has 6 terms and depends on
$(v^{\prime},m_{\nu^{\prime}_{i}},y_{i})$ as well as the mixing angles and
phases in $y_{ij}/y_{i}$.

The orange and blue contours of Fig. 2 show the direct and seesaw
contributions from one $\nu^{\prime}_{i}$ in the
$(y_{i},m_{\nu^{\prime}_{i}})$ plane. For the direct contribution, the scale
$v^{\prime}$ of parity breaking is taken to be $8\times 10^{7}$ GeV, as
required for $e^{\prime}$ to account for the observed dark matter via
freezeout. The mass of $\nu^{\prime}_{i}$ is directly proportional to the
direct neutrino mass contribution

$\displaystyle m_{\nu_{i}^{\prime}}=11\,{\rm GeV}\left(\frac{m_{\nu_{i}}^{\rm
dir}}{0.05\,\rm eV}\right)\left(\frac{v^{\prime}}{8\times 10^{7}\,{\rm
GeV}}\right)^{2}\,.$ (23)

In the orange (blue) shaded regions the direct (seesaw) contribution is larger
than 0.05 eV, so a realistic spectrum requires a cancellation among the 6
contributions. Deeper into the shaded region the required cancellation becomes
stronger. Deep in the unshaded region, both direct and seesaw contributions
are too small to affect current data. There are three copies of this figure,
one for each $\nu^{\prime}_{i}$. Assuming no strong cancellations, at least 2
of the 6 contributions must be close to the edge of the corresponding shaded
region to account for the atmospheric and solar oscillations. Two relevant
contributions could come from a single $\nu^{\prime}_{i}$ if the values of
$(y_{i},m_{\nu^{\prime}_{i}})$ lie close to to the right-hand vertex of the
unshaded triangle, where the 0.05 eV contours of the direct and seesaw
contributions intersect. This intersection point gives the maximum value of
$y_{i}$ possible without cancellations, $y_{i}\lesssim
m_{\nu^{\prime}_{i}}/v^{\prime}$, and similarly the maximum value of
$|y_{ij}|$ for any $j$

$\displaystyle|y_{ij}|\lesssim\frac{m_{\nu^{\prime}_{i}}}{v^{\prime}}\approx
1\times 10^{-7}\left(\frac{m_{\nu}}{0.05\,\rm
eV}\right)\left(\frac{v^{\prime}}{8\times 10^{7}\,\rm GeV}\right)\,.$ (24)

Figure 2: Contours for the direct (orange) and seesaw (blue) contributions to
the light neutrino mass matrix, from the heavy mirror neutrino mass eigenstate
$\nu^{\prime}_{i}$, as defined in Eq. (22), with $v^{\prime}=8\times 10^{7}$
GeV. The orange (blue) shaded regions give $m_{\nu_{i}}^{\rm
dir}(m_{\nu_{i}}^{\rm ss})$ larger than 0.05 eV, and in many theories require
unnatural cancellations. The red line has $m_{\nu_{i}}^{\rm
dir}=m_{\nu_{i}}^{\rm ss}$, giving massless light neutrinos at tree level. The
Singlet Model is constrained to live on this line for all three
$\nu^{\prime}_{i}$. In this theory, radiative neutrino masses of 0.05 eV can
result along the solid part of the red line.

Fig. 2 shows that the $\nu^{\prime}_{i}$ have masses far below $v^{\prime}$,
so that decays via virtual $W^{\prime}$ are negligible and $\nu^{\prime}_{i}$
have long lifetimes. Decays proceed dominantly through the Yukawa coupling
$y_{i}$ of (21), which induces a small mixing angle between $\nu^{\prime}_{i}$
and $\tilde{\nu}_{i}$ of

$\displaystyle\theta_{i}\simeq\frac{y_{i}v}{m_{\nu^{\prime}_{i}}}.$ (25)

For $m_{\nu^{\prime}_{i}}\mathop{}_{\textstyle\sim}^{\textstyle<}m_{W}+m_{e}$,
the dominant decay modes are the three-body beta decays
$\nu^{\prime}_{i}\rightarrow\tilde{e}_{i}\,u_{j}\,\bar{d}_{j},\,\tilde{e}_{i}\,\nu_{j}\,\bar{e}_{j}$,
with $j$ running over two generations of quarks and three generations of
leptons, giving

$\displaystyle\Gamma_{\nu^{\prime}_{i}}\;\simeq\;9\frac{1}{192\pi^{3}}\;\theta_{i}^{2}\;\frac{m_{\nu^{\prime}_{i}}^{5}}{8v^{4}}\,.$
(26)

Similarly, there are $Z$-mediated decay modes. For
$m_{\nu^{\prime}_{i}}\gtrsim m_{W}+m_{e}$, the two-body decay
$\nu^{\prime}_{i}\rightarrow\tilde{\ell}_{i}H$ becomes kinematically allowed
and dominates, with a decay rate

$\displaystyle\Gamma_{\nu^{\prime}_{i}}\;\simeq\;\frac{1}{8\pi}\,y_{i}^{2}\;m_{\nu^{\prime}_{i}}\,.$
(27)

####  3.4.1 Neutrino masses constrained by lepton symmetries

The general form of neutrino masses, (19), has two important limits, one where
$m_{\nu^{\prime}_{i}}=0$ and one where $y_{ij}=0$. This changes the physics
drastically, and Fig. 2 is no longer relevant.

Imposing $B-L$ number, with opposite signs in the SM and mirror sectors, sets
$m_{\nu^{\prime}_{i}}=0$ so that the observed neutrinos are Dirac. The heavy,
unstable $\nu^{\prime}$ states do not exist, changing the cosmological
evolution of the universe. Unlike the charged fermions, where the right-handed
states are part of the SM sector, neutrinos have their right-handed components
coming from the mirror sector, allowing an alternative view of why the
neutrinos are so light. The observed neutrino masses result from Yukawa
couplings $y_{ij}$ of order $10^{-12}$. However, unlike for charged fermions,
these Yukawa interactions arise from dimension 5 operators above $v^{\prime}$,
shown in (3). With $v^{\prime}=10^{8}$ GeV and $M_{D}$ the Planck scale, this
requires $\xi_{ij}$ of order $10^{-2}$, comparable to the largest Yukawa
coupling in the charged lepton sector.

Alternatively, separate lepton parities can be imposed in the SM and mirror
sectors, removing the neutrino portal operator of (3) and forcing $y_{ij}=0$.
The neutrino sector then conserves individual lepton parities for each of the
three generations, and the SM and mirror neutrino masses are related by
$m_{\nu_{i}}=(v^{2}/v^{\prime 2})m_{\nu^{\prime}_{i}}$. Lepton mixing angles
arise because in this basis the Yukawa coupling matrix for the charged leptons
is not diagonal. For $v^{\prime}=10^{8}$ GeV, the lightest mirror neutrino,
$\nu^{\prime}_{1}$, is stable as it is the lightest mirror fermion. The two
heavier mirror neutrinos decay via a loop diagram involving a virtual
$W^{\prime}$, $\nu^{\prime}_{2,3}\rightarrow\nu^{\prime}_{1}+\gamma^{\prime}$.
The decay rate is highly suppressed, giving lifetimes of order $10^{9}$ sec.

####  3.4.2 The singlet model for neutrino masses

In many theories, significant cancellations between different contributions to
the light neutrino mass matrix can only arise from fine tuning. In these
theories the shaded regions of Fig. 2 are fine-tuned, and (24) are naturalness
constraints. However, cancellations cannot be avoided in one of the simplest
UV completions of the dimension 5 operators for neutrino masses Hall:2019qwx ,
which we call the Singlet Model.

Consider a theory where the dimension 5 operators for neutrino masses are
generated by the exchange of three gauge-singlet Weyl fermions $S_{i}$, that
are parity even: $S_{i}\leftrightarrow S_{i}$, via the interactions

$\displaystyle{\cal
L}(S_{i})=S_{i}\,x_{ij}\,(\ell_{j}H+\ell^{\prime}_{j}H^{\prime})+\frac{1}{2}M_{S_{i}}\,S_{i}S_{i}+{\rm
h.c.}$ (28)

Introducing a hatted basis, the Yukawa interactions of $S_{i}$ can be written
as

$\displaystyle{\cal
L}_{x_{i}}=x_{i}\;S_{i}\,(\hat{\ell}_{i}\,H+\hat{\ell^{\prime}}_{i}\,H^{\prime})+{\rm
h.c.},\qquad\hat{\ell_{i}}=x_{ij}\,\ell_{j}/x_{i},\qquad
x_{i}^{2}\equiv\sum_{j}|x_{ij}|^{2},$ (29)

in analogy with (21). The EFT below $M_{S_{i}}$ is

$\displaystyle{\cal
L}_{EFT}(S_{i})\;=\;\frac{1}{2}\,\sum_{i}\frac{x_{i}^{2}}{M_{S_{i}}}\;(\hat{\ell}_{i}\,H+\hat{\ell^{\prime}}_{i}\,H^{\prime})^{2}+{\rm
h.c.}.$ (30)

Note that the $\hat{\ell}_{i}$ are not orthogonal.

Inserting the $H^{\prime}$ vev, and comparing with Eqns. (19) and (21), the
singlet model leads to a correlation between
$m_{\nu^{\prime}_{i}}=x_{i}^{2}v^{\prime 2}/M_{S_{i}}$ and
$y_{i}=x_{i}^{2}v^{\prime}/M_{S_{i}}$. Thus it gives
$y_{i}=m_{\nu^{\prime}_{i}}/v^{\prime}$, leading to

$\displaystyle m_{\nu_{i}}^{\rm dir}=m_{\nu_{i}}^{\rm
ss}=\frac{y_{i}^{2}v^{2}}{m_{\nu^{\prime}_{i}}}=\frac{x_{i}^{2}v^{2}}{M_{S_{i}}}$
(31)

giving the red line in Fig. 2. The light neutrinos are massless at tree level
because each $S_{i}$ couples to only one combination of mirror neutrinos and
neutrinos, $v^{\prime}\hat{\nu}^{\prime}_{i}+v\,\hat{\nu}_{i}$, leaving the
orthogonal combination, $v^{\prime}\hat{\nu}_{i}-v\,\hat{\nu}^{\prime}_{i}$,
without a mass term. This single coupling of $S_{i}$ leads to the combination
$(\hat{\ell}_{i}\,H+\hat{\ell^{\prime}}_{i}\,H^{\prime})^{2}$ appearing in
(30), implying that the coefficient of the Weinberg operators involving
$\hat{\ell}_{i}\hat{\ell}_{i}$ and
$\hat{\ell}^{\prime}_{i}\hat{\ell}^{\prime}_{i}$ are correlated with that of
the portal operator involving $\hat{\ell}^{\prime}_{i}\hat{\ell}_{i}$.
However, this correlation is lost under renormalization group scaling induced
by one-loop electroweak radiative corrections. Hence the light neutrino mass
arising from the interaction of $S_{i}$ is estimated to be

$\displaystyle{\cal
L}_{\nu}(S_{i})\;\simeq\;\frac{1}{2}\;\frac{y_{i}^{2}\,v^{2}}{m_{\nu^{\prime}_{i}}}\;\frac{g_{2}^{2}}{16\pi^{2}}(2L_{+}+L_{-})\;\hat{\nu}_{i}\hat{\nu}_{i},\qquad
L_{+}=\ln\frac{M_{S_{i}}}{v^{\prime}},\;L_{-}=\ln\frac{v^{\prime}}{m_{\nu^{\prime}_{i}}}.$
(32)

The two logs correspond to running in the EFTs above and below $v^{\prime}$.
Since $L_{-}$ is always a large log, the loop factor is expected to be in the
range of 3-30, giving the extent of the solid section of the red line in Fig.
2.

##  4 Mirror Cosmology

After inflation, the maximum temperature of the universe is taken to be less
than $v^{\prime}$, so that the spontaneous breaking of parity does not lead to
unacceptable domain wall densities, but greater than the freeze-out
temperature of mirror electrons, allowing a reliable computation of the relic
$e^{\prime}$ density.

###  4.1 Portals and Decoupling Temperatures

Thermalization of the SM and SM′ sectors can be achieved through the Higgs,
kinetic mixing, neutrino, and heavy axion portals arising from the operators
in Eqns. (3) and (5). We shall assume $\epsilon$ is negligibly small so that
kinetic mixing is unimportant from here on.111The leading non-zero radiative
correction to $\epsilon$ arises at the five-loop level from the interactions
between the axion field $P$, the PQ quarks $\Psi,\Psi^{\prime}$, and Higgses.
Its value is of order $\epsilon_{\rm rad}\sim
1/(16\pi^{2})^{5}|\xi|^{4}|y_{D}|^{4}N_{c}^{2}$. For
$|y_{\Psi}|^{4}|y_{D}|^{4}\lesssim 0.1$, Rutherford scattering of $e^{\prime}$
DM with nuclei in the LZ detector is below detection threshold for 15 ton-
years of exposure Dunsky:2018mqs .

The Higgs and neutrino portals give rise to similar decoupling temperatures
which increase with $v^{\prime}$ and are independent of $f_{a}$. On the other
hand, the axion portal gives rise to a decoupling temperature which increases
with $f_{a}$ but is independent of $v^{\prime}$ for temperatures above the
QCD′ phase transition:

$\displaystyle T_{\rm dec}\simeq\begin{dcases}5\times 10^{4}\,{\rm
GeV}\left(\frac{v^{\prime}}{1\times 10^{8}\,{\rm
GeV}}\right)^{4/3}\quad&\text{Higgs and Neutrino
Portals}\,\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Dunsky:2018mqs}{\@@citephrase{(}}{\@@citephrase{)}}}}\\\
5\times 10^{4}\,{\rm GeV}\left(\frac{f_{a}}{1\times 10^{9}\,{\rm
GeV}}\right)^{2.35}\quad&\text{Axion Portal~{}
\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Salvio:2013iaa,DEramo:2021psx,DEramo:2021lgb}{\@@citephrase{(}}{\@@citephrase{)}}}}}\\\
\end{dcases}$ (33)

As discussed in Sec. 4.2, thermally producing $e^{\prime}$ with the observed
dark matter abundance requires $v^{\prime}\approx 10^{8}$ GeV. From Eq. (33),
we thus see that for $f_{a}\lesssim 10^{9}$ GeV, the axion portal dominates
over the Higgs and Neutrino portals in keeping the two sectors in thermal
equilibrium. Moreover, as discussed in Sec. 5.1, to ensure that the axion does
not generate significant dark radiation around or after neutrino decoupling
requires $f_{a}\lesssim 10^{5}$ GeV. Consequently, in our cosmology, the axion
is the dominant portal responsible for thermal coupling between the two
sectors.

Specifically, the axion keeps the two sectors in equilibrium both above and
below the mirror QCD phase transition temperature,

$\displaystyle T_{\rm QCD^{\prime}}\;\simeq\;1.26\,\Lambda_{\rm
QCD^{\prime}}\;\simeq\;34{\,\rm GeV}\left(\frac{v^{\prime}}{8\times
10^{7}\,\rm GeV}\right)^{4/11},$ (34)

where we used the relation between $T_{\rm QCD^{\prime}}$ and $\Lambda_{\rm
QCD^{\prime}}$ computed in Borsanyi:2012ve . Above $T_{\rm QCD^{\prime}}$,
$a\leftrightarrow gg$ and $a\leftrightarrow g^{\prime}g^{\prime}$ scatterings
keep the two sectors in equilibrium. Below the QCD phase transition
temperature, the two sectors are kept in thermal equilibrium through Primakoff
scatterings, $\gamma e\leftrightarrow ae$ and
$\gamma^{\prime}e^{\prime}\leftrightarrow ae^{\prime}$, until the $e^{\prime}$
number density drops exponentially at temperatures below $m_{e^{\prime}}$.

The first interaction, that between axions and mirror gluons, can be
calculated precisely from the standard axion-gluon scattering rate as computed
in Salvio:2013iaa ; DEramo:2021lgb  with the mapping $g_{3}\leftrightarrow
g_{3}^{\prime}$,

$\displaystyle\Gamma_{ag^{\prime}\leftarrow
g^{\prime}g^{\prime}}\simeq\frac{16}{\pi}\left(\frac{g_{3}^{\prime
2}}{32\pi^{2}}\right)^{2}\frac{T^{3}}{f_{a}^{2}}\mathcal{F}_{g^{\prime}}(g_{3}^{\prime})\,.$
(35)

Here, $\mathcal{F}_{g^{\prime}}$ is a function of $g_{3}^{\prime}$ and is
numerically computed in Salvio:2013iaa ; DEramo:2021lgb  but takes the form
$\mathcal{F}_{g^{\prime}}\approx 2g_{3^{\prime}}^{2}\ln 1.5/g_{3}^{\prime}$
for $g_{3^{\prime}}\lesssim 1$. 222 Note that Eq. (35) is strictly valid in
the limit $m_{a}\ll T$, which is always the case for temperatures above
$T_{\rm QCD^{\prime}}$. When using Eq. (35), we use the temperature-dependent
value of $g_{3}^{\prime}$ at the renormalized energy scale $\mu=T$.

Similarly, the mirror Primakoff scattering rate is given by the Standard Model
rate Bolz:2000fu ; Cadamuro:2011fd  with the replacement of the electron
number density with the mirror electron density, $n_{e}\rightarrow
n_{e^{\prime}}$, as given by

$\displaystyle\Gamma_{ae^{\prime}\leftrightarrow
e^{\prime}\gamma^{\prime}}\simeq\alpha^{\prime}\left(\frac{g_{\gamma^{\prime}}}{f_{a}}\right)^{2}\frac{\pi^{2}}{36\zeta(3)}\left(\ln\frac{T^{2}}{m_{\gamma}^{\prime
2}}+0.82\right)n_{e^{\prime}}\,,$ (36)

where $g_{\gamma^{\prime}}\simeq(E/N)\alpha^{\prime}/2\pi$ and
$m_{\gamma}^{\prime 2}\approx
4\pi\alpha^{\prime}n_{e^{\prime}}/m_{e^{\prime}}$ is the mirror photon plasma
mass. Because the renormalization in the electromagnetic gauge couplings are
small, we take $\alpha^{\prime}=\alpha$.

We define the temperature when the SM and SM’ sectors decouple by the
temperature at which the rate $\Gamma_{a\leftrightarrow\rm{SM^{\prime}}}=3H$,
where $H$ is Hubble and $\Gamma_{a\leftrightarrow\rm{SM}^{\prime}}$ is the sum
of the rates given in Eqns. (35) and (36). The axion decoupling temperature
from the mirror Standard Model is shown in Fig. 3 as a function of $f_{a}$ for
fixed $v^{\prime}=8\times 10^{7}$ GeV ($\Lambda_{\rm QCD^{\prime}}\approx 27$
GeV, $m_{e}^{\prime}\simeq 225$ GeV), which is the required $v^{\prime}$ to
freeze-out mirror electrons as dark matter (see Sec. 4.2). The solid orange
curve shows the axion decoupling temperature from mirror gluon interactions,
Eq. (35), which flatlines at $T_{\rm QCD^{\prime}}\approx 34$ GeV when
$g_{3^{\prime}}$ becomes non-perturbative. Note that since $u^{\prime}$ and
$d^{\prime}$ are much heavier than $\Lambda_{\rm QCD^{\prime}}$, mirror pions
are heavy and thus severely Boltzmann suppressed in number density immediately
after the mirror QCD phase transition. Thus, we do not expect mirror quark
bound states to significantly reduce the decoupling temperature below $T_{\rm
QCD^{\prime}}$ even for $f_{a}$ below $\sim 10^{8}$ GeV. Consequently, we
extrapolate the decoupling temperature due to mirror gluons to lower $f_{a}$
by the dashed horizontal orange line. Similarly, the blue contour of Fig. 3
shows the axion decoupling temperature from mirror Primakoff interactions. We
see that for $f_{a}\lesssim 10^{5}$ GeV, the Primakoff interactions can
dominate, keeping the two sectors in equilibrium to temperatures slightly
below $\Lambda_{\rm QCD^{\prime}}$. Since $f_{a}\gtrsim 10^{5}$ GeV produces
long-lived axions that decay and generate too large $|\Delta N_{\rm eff}|$ as
discussed in Sec. 5.1, we will generally be in the low $f_{a}$ region where
Primakoff interactions set the decoupling temperature between the two sectors
to $T_{\rm decouping}\approx 20-30$ GeV.

Figure 3: Decoupling temperature of SM and mirror sectors from the axion
portal as a function of $f_{a}$. The orange curve shows the decoupling
temperature from axion and mirror gluon interactions which become ineffective
around the mirror QCD phase transition at $T_{\rm QCD^{\prime}}\approx 34$
GeV. The blue curve shows the decoupling temperature from axion and mirror
electron scatterings (Primakoff) which become ineffective at temperatures
below $m_{e}^{\prime}\approx 225$ GeV . For $f_{a}\lesssim 10^{5}$ GeV, the
decoupling temperature from either interaction are comparable and occurs
around $20-30$ GeV.

###  4.2 $e^{\prime}$ Dark Matter

The stability of $e^{\prime}$ and $u^{\prime}$ guarantees the existence of
dark matter in the theory. The freeze-out chronology is similar to
Dunsky:2019upk , except for the existence of the axion, which maintains
thermal equilibrium between the two sectors to far lower temperatures as
discussed in the previous section. An overview of the mirror sector cosmology
is shown on the temperature axis on the right-hand side of Fig. 1, for
$v^{\prime}=8\times 10^{7}$ GeV, which we justify below. For a sufficiently
high reheating temperature after inflation, the following events occur in
succession:

  1. 1.

$t^{\prime}$, $b^{\prime}$, $\tau^{\prime}$, $c^{\prime}$, $\mu^{\prime}$,
$s^{\prime}$ and $d^{\prime}$ freeze-out, and $t^{\prime}$, $b^{\prime}$,
$\tau^{\prime}$ and $c^{\prime}$ subsequently decay. The mirror neutrinos
$\nu^{\prime}$ decouple while relativistic. During this early era, axions keep
the SM and SM’ sectors in equilibrium since axion-gluon and axion-mirror gluon
scatterings are in equilibrium.

  2. 2.

The mirror QCD phase transition occurs at $T_{\rm QCD^{\prime}}\simeq 34$ GeV.
$u^{\prime}$, $d^{\prime}$ and the residual frozen-out $s^{\prime}$ confine to
mirror hadrons which quickly annihilate to $\gamma^{\prime}$ or decay to the
lightest stable hadron, $u^{\prime}u^{\prime}u^{\prime}$ Kang:2006yd ;
Harigaya:2016nlg ; DeLuca:2018mzn . Any residual quarks containing a PQ
component also confine and efficiently annihilate, reducing their abundance
far below that of dark matter and thus cosmologically harmless when they
eventually decay.

  3. 3.

Mirror glueballs ($S^{\prime}$) also form at the QCD′ phase transition and
carry most of the energy and entropy of the mirror gluons. $S^{\prime}$ decays
and inverse decays to $a$ are much faster than Hubble, keeping $S^{\prime}$ in
equilibrium to temperatures $T\ll m_{S^{\prime}}$. The energy and entropy of
$S^{\prime}$ are thus transferred to the Standard and mirror baths which are
still in equilibrium with $a$ due to axion-gluon and axion-mirror Primakoff
interactions, respectively.

  4. 4.

$e^{\prime}$ freezes-out. The exponential reduction in the number density of
$e^{\prime}$ causes the $\gamma^{\prime}$ to decouple from $a$ due to the
reduction in the axion-mirror Primakoff rate. This marks the decoupling of the
Standard and Mirror baths, around $T\sim m_{e^{\prime}}/10\sim(20-30)$ GeV
(see the blue curve of Fig. 3).

  5. 5.

The axion remains in equilibrium with the Standard Model bath due to
interactions with gluons and then pions until $T\sim 10$ MeV.

  6. 6.

$\nu^{\prime}$ and $a^{\prime}$ decay.

Analytically, we can expect that the critical density of $e^{\prime}$ will
scale as Kolb:1990vq

$\displaystyle\Omega_{\rm DM}h^{2}\approx\frac{x_{\rm FO}}{M_{\rm
Pl}\,\sigma_{e^{\prime}\bar{e^{\prime}}\leftrightarrow\gamma^{\prime}\gamma^{\prime}}}\frac{\sqrt{g_{*}(T_{\rm
FO}})}{g_{*S}(T_{\rm FO})}\frac{1}{\rm eV}\simeq
0.12\left(\frac{v^{\prime}}{8\times 10^{7}\,\rm GeV}\right)^{2}$ (37)

where $x_{\rm FO}\equiv m_{e^{\prime}}/{T_{\rm FO}}\simeq 23$ is the ratio of
the mirror electron mass to its freeze-out temperature, and
$\sigma_{e^{\prime}\bar{e^{\prime}}\leftrightarrow\gamma^{\prime}\gamma^{\prime}}\simeq\pi\alpha/m_{e}^{\prime
2}$ is the mirror electron annihilation cross-section into mirror photons.
Because $\Omega_{\rm DM}\propto
1/\sigma_{e^{\prime}\bar{e^{\prime}}\leftrightarrow\gamma^{\prime}\gamma^{\prime}}\propto
m_{e^{\prime}}^{2}$, we expect $\Omega_{\rm DM}\propto v^{\prime 2}$ as shown
in Eq. (37), with $v^{\prime}\simeq 10^{8}$ GeV giving $e^{\prime}$ the
observed dark matter abundance.

In this work, we numerically solve a system of Boltzmann equations describing
the cosmology of $e^{\prime}$ and $u^{\prime}$ freeze-out in order to
precisely determine the $v^{\prime}$ that gives the correct $e^{\prime}$ dark
matter abundance. We find that Eq. (37) is an excellent approximation and that
the relic abundance of $e^{\prime}$ matches that of the observed dark matter
abundance for $v^{\prime}=8\times 10^{7}$ GeV. This justifies using the value
of $v^{\prime}$ in Fig. 1. In addition, we find that the $u^{\prime}$
abundance is subdominant compared to $e^{\prime}$ due to its bound-state
formation and thus enhanced annihilations at $T^{\prime}_{\rm QCD}$. We use
the cross-sections as given in Appendix A of Dunsky:2019upk  to compute this.

Note that in solving the Boltzmann equations, we require that the
$\nu^{\prime}$ do not dominate the energy density of the Universe before
decaying. This condition requires the neutrino Yukawa couplings, $|y_{ij}|$,
to be sufficiently large that the $\nu^{\prime}$ lifetime (26) is not
exceedingly long. Specifically, to avoid $\nu^{\prime}$ matter domination, the
temperature at which the $\nu^{\prime}$ decay
($T_{\Gamma,\nu^{\prime}}\sim\sqrt{M_{\rm Pl}\Gamma_{\nu^{\prime}}}$) must be
greater than the temperature at which they would begin dominating the energy
density of the universe ($T_{{\rm MD},\nu^{\prime}}\sim m_{\nu^{\prime}}Y_{\rm
therm}$),

$\displaystyle D^{-1}\equiv\frac{T_{\Gamma,\nu^{\prime}}}{T_{{\rm
MD},\nu^{\prime}}}=\frac{\sqrt{\Gamma_{\nu^{\prime}}M_{\rm
Pl}\sqrt{\frac{10}{8\pi^{3}g_{*}}}}}{m_{\nu^{\prime}}Y_{\rm
therm}}>1\qquad\text{(No $e^{\prime}$ dilution from $\nu^{\prime}$ decay)}$
(38)

for each of the three $\nu^{\prime}_{i}$. If $D$ is larger than unity, the
$e^{\prime}$ abundance is diluted by $\nu^{\prime}$ decay by a factor $D$; for
$e^{\prime}$ to account for all dark matter requires
$v^{\prime}\approx\sqrt{D}\,10^{8}$ GeV.

In Fig. 4, the region where $\nu^{\prime}_{i}$ decay leads to dilution, $D>1,$
is shaded in orange in the $(y_{i},m_{\nu^{\prime}_{i}})$ plane. Here an
initial thermal abundance of $\nu^{\prime}$ is assumed, as results from
virtual $W^{\prime}$ exchange if the reheat temperature of the universe after
inflation is above $10^{6}$ GeV, for $v^{\prime}$ of $10^{8}$ GeV. Thus,
$e^{\prime}$ freeze-out yields the observed dark matter abundance above and to
the right of this region for $v^{\prime}=10^{8}$ GeV. In the orange region,
$\nu^{\prime}_{i}$ come to dominate the energy density of the universe before
they decay; the $e^{\prime}$ abundance is diluted, requiring larger values of
$v^{\prime}$ to explain the observed dark matter, as studied in section 6.

Figure 4:  $\nu^{\prime}$ masses and Yukawa couplings that yield a successful
cosmology, with the observed dark matter resulting from $e^{\prime}$ freeze-
out with $v^{\prime}=8\times 10^{7}$ GeV. The orange region is excluded from
$\nu^{\prime}$ generating entropy upon its decay and diluting the $e^{\prime}$
abundance below that of dark matter while the green region is excluded from
$\nu^{\prime}$ decaying below $4$ MeV and disrupting nuclear abundances as
predicted by standard BBN. The blue region is excluded by naturalness of the
$\nu$ masses.

If $\nu^{\prime}_{i}$ have significant abundances, they cannot decay during or
after Big Bang Nucleosynthesis (BBN) as this leads to unacceptable nuclear
abundances Kawasaki:2004qu . The green region of Fig. 4 shows the parameter
space where $\nu^{\prime}_{i}$ decay at temperatures below $T_{\rm BBN}=4$
MeV. As the reheat temperature of the universe is lowered below $10^{6}$ GeV,
so that the $\nu^{\prime}_{i}$ abundances drop below the thermal one; the
orange region shrinks rapidly, but the green region is very robust because the
strong BBN bound. Hence, for a wide range of reheat temperatures, a consistent
cosmology for $e^{\prime}$ dark matter with $v^{\prime}=8\times 10^{7}$ GeV
results only in the region above and to the right of the green region.

On the solid blue contours of Fig. 4, the direct and seesaw contributions to
light neutrino masses, as given in (22), are 0.05 eV, so that natural neutrino
masses arise in the region between these contours. In the blue shaded regions,
unnatural fine tunings must be made in the light neutrino mass matrix; for
example a cancellation of a direct mass from one $\nu^{\prime}$ with the
seesaw contribution generated by another $\nu^{\prime}$. The beginning of the
blue shading marks the parameter space where the cancellations require a fine
tuning by a factor of three.

In summary, for this cosmology to be successful and yield $e^{\prime}$ dark
matter with $v^{\prime}=8\times 10^{7}$ GeV, all three $\nu^{\prime}_{i}$ must
lie in the unshaded region of Fig. 4: in the orange region $e^{\prime}$ are
diluted by a thermal abundance of $\nu^{\prime}_{i}$, in the green region late
decaying $\nu^{\prime}_{i}$ adversely affect BBN for a wide range of
$\nu^{\prime}_{i}$ abundances, and in the blue region fine-tuning occurs in
the neutrino mass matrix. This implies $m_{\nu^{\prime}_{i}}\simeq(0.8-4.0)$
GeV and $y_{i}\simeq(0.8-4.0)\times 10^{-8}$. The mirror neutrinos decay well
after the electroweak phase transition, so they do not give a baryon asymmetry
via leptogenesis. The light neutrino mass matrix has relevant contributions
from three seesaw terms and three direct terms. The direct masses are all
larger than about $0.05\,\rm{eV}/10$ and the seesaw masses are all larger than
about $0.05\,\rm{eV}/100$. If the reheating temperature of the universe is
below $10^{6}$ GeV, the constraint given by the orange region may be avoided.

While the BBN constraint is quite robust, it is possible to evade if the
abundance of $\nu^{\prime}_{i}$ is sufficiently low. This requires the reheat
temperature after inflation to be less than about 100 GeV to suppress freeze-
in production via virtual $W^{\prime}$. This is still consistent with
$e^{\prime}$ freeze-out and with gravity wave signals from the QCD′ phase
transition. It also requires
$y_{i}\mathop{}_{\textstyle\sim}^{\textstyle<}10^{-15}$ to suppress freeze-in
production from the Higgs portal. Thus one or more $\nu^{\prime}_{i}$ could
lie to the far left of an extension of Fig. 4, with extremely low $y_{i}$ and
lifetimes longer than $10^{9}$s.

Of particular interest is the possibility that the neutrinos are Dirac, as
discussed in section 3.4.1. In this case the heavy $\nu^{\prime}$ states are
absent and hence cannot dilute $e^{\prime}$ or distort element abundances
during BBN; Fig. 4 is not applicable. With $v^{\prime}=8\times 10^{7}$ GeV,
mirror electrons account for the observed dark matter for any values of the
neutrino masses and mixings consistent with current data.

Alternatively, neutrino masses may be constrained by separate lepton parities
in the SM and mirror sectors, also discussed in section 3.4.1, so that
$y_{ij}=0$. This case gives the $y_{i}=0$ limit of Fig. 4. Since there are
only direct contributions to the light neutrino masses, the two heavier mirror
neutrinos must have masses of order 10 GeV. They have radiative decays, via a
virtual loop of $W^{\prime}$ and $e^{\prime}$, to a mirror photon and the
lightest mirror neutrino, which is stable, with lifetimes of order $10^{9}$s.
These decays do not affect the BBN abundances or CMB data, but we must require
that $e^{\prime}$ dark matter is not diluted and the contribution of the
lightest $\nu^{\prime}$ to dark matter is sub-dominant. This will certainly be
the case if the freeze-in abundance of $\nu^{\prime}$ via virtual $W^{\prime}$
gives a yield $Y_{\nu^{\prime}}<10^{-10}$, which results if the reheat
temperature after inflation is less than about 1 TeV.

##  5 Signals

The mirror world, with $v^{\prime}$ determined to be near $10^{8}$ GeV by the
$e^{\prime}$ freeze-out abundance, has several observational signals that
depend on few free parameters. These signals do not depend on whether
$v^{\prime}$ is set by the Higgs Parity mechanism or by soft $Z_{2}$ breaking
Higgs interactions. Even though this scale is far above that of collider
energies – so that mirror particles are not directly accessible – a
combination of observations could provide strong evidence for the theory. We
first discuss signals from dark radiation in Sec. 5.1, and the rare kaon decay
$K\rightarrow\pi a$, in Sec. 5.2, which depend on $f_{a}$. Dark radiation is
also sensitive to $E/N$, the ratio of PQ anomalies. Next, we discuss signals
from self-interactions of $e^{\prime}$ in galactic halos in Sec. 5.3 and
gravity waves in Sec. 5.4. With $v^{\prime}$ fixed, these signals do not
depend on any unknown parameters of the theory, though there are uncertainties
in the computed signals.

###  5.1 Dark Radiation, $N_{\rm eff}$

There are two effects that lead to a change in the dark radiation, $\Delta
N_{\rm eff}$: 1) the decay of the axion around the neutrino decoupling era
heats up the photon bath relative to the neutrino bath which generates a
negative $\Delta N_{\rm eff}$, and 2) the relic mirror photons act as dark
radiation which generate a positive $\Delta N_{\rm eff}$. Quantitatively,
these effect are described by the effective number of neutrinos,

$\displaystyle N_{\rm
eff}=\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\frac{\rho_{\nu}+\rho_{\gamma^{\prime}}}{\rho_{\gamma}},$
(39)

where the expected Standard Model value is $N_{\rm eff}^{(\rm SM)}=3.044$
Bennett:2020zkv ; Froustey:2020mcq ; particle2020review . In previous work
Dunsky:2022uoq , we numerically computed the value of $N_{\rm eff}$ for
general heavy QCD axion theories with and without frozen-in dark photons. Here
we extend that work to incorporate the contribution to $N_{\rm eff}$ from
mirror photons frozen-out at the two-sector decoupling temperature of around
$\sim 20\mathchar 45\relax 30$ GeV. Fig. 5 shows the constraints in the
$(m_{a},f_{a})$ plane for $E/N=2/3$, where $E/N$ is the ratio of the electric
to color anomaly of the PQ quarks when these quarks possess down-quark like
quantum numbers.

Figure 5: Left: Contours of $N_{\rm eff}$ in the $m_{a}-f_{a}$ plane for
$E/N=2/3$ arising from axion decays and the relic mirror photon abundance. The
dark orange region is excluded at 95$\%$ confidence by Planck, while the light
orange shows the future reach of CMB-S4 experiment at 95$\%$ confidence. The
green region shows the parameter space where the mirror photons thermalize and
the freeze-in picture we use to compute $N_{\rm eff}$ breaks down. The red
dashed line shows the constraint (12) on $m_{a}f_{a}$ for $v^{\prime}=8\times
10^{7}$ GeV. Right: Same as left but zoomed in on the allowed region and
showing constraints from rare kaon decays in purple.

Figure 6: Same as Fig. 5, but for $E/N=0$. For $E/N=0$, the $a\rightarrow
2\gamma^{\prime}$ branching ratio is negligible so that $\Delta N_{\rm
eff}\simeq 0$ when the axion decays well before neutrino decoupling as shown
by the white region in the large $m_{a}$ and low $f_{a}$ region of the left
panel.

The blue contours of Fig. 5 indicate the value of $N_{\rm eff}$ in the
$(m_{a},f_{a})$ plane. The dark orange region shows the current $2\sigma$
limit on $N_{\rm eff}$ from Planck Planck:2018vyg  while the light orange
region shows the future reach of CMB-S4 CMB-S4:2016ple . To get the right
$e^{\prime}$ dark matter abundance, the axion mass must lie on the red strip
as given by Eq. (12) for $v^{\prime}=8\times 10^{7}$ GeV.

We see that most of the axion parameter space is either ruled out by Planck or
will be probed by CMB-S4. The intersection of the presently allowed region
(light orange) and the locus of points on which the axion can live (dashed red
line), powerfully constrains $m_{a}$ to be between $6.5-13.5$ MeV and $f_{a}$
between $(3-7)\times 10^{4}$ GeV. The small white strip in the $(m_{a},f_{a})$
plane where the $N_{\rm eff}$ signal is so small that CMB-S4 cannot probe may
still be able to be probed by rare kaon decays as discussed in the next
section.

The left panel of Fig. 6 shows the same $N_{\rm eff}$ as Fig. 5 but for
$E/N=0$. For this charge-to-color anomaly ratio, the axion coupling to mirror
photons vanishes. Consequently, the axion decay to mirror photons does not
occur and the only positive contribution to $N_{\rm eff}$ originates in the
$\gamma^{\prime}$ frozen-out at $T\sim 30$ GeV. This contribution is small and
for computational simplicity, we neglect it in computing Fig. 6. Compared to
the case for $E/N=2/3$, the region of large $m_{a}$ and small $f_{a}$ opens up
since in this region, the axion decays far before neutrino decoupling and only
into Standard Model particles.

In general, the larger $E/N$ is, the greater the axion branching ratio into
mirror photons and the larger the positive contribution to $N_{\rm eff}$
becomes. For $E/N\gtrsim 1$, $g_{\gamma^{\prime}}>g_{\gamma}$ and the region
of large $m_{a}$ and low $f_{a}$ becomes increasingly ruled out by too large
$N_{\rm eff}$ from $a\rightarrow 2\gamma^{\prime}$ decays. Consequently, we
focus on the case $E/N=2/3$ (PQ quarks with down-like charges) and $E/N=0$
(KSVZ-like axion). The GUT-motivated value of $E/N=8/3$ generates such a large
positive $\Delta N_{\rm eff}$ that Planck excludes the entire parameter space
for $v^{\prime}=8\times 10^{7}$ GeV.

###  5.2 Axion Signal in Kaon Decays

Because the axion mixes with the neutral pseudo-scalar mesons $\pi_{0}$,
$\eta$, and $\eta^{\prime}$, Standard-Model hadronic decays involving such
mesons necessarily include axion decays as well. For sufficiently low $f_{a}$,
the decay $K\rightarrow\pi a$, which arises because of $a-\pi_{0}$ mixing, can
strongly constrain $f_{a}$ from bounds on rare kaon decays. For example, in
the Standard Model, the decay $K^{+}$ to $\pi^{+}$ and missing energy is
strongly suppressed as it is a flavor changing neutral current process. The
experimental limit on the branching ratio
Br($K^{+}\rightarrow\pi^{+}\nu\bar{\nu})$ is approximately $1\times 10^{-10}$,
which can be translated to a similar bound on the branching ratio
Br($K^{+}\rightarrow\pi^{+}a)$ due to its experimentally similar signature of
$\pi^{+}$ and missing energy NA62:2021zjw ; Bauer:2021wjo ; Bauer:2021mvw ;
Goudzovski:2022vbt . For $m_{a}\lesssim 100$ MeV, this translates to a bound
on $f_{a}\gtrsim(5-7)\times 10^{4}$ GeV.

The right panel of Fig. 5 shows a zoomed-in region of the left panel but now
overlaid with the kaon bounds. The lighter purple region shows the parameter
in the $(m_{a},f_{a})$ plane where the axion generates a larger than observed
branching ratio for the decay of $K^{+}$ to $\pi^{+}$ and missing energy in
the NA62 experiment, according to Fig. 7 of Goudzovski:2022vbt . We caution
there is some uncertainty in the matrix element associated with the octet
operator that enhances the $K^{+}\rightarrow\pi^{+}a$ decay rate and show a
more conservative NA62 bound in the darker purple region that is weaker by a
factor of $\sim 30\%$ Co:2022aav , in accordance with the uncertainty of
${\cal O}(m_{K}^{4}/m_{\rho}^{4})$ in the matrix element computed at NLO in
Cirigliano:2011ny . Future constraints from the NA62 Run 2 and the KLEVER
experiment will improve the limit on $f_{a}$ by a factor of $2$ to
$2\sqrt{2}$, respectively Goudzovski:2022vbt . The dot-dashed purple lines
shows the conservative region enhanced by these factors. Fig. 5 demonstrates
that even with some uncertainty in Br$(K\rightarrow a\pi)$, the axion bounds
from rare kaon decays is complementary to the $N_{\rm eff}$ bounds, probing
regions of lower $f_{a}$ that are not constrained by dark radiation. In
particular, the conservative estimate of the bound from $K$ decays implies
that $N_{\rm eff}$ is significantly below the SM value; CMB-S4 or the NA62 Run
2 and KLEVER experiments will either discover a signal or exclude this minimal
mirror axion theory.

###  5.3 Signals from Self-interactions of $e^{\prime}$ Dark Matter

Mirror electrons and positrons form the dominant component of dark matter, and
have self interactions via mirror QED with
$\alpha^{\prime}(\mu=m_{e}^{\prime})\simeq 0.0076$. Galactic halos have a
triaxial structure and self-interactions reduce the anisotropy of the DM
velocity distribution, reducing the ellipticity of the halos. The measured
ellipticity of NGC720 has led to constraints on self-interactions Feng:2009mn
. Remarkably, for $m_{e^{\prime}}=225$ GeV, our predicted value for
$\alpha^{\prime}$ is right on the upper limit Agrawal:2016quu . As stressed in
Agrawal:2016quu , these constraints are not strict bounds. One caveat is that
the constraint results from calculations of the time scale for the halo
velocity distribution to become isotropic, whereas the measured ellipticity is
probing the mass density distribution, which could evolve on a longer time
scale. $N$-body simulations could lead to more precise bounds. Also, the
constraint relies on observations of a single galaxy, NGC720, which might have
special features, for example low ellipticity; a more robust result would
require analysis of further galaxies. Finally, isotropization effects from
large self-interactions could be masked by galaxy mergers, which increase
ellipticities.

As mirror QED is unbroken, the effects on large scale structure from mirror
electron self-interactions results from frequent small angle scattering,
rather than rare large angle scattering. Using the results of Agrawal:2016quu
, the momentum-transfer cross-section for our $e^{\prime}$ dark matter is

$\displaystyle\frac{\sigma_{T}}{m_{e^{\prime}}}\simeq
2\,\frac{\rm{cm}^{2}}{\rm{g}}\left(\frac{300\;\rm{km\;s}^{-1}}{v}\right)^{4},$
(40)

and is highly dependent on the speed $v$. Typical speeds in (dwarf galaxies,
galaxies, galaxy clusters) are $(30,300,1000)\rm{km\;s}^{-1}$. Numerical
N-body simulations, including only dark matter, show that on galactic scales a
velocity-independent $\sigma_{T}/m$ of $1\rm{g}/\rm{cm}^{2}$ leads to
significant departures from CDM simulations Fischer:2022rko . These departures
are seen for several observables, including the matter power spectrum, halo
and subhalo mass functions, the number of satellites, halo density and
circular velocity profiles, and halo shape distributions. In particular, for
the number of satellites, the density profiles and the shape distributions,
the departures from CDM are larger when frequent small-angle scattering
dominates over rare large-angle scattering. This suggests that future
numerical simulations, with $\sigma_{T}\propto 1/v^{4}$, would identify
signals for mirror electron dark matter that could be searched for in the
data.

###  5.4 Gravitational Waves from $SU(3)^{\prime}$ Phase Transition

For $v^{\prime}\gtrsim 10^{6}$ GeV, the mirror QCD scale is below the mass of
$u^{\prime}$, the lightest mirror quark. Consequently, for $v^{\prime}=8\times
10^{7}$ GeV, the mirror QCD phase transition is first order Yaffe:1982qf ;
Svetitsky:1982gs . The first-order phase transition creates vacuum bubbles
that expand and relativistically collide, generating a gravitational wave
spectrum Witten:1984rs . The gravitational waves arise from bubble collisions,
plasma turbulence, and plasma sound waves caprini2018cosmological , so that
the spectrum maybe written as

$\displaystyle\Omega_{\rm GW}h^{2}=\left(\Omega_{\rm GW}^{(\rm
col)}+\Omega_{\rm GW}^{(\rm turb)}+\Omega_{\rm GW}^{(\rm
sound)}\right)h^{2}\,.$ (41)

Two key input quantities that determine the amplitude of the spectrum are the
ratio of vacuum energy density to the background radiation density at bubble
nucleation,

$\displaystyle\alpha\equiv\frac{\rho_{\rm vac}}{\rho_{\rm rad}}$ (42)

and the ratio of vacuum energy density to the component of the plasma density
that strongly couples to the expanding bubble,

$\displaystyle\alpha^{\prime}\equiv\frac{\rho_{\rm vac}}{\rho_{\rm
coupled}}\,.$ (43)

These energy density ratios affect the bubble wall velocity
huber2008gravitational ,

$\displaystyle v_{b}\simeq\frac{\sqrt{\frac{1}{3}}+\sqrt{\alpha^{\prime
2}+\frac{2\alpha^{\prime}}{3}}}{1+\alpha^{\prime}}\,,$ (44)

the efficiency of the transfer of vacuum energy into bubble kinetic energy
huber2008gravitational

$\displaystyle\kappa\equiv\frac{\rho_{\rm kin}}{\rho_{\rm
vac}}\simeq\frac{1}{1+0.715\alpha^{\prime}}\left(0.715\alpha^{\prime}+\frac{4}{27}\sqrt{\frac{3\alpha^{\prime}}{2}}\right)\,$
(45)

and the fractional kinetic energy density of the bubbles to the total energy
density of the universe Caprini:2010xv

$\displaystyle\frac{\rho_{\rm kin}}{\rho_{\rm rad}+\rho_{\rm
vac}}=\kappa\frac{\alpha}{1+\alpha}\,.$ (46)

Assuming that the energy density of the mirror gluon gas at $T_{\rm
QCD^{\prime}}$ ($\rho_{g}^{\prime}\simeq 1.07T_{\rm QCD^{\prime}}^{4}$
Borsanyi:2012ve ) is comparable to the vacuum energy, then $\alpha\approx
0.03$. The precise value of $\alpha^{\prime}$ is not as important since the
spectrum is fairly insensitive to $\alpha^{\prime}$ as long as it is
$\mathcal{O}(1)$. For a fiducial value, we take $\alpha^{\prime}\simeq 0.343$
Morgante:2022zvc . Last, we assume that the efficiency factor is the same for
collisions, turbulence, and sound waves. Other numerical fits of the
efficiency factor, such as in Espinosa:2010hh , give very similar values to
Eq. (45) in the $\alpha^{\prime}\sim 1$ limit.

The gravitational-wave spectra due to bubble collisions huber2008gravitational
; Caprini:2010xv , turbulence Caprini:2009yp ; Caprini:2010xv , and sound
waves Hindmarsh:2015qta ; Hindmarsh:2017gnf ; Caprini:2019egz  are

$\displaystyle\Omega_{\rm GW}^{(\rm col)}h^{2}\simeq\;$
$\displaystyle\Omega_{\rm
rad}h^{2}\left(\frac{\kappa\alpha}{1+\alpha}\right)^{2}\left(\frac{H}{\beta}\right)^{2}S_{c}(f)$
(47) $\displaystyle\Omega_{\rm GW}^{(\rm turb)}h^{2}\simeq\;$
$\displaystyle\Omega_{\rm
rad}h^{2}\left(\frac{\kappa\alpha}{1+\alpha}\right)^{3/2}\left(\frac{H}{\beta}\right)S_{t}(f)$
(48) $\displaystyle\Omega_{\rm GW}^{(\rm sound)}h^{2}\simeq\;$
$\displaystyle\Omega_{\rm
rad}h^{2}\left(\frac{\kappa\alpha}{1+\alpha}\right)^{3/2}\left(\frac{H}{\beta}\right)^{2}S_{s}(f)$
(49)

where $\Omega_{\rm rad}h^{2}\simeq 4.16\times 10^{-5}$ is the fractional
density of radiation today and $\beta/H$ quantifies the duration of the phase
transition. Note that the scaling of the efficiency in $\Omega_{\rm GW}^{(\rm
sound)}$ changes for long duration sources which are not applicable in the
large $\beta/H$ limit that we take here Caprini:2019egz . The frequency
dependence and effect of the bubble wall velocity is encoded in the spectral
shape functions $S(f)$, which are given for bubble collisions, turbulence, and
sound waves by

$\displaystyle S_{c}$
$\displaystyle=\frac{2}{3\pi^{2}}\frac{v_{b}^{3}}{0.42+v_{b}^{2}}\frac{(f/f_{c})^{3}}{1/3+(f/f_{c})^{4}}$
(50) $\displaystyle S_{t}$
$\displaystyle=\frac{2}{\pi^{2}}v_{b}^{2}\frac{(f/f_{t})^{3}(1+f/f_{t})^{-11/3}}{(f/f_{t})(\beta/H)+v_{b}/4\pi^{2}}$
(51) $\displaystyle S_{s}$ $\displaystyle=2.4\times
10^{-2}\,v_{b}^{2}(f/f_{s})^{3}\left(\frac{7}{4+3(f/f_{s})^{2}}\right)^{7/2}\,,$
(52)

with reference frequencies

$\displaystyle f_{c}$ $\displaystyle=3.5\times 10^{-6}\,{\rm
Hz}\left(\frac{\beta}{H}\right)\left(\frac{T_{\rm QCD^{\prime}}}{35\,\rm
GeV}\right)\left(\frac{g_{*}}{100}\right)^{1/6}\left(\frac{1}{1.8-0.1v_{b}+v_{b}^{2}}\right)$
(53) $\displaystyle f_{t}$ $\displaystyle=8.5\times 10^{-6}\,{\rm
Hz}\left(\frac{\beta}{H}\right)\left(\frac{T_{\rm QCD^{\prime}}}{35\,\rm
GeV}\right)\left(\frac{g_{*}}{100}\right)^{1/6}\left(\frac{1}{v_{b}}\right)$
(54) $\displaystyle f_{s}$ $\displaystyle=3.0\times 10^{-6}\,{\rm
Hz}\left(\frac{\beta}{H}\right)\left(\frac{T_{\rm QCD^{\prime}}}{35\,\rm
GeV}\right)\left(\frac{g_{*}}{100}\right)^{1/6}\left(\frac{1}{v_{b}}\right)\,.$
(55)

In this work, the temperature of the phase transition, $T_{\rm QCD^{\prime}}$,
is determined by (34) with $v^{\prime}\simeq 8\times 10^{7}$ GeV in order for
the relic freeze-out abundance of $e^{\prime}$ to match the observed dark
matter abundance. As given by Eq. (11), this sets $\Lambda_{\rm
QCD^{\prime}}\simeq 27$ GeV, or equivalently, $T_{\rm QCD^{\prime}}\simeq 34$
GeV. The duration of the phase transition, $\beta^{-1}$, is the main unknown
parameter that enters into the gravitational wave spectrum. Typically, the
duration of the phase transition is expected to be some fraction of Hubble,
$H$, with $\beta/H$ between $100-10^{5}$. Here 100 is the maximal possible
value Hogan:1983ixn , while $10^{5}$ is a value suggested by a holographic
computation Morgante:2022zvc . Note the larger value of $\beta/H$ (i.e., a
shorter duration phase transition) gives rise to a smaller amplitude in
$\Omega_{\rm GW}h^{2}$ but also a higher peak frequency of $f_{c}$ or $f_{t}$.

$\Omega_{\rm GW}h^{2}$ is shown in Fig. 7 for $\Lambda_{\rm
QCD^{\prime}}\simeq 27$ GeV ($T_{\rm QCD^{\prime}}\simeq 34$ GeV) and a range
of potential values for $\beta/H$. We see that future gravitational wave
detectors like LISA, BBO, and DECIGO amaro2017laser ; harry2006laser ;
Kawamura:2020pcg ; Schmitz:2020syl  may be able to observe a gravitational
wave signal from the mirror QCD phase transition if $\beta/H$ is not too
large. Once $\beta/H$ is known more precisely from simulations and theory, the
amplitude and peak frequency of the gravitational wave signal will be
predicted which can provide strong evidence for the simple model of high
quality axions in this work.

Figure 7: Gravitational wave signal arising from the first-order mirror QCD
phase transition. The solid, dashed, dotted, and dot-dashed curves show
$\Omega_{\rm GW}h^{2}$ for different assumptions of $\beta/H$. For
$\beta/H\lesssim 10^{4}$, the gravitational wave signal may be detected by the
BBO and DECIGO experiments.

##  6 Increasing the Mirror Electroweak Scale

Figure 8: Generalization of Fig. 4 but for $v^{\prime}\geq 10^{8}$ GeV. The
solid, dashed, and dotted contours correspond to dilution factors of
$1,10^{2}$ and $10^{4}$ generated by the decay of $\nu^{\prime}$. For
$v^{\prime}=10^{8}$ GeV, no dilution is required and the allowed parameter
space must lie to the right of the solid orange $D=1$ contour, and for natural
$\nu$ masses, within the solid blue triangle. For
$v^{\prime}=10^{9}\,(10^{10})$ GeV, $D=10^{2}\,(10^{4})$ is required and at
least one neutrino mass must lie on the dashed (dotted) orange contour and for
natural $\nu$ masses, within the dashed (dotted) blue triangle. The blue
solid, dashed, and dotted seesaw contributions are artificially displaced
slightly to make the allowed natural triangular region for each $v^{\prime}$
easier to see. The dashed red contour shows the allowed
$(m_{\nu^{\prime}_{i}},y_{i})$ for the singlet model as discussed in Sec.
3.4.2.

In this paper, we have discussed the cosmology and phenomenological signals of
a heavy QCD axion arising from a mirror world. To recap, the logic proceeded
as follows: matching the freeze-out abundance of mirror electrons to the
observed dark matter abundance fixed $v^{\prime}\simeq 8\times 10^{7}$ GeV.
Knowing $v^{\prime}$, we determined the mirror QCD scale from Eq. (11), the
axion mass as a function of $f_{a}$ (12), the mirror neutrino masses (23), and
the gravity wave spectrum (Fig. 7). A question naturally arises: if we relax
the constraint that the mirror electron is set by pure freeze-out, can
$v^{\prime}$ then be increased above $10^{8}$ GeV? If so, how does this change
the observable signals?

First, with larger $v^{\prime}$, the parameter space where the neutrino masses
are not fined-tuned grows beyond the triangular region of Fig. 2. In general,
each $\nu^{\prime}_{i}$ (with $i=1,2,3$) leads to two contributions to the
light neutrino mass matrix. As shown in (22) there is a direct mass term
$(v/v^{\prime})^{2}m_{\nu^{\prime},i}$ for $\nu_{i}$ and a seesaw term
$y_{i}^{2}v^{2}/m_{\nu^{\prime},i}$ for the active neutrino $\tilde{\nu}_{i}$,
defined by (21). Each term can naturally be at most $\sim 0.05$ eV (that is,
without any fine-tuned cancellation). The horizontal ‘Direct’ branch contours
in Fig. 8 indicate where the direct mass term is 0.05 eV for
$v^{\prime}=10^{8}$ GeV (solid), $10^{9}$ GeV (dashed), and $10^{10}$ GeV
(dotted), which grows quadratically with $v^{\prime}$. Similarly, the diagonal
‘Seesaw’ branch contours in Fig. 8 indicate where the seesaw mass term is 0.05
eV for the same three values of $v^{\prime}$. Note it is possible that one
light neutrino is much less than $0.05$ eV, in which case the allowed
parameter space lies within the triangle bounded by the Direct and Seesaw
contours. However, at least one light neutrino should have a mass comparable
to $0.05$ eV in which case at least one $m_{\nu^{\prime},i}$ lives
approximately on the Direct or Seesaw contour.

Second, when $v^{\prime}$ increases, $m_{e}^{\prime}$ increases so that
matching the relic yield of $e^{\prime}$ to the observable dark matter energy
density requires that $e^{\prime}$ possesses a yield lower than that from
freeze-out. One way to realize this is to dilute the $e^{\prime}$ abundance
after it freezes-out, which can be accomplished naturally by long-lived
$\nu^{\prime}$. The orange contours of Fig. 8 show different values of the
dilution factor (see Eq. (38)) in the $y_{\nu}-m_{\nu}^{\prime}$ plane. To the
right of the solid orange ($D=1$) line, $y_{\nu}$ is sufficiently large that
$\nu^{\prime}$ decays early enough to avoid generating any dilution. This is
the region for $v^{\prime}=10^{8}$ GeV, where no dilution of $e^{\prime}$ is
necessary, and is shown enlarged in Fig. 4. To the left of this line,
$e^{\prime}$ dilution from $\nu^{\prime}$ decay occurs. Since the frozen-out
density of $e^{\prime}$ scales as $\Omega_{\rm e^{\prime}}\propto
m_{e^{\prime}}^{2}\propto v^{\prime 2}$, and since no dilution of $e^{\prime}$
is required for $v^{\prime}=10^{8}$ GeV, it follows that the required
$e^{\prime}$ dilution as a function of $v^{\prime}$ is
$D\approx(v^{\prime}/10^{8}\,{\rm GeV})^{2}$. The dashed and dotted orange
contours show the parameter space where $D=10^{2}$ and $10^{4}$, respectively.
To achieve the correct dilution for $v^{\prime}=10^{9}$ GeV and $10^{10}$ GeV
requires that one of the $\nu^{\prime}_{i}$ lives on the dashed or dotted
orange contours, respectively, and none live to the left.

This case of $e^{\prime}$ dark matter in the mirror world, via freeze-out and
dilution from $\nu^{\prime}$ decays, was considered in Dunsky:2019upk . In
that paper there was no axion, so $y_{D}$ of (8) was absent and $v^{\prime}$
was determined by the Higgs Parity mechanism, with no soft breaking of the
$Z_{2}$ and with only SM interactions responsible for the radiative generation
of the SM quartic. Fig. 10 of Dunsky:2019upk  then shows that $e^{\prime}$
dark matter requires
$v^{\prime}\mathop{}_{\textstyle\sim}^{\textstyle<}3\times 10^{10}$ GeV. Such
$v^{\prime}$ do not result from the central experimental values of
$m_{t}=(172.56\pm 0.48)$ GeV and $\alpha_{s}(M_{Z})=0.1179\pm 0.0009$. For
example, for central $\alpha_{s}$ they require $m_{t}$ to be at least
2$\sigma$ high; alternatively if $m_{t}$ is 1$\sigma$ high then $\alpha_{s}$
must be more than 1$\sigma$ low. More precise measurements of $m_{t}$ and
$\alpha_{s}$ will determine if this scheme of Dunsky:2019upk  is allowed.

Third, to get the necessary dilution, $\nu^{\prime}$ cannot decay after BBN
without distorting heavy element abundances. The green region of Fig. 8 shows
the parameter space where $\nu^{\prime}$ decays at temperatures below $T_{\rm
BBN}=4$ MeV. Fig. 8 demonstrates that the necessary dilution for $e^{\prime}$
dark matter, $D=10^{2}$ for $v^{\prime}=10^{9}$ GeV ($D=10^{4}$ for
$v^{\prime}=10^{10}$ GeV), requires one $\nu^{\prime}_{i}$ to live on the
short segment of the orange dashed (dotted) line in the unshaded region. The
length of this line segment is constrained by avoiding fine-tuning in the
neutrino sector and requiring $\nu^{\prime}$ to decay before BBN. This
$\nu^{\prime}_{i}$, with a mass of order $10^{3}$ GeV ($10^{5}$ GeV), gives a
relevant direct mass contribution to the corresponding light neutrino. The
other two $\nu^{\prime}_{i}$ can lie anywhere to the right of the orange
dashed (dotted) line (although still within the blue dashed or dotted lines
for naturalness), and at least one of them must give another relevant
contribution to the light neutrino masses.

Fourth, we emphasize that while allowing $\nu^{\prime}$ to dilute the
$e^{\prime}$ abundance makes $v^{\prime}>10^{8}$ GeV feasible, it does so at
the cost of diminishing some cosmological signals. For example, the
gravitational wave spectrum shown in Fig. 7 is reduced by a factor of
$D^{4/3}$, which can weaken the signal below the detection threshold of LISA,
BBO and DECIGO except for the largest $\beta/H$. Similarly, the dilution
generated by $\nu^{\prime}$, which is generated close to BBN, erases most of
the $\Delta N_{\rm eff}$ generated by axion decays or residual dark photon
radiation shown in Fig. 5.

Fifth, increasing $v^{\prime}$ increases $\Lambda_{\rm QCD}^{\prime}$ so that
the axion mass is also heavier according to Eq. (12). The three dashed red
lines of Fig. 9 shows the $m_{a}-f_{a}$ correlation for
$v^{\prime}=(10^{8},10^{9},10^{10})$ GeV. The orange regions shows the
constraints from SN 1987A summarized in Kelly:2020dda  with the darker region
Chang:2018rso  more conservative for low $f_{a}$ than the lighter region
Ertas:2020xcc . The purple region shows the bound from searches for the rare
kaon decay $K^{+}\rightarrow\pi^{+}a$.

Sixth, since the constraints from $\Delta N_{\rm eff}$ become less powerful
with increasing $v^{\prime}$, larger values of the PQ anomaly ratio, $E/N$,
are allowed. Thus, unlike at $v^{\prime}=10^{8}$ GeV, it is possible for the
SM to be unified, for example into an $SU(5)$ theory with $E/N=8/3$. This
introduces two more observational signals that have event rates with very
different dependencies on the unified scale $v_{G}$. The proton decay event
rate scales like $1/v_{G}^{4}$, while the event rate for the direct detection
of $e^{\prime}$ dark matter from kinetic mixing scales as a positive of
$v_{G}$ (because $\epsilon$ is generated by inserting unified symmetry
breaking into a higher dimension operator). Results for these signals are
shown in Figure 5 of Dunsky:2019upk . With $\epsilon$ arising from a dimension
8 operator, proton decay will be discovered at the Hyper-Kamiokande experiment
and dark matter at the LZ experiment.

Finally, what if the neutrinos are Dirac, or a lepton parity sets $y_{i}=0$,
as discussed in section 3.4.1? In the case of Dirac neutrinos, $\nu^{\prime}$
become $\nu_{R}$, degenerate with $\nu_{L}$. Since there are no heavy
$\nu^{\prime}$ states, it is not possible to dilute the mirror electron
abundance by $\nu^{\prime}$ decay; dark matter requires $v^{\prime}=10^{8}$
GeV. For the case with $y_{i}=0$ from lepton parity, the two heavier
$\nu^{\prime}$ decay radiatively to the mirror photon and the lightest
$\nu^{\prime}$, which is stable. To give $e^{\prime}$ dilution and successful
nucleosynthesis, these decays must occur before 4 MeV, requiring much larger
$m^{\prime}_{\nu}$ and therefore larger $v^{\prime}$ of at least $10^{10}$
GeV. Any further increase in $v^{\prime}$ makes $\nu^{\prime}$ sufficiently
heavy that the mirror beta decay $\nu^{\prime}\rightarrow
e^{\prime}u^{\prime}\bar{d}^{\prime}$ opens up. The resulting rapid decays of
$\nu^{\prime}$ imply that, at these large values of $v^{\prime}$, the freeze-
out abundance of $e^{\prime}$ cannot be sufficiently diluted. Hence, this
scheme works, if at all, only for a narrow range of $v^{\prime}$ near $3\times
10^{10}$ GeV.

Figure 9:  For $v^{\prime}>10^{8}$ GeV, dilution of $e^{\prime}$ dark matter
by $\nu^{\prime}$ decays is required, significantly altering the $N_{\rm eff}$
bounds of Fig. 5. Nevertheless, supernova (orange) and rare kaon decay
(purple) constraints remain. The dashed red contours show the allowed axion
mass for $v^{\prime}=10^{8},10^{9}$ and $10^{10}$ GeV.

##  7 Summary

In this paper, we studied a mirror-world scenario with a PQ symmetry, where
the strong CP problem is solved by an axion. The mirror QCD scale is much
larger than the SM QCD scale because of the large mirror electroweak scale,
$v^{\prime}\gg v$. The axion mass is dominantly given by mirror QCD dynamics
and is larger than the standard case. The axion decay constant may be much
below $10^{9}$ GeV. Because of the larger axion mass and the smaller decay
constant, the PQ symmetry may be easily understood as an accidental symmetry.

The mirror electron is absolutely stable because of mirror electromagnetic
charge conservation and is a natural dark matter candidate. Self-interactions
from mirror QED affect large scale structure, such as galactic density
profiles and shape distributions, and the number of galactic satellites. The
observed dark matter abundance is explained by freeze out of the annihilation
of mirror electrons if the mirror electroweak scale $v^{\prime}$ is around
$10^{8}$ GeV. The Higgs Parity mechanism can set such a value of $v^{\prime}$
if the coupling of the Higgs with the KSVZ fermions is $O(1)$. The mirror QCD
scale is uniquely predicted, and the axion mass is predicted as a function of
the decay constant, as shown in Eq. (12). The mirror QCD phase transition is
first order and produces primordial gravitational waves. Since the unique
parameter of the phase transition, namely, the mirror QCD scale, is fixed the
spectrum of the gravitational waves can in principle be uniquely predicted.
Practically, the mirror QCD phase transition is induced by strong dynamics so
we cannot precisely compute the spectrum. We show the spectrum as a function
of parameters that characterize the phase transition in Fig. 7.

The heavy axion plays an important role in the cosmology of the mirror world.
Without the axion, the mirror QCD entropy would be transferred into mirror
photons, producing too much dark radiation. The axion instead maintains the
thermal equilibrium of the SM and mirror sectors and the mirror QCD entropy is
shared with the SM bath. Still, the radiation energy density of the universe
is modified by the presence of mirror photons and by the decay of the axion
into SM and mirror photons after neutrinos decouple. The constraint on the
parameter space is shown in Figs. 5 and 6. Future observations of the cosmic
microwave background and rare kaon decays will probe all of the remaining
parameter space.

The mirror neutrinos also affect the cosmological evolution. If they obtain
Majorana masses from dimension-five operators, they are massive and may decay
after they dominate the universe, diluting mirror electron dark matter and
gravitational waves. This is avoided if the observed neutrino masses are
explained by comparable contributions from the dimension-five operators and
the see-saw mechanism from the mirror neutrinos, as shown in Fig. 8. Such
comparable contributions are a prediction of the minimal UV completion of the
dimension-five operators. If significant entropy production does occur,
$v^{\prime}>10^{8}$ GeV becomes also viable. If lepton symmetry is preserved,
the mirror neutrinos obtain Dirac masses with the SM neutrinos. They are
nearly massless and behave as dark radiation.

###### Acknowledgements.

DD is supported by the James Arthur Postdoctoral Fellowship. The work of LH
was supported by the Director, Office of Science, Office of High Energy
Physics of the U.S. Department of Energy under the Contract No. DE-
AC02-05CH11231 and by the NSF grants PHY-1915314 and PHY-2210390.

##  Appendix A Accidental PQ symmetry

In this appendix, we discuss a UV completion where the PQ symmetry
accidentally arises as a result of other exact symmetry. We impose
$Z_{2}\times Z_{3}$ symmetry and introduce a complex scalar field $P$ and
three pairs of quarks $D_{i}$, $\bar{D}_{i}$ and their mirror partners
$D^{\prime}_{i}$, $\bar{D}^{\prime}_{i}$. Under the $Z_{2}\times Z_{3}$
symmetry, they transform as

$\displaystyle Z_{2}$ $\displaystyle:P\rightarrow-P,~{}~{}D_{i}\rightarrow-
D_{i}^{\prime},~{}~{}\bar{D}_{i}\rightarrow\bar{D}_{i}^{\prime}$
$\displaystyle Z_{3}$ $\displaystyle:P\rightarrow e^{i2\pi/3}\times
P,~{}~{}D_{i}\rightarrow e^{-i2\pi/3}\times
D_{i},~{}~{}D^{\prime}_{i}\rightarrow e^{-i2\pi/3}\times D^{\prime}_{i}$ (56)

$Z_{3}\times Z_{2}$ does not have QCD or gravitational anomaly. Instead of
adding three pairs of fermions in the fundamental representation of $SU(3)$,
we may consider higher representations of $SU(3)\times SU(2)$, such as $({\bf
8},{\bf 1})$ and $({\bf 3},{\bf 3})+({\bf\bar{3}},{\bf 3})$. For the latter
case, $E/N$ becomes large so we need dilution after the mirror photon
decouples or more fermions to reduce the electromagnetic anomaly coefficient.

The renormalizable couplings among them allowed by the symmetry are

$\displaystyle\xi_{i}P\left(D_{i}\bar{D}_{i}+D^{\prime}_{i}\bar{D}^{\prime}_{i}\right)-\lambda\left(|P|^{2}-\frac{9}{2}f_{a}^{2}\right)^{2},$
(57)

which preserve an accidental $U(1)_{\rm PQ}$ symmetry. Here the factor of $9$
comes from the square of the domain wall number. The leading explicit PQ-
breaking term is $P^{6}/M^{2}$, where $M$ is the cutoff scale. As long as

$\displaystyle m_{a}\gtrsim 10^{5}\frac{f_{a}^{2}}{M}\sim 0.1~{}{\rm
MeV}\left(\frac{f_{a}}{10^{5}~{}{\rm GeV}}\right)^{2}\frac{M_{\rm pl}}{M},$
(58)

the strong CP problem is solved.

It is intriguing to identify the $Z_{2}$ symmetry above with the mirror
symmetry. However, for that case, the leading explicit PQ breaking term is
$P^{3}(|H|^{2}-|H^{\prime}|^{2})/M$, and the strong CP problem cannot be
solved. This can be avoided by embedding the theory into solutions to the
electroweak hierarchy problem where $|H|^{2}$ and $|H^{\prime}|^{2}$ are
controlled by additional structure, such as supersymmetry or compositeness of
the Higgs. Another possibility is to generalize the model into the ones with
$Z_{2n+1}\times Z_{2}$ symmetry with $n>1$ and forbid explicit PQ breaking up
to $P^{(2n+1)}(|H|^{2}-|H^{\prime}|^{2})$. This can be achieved by introducing
$2n+1$ pairs of fermions.

## References

  * (1) T. D. Lee and C.-N. Yang, _Question of Parity Conservation in Weak Interactions_ , _Phys. Rev._ 104 (1956) 254.
  * (2) I. Y. Kobzarev, L. B. Okun and I. Y. Pomeranchuk, _On the possibility of experimental observation of mirror particles_ , _Sov. J. Nucl. Phys._ 3 (1966) 837.
  * (3) L. B. Okun, _Mirror particles and mirror matter: 50 years of speculations and search_ , _Phys. Usp._ 50 (2007) 380 [hep-ph/0606202].
  * (4) R. Foot, _Mirror dark matter: Cosmology, galaxy structure and direct detection_ , _Int. J. Mod. Phys. A_ 29 (2014) 1430013 [1401.3965].
  * (5) B. Holdom, _Two U(1)’s and Epsilon Charge Shifts_ , _Phys. Lett. B_ 166 (1986) 196.
  * (6) H. Goldberg and L. J. Hall, _A New Candidate for Dark Matter_ , _Phys. Lett. B_ 174 (1986) 151.
  * (7) F.-Y. Cyr-Racine and K. Sigurdson, _Cosmology of atomic dark matter_ , _Phys. Rev. D_ 87 (2013) 103515 [1209.5752].
  * (8) S. Bansal, J. Barron, D. Curtin and Y. Tsai, _Precision Cosmological Constraints on Atomic Dark Matter_ , 2212.02487.
  * (9) R. Foot, H. Lew and R. R. Volkas, _Unbroken versus broken mirror world: A Tale of two vacua_ , _JHEP_ 07 (2000) 032 [hep-ph/0006027].
  * (10) R. Barbieri, T. Gregoire and L. J. Hall, _Mirror world at the large hadron collider_ , hep-ph/0509242.
  * (11) L. J. Hall and K. Harigaya, _Implications of Higgs Discovery for the Strong CP Problem and Unification_ , _JHEP_ 10 (2018) 130 [1803.08119].
  * (12) D. Dunsky, L. J. Hall and K. Harigaya, _Dark Matter, Dark Radiation and Gravitational Waves from Mirror Higgs Parity_ , _JHEP_ 02 (2020) 078 [1908.02756].
  * (13) R. D. Peccei and H. R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * (14) R. D. Peccei and H. R. Quinn, _Constraints Imposed by CP Conservation in the Presence of Instantons_ , _Phys. Rev. D_ 16 (1977) 1791.
  * (15) Z. Berezhiani, L. Gianfagna and M. Giannotti, _Strong CP problem and mirror world: The Weinberg-Wilczek axion revisited_ , _Phys. Lett. B_ 500 (2001) 286 [hep-ph/0009290].
  * (16) V. A. Rubakov, _Grand unification and heavy axion_ , _JETP Lett._ 65 (1997) 621 [hep-ph/9703409].
  * (17) H. Fukuda, K. Harigaya, M. Ibe and T. T. Yanagida, _Model of visible QCD axion_ , _Phys. Rev. D_ 92 (2015) 015021 [1504.06084].
  * (18) A. Hook, S. Kumar, Z. Liu and R. Sundrum, _High Quality QCD Axion and the LHC_ , _Phys. Rev. Lett._ 124 (2020) 221801 [1911.12364].
  * (19) Z. Chacko, H.-S. Goh and R. Harnik, _The Twin Higgs: Natural electroweak breaking from mirror symmetry_ , _Phys. Rev. Lett._ 96 (2006) 231802 [hep-ph/0506256].
  * (20) R. Barbieri, L. J. Hall and K. Harigaya, _Minimal Mirror Twin Higgs_ , _JHEP_ 11 (2016) 172 [1609.05589].
  * (21) H. Fukuda, M. Ibe and T. T. Yanagida, _Dark Matter Candidates in a Visible Heavy QCD Axion Model_ , _Phys. Rev. D_ 95 (2017) 095017 [1702.00227].
  * (22) R. Barbieri, L. J. Hall and K. Harigaya, _Effective Theory of Flavor for Minimal Mirror Twin Higgs_ , _JHEP_ 10 (2017) 015 [1706.05548].
  * (23) D. Buttazzo, G. Degrassi, P. P. Giardino, G. F. Giudice, F. Sala, A. Salvio et al., _Investigating the near-criticality of the Higgs boson_ , _JHEP_ 12 (2013) 089 [1307.3536].
  * (24) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * (25) W. E. Caswell, _Asymptotic Behavior of Nonabelian Gauge Theories to Two Loop Order_ , _Phys. Rev. Lett._ 33 (1974) 244.
  * (26) D. R. T. Jones, _Two Loop Diagrams in Yang-Mills Theory_ , _Nucl. Phys. B_ 75 (1974) 531.
  * (27) M. E. Machacek and M. T. Vaughn, _Two Loop Renormalization Group Equations in a General Quantum Field Theory. 1. Wave Function Renormalization_ , _Nucl. Phys. B_ 222 (1983) 83.
  * (28) S. Durr, Z. Fodor, C. Hoelbling and T. Kurth, _Precision study of the SU(3) topological susceptibility in the continuum_ , _JHEP_ 04 (2007) 055 [hep-lat/0612021].
  * (29) Flavour Lattice Averaging Group (FLAG) collaboration, _FLAG Review 2021_ , _Eur. Phys. J. C_ 82 (2022) 869 [2111.09849].
  * (30) M. Kamionkowski and J. March-Russell, _Planck scale physics and the Peccei-Quinn mechanism_ , _Phys. Lett. B_ 282 (1992) 137 [hep-th/9202003].
  * (31) S. M. Barr and D. Seckel, _Planck scale corrections to axion models_ , _Phys. Rev. D_ 46 (1992) 539.
  * (32) S. Ghigna, M. Lusignoli and M. Roncadelli, _Instability of the invisible axion_ , _Phys. Lett. B_ 283 (1992) 278.
  * (33) 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].
  * (34) D. Aloni, Y. Soreq and M. Williams, _Coupling QCD-Scale Axionlike Particles to Gluons_ , _Phys. Rev. Lett._ 123 (2019) 031803 [1811.03474].
  * (35) L. J. Hall and K. Harigaya, _Higgs Parity Grand Unification_ , _JHEP_ 11 (2019) 033 [1905.12722].
  * (36) D. Dunsky, L. J. Hall and K. Harigaya, _CHAMP Cosmic Rays_ , _JCAP_ 07 (2019) 015 [1812.11116].
  * (37) A. Salvio, A. Strumia and W. Xue, _Thermal axion production_ , _JCAP_ 01 (2014) 011 [1310.6982].
  * (38) F. D’Eramo, F. Hajkarim and S. Yun, _Thermal axion production at low temperatures: a smooth treatment of the QCD phase transition_ , 2108.04259.
  * (39) F. D’Eramo, F. Hajkarim and S. Yun, _Thermal QCD Axions across Thresholds_ , 2108.05371.
  * (40) S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, _Precision SU(3) lattice thermodynamics for a large temperature range_ , _JHEP_ 07 (2012) 056 [1204.6184].
  * (41) M. Bolz, A. Brandenburg and W. Buchmuller, _Thermal production of gravitinos_ , _Nucl. Phys. B_ 606 (2001) 518 [hep-ph/0012052].
  * (42) D. Cadamuro and J. Redondo, _Cosmological bounds on pseudo Nambu-Goldstone bosons_ , _JCAP_ 02 (2012) 032 [1110.2895].
  * (43) J. Kang, M. A. Luty and S. Nasri, _The Relic abundance of long-lived heavy colored particles_ , _JHEP_ 09 (2008) 086 [hep-ph/0611322].
  * (44) K. Harigaya, M. Ibe, K. Kaneta, W. Nakano and M. Suzuki, _Thermal Relic Dark Matter Beyond the Unitarity Limit_ , _JHEP_ 08 (2016) 151 [1606.00159].
  * (45) V. De Luca, A. Mitridate, M. Redi, J. Smirnov and A. Strumia, _Colored Dark Matter_ , _Phys. Rev. D_ 97 (2018) 115024 [1801.01135].
  * (46) E. W. Kolb and M. S. Turner, _The Early Universe_ , vol. 69. 1990, 10.1201/9780429492860.
  * (47) 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].
  * (48) J. J. Bennett, G. Buldgen, P. F. De Salas, M. Drewes, S. Gariazzo, S. Pastor et al., _Towards a precision calculation of $N_{\rm eff}$ in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED_, _JCAP_ 04 (2021) 073 [2012.02726].
  * (49) J. Froustey, C. Pitrou and M. C. Volpe, _Neutrino decoupling including flavour oscillations and primordial nucleosynthesis_ , _JCAP_ 12 (2020) 015 [2008.01074].
  * (50) P. D. Group, P. Zyla, R. Barnett, J. Beringer, O. Dahl, D. Dwyer et al., _Review of particle physics_ , _Progress of Theoretical and Experimental Physics_ 2020 (2020) 083C01.
  * (51) D. I. Dunsky, L. J. Hall and K. Harigaya, _Dark Radiation Constraints on Heavy QCD Axions_ , 2205.11540.
  * (52) Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * (53) CMB-S4 collaboration, _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * (54) NA62 collaboration, _Measurement of the very rare K +→${\pi}^{+}\nu\overline{\nu}$ decay_, _JHEP_ 06 (2021) 093 [2103.15389].
  * (55) M. Bauer, M. Neubert, S. Renner, M. Schnubel and A. Thamm, _Consistent Treatment of Axions in the Weak Chiral Lagrangian_ , _Phys. Rev. Lett._ 127 (2021) 081803 [2102.13112].
  * (56) M. Bauer, M. Neubert, S. Renner, M. Schnubel and A. Thamm, _Flavor probes of axion-like particles_ , _JHEP_ 09 (2022) 056 [2110.10698].
  * (57) E. Goudzovski et al., _New Physics Searches at Kaon and Hyperon Factories_ , 2201.07805.
  * (58) R. T. Co, T. Gherghetta and K. Harigaya, _Axiogenesis with a heavy QCD axion_ , _JHEP_ 10 (2022) 121 [2206.00678].
  * (59) V. Cirigliano, G. Ecker, H. Neufeld, A. Pich and J. Portoles, _Kaon Decays in the Standard Model_ , _Rev. Mod. Phys._ 84 (2012) 399 [1107.6001].
  * (60) J. L. Feng, M. Kaplinghat, H. Tu and H.-B. Yu, _Hidden Charged Dark Matter_ , _JCAP_ 07 (2009) 004 [0905.3039].
  * (61) P. Agrawal, F.-Y. Cyr-Racine, L. Randall and J. Scholtz, _Make Dark Matter Charged Again_ , _JCAP_ 05 (2017) 022 [1610.04611].
  * (62) M. S. Fischer, M. Brüggen, K. Schmidt-Hoberg, K. Dolag, F. Kahlhoefer, A. Ragagnin et al., _Cosmological simulations with rare and frequent dark matter self-interactions_ , _Mon. Not. Roy. Astron. Soc._ 516 (2022) 1923 [2205.02243].
  * (63) L. G. Yaffe and B. Svetitsky, _First Order Phase Transition in the SU(3) Gauge Theory at Finite Temperature_ , _Phys. Rev. D_ 26 (1982) 963.
  * (64) B. Svetitsky and L. G. Yaffe, _Critical Behavior at Finite Temperature Confinement Transitions_ , _Nucl. Phys. B_ 210 (1982) 423.
  * (65) E. Witten, _Cosmic Separation of Phases_ , _Phys. Rev. D_ 30 (1984) 272.
  * (66) C. Caprini and D. G. Figueroa, _Cosmological backgrounds of gravitational waves_ , _Classical and Quantum Gravity_ 35 (2018) 163001.
  * (67) S. J. Huber and T. Konstandin, _Gravitational wave production by collisions: more bubbles_ , _Journal of Cosmology and Astroparticle Physics_ 2008 (2008) 022.
  * (68) C. Caprini, R. Durrer and X. Siemens, _Detection of gravitational waves from the QCD phase transition with pulsar timing arrays_ , _Phys. Rev. D_ 82 (2010) 063511 [1007.1218].
  * (69) E. Morgante, N. Ramberg and P. Schwaller, _Echo of the Dark: Gravitational Waves from Dark SU(3) Yang-Mills Theory_ , 2210.11821.
  * (70) J. R. Espinosa, T. Konstandin, J. M. No and G. Servant, _Energy Budget of Cosmological First-order Phase Transitions_ , _JCAP_ 06 (2010) 028 [1004.4187].
  * (71) C. Caprini, R. Durrer and G. Servant, _The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition_ , _JCAP_ 12 (2009) 024 [0909.0622].
  * (72) M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, _Numerical simulations of acoustically generated gravitational waves at a first order phase transition_ , _Phys. Rev. D_ 92 (2015) 123009 [1504.03291].
  * (73) M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, _Shape of the acoustic gravitational wave power spectrum from a first order phase transition_ , _Phys. Rev. D_ 96 (2017) 103520 [1704.05871].
  * (74) C. Caprini et al., _Detecting gravitational waves from cosmological phase transitions with LISA: an update_ , _JCAP_ 03 (2020) 024 [1910.13125].
  * (75) C. J. Hogan, _NUCLEATION OF COSMOLOGICAL PHASE TRANSITIONS_ , _Phys. Lett. B_ 133 (1983) 172.
  * (76) P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender et al., _Laser interferometer space antenna_ , _arXiv preprint arXiv:1702.00786_ (2017) .
  * (77) G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner and E. S. Phinney, _Laser interferometry for the big bang observer_ , _Classical and Quantum Gravity_ 23 (2006) 4887.
  * (78) S. Kawamura et al., _Current status of space gravitational wave antenna DECIGO and B-DECIGO_ , _PTEP_ 2021 (2021) 05A105 [2006.13545].
  * (79) K. Schmitz, _New Sensitivity Curves for Gravitational-Wave Signals from Cosmological Phase Transitions_ , _JHEP_ 01 (2021) 097 [2002.04615].
  * (80) K. J. Kelly, S. Kumar and Z. Liu, _Heavy axion opportunities at the DUNE near detector_ , _Phys. Rev. D_ 103 (2021) 095002 [2011.05995].
  * (81) J. H. Chang, R. Essig and S. D. McDermott, _Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle_ , _JHEP_ 09 (2018) 051 [1803.00993].
  * (82) F. Ertas and F. Kahlhoefer, _On the interplay between astrophysical and laboratory probes of MeV-scale axion-like particles_ , _JHEP_ 07 (2020) 050 [2004.01193].