The Supersymmetric R-Parity Violating Dine-Fischler-Srednicki-Zhitnitsky Axion Model

We propose a complete R-parity violating supersymmetric model with baryon triality that contains a Dine-Fischler-Srednicki-Zhitnitsky axion chiral superfield. We parametrize supersymmetry breaking with soft terms, and determine under which conditions the model is cosmologically viable. As expected we always find a region of parameter space in which the axion is a cold dark matter candidate. The mass of the axino, the fermionic partner of the axion, is controlled by a Yukawa coupling. When heavy [${\mathcal O}$(TeV)], the axino decays early and poses no cosmological problems. When light [$\mathcal{O}$(keV)] it can be long lived and a warm dark matter candidate. We concentrate on the latter case and study in detail the decay modes of the axino. We find that constraints from astrophysical X- and gamma rays on the decay into photon and neutrino can set new bounds on some trilinear supersymmetric R-parity violating Yukawa couplings. In some corners of the parameter space the decays of a relic axino can also explain a putative 3.5 keV line.


Introduction
Supersymmetry (SUSY) is the best candidate to explain the hierarchy problem of the Standard Model (SM) [1,2], and thus a leading candidate for the discovery of new physics effects at the LHC. Most supersymmetric models which have been searched for at the LHC make the assumption of conserved R-parity [3]. A very restrictive and widely considered version with universal boundary conditions at the unification scale, the constrained minimal supersymmetric Standard Model (CMSSM) [4], is now in considerable tension with collider data [5]. The R-parity violating (RpV) CMSSM is still very much allowed [6].
R-parity is a discrete multiplicative symmetry invoked in the CMSSM to forbid baryon and lepton number violating operators which together lead to rapid proton decay. The bonus of imposing R-parity is that the lightest supersymmetric particle (LSP), typically the neutralino, is stable and provides a dark matter (DM) candidate in the form of a weakly interacting massive particle (WIMP). This is the most popular and most searched for DM candidate. However, to-date no neutralino dark matter has been found [7][8][9]. It is thus prudent to investigate other DM candidates, also within alternative supersymmetric models.
R-parity does not forbid some dimension-5 operators dangerous for proton decay [10]. This issue can be resolved by imposing a Z 6 discrete symmetry, known as proton hexality (P 6 ), which leads to the same renormalizable superpotential as R-parity [11,12], but is incompatible with a GUT symmetry. Alternatively, one can impose a Z 3 -symmetry known as baryon-triality (B 3 ) [11,[13][14][15][16][17][18]. The latter allows for lepton number violating operators in the superpotential thanks to which the neutrinos acquire Majorana masses, without introducing a new very high energy Majorana mass scale [15][16][17][19][20][21][22][23][24]. This is a virtue of B 3 models. 1 A possible handicap is that the LSP is unstable and is not a dark matter candidate. This can be naturally resolved when taking into account the strong CP problem, as we discuss in the following.
Besides the hierarchy problem, every complete model should also address the strong CP problem [30], which plagues the SM, as well as its supersymmetric generalizations. In its simplest forms this necessitates introducing the pseudo scalar axion field [31,32]. In the supersymmetric versions, the axion is part of a chiral supermultiplet and is accompanied by another scalar, the saxion, and a spin-1/2 fermion, the axino.
In this paper we propose to study a B 3 model with the inclusion of an axion supermultiplet of the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) type [33,34]. The axion contributes to the DM energy density in the form of cold DM [35,36]. Depending on the value of its decay constant, f a , it constitutes all or a fraction of the DM. The gravitino is heavy [O(TeV)], decays early in the history of the universe, and does not pose cosmological problems. The axino mass is proportional to λ χ v χ [see Eq. (2.18)], where v χ is the vacuum expectation value (VEV) of a scalar field of the order of the soft masses [O(TeV)], and λ χ is a dimensionless Yukawa coupling in the superpotential. We consider two cases: (i) a heavy axino (∼ TeV, which requires λ χ ∼ 1), and (ii) a light one (∼ keV, requiring λ χ 1). In case (i) the axino decays early and is not a DM candidate. In the more interesting case (ii), it is the LSP, its lifetime is longer than the age of the universe, and contributes as warm DM [37]. We study in detail its three decay channels: into an axion and a neutrino, into three neutrinos, into a neutrino and a photon, the last one subject to constraints from Xand gamma-ray data.
In 2014 a potential anomaly was observed in X-ray data coming from various galaxy clusters and the Andromeda galaxy [38,39]. There have been several papers discussing this in terms of an axino in R-parity violating supersymmetry [40][41][42]. Although we were able to show that this explanation is not probable [43].
The outline of the paper is as follows. In Sec. 2, we introduce the model, describe its parameters and mass spectrum. In Sec. 3 we compute in detail the axino decay rates and branching fractions. In Sec. 4 we consider cosmological and astrophysical constraints on the model, with a focus on the more interesting case of a light axino. We conclude with a discussion of our results and point out possible future directions of investigation.
In Appendix A we include some general comments on the axino mass, we discuss its dependence on the SUSY breaking scale and point out that DFSZ models accommodate more easily a light axino, as opposed to modelsàla Kim-Shifman-Vainshtein-Zakharov (KSVZ) [44,45]. In Appendix B we furthermore give details on the derivation of the mixings of the axino mass eigenstates.

