Singlet fermion dark matter and Dirac neutrinos from Peccei-Quinn symmetry

The Peccei-Quinn (PQ) mechanism not only acts as an explanation for the absence of strong CP violation but also can play a main role in the solution to other open questions in particle physics and cosmology. Here we consider a model that identifies the PQ symmetry as a common thread in the solution to the strong CP problem, the generation of radiative Dirac neutrino masses and the origin of a multicomponent dark sector. Specifically, scotogenic neutrino masses arise at one loop level with the lightest fermionic mediator field acting as the second dark matter (DM) candidate thanks to the residual $Z_2$ symmetry resulting from the PQ symmetry breaking. We perform a phenomenological analysis addressing the constraints coming from the direct searches of DM, neutrino oscillation data and charged lepton flavor violating (LFV) processes. We find that the model can be partially probed in future facilities searching for WIMPs and axions, and accommodates rates for rare leptonic decays that are within the expected sensitivity of upcoming LFV experiments.

a common thread in the solution to the strong CP problem, the generation of radiative Dirac neutrino masses and the origin of a multicomponent dark sector. Specifically, scotogenic neutrino masses arise at one loop level with the lightest fermionic mediator field acting as the second dark matter (DM) candidate thanks to the residual Z 2 symmetry resulting from the PQ symmetry breaking. We perform a phenomenological analysis addressing the constraints coming from the direct searches of DM, neutrino oscillation data and charged lepton flavor violating (LFV) processes. We find that the model can be partially probed in future facilities searching for WIMPs and axions, and accommodates rates for rare leptonic decays that are within the expected sensitivity of upcoming LFV experiments.

