A Heavy QCD Axion and the Mirror World

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

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

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

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

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

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.

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

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

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

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

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.

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

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 .

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

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

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

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.

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

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

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

The decay rate into $2a$ is given by

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

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

The light neutrino mass matrix is therefore

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

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}$

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

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$

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

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

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

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.¹¹1The 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:

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,

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}$,

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$. ²²2 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

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.

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

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}$),

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.

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.