Fully Testable Axion Dark Matter within a Minimal SU (5) GUT

: We present a minimal Grand Uniﬁed Theory model, based on SU (5) gauge symmetry and a global U (1) Peccei-Quinn symmetry, that predicts the existence of an ultralight axion dark matter within a narrow mass range of m a ∈ [0 . 1 , 4 . 7] neV. This mass window is determined through an interplay between gauge coupling uniﬁcation constraints, partial proton decay lifetime limits, and the need to reproduce experimentally observed fermion mass spectrum. The entire parameter space of the proposed model will be probed through a synergy between several low-energy experiments that look for proton decay (Hyper-Kamiokande), axion dark matter through axion-photon coupling (ABRACADABRA and DMRadio-GUT), and nucleon electric dipole moments (CASPEr Electric).

A Renormalization group running of the gauge couplings 27 1 Introduction The Standard Model (SM) of elementary particle physics has performed exquisitely in explaining a multitude of experimental observations.There are, however, several important questions that evidently require physics beyond the SM in order to be fully addressed.For example, one of the most important discoveries in particle physics is the observation of nonzero neutrino masses, whereas neutrinos are strictly massless within the SM framework.Furthermore, it is well established that approximately 26% of the total energy density of the universe is in the form of the so-called dark matter that cannot be of the SM origin.This is especially puzzling as the stable SM matter only represents about 5% of the energy density of the universe.Also, the strong CP problem -why the QCD θ parameter takes the value 10 −10 or less -is still an open issue within the SM.
It might be that all these issues are related.In fact, the unified gauge theory [1][2][3][4][5][6] formulation of the elementary particle interactions is a very popular and successful tool for tackling the aforementioned shortcomings of the SM.The simplest possible scenario, among various possible choices of the Grand Unified Theory (GUT) groups, is the Georgi-Glashow model [3] that embeds the entire SM gauge group within an SU (5).In that construction, one 5-dimensional and one 10-dimensional representation of SU (5) comprise all the fermions of a single SM family.The SU (5) symmetry is broken down to the SM gauge group when a real Higgs in the adjoint representation acquires a vacuum expectation value (VEV).The SM symmetry is subsequently broken to SU (3) × U (1) em by the VEV of the SM Higgs doublet that resides within a fundamental representation.The Georgi-Glashow model, however, is incomplete since (i) it fails to achieve gauge coupling unification, (ii) it predicts wrong mass relations between down-type quarks and charged leptons, and (iii) neutrinos remain massless as in the SM.On top of that, the Georgi-Glashow model does not address the strong CP problem, nor does it include a dark matter candidate.
The most compelling new physics resolution of the strong CP problem is given in terms of the Peccei-Quinn (PQ) symmetry [7,8].In the PQ framework, a global U (1) PQ symmetry is spontaneously broken by a complex scalar leading to a nearly massless pseudoscalar particle [9][10][11][12][13][14], namely the "axion", which can, in turn, serve as a cold dark matter candidate [15][16][17].Intriguingly, as first shown in Ref. [18], the axion can be embedded within the scalar representation that breaks the GUT symmetry.The model presented in Ref. [18] did not, however, address several important GUT issues, such as neutrino mass generation and gauge coupling unification.For a sample of models that pursue this particular approach, but with a more realistic agenda, see Refs.[19][20][21][22].See also Ref. [23] for a discussion of other light axion-like particles and their connections to GUTs.
Our primary interest in this manuscript is to combine the PQ symmetry with a simple, yet realistic, SU (5) GUT scenario [24,25] and to investigate the main predictions of such a setup.The SU (5) proposal [24,25] in question extends the particle content of the Georgi-Glashow model by a 35-dimensional Higgs representation and a 15-dimensional vectorlike fermion representation.Remarkably, within that scenario, the observed mismatch between the down-type quarks and charged leptons is intrinsically connected to the neutrino mass generation.More specifically, the difference between the down-type quark and charged lepton mass matrices is given by a rank-one matrix.This stipulates that the down-type quarks and charged leptons have similar, yet, different masses, in accordance with experimental observations.The neutrino mass matrix, on the other hand, is made out of a sum of two rank-one matrices that are transpose of each other.This, in turn, dictates that one of the neutrinos is strictly a massless particle.Moreover, since the model relates these three rank-one matrices, the neutrino masses consequentially mirror the mismatch between the down-type quark and charged lepton masses and are thus of the normal hierarchy.
We extend the minimal realistic SU (5) proposal [24,25] with a PQ symmetry to address the strong CP problem as well as the origin of dark matter and show that such a simple extension still preserves the most prominent features of the original model.Our detailed study reveals that the proposed setup is highly predictive, and that the entire parameter space of the theory will be fully tested in the near future through a combination of several experimental efforts.These comprise the proton decay experiment Hyper-Kamiokande as well as the axion dark matter experiments ABRACADABRA, DMRadio-GUT, and CASPEr Electric.
The manuscript is organized as follows.In Sec. 2 we introduce the particle content and symmetries of the model.The details of the PQ symmetry implementation and the nature of the axion dark matter are discussed in detail in Sec. 3. A numerical study of the model is performed in Sec. 4, where we also present the most relevant experimental predictions.We briefly conclude in Sec. 5.
for convenience, refer to a given irrep/multiplet by using either its dimensionality with respect to the appropriate gauge group or the associated symbol.
Beside the non-trivial assignment under the Lorentz symmetry, the aforementioned SU (5) irreps carry the PQ U (1) PQ charges that are presented in Table II.
Before we write down and discuss relevant parts of the model Lagrangian, we briefly justify the proposed particle content.
• 24 H breaks the SU (5) × U (1) PQ symmetry.It furthermore provides axion dark matter (DM), helps to generate unification of the SM gauge coupling constants, and facilitates a process of creation of the experimentally observed mismatch between the down-type quark and charged lepton masses.
• 5 H 1 and 5 H 2 jointly break the SM gauge symmetry down to SU (3)×U (1) em .5 H 2 also provides the up-type quark masses through its vacuum expectation value (VEV), whereas 5 H 1 and 5 H 2 together play an indispensable role in three different mechanisms that create phenomenologically viable masses for the downtype quarks, charged leptons, and neutrinos.
• 35 H is essential for neutrino mass generation.It also helps to provide the gauge coupling unification at scales compatible with the existing limits on partial proton decay lifetimes.
• 15 F and 15 F participate in the neutrino mass generation mechanism.In addition to that, these SU (5) irreps are, together with 24 H and 5 H 1 , instrumental in addressing the observed mismatch between the down-type quark and charged lepton masses.