The Lagrangian
The building blocks of the minimal supersymmetric DSFZ axion model are the particle content of the MSSM plus three superfields, which are necessary to construct a self-consistent axion sector as discussed in Ref. [46]. The particle content, the quantum numbers with respect to the SM gauge sector SU (3) c ×SU (2) L ×U (1) Y , and the charges under the global Peccei-Quinn (PQ) symmetry [47], U (1) PQ , are summarized in Table 1. The renormalizable superpotential is given by

2)
(2.4) Here and in the following we suppress generation, as well as isospin and color indices. Superfields are denoted by a hat superscript. We have imposed the discrete Z 3 symmetry [11] known as baryon triality, B 3 . This leaves only the RpV operators 1 2 λLLÊ, λ LQD , and c 2ÂLĤu . We have generalized the usual bilinear operators µH u H d and κ i L i H u to obtain PQ invariance. The PQ charges satisfy from the terms in Eq. (2.3), and similar relations for the other fields that are readily obtained from the terms in the rest of the superpotential. The PQ symmetry is spontaneously broken at the scale f a , due to scalar potential resulting from W PQ , and the scalar parts ofÂ andÂ get vacuum expectation values (VEVs) v A , vĀ, respectively, and we denote their ratio as Below the PQ breaking scale, we generate an effective µ-term, µ effĤuĤd , and effective bilinear RpV κ-terms, κ eff,iLiĤu , with For successful spontaneous breaking of the electroweak (EW) symmetry, µ eff should be of the order of the EW scale, while κ eff is constrained by neutrino physics, as discussed below, and should be ≤ O(MeV). Since f a > 10 9 GeV from astrophysical bounds on axions [48], the couplings c 1 and c 2 must be very small, roughly c 1 < 10 −6 and c 2 < 10 −11 . They are, nonetheless, radiatively stable. We have assumed that W PQ respects an R-symmetry, under which the fieldχ has charge 2. This forbids the superpotential termsχ 2 andχ 3 . See also Ref. [25].
The soft-supersymmetric terms of this model consist of scalar squared masses, gaugino masses, and the counterparts of the superpotential couplings. The full soft Lagrangian reads with f ∈ { e, , d, u, q}. We assume that SUSY breaking violates the R-symmetry of the axion sector, hence we include the soft terms T λχ χAĀ and L V χ. Note that, even if set to zero at some scale, these terms would be generated at the two-loop level because of the small coupling to the MSSM sector, where no R-symmetry is present. At low energy the Lagrangian contains interactions of the axion with the gauge fields. These are best understood by performing a field rotation 2 that leads to a basis in which all the matter fields are invariant under a PQ transformation [50]. We have Here a is the axion field which in our model is given by a 1 √ 2 (σ A − σĀ) to a good approximation 3 , cf. Eq. (2.6). B µν , W a µν , G a µν are the SM field strength tensors, here a is an adjoint gauge group index, with X aµν ≡ µνρσ X a ρσ , and the coefficients are given in terms of the PQ charges In the supersymmetric limit the couplings of the axino to gauginos and gauge fields are related to those in Eq. (2.11). We are, however, interested in the case of broken SUSY.
In general, such axino couplings can be calculated explicitly, by computing triangle loop diagrams, once the full Lagrangian is specified, as it is in our model. Later in the paper we use the full Lagrangian to compute the one-loop decay of the axino into a photon and a neutrino, including the couplings.

