Filtered asymmetric dark matter during the Peccei-Quinn phase transition

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.


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 [1]. 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 [2], though direct detection searches disfavor some of these scenarios and lead us to consider other possibilities such as non-thermal DM models [3] and scenarios in which DM is produced in the decay of thermally-decoupled heavy particles [4,5]. 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 [4,6].
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 [7,8]. 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 thermallyproduced DM and a possibility to produce DM masses above 100 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 PeV, the Peccei-Quinn (PQ) symmetry U(1) PQ is spontaneously broken according to QCD axion models, resolving the strong CP problem [9,10]. 1 Provided that the PQ PT is first order, the highly massive DM candidates can be naturally generated through the bubble filteringout 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 [13,14] 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 [15][16][17]. 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 [4,18]. 2 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 TeV [22][23][24][25]. 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.

The Model
Due to the anomalous U(1) axial symmetry of strong interactions and QCD vacuum structure, an effective CP-violating term associated with θ-vacuum is allowed in the Lagrangian. This term contributes to the neutron electric dipole moment, d n ≈ 3.6 × 10 −16θ e cm [26], which is experimentally constrained |d n | < 2.9 × 10 −26 e cm (90% C.L.) [27] and therebȳ θ 10 −10 . This poses the strong CP problem in that there is no reason in the SM whyθ should be very small.
An interesting solution to this problem is based on a global chiral U(1) PQ symmetry proposed by Peccei and Quinn [9]. In fact, at QCD scales, the interaction of a pseudo-scalar field, the axion a, is effectively added to the Lagrangian, (a/f a +θ)G G, where f a is the 1 Considering astrophysical constraints, the symmetry breaking scale is bounded between few 10 8 GeV and few 10 17 GeV [11,12]. 2 For baryogenesis scenarios fulfilled via lepton-number violation at a first order PT, see [17,[19][20][21].
axion decay constant, G µν is the gluon field strength and G denotes its dual. Therefore, the CP-violatingθ-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 U(1) 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 [28,29] which contain extra heavy quarks and PQ scalar fields, carrying the PQ charge, and DFSZ models [13,14] 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 2) Φ is the PQ scalar field singlet under the SM gauge symmetry group SU(2) L × U(1) Y , H u and H d denote two SU(2) L doublets, and the interaction of singlet RH neutrinos with Φ violates the lepton number by two units ∆L = 2 in case L(N R i ) = 1. 3 We also consider a discrete 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 where X X is the PQ charge of a given field. Moreover, imposing the orthogonality between PQ and corresponding hypercharge currents, −X Hu v 2 where v is the Electroweak (EW) vev, and choosing X Φ = 1, we may determine all charges through  Table 1. Charge assignment and group representation of the fields in the model v a v φ . 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 Φ, such that Φ = 1/ √ 2(φ + f a ) exp(ia/f a ) and the axion field is shifted by the PQ transformation.

