Axion dark matter in a $3-3-1$ model

Slightly extending a right-handed neutrino version of the $3-3-1$ model, we show that it is not only possible to solve the strong CP problem but also to give the total dark matter abundance reported by the Planck collaboration. Specifically, we consider the possibility of introducing a $3-3-1$ scalar singlet to implement a gravity stable Peccei-Quinn mechanism in this model. Remarkably, for allowed regions of the parameter space, the arising axions with masses $m_a\approx$ meV can both make up the total dark matter relic density through non-thermal production mechanisms and be very close to the region to be explored by the IAXO helioscope.


I. INTRODUCTION
The impressive observation that almost thirty percent of the energy content of the Universe is due to dark matter (DM) is challenging our understanding of particle physics and cosmology. For a historical review see Ref. [1]. Much effort have been done in order to unravel the nature of DM. Experiments designed to detect weakly interacting massive particles (WIMPs), the, so far, DM candidate paradigm, have failed in providing positive results [2,3].
At the same time, the Large Hadron Collider (LHC) has not been able to produce any signal of a DM candidate, as is the case of the lightest supersymmetric partners of the standard model (SM) neutral particles (gauge or scalar), called neutralinos, or gravitinos (partners of the graviton) [4].
As a consequence of these negative results, it is noticeable the growing interest in studying axions and axionlike particles (ALPs) because they are well motivated alternatives to WIMPs. Moreover, they can be linked to solutions of still intriguing astrophysical phenomena [5]: (i) ALPs may be the explanation to the TeV photon cosmic transparency if there are gamma ray ←→ ALP oscillations. If so, gamma rays could be converted to ALPs due to the magnetic fields near active galactic nuclei, for instance, traveling "freely" for a long distance to our galaxy and then reconverted into gamma rays in the galactic magnetic fields; (ii) also, ALPs may explain the anomalous energy loss of white dwarfs because from the luminosity of this kind of stars it is inferred that a new energy loss mechanism is needed. In the present scenario, this mechanism could be related to axions or ALPs bremsstrahlung if they directly couple to electrons. All these astrophysical processes constrain the relevant parameters describing axions and ALPs physics. In fact, besides these theoretical arguments for considering axions and/or ALPs, there is also much experimental effort searching for this kind of particles [6]. A variety of experiments have been designed and, in general, they are classified as haloscopes, helioscopes and light-shining through a wall, and most of them are based on the conversion of axions or ALPs into gamma rays in the presence of strong magnetic fields [7].
The axion field was initially introduced as a dynamical solution for the so-called strong CP problem. This problem comes from the extra term which has to be added to the QCD Lagrangian due to the nontrivial structure of the QCD vacuum: where G a µν is the gluon field strength andG a µν its dual. This θ-term violates P, T and CP symmetries and, hence, it induces a neutron electric dipole moment (NEDM). In order to be in agreement with experimental NEDM data the value of the θ parameter must be θ 0.7 × 10 −11 [8]. The strong CP problem is, then, to explain why this parameter is so small. After including weak interactions, the coefficient of the GG term changes toθ = θ − arg det M q , where M q is the quark mass matrix. The Peccei-Quinn (PQ) solution to this problem is implemented by introducing a global U(1) symmetry that must be spontaneously broken and afflicted by a color anomaly. The axion is then the Nambu-Goldstone boson associated to the breaking of that U(1) symmetry, which is now known as the U(1) PQ symmetry. After including the axion field, a(x), the total Lagrangian has a term proportional to the color anomaly N C : wheref a /N C ≡ f a is the axion-decay constant and it is related to the magnitude of the vacuum expectation value (VEV) that breaks the U(1) PQ symmetry. We also have that the divergence of the PQ current, ∂ µ J µ PQ , is N C g 2 s 32π 2 G a µνGa µν = 0. Hence, the CP violating term GG is now proportional to (θ + N C a(x)/f a ) and it is shown that a(x) = −f aθ /N C minimizes the axion effective potential so that, when the axion field is redefined, a(x) → a(x) − a(x) , the CP violating term GG is no longer present in the Lagrangian, solving in this way the strong CP problem. Although the axion is massless at tree level, it is, in fact, a pseudo-Nambu-Goldstone boson since it gains a mass due to nonperturbative QCD effects related to the U(1) PQ color anomaly. The axion mass and all its couplings are governed by the value of f a . The original conception of the axion was ruled out long ago because f a was thought to be near the electroweak scale, implying in a "visible" axion, in contradiction with laboratory and astrophysical constraints. Few years after the PQ proposal it was realized that for large enough values of f a the axion could be a cold dark matter candidate [9][10][11].
In fact, for high symmetry breaking scales, the axion is a nonbaryonic extremely weaklyinteracting massive particle, stable on cosmological time scales, which makes it a candidate to dark matter. Later in the text we discuss the constraints on f a coming from NEDM, "invisibility" of the axion, and astrophysical data.
In order to consider the axion a viable DM candidate we must deal with its relic abundance which strongly depends on the history of the Universe. In particular, the cosmological scenario for the axion production changes significantly if the PQ symmetry is broken before or after the inflationary expansion of the Universe. The main issue related to the order of these events concerns the axion-production mechanisms. There are production mechanisms due to topological defects, like axionic strings and domain walls, that are comparable to the vacuum misalignment one. Hence, on one hand, if the PQ-symmetry breaking occurs before inflation, inflation will erase these topological defects. On the other hand, if the PQ-symmetry breaking happens after inflation, it is expected an additional number of axions to be produced due to the decay of the topological defects, affecting directly the relic abundance estimative. In this work we consider axions as DM candidates in the so-called post-inflationary scenario, when the reheating temperature, T R , is high enough to restore the PQ symmetry, T R > T C ∼ f a , which will be broken at a later time, when the temperature of the Universe falls below the critical temperature T C .
As we can see, axions present some features with relevant implications not only in particle physics but also in cosmology and it is also a strong indication that physics beyond the SM is in order. In this vein a large variety of models, extensions of the SM, has been proposed.
Most of them claim for very appealing achievements relating the DM solution to another yet unsolved issue in particle physics [12][13][14], as it is the case of the lightness of the active neutrino masses, the smallness of the strong CP violation, or the hierarchy problem, for instance.
Among others, a way of introducing new physics is to consider a model with a larger symmetry group. In particular, there is a class of models based on the SU(3) C ⊗ SU(3) L ⊗ U(1) X gauge group (the so called 3 − 3 − 1 models, for shortness), which are interesting extensions of the SM. In general, these 3 − 3 − 1 models bring welcome features which we review very shortly here. We can take advantage of the larger group representation to choose the matter content in order to introduce new degrees of freedom which are appropriate to implement, for instance, a mechanism to generate tiny active neutrino masses, in the lepton sector. The quark sector will also have new degrees of freedom and, depending on the particular representation, the model can have quarks with exotic electric charges or not.
The issue of the chiral anomaly cancellation is solved provided we have the same number of triplets and anti-triplets, including color counting. Then, considering that we have the same number of lepton and quark families, say n f , we find that n f must be three or a multiple of three. However, from the QCD asymptotic freedom we find that the number of families must be just three in order to get the correct, negative, sign of the renormalization group β function. Note that, contrarily to the SM, the total number of families must be considered altogether in order to get the model anomaly free. Hence, the number of families and the number of colors are related to each other by the anomaly cancellation condition. This fact is a direct consequence of the 3 − 3 − 1 gauge invariance and it can be seen as a hint to the solution to the family replication issue. We can still mention other interesting features: (i) the electric charge quantization does not depend if neutrinos are Majorana or Dirac fermions [15]; (ii) the model described in Refs. [16][17][18] presents the relation t 2 = (g /g) 2 = sin 2 θ W /(1 − 4 sin 2 θ W ), which relates the U(1) X and the SU(3) L coupling constants, g and g, respectively, to the electroweak θ W angle. This relations shows a Landau-like pole at some O(TeV), energy scale, µ, for which sin 2 θ W (µ) = 1/4 [19], and it would be an explanation to the observed value sin 2 θ W (M Z ) < 1/4. (iii) The Peccei-Quinn symmetry, usually introduced to solve the strong CP problem, can be introduced in a natural way [20]. In this work we consider a version of a 3 − 3 − 1 model where a gravity stable PQ mechanism can be implemented. We analyze the conditions under which the axion, resulting from the spontaneous breaking of the PQ symmetry in this model, can be considered a dark matter candidate.
This work is organized as follows. In Sec. II, we present the general features of the 3 − 3 − 1 model, including its matter content, Yukawa interactions and scalar potential. In Sec. III we show the main steps to make the axion invisible and the PQ mechanism stable against gravitational effects. We also show the axion effective potential from which its mass is derived. In Sec. IV we consider the axion production mechanisms in order to compute its abundance in the Universe. Results for the vacuum misalignment and decay of the string and string-wall system mechanisms are given. In Sec. V we confront the predictions from the previous section with the observational constraints, coming mainly from the Planckcollaboration results for the DM abundance, the NEDM data and direct axion searches, in order to constrain the parameter space of the model. Section VI is devoted to our final discussions and conclusions.
In order to generate the fermion and boson masses, the SU(3) C ⊗ SU(3) L ⊗ U(1) X symmetry must be spontaneously broken to the electromagnetic group, i.e., to the U (1) Q symmetry, where Q is the electric charge. To do this, it is necessary to introduce, at least, three SU(3) L triplets, η, ρ, χ, as shown in Ref. [30], which are given by Once these fermionic and bosonic fields are introduced in the model, we can write the most general Yukawa Lagrangian, invariant under the local gauge group, as follows with where ijk is the Levi-Civita symbol and a , i, j, k = 1, 2, 3 and a, b, s, t are in the same range as in Eq. (3). It is also straightforward to write down the most general scalar potential consistent with gauge invariance and renormalizability as with We have divided the total scalar potential V (η, ρ, χ) in two pieces, (4,5)R , and all the other fields even by the symmetry), and V Z 2 (η, ρ, χ), which breaks Z 2 . This discrete symmetry is motivated by the implementation of the PQ mechanism as shown below.
It is well known that the minimal vacuum structure needed to give masses to all the particles in the model is which correctly reduces the SU (3) C ⊗ SU (3) L ⊗ U (1) X symmetry to the U (1) Q one. In principle, the remaining neutral scalars, η 0 3 and χ 0 1 , can also gain VEVs. However, in this case, dangerous Nambu-Goldstone bosons can arise in the physical spectrum, as shown in Ref. [37]. In this paper, we are going to consider only the minimal vacuum structure given in Eq. (14).