The Mass Spectrum
We now turn our attention to the spin-1/2 neutral fermion mass spectrum of this model. The 10 × 10 mass matrix in the basis (λ B , W 0 , H 0 u , H 0 d , ν L,i , A, Ā , χ), which is a generalization of the MSSM neutralino mass matrix, reads (2.14) M χ 0 was obtained by first implementing our model with the computer program SARAH [51][52][53], then setting the sneutrino VEVs to zero, which corresponds to choosing a specific basis [54,55]. This is possible, since the superfieldsĤ d andL i have the same quantum numbers in RpV models. The choice of basis is scale dependent, in the sense that sneutrino VEVs are generated by renormalization group (RG) running, even if set to zero at a given scale. Hence, one needs to specify at what scale the basis is chosen. For later convenience, when we compute the axino decays, we choose this scale to be the axino mass. Neglecting for a moment the axino sector, which couples very weakly to the rest, and Takagi diagonalizing the upper-left 7 × 7 block [56], one finds [20,57] two massless neutrinos and a massive one with Requiring m ν ≈ 0.1 eV, with tan β ≈ 1 and µ eff ≈ 1 TeV, one finds that κ eff,i is of order MeV, and correspondingly smaller for larger tan β. The two neutrinos which are massless at tree level acquire a small mass at one loop [19,23,24], which we have not included in the 10x10 mass matrix above.
In addition to the neutrinos, there are four eigenstates, mainly built from the MSSM neutralinos (λ B , W 0 , H 0 u , H 0 d ), with masses between 100 and 1000 GeV, and three axino eigenstates with masses that can be calculated analytically in the limit v A = vĀ = fa √ 2 , corresponding to tan β = 1, and c 1,2 → 0. One finds for the axinos [46] The lightest of these three states corresponds to the fermionic component of the linear combination of superfields 1 √ 2 (Ā − A). It is interpreted as the axino with a mass [46] 4 For the more general case of tan β = 1, we can not give analytical expressions. Instead, we find at the lowest order in We also give the power of the next term in the expansion, where f a enters. This is the main contribution to the axino mass, we neglect small corrections proportional to c 1 and c 2 . We expect v χ to be at the soft SUSY breaking scale, O(TeV). The Yukawa coupling λ χ then is our main parameter to control m a . λ χ is radiatively stable: if we set it small at tree level, it remains small. See also Ref. [58] for a quantitative treatment of Yukawa RGEs inn supersymmetric models. The most interesting case we consider later is that of a light axino, with m a = O(keV), which means λ χ must be very small.
In order to determine the interactions of the axino, in particular the potential decay modes, we would like to estimate the mixing between the axino and the neutrinos. To do this one first diagonalizes the upper 7x7 block, neglecting the axinos which are very weakly coupled. The massive neutrino eigenstate will be mostly a combination of the gauge eigenstates ν L,i , with small components of bino, wino and higgsinos. Next include the lower 3x3 axino block. The off-diagonal entries between the neutrino and axino blocks are Now rotate the 3x3 axino block to the mass eigenstates. The axino is a linear combination of A andĀ, the other two heavier fermions are a combination of A,Ā and χ. The important point is that the entries in the 3x3 Takagi diagonalization axino matrix are of order one. This implies that the pieces in the off diagonal blocks remain of order κ eff v u /f a after rotating the axino block, and we obtain the form of the matrix From this we can read off the axino-neutrino mixing For the purpose of this estimate and those in Sec The axino-higgsino mixing can be estimated by observing that the higgsino mass is of order µ eff , and the off-diagonal element between axino and higgsino is of order , we can estimate the bino-higgsino mixing as g 1 v M SUSY . Multiplying this times x a, H we obtain the axino-bino mixing Next we briefly consider the scalar sector. It is instructive to start by considering the limit c 1,2 → 0, in which there is no mixing between the MSSM sector and the axion sector, and the axion block of the scalar squared mass matrix can easily be diagonalized. We find a massless axion and a saxion with a mass squared ∼ m 2 A or m 2Ā , the parameters in the soft Lagrangian [m A , mĀ = O(TeV), cf. Eq. (2.10)]. They are, respectively, the CP-odd and the CP-even mass eigenstates, corresponding to the combination 1 √ 2 (Ā − A). The result is not appreciably altered once we turn on the small couplings c 1 and c 2 . The other four real scalar degrees of freedom in this sector are heavy, with masses of order λ χ f a . The masses of the scalars in the MSSM sector are to a very good approximation those already well studied in the RpV SUSY literature [19,58].

Light Axino Decay Modes
In this section we consider an axino mass lower than twice the electron mass, in which case the axino has only three three decay modes: 1. into a neutrino and an axion (via a dimension 5 operator),ã → ν i + a; 2. into three neutrinos (tree-level),ã → ν i + ν j + ν k ; 3. into a neutrino and a photon (one-loop),ã → ν i + γ.
We discuss them in turn.

Decay into Neutrino and Axion
The dominant operator responsible for the decayã → a + ν is dimension 5: This is in the basis of the mass matrix in Eq. (2.20), before we diagonalize to the final mass eigenstate. Here a is the axion, ψã is the four-component spinor denoting the axino mass eigenstate and the neutrino arises due to mixing with the axino, once we diagonalize.
A simple way to understand the origin of this operator is to consider the kinetic term ψãγ µ ∂ µ ψã. After the chiral rotation ψã → e iγ 5 a/fa ψã, we obtain the above operator, when we expand the exponential. The diagram for the decay is shown in Fig. 1a, where we include the axino-neutrino mixing [see Eq. (2.21)] in the final state. The partial decay width into this channel is

