Axion and FIMP Dark Matter in a $U(1)$ extension of the Standard Model

In the Standard Model a Dark Matter candidate is missing, but it is relatively simple to enlarge the model including one or more suitable particles. We consider in this paper one such extension, inspired by simplicity and by the goal to solve more than just the Dark Matter issue. Indeed we consider a local $U(1) $ extension of the SM providing an axion particle to solve the strong CP problem and including RH neutrinos with appropriate mass terms. One of the latter is decoupled from the SM leptons and can constitute stable sterile neutrino DM. In this setting, the PQ symmetry arises only as an accidental symmetry but its breaking by higher order operators is sufficiently suppressed to avoid introducing a large $ \theta $ contribution. The axion decay constant and the RH neutrino masses are related to the same v.e.v.s and the PQ scale and both DM densities are determined by the parameters of the axion and scalar sector. The model predicts in general a mixed Dark Matter scenario with both axion and sterile neutrino DM and is characterised by a reduced density and observational signals from each single component.


I. INTRODUCTION
Till date, the Standard Model (SM) has been a very successful theory in describing nature and with the recent discovery of the Higgs boson it is now a complete model that can be extended to higher scale [1]. Nevertheless, the SM cannot be the full theory as it does not address a few important issues, e.g. the absence of a Dark Matter candidate, whose presence is seen in many astrophysical and cosmological observations [2], or the generation of neutrino masses and mixings [3][4][5]. In this paper we add to the Standard Model a new dark sector in order to give a concerted solution to these issues as well as to the strong CP problem, i.e. the very long standing puzzle connected with the suppression of the following term in the SM lagrangian whereG µν = µναβ G αβ and m q is the mass matrix for the SM quarks. This term breaks explicitly the Parity symmetry and generates a non-vanishing electric dipole moment for the neutron of the order d n ∼ 10 −16θ e cm. The present experimental bound d n < 10 −26 e cm, puts a very strong bound on the θ−term i.eθ < 10 −10 [6]. Such small value for the parameterθ is unnatural, asθ is an angular variable, including both the vacuum and the EW contributions, and is expected to be of order O(1) 1 .
One of the most famous and still viable solutions of this problem is the introduction of an axial global symmetry namely the Peccei-Quinn symmetry [10][11][12][13] broken spontaneously at a high scale F a . In that context the θ-term becomes dynamical and it relaxes to the minimum of potential (atθ ∼ 0) in the course of the cosmological evolution. In order to achieve a dynamical θ we need a new pseudoscalar field and therefore we have to extend the SM to contain an additional extra CP-odd degree of freedom, the axion field. But as global symmetries are not respected by gravity [14][15][16][17][18][19], we expect such a PQ symmetry to be only accidental and to be broken by non-renormalisable operators at the Planck scale.
It is well-known that the axion can be a good Dark Matter candidate [20,21], solving another of the open issues of the SM, and that in this case F a ∼ 10 10−11 GeV has to be at an intermediate scale between the EW scale and the Planck scale. Such a mass scale can also fit with the mass of the RH neutrinos and so the light neutrino masses could be sourced by the seesaw mechanism to obtain active neutrino masses in the sub-eV scale. The neutrino masses have been addressed in the context of axion models in [22,23].
In the present work, we propose a new realisation of a gauged U (1) extension of the SM which contains an accidental PQ symmetry and an axion field. We add to the SM gauge group an additional local gauge group U (1) X and a discrete symmetry Z 2 . The gauge symmetry is broken at a high scale by two complex scalar v.e.v.s such that one pseudoscalar field remains physical and plays the role of the axion. Moreover, apart for the new gauge boson and the new scalar sector, we enlarge the particle spectrum by three RH neutrinos and two sets of Dirac fermions, the latter charged under the color gauge group as well as the U (1) X . We investigate the anomaly constraints of the U (1) X and we solve those conditions fixing the new charges such that the operators providing the RH neutrino masses are present in the lower orders in the same v.e.v.s.
The charges of the U (1) X field are also such that a global PQ symmetry will accidentally emerge from this local U (1) X . Moreover the global U (1) symmetry has a color anomaly due to the presence of the additional colored fermions and generates the θ−term for the axion field after the exotic quarks are integrated out. After the QCD phase transition a potential for the axion arises from QCD instantons [24]. Due to gravitational effects not respecting the global symmetry, we expect higher dimensional operators at the Planck scale which violate the global PQ-symmetry but conserve the gauged (U (1) X ) symmetry. Because of this extra contribution from higher dimensional operator, the axion minimum shifts from the valueθ = 0. Such shift is restricted by the bound on the neutron electric dipole moment and so also the range of the scalar v.e.v.s will be bounded from above. As we will see, reaching the value of F a compatible with the full DM density via the misalignment mechanism [20,21] is not always possible in the model, but luckily, if we assume the RH neutrinos mass to be at the GeV scale and the U (1) X gauge coupling very small, we can consider one of the RH neutrinos as a feebly interacting massive particle (FIMP) DM component. We will explore in the following the full parameter space of the model and see that in most regions a mixed DM density arises, lowering possible axion detection signals as the FIMP fills the gap to the observed density.
The paper is organised in the following way. In Section II, we have addressed the model in detail including the discussion of the anomaly constraints, in Section III we show how we can fit the present neutrino data within the model with only two RH neutrinos. In Section IV we have a look on the axion field in the model and compute its relic density by misalignment. In Section V we consider the lightest RH neutrino as a FIMP DM candidate and show that it can fill the gap to reach the observed DM density. Finally we conclude in Section VI.

