††institutetext: School of Particles and Accelerators, Institute Research in
Fundamental Sciences (IPM), P. O. Box 19395-5531, Tehran, Iran

# Filtered asymmetric dark matter during the Peccei-Quinn phase transition

M. Ahmadvand<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a bubble filtering-out mechanism for an asymmetric
dark matter scenario during the Peccei-Quinn (PQ) phase transition. Based on a
QCD axion model, extended by extra chiral neutrinos, we show that the PQ phase
transition can be first order in the parameter space of the model and
regarding the PQ symmetry breaking scale, the mechanism can generate PeV-scale
heavy neutrinos as a dark matter candidate. Considering a CP-violating source,
during the phase transition, discriminating between the neutrino and
antineutrino number density, we find the observed dark matter relic abundance,
such that the setup can be applied to the first order phase transition with
different strengths. We then calculate effective couplings of the QCD axion
addressing the strong CP problem within the model. We also study the energy
density spectrum of gravitational waves generated from the first order phase
transition and show that the signals can be detected by future ground-based
detectors such as Einstein Telescope. In particular, for a visible heavy axion
case of the model, it is shown that gravitational waves can be probed by
DECIGO and BBO interferometers. Furthermore, we discuss the dark matter-
standard model neutrino annihilation process as a source for the creation of
PeV-scale neutrinos.

###### Keywords:

Dark matter, Right-handed neutrinos, Peccei-Quinn phase transition, Axion

##  1 Introduction

Astrophysical and cosmological observations including measurements of galactic
rotation curves, gravitational lensings, and anisotropies in the Cosmic
Microwave Background (CMB) support the existence of Dark Matter (DM) forming
around 26% of the Universe energy density Bertone:2016nfn . However, the
nature of DM in fundamental physics is still unknown.

To explain the relic abundance of the cold DM from the early Universe, various
mechanisms and candidates have been suggested. Due to the coincidence between
the weak scale and massive particles with a mass range from GeV to TeV, and
also relying on the cosmological expansion and thermal freeze-out mechanism,
Weakly Interacting Massive Particle (WIMP) paradigm has been an attractive
scenario Bertone:2004pz , though direct detection searches disfavor some of
these scenarios and lead us to consider other possibilities such as non-
thermal DM models Allahverdi:2010rh  and scenarios in which DM is produced in
the decay of thermally-decoupled heavy particles Ahmadvand:2020izy ;
Kolb:1998ki . Another remarkable fact about the DM energy density is its
closeness to the abundance of asymmetric baryonic matter. This may imply a
common asymmetric origin for dark and visible matter and is the motivation for
asymmetric DM scenarios Ahmadvand:2020izy ; Kaplan:2009ag .

In this paper, we use the so-called filtered DM mechanism through which DM
dynamically acquires a mass and its number density is abruptly frozen out
Baker:2019ndr ; Chway:2019kft . Based on first order Phase Transitions (PTs)
and DM interactions with bubbles, this mechanism provides a framework to evade
the Griest-Kamionkowski (GK) bound on the mass of thermally-produced DM and a
possibility to produce DM masses above $100\,\mathrm{TeV}$.

On the other hand, at these high energy scales, one of the well-motivated
cosmological PTs could have taken place. In the early Universe at temperatures
around $100\,\mathrm{PeV}$, the Peccei-Quinn (PQ) symmetry
$\mathrm{U(1)}_{\mathrm{PQ}}$ is spontaneously broken according to QCD axion
models, resolving the strong CP problem Peccei:1977hh ; Peccei:2006as
.111Considering astrophysical constraints, the symmetry breaking scale is
bounded between few $10^{8}\,\mathrm{GeV}$ and few $10^{17}\,\mathrm{GeV}$
Raffelt:2006cw ; Arvanitaki:2009fg . Provided that the PQ PT is first order,
the highly massive DM candidates can be naturally generated through the bubble
filtering-out mechanism. Moreover, in this work according to this setup, we
propose an asymmetric DM scenario during the PQ PT.

We use a Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) axion model Dine:1981rt ;
Zhitnitsky:1980tq  supplemented with chiral neutrinos, where one of the
flavors can play the role of DM and we focus on this flavor in our
discussions. Taking loop quantum effects into account at finite temperature,
within the parameter space of the theory we show that the PQ PT can be first
order and consequently bubbles of the broken phase can nucleate and expand.
DMs whose kinetic energy is greater than their masses in the broken phase
enter the bubbles. In addition, Right-Handed (RH) neutrino interactions with
bubbles violate the lepton number. We also consider the possibility of a CP
violation source varying during the PT Baldes:2016gaf ; Cline:2020jre ;
Long:2017rdo . Due to CP-violating effects, we obtain the net dark number
density and hence the observed relic abundance can be found and survive after
the PT. Moreover, the baryon asymmetry can be fulfilled due to asymmetric
visible neutrinos via leptogenesis scenarios Ahmadvand:2020izy ;
Davidson:2008bu .222For baryogenesis scenarios fulfilled via lepton-number
violation at a first order PT, see Long:2017rdo ; Pascoli:2016gkf ;
Cohen:1990py ; Cohen:1990it .

We then explore phenomenological consequences of the model. As a solution to
the strong CP problem, we study effective interactions of the QCD axion in the
model. Regarding the generation of Gravitational Waves (GWs) from first order
PTs, GW direct detection experiments are powerful tools to probe early
Universe events. By obtaining the bubble profile and other required
quantities, including the vacuum energy and duration of the PT, we find the
energy density spectrum of GWs produced from the first order PQ PT.
Furthermore, we discuss the detectability of possible GWs from the PT in the
case of visible QCD axions with PQ symmetry breaking scale around
$100\,\mathrm{TeV}$ Rubakov:1997vp ; Berezhiani:2000gh ; Hook:2014cda ;
Fukuda:2015ana . Finally, we investigate some of effective DM interactions
with Standard Model (SM) particles.

In Section 2, we introduce the model and discuss the PQ PT. In section 3, the
DM relic abundance is obtained within the model. We study phenomenological
aspects of the model in Section 4 and conclude in Section 5.

##  2 The Model

Due to the anomalous U(1) axial symmetry of strong interactions and QCD vacuum
structure, an effective CP-violating term associated with $\theta$-vacuum is
allowed in the Lagrangian. This term contributes to the neutron electric
dipole moment, $d_{n}\approx 3.6\times 10^{-16}\bar{\theta}\,e\,\mathrm{cm}$
Crewther:1979pi , which is experimentally constrained $|d_{n}|<2.9\times
10^{-26}\,e\,\mathrm{cm}~{}(90\%~{}\mathrm{C.L.})$ Baker:2006ts  and thereby
$\bar{\theta}\lesssim 10^{-10}$. This poses the strong CP problem in that
there is no reason in the SM why $\bar{\theta}$ should be very small.

An interesting solution to this problem is based on a global chiral
$\mathrm{U(1)}_{\mathrm{PQ}}$ symmetry proposed by Peccei and Quinn
Peccei:1977hh . In fact, at QCD scales, the interaction of a pseudo-scalar
field, the axion $a$, is effectively added to the Lagrangian,
$(a/f_{a}+\bar{\theta})G\widetilde{G}$, where $f_{a}$ is the axion decay
constant, $G_{\mu\nu}$ is the gluon field strength and $\widetilde{G}$ denotes
its dual. Therefore, the CP-violating $\bar{\theta}$-term can be cancelled via
the vacuum expectation value (vev) of axion at the minimum of its potential,
addressing the strong CP problem. The UV completion of such non-renormalizable
interactions can be constructed in a model invariant under
$\mathrm{U(1)}_{\mathrm{PQ}}$. At some high energy scale, the symmetry is
spontaneously broken and the resulting pseudo-Nambu Goldstone boson would be
matched with the axion so that at low energies the effective interaction can
be induced due to the chiral anomaly and QCD instantons.

In general, considering astrophysical bounds on the PQ symmetry breaking
scale, two types of axion models can be categorized: Kim-Shifman-Vainshtein-
Zakharov (KSVZ) models Kim:1979if ; Shifman:1979if  which contain extra heavy
quarks and PQ scalar fields, carrying the PQ charge, and DFSZ models
Dine:1981rt ; Zhitnitsky:1980tq  in which beside the PQ scalar field, an
additional Higgs field is introduced. We here employ a DFSZ type axion model
extended by three RH neutrinos, where one of the flavors can be regarded as a
DM candidate.

The Lagrangian of the model is given by

$\mathcal{L}\supset|\partial_{\mu}\Phi|^{2}+V(\Phi,H_{u},H_{d})+\overline{N}_{R_{i}}i\partial_{\mu}N_{R_{i}}+y_{ij}\,\Phi\overline{N}_{R_{i}}N_{R_{j}}^{c}$
(1)

where