Decay into Three Neutrinos
This tree-level decay proceeds via the axino mixing with the neutrino and the exchange of a Z boson between the fermionic currents, as shown in Fig. 1b. Assuming for simplicity that the neutrino admixture of the axino is measured by x a,ν for all neutrino flavors, the decay width into all possible flavors of final neutrinos reads [59] (3.5) Γã →3ν and Γã →νa give the main contribution to the total decay width in most of the parameter space. We find that for 10 10 GeV < f a < 10 12 GeV and 1 keV < m a < 100 keV, the axino lifetime ranges from 10 23 to 10 41 s, which is longer than the age of the universe (∼ 4 · 10 17 s). Figure 1: Left: Feynman diagram for the axino decay into an axion and a neutrino, as described by the five-dimensional operator of Eq. (3.1). Right: One of the Feynman diagrams describing the decay of the axino into three neutrinos. In our convention the ⊗ on a fermion line stands for a mixing between fermionic eigenstates.

Decay into Photon Plus Neutrino
The one-loop diagrams for the decay a(p) → ν(k 1 ) + γ(k 2 ) have an amplitude with the following structure [60] dictated by gauge invariance Here η ν and η a are the signs of the mass eigenvalues of the neutrino and the axino, respectively, k 2 the photon momentum, and its polarization vector. g aνγ is a function with the details of the loop integrals, and has the dimension of inverse mass. The decay rate then has the form We are neglecting the final state neutrino mass. The corresponding Feynman diagrams are shown in Figs. 2 and 3. Here we are interested in an estimate, and thus wish to determine which diagram(s) is (are) dominant. Thus we first consider the individual amplitudes squared, in turn, as well as the corresponding decay width, without including interference terms. If not otherwise specified, we make use of the formulae provided in Ref. [61] to compute the loop integrals. We do not include diagrams with sneutrino VEVs insertions (see e.g. Ref. [62]), since, as previously mentioned, we work in the basis, defined at the axino mass scale, where such VEVs are zero. The first diagrams we consider are those depicted in Fig. 2a , 2b. The corresponding decay width has been computed in the analogous case of a sterile neutrino decaying into a neutrino and a photon [63] (see also Ref. [64] for a review), Here, the first subscript on Γ indicates the state the incoming axino mixes with, while the subscripts in parentheses denote the particles running in the loop. We use this notation for Figure 2: Diagrams with no dependence on the the RpV trilinear couplings λ ( ) . The cross on a fermion line indicates a mass insertion needed for the chirality flip, while the ⊗ on a scalar or fermion line is a (dimensionless) mixing. We write the coupling constants in green at the vertices.
the rest of the section. Let us anticipate that in most of the parameter space we consider, Γ ν(W l) gives the dominant contribution. It is therefore useful to compare the other partial widths to Γ ν(W l) . Incidentally, a decay width of Γã = 3·10 −58 GeV corresponds to a lifetime τã = 2 · 10 23 s and the lifetime of the universe is about τ univ. ≈ 4 · 10 17 s. The main contribution to the two diagrams in Fig. 2c, 2d comes from the heaviest lepton that can run in the loop, the τ . We have The m 2 τ in the numerator comes from the mass insertion that flips the chirality in the fermion line in the loop, while Bµ/m 2 H ± is an estimate of the mixing between the charged scalar Higgses, Bµ being the bilinear parameter in the Higgs potential, see Eq. (2.10). Taking Bµ ∼ m 2 H ± as an estimate, we find the ratio which is highly suppressed, and we neglect the contributions from the corresponding diagrams in the following. Next we consider the diagrams of Fig. 2e, 2f. We assume theτ to be the lightest charged slepton and give the leading contribution. We find Here, xτ LR is the mixing between the left and right staus, µ eff is the Higgsino mass. Note that the expression is finite in the limit µ eff = mτ . Taking xτ LR ≈ 1, we find the ratio . This is also highly suppressed and we neglect it in the following. Next, we examine the diagrams of Fig. 3, which involve at least one trilinear RpV coupling, λ or λ . To avoid cluttering the equations, we use λ to denote either. For the diagrams of Fig. 3a , 3b we find Here, f ( f ) denotes either a down-type quark (squark), in the case of a λ coupling, or a lepton (slepton), in the case of a λ coupling. Q f is the electric charge of f in units of e, and g 1 = e/ cos θ W is the hypercharge gauge coupling. The dominant contributions come from loops of bottom-sbottom, or tau-stau; the two are of the same order of magnitude, under the assumption that sleptons and squarks have comparable masses, as the smaller mass of τ is compensated by its larger electric charge. The axino-bino mixing is estimated Figure 3: Diagrams with an explicit dependence on the RpV trilinear couplings λ ( ) . The notation is as in Fig. 2. in Eq. (2.23). Comparing this contribution to the one of Eq. (3.8), we find where m b is the mass of the bottom quark. For sfermions lighter than a TeV, Γ B(ff ) would give a contribution larger than Γ ν(W l) . Here we restrict ourselves to m f > 1 TeV, motivated by recent analyses [6,65], which disfavor lighter squarks and sleptons. We comment further on this contribution below.
For the diagrams of Fig. 3c , 3d we find (note there is no mass insertion in the fermion line in the loop) Here, f ( f ) can be either a down-type quark (squark), in the case of a λ coupling, or a lepton (slepton), in the case of a λ coupling. Y is either the Y d or the Y e Yukawa, and the axino-higgsino mixing is given in Eq. (2.22). As in the previous case, the main contributions to these diagrams are from bottom-sbottom and tau-stau loops, respectively. Γ H(ff ) is suppressed compared to Γ B(ff ) by a factor m 2 a /m 2 f . We have The last two diagrams to evaluate are those of Fig. 3e, 3f. The two RpV couplings appearing in each of these diagrams have no restrictions on their flavor index structure, contrary to those of Figs. 3a to 3d, which are diagonal in the singlet-doublet mixings. For the sake of our discussion we also do not distinguish between the two potentially different λ here, as we consider only one RpV coupling to be different from zero at a time. We find .
From our estimates we see that Γ B(ff ) and Γ ν(W l) give the dominant contributions, the other diagrams being always negligible. Thus the axion decay width into photon and neutrino is given to a good approximation by Γ B(ff ) + Γ ν(W l) . In Fig. 4 we plot the branching ratios corresponding to the three decay channels of our light axino. It is worth examining Γ B(ff ) in more detail. First, we reintroduce the generation indices in the superpotential terms As the main contributions to Γ B(ff ) are from tau-stau and bottom-sbottom loops, the trilinear couplings of interest are those shown in Table 2. The first index of λ ijk and λ ijk  refers, in our case, to the neutrino in the lepton doublet. We see in Fig. 5 that when we consider the decay a → ν 2 γ, with λ 233 and λ 233 saturating the upper bounds in Table 2, the diagrams giving Γ B(ff ) provide the dominant contribution. We can use this fact to set new bounds on these trilinear couplings by considering constraints from X and gamma rays [66,67] on the decays a → ν 2,3 γ. We proceed as follows. First, we assume that the entire dark matter is constituted by the axino, Ω a h 2 = Ω DM . This fixes f a for any given value of m a . Then the excluded region corresponds to where τ bound is given, for instance, in Fig. 9 of Ref. [67] as a function of the mass of the decaying particle, the axino in our case. In Fig. 6 we see that, for m a = 100 keV, we can set bounds on the trilinears λ 233 , λ 233 , λ 333 of order 10 −2 . These are to be compared with the bounds in Table 2. However, it is important to keep in mind that these new bounds, contrary to those in Table 2, rely on the presence of the axino and on the assumption that it is the whole dark matter.
In Fig. 8 Table 2: Upper bounds on the magnitude of R-parity violating couplings at the 2σ confidence level, taken from Ref. [6]. The constraints arise from indirect decays. The concrete processes are described in detail in Ref. [68].   Table 2. excluded regions in purple correspond to Here, Γ a→νγ = max[Γ B(ff ) , Γ ν(W l) ], and we multiply by Ω a /Ω DM to account for the region of parameter space where the axino is under abundant. We discuss the axino relic abundance in Sec. 4.