II. THE MODEL
In the present work, we introduce a U (1) X gauge extension of the SM, together with an extended particle content with three right-handed neutrinos N 1 and N 2 , N 3 with U (1) X charges n e and n, two sets of exotic quarks ψ i L , ψ i R , and χ i L , χ i R with charges under U (1) X as α L , α R , β L and β R , respectively, as well as two singlet scalars φ 1 , φ 2 with U (1) X charges α L − α R , β L − β R . Here we consider different charges in the RH neutrino sector such that one of the RH neutrinos does not couple to the light leptons and remains as a possible DM candidate. Moreover the charges of the scalars are fixed to allow a Yukawa coupling with the two sets of fermions generating the fermion masses at the PQ scale F a . The Standard Model leptons have the same charge as the corresponding RH neutrino, completing a vectorial representation of the U (1) X symmetry and similarly the SM quarks all have the same vectorial charge m.
Among the additional particles, only the two sets of exotic quarks transform non trivially under SU (3) c gauge group, while all the rest are SM singlets. We introduce as well a Z 2 symmetry which ensures the lightest RH neutrino N 1 as stable dark matter candidate and forbids the mixing among the extra fermions with the same QCD charge. In Table (I, II), we show all the SM and beyond SM particles with their corresponding charges under the complete gauge group

A. Gauge Anomaly Cancellation
In order to fix the U (1) X charges of the new particles and the SM particles, we need to consider the gauge anomaly involving the new local gauge group and the SM gauge group, i.e. the following gauge combinations: 1) Y and [Gravity] 2 ×U (1) X . As we have chosen the RH neutrinos and the SM leptons to have vectorial charges under the new symmetry, they do not contribute to most of the anomaly constraints. The following non-trivial conditions on the U (1) X charges arise: where n χ = Nχ N ψ .
same as the previous constraint.
• the [U (1) X ] 3 anomaly gives where in the last line we have used Eq. (2) to simplify the equation to a second order equation in the charges.
Now let us consider more in detail the constraints on the exotic quark charges. We can rewrite Eq. (4) to depend only on the charge difference ∆β = β R − β L and the charge ratios So here we have a quadratic equation in both y and z as long as n χ = 1. The special case N ψ = N χ is not interesting since then we have opposite charges for the new scalar fields and then the PQ symmetry coincides with the U (1) X . The solutions of the above quadratic equation for y are given as, and they simplify to a single non-vanishing rational root for z = ±1, i.e.
The other solution y + = 0 just corresponds to vectorial charges for the exotic fermions and therefore vanishing X charge for the exotic scalar fields.
A physical solution has to be real, so the discriminant must be positive, giving the bound The above equation is quadratic in z and has always real solutions, given by Depending on the sign of 4n 2 χ − 1, different ranges for z are allowed , i.e.
Without loss of generality, we assume that n χ > 1 because the case n χ < 1 is the same just exchanging the role of the two sets of exotic fermions, i.e. α i ↔ β i , N ψ ↔ N χ , n χ ↔ 1/n χ . We can now determine the different values of α L,R and β L,R fixing the relative number of exotic quarks n χ and the ratio z. For the case |z| = 1, the non trivial solutions are unique and rational and we concentrate on this particular case. We show in Table III the U (1) X the solutions of the anomaly conditions for different values of n χ , z = ±1 and one can easily determine the charges for other cases as well as a function of one reference charge α R using eqs. (3), (8) and the definitions of y, z. Note that the value of n χ has to be large in order to obtain very different charges for the scalar fields and so suppress all the possible higher dimensional operators in the scalar potential. Now we are left to fix the charges of the SM fields. In order to be able to write down the Dirac and Majorana mass terms associated with the left and right handed neutrinos as couplings with the exotic scalars, we fix n e and n in terms of the charges of the exotic quarks as n e = 1+nχ 2 (β R − β L ) = 1+nχ 2 y α R and n = 1−nχ 2 (β R − β L ) = 1−nχ 2 y α R . Finally from eq. (5), the quark charge has to be m = − ne+2n 9 = nχ−3 18 y α R .  So far the total number of exotic quarks N ψ + N χ = N ψ (1 + n χ ) has not been fixed and the U (1) X charges only depend on the relative number n χ , but we will see later that in order to suppress higher order operators in the axion potential this number has to be large. Of course the new exotic quarks will affect the running of the SU (3) coupling above their mass scale and contribute to the beta function as the SM quarks, i.e. at one loop we have so that to keep asymptotic freedom only up to 10 exotic quarks are allowed. We will discuss later if this constraint can be fulfilled.
With the charge assignments in Table III, the complete Lagrangian for the particles in the Table I and Table II is given by where D µ corresponds to the co-variant derivative and depends on the different charge assignments of the fields. L SM represents the Lagrangian for the SM particles, the next terms are the kinetic terms for right handed neutrinos, extra fermions and singlet scalars, respectively. L Y uk N contains the Yukawa couplings for leptons and right handed neutrinos including operators to dimension 5 suppressed by the Planck scale. These are as follows We see here clearly that the N 1 RH neutrino is singled out by having a different U (1) X and Z 2 charges and cannot mix with the other two, nor with the light neutrinos. The last terms are higher dimensional operators for the right handed neutrinos which generate a mass once φ 1,2 develop vevs and will be discussed in the next section. The term L Y uk BSM represents the Yukawa sector for the additional exotic quarks and takes the following form Here also we can forbid the cross terms between the two sectors assuming opposite Z 2 parity 2 . These Yukawa terms are important because they are compatible with an accidental global Peccei-Quinn symmetry and by assigning a Peccei-Quinn charge to the exotic fermions and scalars, we can identify one of the pseudoscalars with the axion and generate the dimension five operator coupling of such an axion with the gluons. So our model is a generalisation of the KSVZ type of axion models with multiple exotic quarks [25,26].
In order to break the PQ symmetry as well as the U (1) X symmetry, the scalar fields have to take a vacuum expectation value. The complete potential V(φ h , φ i ) up to mass dimension four has the following form, In principle also higher dimensional operators could appear here, but mostly they depend just on |φ i | 2 due to the different charge assignments. These terms would just change slightly the vacuum expectation values v i , but remain invariant under the phase shift of the scalar fields and therefore do not depend on the axion. Other higher dimensional operators affecting the axion have to be suppressed in order not to spoil the solution of the strong CP problem and they will be discussed in Section IV. After spontaneous symmetry breaking, the scalars can be written as We postpone the discussion of the CP-odd component of the singlet scalars till we reach Section IV. The tad-pole conditions for the potential mentioned above are as follows