$\displaystyle V(\Phi,H_{u},H_{d})$
$\displaystyle=\lambda_{\phi}\left(|\Phi|^{2}-v_{\phi}^{2}/2\right)^{2}+\left|H_{d}\right|^{2}\left(\kappa_{d}|\Phi|^{2}-\mu_{d}^{2}\right)+\left|H_{u}\right|^{2}\left(\kappa_{u}|\Phi|^{2}+\mu_{u}^{2}\right)$
(2)
$\displaystyle-\left(\kappa\Phi^{\dagger}H_{u}H_{d}+h.c.\right)+\lambda_{d}\left|H_{d}\right|^{4}+\lambda_{u}\left|H_{u}\right|^{4}+\lambda_{1}\left|H_{u}H_{d}\right|^{2}+\lambda_{2}\left|H_{u}\right|^{2}\left|H_{d}\right|^{2}$

$\Phi$ is the PQ scalar field singlet under the SM gauge symmetry group
$\mathrm{SU(2)}_{L}\times\mathrm{U(1)}_{Y}$, $H_{u}$ and $H_{d}$ denote two
$\mathrm{SU(2)}_{L}$ doublets, and the interaction of singlet RH neutrinos
with $\Phi$ violates the lepton number by two units $\Delta L=2$ in case
$L(N_{R_{i}})=1$.333We assume the coupling of such interactions
$\Phi\overline{N}_{L}N_{R}$ is very small compared to that of L-violating
Yukawa interactions. We also consider a discrete $\mathbb{Z}_{2}$ symmetry
under which all fields are even except for the dark neutrino flavor which is
odd so that the massive DM gets thermally-decoupled after the PT. We assume
that $H_{u}$ interacts with $u$ quarks and the rest of the SM as well as two
visible RH neutrinos are coupled to $H_{d}$

$y_{u}\overline{Q}_{L}H_{u}u_{R}+y_{d}\overline{Q}_{L}H_{d}d_{R}+y_{l}\overline{\Psi}_{L_{l}}H_{d}l_{R}+y_{N_{\alpha}}\overline{\Psi}_{L_{l}}\widetilde{H}_{d}N_{R_{\alpha}}+h.c..$
(3)

The Lagrangian should be invariant under $\mathrm{U(1)}_{\mathrm{PQ}}$. PQ
invariance of $\Phi^{\dagger}H_{u}H_{d}$ implies
$X_{H_{u}}+X_{H_{d}}=X_{\Phi}$ where $X_{X}$ is the PQ charge of a given
field. Moreover, imposing the orthogonality between PQ and corresponding
hypercharge currents, $-X_{H_{u}}v_{u}^{2}+X_{H_{d}}v_{d}^{2}=0$
DiLuzio:2020wdo . Defining $v_{u}=v\sin\theta$ and $v_{d}=v\cos\theta$, where
$v$ is the Electroweak (EW) vev, and choosing $X_{\Phi}=1$, we may determine
all charges through

$\displaystyle
X_{H_{d}}=\sin^{2}\theta,~{}~{}~{}~{}~{}X_{H_{u}}=\cos^{2}\theta,$ (4)

as shown in Table 1. The SM Higgs is indeed
$H=H_{d}\cos\theta+\widetilde{H}_{u}\sin\theta$, where
$\widetilde{H}_{u}=i\sigma_{2}H^{*}_{u}$, and $v^{2}=v_{d}^{2}+v_{u}^{2}$.
Given that
$v_{a}^{2}=\sum_{i}X_{i}v_{i}^{2}=v_{\phi}^{2}+v^{2}\sin^{2}(2\theta)/4$ and
$v_{\phi}\gg v$, we have $v_{a}\simeq v_{\phi}$. We will also show in Section
4 that $v_{a}=f_{a}$ in the model. The PQ symmetry is spontaneously broken by
the vev of $\Phi$, such that $\Phi=1/\sqrt{2}(\phi+f_{a})\exp(ia/f_{a})$ and
the axion field is shifted by the PQ transformation.

Field | $\mathrm{SU(3)}_{c}$ | $\mathrm{SU(2)}_{L}$ | $\mathrm{U(1)}_{Y}$ | $\mathrm{U(1)}_{\mathrm{PQ}}$
---|---|---|---|---
$\Phi$ | $\mathbf{1}$ | $\mathbf{1}$ | 0 | 1
$H_{u}$ | $\mathbf{1}$ | $\mathbf{2}$ | $-\frac{1}{2}$ | $\cos^{2}\theta$
$H_{d}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\frac{1}{2}$ | $\sin^{2}\theta$
$Q_{L}$ | $\mathbf{3}$ | $\mathbf{2}$ | $\frac{1}{6}$ | $\cos^{2}\theta$
$u_{R}$ | $\mathbf{3}$ | $\mathbf{1}$ | $\frac{2}{3}$ | 0
$d_{R}$ | $\mathbf{3}$ | $\mathbf{1}$ | $-\frac{1}{3}$ | 2$\cos^{2}\theta$-1
$\Psi_{L}$ | $\mathbf{1}$ | $\mathbf{2}$ | $-\frac{1}{2}$ | $\frac{1}{2}-\sin^{2}\theta$
$e_{R}$ | $\mathbf{1}$ | $\mathbf{1}$ | -1 | $\frac{1}{2}-2\sin^{2}\theta$
$N_{R}$ | $\mathbf{1}$ | $\mathbf{1}$ | 0 | $\frac{1}{2}$
Table 1: Charge assignment and group representation of the fields in the model

###  2.1 First order Peccei-Quinn phase transition

In order to describe the PT, we first express the potential in terms of
electrically neutral components, $H_{u}=h_{u}/\sqrt{2}$ and
$H_{d}=h_{d}/\sqrt{2}$ at the tree-level

$\displaystyle V$
$\displaystyle=\frac{\lambda_{\phi}}{4}\left(\phi^{2}-f_{a}^{2}\right)^{2}+\frac{1}{2}h_{d}^{2}\left(\frac{\kappa_{d}}{2}\phi^{2}-\mu_{d}^{2}\right)+\frac{1}{2}h_{u}^{2}\left(\frac{\kappa_{u}}{2}\phi^{2}+\mu_{u}^{2}\right)-\frac{\kappa}{\sqrt{2}}\phi
h_{u}h_{d}$ (5)
$\displaystyle+\frac{\lambda_{d}}{4}h_{d}^{4}+\frac{\lambda_{u}}{4}h_{u}^{4}+\frac{\lambda_{1}+\lambda_{2}}{4}h_{u}^{2}h_{d}^{2}.$

We assume the case in which $h_{u}$ vanishes during the PT, hence we will not
consider its dynamics during the PT, though in the loop level its loop effects
appear in the potential for $\phi$ and $h_{d}$. Consequently, there are two
minima: 1- ($\phi=f_{a},~{}~{}h_{d}=h_{u}=0$) and 2-
($h_{d}=\mu_{d}^{2}/\lambda_{d},~{}~{}\phi=h_{u}=0$). For this PT, we study
the first direction.444About the possibility of a supercooling period during
the PT in axion models see VonHarling:2019rgb ; DelleRose:2019pgi ;
Ghoshal:2020vud . To have such a minimum, the potential should satisfy the
condition $V(\mathrm{min}1)<V(\mathrm{min}2)$ from which
$\lambda_{\phi}>\mu_{d}^{4}/(f_{a}^{4}\lambda_{d})$. Moreover, we are
interested in the case in which $m_{N_{i}}>m_{\phi}$ in the broken phase, thus
$y_{i}>\sqrt{2\lambda_{\phi}}$, where $y_{i}$ stands for the Yukawa couplings
of heavy neutrinos in the diagonal form.

To account for quantum and thermal effects on the potential, we obtain the
one-loop effective potential at finite temperature containing the resummed
daisy (ring) diagrams,

$V_{t}=V+V_{\mathrm{CW}}+V_{\mathrm{th}}+V_{r}$ (6)

where at zero temperature, the one-loop Coleman-Weinberg quantum correction
can be written as Coleman:1973jx

$V_{\mathrm{CW}}\left(\phi\right)=\sum_{i}(-1)^{F_{b/f}}g_{i}\frac{m_{i}^{4}\left(\phi\right)}{64\pi^{2}}\left[\ln\left(\frac{m_{i}^{2}\left(\phi\right)}{\Lambda^{2}}\right)-c_{i}\right].$
(7)

Here we take $\Lambda=f_{a}$, $F_{b/f}=1(0)$ for fermions (bosons), $g_{i}$ is
the number of degrees of freedom for a given field, and $c_{i}=3/2(5/2)$ for
scalars and fermions (vectors). Also, the one-loop thermal correction is given
by Quiros:1999jp ; Curtin:2016urg

$V_{T}\left(\phi,T\right)=\sum_{i}(-1)^{F_{b/f}}g_{i}\frac{T^{4}}{2\pi^{2}}J_{b/f}\left[\frac{m_{i}^{2}\left(\phi\right)}{T^{2}}\right]$
(8)