Comment on the 3.5 keV Line.
This model can fit the 3.5 keV line observed in galaxy clusters [38][39][40]42].  Figure 6: We show regions excluded by X and gamma-ray constraints [66,67]. Here we use Eq. (3.20) and assume the axino to contribute to the totality of dark matter. The dark (light) gray region is excluded when m a = 10 keV (100 keV). In λ i33 , the index i can be either 2 or 3, in which case the dominant contribution to the decay a → ν i γ is from Γ B(ff ) .

The Gravitino
We assume a supergravity completion of our model, with SUSY broken at a high scale ( √ F ∼ 3 · 10 10 GeV) and the breaking effects mediated via Planck suppressed operators, so that the soft scale is M SUSY ∼ F M Pl = O(TeV). The gravitino is heavy, with m 3/2 ∼ M SUSY , decays before BBN, and poses no cosmological problems.

The Saxion
In the context of supersymmetric axion models one usually has to worry about the saxion: if it is too light it can pose serious issues (see e.g. Ref. [78]). We briefly review why. The saxion is a pseudomodulus, meaning that its potential is quite flat. During inflation we expect it to sit at a distance of order f a from the minimum of its zero-temperature potential. After inflation, when the Hubble parameter H becomes comparable to the saxion mass, m s , the saxion starts to oscillate about its minimum. At that time the energy density in the oscillating field is ρ s ∼ m 2 s f 2 a . The energy density of radiation is ρ r ∼ (T osc s ) 4 , so that ρs ρr osc Pl . The energy density of the oscillating saxion scales as T 3 , behaving like matter, which implies that in the absence of decays it would come to dominate the energy density of the universe at a temperature It is important to determine whether it decays before or after domination. To proceed with this simple estimate we parametrize the decay rate as Γ s ∼ 1 16π , neglecting the masses of the decay products. The factor of f 2 a at the denominator is due to the fact that the saxion couplings to matter are suppressed by f a . Then the decay temperature is obtained from H(T dec s ) ∼ Γ s : From these estimates we find