I. INTRODUCTION
The apparent non-observation of CP violation in the QCD Lagrangian represents one of the most active subjects in high energy physics, both theoretical and experimentally speaking. In the theory side, the absence of strong CP violation can be dynamically explained invoking the Peccei-Quinn (PQ) mechanism [1], which considers the spontaneous breaking of an anomalous global U (1) symmetry with the associated pseudo-Nambu-Goldstone boson, the (QCD) axion [2,3]. The axion itself turns to be a promising candidate for making up the dark matter (DM) of the Universe thanks to a variety of production mechanisms [4], for instance via the vacuum misalignment mechanism [5][6][7]. Besides this, it is remarkable that the physics behind the PQ mechanism can be also used to address other open questions in particle physics and cosmology such as neutrino masses [8][9][10][11][12][13][14][15], baryon asymmetry [16][17][18][19] and inflation [20][21][22][23].
The recent analysis [24] considering the PQ mechanism as the responsible for the massiveness of neutrinos revealed that it is also possible to consistently accommodate radiative Dirac neutrino masses 1 with a viable WIMP DM candidate, thus providing a set of multicomponent scotogenic models with Dirac neutrinos. Concretely, in these scenarios one-loop Dirac neutrino masses are generated through the d = 5 effective operatorLHN R S [40,41] once the axionic field S develops a vacuum expectation value, with the contributions arising from the tree-level realizations of such an operator being forbidden thanks to the charge assignment. As a further consequence of the PQ symmetry, the residual discrete symmetry that is left over renders stable the lightest particle mediating the neutrino masses, and since such a particle must be electrically neutral, it turns out that the setup also accommodates a second DM species [26,29,[42][43][44][45].
In this work we perform a phenomenological analysis of the T3-1-A-I model introduced in Ref. [24]. In order to determine the viable parameter space of the model we take into account the constraints coming from direct detection experiments, lepton flavor violation (LFV) processes, DM relic density and neutrino physics. We find that for a wide and typical range of the parameter values, the model easily satisfies these constraints and, additionally, future experiments will be able to test a considerable portion of the parameter space.
The layout of this paper is organized as follows. The main features of the model are presented in Section II. In Section III we determine the elastic scattering cross section between the WIMP particle and nucleons, and estimate the expected number of events in current and future direct detection experiments. Section IV is dedicated to a numerical analysis addressing the DM and LFV phenomenology. Finally, we conclude in Section V.

II. THE MODEL
As usual in models with massive Dirac neutrinos, this model extends the SM with three singlet Weyl fermions N Rβ (β = 1, 2, 3) −the right-handed partners of the SM neutrinos. The one-loop neutrino mass generation additionally demands [24] the introduction of one SU (2) L doublet scalar H 2 , one singlet scalar Φ and two singlet Dirac fermions ψ i (i = 1, 2). As the last piece we have a chiral exotic down-type quark D which is added in order to realize the hadronic KSVZ-type axion model [46,47]. In Table I  The relevant part of the scalar potential can be expressed as where the coupling constants associated to the quartic interaction terms |Φ| 2 |S| 2 , |H 2 | 2 |S| 2 , |H 2 | 2 |Φ| 2 and |H 1 | 2 |S| 2 have been assumed to be small. Since the (H † 1 H 2 ) 2 + (H † 2 H 1 ) 2 term is forbidden, it follows that the neutral component of H 2 remains as a complex field and does not get splitted into a CP-even and a CP-odd field. Nevertheless, it does get mixed with Φ through the term proportional to λ (since both scalar fields do not acquire a nonzero vacuum expectation value). We parametrize the scalar fields as where ρ stands for the radial component of the field S whose mass is set by the scale of the PQ symmetry breaking v S , whereas the angular part of S corresponds to the QCD axion a, h is the SM Higgs boson and v SM = 246.22 GeV.
In the basis (H 0 , Φ), the mass matrix for the Z 2 -odd neutral scalars leads to the mass eigenstates S 1,2 (both with two degrees of freedom since they are complex) via the transformation where we have defined C H1 = cos ϕ, C H2 = sin ϕ, C Φ1 = − sin ϕ and C Φ2 = cos ϕ. The mixing angle is defined through the expression 3 being m S 1,2 the eigenvalues of M 0 . In our analysis we will take m S 1 < m S 2 , which implies that for ϕ → 0 the heavier scalar is mainly singlet, whereas for ϕ → ± π/2 the heavier one is mainly doublet. For the charged scalar H ± we find that its mass is given by m 2 H ± = µ 2 2 + 1 2 λ 3 v 2 SM , just as happens in the inert doublet model.

A. Neutrino masses and charged LFV
The new Yukawa interactions involving the SM leptons are given by where y iβ and h βi are 2 × 3 and 3 × 2 Yukawa matrices, respectively. After the spontaneous symmetry breaking, scotogenic Dirac neutrino masses are generated at one loop [40,41,[48][49][50][51][52][53][54][55][56][57] as 3 A tiny value for λ is not only necessary to reproduce the observed neutrino phenomenology but also to have the complex scalars S1,2 at or below the TeV mass scale. This is in consonance with the requirement of demanding a tiny value for the scalar couplings between the axion field S and the other scalar fields, as happens in most of the axion models.
is illustrated in Fig. 1. The effective mass matrix can be written as where and F (x) = x ln(x)/(x − 1). The Dirac neutrino mass matrix M ν can be diagonalized through with h 11 = h 12 = 0, whereas in the inverted hierarchy (IH) case with h 31 = h 32 = 0. Notice that one of the right-handed neutrinos becomes decoupled because we are considering a scenario with the minimal set of singlet fermions 4 .
Although Diracness of neutrino masses is compatible with the conservation of the total lepton number, family lepton number violation is unavoidable due to neutrino oscillations. In this model, LFV processes involving charged leptons such as α → β γ, α → 3 β and µ − e conversion in nuclei are induced at one-loop level, involving only y iβ Yukawa interactions and mediated by the H ± charged scalar and the ψ i neutral fermions. The decay rate for the α → β γ processes (see the top left panel of Fig. 2) neglecting the lepton masses at the final states is given by Feynman diagrams for the α → β γ and α → β¯ β β charged LFV processes present in the model.
Concerning the α → β¯ β β processes, there are two kinds of diagrams (see Fig. 2): the γand Zpenguin diagrams (top-right panel) and the box diagrams (bottom panel).
The contribution from the Higgs-penguin diagrams is suppressed for the first two charged leptons generations due to their small Yukawa couplings 5 . It follows that the α → β¯ β β processes contain four kind of contributions: the photonic monopole, photonic dipole, Z-penguin and box diagrams 6 . In contrast, the photonic dipole contribution is the only one present in the α → β γ processes.
Finally, the µ − e conversion diagrams are obtained when the pair of lepton lines attached to the photon and Z boson in the penguin diagrams (see top panels of Fig. 2) are replaced by a pair of light quark lines 7 . For the µ − e conversion in nuclei there are no box diagrams since the Z 2 -odd particles do not couple to quarks at tree level. Accordingly, the photonic non-dipole and dipole terms along with the Z−penguin one are the only terms that contribute to the µ − e conversion processes. In this work we calculate the rates for α → β¯ β β and µ − e through the chain SARAH [60,61], SPheno [62,63] and FlavorKit [64]. 5 The contribution of those processes involving tau leptons is not negligible, but the corresponding limits are less restrictive. 6 Notice that the Yukawa interactions also lead to neutrino three-body decays ν h → ν lνl ν l , with mν h > mν l , via a box diagram similar to the bottom panel in figure 2, with the charged leptons replaced by neutrinos and the charged scalar H ± by the neutral one H 0 . Since the decay rate for these processes [59] is proportional to the ratio m 5 ν h /m 4 S 1 , the expected lifetime is several orders of magnitude larger than the age of the Universe. 7 Higgs-penguin diagrams are again suppressed, in this case by the Yukawa couplings to light quarks.

B. Two-component dark matter
The natural DM candidate in models featuring a PQ symmetry is the axion itself since the associated energy density decreases as non-relativistic matter does [4,[65][66][67]. The amount of axion relic abundance depends on whether the PQ symmetry is broken before or after inflation.
On the one hand, if PQ symmetry is broken after the inflationary epoch the axion field would be randomly distributed over the whole range (0, 2πv S ), meaning that the initial misalignment angle θ a takes different values in different patches of the Universe resulting in the average θ a = π/ √ 3.
In this case, topological defects such as string axions and domain walls [68][69][70][71] also contribute to the axion relic abundance. On the other hand, when PQ symmetry is broken before the inflationary epoch and is not restored during the reheating phase the axion field is uniform over the observable Universe, meaning that the initial misalignment angle takes a single value in the interval (0, 2π].
For simplicity purposes, we assume that the reheating temperature after inflation is below the PQ symmetry breaking scale, in which case the axion abundance is settled to [6,72] It follows that the axion can be the main DM constituent if v S ∼ 10 12 GeV for a no fine-tuned θ a (that is θ a ≈ O(1)). Under this premise, the axion window becomes m a ∼ (1 − 10) µeV.
Nevertheless, the axion can give a subdominant contribution to the relic DM abundance for lower values of v S , thus allowing for a multicomponent DM scenario.
In addition to the axion, this model brings along with a second DM candidate since the remnant Z 2 symmetry renders stable the lightest particle charged under it. The case of S 1 being that candidate turns to be ruled out since direct detection searches have excluded models where the DM candidate has a direct coupling to the Z gauge boson. Therefore, ψ 1 ≡ ψ becomes the second viable DM candidate of the model. According to Eq. (7), ψ only interacts with the SM particles via the Yukawa interactions y and h, and since these must be non-negligible in order to explain the neutrino oscillation parameters ψ necessarily reaches thermal equilibrium with the SM plasma.
The ψ relic abundance is determined by the cross sections of the annihilation and co-annihilation processesψψ →¯ ,νν and H + ψ → γ, respectively. Let us to stress that the h interactions can actually take large values because they are not taking part of the LFV processes, which means that ψ may feature a large annihilation cross section. If the fermion DM and scalar mediator masses are assumed to be sufficiently non-degenerate, co-annihilations can be neglected, and thus the relic abundance simply depends on the ψ-annihilation cross section. In our numerical analysis, nonetheless, we use micrOMEGAs [73] in order to take into account all the relevant processes that contribute to the setting of the relic abundance of ψ.

III. DIRECT DETECTION OF FERMION DARK MATTER
Being a SM singlet that couples to leptons and Z 2 -odd scalars, ψ does not have tree-level interactions with the SM quarks. However, the interactions involving Z 2 -odd particles allow us to construct effective interactions at one-loop level between singlet fermion DM and quarks.
In the basis of mass eigenstates, the relevant interaction terms involved in the direct detection where a sum over repeated latin and greek indices is implied. The new coefficients in this expression are defined as The interplay of the above interactions with the gauge and scalar interactions lead to the effective one-loop couplings between ψ and the Higgs, photon and Z bosons as shown in Fig. 3.
The differential spin-independent cross section for the fermion DM particle being scattered by a target nucleus of mass m T and atomic and mass numbers A and Z can be expressed as [74] dσ Here v rel is the relative velocity between ψ and the nucleus, and E R denotes the recoil energy of the nucleus due to the interaction. The maximum value of The effective couplings C γ M (E) and C γ R correspond to the Wilson coefficients describing the interaction with the photon field as a result of the magnetic (electric) dipole moment and the charge radius of the fermion DM. In this model these quantities are given by where explicit expressions for the one-loop functions g M 1 and g R1 are reported in Appendix A. The vanishing of C γ E has its origin in the absence of a coupling between the right-handed electron and ψ as can be seen from Eq. (14). The interaction with nucleons is described by the effective scalar and vector couplings f where m N is the nucleon mass. On the other hand, the vector coupling can be expressed in terms of the ψ-quark vector couplings C q V as In this model, the ψ effective interactions with quarks and gluons can be cast as being C ψh(Z) the effective scalar (vector) coupling between ψ and the Higgs (Z) boson. These are given by As indicated above, a sum over repeated indices is implied and the definitions of the one-loop functionsg h1 , g Z1 ,g Z1 and g Z2 are reported in Appendix A. By last, the recoil-energy dependent nuclear form factor F 2 c (E R ) in Eq. (16) reads [77,78] F 2 c (E R ) = 3 where j 1 is the spherical Bessel function of the first kind, q = √ 2m T E R and R = c 2 + 7 3 π 2 a 2 − 5s 2 with c = (1.23A 1/3 − 0.6) fm, a = 0.52 fm and s = 0.9 fm.
In order to estimate the expected number of ψ-nuclei scattering events in a direct detection experiment like XENON1T [79], we calculate the differential event rate per unit of detector mass through [74] yield a recoil with energy E R , which can be determined from and f ⊕ (v rel ) stands for the DM velocity distribution measured with respect to the lab frame. With respect to the galactic frame, this distribution is assumed to follow a Maxwell-Boltzmann one, i.e.
where the maximum speed is equal to the galaxy escape velocity, v esc , and In this way, if v E is the velocity of the Earth with respect to the galactic frame, then 9 For the numerical analysis we took the values used by the XENON1T collaboration [80], namely, v 0 = 220 km/s, v esc = 544 km/s and v E = 232 km/s.
In the case of the direct detection experiment XENON1T, the number of expected events, N events , can be determined as [81] Here ω exp = 278.8 days×1. 30 [80]. ν(E R ), for its part, was calculated as 9 Given the functional dependence of dσ ψT dE R with v rel , the integrals must be calculated when determining dR dE R . Analytical expressions for these integrals can be found in Appendix C of Ref. [74].
From the most recent data reported by XENON1T, and with the aid of a Test Statistic (TS), we can obtain an upper bound for N events . Closely following Ref. [86], we take with and L BG ≡ L(0). It follows that by demanding TS(m ψ ) > 2.71, limits for N events are obtained at 90% CL. For N obs = 14 (number of observed events) and N BG = 7.36(61) (number of background events) [80], the expected number of events must fulfill N events 19.5 [74].

IV. NUMERICAL RESULTS
In order to study the fermion DM phenomenology and take into account the constraints associated with charged LFV processes, we have implemented the model in SARAH [60,61] to calculate, via SPheno [62,63] and FlavorKit [64], the flavor observables. In addition, we have used micrOMEGAs [73] to calculate the ψ relic abundance. We have performed a random scan over the relevant free parameters of the model as shown in Table II and assumed λ 2 = λ Φ = λ hΦ = λ S = 0.01. Moreover, the mass of the exotic quark, M D , has been set to ∼ 10 TeV along with f β = 0.1 in order to avoid  Fig. 4, we present N events as a function of the scalar mixing angle (| sin ϕ|). All the points satisfy the current LFV constraints and the current limit imposed by the XENON1T collaboration. The prospect limits expected by XENONnT are indicated by the horizontal dashed line [80]. Notice that a large fraction of the parameter space (red points) will be explored in the next years (those featuring a small mixing angle are beyond the projected sensitivity, although such a small mixing angle is favoured by the tiny neutrino masses). From Fig. 4 we also notice that for either ϕ → 0 or ϕ → ± π/2, the number of events decreases rapidly, being steeper for |ϕ| → π/2. This happens because when ϕ → 0, the effective coupling between ψ and the Higgs (C ψh ) becomes more suppressed than the case |ϕ| → π/2 (see Eq. (26)). However, the coupling between ψ and the Z boson (C ψZ ) remains unaltered for both In Fig. 5 we display the mass splitting ∆M ≡ m H ± − m ψ as a function of the scalar mixing angle. As in Fig. 4, the gray points are allowed by the current experimental searches, whereas the red ones represent the region that will be explored in the next years. Notice that in the case of maximal mixing (|ϕ| ∼ π/4), where the number of events reaches its maximum value, there are some points localized in the allowed region. For these points the mass splitting between the charged scalar and ψ is small, ∆M 50 GeV, which means that the co-annihilation processes are relevant.
Conversely, for either ϕ → 0 or ϕ → ± π/2 where the number of events drops sharply, the mass splitting can take any value in the allowed range determined in the scan.
We now turn to discuss the axion phenomenology. The contribution of ψ to the total DM relic density R ψ = Ω ψ /(Ω ψ + Ω a ) as a function of the axion mass is shown in Fig. 6. Each point reproduces the observed DM relic density Ωh 2 = 0.120 ± 0.001 at 3σ [112] and satisfies the direct detection constraints on ψ as well as the charged LFV bounds. The color code represents the corresponding axion-photon coupling, g aγ ∼ α 2πv S , which plays a main role in axion searches through helioscope and haloscope experiments [113][114][115]. For v S 10 10 GeV the main contribution FIG. 5. The mass splitting ∆M ≡ m H ± − m ψ as a function of the scalar mixing angle. All the points satisfy the current limit imposed XENON1T and the red ones are those that will be explored by XENONnT.
FIG. 6. ψ contribution to the total DM relic density as a function of the axion mass m a and the axionphoton coupling, g aγ ∼ α 2πv S , for different values of the initial misalignment angle θ a . The vertical dashed lines represent the regions that will be explored by the cavity haloscope experiments ADMX [104,105], CULTASK [106,107] and MADMAX [108][109][110][111].
to the DM relic abundance comes from the fermion DM candidate, with the corresponding axion mass window laying outside the experimental searches. However, for increasing values of v S the mixed fermion-axion DM scenario becomes more relevant and the axion can account for a fraction or the whole of the DM relic abundance. In this case, a large fraction of the axion mass window (with the axion-photon coupling taking values in the range g aγ ∼ (5×10 −16 −10 −13 ) GeV −1 ) can be explored by several haloscope experiments [114,115]: ADMX for m a ∼ (2.5 − 3.5) µeV, CULTASK for m a ∼ (3.5 − 12) µeV and MADMAX for m a ∼ (0.04 − 0.4) meV. Let us notice, however, that some regions are beyond the reach of the projected sensitivity of the experiments. Nonetheless, by enlarging the particle content or changing the PQ charge assignment on the current fields of the model, the chiral anomaly coefficient in the g aγ coupling can be modified in such a way that the entire region planned to search QCD axions in KSVZ and DFSZ models becomes experimentally accessible [114][115][116].
Regarding the charged LFV processes, we focus in the observables involving muons in the initial state. In Fig. 7 are displayed the branching ratios B(µ → eγ) (cyan points), B(µ → 3e) (magenta points) and the rate for the µ → e conversion in titanium R µe (blue points) as a function of the Yukawa coupling y 21 . The left (right) panel stands for the normal (inverted) neutrino mass hierarchy, with the dotted, solid and dashed horizontal lines representing the projected sensitivity of LFV experiments for B(µ → eγ), B(µ → 3e) and R µe , respectively. It follows that for both neutrino mass hierarchies the current bounds demand |y 21 | 0.1 whereas the future searches will explore values as low as 0.01, with the conversion in nuclei being the most relevant process (we have found similar results for the other Yukawa couplings y iβ ). In addition to this, we notice that the observables associated with the lepton conversion in nuclei and the process with three electrons in the final state are strongly correlated. This is because the conversion process in nuclei does not involve box diagrams whereas the corresponding contribution in µ → 3e becomes suppressed by a factor |y iβ | 4 and therefore the penguin diagrams give the dominant contribution for both charged LFV processes. This strong correlation, along with the correlation from B(µ → eγ), is shown in Fig. 8. On the other hand, we observed that the Yukawa couplings h βi can take values along the whole range considered in the scan. This result along with the ones for y iβ discussed above translate to that the scalar coupling |λ | lies below 10 −7 in order to reproduce the observed neutrino mass scale (see Eq. (9)).
A final comment is in order concerning the exotic quark D. Since it couples to the SM sector through the Yukawa term, q L H 2 D R , it can decay into a S 1,2 scalar and a SM quark, thus avoiding potential issues that arise when an exotic quark is considered cosmologically stable [117]. Furthermore, such an exotic quark can be produced at colliders via quark and gluon fusion, leading to the specific signature of jets plus missing energy. From the analysis presented in Ref. [44], which deals with a scenario similar to the one studied here, the LHC searches for exotic quarks imply that M D 700 GeV, which is well below the value considered in our analysis.

V. CONCLUSIONS
In the class of models known as scotogenic models the generation of radiative neutrino masses is associated with the existence of a discrete symmetry that forbids the tree-level contribution and stabilizes the DM particle. On these lines, the PQ symmetry can also be invoked to simultaneously provide a solution to the strong CP problem, radiatively induce neutrino masses and stabilize the particles lying in the dark sector. In this work we considered a multicomponent scotogenic model with Dirac neutrinos where the dark sector is composed by the QCD axion and a SM singlet fermion, the latter stabilized by the PQ symmetry. We computed the expected number of fermion-DM nuclei scatterings in XENON1T and identified regions of the parameter space compatible with observed DM abundance, direct searches of singlet fermion and axion DM, neutrino oscillation data and the upper bounds on charged LFV processes. Furthermore, we find that for some choices of parameters both the singlet fermion and the axion present detection rates are within the expected sensitivity of XENONnT and haloscope experiments such as ADMX, CULTASK, and MADMAX, respectively, and that the LFV processes involving muons in the initial state may be probed in upcoming rare leptonic decay experiments.