B. Neutral Scalars
Once the U (1) X , the PQ and electro-weak symmetry break completely, using the above tadpole conditions, we can write down the mass matrix for the neutral scalars in the basis (H, H 1 , H 2 ) in the following way, From Eq. (18), one can see that to get the mass eigenstates for the neutral scalars (h, h 1 , h 2 ), then one needs to diagonalize the above mass matrix. As we know by the Euler angles (θ 01 , θ 02 , θ 12 ) we can make the 3 × 3 matrix in the diagonal form and the two basis are related in the following way, where c ij = cos θ ij and s ij = sin θ ij (i,j = 0, 1, 2). The v.e.v. of the exotic scalars v 1,2 determine the PQ scale and have therefore to take a high value way above the EW scale. This implies a certain tuning of the mixed coupling λ hφ i . If we consider that quartic couplings λ 12 , λ hφ i , λ φ i ∼ λ h (i = 1, 2) then the mixing angles between the SM Higgs and the beyond SM Higgses become very small, so θ 01 , θ 02 ∼ 10 −9 . For this choice of parameters, the exotic Higgs fields are very heavy M h 1,2 ∼ 10 7 TeV. The other possibility would be to consider h 1,2 in the mass range around TeV scale and mixing with the SM Higgs O(0.1), in that case, we need quartic couplings associated with the extra singlets scalars very small (∼ 10 −7 ). In this work, we consider h 1,2 mass in between these two extreme cases i.e. in the range TeV to a few hundreds TeV.

III. NEUTRINO MASSES
In the present model we can generate the light neutrino masses below the eV scale by the Type I seesaw mechanism. Once the singlet scalars and the SM Higgs doublet take v.e.v.s (as given in Eq. (16)), the Dirac mass matrix in the basis (ν e , ν µ , ν τ ) and (N 2 , N 3 ) takes the following form, where f = e, µ, τ and i = 2, 3. We can absorb the phases of the first column in the lepton fields, but then the phases of the second column cannot be removed. The elements in the first row are coming from higher dimensional operators and are therefore suppressed by v 1 /M P L . The right handed neutrino mass matrix in the basis (N 2 , N 3 ) takes instead the following form, where 3) as obtained from Eq. (13). The elements of the right handed neutrino mass matrix are coming from the higher dimensional operators and give naturally RH neutrino masses at the EW scale. The light neutrino masses arise by the Type I seesaw mechanism as The light neutrino mass matrix can be diagonalised by the PMNS matrix consist of three mixing angles θ 12 , θ 13 , θ 23 and one phase δ. The mixing angles and the mass square differences of the light neutrinos have been measured very precisely by the oscillation experiments [5,27]. The present day 3σ bounds on the oscillation parameters are as follows, • bound on mass squared differences in 3σ range [5,27]  As in our setting only two of the RH neutrinos participate in the see-saw, one of the light neutrinos remains massless, so that the sum of the neutrino masses is of the order of 0.06 eV, consistent with the present cosmological bounds [28,29].
As the oscillation experiments suggest sub-eV scale neutrino mass, we choose our parameters accordingly and the masses of the right handed neutrinos and the Dirac mass matrix elements can be contained in the ranges, In Fig. (1), we show the scatter plots among the different elements of Majorana mass matrix (m R ) and Dirac mass matrix (m d ) after fitting the neutrino oscillation data. As can be seen, there are more than five free parameters associated with the neutrino mass so we can easily satisfy the neutrino oscillation data both for normal and inverted hierarchy. In this figure we show the normal hierarchy case.
As each element of the neutrino mass matrix m ν has to be smaller than 10 −10 GeV, the elements of the Dirac mass matrix have to be small, around 10 −3 GeV, while the remaining suppression comes from the inverse Majorana mass matrix. In the RP, we see a nice correlation among the red points which is due to the bound on the mass square differences.
In the LP and RP of Fig. (2), we show the variation among the different elements of the Dirac mass matrix. Different colours correspond to different elements of the Dirac mass matrix and have been explained in the caption of the figure. From both the figure we can see that individual elements can reach values up to 10 −3 GeV but only for the last two generations of light neutrinos. The elements related to the electron neutrino are smaller.
Since the operators giving rise to the RH neutrino mass matrix are of dimension 5, even if the scalar v.e.v.s are large, the mass matrix entries remain in the 10-100 GeV range and so without ad-hoc cancellations the mass eigenvalues are in the same range.