where

$J_{b/f}\left(y^{2}\right)=\int_{0}^{\infty}dx\,x^{2}\ln\left[1\mp
e^{-\sqrt{x^{2}+y^{2}}}\right]$ (9)

and the useful form of these thermal functions in the high temperature limit
is given by

$J_{b}\left(\frac{m^{2}}{T^{2}}\right)=-\frac{\pi^{4}}{45}+\frac{\pi^{2}}{12}\left(\frac{m}{T}\right)^{2}-\frac{\pi}{6}\left(\frac{m^{2}}{T^{2}}\right)^{3/2}-\frac{1}{32}\left(\frac{m}{T}\right)^{4}\ln\left(\frac{m^{2}}{a_{b}T^{2}}\right)+\cdots$
(10)
$J_{f}\left(\frac{m^{2}}{T^{2}}\right)=\frac{7\pi^{4}}{360}-\frac{\pi^{2}}{24}\left(\frac{m}{T}\right)^{2}-\frac{1}{32}\left(\frac{m}{T}\right)^{4}\ln\left(\frac{m^{2}}{a_{f}T^{2}}\right)+\cdots$
(11)

where $\ln(a_{b})=5.4076$ and $\ln(a_{f})=2.6351$. The leading order of multi-
loop corrections is included in the so-called daisy diagram whose contribution
is given by

$V_{r}\left(\phi,T\right)=-\sum_{i\in
b}\frac{g_{i}T}{12\pi}\left(\left[m_{i}^{2}(\phi)+\Pi_{i}(T)\right]^{3/2}-\left[m_{i}^{2}(\phi)\right]^{3/2}\right)$
(12)

where the thermal mass squared of the scalars are

$\Pi_{\phi}(T)=\frac{T^{2}}{6}\left(\kappa_{d}+\kappa_{u}+2\lambda_{\phi}\right),$
(13) $\Pi_{h_{d}}(T)=\frac{T^{2}}{48}\left(9g^{2}+3g^{\prime
2}+\frac{12\lambda_{t}^{2}}{\cos^{2}\theta_{\mathrm{w}}}+24\lambda_{d}+4\kappa_{d}+8\lambda_{1}+8\lambda_{2}\right),$
(14) $\Pi_{h_{u}}(T)=\frac{T^{2}}{48}\left(9g^{2}+3g^{\prime
2}+24\lambda_{u}+4\kappa_{u}+8\lambda_{1}+8\lambda_{2}\right).$ (15)

where $\lambda_{t}\simeq 1$, $g\simeq 0.42$, $g^{\prime}\simeq 0.12$, and
$\cos^{2}\theta_{\mathrm{w}}\simeq 0.77$. We calculate the effective potential
accordingly with some representative values of parameters, e.g., for
$y_{1}=\sqrt{2}$, $y_{2}=\sqrt{2}/2$, $y_{3}=\sqrt{2}\times 10^{-3}$,
$\kappa_{u}=3$, $\lambda_{u}=\lambda_{d}=\lambda_{1}=\lambda_{2}=10^{-2}$ and
$\kappa_{d}<\lambda_{\phi}=10^{-8}$. As shown in Fig. (1), for
$f_{a}=100\,\mathrm{PeV}$, we find the PQ PT is of first order and the two
degenerate phases are separated through a barrier at the critical temperature
$T_{c}=16.88\,\mathrm{PeV}$. Based on these parameter values, the first order
PT is mainly induced by loop correction terms of $h_{u}$. Despite the effect
of daisy diagram corrections which reduce the altitude of the potential energy
barrier and weaken the strength of the PT, the barrier is generated due to the
thermal cubic term, Eq. (10).

Figure 1: For $f_{a}=100\,\mathrm{PeV}$, we display the total potential at
three different temperatures within the aforementioned parameter space of the
model.

At some lower temperature, bubbles of the new phase nucleate, expand and
eventually collide with each other until the Universe is completely
transmitted into the broken phase. Indeed, bubbles are non-perturbative
solutions of the following three-dimensional Euclidean bounce action
quantifying the tunneling process

$S_{3}(T)=\int_{0}^{\infty}4\pi
r^{2}~{}dr\Big{[}\frac{1}{2}\Big{(}\frac{d\phi}{dr}\Big{)}^{2}+V_{t}(\phi,T)\Big{]}.$
(16)

By extremizing the action, we obtain the bounce equation and according to the
following boundary condition, the equation can be numerically solved

$\frac{d^{2}\phi}{dr^{2}}+\frac{2}{r}\frac{d\phi}{dr}=\frac{\partial
V_{t}}{\partial\phi},~{}~{}~{}~{}~{}~{}~{}\left.\frac{d\phi}{dr}\right|_{r=0}=0,~{}~{}~{}~{}\phi(\infty)=0.$
(17)

Using the any bubble code Masoumi:2016wot , as can be seen from Fig. (2), we
obtain the bubble profile connecting the two phases.

Figure 2: The bubble profile as the solution of the bounce equation is shown.
$l_{w}$ denotes the wall width.

Bubbles nucleate when the bubble formation probability, proportional to
$\exp\left(-S_{3}(T)/T\right)$, is of the order of one. Therefore, we can find
the nucleation temperature, $T_{n}$, by the following relation Linde:1980tt

$\frac{S_{3}(T_{n})}{T_{n}}=4\ln\Big{(}\frac{T_{n}}{H_{*}}\Big{)},$ (18)

where $H_{*}\simeq T_{n}^{2}/M_{\mathrm{pl}}$ is the Hubble parameter and
$M_{\mathrm{pl}}\simeq 2.43\times 10^{18}\,\mathrm{GeV}$ is the reduced Planck
mass. As a result, using the obtained solution, from Eq. (18), we find the
nucleation temperature of the PT, $T_{n}\simeq 9.19\,\mathrm{PeV}$.

In the next section, relying on the first order PQ PT and the obtained
parameters, we proceed to calculate the DM relic abundance.

##  3 Dark matter relic abundance

In this section, using the DM filtration mechanism, we will find the net
number density of DM by solving the Boltzmann equation during the PQ PT.
According to the first order PT breaking the PQ symmetry, DM acquires a mass
$m_{N}^{\mathrm{in}}$ ($N\equiv N_{\mathrm{DM}}$) and its interactions can be
put out of equilibrium inside the bubbles. DM particles and antiparticles that
have enough kinetic energy, greater than $m_{N}^{\mathrm{in}}$, enter the
bubbles and otherwise reflected DMs remain in equilibrium with the thermal
bath. Furthermore, we are interest in the possibility of a varying CP-
violating source, similar to the varying bubble profile, during the PT and
because of the L-violating interaction of DMs with bubbles, this process gives
rise to a difference between the number density of DM particles and
antiparticles, so that the net abundance survives after the PT.

We assume the wall is planar, perpendicular to the $z$-axis, and moving with
the bubble wall velocity $v_{w}$ in the negative $z$ direction. Also, assuming
the system has reached a steady state and is translation invariant in $x$ and
$y$, one can write the Liouville operator of phase space distribution function
$f_{N}$ for DMs in the wall frame as

$\mathbf{L}[f_{N}]=\frac{df_{N}}{dt}=\frac{p_{z}}{E}\frac{\partial
f_{N}}{\partial z}+F\frac{\partial f_{N}}{\partial p_{z}}.$ (19)

$F$ is the semiclassical force which is given by Cline:2020jre ; Cline:2000nw
; Fromme:2006wx

$F=-\frac{m_{N}(z)}{E}\frac{\partial m_{N}}{\partial z}\pm
s\frac{m_{N}(z)}{EE_{z}}\left(\frac{\partial m_{N}}{\partial
z}\frac{\partial\theta}{\partial
z}+\frac{m_{N}(z)}{2}\frac{\partial^{2}\theta}{\partial
z^{2}}-\frac{m_{N}^{2}(z)}{2E^{2}}\frac{\partial m_{N}}{\partial
z}\frac{\partial\theta}{\partial z}\right)$ (20)

where

$E_{z}^{2}=E^{2}-\mathbf{p}^{2}_{\|}=m_{N}^{2}+p_{z}^{2}$ (21)

and $\mathbf{p}_{\|}$ denotes the momentum parallel to the wall. The positive
sign $(+)$ is assigned to particles and $(-)$ to antiparticles. Also, boosting
to the frame in which $\mathbf{p}_{\|}=0$ for helicity eigenstates, $s=\pm 1$
Cline:2020jre . The last three terms of (20) are CP-violating sources during
the PT and originated from the complex mass term

$\hat{m}_{N}(z)=y_{i}\phi(z)/\sqrt{2}\exp[i\theta(z)]=m_{N}(z)\exp[i\theta(z)].$
(22)

Based on the bubble profile found in the previous section, we solve in the
following the Boltzmann equation in $T_{n}$ units and hence we can model the
$z$-dependent solution by