First order Peccei-Quinn phase transition
In order to describe the PT, we first express the potential in terms of electrically neutral components, 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 φ and h d . Consequently, there are two minima: For this PT, we study the first direction. 4 To have such a minimum, the potential should satisfy the condition Moreover, we are interested in the case in which m N i > m φ in the broken phase, thus y i > 2λ φ , 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, where at zero temperature, the one-loop Coleman-Weinberg quantum correction can be written as [34] V Here we take Λ = 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 [35,36] and the useful form of these thermal functions in the high temperature limit is given by 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 where the thermal mass squared of the scalars are where λ t 1, g 0.42, g 0.12, and cos 2 θ w 0.77. We calculate the effective potential accordingly with some representative values of parameters, e.g., for y 1 = √ 2, As shown in Fig. (1), for f a = 100 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 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. (2.10).
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  By extremizing the action, we obtain the bounce equation and according to the following boundary condition, the equation can be numerically solved Using the any bubble code [37], as can be seen from Fig. (2), we obtain the bubble profile connecting the two phases. Bubbles nucleate when the bubble formation probability, proportional to exp (−S 3 (T )/T ), is of the order of one. Therefore, we can find the nucleation temperature, T n , by the following relation [38] where H * T 2 n /M pl is the Hubble parameter and M pl 2.43 × 10 18 GeV is the reduced Planck mass. As a result, using the obtained solution, from Eq. (2.18), we find the nucleation temperature of the PT, T n 9.19 PeV.
In the next section, relying on the first order PQ PT and the obtained parameters, we proceed to calculate the DM relic abundance.

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 in N (N ≡ N 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 in N , 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 F is the semiclassical force which is given by [16,39,40] and 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 p = 0 for helicity eigenstates, s = ±1 [16]. The last three terms of (3.2) are CP-violating sources during the PT and originated from the complex mass term 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 where A = 7.2, and B = 3.9. Moreover, we consider the CP violating phase as [41] θ(z) = arctan ∆θ 2 1 + tanh( z l w ) . (3.6) Using the ansatz f N = A(z, p z ) exp(−E p /T ) for the distribution function [7], we describe the deviation from equilibrium by 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 Integrating over p x and p y , for RH chirality the Liouville operator would be where because of the friction effects produced from reflected and penetrated (anti)particles, we used non-relativistic bubble wall velocities, so that where the integration over the collision term is given by [7] g N dp x dp y (2π where σ is the spin-averaged cross section of such processes N R (p p )+N R (q p ) → φ(k p )+φ(l p ) To solve the Boltzmann equation, we use the method explained in [7] and rewrite the equation as (3.11) and λ parameterizes a given curve on the (z, p z ) plane. Therefore, from Eqs. (3.8, 3.9), imposing appropriate boundary conditions, we first obtain z(λ) and p z (λ), then we find A ± 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   figure), the factor A + for DM particles is displayed. Since to obtain the observed relic abundance, it requires ∆θ 1 and hence since A + and A − are displayed almost the same, A + is only shown.
ensuring the equilibrium phase space distribution function far away from the bubble wall. Using the obtained curves of (z → −∞, p z > m in N ), for (anti)particles deep inside the bubble (z → ∞) with p z < 0 and similar dynamics, we set the boundary condition as [7] lim z→∞ A ± (p z ) = lim z→∞ A ± (−p z ). (3.13) Interpolating between different solutions obtained for several curves, we find A ± as displayed in Fig. (3). Finally, integrating over p z at deep inside the bubble, we obtain the net number density ∆n N /T 3 n in T n units, where ∆n N ≡ 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.
The DM relic abundance can be found in the asymmetric case by the following relation 5 where g * 0 = 3.9, T 0 0.235 meV, g * ∼ 110 and the present Hubble parameter H 0 = 100 h km sec −1 Mpc −1 . As a result, the observed relic abundance, h 2 Ω DM 0.12 [42], can be obtained for, e.g., y 1 = √ 2 or m in N 66.3 PeV and from the attained ∆n N /T 3 n which is determined by ∆θ ∼ 10 −12 .
Furthermore, analogous to the analytic approach of [7], estimated the DM number density proportional to its equilibrium abundance inside the bubble n eq /4 where here, considering the net number density as ∆n N ∼ δθ n eq , with T n 9.19 PeV and m in N 66.3 PeV, one can obtain the observed relic abundance for δθ 0.19 × 10 −12 .

Phenomenology
In this section we study phenomenological consequences of the model.

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θ term. Along this direction, the anomalous PQ current of the axial U (1) PQ symmetry is given as follows where g s is the SU(3) c strong coupling, e is the electric charge, F µν is the photon field strength, N and 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 n g is the number of SM fermion f generation, f a = v a /(2N ), m u and m d are u and d quark masses, respectively. The cosmological Domain Wall (DW) problem [43] can be avoided if N DW ≡ 2N = 1 [30]. As explained in Section 2, X Hu + X H d = X Φ . Therefore, according to Table 1

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 [44][45][46][47][48].
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 [49][50][51][52]. 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. 6 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, and v t (T n ) is the true vacuum at T n . There is a critical value of α, denoted by α ∞ , such that for α > α ∞ bubbles can run away [53], where n i (n i /2) is the number of degrees of freedom for boson (fermion) species, and ∆m 2 i 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 β H * = T n d dT Based on the analysis described in Section 2.1, we found T n 9.19 PeV. At this temperature we obtain α 0.55, α ∞ 0.6. In addition, from Eq. (4.9) and the obtained bounce action we find β/H * 246.
In the non-runaway case where α < α ∞ , dominant contributions to GWs come from sound waves and MHD turbulence, i.e., h 2 Ω(f ) h 2 Ω sw + h 2 Ω tu where [52,54] where the spectral shapes are given by [53], with h * = 16.5 × 10 −2 [Hz] T n PeV g * 100 1 6 . (4.14) The red-shifted peak frequencies in the spectral shapes are given by 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 c s = 1/ √ 3, v w = c s which is the case of transition from subsonic to supersonic deflagrations, and the large velocity limit, v w → 1, this efficiency factor is expressed respectively as [55] Moreover, for subsonic deflagrations v w c s , Jouguet detonations, and detonations with v w v J , the efficiency factor is given by the following fits [55] κ The fraction of plasma motion in turbulence, ε = κ tu /κ v , can be set as ε = 0.05 [52], so that the main source contributing to the GW energy density from the plasma motion is the source of sound waves, κ sw = (1 − ε)κ v . Eventually, substituting key parameters in Eqs. (4.10, 4.11), 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) [56,57]  For v w = 0.9, we use Eq. (4.22) for the efficiency factor. The dashed curve is the approximate prospective sensitivity curve of ET [56,57]. We show that the GW signals may be within the reach of this planned GW detector.

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 ∼ 1 − 100 TeV [25].
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 Z N mirror symmetry with N mirror worlds, so that with additional confining sectors the axion potential would change, (a/f a +θ eff )G k G k [22][23][24][25]. Consequently, in this case the axion mass would be [25] where m π and f π are the pion mass and its decay constant in the mirror sectors and z = m u /m d . For example taking z = z ∼ 0.5, m π = 135 GeV, f π = 93 GeV and f a = 100 TeV, the axion mass will be m a ∼ 60 MeV. Connecting the mirror worlds by the axion [58], in the UV axion model the PQ scalar can be transformed as φ → exp(2πi/N)φ. Thus, considering f a = 100 TeV, the PT parameters can be retained at T n 9.19 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 [59,60].
We leave further investigations on these models for another future work.

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 → N e scattering (Fig(6)) and explore the electron recoil energy. DMs move at non-relativistic velocities, v D ∼ 10 −3 , in typical DM halos [61]. Considering these speeds for the PeV-scale DM, one can calculate the transfer recoil energy as follows [62] E hence in the lab frame for v e = 0 and m e 0.5 MeV, the electron recoil energy would be E r eV. Therefore, such non-relativistic massive DMs cannot explain the excess events in a 1 − 7 keV recoil energy reported by the XENON1T experiment [63]. 7 However, the heavy DMs can be considered as a source for producing ultra-high energy neutrinos [65]. To show this process according to the model (Fig. (6)), we consider the annihilation of DM to SM neutrinos N (p) N (q) → v l (k) v l (g). In the center of mass frame, for the non-relativistic DM, v D ∼ 10 −3 , we have |p| = |q| ∼ m N v D and hence the energy of incoming particles E ∼ m N . Also, for m v l ∼ 0 and due to the energy and momentum conservation, E v l ∼ |k| ∼ m N . Finally, calculating the amplitude of the process, M, we can find the cross section as follows where Y ≡ y 1 y α y 2 Nα and we take m N = y 1 f a / √ 2 m φ . Therefore, as obtained from Eq. (4.25), for E v l ∼ m N = 1 PeV, the cross section is compatible with the bound on DM annihilation processes [66].

Conclusion
The bubble filtration mechanism provides a setup which allows DM masses above 100 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 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 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 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.