IV. AXION
In the model section we have given a detailed discussion of the scalar particle content, here we want to focus on the CP odd part of the two singlet scalars φ 1 and φ 2 and show that one combination of those plays the role of the QCD axion.
If we write the scalar fields as where v 1 , v 2 are vacuum expectation values of the two fields, we can identify the fields a 1 , a 2 with the massless Goldstone bosons. The kinetic term for the above two scalars with the extra gauge boson Z BL takes the following form, where F µν X is the field strength tensor for the U (1) X gauge boson Z X and the covariant derivative is written as, where g X and q i are the U (1) X gauge coupling and the charges of the corresponding fields, respectively q 1 = α L − α R , q 2 = β L − β R as in Table III. The value of these charges can be rescaled by fixing the free parameter α R and if not otherwise specified, we choose it such to have q 2 = 1. Using the expression Eq. (25), Eq. (26) takes the following form, where Here we see that one combination of the pseudoscalars a i mixes with the Z X gauge boson and can be absorbed into that field as the longitudinal component of the massive gauge boson, so define the new pseudoscalars as, where in the last expressions we have used Eq. (2) and q 1 q 2 = α L −α R β L −β R = − Nχ N ψ = −n χ . So the B scalar is absorbed by the massive gauge boson Z X , while the extra Goldstone boson A remains massless at the scale of U (1) X breaking.
One can easily notice that the complete Lagrangian (see Eq. (12)) has an accidental (axial for the exotic sector) U (1) global symmetry with the charge assignments shown in Table II. We call this accidental symmetry the global Peccei-Quinn (PQ) symmetry. With a field dependent axial transformation for the exotic quarks, we can eliminate the pseudoscalar fields from the Yukawa couplings and generate the axion field coupling with the gluons as follows, from Eq. (30), one can see that the decay constant corresponding to the CP odd particle A is for N ψ = 1, and then θ a = A Fa varies in the range 0 to 2π. The PQ symmetry is broken not only by the spontaneous breaking of the U (1) X , but also by instanton effects, such that after the QCD phase transition, the axion has the usual periodic potential where M a is the axion mass and given by where m u , m d are the quark masses, while m π = 139.57 MeV is the pion mass and f π ∼ 130 MeV the pion decay constant. So a v.e.v. of the axion A can cancel theθ term and we can define a field with vanishing v.e.v. asθ a = θ a +θ.