$\phi(z)=\frac{A}{2}\left(1+\tanh[Bz]\right)$ (23)

where $A=7.2$, and $B=3.9$. Moreover, we consider the CP violating phase as
Bruggisser:2017lhc

$\theta(z)=\arctan\left(\frac{\Delta\theta}{2}\left[1+\tanh(\frac{z}{l_{w}})\right]\right).$
(24)

Using the ansatz $f_{N}=\mathcal{A}(z,p_{z})\exp(-E^{p}/T)$ for the
distribution function Baker:2019ndr , we describe the deviation from
equilibrium by $\mathcal{A}$ which also includes the chemical potential. The
energy of a particle in the plasma frame $E^{p}$ is related to the one in the
wall frame as

$E_{p}^{p}=\gamma_{w}(E-v_{w}p_{z}),~{}~{}~{}~{}\gamma_{w}=\frac{1}{\sqrt{1-v_{w}^{2}}}.$
(25)

Integrating over $p_{x}$ and $p_{y}$, for RH chirality the Liouville operator
would be

$\displaystyle g_{N}$
$\displaystyle\int\frac{dp_{x}dp_{y}}{(2\pi)^{2}}\mathbf{L}\left[f_{N}\right]\sim$
(26) $\displaystyle\sim$
$\displaystyle\left[\left(\frac{p_{z}}{m_{N}}\frac{\partial}{\partial
z}-\left(\frac{\partial m_{N}}{\partial
z}\pm\left(\frac{1}{2m_{N}}\frac{\partial m_{N}}{\partial
z}\frac{\partial\theta}{\partial
z}+\frac{1}{2}\frac{\partial^{2}\theta}{\partial
z^{2}}\right)\right)\frac{\partial}{\partial p_{z}}-\left(\frac{\partial
m_{N}}{\partial
z}\right)\frac{v_{w}}{T}\right)\mathcal{A}_{\pm}\left(z,p_{z}\right)\right]$
$\displaystyle\times$
$\displaystyle\frac{g_{N}m_{N}T}{2\pi}\exp\left(\frac{v_{w}p_{z}-\sqrt{m_{N}^{2}+p_{z}^{2}}}{T}\right)$

where because of the friction effects produced from reflected and penetrated
(anti)particles, we used non-relativistic bubble wall velocities, so that
$E^{p}\simeq E-v_{w}p_{z}$. According to the Boltzmann equation, we have
$\mathbf{L}\left[f_{N}\right]=\mathbf{C}\left[f_{N}\right]$, where the
integration over the collision term is given by Baker:2019ndr

$g_{N}\int\frac{dp_{x}dp_{y}}{(2\pi)^{2}}\mathbf{C}\left[f_{N}\right]=-g_{N}^{2}\left[\mathcal{A}_{\pm}\left(z,p_{z}\right)-1\right]\int\frac{dp_{x}dp_{y}}{(2\pi)^{2}2E_{p}^{p}}\frac{d^{3}q}{(2\pi)^{3}2E_{q}^{p}}4\tilde{F}\sigma\exp\left(-\frac{E_{p}^{p}+E_{q}^{p}}{T}\right)$
(27)

where $\sigma$ is the spin-averaged cross section of such processes
$N_{R}(p^{p})+N_{R}(q^{p})\rightarrow\phi(k^{p})+\phi(l^{p})$ and
$\tilde{F}=\sqrt{((p+q)^{2}-2m_{N}^{2})^{2}-4m_{N}^{4}}/2$. To solve the
Boltzmann equation, we use the method explained in Baker:2019ndr  and rewrite
the equation as

$a_{\pm}\left(z,p_{z}\right)\frac{\partial\mathcal{A}_{\pm}}{\partial
z}+b_{\pm}\left(z,p_{z}\right)\frac{\partial\mathcal{A}_{\pm}}{\partial
p_{z}}=c_{\pm}\left(\mathcal{A}_{\pm},z,p_{z}\right)$ (28)

where

$\frac{dz(\lambda)}{d\lambda}=a_{\pm}\left(z,p_{z}\right),~{}~{}~{}~{}~{}\frac{dp_{z}(\lambda)}{d\lambda}=b_{\pm}\left(z,p_{z}\right),~{}~{}~{}~{}~{}\frac{d\mathcal{A}_{\pm}\left(z(\lambda),p_{z}(\lambda)\right)}{d\lambda}=c_{\pm}\left(\mathcal{A}_{\pm}(\lambda),z(\lambda),p_{z}(\lambda)\right)$
(29)

and $\lambda$ parameterizes a given curve on the ($z$, $p_{z}$) plane.
Therefore, from Eqs. (26, 27), imposing appropriate boundary conditions, we
first obtain $z(\lambda)$ and $p_{z}(\lambda)$, then we find
$\mathcal{A}_{\pm}$ for a region of $z$ and $p_{z}$ values.

For (anti)particles outside the bubble traveling with $p_{z}>0$, the boundary
condition is set as

$\lim_{z\rightarrow-\infty}\mathcal{A}_{\pm}\rightarrow 1,$ (30)

ensuring the equilibrium phase space distribution function far away from the
bubble wall. Using the obtained curves of ($z\rightarrow-\infty$,
$p_{z}>m_{N}^{\mathrm{in}}$), for (anti)particles deep inside the bubble
($z\rightarrow\infty$) with $p_{z}<0$ and similar dynamics, we set the
boundary condition as Baker:2019ndr

$\lim_{z\rightarrow\infty}\mathcal{A}_{\pm}(p_{z})=\lim_{z\rightarrow\infty}\mathcal{A}_{\pm}(-p_{z}).$
(31)

Interpolating between different solutions obtained for several curves, we find
$\mathcal{A}_{\pm}$ as displayed in Fig. (3). Finally, integrating over
$p_{z}$ at deep inside the bubble, we obtain the net number density $\Delta
n_{N}/T_{n}^{3}$ in $T_{n}$ units, where $\Delta n_{N}\equiv
n_{N}^{+}-n_{N}^{-}$. We take $v_{w}=0.01$ in the following calculations,
however, for other values, for example $v_{w}=0.1$, the final result can be
equivalently obtained by slightly different values of parameters.

Figure 3: For two bubble velocities, $v_{w}=0.01$ (left figure) and
$v_{w}=0.1$ (right figure), the factor $\mathcal{A}_{+}$ for DM particles is
displayed. Since to obtain the observed relic abundance, it requires
$\Delta\theta\ll 1$ and hence since $\mathcal{A}_{+}$ and $\mathcal{A}_{-}$
are displayed almost the same, $\mathcal{A}_{+}$ is only shown.

The DM relic abundance can be found in the asymmetric case by the following
relation555During the PT, L-violating processes with CP-violating effects,
present during this period, occur inside the bubbles, thus the net number
density is used in the DM relic abundance relation.

$h^{2}\Omega_{\mathrm{DM}}\simeq\frac{m^{\mathrm{in}}_{N}\Delta
n_{N}~{}g_{*0}T_{0}^{3}}{3M_{\mathrm{pl}}^{2}\left(H_{0}/h\right)^{2}g_{*}T_{n}^{3}}$
(32)

where $g_{*0}=3.9$, $T_{0}\simeq 0.235\,\mathrm{meV}$, $g_{*}\sim 110$ and the
present Hubble parameter
$H_{0}=100\,h~{}\mathrm{km}\,\mathrm{sec}^{-1}\,\mathrm{Mpc}^{-1}$. As a
result, the observed relic abundance, $h^{2}\Omega_{\mathrm{DM}}\simeq 0.12$
Planck:2018vyg , can be obtained for, e.g., $y_{1}=\sqrt{2}$ or
$m^{\mathrm{in}}_{N}\simeq 66.3\,\mathrm{PeV}$ and from the attained $\Delta
n_{N}/T_{n}^{3}$ which is determined by $\Delta\theta\sim 10^{-12}$.

Furthermore, analogous to the analytic approach of Baker:2019ndr , estimated
the DM number density proportional to its equilibrium abundance inside the
bubble $n_{\mathrm{eq}}/4$ where

$n_{\mathrm{eq}}=\frac{g_{N}\left(m_{N}^{\mathrm{in}}T_{n}\right)^{3/2}}{(2\pi)^{3/2}}e^{-m_{N}^{\mathrm{in}}/T_{n}},$
(33)

here, considering the net number density as $\Delta
n_{N}\sim\delta\theta\,n_{\mathrm{eq}}$, with $T_{n}\simeq 9.19\,\mathrm{PeV}$
and $m^{\mathrm{in}}_{N}\simeq 66.3\,\mathrm{PeV}$, one can obtain the
observed relic abundance for $\delta\theta\simeq 0.19\times 10^{-12}$.

##  4 Phenomenology

In this section we study phenomenological consequences of the model.

###  4.1 Axion couplings