(4.4)
In our model the saxion mass is of order TeV, so it decays way before it comes to dominate the energy density of the universe. 5 Note the temperature when it decays is of order GeV, safely above BBN.

Axion Dark Matter
In this model the axion is a dark matter candidate. We do not review here the calculations of its relic abundance, but refer the reader to Ref. [81,82]. There are two important mechanisms to produce axion cold DM. One is the misalignment mechanism, which results in the abundance Ω a,mis h 2 = 0.18 θ 2 1 f a 10 12 GeV with θ 1 the initial misalignment angle. The other contribution comes from axions emitted by global axionic strings, which also contribute as non-relativistic axions: Here ξ is the length parameter which represents the average number of strings in a horizon volume. It is determined by numerical simulations [83], ξ 1.0. The parameterr is related to the average energy of the axions emitted in string decays. Its value has been the subject of a long debate. Two scenarios have been put forth, one predictsr ≈ 1, the otherr ≈ 70.
See Ref. [84] and references therein for details, and Ref. [85] for recent considerations on the issue. While the misalignment mechanism produces axions when the temperature is T ∼ Λ QCD , the strings form at a much higher temperature T ∼ f a . If the PQ phase transition happened before the end of inflation, or equivalently if the reheating temperature (T RH ) were lower than f a , then the strings would be inflated away and we would no longer have a contribution corresponding to Eq. (4.6). In this scenario the axion abundance is from Eq. (4.5), and θ 1 can take any value between −π and π. In the scenario with the PQ phase transition post inflation (T RH > f a ), one has to average θ 2 1 over many QCD horizons, θ 2 1 = π 2 /3, and the contribution from string decays can be comparable to that from misalignment or larger, depending on the value ofr.

A Heavy Axino
For λ χ ∼ 1, we have m a ∼ v χ ∼ M SUSY , an axino mass of order TeV. Such an axino has many open decay channels into MSSM particles, with a total width that can be estimated as Γ a ∼ 1 8π m 3 a f 2 a . It decays safely before BBN, and does not pose cosmological issues.