A. Gravitational effects in the axion potential
It is well-known that Planck scale physics does not conserve global symmetries. Wormholes [14,16,18] are one example in support of this and black holes also can absorb charges without effect, as predicted by the no-hair theorem [30]. More explicitly, if in a scattering process a virtual or non-virtual black hole is formed from the initial state with definite charge, such global charge is not conserved in the final state due to Hawking evaporation [31].
In this work, we assume a gauged U (1) X symmetry as underlying the accidental PQ symmetry U (1) P Q . Therefore, at the Planck scale we can expect only higher dimensional operators which conserve the gauge symmetry but violate the global U (1) P Q symmetry. Thanks to our choice of U (1) X charges, many operators are forbidden by the local gauge symmetry and so the lowest dimensional operator of this kind is of high dimension [32][33][34], where g = |g|e iδ is a complex coupling 3 . M P L = 1.22 × 10 19 GeV is the Planck mass. Note that this operator is invariant with respect of the U (1) X local symmetry due to the anomaly condition in eq. (3), but for N ψ = N χ it breaks the PQ symmetry. For large number of fields N ψ and/or N χ , it is strongly suppressed and does not affect the axion potential too much. Indeed after spontaneous symmetry breaking, it gives an extra contribution to the axion potential as follows [32], where M g a is the quantum-gravitational induced axion mass and have the following expression, and p = −N ψ + N χ = N ψ (n χ − 1). So the axion potential becomes where r g = (M g a ) 2 M 2 a 1. Therefore the minimum of the axion potential (as shown in Eq. (32)) is shifted fromθ a = 0 due to the presence of the phase δ in the second term. The axion can still address the strong CP-problem if the shift lies below the upper bound, ∆θ ≤ 10 −10 . By minimising the complete potential of the axion field in Eq. (37), we obtain the following value of θ a (defined it as ∆θ), ∆θ = r g |p sin δ| 1 + p 4 r 2 g + 2p 2 r g cos δ 3 By assuming this coupling complex, we are not explicitly assuming that Planck scale contribution is complex.
The phase can also come from the chiral rotation of the fermions to diagonalise the mass matrix (m q ) and in that case this angle is δ arg(det[m q ]) If we assume |p sin δ| ∼ 1, p 2 r g 1 since for very high v.e.v. of the singlet scalars (v i ∼ 10 10 GeV, i=1, 2) the parameter r g 1, Eq. (38) takes the following form So taking the values N ψ = 1 and N χ = n χ 1, we have a stronger dependence on v 2 . In the Section IV B, we show the variation of ∆θ in the v 1 − v 2 plane.
For the special case v 1 = v 2 and |g| ∼ 1 the expression gives Using Stirling's formula for n χ 1, this can be reduced to We see that we need n χ to be sufficiently large to give a strong suppression of the shift in the vacuum value of the θ angle, as the dominant behaviour is determined by the exponential of n χ (1 + ln Fa M P L ) < 0. Considering n χ = 9, the largest number of fields compatible with asymptotic freedom in QCD up to the high scale, we obtain ∆θ ∼ 29.41 × 10 −10 F a 10 10 GeV 10 . ( so we see that we need unfortunately to introduce of the order of 10 exotic fermions or more to satisfy the present experimental constraints without relying on fine-tuning or reducing the axion decay constant around or below the astrophysical bounds from the cooling of the 1987A Supernovae [35], giving F a ≥ 2 × 10 8 GeV. We have checked numerically that n χ = 8 is ruled out from the θ−parameter bound. In Figure 3 we show the allowed parameter space in the v 1 versus v 2 plane for n χ = 9, 10, 11, 12 obtained numerically with the full expression in eq. (38). As we can easily see from the figure, increasing the n χ value, enlarges the allowed range of v 2 . This is because larger n χ gives a larger Planck scale suppression. We can also see a turnover of the blue region near v 2 ∼ 2 × 10 10 which is due to the lower limit for the axion relic density we have considered in the present work i.e. Ω a h 2 ∼ 10 −4 . Moreover, the right handed neutrinos mass also depends on the v.e.v.s v 1,2 , so the cyan lines represent the different values of right handed neutrino mass for y ij = 10 −3 (i, j = 2, 3). In the next sections we will investigate more in detail these cases not too far from the asymptotic freedom case with respect to the DM density. Cyanide lines represent the different mass eigenvalues of the right handed neutrino mass matrix (as defined in Eq. (22)) when y ij = 10 −3 (i, j = 2, 3). The plots from top left to right bottom correspond to N ψ = 1 and n χ = 9, 10, 11 and 12 , respectively.

B. Axion as Dark Matter
We consider here a sub-eV scale axion and it can be a very good cold dark matter candidate. The axion contribution to the cold dark matter density from the misalignment mechanism takes the following form [36], where the initial value θ i ∼ 1, as averaged over many domains if the PQ symmetry is broken after inflation or if the initial condition is not fine-tuned for the case of PQ symmetry breaking before inflation. In Fig. (3), we see that the bound from ∆θ constrains the value F a below 10 12 GeV, so that for most of the parameter space, the axion density does not constitute the full dark matter relic density of the universe, apart for an unnaturally large θ i . Then we need an additional Dark Matter component or an additional axion production mechanism. Topological defects like axion strings could contribute to DM production, in the case when the PQ and U (1) X are broken after the inflationary epoch, but the predictions are still uncertain [37][38][39][40]. Generically that additional production could move the region of axion DM to a lower value of F a by about one order of magnitude [40]. But in our case the strings formed are also U (1) X strings, which allows them to also decay in other particles than axions and could reduce axion production. We therefore consider the conservative case of main production via misalignment. Note that our model has no domain wall problem for N ψ = 1, even if N χ is large, as the smallest of the two plays the role of the overall domain wall number in eq. (31) and we do not expect axion production from domain walls.
In Fig. (4), we show the scatter plots in the F a − v 1 plane after satisfying the bound on the θ parameter, ∆θ < 10 −10 and showing in color the axion DM density via misalignment from Eq. (43). The four plots from top-left to bottom-right correspond to N ψ = 1 and n χ = 9, 10, 11 and 12, respectively. For low values of n χ , the constraint from the neutron dipole moment leaves open only a parameter region with a too small contribution to the DM relic density. For n χ = 10, 11 the maximal value of F a reaches 4 × 10 10 , 2 × 10 11 GeV respectively.
For n χ = 12 or larger, we can reach a maximal value F a = 10 12 GeV, so the axion can provide the full DM density in the Universe. But still in a large region of the parameter space the density is too small and therefore we need another DM component. In this model such component can be the third right-handed neutrino N 1 , which is stable thanks to the Z 2 symmetry, and will be discussed in section V C. Axion coupling with SU (2) L and U (1) Y gauge bosons In the present work, in order to single out one of the RH neutrinos, the first generation of leptons are charged under U (1) P Q symmetry. Because of these non-trivial charges and the fact that the EW couplings are non-vectorial, we obtain non-vanishing couplings of the axion field with the SU (2) L gauge bosons and the U (1) Y gauge boson. As defined in Table (I, II), the PQ charge associated with L e and e R are equal to −2q a . The W i W i A (i = 1, 2, 3) coupling in the present model due to the nontrivial contribution from the SU (2) L × SU (2) L × U (1) P Q anomaly term can be expressed as In the same way, we also obtain a non-trivial contribution from the U (1) Y × U (1) Y × U (1) P Q anomaly term which reads as, where we have used Y Le = − 1 2 and Y e = −1 . Combining the W 3 and B Y we obtain a vanishing coupling to the photon due to an exact cancellation. Indeed in this model the exotic fermions are not charged under the EM interaction, so they also do not contribute, and the leptons have vectorial couplings with respect of the electromagnetic U (1). Nevertheless a non-vanishing coupling arises from the pion-axion mixing giving [24] Therefore, the axion has substantial interaction with the photons and with the electroweak gauge bosons which can be measured at ongoing and future experiments [42]. In our model the density of the axion can be lower than the total DM density, but the axion-photon coupling can be larger due to the constraint on F a from ∆θ, so that it may be still detectable. Unfortunately the mass range is shifted accordingly, leading to a larger mass than those recently explored by ADMX [41], but those mass range could be accessible in the future [42].