The Lagrangian of the QCD axion as a dynamical solution for the strong CP
problem is given such that the axion effective interaction with gluons and its
vev give rise to cancellation of the $\bar{\theta}$ term. Along this
direction, the anomalous PQ current of the axial $U(1)_{\mathrm{PQ}}$ symmetry
is given as follows

$\partial^{\mu}J_{\mu}^{\mathrm{PQ}}=\frac{g_{s}^{2}\mathcal{N}}{16\pi^{2}}G\widetilde{G}+\frac{e^{2}\mathcal{E}}{16\pi^{2}}F\widetilde{F}$
(34)

where $g_{s}$ is the $\mathrm{SU(3)}_{c}$ strong coupling, $e$ is the electric
charge, $F_{\mu\nu}$ is the photon field strength, $\mathcal{N}$ and
$\mathcal{E}$ are the QCD and electromagnetic anomaly coefficients,
respectively. Due to the spontaneous PQ symmetry breaking, the effective
interactions of the axion are expressed as

$\mathcal{L}_{a}\supset\frac{g_{s}^{2}}{32\pi^{2}f_{a}}aG\widetilde{G}+\frac{g_{a\gamma}}{4}aF\widetilde{F}-\frac{c^{0}_{f}}{2f_{a}}\partial_{\mu}a\bar{f}\gamma^{\mu}\gamma_{5}f-\frac{c^{0}_{N_{R}}}{2f_{a}}\partial_{\mu}a\overline{N}_{R}\gamma^{\mu}N_{R}$
(35)

where

$g_{a\gamma}=\frac{e^{2}}{8\pi^{2}f_{a}}\frac{\mathcal{E}}{\mathcal{N}}-\frac{e^{2}}{8\pi^{2}f_{a}}\left(\frac{2}{3}\frac{4m_{d}+m_{u}}{m_{d}+m_{u}}\right)$
(36)
$c^{0}_{f}=\frac{X_{f_{R}}-X_{f_{L}}}{2\mathcal{N}}=\frac{X_{H_{f}}}{2\mathcal{N}},~{}~{}~{}~{}~{}~{}~{}c^{0}_{N_{R}}=\frac{X_{\Phi}}{2\mathcal{N}}$
(37) $\mathcal{N}=\frac{n_{g}}{2}\left(X_{H_{u}}+X_{H_{d}}\right)$ (38)
$\mathcal{E}=n_{g}\left(3\left(\frac{2}{3}\right)^{2}X_{H_{u}}+3\left(-\frac{1}{3}\right)^{2}X_{H_{d}}+\left(-1\right)^{2}X_{H_{d}}\right)$
(39)

$n_{g}$ is the number of SM fermion $f$ generation,
$f_{a}=v_{a}/(2\mathcal{N})$, $m_{u}$ and $m_{d}$ are $u$ and $d$ quark
masses, respectively. The cosmological Domain Wall (DW) problem
Zeldovich:1974uw  can be avoided if $N_{\mathrm{DW}}\equiv 2\mathcal{N}=1$
DiLuzio:2020wdo . As explained in Section 2, $X_{H_{u}}+X_{H_{d}}=X_{\Phi}$.
Therefore, according to Table 1, $\mathcal{N}=1/2$ and thereby $f_{a}=v_{a}$
and $N_{\mathrm{DW}}=1$ can be fulfilled if one family is charged under the PQ
symmetry.

Moreover, depending on the fermion type, $\tan\theta\in[0.25,170]$
DiLuzio:2020wdo , thus we obtain $\mathcal{E}=4/3$,
$c_{u}^{0}=\cos^{2}\theta$, $c_{d}^{0}=\sin^{2}\theta$, $c_{N_{R}}^{0}=1$ from
which other couplings including the axion-pion coupling can de determined.

###  4.2 Gravitational wave signals

Regarding growing efforts and progresses on the GW direct detection,
cosmological PTs, specially first order PTs as a source for the GW radiation,
are important events that should be studied Mazumdar:2018dfl ;
Ahmadvand:2017xrw ; Ahmadvand:2017tue ; Abedi:2019msi ; Ahmadvand:2020fqv .

During cosmological first order PTs and bubble evolution processes, three
sources for the GW generation have been proposed: bubble collision, sound
waves and magnetohydrodynamic (MHD) turbulence Kosowsky:1992rz ;
Kamionkowski:1993fg ; Kosowsky:2001xp ; Hindmarsh:2015qta . The GW energy
density spectrum can be characterized by some of important PT quantities
including the bubble wall velocity, the released latent heat, and duration of
the PT, calculated at the nucleation temperature.666Since we do not
differentiate between the temperature at which GWs are generated and the
nucleation temperature in this prompt PT, we calculate the parameters at
$T_{n}$.

The parameter which is associated with the latent heat and appears in the GW
energy density computation is the ratio of the vacuum energy density to the
thermal energy density,

$\alpha=\frac{\epsilon}{\frac{\pi^{2}}{30}g_{*}T_{n}^{4}},~{}~{}~{}~{}~{}~{}\epsilon=\left(\Delta
V_{t}(T)-T\frac{d\Delta V_{t}(T)}{dT}\right)\Bigg{|}_{T=T_{n}},$ (40)

where $\Delta V_{t}(T_{n})=V_{t}[0,T_{n}]-V_{t}[v_{t}(T_{n}),T_{n}]$ and
$v_{t}(T_{n})$ is the true vacuum at $T_{n}$. There is a critical value of
$\alpha$, denoted by $\alpha_{\infty}$, such that for $\alpha>\alpha_{\infty}$
bubbles can run away Caprini:2015zlo , where

$\alpha_{\infty}=\frac{30}{24\pi^{2}}\frac{\sum_{i}n_{i}\Delta
m_{i}^{2}}{g_{*}T_{n}^{2}},$ (41)

$n_{i}\,(n_{i}/2)$ is the number of degrees of freedom for boson (fermion)
species, and $\Delta m_{i}^{2}$ is the squared mass difference of particles
between the broken and symmetric phase at the nucleation temperature. Another
key parameter related to the inverse of PT duration is calculated by

$\frac{\beta}{H_{*}}=T_{n}\frac{d}{dT}\left(\frac{S_{3}(T)}{T}\right)\Bigg{|}_{T_{n}}.$
(42)

Based on the analysis described in Section 2.1, we found $T_{n}\simeq
9.19\,\mathrm{PeV}$. At this temperature we obtain $\alpha\simeq 0.55$,
$\alpha_{\infty}\simeq 0.6$. In addition, from Eq. (42) and the obtained
bounce action we find $\beta/H_{*}\simeq 246$.

In the non-runaway case where $\alpha<\alpha_{\infty}$, dominant contributions
to GWs come from sound waves and MHD turbulence, i.e., $h^{2}\Omega(f)\simeq
h^{2}\Omega_{\mathrm{sw}}+h^{2}\Omega_{\mathrm{tu}}$ where Hindmarsh:2015qta ;
Caprini:2009yp

$h^{2}\Omega_{\mathrm{sw}}(f)=2.65\times
10^{-6}\Big{(}\frac{H_{*}}{\beta}\Big{)}\Big{(}\frac{\kappa_{\mathrm{sw}}\alpha}{1+\alpha}\Big{)}^{2}\Big{(}\frac{100}{g_{*}}\Big{)}^{\frac{1}{3}}v_{w}~{}S_{\mathrm{sw}}(f),$
(43) $h^{2}\Omega_{\mathrm{tu}}(f)=3.35\times
10^{-4}\Big{(}\frac{H_{*}}{\beta}\Big{)}\Big{(}\frac{\kappa_{\mathrm{tu}}\alpha}{1+\alpha}\Big{)}^{\frac{3}{2}}\Big{(}\frac{100}{g_{*}}\Big{)}^{\frac{1}{3}}v_{w}~{}S_{\mathrm{tu}}(f),$
(44)

where the spectral shapes are given by Caprini:2015zlo ,

$\displaystyle S_{\mathrm{sw}}(f)$ $\displaystyle=$
$\displaystyle\Big{(}\frac{f}{f_{\mathrm{sw}}}\Big{)}^{3}\Big{(}\frac{7}{4+3(\frac{f}{f_{\mathrm{sw}}})^{2}}\Big{)}^{\frac{7}{2}},$
(45) $\displaystyle S_{\mathrm{tu}}(f)$ $\displaystyle=$
$\displaystyle\frac{(\frac{f}{f_{\mathrm{tu}}})^{3}}{(1+\frac{f}{f_{\mathrm{tu}}})^{\frac{11}{3}}(1+\frac{8\pi
f}{h_{*}})},$ (46)

with

$h_{*}=16.5\times
10^{-2}[\mathrm{Hz}]\Big{(}\frac{T_{n}}{\mathrm{PeV}}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}}.$
(47)

The red-shifted peak frequencies in the spectral shapes are given by