III. IMPLEMENTING A GRAVITY STABLE PQ MECHANISM
The key ingredient to implement the PQ mechanism is the invariance of the entire Lagrangian under a global U (1) symmetry, called U (1) PQ , which must be both afflicted by a color anomaly and spontaneously broken [38][39][40][41]. In general, the implementation of the PQ mechanism in the 3 − 3 − 1 models is relatively straightforward [20,28]. In particular, in Ref. [20] a gravitationally stable PQ mechanism for the model considered here is successfully implemented. We are going to review its main results for completeness.
First of all, we search for all U (1) symmetries of the Lagrangian given in Eqs. (7) and (11). Doing so, we find only two symmetries, U(1) X and U(1) B , which clearly do not satisfy the two minimal conditions required for the U (1) PQ symmetry. See Table I for the quantum number assignments of the fields for these symmetries. In other words, the U(1) PQ is not naturally allowed by the gauge symmetry. However, if the Lagrangian is slightly modified Table I: The U(1) symmetries of the Lagrangian given by Eqs. (7) and (11).
all terms in V Z 2 (η, ρ, χ) are forbidden. In addition, the Yukawa Lagrangian interactions given in Eqs. (8)(9)(10) are slightly modified to Consequently, with the imposition of this Z 2 symmetry a U(1) PQ symmetry is automatically introduced with the charges given in Table II. As η, ρ, χ get VEVs, an axion appears in the physical spectrum. However, it is a visible axion because the U(1) PQ symmetry is actually broken by v ρ 0 2 , which is upper bounded by the value of v SM 246 GeV, as shown in Refs. [20,37]. Hence, this scenario is ruled  [42]. Nevertheless, a singlet scalar, φ ∼ (1, 1, 0), can be introduced in order to make the axion invisible. Its role is to break the PQ symmetry at an energy scale much larger than the electroweak one. This field does not couple directly to quarks and leptons, however it couples to the scalar triplets, η, ρ and χ, through Hermitian terms and the non-Hermitian term λ PQ ijk η i ρ j χ k φ, from which it gets a PQ charge equal to 6, cf. Table II. Notice that this term is allowed as long as the φ field is odd under the Z 2 symmetry, i.e., Z 2 (φ) = −φ.
Although the Z 2 discrete symmetry apparently introduces the PQ mechanism in the model, there are two issues with it. First, the Z 2 and gauge symmetries allow some renormalizable terms in the scalar potential, such as φ 2 , φ 4 , ρ † ρφ 2 , η † ηφ 2 , χ † χφ 2 , that explicitly violate the PQ symmetry in an order low enough to make the PQ mechanism ineffective. Second, since the PQ symmetry is global, it is expected to be broken by gravitational effects [43,44]. Thus, a mechanism to stabilize the axion solution has to be introduced.
As usual, the entire Lagrangian is considered to be invariant under a Z D discrete gauge symmetry (anomaly free) [20,28,[45][46][47][48] and, in addition, this symmetry is supposed to induce the U(1) PQ symmetry. For Z D≥10 it is found that all effective operators of the form

Pl
(where N ≥ D is a positive integer and M Pl is the reduced Planck mass) that can jeopardize the PQ mechanism are suppressed. In particular, in Ref. [20] two different symmetries, Z 10 and Z 11 , were found to stabilize the PQ mechanism for the Lagrangian given by Eqs. (11,(15)(16)(17). The specific charge assignments for these symmetries are shown in Table III. Note that the term λ 15 ijk η i ρ j χ k in the scalar potential is prohibited by both of these discrete symmetries and it must be removed from the entire Lagrangian.
We remark that both the Z 10 and Z 11 discrete symmetries in Table III are anomaly free.
This type of discrete symmetry is known as gauge discrete Z N symmetry and it is assumed to be a remnant of a gauge (local) symmetry valid at very high energies, [45]. The anomalyfree conditions are necessary in order to truly protect the PQ mechanism against gravity effects [46,[49][50][51], Specifically, these discrete symmetries satisfy respectively. Other anomalies, such as Z 3 N , do not give useful low energy constraints because these depend on some arbitrary choices concerning to the full theory. In particular, the Z 3 N anomaly depends on the fermions which get masses at very high energy and are integrated out in the low-energy Lagrangian. All the details of these anomaly conditions applied to the 3 − 3 − 1 model can be found in Ref. [20].
In both cases, the axion, a (x), is the phase of the φ field, i.e., φ (x) ∝ exp ia (x) /f a , which impliesf a ≈ v φ . As it is well known, to make the axion compatible with astrophysical and cosmological considerations, the axion-decay constant f a (related tof a by f a =f a /N C =f a /N DW , with N DW being the number of domain walls in the theory. In this model we have N C = N DW = 3), must be in the range 10 9 GeV f a 10 12 GeV (we are assuming a post-inflationary PQ symmetry breaking scenario). Note that this high , justifies the approximation in the form of axion eigenstate. It is also important to remember that in this model v 2  Table III. These two effects induce an effective potential for the axion, V eff , from which it is possible to determine the axion mass.
In more detail, as the U(1) PQ symmetry is anomalous, we will have a V PQ term in the effective potential, which can be written as where m π 135 MeV and f π 92 MeV are the mass and decay constant of the neutral pion, respectively; m u and m d are the masses of the up and down quarks. Note that V PQ has a minimum when a (x) /f a = 0, which solves the strong CP problem in the usual way.
However, because of the PQ symmetry is also explicitly broken by gravity effects, the effective potential gets another term, V gravity , which reads where N = 10, 11 for Z 10 and Z 11 , respectively. The phase δ D inside the trigonometric function can be written as where δ is the phase of the g coupling constant andθ is the parameter which couples to the gluonic field strength and its dual. This extra term in the scalar potential, Eq. where V eff has a minimum. Expanding V eff = V PQ + V gravity in powers of a(x) fa , we find that in the minimum, the axion VEV satisfies where we have used v φ ≈f a = N DW f a . Note that for |g| = 0 (or for δ D = 0) we have that a(x) fa = 0 in the minimum, as it should be to solve the strong CP problem. However, in the general case, the value of a(x) fa does not satisfy the NEDM constraint [8], which imposes In addition, V gravity brings a mass contribution for the axion, m a, gravity . From Eq.
This contribution can, in general, be much larger than the well-known axion-mass term coming from the QCD nonperturbative terms, Eq. (18), Thus, in order to maintain the axion mass stable, we are going to look for values of the parameters |g|, f a and δ D for N = 10, 11 that both satisfy the NEDM constraint and leave the axion mass stable (m a, QCD m a, gravity ).
Before closing this section, it is important to remark that although the 3 − 3 − 1 model considered in this paper has additional contributions to CP -violating processes that in principle can contribute to the NEDM, these do not require tuning the model parameters at the same order of the θ parameter as it was correctly estimated in Ref. [20]. Roughly

TER
For the postinflationary f a values considered here, cold dark matter in the form of axions can be produced by three different processes: the misalignment mechanism [52], where the axion field oscillates about the minimum of its potential, trying to decrease the energy after the breaking of the PQ symmetry; and the decay of one-dimensional (global strings [53]) and two-dimensional (domain walls [54]) topological defects, which appear after breaking this symmetry. Now, we will briefly review the general expressions for the axion relic density in these three mechanisms following Ref. [55].

A. Misalignment mechanism
The equation of motion for the axion field a in a homogeneous and isotropic Universe, is of the type of a damped harmonic oscillator with a natural frequency equal to the axion mass. In this case, taking into account nonperturbative effects of QCD at finite temperature and considering the interacting instanton liquid model (IILM) [56], the axion mass depends on the temperature as [57] where the values of the parameters are c T = 1.68 × 10 −7 , n = 6.68 and Λ QCD = 400 MeV [57]. This dependence, is valid in the regime where the axion mass at temperature T is less than its value at temperature zero, given by m a (0) 2 = c 0 , where c 0 = 1.46 × 10 −3 , which leads to a minimum temperature ∼ 103 MeV for the validity of the fit. The temperature T osc at which the axion field begins to oscillate is given by [55] T osc = 2.29 GeV g * (T osc ) 80 where g * (T osc ) is the number of relativistic degrees of freedom at temperature T osc . Eq. (26) is valid for temperatures greater than 103 MeV, where Eq. (25) holds, and it is also assumed a not too strong dependence on the temperature of g * , which, for the range 10 9 GeV < f a < 10 12 GeV analyzed in this work, varies between 80 and 85 [58], what would change the abundance of axion dark matter by a factor of ≈ 1.02. Once the adiabatic condition is satisfied, both the entropy and the number of axions with momentum zero per comoving volume are conserved [9], and it is possible to obtain the dark matter abundance [55] Ω a,mis h 2 = 4.63 × 10 −3 f a 10 10 GeV 6+n 4+n , where g * (T osc ) = 80 and Λ QCD = 400 MeV have been used.

B. Decay of global strings
Global strings are the first of the topological defects that appear after the breaking of the U(1) PQ symmetry at T v φ because the field φ (with PQ charge equal to 6 in the 3 − 3 − 1 model considered here) acquires a VEV | φ | = v φ [55,59]. Actually, the breaking of the PQ symmetry leads to the formation of a densely knotted network of cosmic axion strings, which oscillate under their own tension, losing their energy by radiating axions [60]. The radiation process lasts from the PQ-symmetry breaking time to the QCD phase transition time. Using results of numerical studies which provide the time dependence of ρ string (energy density of strings) and ρ a,string (energy density of axions produced by the string decays), it is possible to obtain the nowadays abundance of radiated axions [61,62], with α = (7.3 ± 3.9) × 10 −3 , g * (T osc ) = 80 and Λ QCD = 400 MeV. N DW = 3 is the number of domain walls in this model, and n = 6.68 is the same parameter that appears in Eq. (25).

C. Decay of string-wall systems
In that these domain walls decay at a certain time after being formed [63]. Actually, the domain walls bounded by strings begin to oscillate and eventually, when their tensions are greater than the tensions of the strings, their annihilations lead to axion production [64,65].
The energy density of domain walls can overclose the Universe due to its dependence on the inverse of the square of the scale factor, R, which decreases at a slower rate than the corresponding to matter, ρ ∼ R −3 , and radiation, ρ ∼ R −4 . In our case, this problem is solved by the introduction of a Planck-suppressed operator in the effective potential for the axion field a, parametrized as in Eq. (19).
The current axion abundance is given by the expression [55,66]: where N −4 , and β = 1.65 ± 0.47 is a parameter obtained from numerical simulations. Finally, we will refer to the case p = 1 as the exact scaling, and p = 1 as the deviation from scaling. From here on, we use p = 0.926 for the deviation from scaling case, since it is the suggested value by numerical simulations [55].
In order to conclude this section, we have seen that axions can be produced by three different non-thermal mechanisms, which leads to the result that the total abundance of axions in the Universe can be written as the sum of all these contributions, Eqs. (27), (28) and (29), i.e., The total dark matter abundance due to axions is upper bounded by the observational constraint on the current relic density Ω Planck DM h 2 = 0.1197 ± 0.0066 (at 3σ) as reported by the Planck Collaboration [67]. In the next section, we will analyze the behavior of each contribution to the total abundance, in order to establish a suitable region of parameters for the model analyzed in this work.

MATTER
In general, the total dark matter relic density due to axions in this 3 − 3 − 1 model depends on f a , g, N DW and Z N . The dependence on f a , g, and N DW is direct because Ω a,mis , Ω a,string and Ω a,wall explicitly depend on these parameters. Nevertheless, the dependence on Z N is indirect. Roughly speaking, this discrete symmetry constrains the order of the dominant gravity-induced operator gφ N /M N −4 Pl . In other words, the discrete symmetry sets the exponent N which directly affects the total dark matter due to axions. Actually, we have two discrete symmetries, Z 10 and Z 11 (see Table III), that stabilize the PQ mechanism, which implies that there are two cases to be considered, N = 10 and N = 11. On the other hand, the domain wall parameter, N DW , is set to be equal to 3 by the PQ symmetry and the matter content in the model. Thus, we are interested in knowing if the model with Z 10 or/and Z 11 symmetry provides the total dark matter reported by the Planck collaboration [67] when f a , g, take their allowed values, without conflicting with the constraints on the axion phenomenology.
In order to do that, it is convenient, first, to study separately the behavior of the three axion production mechanisms which results are shown in Fig. 1. Specifically, the cyan and black lines show the axion abundances produced by misalignment and global string decay mechanisms, respectively. On the other hand, the blue lines show the abundance of axion dark matter due to the decay of domain wall systems for N = 10 and N = 11, calculated for the coupling constant value |g| = 1. Two shaded regions are also shown: the light red one corresponds to the exclusion region coming from the constraint of the over closure of the Universe [67], and the yellow region gives the possible interval for the axion decay constant f a , for which no over abundance of axions from decay of global strings or domain walls is produced. Finally, the dark green line corresponds to the total abundance of axions, Ω a h 2 , as given by Eq. cause Ω a,string has an extra N 2 DW = 9 global factor. Indeed, the misalignment mechanism contributes at most by ≈ 7% for the total dark matter density. In contrast, Ω a,wall is decreasing with f a and thus it dominates Ω a for the smaller values of f a (3.6 × 10 9 GeV f a 5.3 × 10 9 GeV). That can be understood realizing that the domain-wall time decay is larger for smaller f a values, making the domain wall more stable and, in this way, explaining why this mechanism contributes more for the axion relic density when f a is smaller.
The opposite behavior of Ω a,string and Ω a,wall allow to set an upper and lower bound on f a . For |g| = 1, f a is constrained to be 3.6 × 10 9 GeV < f a < 1.7 × 10 10 GeV in order to satisfy Ω a,wall h 2 and Ω a,string h 2 Ω Planck DM h 2 [67]. Actually, the interval of allowed f a values is slightly thinner because all of the three axion production mechanisms contribute simultaneously. Also, note that the f a upper bound above is independent on the value of N and on the value of |g|, as can be seen from Eq. (28). In contrast, the lower bound is only valid for the case of N = 10. Actually, the case of Z 11 is completely ruled out and, for this reason, our analysis will be concerned exclusively with the Z 10 symmetry case. Once we have gained a general knowledge about the behavior of Ω a h 2 as function of f a for |g| = 1, we can go further studying the parameter space for the Z 10 case, allowed by the axion phenomenology. In particular, in Fig. 2  The blue curves correspond to the regions where the total axion dark matter abundance is equal to Ω Planck DM h 2 , taking into account the uncertainties in the parameters α and β in Eqs. (28) and (29). Notice that for a given value of f a , |g| is lower bounded by these lines. Larger values of |g| imply Ω a h 2 < Ω Planck DM h 2 . The light blue shaded region is ruled out by the over closure of the Universe for the case of the parameter β = 2.12 in Eq. (29) and for the α = 7.3 − 3.9 = 3.4 factor in Eq. (28). From the remaining region, it is possible to exclude another large part applying the axion mass stability condition, m a,QCD > m a,gravity [see the photons, g aγγ , depends on the f a decay constant, the electromagnetic and color anomaly coefficients, E and N C , respectively. It is known that these anomaly coefficients are completely determined by the fermion content and the U(1) PQ charges of the model, cf. Table (II). Standard calculations for anomaly coefficients [41,68]  Moreover, it is notable that for the range with larger masses (blue line), the axion parameters of this 3 − 3 − 1 model are very close to the projected region which is going to be explored by the IAXO experiment [66,69]. ALPS-II [70], of the helioscope IAXO [69], of the haloscopes ADMX and ADMX-HF [71,72]. The yellow band corresponds to the generic prediction for axion models in QCD. In addition, the two (one) thick red (blue) lines stand for the predicted mass ranges and coupling to photons in this model, for |g| = 0.1 (|g| = 1), where axions make up the total DM relic density.

VI. CONCLUSIONS
In this work, we consider a version of an alternative electroweak model based on the SU(3) L ⊗ U(1) X gauge symmetry, the so called 3 − 3 − 1 models, when the color gauge group is added. For this version, which includes right-handed neutrinos, it is shown in Ref. [20] that the PQ mechanism for the solution of the strong CP problem can be implemented. In this implementation, the axion, the pseudo Nambu-Goldstone boson that emerges from the PQ-symmetry breaking, is made invisible by the introduction of the scalar singlet φ ∼ (1, 1, 1) whose VEV, v φ ≈f a , is much larger than v SM , and any other VEV in the model. Moreover, the axion is also protected against gravitational effects, that could destabilize its mass, by a discrete Z N symmetry, with N = 10, 11.
Once we have set this consistent scenario, we investigate the capabilities of this axion, produced in the framework of this particular 3 − 3 − 1 model, to be a postinflationary cold dark matter candidate. We started focusing in the axion-production mechanisms. As it was explained in the previous section, from Fig. 1 we see that the vacuum misalignment mechanism does not dominate the DM relic abundance, and, if it was the only production mechanism in action, an upper bound for f a could be set by imposing that it should account for all the DM abundance, i.e., Ω a,mis h 2 = Ω Planck DM h 2 , and we would find the corresponding value f a ≈ 1.5 × 10 11 GeV, for the parameters determined by the model, in this case N DW = 3. However, there are two other more efficient mechanisms due to the decay of topological defects: cosmic strings and domain walls. As the curves for Ω a,string h 2 and Ω a,wall h 2 grow in opposite directions, relatively to the f a values, we can determine an upper bound and a lower bound for f a by imposing the total Ω a h 2 matches the observed Planck results. This is the case when we add up all the contributions for N = 10, and we find 3.6 × 10 9 GeV < f a < 1.7 × 10 10 GeV. However, we would like to stress that this is not the case for N = 11. For N = 11 there is no value of f a for which the addition of the partial abundances lies below the observed result. It means that the Z 11 , which possesses the good quality of stabilizing the axion, is not appropriate for the axion-production issue since it makes the domain wall mechanism too efficient and overpopulates the Universe.
As it can be seen from Fig. 1, for any fixed allowed value of |g|, there are two values of f a that are in agreement with the value of Ω Planck DM h 2 . In fact they are regions, if we take into account the uncertainties following the discussion in the previous section for Fig. 2.
Outside these regions, the axion abundance will be a fraction of Ω Planck DM h 2 . See the solid dark green curve in Fig. 1 for |g| = 1. If this happens to be the case, i.e., if these predicted regions are somehow excluded, by future experimental data for the axion mass value, for instance, then, another kind of DM will be needed. We have also found special values for δ D , (0.4 − 4.1) × 10 −5 , by requiring the minimal compatible intersection region between the curves that obey the NEDM and Ω Planck DM h 2 constraints. This value was obtained considering the maximum value of |g|, i.e., |g| = √ 4π, cf. Fig. 2(a). For lower values of |g|, higher tuning on δ D is required. However, it seems unnatural to require severe levels of tuning on δ D , since for this quantity a tiny value is the result of the difference between two terms that have completely different origins.
Regarding the capabilities of detecting the axion dark matter, Fig. 3 shows the sensitivities of several experiments in the m a − g aγγ plane. In this plot, the thick blue and red lines are the regions where the axion abundance is responsible for all the observed DM. These lines were obtained by using |g| of order one. Moreover, the blue region corresponding to masses of the order of meV and g aγγ ≈ 10 −12 GeV −1 , lies very close to the projected IAXO sensitivity, so that it will be reachable in the near future.
Looking back to our results we can conclude that this version of the 3 − 3 − 1 model, concerning the axion DM issue and the strong CP problem, is phenomenologically consistent.
This model, besides its good qualities presented in the introduction, also possesses new degrees of freedom that are not yet experimentally probed. For instance, the model has charged and neutral scalars (besides the Higgs), extra vector bosons and extra quarks, that are expected to be heavy, and could, in principle, be searched at colliders. See Refs. [73,74] for recent studies concerning the 3 − 3 − 1 model phenomenology, in general, at the LHC.