A Light Axino
The axino, if light, can also be a dark matter candidate [37,[86][87][88][89][90][91]. It was pointed out in Ref. [37] that for T RH > f a a stable axino should not exceed the mass of 1 or 2 keV, otherwise it would result in overabundant DM. More recently Bae et al. [82,92] showed that the axino production rate is suppressed if M Φ f a , where M Φ is the mass of the heaviest PQ-charged and gauge-charged matter supermultiplet in the model. This suppression is significant in DFSZ models, where M Φ corresponds to the Higgsino mass. The consequence of such a suppression is that the upper bound of 2 keV on the axino mass quoted in Ref. [37] is relaxed. In our model the axino is not stable. However, if its mass is below the electron mass the only decay modes are 1. into a neutrino and an axion (tree-level, with dimension 5 operator),ã → ν i + a; 2. into three neutrinos (tree-level),ã → ν i + ν j + ν k ; 3. into a neutrino and a photon (one-loop),ã → ν i + γ.
The resulting lifetime, as we show in the following subsections, is longer than the age of the universe, thus making the axino a dark matter candidate, with relic abundance Ω a h 2 2.1 × 10 −5 m a 1 keV (4.7) This does not depend on T RH , as long as T RH > 10 4 GeV [82,92]. For lower T RH and different cosmological histories the axion abundance can be different [93], but we will not consider such scenarios in this work. We learn that, for f a = 10 12 GeV, the axino can have a mass as large as a few MeV before it becomes overabundant. The important feature of Eq. (4.7) is that the abundance is inversely proportional to f 2 a . It is easy to understand why. Figure 7: We show, in the plane f a vs the axino mass m a , the regions where the combination of axion and the axino dark matter exceeds the observed abundance, (Ω a + Ω a )h 2 > Ω DM h 2 = 0.12 [94]. In the grey region the main contribution is from the axion, in the red from the axino [Eq. (4.7)]. Along the black solid line the axion alone results in the correct abundance, while along the black dashed line the axino alone accounts for Ω DM . In the left panel the reheat temperature, T RH , is below the PQ phase transition scale, f a , and the only contribution to axion dark matter is from the misalignment mechanism, Eq. (4.5).
Here we take the initial misalignment angle θ 1 = 1. In the white region, the combination of axion plus axino results in a total abundance below Ω DM . In the right panel T RH is above f a , and we also have a contribution to axion dark matter from the decay of axionic strings, Eq. (4.6). We show the excluded region for two different values of the parameterr. Note that forr = 70, axion plus axino are overabundant on the entire plane. Furthermore, the lifetime of the axino is longer than the age of the universe in the entire region plotted.
The axino never reaches thermal equilibrium in the early universe because its interactions are very weak. It is produced in scattering processes (listed e.g. in Table 1 of Ref. [82]) with a cross section proportional to 1/f 2 a . The larger the cross section the more the axino is produced, hence the 1/f 2 a dependence in Eq. (4.7). In Fig. 7 we show the allowed region of parameter space, in the plane f a vs m a , where we can have both axion and axino dark matter.

Summary
We have presented for the first time a complete RpV SUSY model with baryon triality and with a DFSZ axion superfield. We have studied the mass spectrum and then investigated some cosmological bounds. The axion is a good dark matter candidate when its decay constant f a = O(10 11 ) GeV. For lower values of f a the axino, with a mass roughly between 1 and 100 keV, can be the dominant dark matter candidate. We have looked at the possible decay modes of the light axino in detail. For such a light axino, its lifetime is longer than the age of the universe, but its decay into photon and neutrino still leads to interesting phenomenology. We have shown that X-and gamma-ray constraints on this decay give new bounds on some trilinear RpV couplings. These are model dependent and rely on the axino constituting the whole dark matter in the universe. We have also shown that in this model there is a corner of parameter space where a 7 keV axino could fit the 3.5 keV line observed in galaxy clusters. This model could be explored in further details in relation to neutrino physics or possible collider signatures. We leave it to future work. rather than in the KSVZ models [44,45]. We then discuss how the mass is related to the PQ and SUSY breaking energy scales.

A.1 KSVZ vs DFSZ
Axion models fall into two broad categories denoted as KSVZ and DFSZ. Both contain a complex scalar field, A(x), charged under a global U (1) PQ , PQ symmetry [47]. The axion field, a(x), is identified with the phase of A(x). 6 In the KSVZ model, A(x) couples to exotic heavy quark fields, Q(x), that are charged under color and under the U (1) PQ (and possibly U (1) EM ), while the rest of the SM fields do not carry any PQ charge. The coupling f Q AQQ generates, via a triangle diagram, the interaction term aG µν G µν which is crucial in solving the strong CP problem. Here G µν (x) is the gluon field strength tensor and G µν (x) its dual. In the supersymmetric version, the KSVZ model is defined by the superpotential W KSVZ = f QÂQQ , where we use hats to denote chiral superfields and Q,Q are two distinct superfields in the 3 and3 representations of SU (3) c , respectively, but with equal PQ charge. Embedding the model in supergravity results in a one-loop contribution to the axino mass of order [37,99] m a ∼ 10 GeV where m 3/2 is the gravitino mass. Moreover, even if supergravity effects are neglected, there is an irreducible two-loop contribution involving the gluino [100] m a ∼ 0.3 GeV where m g is the gluino mass. This implies that even if the axino mass is tuned to keV values at tree level, loop corrections tend to raise it to GeV values. Thus a KSVZ axino prefers to be heavier than a few keV. In DFSZ models the field A(x) couples to two Higgs doublets, H u and H d , which are also charged under the U (1) PQ , as well as the other SM fields are. Thus the supersymmetric version of the DFSZ axion model is more economical than the KSVZ counterpart, as SUSY already requires at least two Higgs doublets. The SUSY coupling of axions to Higgses can be written as [37] W DFSZ = c 1ÂĤuĤd . The field A has to get a large vacuum expectation value (VEV) A ∼ f a , with f a > 10 9 GeV, otherwise the corresponding axion would be excluded by supernova constraints [48]. In turn, this implies a tiny coupling c 1 . Note that the corresponding coupling in non supersymmetric DFSZ models also has to be very small, for the axion to be invisible. In the SUSY context it has been proposed that the operator c 1ÂĤuĤd could be replaced [101] by a non-renormalizable one g M PlÂ 2Ĥ uĤd . With the latter operator one easily obtains a µ-term at the TeV scale, while with the former we are forced to take very small values the coupling c 1 . We do not address the µ-problem, but we point out that Ref. [102] showed that the operator we consider can be derived consistently within a string theory framework. The fact that c 1 is tiny helps with the axino mass: once we set it to keV values at tree level we are guaranteed that radiative corrections, proportional to c 1 , will be negligible as opposed to the KSVZ case we discussed above.

A.2 Low-scale vs High-scale SUSY Breaking
In a previous paper [46], we indicated a simple way to understand that in models where the SUSY breaking scale, M SB , is lower than the PQ scale, f a , the axino mass is of order O M 2 SUSY fa , with M SUSY the scale of the soft supersymmetry breaking terms, while in models with M SB > f a the axino mass is typically of order M SUSY . The former models, with low M SB are representative of global SUSY, for which the best known framework is gauge mediation of supersymmetry breaking [103]. The latter ones, with high M SB usually fall in the scheme of gravity mediation of supersymmetry breaking, in local supersymmetry, and the scale of the soft terms is that of the gravitino mass, M SUSY ∼ m 3/2 .
In light of these considerations one would think that a light axino is more natural in the context of gauge mediation. However it turns out that this is difficult to accommodate. The problem is with the saxion. Consider a model of minimal gauge mediation (see e.g. Ref. [103]), where the messengers do not carry any PQ charge and communicate with the visible sector only via gauge interactions. The leading contribution to the mass of the sfermions of the minimal supersymmetric Standard Model (MSSM) is generated at two loops. However the saxion is a gauge singlet and its mass can only be generated at three loops, as shown schematically in Fig. 9. The coupling c 1 that appears twice in the diagram is the one we discussed in the previous subsection, and the same that appears in the superpotential of Eq. (2.1). The saxion squared mass, m 2 s is suppressed by a factor of ∼ This light saxion could pose serious cosmological problems [104] as it would come to dominate the energy density of the universe for a long time before it decays. We reviewed some of the issues associated with saxion cosmology in Sec. 4. One way out of this problem would be to make the saxion heavier. This could be achieved by either coupling the axion superfield to the messengers or to fields in the hidden sector responsible for SUSY breaking. Apart from the extra field content that the procedure would bring into the model there is another issue that seems difficult to overcome: the same couplings needed to make the saxion heavier would very likely produce a heavy axino [105]. Perhaps there is a clever way to arrange for a heavy saxion and a light axino in gauge mediation, but in light of our considerations it seems that such a model would have to be quite complicated. We do not investigate this aspect further in this work. Rather we choose a supergravity (SUGRA) model. We showed in Ref. [46] that in SUGRA the axino would typically have a mass comparable to the gravitino, thus in the TeV range. However, we can adjust a single parameter, λ χ , to lower the axino mass down to the keV range. Doing so does not affect other masses, such as those of saxion and gravitino for instance. Also, as we argued above, since the axino we consider is of the DFSZ type, once we set its mass at tree level it will not be affected appreciably by loop corrections. Therefore the SUGRA model can be kept minimal, as opposed to a possible model with gauge mediation, at the expense of some tuning, needed to lower the mass of the axino. In the model presented in this paper the axino is much lighter than the gravitino. This represents an explicit exception to the generic argument that the axino should be heavier than the gravitino, put forward in Ref. [106].

B Axino-Gaugino Mixing
Most of the amplitudes involved in the radiative decay of the axino bear a strong dependence on its mixing with some of the other neutralino mass eigenstates. In this work we rely on well-motivated estimates for these quantities, but in order to obtain their precise values one should diagonalize the neutralino mass matrix M χ 0 . While an exact analytical algorithm exists for the case of the MSSM [107], adding one more mass eigenstate, as for the case of the Next-to-Minimal Supersymmetric Standard Model (NMSSM), does not allow for a closed-form solution. In this section we show how the approximate blockdiagonalization procedure of Ref. [108] for the NMSSM 5x5 neutralino mass matrix can be used to obtain a numerical value for the axino-gaugino mixing. The only requirement is a small mixing between the the singlet and the doublet states, which in our case corresponds to the mixing between any of the two Higgs doublets and the axion field A. From M χ 0 in Eq. (2.14) we see that such mixing is indeed small, since Recalling that the axino mass eigenstate consists of a linear combination of the fermionic components of the fields A andĀ only, we define the following neutralino mass 5x5 sub-matrix: where the axino mass m a is the lightest eigenvalue which results from the diagonalization of the 3x3 lower-right block of the neutralino mass matrix M χ 0 in Eq. (2.14). The diagonalization of M 5 proceeds in two steps: first, the 4x4 matrix V rotates only the MSSM upper-left 4x4 block into a diagonal form; next a second matrix block-diagonalizes the full resulting 5x5 matrix, in a way that the 4x4 MSSM sub-matrix M 4 is still diagonal up to corrections of the second order in the singlet-doublet mixing parameter µ eff fa . The 5x5 unitary matrix U which combines these two steps is The axino-gaugino mixings can then be read off from the off-diagonal entry 4 j=1 V ij Λ j of U , with j = B 0 , W 0 , H 0 u , H 0 d . In Fig. 10 we show two examples of how a more careful analysis can change significantly such quantities, and thus ultimately the phenomenology of the model. This possibly constitutes caveats to our above reasoning. In the left panel we observe that a value of tan β close to 1 strongly suppresses the value of x a, B , such that our previous estimates are off by several orders of magnitude. In the right panel instead we see how allowing for larger values of µ eff might lower x a, H , if M 1 < µ.