$f_{\mathrm{sw}}=19\times
10^{-2}[\mathrm{Hz}]\Big{(}\frac{1}{v_{w}}\Big{)}\Big{(}\frac{\beta}{H_{*}}\Big{)}\Big{(}\frac{T_{n}}{\mathrm{PeV}}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}},$
(48) $f_{\mathrm{tu}}=27\times
10^{-2}[\mathrm{Hz}]\Big{(}\frac{1}{v_{w}}\Big{)}\Big{(}\frac{\beta}{H_{*}}\Big{)}\Big{(}\frac{T_{n}}{\mathrm{PeV}}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}}.$
(49)

In this case, the bubble wall velocity reaches to a subluminal value and
depending on the velocity, the efficiency factor for the conversion of the
latent heat to the plasma motion would be different. We study the effect of
various bubble wall velocities and combustion modes on the GW signals. For
small velocities, $v_{w}\ll c_{s}=1/\sqrt{3}$, $v_{w}=c_{s}$ which is the case
of transition from subsonic to supersonic deflagrations, and the large
velocity limit, $v_{w}\rightarrow 1$, this efficiency factor is expressed
respectively as Espinosa:2010hh

$\kappa^{A}_{v}=\frac{6.9\alpha~{}v_{w}^{6/5}}{1.36-0.037\sqrt{\alpha}+\alpha},~{}~{}~{}~{}~{}v_{w}\ll
c_{s},$ (50)
$\kappa^{B}_{v}=\frac{\alpha^{2/5}}{0.017+(0.997+\alpha)^{2/5}},~{}~{}~{}~{}~{}v_{w}=c_{s},$
(51)
$\kappa^{C}_{v}=\frac{\alpha}{0.73+0.083\sqrt{\alpha}+\alpha},~{}~{}~{}~{}~{}v_{w}\rightarrow
1.$ (52)

Moreover, for subsonic deflagrations $v_{w}\lesssim c_{s}$, Jouguet
detonations, and detonations with $v_{w}\gtrsim v_{J}$, the efficiency factor
is given by the following fits Espinosa:2010hh

$\kappa_{v}\left(v_{w}\lesssim
c_{s}\right)=\frac{c_{s}^{11/5}\kappa_{v}^{A}\kappa_{v}^{B}}{\left(c_{s}^{11/5}-v_{w}^{11/5}\right)\kappa_{v}^{B}+v_{w}\,c_{s}^{6/5}\kappa_{v}^{A}},$
(53) $\kappa^{J}_{v}=\frac{\sqrt{\alpha}}{0.135+\sqrt{0.98+\alpha}},\quad
v_{w}=v_{J}=\frac{\sqrt{2\alpha/3+\alpha^{2}}+\sqrt{1/3}}{1+\alpha},$ (54)
$\kappa_{v}\left(v_{w}\gtrsim
v_{J}\right)=\frac{\left(v_{J}-1\right)^{3}\left(v_{J}/v_{w}\right)^{5/2}\kappa_{v}^{J}\kappa_{v}^{C}}{\left[\left(v_{J}-1\right)^{3}-\left(v_{w}-1\right)^{3}\right]v_{J}^{5/2}\kappa_{v}^{J}+\left(v_{w}-1\right)^{3}\kappa_{v}^{C}}.$
(55)

The fraction of plasma motion in turbulence,
$\varepsilon=\kappa_{\mathrm{tu}}/\kappa_{v}$, can be set as
$\varepsilon=0.05$ Hindmarsh:2015qta , so that the main source contributing to
the GW energy density from the plasma motion is the source of sound waves,
$\kappa_{\mathrm{sw}}=(1-\varepsilon)\kappa_{v}$.

Eventually, substituting key parameters in Eqs. (43, 44), we find the GW
spectrum for different cases of the bubble wall velocity. As shown in Fig.
(4), the GWs of the first order PQ PT can be within the reach of ongoing third
generation ground-based detectors such as Einstein Telescope (ET) Hild:2008ng
; Sathyaprakash:2012jk

Figure 4: GW energy density spectra of the first order PQ PT with
$f_{a}=100\,\mathrm{PeV}$ for different bubble wall velocities and combustion
regimes are displayed. For $v_{w}=0.9$, we use Eq. (55) for the efficiency
factor. The dashed curve is the approximate prospective sensitivity curve of
ET Hild:2008ng ; Sathyaprakash:2012jk . We show that the GW signals may be
within the reach of this planned GW detector.

####  4.2.1 Gravitational wave signals within a heavy axion case

It is also interesting to study the possibility of the so-called visible axion
models in which the PQ symmetry breaking scale can be $f_{a}\sim
1-100\,\mathrm{TeV}$ Fukuda:2015ana .

We here only explore the detectability of possible GWs from the first order PQ
PT in a heavy axion model. Such class of models can be justified by
considering an additional $\mathbb{Z}_{\mathbf{N}}$ mirror symmetry with
$\mathbf{N}$ mirror worlds, so that with additional confining sectors the
axion potential would change,
$(a/f_{a}+\bar{\theta}_{\mathrm{eff}})G_{k}\widetilde{G}_{k}$ Rubakov:1997vp ;
Berezhiani:2000gh ; Hook:2014cda ; Fukuda:2015ana . Consequently, in this case
the axion mass would be Fukuda:2015ana

$m_{a}\simeq\frac{\sqrt{z^{\prime}}}{1+z^{\prime}}\frac{f_{\pi^{\prime}}m_{\pi^{\prime}}}{f_{a}}$
(56)

where $m_{\pi^{\prime}}$ and $f_{\pi^{\prime}}$ are the pion mass and its
decay constant in the mirror sectors and
$z^{\prime}=m_{u^{\prime}}/m_{d^{\prime}}$. For example taking
$z^{\prime}=z\sim 0.5$, $m_{\pi^{\prime}}=135\,\mathrm{GeV}$,
$f_{\pi^{\prime}}=93\,\mathrm{GeV}$ and $f_{a}=100\,\mathrm{TeV}$, the axion
mass will be $m_{a}\sim 60\,\mathrm{MeV}$.

Connecting the mirror worlds by the axion DiLuzio:2021pxd , in the UV axion
model the PQ scalar can be transformed as $\phi\rightarrow\exp(2\pi
i/\mathbf{N})\phi$. Thus, considering $f_{a}=100\,\mathrm{TeV}$, the PT
parameters can be retained at $T_{n}\simeq 9.19\,\mathrm{TeV}$. As a result,
as displayed in Fig. (5), in this case the GW signals of the PQ PT can be
probed by the future space-based interferometers including DECIGO and BBO
Yagi:2011wg ; Kudoh:2005as .

We leave further investigations on these models for another future work.

Figure 5: GW energy density spectra of the first order PQ PT with
$f_{a}=100\,\mathrm{TeV}$ for different bubble wall velocities are displayed.
Dashed curves are the expected sensitivity curves of proposed space-based
DECIGO and BBO detectors Yagi:2011wg ; Kudoh:2005as , which may probe the
generated GW signals.

###  4.3 Effective dark matter interactions

Depending on DM interactions with the visible sector, a model can be also
tested by direct detection experiments and indirect measurements. In this
section we illustrate some of these processes, particularly focusing on the
DM-SM neutrino interaction. However, we first consider elastic
$N\,e\rightarrow N^{\prime}\,e^{\prime}$ scattering (Fig(̇6)) and explore the
electron recoil energy.

Figure 6: Left: The DM annihilation to SM neutrinos is shown, where $h$ is the
Higgs field. Right: The DM-SM fermion scattering is shown.

DMs move at non-relativistic velocities, $v_{D}\sim 10^{-3}$, in typical DM
halos Smith:2006ym . Considering these speeds for the PeV-scale DM, one can
calculate the transfer recoil energy as follows Kannike:2020agf

$E_{r}\equiv E_{e}-E_{e^{\prime}}\simeq 2m_{e}v_{D}(v_{D}-v_{e}),$ (57)

hence in the lab frame for $v_{e}=0$ and $m_{e}\simeq 0.5\,\mathrm{MeV}$, the
electron recoil energy would be $E_{r}\simeq\mathrm{eV}$. Therefore, such non-
relativistic massive DMs cannot explain the excess events in a
$1-7\,\mathrm{keV}$ recoil energy reported by the XENON1T experiment
XENON:2020rca .777For an explanation for these events, see for example
Shakeri:2020wvk .

However, the heavy DMs can be considered as a source for producing ultra-high
energy neutrinos Arguelles:2019ouk . To show this process according to the
model (Fig. (6)), we consider the annihilation of DM to SM neutrinos
$N(p)\,N(q)\rightarrow v_{l}(k)\,v_{l}(g)$. In the center of mass frame, for
the non-relativistic DM, $v_{D}\sim 10^{-3}$, we have
$|\mathbf{p}|=|\mathbf{q}|\sim m_{N}v_{D}$ and hence the energy of incoming
particles $E\sim m_{N}$. Also, for $m_{v_{l}}\sim 0$ and due to the energy and
momentum conservation, $E_{v_{l}}\sim|\mathbf{k}|\sim m_{N}$. Finally,
calculating the amplitude of the process, $\mathcal{M}$, we can find the cross
section as follows