Scalar sector
There are several parts of the scalar sector of the model that need to be discussed in detail.The SU (5) × U (1) PQ symmetry breaking is due to The VEV of φ α β that does the SU (5) symmetry breaking reads where we assume that the VEV of the electrically neutral component of the SU (2) triplet φ 1 (∈ 24 H ) is negligible.The squares of masses of multiplets in 24 H , as generated via Eqs.(2.1) and (2.2), are (2.10) These results are summarized in Table III for convenience.We again emphasize that 24 H also breaks the PQ symmetry while we currently discuss solely the SU (5) symmetry breaking.(Hence the omission of an overall phase in Eq. (2.2).The exact role of that phase will be discussed in Sec.Table III: Mass-squared spectrum of a complex irrep 24 H ≡ φ.
The potential given by Eq. (2.1) dictates that the imaginary part of φ 0 (∈ 24 H ) is massless.In fact, the axion is mostly composed of that particular state, as we show later on.The real components of φ 3 (∈ 24 H ) and φ 3 (∈ 24 H ), on the other hand, provide the necessary degrees of freedom for the proton decay mediating gauge bosons in 24 V to obtain a mass M GUT , where Here, M GUT is also the scale of gauge coupling unification, and α GUT is the corresponding SU (5) gauge coupling.
The scalar fields in the fundamental irreps of SU (5) couple via where we suppress SU (5) indices.The doublet-triplet spitting, i.e., breaking of the mass degeneracy between Ξ a and ω a multiplets, is accomplished via the following additional terms in the scalar potential: The mass-squared matrices of ω a and Ξ a multiplets, in the Λ 1 -Λ 2 basis, are . (2.15) Clearly, the required doublet-triplet splitting can be obtained by an appropriate choice of the model parameters.The linear combinations of ω 1 and ω 2 will consequently yield mass eigenstates we denote T 1 and T 2 in the rest of the manuscript.Also, Ξ 1 and Ξ 2 will produce mass eigenstates H 1 and H 2 , where H 1 is identified with the SM Higgs with 125 GeV mass.Finally, the VEVs of 5 Ha that break The lepton number conservation is violated via a single term in the Lagrangian that reads The neutrino masses will thus be directly proportional to the dimensionless parameter λ of Eq. (2.16).
The masses of the SM gauge group multiplets in 35 H are determined by the following SU (5) contractions (2.17) The contractions of Eq. (2.17) yield (2.21) These, in turn, produce a single mass-squared relation that reads The mass spectrum given in Table III and the mass relation presented in Eq. (2.22) are necessary input for the gauge coupling unification analysis.

Fermion sector
The Yukawa sector of the model is where the PQ charge assignment of Table II and the SU (5) indices are all implicitly understood.The Yukawa matrix elements of the model are and y, where we have used the freedom to rotate irreps in the SU (5) group space to reach this particular Yukawa coupling basis.The model accordingly has nineteen real parameters and fifteen phases in the Yukawa sector to accommodate all of the masses and mixing parameters of the SM fermions as well as the masses of fermions in the 15 F -15 F vectorlike pair.
The PQ charge assignment forbids a bare-mass term for the 15 F -15 F pair.The masses of the associated SM gauge group multiplets are thus generated solely through the last term of Eq. (2.23), which reads where the overall phase of 24 H , once again, is not shown for simplicity.We subsequently define It is important to point out that, apart from different Clebsch-Gordan coefficients, all submultiplets within 15 F have a common mass scale.(Even though Σ 1 and Σ 3 mix with the fermions in 5 F i and 10 F i , this does not affect equalities in Eqs.(2.25) and (2.26).)We will show, later on, that the product yv φ is rather constrained by a requirement for the model to simultaneously generate large enough unification and neutrino mass scales.
The masses of the SM fermions are obtained after the breaking of the SM gauge group down to SU (3) × U (1) em as follows.The down-type quark sector 4 × 4 mass matrix can be written as where we introduce v φ = − 1 4 5 3 v φ .This matrix can be transformed into a blockdiagonal form comprising a 3 × 3 part denoted M d and a mass parameter M H as follows where unitary matrices X and Y take the form with ) Here, .This is possible due to the fact that Σ 3 ∈ 15 F and Q i ∈ 10 F i transform in the exact same way under the SM gauge group [26].
The charged fermion mass matrices of the model can be succinctly written as ) where GeV.We note the two most prominent features of the charged fermion sector.First, M u can be treated as a symmetric matrix in the flavor space.Second, a mismatch between the charged lepton and down-type quark mass matrices is proportional to a rank-one matrix Y c Y a .We again stress that we work in the basis where , where i = 1, 2, 3, are the masses of the SM charged leptons.
The neutrino mass in this model is generated by utilizing the Yukawa couplings Y a and Y b that appear in Eq. (2.23) and the lepton number violating term of Eq. (2.16).Completion of the neutrino mass loop requires, in addition to the SM fields, the presence of (1, 3, 1) + (1, 3, −1)(⊂ 15 F + 15 F ) vectorlike fermions and the scalar quadruplet (1, 4, −3/2)(⊂ 35 H ). The corresponding Feynman diagram illustrating the neutrino mass generation mechanism is shown in Fig. 1.This particular one-loop mechanism to generate neutrino masses has been introduced in Ref. [27,28].
The neutrino mass matrix elements (M ν ) ij , at the leading order, read where m 2 and m 3 are neutrino mass eigenstates and N is a unitary matrix.Note that one of the neutrinos is a strictly massless particle due to the fact that M ν is constructed out of two rank-one matrices with elements Y a i Y b j and Y b i Y a j .This is accordingly encoded in the right-hand side of Eq. (2.37).
Since the charged lepton mass matrix in Eq. (2.36) is already in a diagonal form, we can write that where V PMNS is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) unitary mixing matrix, that is defined as ).Here we use the PDG parametrization [29] for the R 23 , U 13 , and R 12 matrices.Note that there is only one Majorana phase β ν appearing in Q due to the fact that one of the neutrinos is massless.
One especially convenient feature of the neutrino sector is that the matrices Y a and Y b can be expressed in terms of the PMNS matrix parameters and phases η ν i , i = 1, 2, 3. Using the parametrization mentioned in Refs.[30,31] we can write the two Yukawa coupling vectors Y a and Y b as where N ij denotes the ij-th element of the unitary matrix N , r 2 = m 2 /m 0 , and r 3 = m 3 /m 0 .Moreover, ξ is a dimensionless scaling parameter that needs to be introduced if one is to accurately scan over all possible phenomenologically viable entries in Y a and Y b that accommodate experimental observables in the neutrino sector with utmost certainty.Eq. (2.39) is applicable solely to the normal neutrino mass hierarchy scenario since that is one of the model predictions, as we will discuss later.
3 Peccei-Quinn symmetry and axion dark matter We discuss the implementation of the PQ symmetry within our setup and elaborate on the model's main ingredients and experimental detection prospects in the following.
In the "invisible axion" models [11][12][13][14] the PQ symmetry is broken by a scalar field that carries a non-trivial PQ charge, where the scalar is a singlet under the SM.We embed this scalar within the 24-dimensional Higgs irrep that is charged under the U (1) PQ symmetry, as shown in Table II.Consequently, our setup unifies the GUT and PQ breaking scales.The VEV of 24 H ≡ φ α β can be written as [18] where the pseudoscalar part, i.e., field a φ (x), essentially remains massless, whereas the radial mode acquires a mass of the order of the GUT scale while the global U (1) PQ symmetry is spontaneously broken with order parameter v φ .To correctly identify the massless axion, one also needs to include all other Higgses that carry PQ charges and participate in symmetry breaking.The non-Hermitian operators that are responsible for the breaking of the rephasing symmetry of the three scalar fields are given by the terms in the second line of Eq. (2.13).The VEVs of neutral components of the SU (2) doublets can be re-written as where we take all VEVs to be real, and, as mentioned before, we neglect the VEV of the SU (2) triplet in 24 H .With these assumptions, the axion field is identified as [42], where x i denotes the PQ charge of the corresponding i-th scalar (and 16 GeV and v Λa ∼ 10 2 GeV, the axion mostly resides in 24 H with a ≈ a φ .
The axion field must also be orthogonal to the Goldstone field eaten up by the Z-boson.This translates into the following condition which, in our benchmark charge assignment, fixes tan β = 1.Here, we do not present the expression of the SM Higgs mass eigenstate, which can be obtained via the diagonalization of the 4 × 4 mass matrix of the CP-even states.The heaviest one is expected to reside at the GUT scale, and the lightest one is the SM Higgs boson.Depending on the chosen hierarchy, the remaining two eigenstates -one coming from the triplet and the other from the pair of doublets -can live anywhere in between the electroweak and GUT scales.Now, performing a field-dependent axial transformation that is anomalous under QCD, the axion can be disentangled from the Yukawa interactions.This transformation generates the effective anomalous interactions of the following types: Here, G (F ) is the gluon (photon) field strength tensor, G ( F ) is its dual, and f a is the axion decay constant.The effective operator of the form aG G is the key to the PQ solution to the strong CP problem.Since these sub-multiplets carry color and electromagnetic charges, the PQ current has both QCD and electromagnetic anomalies, with the corresponding anomaly coefficients [43], where sums are taken over all fermions, which we generically denote by ψ.Using well-known formulas, we obtain |N | ≡ N = 13/2 and |E| ≡ Ê = 52/3 while the domain-wall number, which is relevant for cosmology, is N DW = 2 N = 13.Subsequently, we find the axion decay constant to be Since the decay constant is of the order of the GUT scale, i.e., f a ∼ M GUT , we refer to the axion as the "GUT axion".Once strong interactions confine, non-perturbative QCD effects generate a potential that gives rise to a tiny axion mass [44,45] m a = 5.7 neV 10 15 GeV f a = 5.7 neV This shows that the axion mass is predicted if the grand unification scale M GUT is known.We accordingly compute the predicted range of the GUT scale within our model in Sec. 4 by taking into account all relevant constraints.Since the non-observation of proton decay requires the GUT scale to be large, the axion mass is expected to be around the neV scale within our setup.An axion in this mass range is extremely weakly coupled to the SM particles due to an extremely large decay constant.Remarkably, an axion with neV mass can serve as an excellent dark matter candidate and can be searched for efficiently in direct detection experiments [50] hunting for ultra-light axions.
Next, we consider the most relevant axion couplings for experimental sensitivity.In the low-energy effective Lagrangian for the axion, it is sometimes convenient to eliminate the axion coupling to the gluons via a field-dependent axial transformation of the SM quarks.After making such a rotation, the axion coupling to the photons is given by [45], where the model-dependent quantity, apart from f a (see Eq. (3.9)), in our case, is given by Ê/ N = 8/3.In fact, the dark matter experiment ABRACADABRA [46] has a great potential to look for an axion dark matter in the mass range of interest.As Expected reach in the m a vs. g aγγ plane for the broadband (Broad) and resonant (Res) strategies of the ABRACADABRA (ABD) experiment [46].The blue line (that lies on the QCD line) corresponds to the prediction of our model.The projected 3 σ sensitivity of DMRadio-GUT [47,48] is also presented in the green shaded region.Furthermore, the expected theoretical reach using the optomechanical cavity method [49] is shown with solid black lines.See text for more details.
shown in Fig. 2, a major part of the parameter space of our theory will be probed by this dark matter direct detection experiment.Fig. 2 is obtained by varying the model parameters while imposing all relevant constraints.The details of our numerical procedure are relegated to Sec. 4.
Another axion dark matter experiment, the DMRadio-GUT [47,48], will also be sensitive in detecting axions with GUT scale decay constant f a (∼ 10 16 GeV).DMRadio-GUT will be far more sensitive compared to its previous two phases, i.e., DMRadio-50L and DMRadio-m 3 , since it will have a factor of three enhancement in the field and a factor of ten enhancement in volume relative to DMRadio-m 3 .The projected 3 σ sensitivity of DMRadio-GUT is also presented in Fig. 2 by a green shaded region, which will probe a significant portion of the parameter space.Yet another proposal utilizing an optomechanical cavity [49] filled with superfluid helium is shown to be highly promising in detecting ultra-light axion dark matter.This proposed experimental method, with a cavity size of order O(10 m) is expected to be sensitive to axion-photon couplings for axions with the GUT scale size decay constant.In Fig. 2, the corresponding theoretical reach is shown with solid black lines.The ABRACADABRA experiment will be sensitive to axion masses as low as m a ∼ 2 neV, whereas the sensitivity of DMRadio-GUT and optomechanical cavity is about m a ∼ 0.4 neV and m a ∼ 0.1 neV, respectively.A combination of all these axion dark matter experiments will eventually probe the entire parameter space of  The blue band (that lies on the QCD line) corresponds to the prediction of our model; see text for details.The shaded regions show the sensitivity projections of CASPEr Electric [51,52] in its various phases.Moreover, the ultimate sensitivity limit is given by the nuclear spin noise.
the proposed model.Intriguingly, ultra-light axion dark matter can also be efficiently searched for via oscillating nucleon electric dipole moments (EDM).As already stated, the QCD axion solves the strong CP problem by promoting the θ parameter into the dynamical axion field.Consequently, the effective θ angle gives rise to an EDM for nucleons sourced by the axion.Owing to the dynamical nature of the axion, this EDM will change in time, giving rise to unique signals.In the effective Lagrangian, the coupling of the axion to nucleon n takes the following form, The nucleon electric dipole moment generated through the above operator is given by d n = g aD a.The classical field that describes the axion field can be written as a = a 0 cos(m a t).The amplitude, a 0 , is determined from the local dark matter density, namely, ρ DM = 1 2 m 2 a a 2 0 , which assumes the axion comprises 100% of dark matter within our setup.The nucleon electric dipole moment is then determined by the dark matter energy density, d n = √ 2g aD √ ρ DM cos(m a t)/m a .Moreover, the nucleon electric dipole moment can also be expressed in terms of the axion decay constant.In terms of our model parameters, it can be re-written in the following form [53]: with roughly a 40% uncertainty [54], where the decay constant is given in Eq. (3.9).(See also Refs.[55][56][57].)The corresponding coupling as a function of the axion mass is shown in Fig. 3.As can be seen from this figure, excitingly, the CASPEr Electric [51,52] experiment alone will probe almost the entire parameter space of our model.The width of the band corresponds to the calculation uncertainty as mentioned before.Fig. 3 is also obtained by varying model parameters after one imposes all the relevant constraints.The exact details will be discussed later in the text.
Since the axion is ultra-light in our setup, it can constitute the entirety of the dark matter.It is important to point out that the breaking of the GUT symmetry to that of the SM gauge group SU (5) × U (1) PQ → SU (3) × SU (2) × U (1) leads to an overproduction of super-heavy monopoles that must be inflated away.As discussed above, spontaneous breaking of the PQ symmetry leads to N DW distinct degenerate vacua, giving rise to a domain-wall problem, which also requires dilution to be consistent with cosmology.Both of these problems, along with the horizon and flatness problems, can be elegantly solved via inflation taking place after the GUT symmetry breaking.We, however, do not specify the details of the inflationary dynamics, which is beyond the scope of this work.The amount of axion dark matter produced then depends on whether the PQ symmetry is restored or not after inflation.We assume that the U (1) PQ remains broken during inflation and is never restored afterwards.In such a scenario, the relic abundance of the axion dark matter is given by [58] Ωh 2 ∼ 0.12 5 neV m a 1.17 which shows that the initial value of θ i = a i /f a , where a i is the value of the axion field, needs to be somewhat smaller than unity to be consistent with the observed dark matter relic abundance Ωh 2 ∼ 0.12 ± 0.001 [59].Thus, for θ i ∼ 10 −2 , the axion is composed of all the dark matter.(μ) Moreover, since the partial proton lifetimes are proportional to the fourth power of the GUT scale, our model can be simultaneously probed with axion dark matter and proton decay experiments.

Unification
The renormalization group equations (RGEs) for the gauge couplings can, at the 2-loop level, be written as [60] Note that we neglect the effect of the Yukawa couplings Y a , Y b , and Y c on the running of the gauge couplings.
In order to investigate the viable part of the parameter space giving gauge coupling unification, we freely vary the masses of the fields φ Re , and H 2 between the TeV and the GUT scale, while taking into account the mass spectrum constraints presented in Sec. 2. Some of these states remain light to achieve high scale unification to be compatible with proton decay bounds.To get an understanding of how these states can be light, let us consider the masses of φ Re 1 and φ Im 8 , which are given by Eq. (2.5) and Eq.(2.8), respectively.It can be easily seen that by adjusting the relevant quartic couplings, they can, in principle, live in the low energies (similar arguments are applicable for the rest of the states that reside somewhat below the GUT scale; for details, see Sec. 2).Such adjustments, however, introduce additional fine-turning problems on top of the usual doublet-triplet splitting problem, which we accept.We also ensure that the scalar leptoquark mediated proton decay is sufficiently suppressed by varying the masses of T 1 and T 2 between 3 × 10 11 GeV and the GUT scale.The numerical fit is performed by running the gauge couplings at the 2-loop level from the GUT scale to the Z mass scale at which a χ 2 -function that we define later in detail is minimized.We use the low-scale values g 1 = 0.461425 +0.000044 −0.000043 , g 2 = 0.65184 +0.00018 −0.00017 , and g 3 = 1.2143 +0.0035 −0.0036 [61] as our input, where g i = √ 4πα i .To demonstrate that within our setup, the gauge couplings can indeed unify, Fig. 4 shows one possible particle mass spectrum giving exact gauge coupling unification that is in agreement with the current proton decay constraints and that yields correct neutrino mass scale via Eq.(2.37).

Proton decay
The formulae for the proton decay widths of various decay channels can be found in Refs.[62,63].For example, the decay width for the proton decay channel having a pion and a charged lepton in the final state is given by1 Here, m p = 0.9393 GeV and m π = 0.134 GeV denote the proton and pion masses, respectively, while α = 1, 2 with e + 1 ≡ e + and e + 2 ≡ µ + .The leading log renormalization of the dimension six operators is encoded via A L = 1.2 [65] Table IV: Present experimental bounds on the partial lifetimes τ p as well as future sensitivities for 10 years of runtime, both at 90% confidence level. where2 with γ L(R)i = (23(11)/20, 9/4, 2) [66 -68].Moreover, we take the hadron matrix elements, such as, for example, 13) GeV 2 , from Refs.[69,70].Finally, the c-coefficients of Eq. ( c(e c α , where the unitary matrices U L/R , E L/R , D L/R , and N diagonalize the SM fermion mass matrices through the following transformations The current experimental constraints and future sensitivities for the various partial lifetimes that we use in our numerical analysis are presented in Table IV.For a recent review on the subject, see Ref. [74].

Numerical procedure
We start our numerical analysis by constructing matrices M u , M e , Y a , Y b , and Y c at the GUT scale, as described in the next few paragraphs.
Since the up-type quark mass matrix M u is approximately symmetric, we have that U R = U * L .This allows us to express M u as We furthermore parametrize the up-type quark mixing matrix U L in terms of the down-type quark mixing matrix D L , the Cabibbo-Kobayashi-Maskawa (CKM) matrix V CKM , and five GUT phases , and η u 3 , as In our analysis, we set η u 1 = η u 2 = η u 3 = 0 since these three phases do not affect the proton decay predictions at all.
We set E L = E R = 1 since M e is diagonal and real.This also means that we can simply construct M e via an equality that reads Y a and Y b are constructed via Eq.(2.39) using the neutrino mixing matrix N = diag(e iη ν 1 , e iη ν 2 , e iη ν3 )V * PMNS as an input.Note that V PMNS contains the CP violating phase δ ν as well as the Majorana phase β ν .We furthermore take Y c to be a general complex 1 × 3 matrix through Once the parameter dependence of M u , M e , Y a , Y b , and Y c is properly accounted for, as described above, we can also construct M d and M ν that are given by Eqs.(2.35) and (2.37), respectively.We treat λ in M ν as a free parameter while the two Higgs VEVs that enter M d and M ν are given by v Λ 1 = v Λ 2 = 174/ √ 2 GeV due to the constraint that tan β of Eq. (3.4) is equal to one.
In summary, the free parameters for our numerical analysis are the unification scale M GUT and the corresponding gauge coupling α GUT , the masses of the fields , the quartic Higgs coupling λ, and the scaling parameter ξ.These 24 parameters are fitted to the experimental observables that are the SM gauge couplings g 1 , g 2 , and g 3 , and the down-type quark masses m d , m s , and m b , while requiring that the current proton decay constraints, as given in Table IV, are satisfied.
Note that the charged lepton masses, the up-type quark masses, the neutrino mass squared differences, the CKM mixing parameters, and the known PMNS mixing parameters are all automatically accounted for.
Since there are more parameters than observables, proton decay cannot be predicted sharply in all decay channels as we will discuss in the next section.But, due to the fact that the neutrino mass matrix is connected to the mismatch between the charged lepton and down-type quark mass matrices, our model predicts the PMNS parameters δ ν and β ν to be in relatively narrow intervals.
The gauge couplings are fitted to their low-energy scale values [61] after the 2loop level running from the high scale to the low scale is performed.To simplify the analysis, we do not run the Yukawa parameters from low scale to the GUT scale using RGEs, and the down-type quark and neutrino masses are directly fitted at the high scale using the high scale values provided in Ref. [82].The χ 2 -function is obtained comparing the theoretical prediction p i with the experimental central value e i , normalized with the corresponding experimental standard deviation σ i of the i-th observable via To minimize the χ 2 -function we apply a differential evolution algorithm.This minimization yields a satisfactory benchmark point and thus proves the viability of our model.Then, starting from this benchmark point, a Markov-chain-Monte-Carlo (MCMC) analysis with a flat prior distribution, involving a Metropolis-Hasting algorithm, is performed giving us a total of 6 × 10 6 datapoints.Finally, we use these points to calculate the highest posterior density (HPD) regions of various quantities.
For the numerical analysis, all parameters are freely varied in such a way that the perturbativity of all Yukawa and Higgs couplings is satisfied.In particular, the absolute values of all entries in Y a , Y b , and Y c as well as the absolute value of λ are all required to be less than or equal to 1.To this end, the scaling parameter ξ ensures that the full parameter space is covered with the chosen parametrization of the matrices Y a and Y b .Furthermore, although we fix some model parameters during the fitting/minimization procedure by directly plugging in experimental central values of some observables, we still vary these parameters in the subsequent MCMC analysis.
It is not necessary to fit the up-type quark masses in our numerical analysis, as already mentioned above.We, however, briefly address perturbativity of the top quark Yukawa coupling.The Yukawa couplings, once the GUT symmetry is broken down to the SM gauge group, need to be matched with the Yukawa couplings of the effective theory, which in our case resembles the type-II 2HDM (two Higgs doublet model).Thus, in the initial basis, the up-type quark sector interacts with one Higgs doublet whereas the down-type quark and charged lepton sectors interact with the other Higgs doublet.Consequently, the SM top quark Yukawa coupling reads y SM t = cos β ŷ, where ŷ is the coupling we are interested in while cos β = 1/ √ 2. To investigate the perturbativity of ŷ we need to study the following RGEs [83]   To generate Fig. 5, we solely consider an effective 2HDM scenario.Our scenario, however, is more complex since several scalar multiplets live significantly below the GUT scale.These fields affect the running of the gauge couplings, as can be seen in Fig. 4, and make the couplings substantially larger, at higher energy scales, when compared to the 2HDM case.This simply means that the negative contributions in Eq. (4.14) are even more important than in the 2HDM case and can thus easily result in smaller coupling ŷ at the high energy scale, if compared to naive expectation.This effect has already been pointed out in Ref. [25], where the RGE running of the charged fermion Yukawa couplings has been implemented.

Results
In this section, we present the outcome of our numerical study.We are interested in the full axion mass range, the predictions for partial proton decay lifetimes, and the viable range of the Dirac CP and Majorana phases of the PMNS matrix.
The axion mass m a is connected to the GUT scale M GUT and gauge coupling α GUT via Eq.(3.10).We can therefore obtain the predicted range of the axion mass by maximizing and minimizing Eq. (3.10).We demand viable gauge coupling unification and correct neutrino mass scale while making sure that none of the current proton decay constraints are violated.We find m a ∈ [0.1, 4.7] neV which we present in Figs. 2 and 3.As discussed in more detail in Sec. 3, this already demonstrates that the full parameter space will be probed by two kinds of future axion DM experiments that are sensitive to either the axion to photon coupling or to the nucleon EDM.
To start our numerical analysis, we find a viable benchmark point from a full χ 2 fit.In particular, for the case of normal neutrino mass ordering, we obtain that Y a = −0.120+ i 0.00943, 0.513 + i 0.200, 0.898 ,  63 GeV, M Σ 3 = 10 13.24 GeV, M Φ 1 = 10 11. 63 GeV, M Φ 3 = 10 5. 28 GeV, M Φ 6 = 10 4. 18 GeV, M Φ 10 = 10 11.63 GeV, α −1 GUT = 15.62, and λ = 1.00.If the proton decay pull is neglected, this choice of the input parameters gives χ 2 below 0.01.This is thus a perfect fit for the gauge couplings as well as for the fermion masses and mixings.It is to be pointed out that even though the benchmark point presented here corresponds to a scenario with T 1 and T 2 masses being two orders of magnitude smaller than the GUT scale, these fields can easily reside at the GUT scale without significantly affecting the value of M GUT .
The PMNS Dirac CP phase, for this benchmark point, is given by δ ν = −48.5 • , whereas the PMNS Majorana phase is β ν = −71.3• .We note that for the case of inverted neutrino mass ordering, no good fit-point can be obtained.This is due to the fact that the Yukawa matrix Y a is needed to generate both the viable neutrino masses and the correct mismatch between the charged lepton and down-type quark masses.In the case of inverted ordering, the first two entries in Y a would need to be somewhat larger than the third entry.This is, however, in conflict with the downtype quark mass fit that requires the first entry of Y a to be smaller than the second and third entries.Therefore, our model predicts that the neutrinos have normal mass ordering.
From the aforementioned benchmark point, we start an MCMC analysis with a flat prior.All obtained points are presented in Fig. 6 in a plane of axion mass vs. par- tial proton decay lifetime in the dominant decay channel p → π 0 e + .We also present the future sensitivities of the DM experiments ABRACADABRA, DMRadio-GUT, and CASPEr Electric, as discussed in Sec. 3, as well as the future sensitivity of the proton decay experiment Hyper-Kamiokande, as discussed in Sec.4.2.Fig. 6 nicely visualizes how various parts of the model parameter space can be probed through the synergy between three different kinds of experiments testing (i) the axion to photon coupling, (ii) the nucleon EDM, and (iii) proton decay.For example, if the axion mass is observed to be above 3 neV, proton decay via p → π 0 e + necessarily has to be seen by Hyper-Kamiokande if our model is realized in nature.Moreover, regardless of whether proton decay will be observed by Hyper-Kamiokande, the former two kinds of experiments will be able to cover the entire parameter space of our model.We are also interested in the proton decay predictions of all two-body decay channels within our model.First, we want to obtain the full allowed range for all partial proton lifetimes, which is for the decay channel p → π 0 e + already hinted in Fig. 6.To this end, we vary all the parameters, including the intermediate-scale particle masses in the MCMC analysis.The 1 σ (dark) and 2 σ (light) HPD results of this analysis are shown in Fig. 7.The blue line segments indicate the current experimental bounds, while the red line segments represent the future sensitivities.(See, for example, Table IV.)Fig. 7 shows that a part of the predicted 1 σ HPD interval for the two decay channels p → π 0 e + and p → η 0 e + will be tested by Hyper-Kamiokande.The large uncertainty in these partial lifetime predictions that are coming from the dependence on the fourth power of the GUT scale can be erased by considering ratios of specific decay channels. 4Fig. 8 shows the prediction of such ratios with the dominant decay channel p → π 0 e + in the denominator.Especially interesting is the prediction for the ratio τ (p → η 0 e + )/τ (p → π 0 e + ), since both τ (p → η 0 e + ) and τ (p → π 0 e + ) will be partly tested by Hyper-Kamiokande.This ratio is predicted very sharply.However, such a sharp prediction for this particular ratio is not only specific to our model but a more common feature of models in which gauge boson mediated proton decay is dominant and in which the contribution involving the c-coefficient c(e c , d) dominates the contribution with the c-coefficient c(e, d c ).Nevertheless, these ratios of partial proton lifetimes provide an interesting additional opportunity to probe our model.
In order to understand the dependence of different decay channels on the flavor structure in the fermion mass matrices, we fix the mass scales to the same values as listed below Eq. (4.18) and only vary the parameters in the fermion mass matrices in the MCMC analysis, computing for each point the proton decay prediction for individual decay channels.We visualize the 1 σ (dark) and 2 σ (light) HPD results of this analysis in Fig. 9, where the blue line segments indicate the current experimental bounds at 90% confidence level presented in Table IV.Interestingly, the partial lifetimes, for some channels, are much more sharply predicted than for the others.The sharp prediction for the decay channels with an antineutrino in the final state is a generic feature for models with a (nearly) symmetric up-type quark mass matrix.On the other hand, the fact that the partial lifetime of the decay channel p → π 0 e + has such a sharp prediction is uncommon and represents a nice feature of our model, which also implies that this decay channel is predicted to be the dominant one. 5he interesting result that some decay channels yield much sharper predictions than others can be understood by investigating the freedom in the mixing matrices that are defined in Eq. (4.8).This is demonstrated in the following example, where we compare the predictions for the p → π 0 e + and p → π 0 µ + lifetimes.The relevant c-coefficients of the two decay channels in question read

Figure 1 :
Figure 1: The 1-loop Feynman diagram which is responsible for neutrino mass generation.

Figure 2 :
Figure2: Expected reach in the m a vs. g aγγ plane for the broadband (Broad) and resonant (Res) strategies of the ABRACADABRA (ABD) experiment[46].The blue line (that lies on the QCD line) corresponds to the prediction of our model.The projected 3 σ sensitivity of DMRadio-GUT[47,48] is also presented in the green shaded region.Furthermore, the expected theoretical reach using the optomechanical cavity method[49] is shown with solid black lines.See text for more details.

Figure 3 :
Figure3: Axion coupling to the nucleon EDM operator as a function of the axion mass.The blue band (that lies on the QCD line) corresponds to the prediction of our model; see text for details.The shaded regions show the sensitivity projections of CASPEr Electric[51,52] in its various phases.Moreover, the ultimate sensitivity limit is given by the nuclear spin noise.

Figure 4 :
Figure 4: Example for the choice of the intermediate-scale particle masses giving gauge coupling unification.

. 1 )
Here, b SM i (b SM ij ) are the SM 1-loop (2-loop) gauge coefficients, while b J i (b J ij ) are the 1-loop (2-loop) gauge coefficients of the multiplets J with intermediate-scale masses M J , i.e., M Z < M J < M GUT .These coefficients are listed in Appendix A. Moreover, β Y i are the Yukawa contributions and H is the Heaviside step function defined as

(c 3 ,c 2 ,c 1 )
2HDM =(−7,−3,7) simplicity, we only consider an effective theory of 2HDM.We evolve these RGEs from the M Z scale to the GUT scale, which we choose to be 10 16 GeV.The initial value of the top quark Yukawa is extracted from the value of the top quark mass in the MS scheme, i.e., m t (m t ) = 163 GeV[84], providing us with ŷ = m t /(v cos β), where v = 174.104GeV.(Recall that m t is the scale-dependent mass and not the physical mass of the top-quark[84].)Our result for the running of ŷ is presented in Fig.5.Clearly, the coupling ŷ remains perturbative up to the GUT scale due to an interplay between the gauge and Yukawa coupling contributions of Eq. (4.14).

Figure 5 :
Figure 5: RGE running of the Yukawa coupling that is relevant for to top quark mass generation.See text for details.

Figure 6 :
Figure6: The generated points from the MCMC analysis presented in the m a − τ (p → π 0 e + ) plane.The current Super-Kamiokande bound is represented by a gray box, while the future Hyper-Kamiokande sensitivity is indicated by a blue dotted line.Moreover, the projected sensitivity of various axion DM experiments is also shown: ABRACADABRA (ABD) with a red dotted line, DMRadio-GUT with a green dotted line, CASPEr Electric with a brown dotted line.For details, see the main text.