V. FIMP DARK MATTER
From Fig. (4) we see that the axion cannot be the sole dark matter component for n χ < 12 and neither for the entire region in the v 1 −F a plane if n χ = 12, since the EDMs bound does not allow large enough F a to match the total DM density. Fortunately, we have another DM candidate in the model since we have assigned Z 2 parities to the right handed neutrinos such that N 1 is the lightest Z 2 odd particle and stable, since the Higgs fields do not break this global symmetry. In order for N 1 to be a good DM candidate, it has to be produced in the early Universe in right amount; due to its very weak interactions a natural way for that to happen is through the FIMP mechanism [43]. Indeed the only interactions of the right handed neutrino N 1 are the U (1) X gauge interaction and the non-renormalisable operator generating its Majorana mass as where n e and g X are the charge of N 1 and gauge coupling associated with U (1) X gauge group.
We can produce the N 1 from the decay of scalars contained in φ 1,2 (h 1 , h 2 assuming the EW symmetry is unbroken and v = 0 in the mass matrix in eq. (18)) and of the gauge boson (Z X ). The production of N 1 via the freeze-in mechanism requires feeble couplings and therefore also the gauge coupling g X has to be very small and the mass of the new Z X boson from v 1 , v 2 will be naturally in the TeV range. Consequently, the extra gauge boson will not reach thermal equilibrium in the early Universe and will be itself produced via the FIMP mechanism, similarly to what happens also in the model in [44]. Once Z X is produced sufficiently, it will then decay to two N 1 particles. This kind of dark matter production is similar to the superWIMP mechanism, just with a non thermal density for the decaying particle, we can call it 'superFIMP' production. We can determine the DM relic density once we know the amount of Z X produced and solved the relevant Boltzmann equations.
For a decay process A → BC, where C is the DM candidate feebly coupling and therefore out-of-equilibrium, the relic density obtained from the mother particle decay in equilibrium is simply given by [43] where g A counts the internal d.o.f. of the mother particle, while g s , g ρ (∼ 100) denote the number of relativistic d.o.f. in the entropy and energy density of the Universe. The DM production from the Higgs decays for h i , i = 1, 2 is of this type assuming that the scalars are in equilibrium thanks to the interaction with the SM Higgs and the additional quartic and cubic couplings. Indeed due to the scatterings with SM fields, the exotic Higgs scalars are in equilibrium as long as their quartic couplings are larger than 10 −13 . The RH neutrino energy density produced from these decays is therefore well-described by the formula above taking into account that two neutrinos are produced in each decay as Here the decay of the Higgs into neutrinos is determined by the same operator giving the neutrino mass term and so we have for each Higgs field, assuming small mixings so that H 1,2 are nearly equal to h 1,2 , leading to suppressed by both the Higgs mass and the scalar v.e.v. v 1 = v 2 = F a n 2 χ + 1 . We see therefore that the FIMP production from direct Higgs decay into neutrino is often too suppressed to give a substantial DM component as long as For the production of neutrinos from the gauge boson Z X , we have instead to consider the coupled Boltzmann equations, as also Z X does not equilibrate. The production of Z X itself takes place from the decay of the Higgs fields since it is lighter than the scalar fields due to the suppressed gauge coupling. Otherwise it can also be produced in the scatterings of the fermions of the SM, which are charged under U (1) X , but the scattering processes are negligible compared to the decays as discussed in [44].
The Boltzmann equation to determine the phase-space distribution f Z X (p) of the gauge boson Z X has to be solved to be able to compute its non-thermal average decay rate into N 1 as described in [44][45][46]. The non-thermal average of the decay width of the Z X is indeed given by, In the equation for the yield, both a production and a decay term appear as they may be active at the same time. Using as a time variable z = M h 1 /T , with M h 1 the mass of the lighter Higgs scalar, we have the coupled equations where M P L = 1.22 × 10 19 GeV is the Planck mass and g (z) = gs(z) √ gρ(z) Here we are neglecting the backreaction due to the N 1 density, as it is negligible in the whole regime. Finally, one can solve the Eq. (54) and determine the co-moving number density of N 1 and its relic density by the following expression, In Fig. 5, we show the variation of the co-moving number density of the gauge boson Z X and the relic density of the DM N 1 . The LP shows the curves for three values of g ef f X = q 2 g X , which are 1 × 10 −12 , 3 × 10 −12 and 6 × 10 −12 , respectively. We can notice that for all the three values of g ef f X the co-moving number density of the extra gauge boson is the same. This is because the decay width for the process h i → Z X Z X decay is independent of the gauge coupling as in As the extra gauge boson Z X is produced from the decay of the Higgses h i (i = 1, 2) in thermal equilibrium, we can determine the relic density of Z X on the plateau from Eq. (48). The decay rate is giving corresponding for g ef f X = 3 × 10 −12 and the parameters as in Fig.5 to Ω F IM P Z X h 2 ∼ 610. On the other hand if we compute the Z X relic density using the full Boltzmann equation including the back-reaction as in Eq. (53), then at the plateau the co-moving number density is Y Z X ∼ 7.14632 × 10 −10 , corresponding to Ω BE Z X h 2 = 590 using eq. (54). So we can conclude that the back-reaction in the Z X equation does not have a strong impact and that actually most of the Z X decays happen after the Z X density has frozen-in, at least for small g X . Indeed, contrary to the production, the decay processes depends on the exact value of g X , as Γ Z X ∝ g 2 X . So the larger g X , the faster Z X decays and consequently the faster a N 1 density is produced. As long as the Z X production and decay are well separated in time, we can approximate the final density of N 1 by using a formula analogous to the SuperWIMP mechanism, but for the FIMP Z X density as mother particle density, i.e.
(Ω SF where BR Z X →N 1 N 1 is the branching ratio of decay into two N 1 . So we see that also in this case the N 1 density is independent from the exact value of the gauge coupling g X as long as it is small enough to avoid equilibration and separate the production and decay epochs. Then the approximate SuperFIMP formula above matches the relic density obtained from the full Boltzmann equation. We can also have an analytic estimate of the dependence on the Higgs parameters for the N 1 production as determined by the v.e.v.s of the scalar fields, their masses and v.e.v.s. The last expression is valid for the case v 1 = v 2 = F a n 2 χ + 1, otherwise the largest v.e.v. dominates and the suppression scale becomes n 2 χ v 2 i instead of F 2 a (n 2 χ + 1) 2 . So we see that again the neutrino energy density is suppressed by the PQ scale or the v.e.v. and n χ , but enhanced by a large Higgs and neutrino masses 4 .
The branching fraction of Z X decay into N 1 's is just determined by the charges of the SM fields, n, n e and m, assuming all other exotic states to be too heavy, i.e. we have for large n χ 1.
In the RP of Figure 5, we show the variation of DM relic density from the decay of the gauge boson and the Higgs fields for three different values of DM masses, 1 GeV, 6 GeV and 10 GeV Since the v.e.v. v 1 = 2×10 11 GeV, DM is produced dominantly from the Z X decay and the Higgs contribution is negligible as visible in the figure. For N 1 masses above 100 GeV, the DM becomes overabundant because the coupling is directly proportional to the DM mass. In considering the production of N 1 from the Higgs's decays we have assumed here small mixing among the Higgses without loss of any generality as for larger mixing our conclusion remains same.

A. Scatter Plots
In generating the scatter plots among the different parameters of the present model, we have varied the model parameters in the following range, In the scatterplots we have both contributions from the decay of the scalar particles in equilibrium (h 1 , h 2 ) and from the Z X decay not in equilibrium. For the latter, in order to speed up the computation, we use the approximate expression where ξ is an efficiency factor including the branching ratio, ξ < 1. We estimate the ξ factor from the Fig. (5) obtaining ξ = 1 10 , in good agreement with the branching fraction estimate above. We compare the N 1 energy density with the cosmological 3σ band measured by the Planck collaboration [28,29].
In Fig. (6), we show scatter plots in v 1 − M N 1 (LP) and M Z X − g X (RP) planes after satisfying the bound on the ∆θ parameter as well as the constraint on the DM relic density. Cyan points represent the contribution from the scalar decays, green points correspond to the contribution from the Z X decay and the red points correspond to the total contribution coming from both decays. The contribution from the scalars decay to DM is proportional to, Ω F IM P N 1 , so we can see a linear relation among the cyan points till v 1 < 5 × 10 11 which is consistent with the analytical formula. For v 1 > 5 × 10 11 , the N 1 relic density is controlled by the other v.e.v. v 2 , which is restricted by the bound from the θ parameter. So from that point the dependence on v 1 is weaker and the correlation between the mass and the v.e.v. is broken. We see that The contribution from Z X decay is proportional to the ratio of the dark matter mass and the Z X mass as in the SuperWIMP case, so that there is a sharp correlation among the green points. The red points correspond to the total contribution which is sum of the scalar and gauge boson contribution. Most of the green points are not allowed because the scalars overproduce DM, while most of the cyan region is also ruled out because in those region DM is overproduced by the decay Z X . Therefore, taking into account both contributions, only the region at large v 1 remains. There the value of F a is actually determined by the other (smaller) v.e.v. v 2 so that the axion energy density can be much smaller than the N 1 energy density. Moreover, the DM relic density constraint can be satisfied also for higher values of the DM mass, if the phase space suppression starts to play a role in the decay rate of the Z X boson into neutrinos reducing the BR.
On the other hand,in the RP we show the points satisfying the DM energy density bound in the M Z X − g ef f X plane. The points there follow the very simple correlation M Z X g X v 1 due to the symmetry breaking. Here we see that for the production via Higgs decay the viable region is larger as we vary the h 1,2 in a broad window, but for the most part of the parameter space the red and green point overlap, as the dominant production comes from the Z X decay and it is independent of the exact gauge coupling and gauge boson mass, as given in eq. (59). Nevertheless, the gauge coupling has still to be small enough to avoid equilibration of the gauge boson and the neutrino, g X ≤ 10 −6 , as well as ensure M h 1,2 > M Z X .

VI. CONCLUSION
We have presented a new U (1) X extension of the Standard Model providing solutions and relating three different open questions in the SM. Indeed the model contains an accidental PQ symmetry and an axion field solving the strong CP problem. The same scalar v.e.v.s breaking the PQ symmetry give masses to the RH neutrinos and contributes to the light neutrino masses via the seesaw. The RH neutrino masses are due to non-renormalisable operators, so that they are at the electroweak scale. Only two of the RH neutrinos take part to the seesaw mechanism, so that the lightest active neutrino remains massless. The charges of the exotic fields in the model are fixed by anomaly freedom of the new gauged U (1) X and these constraints also strongly reduce the number of possible non-renormalisable operators, allowing for an accidental PQ symmetry to emerge.
In the model there are naturally two DM candidates, the axion and the only RH neutrino that does not couple to the active neutrinos and is stable due to the Z 2 symmetry. Depending on the field content of the model, the constraints on the scale F a from higher order operators violating the PQ symmetry are such that the axion cannot provide the full DM density and the gap has to be filled by the RH neutrino. Different scenarios arise depending on the number of exotic quarks present in the model encoded in the parameter n χ for N ψ = 1. Such parameter fixes the order of the gravity-induced operators that modify the θ angle breaking the accidental PQ symmetry. Only the value n χ ≤ 9 allows for the asymptotic freedom of QCD to survive to the high scale and in that case the EDM bounds limit the value of F a so that the axion cannot provide a substantial number density for Dark Matter. Then the RH neutrino can give the total DM density for small gauge couplings g X and masses in the tens of GeV with the main production from the Z X gauge boson decay out of equilibrium. For values of n χ ≥ 12 instead also mainly axion DM is possible, with a small RH neutrino component suppressed by a small mass or large v.e.v.. Indeed the parameters determining the two energy densities depend on different combinations of the two vacuum expectation values and while for the axion field the key parameter is F a , determined by the smaller v.e.v., the density of the RH neutrino is sensitive as well to the Higgs masses and the larger v.e.v..
Since the DM density of the axion and the RH neutrino is reduced compared to the total density, we expect a weaker signal in direct detection experiments. Indeed the coupling of the axion to photons is the same as for hadronic KSVZ models with electrically neutral exotic quarks, which is now reached by the ADMX experiment for the mass window 2.6 − 4.2µ eV. In our case though often the axion mass is larger due to the smaller allowed F a , in the range 5 µeV to 28 meV and the energy density smaller than Ω DM . The N 1 instead couples to the SM fields via the U (1) X gauge interaction with a tiny gauge coupling. So in this scenario the new Z X gauge bosons is in the TeV range even if the v.e.v.s of the Higgs fields is of order F a . Unfortunately the gauge coupling is too small to allow for gauge boson production at a collider. On the other hand the direct detection of N 1 Dark Matter via Z X and Higgs exchange is independent of the small gauge coupling but suppressed by the scalar v.e.v.s and masses. The RH neutrino sector at the EW scale can provide additional signatures of the model, but in that case based on the standard seesaw coupling as discussed in [47]. The model may be able also to include leptogenesis at the low scale via the ARS mechanism or via resonant leptogenesis, if these two states N 2,3 are sufficiently degenerate. Finally, the exotic scalar fields mix with the SM Higgs after EW symmetry breaking and they open additional Higgs decays into axions and modify generically the Higgs couplings and possibly as well the Higgs potential [48], but these effects are always suppressed by the PQ scale F a and difficult to measure.

VII. ACKNOWLEDGEMENTS
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska -Curie grant agreement No 860881-HIDDeN. This work used the Scientific Compute Cluster at GWDG, the joint data center of Max Planck Society for the Advancement of Science (MPG) and University of Göttingen.