$\sigma\sim\frac{|\mathbf{k}|\,|\mathcal{M}|^{2}}{64\pi
E^{2}\,|\mathbf{p}|}\sim
10^{-48}\,\mathrm{cm}^{2}\left(\frac{\mathcal{Y}}{10^{-2}}\right)^{2}\left(\frac{m_{N}}{1\,\mathrm{PeV}}\right)^{-2}$
(58)

where $\mathcal{Y}\equiv y_{1}y_{\alpha}y_{N_{\alpha}}^{2}$ and we take
$m_{N}=y_{1}f_{a}/\sqrt{2}\gg m_{\phi}$. Therefore, as obtained from Eq. (58),
for $E_{v_{l}}\sim m_{N}=1\,\mathrm{PeV}$, the cross section is compatible
with the bound on DM annihilation processes Chianese:2021htv .

##  5 Conclusion

The bubble filtration mechanism provides a setup which allows DM masses above
$100\,\mathrm{TeV}$, the GK bound constrained the DM mass within the thermal
freeze-out mechanism. In this work, using the filtering mechanism, we have
presented an asymmetric DM scenario during the PQ PT through which DMs can
naturally acquire these large masses. Based on a QCD axion model extended by
chiral neutrinos, where one of the flavors plays the role of DM, we find the
one-loop finite temperature effective potential and show that the PT can be
first order within the parameter space of the theory. We obtain the profile of
bubbles nucleated during the PT at the PQ symmetry breaking scale around
$100\,\mathrm{PeV}$. Relying on the lepton number violating interaction of RH
neutrinos, we have shown an asymmetry between heavy neutrinos and
antineutrinos, provided that a CP-violating source is imposed during the PT.
Indeed, we solve the Boltzmann equation numerically, and find the net number
density as well as the observed relic abundance of the DM. It is interesting
to note that because of CP violation effects varying during the PT, the
scenario is not restricted to $m_{N}\gg T_{n}$ for obtaining the relic
abundance and also the resulting net abundance remains after the PT.

As for resolving the strong CP problem in the model, we have then obtained the
QCD axion couplings. Furthermore, calculating the vacuum energy and duration
of the PT at the nucleation temperature, we find the energy density spectrum
of GWs generated from the PQ PT for different combustion modes. We show that
the signals can be detected by the future ground-based detectors such as ET.
In particular, we have investigated possible GWs of the first order PQ PT for
a class of heavy axion models at $f_{a}=100\,\mathrm{TeV}$. We show the GW
signals in this case can be explored by DECIGO and BBO detectors. Eventually,
considering the annihilation of DM to SM neutrinos, we compute the cross
section which is consistent with the constraint on DM annihilation processes
and show that the interaction can be regarded as a source for the ultra-high
energy, PeV scale, neutrinos.

###### Acknowledgements.

I would like to thank Soroush Shakeri for helpful comments and discussions.

## References

  * (1) G. Bertone and D. Hooper, “History of dark matter,” Rev. Mod. Phys. 90, no.4, 045002 (2018) [arXiv:1605.04909 [astro-ph.CO]].
  * (2) G. Bertone, D. Hooper and J. Silk, “Particle dark matter: Evidence, candidates and constraints,” Phys. Rept. 405, 279-390 (2005) [arXiv:hep-ph/0404175 [hep-ph]].
  * (3) R. Allahverdi, B. Dutta and K. Sinha, “Cladogenesis: Baryon-Dark Matter Coincidence from Branchings in Moduli Decay,” Phys. Rev. D 83, 083502 (2011) [arXiv:1011.1286 [hep-ph]].
  * (4) M. Ahmadvand, “Matter and dark matter asymmetry from a composite Higgs model,” Eur. Phys. J. C 81, no.4, 358 (2021) [arXiv:2010.10121 [hep-ph]].
  * (5) E. W. Kolb, D. J. H. Chung and A. Riotto, “WIMPzillas!,” AIP Conf. Proc. 484, no.1, 91-105 (1999) [arXiv:hep-ph/9810361 [hep-ph]].
  * (6) D. E. Kaplan, M. A. Luty and K. M. Zurek, “Asymmetric Dark Matter,” Phys. Rev. D 79, 115016 (2009) [arXiv:0901.4117 [hep-ph]].
  * (7) M. J. Baker, J. Kopp and A. J. Long, “Filtered Dark Matter at a First Order Phase Transition,” Phys. Rev. Lett. 125, no.15, 151102 (2020) [arXiv:1912.02830 [hep-ph]].
  * (8) D. Chway, T. H. Jung and C. S. Shin, “Dark matter filtering-out effect during a first-order phase transition,” Phys. Rev. D 101, no.9, 095019 (2020) [arXiv:1912.04238 [hep-ph]].
  * (9) R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38, 1440-1443 (1977).
  * (10) R. D. Peccei, “The Strong CP problem and axions,” Lect. Notes Phys. 741, 3-17 (2008) [arXiv:hep-ph/0607268 [hep-ph]].
  * (11) G. G. Raffelt, “Astrophysical axion bounds,” Lect. Notes Phys. 741, 51-71 (2008) [arXiv:hep-ph/0611350 [hep-ph]].
  * (12) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, “String Axiverse,” Phys. Rev. D 81, 123530 (2010) [arXiv:0905.4720 [hep-th]].
  * (13) M. Dine, W. Fischler and M. Srednicki, “A Simple Solution to the Strong CP Problem with a Harmless Axion,” Phys. Lett. B 104, 199-202 (1981).
  * (14) A. R. Zhitnitsky, “On Possible Suppression of the Axion Hadron Interactions. (In Russian),” Sov. J. Nucl. Phys. 31, 260 (1980).
  * (15) I. Baldes, T. Konstandin and G. Servant, “Flavor Cosmology: Dynamical Yukawas in the Froggatt-Nielsen Mechanism,” JHEP 12, 073 (2016) [arXiv:1608.03254 [hep-ph]].
  * (16) J. M. Cline and K. Kainulainen, “Electroweak baryogenesis at high bubble wall velocities,” Phys. Rev. D 101, no.6, 063525 (2020) [arXiv:2001.00568 [hep-ph]].
  * (17) A. J. Long, A. Tesi and L. T. Wang, “Baryogenesis at a Lepton-Number-Breaking Phase Transition,” JHEP 10, 095 (2017) [arXiv:1703.04902 [hep-ph]].
  * (18) S. Davidson, E. Nardi and Y. Nir, “Leptogenesis,” Phys. Rept. 466, 105-177 (2008) [arXiv:0802.2962 [hep-ph]].
  * (19) S. Pascoli, J. Turner and Y. L. Zhou, “Baryogenesis via leptonic CP-violating phase transition,” Phys. Lett. B 780, 313-318 (2018) [arXiv:1609.07969 [hep-ph]].
  * (20) A. G. Cohen, D. B. Kaplan and A. E. Nelson, “WEAK SCALE BARYOGENESIS,” Phys. Lett. B 245, 561-564 (1990).
  * (21) A. G. Cohen, D. B. Kaplan and A. E. Nelson, “Baryogenesis at the weak phase transition,” Nucl. Phys. B 349, 727-742 (1991).
  * (22) V. A. Rubakov, “Grand unification and heavy axion,” JETP Lett. 65, 621-624 (1997) [arXiv:hep-ph/9703409 [hep-ph]].
  * (23) Z. Berezhiani, L. Gianfagna and M. Giannotti, “Strong CP problem and mirror world: The Weinberg-Wilczek axion revisited,” Phys. Lett. B 500, 286-296 (2001) [arXiv:hep-ph/0009290 [hep-ph]].
  * (24) A. Hook, “Anomalous solutions to the strong CP problem,” Phys. Rev. Lett. 114, no.14, 141801 (2015) [arXiv:1411.3325 [hep-ph]].
  * (25) H. Fukuda, K. Harigaya, M. Ibe and T. T. Yanagida, “Model of visible QCD axion,” Phys. Rev. D 92, no.1, 015021 (2015) [arXiv:1504.06084 [hep-ph]].
  * (26) R. J. Crewther, P. Di Vecchia, G. Veneziano and E. Witten, “Chiral Estimate of the Electric Dipole Moment of the Neutron in Quantum Chromodynamics,” Phys. Lett. B 88, 123 (1979) [erratum: Phys. Lett. B 91, 487 (1980)].
  * (27) C. A. Baker, D. D. Doyle, P. Geltenbort, K. Green, M. G. D. van der Grinten, P. G. Harris, P. Iaydjiev, S. N. Ivanov, D. J. R. May and J. M. Pendlebury, et al. “An Improved experimental limit on the electric dipole moment of the neutron,” Phys. Rev. Lett. 97, 131801 (2006) [arXiv:hep-ex/0602020 [hep-ex]].
  * (28) J. E. Kim, “Weak Interaction Singlet and Strong CP Invariance,” Phys. Rev. Lett. 43, 103 (1979).
  * (29) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, “Can Confinement Ensure Natural CP Invariance of Strong Interactions?,” Nucl. Phys. B 166, 493-506 (1980).
  * (30) L. Di Luzio, M. Giannotti, E. Nardi and L. Visinelli, “The landscape of QCD axion models,” Phys. Rept. 870, 1-117 (2020) [arXiv:2003.01100 [hep-ph]].
  * (31) B. Von Harling, A. Pomarol, O. Pujolàs and F. Rompineve, “Peccei-Quinn Phase Transition at LIGO,” JHEP 04, 195 (2020) [arXiv:1912.07587 [hep-ph]].
  * (32) L. Delle Rose, G. Panico, M. Redi and A. Tesi, “Gravitational Waves from Supercool Axions,” JHEP 04, 025 (2020) [arXiv:1912.06139 [hep-ph]].
  * (33) A. Ghoshal and A. Salvio, “Gravitational waves from fundamental axion dynamics,” JHEP 12, 049 (2020) [arXiv:2007.00005 [hep-ph]].
  * (34) S. R. Coleman and E. J. Weinberg, “Radiative Corrections as the Origin of Spontaneous Symmetry Breaking,” Phys. Rev. D 7, 1888-1910 (1973).
  * (35) M. Quiros, “Finite temperature field theory and phase transitions,” [arXiv:hep-ph/9901312 [hep-ph]].
  * (36) D. Curtin, P. Meade and H. Ramani, “Thermal Resummation and Phase Transitions,” Eur. Phys. J. C 78, no.9, 787 (2018) [arXiv:1612.00466 [hep-ph]].
  * (37) A. Masoumi, K. D. Olum and B. Shlaer, “Efficient numerical solution to vacuum decay with many fields,” JCAP 01, 051 (2017) [arXiv:1610.06594 [gr-qc]].
  * (38) A. D. Linde, “Fate of the False Vacuum at Finite Temperature: Theory and Applications,” Phys. Lett. B 100, 37-40 (1981).
  * (39) J. M. Cline, M. Joyce and K. Kainulainen, “Supersymmetric electroweak baryogenesis,” JHEP 07, 018 (2000) [arXiv:hep-ph/0006119 [hep-ph]].
  * (40) L. Fromme and S. J. Huber, “Top transport in electroweak baryogenesis,” JHEP 03, 049 (2007) [arXiv:hep-ph/0604159 [hep-ph]].
  * (41) S. Bruggisser, T. Konstandin and G. Servant, “CP-violation for Electroweak Baryogenesis from Dynamical CKM Matrix,” JCAP 11, 034 (2017) [arXiv:1706.08534 [hep-ph]].
  * (42) N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020) [arXiv:1807.06209 [astro-ph.CO]].
  * (43) Y. B. Zeldovich, I. Y. Kobzarev and L. B. Okun, “Cosmological Consequences of the Spontaneous Breakdown of Discrete Symmetry,” Zh. Eksp. Teor. Fiz. 67, 3-11 (1974) SLAC-TRANS-0165.
  * (44) A. Mazumdar and G. White, “Review of cosmic phase transitions: their significance and experimental signatures,” Rept. Prog. Phys. 82, no.7, 076901 (2019) [arXiv:1811.01948 [hep-ph]].
  * (45) M. Ahmadvand and K. Bitaghsir Fadafan, “Gravitational waves generated from the cosmological QCD phase transition within AdS/QCD,” Phys. Lett. B 772, 747-751 (2017) [arXiv:1703.02801 [hep-th]].
  * (46) M. Ahmadvand and K. Bitaghsir Fadafan, “The cosmic QCD phase transition with dense matter and its gravitational waves from holography,” Phys. Lett. B 779, 1-8 (2018) [arXiv:1707.05068 [hep-th]].
  * (47) H. Abedi, M. Ahmadvand and S. S. Gousheh, “Electroweak phase transition in the presence of hypermagnetic field and the generation of gravitational waves,” [arXiv:1901.05912 [hep-ph]].
  * (48) M. Ahmadvand, K. Bitaghsir Fadafan and S. Rezapour, “Gravitational waves of a first-order QCD phase transition at finite coupling from holography,” [arXiv:2006.04265 [hep-th]].
  * (49) A. Kosowsky, M. S. Turner and R. Watkins, “Gravitational waves from first order cosmological phase transitions,” Phys. Rev. Lett. 69, 2026-2029 (1992).
  * (50) M. Kamionkowski, A. Kosowsky and M. S. Turner, “Gravitational radiation from first order phase transitions,” Phys. Rev. D 49, 2837 (1994) [astro-ph/9310044].
  * (51) A. Kosowsky, A. Mack and T. Kahniashvili, “Gravitational radiation from cosmological turbulence,” Phys. Rev. D 66, 024030 (2002) [arXiv:astro-ph/0111483 [astro-ph]].
  * (52) 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, no. 12, 123009 (2015) [arXiv:1504.03291 [astro-ph.CO]].
  * (53) C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau, P. Schwaller and G. Servant, et al. “Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions,” JCAP 04, 001 (2016) [arXiv:1512.06239 [astro-ph.CO]].
  * (54) 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, 024 (2009) [arXiv:0909.0622 [astro-ph.CO]].
  * (55) J. R. Espinosa, T. Konstandin, J. M. No and G. Servant, “Energy Budget of Cosmological First-order Phase Transitions,” JCAP 06, 028 (2010) [arXiv:1004.4187 [hep-ph]].
  * (56) S. Hild, S. Chelkowski and A. Freise, “Pushing towards the ET sensitivity using ’conventional’ technology,” [arXiv:0810.0604 [gr-qc]].
  * (57) B. Sathyaprakash, M. Abernathy, F. Acernese, P. Ajith, B. Allen, P. Amaro-Seoane, N. Andersson, S. Aoudia, K. Arun and P. Astone, et al. “Scientific Objectives of Einstein Telescope,” Class. Quant. Grav. 29, 124013 (2012) [erratum: Class. Quant. Grav. 30, 079501 (2013)] [arXiv:1206.0331 [gr-qc]].
  * (58) L. Di Luzio, B. Gavela, P. Quilez and A. Ringwald, “An even lighter QCD axion,” JHEP 05, 184 (2021) [arXiv:2102.00012 [hep-ph]].
  * (59) H. Kudoh, A. Taruya, T. Hiramatsu and Y. Himemoto, “Detecting a gravitational-wave background with next-generation space interferometers,” Phys. Rev. D 73, 064006 (2006) [arXiv:gr-qc/0511145 [gr-qc]].
  * (60) K. Yagi and N. Seto, “Detector configuration of DECIGO/BBO and identification of cosmological neutron-star binaries,” Phys. Rev. D 83, 044011 (2011) [erratum: Phys. Rev. D 95, no.10, 109901 (2017)] [arXiv:1101.3940 [astro-ph.CO]].
  * (61) M. C. Smith, G. R. Ruchti, A. Helmi, R. F. G. Wyse, J. P. Fulbright, K. C. Freeman, J. F. Navarro, G. M. Seabroke, M. Steinmetz and M. Williams, et al. “The RAVE Survey: Constraining the Local Galactic Escape Speed,” Mon. Not. Roy. Astron. Soc. 379, 755-772 (2007) [arXiv:astro-ph/0611671 [astro-ph]].
  * (62) K. Kannike, M. Raidal, H. Veermäe, A. Strumia and D. Teresi, “Dark Matter and the XENON1T electron recoil excess,” Phys. Rev. D 102, no.9, 095002 (2020) [arXiv:2006.10735 [hep-ph]].
  * (63) E. Aprile et al. [XENON], “Excess electronic recoil events in XENON1T,” Phys. Rev. D 102, no.7, 072004 (2020) [arXiv:2006.09721 [hep-ex]].
  * (64) S. Shakeri, F. Hajkarim and S. S. Xue, “Shedding New Light on Sterile Neutrinos from XENON1T Experiment,” JHEP 12, 194 (2020) [arXiv:2008.05029 [hep-ph]].
  * (65) C. A. Argüelles, A. Diaz, A. Kheirandish, A. Olivares-Del-Campo, I. Safa and A. C. Vincent, “Dark Matter Annihilation to Neutrinos,” [arXiv:1912.09486 [hep-ph]].
  * (66) M. Chianese, D. F. G. Fiorillo, R. Hajjar, G. Miele, S. Morisi and N. Saviano, “Heavy decaying dark matter at future neutrino radio telescopes,” JCAP 05, 074 (2021) [arXiv:2103.03254 [hep-ph]].