Coupled Boltzmann computation of mixed axion neutralino dark matter in the SUSY DFSZ axion model

The supersymmetrized DFSZ axion model is highly motivated not only because it offers solutions to both the gauge hierarchy and strong CP problems, but also because it provides a solution to the SUSY mu problem which naturally allows for a Little Hierarchy. We compute the expected mixed axion-neutralino dark matter abundance for the SUSY DFSZ axion model in two benchmark cases-- a natural SUSY model with a standard neutralino underabundance (SUA) and an mSUGRA/CMSSM model with a standard overabundance (SOA). Our computation implements coupled Boltzmann equations which track the radiation density along with neutralino, axion (produced thermally (TH) and via coherent oscillations (CO)), saxion (TH- and CO-produced), axino and gravitino densities. In the SUSY DFSZ model, axions, axinos and saxions go through the process of freeze-in-- in contrast to freeze-out or out-of-equilibrium production as in the SUSY KSVZ model-- resulting in thermal yields which are largely independent of the re-heat temperature. We find the SUA case with suppressed saxion-axion couplings (\xi=0) only admits solutions for PQ breaking scale f_a~<5\times 10^{12} GeV where the bulk of parameter space tends to be axion-dominated. For SUA with allowed saxion-axion couplings (\xi =1), then f_a values up to ~ 2\times 10^{14} GeV are allowed. For the SOA case, almost all of SUSY DFSZ parameter space is disallowed by a combination of overproduction of dark matter, overproduction of dark radiation or violation of BBN constraints. An exception occurs at very large f_a~ 10^{15}-10^{16} GeV where large entropy dilution from CO-produced saxions leads to allowed models.


Introduction
The recent discovery of the Higgs boson [1,2] with mass m h = 125.5 ± 0.5 GeV confirms the particle content of the Standard Model (SM) but carries with it a puzzle: why is the Higgs mass so light? Radiative corrections to the Higgs mass are of the form where c i is a loop dependent factor with |c i | ∼ 1 and Λ is the cutoff scale below which the SM ought to be valid. Setting δm 2 h = m 2 h and using e.g. c i = 1 tells us that Λ 1 TeV, i.e. that we expect new physics starting near the TeV scale. Yet so far, LHC data are in strong agreement with the SM.
The introduction of supersymmetry (SUSY) into the theory tames the quadratic divergences, and furthermore relates the Higgs mass to the Z mass, predicting m h 135 GeV within the Minimal Supersymmetric Standard model, or MSSM [3]. Comparing the measured value of m h to the theory prediction, one finds that the Higgs mass falls squarely within the narrow window predicted by SUSY.
A further problem with the SM arises in the QCD sector, where naively the U (2) L × U (2) R chiral symmetry of the light quark sector implies the existence of four-and not three-light pions. 't Hooft resolved this problem [4] via discovery of the QCD θ vacuum where the anticipated U (1) A symmetry is not respected [5]. A consequence of 't Hooft's solution is that the QCD Lagrangian contains a CP -violating term whereθ ≡ θ + arg det(M ), with M being the quark mass matrix. Measurements of the neutron EDM implyθ 10 −10 , so the term is somehow minuscule. An elegant resolution of this "strong CP " problem involves the introduction of a spontaneously broken global Peccei-Quinn (PQ) symmetry [6] and its concommitant axion field a [7]. For realistic models [8,9], the scale of PQ symmetry breaking f a is required to be f a 10 9 GeV lest red giant stars cool too quickly [10]. By enlarging the SM to include the PQ axion, Eq. 1.1 implies a Higgs mass m h ∼ f a . To solve the strong CP problem while simultaneously taming the Higgs mass, it seems both SUSY and PQ are required. In this case, the axion comprises but one element of an axion superfield given by where now the θ are spinorial Grassmann co-ordinates and F a is the axion auxiliary field. Here, s is the R-parity even spin-0 saxion field andã is the R-parity-odd spin-1 2 axino field. In gravity-mediated SUSY breaking models, then m s is a soft SUSY breaking term which is expected to be ∼ m 3/2 and the (more model-dependent) axino mass mã is also expected to be of order m 3/2 [11]. Here, the gravitino mass m 3/2 generated via the super-Higgs mechanism is expected to be of order the weak scale ∼ 1 TeV while the visible sector sparticle masses are also expected to be of order m 3/2 [12]. Lack of a SUSY signal at LHC8, along with a decoupling solution [13] to the SUSY flavor, CP , proton decay and gravitino problems all suggest m 3/2 to be more like ∼ 10 − 20 TeV. Meanwhile, SUSY electroweak naturalness requires the superpotential µ-term to be ∼ 100 − 200 GeV [14]. In such a case, one would expect the lightest neutralino to be the stable LSP and to be a higgsino-like WIMP dark matter candidate. However, in this case dark matter would be composed of an axion-neutralino admixture, i.e. two dark matter particles! A further problem with SUSY models is the so-called µ-problem. The superpotential Higgs/higgsino mass term µ is supersymmetric so that one expects it naively to have values of order the GUT or reduced Planck scales. But since it gives mass to the newly discovered Higgs boson (along with W ± and Z 0 ), phenomenology dictates it to be of order the weak scale. An elegant solution occurs within the context of the SUSY DFSZ axion model [9]. In this case, the SM Higgs doublets carry PQ charge so that the µ term is in fact forbidden. But there may exist non-renormalizable couplings of the Higgs doublets to a PQ-charged superfield S: where n is an integer ≥ 1. In this Kim-Nilles solution to the SUSY µ problem [15], under PQ symmetry breaking S receives a vev S ∼ f a so that an effective µ term is generated with This mechanism allows for µ m 3/2 since the µ-term arises from PQ symmetry breaking whilst m 3/2 might arise from hidden sector SUSY breaking. 1 For n = 1 and λ ∼ 1, µ ∼ 100 GeV requires f a ∼ 10 10 GeV while n > 1 allows for much larger values of f a . Alternatively, the Giudice-Masiero solution [16] to the µ-problem favors µ ∼ m 3/2 wherein tension then arises between SUSY naturalness and LHC sparticle mass bounds.
In the supersymmetric DFSZ model, the axion domain wall number is N DW = 6 since the quark doublet superfields carry PQ charge. As a result, the PQ symmetry must be broken before or during inflation 2 in order to avoid the overclosure of the universe through the production of stable domain walls [17]. In this case the axion misalignment angle (θ i ) is constant in our patch of the universe and the relic density from coherent oscillations of the axion is given by (ignoring possible entropy dilution effects) [18,19,20]: where f (θ i ) = ln e/(1 − θ 2 i /π 2 7/6 . It has been pointed out recently [21,22,23] that the measurement of a large tensorto-scalar ratio (r 0.2) by the BICEP2 collaboration [24] provides strong constraints on axion models. In particular, the breaking of the PQ symmetry before inflation (as required in the DFSZ model) would lead to too large isocurvature perturbations thus excluding this possibility. However, simple extensions of the PQ breaking sector are possible that can significantly affect the inflationary cosmology. One possible extension is to introduce an inflaton-dependent interaction that explicitly breaks the PQ symmetry. In this case the axion becomes massive during inflation and isocurvature perturbations do not develop [21]. Another possibility is to consider the case where the PQ breaking scale during inflation is larger than in the current universe, so isocurvature perturbations are suppressed. This scenario can be realized through the D-term interaction of the anomalous U (1) gauge symmetry in the PQ sector [25] or from Planck-suppressed interactions between the axion superfield and the inflaton superfield in the Kähler potential [26]. In the following discussion, we may assume that the isocurvature perturbation is suppressed by one of these extended PQ breaking scenarios, so the SUSY DFSZ model can be made compatible with the BICEP2 measurement. Alternatively, it remains to be seen whether the BICEP2 result is verified by further measurements at different frequency values [27,28].
Besides being produced through coherent oscillations, axions are also produced through thermal scatterings in the early universe. In this case, however, they are relativistic and constitute dark radiation, contributing to the number of relativistic degrees of freedom in the early universe. The amount of dark radiation produced during big bang nucleosynthesis (BBN) or during matter-radiation decoupling is usually parametrized by the number of effective neutrinos, which is conservatively constrained by BBN and CMB data to be N eff < 4.6 (or ∆N eff < 1.6). 3 However, as discussed in Sec. 2, the thermal production of axions is suppressed at temperatures below the Higgs/higgsino masses, resulting in a negligible contribution to ∆N eff . Nonetheless, relativistic axions may also be produced from saxion decays. The s → aa branching ratio is controlled by the axion-saxion effective coupling [11]: where ξ is a model dependent parameter, which can be small (or even zero) or as large as 1. Since the saxion decays strongly depend on ξ, we discuss the ξ = 0 and ξ = 1 limiting cases separately in Sec. 3. As already mentioned above, the total DM abundance in the SUSY DFSZ scenario receives contribution from both CO axions and relic neutralinos. The relic abundance of neutralinos in the SUSY DFSZ model was first considered in Ref. [30,31,32,33]. Neutralinos are produced through the usual freeze-out mechanism as well as through injection from saxion and axino decays. Therefore, in order to compute the final neutralino relic abundance it is necessary to determine the axino and saxion production rates and decay widths in the early universe. The axion multiplet couples to the MSSM primarily through its coupling with the Higgs supermultiplets, generated after breaking of the PQ symmetry as [33]: where 1 + Bθ 2 is a SUSY breaking spurion field and c H is the PQ charge of the Higgs bilinear operator H u H d . Since Eq. 1.8 generates tree level interactions of the type AH u H d , the thermal production of saxions, axions and axinos happen through the freeze-in mechanism [34]. In this case the production is maximal at T ∼ m s,ã , leading to thermal yields which are largely independent of the re-heat temperature (T R ) [31,35]. As discussed in Sec. 2, in some regions of parameter space the thermal production and decay of axinos and saxions are competing processes and cannot be treated separately. As a result, the sudden decay approximation is no longer valid and a precise calculation of the neutralino relic abundance (which receives contributions from axino and saxion decays) requires the numerical integration of the Boltzmann equations.
In the present work, we continue to refine the calculation of mixed axion-neutralino CDM in the SUSY DFSZ model. Here we compute the evolution of the axion, axino, saxion, neutralino and gravitino relic abundances using the appropriate system of coupled Boltzmann equations. In Ref's. [36,37], a similar calculation was performed in the SUSY KSVZ scenario that allowed for a more precise computation of the dark matter relic abundance; this method included the effects of the temperature-dependence of the neutralino annihilation cross section ( σv (T )) and the non-thermal production of neutralinos in models with large entropy injection from saxion decays. Here we apply a similar formalism to the SUSY DFSZ model, using the axino/saxion thermal production rates and decay rates computed in previous works [31,33,35]. This approach allows for • correct calculation of axino and saxion thermal yields for small f a values, • inclusion of temperature-dependent σv (T ) such as occurs for bino-like CDM with mainly p-wave annihilation, • inclusion of non-sudden axino/saxion decays and • accurate calculation of entropy production and injection in the early universe.
Furthermore, we also scan the parameter space of the SUSY DFSZ model and identify the regions consistent with dark matter, BBN and dark radiation constraints. The remainder of this paper is organized as follows. In Sec. 2, we discuss the set of coupled Boltzmann equations used to compute our numerical results. In Sec. 3, we present the two benchmark models used in our analysis and discuss the behavior of the dark matter relic abundance in these models as a function of the PQ parameters. In order to keep our results general, we scan over the most relevant PQ parameters and numerically solve the Boltzmann equations for each point. We also discuss the BBN and ∆N eff constraints in these models. Finally, in Sec. 4, we present a brief summary and conclusions.

Coupled Boltzmann equations
Our goal is to numerically solve the coupled Boltzmann equations which track the number and energy densities of neutralinos Z 1 , gravitinos G, saxions s, axinosã, axions a and radiation as a function of time starting at the re-heat temperature T = T R at the end of inflation until today. For axions and saxions, we include separately the thermally produced (TP) and coherent oscillating (CO) components. The simplified set of Boltzmann equations for the SUSY KSVZ model as well as the method for their numerical solution were presented in detail in Ref's. [36,37]. In this section we discuss the main differences between the KSVZ and DFSZ scenarios and how the simplified Boltzmann equations derived for the KSVZ case must be generalized in order to allow for a proper computation of the relic abundances in the DFSZ model.
In the KSVZ model considered in Ref. [36], the thermal production of saxions, axions and axinos is maximal at T ∼ T R (for re-heat temperatures below the decoupling temperature of saxions and axinos), resulting in a thermal yield proportional to the reheat temperature [38]. Also, since the axino/saxion decay widths are suppressed by the loop factor as well as by the PQ scale, their decays tend to take place at temperatures T m, where m is the axino or saxion mass. Hence the thermal production and decay processes can be safely treated as taking place at distinct time scales. Furthermore, the inverse decay process (a + b →ã, s) is always Boltzmann-suppressed when the decay term becomes sizable (Γ ∼ H), thus we can easily neglect the inverse decay contributions. In the DFSZ scenario, however, the situation can be drastically different. Here, the tree-level couplings between the axion supermultiplet and the Higgs superfields (Eq. 1.8) modify the thermal scatterings of saxions, axions and axinos and can significantly enhance their decay widths. As computed in Ref. [35], the thermally averaged cross-section for axino (or saxion) production is proportional to: where M is the characteristic scale of the scattering (either the higgsino or saxion/axino mass) and K 2 is the modified Bessel function. As we can see from the above expression (unlike the KSVZ case) the production is maximal at T M/3 T R . Hence most of the thermal production of axinos and saxions takes place at T ∼ M , resulting in thermal yields which are independent of T R . This behavior is similar to the freeze-in mechanism [34], where a weakly interacting (and decoupled) dark matter particle becomes increasingly coupled to the thermal bath as the universe cools down. However, in the current scenario, the "frozenin" species (axinos and saxions) are not stable and their decays will only contribute to the dark matter (neutralino or axion) relic abundance if they take place after neutralino freezeout.
The coupling in Eq. 1.8 can also enhance the axino/saxion decay width for large µ values, since the coupling to Higgs/higgsinos is proportional to µ/f a . As a result, saxions and axinos may decay at much earlier times (larger temperatures) when compared to the KSVZ scenario. If their decay temperatures are of order of their masses, then inverse decay processes such as Z 1 + h →ã or h + h → s can no longer be neglected. In fact, in Ref. [33] it was shown that the decay temperatures can indeed be larger than the axino or saxion mass, so the inverse decay process can be significant. The main effect of including the inverse decay process is to delay the axino/saxion decay. This is an important effect which cannot be accounted for in the sudden decay approximation and requires the numerical solution of the Boltzmann equations. We point out, however, that the inverse decay process is only relevant for T decay M , since if the decay happens at lower temperatures, the inverse decay process is Boltzmann-suppressed. As a result, the inverse decay process will only be relevant for the cases where axinos and saxions decay before neutralino freeze-out (since T f r ∼ m Z 1 /20 M ) and we do not expect it to modify the neutralino relic abundance. Nevertheless, it is essential to include inverse decays in the Boltzmann equations for consistency.
As discussed above, inverse decay processes were not relevant for the KSVZ case and were neglected in Ref's [36,37]. With the addition of the inverse decay process, the Boltzmann equations for the number (n i ) and energy (ρ i ) densities of a thermal species i (= a, s orã) reads: 4 ib andn i is the equilibrium density of particle species i. 5 It is also convenient to use the above results to obtain a simpler equation for ρ i /n i : where P i is the pressure density (P i 0 (ρ i /3) for non-relativistic (relativistic) particles).
As discussed in Ref. [36], we track separately the CO-produced and the TP components of the axion and saxion fields since we assume the CO components do not have scattering contributions. Under this approximation, the equations for the CO-produced fields (axions and saxions) read: The generalization of the Boltzmann equations to include decays to n-body final states (n > 2) is straightforward. For the cases where 3-body decays are relevant (such as gravitino decays), we use the appropriate generalized equations. 5 In general the leading thermal scatterings for saxions, axions and axinos involve single production of these, resulting in a contribution to the Boltzmann equations linear in their densities [35]: (ninj −ninj) σv where i = {s, a,ã} and j = MSSM. We simply assume nj = ni andnj =ni for scattering contributions for saxion, axion and axino to avoid the complication that comes from dealing with thousands of processes. This assumption does not make a significant difference in the final result since the threshold scale of the process does not change substantially.
The amplitude of the coherent oscillations is defined by the initial field values, which for the case of PQ breaking before the end of inflation is a free parameter for both the axion and saxion fields. We parametrize the initial field values by θ i = a 0 /f a and θ s = s 0 /f a . Finally, we must supplement the above set of simplified Boltzmann equations with an equation for the entropy of the thermal bath: where R is the scale factor and BR(i, X) is the fraction of energy injected in the thermal bath from i decays. In order to solve the above equations, it is necessary to compute the values of the decay widths and σv appearing in Eqs. 2.2, 2.4 and 2.6. Since these have been presented in previous works, we just refer the reader to the relevant references. The value of σv for thermal axino production is given in Ref's [31,35], while σv for neutralino annihilation is extracted from IsaReD [39]. For thermal saxion and axion production, it is reasonable to expect annihilation/production rates similar to axino's, since supersymmetry assures the same dimensionless couplings. Hence we apply the result for axino thermal production from Ref's [31,35] to saxions and axions. For the gravitino thermal production we use the result in Ref. [40]. The necessary saxion and axino partial widths and branching fractions can be found in Ref. [33], while the gravitino widths are computed in Ref. [41].
In order to illustrate the effects discussed above, in Fig. 1 we show a specific solution of the Boltzmann equations, where we take a MSSM model with µ = 2.5 TeV (the SOA benchmark defined in Sec. 3) and T R = 10 7 GeV, f a = 10 10 GeV, m G = 10 TeV, mã = 1 TeV and m s = 500 GeV. We also take the saxion and axion mis-alignment angles (θ s and θ i ) equal to 1 and ξ = 1 so s → aa and s →ãã decays are turned on. The figure shows the evolution of the yields (n i /s) of both the thermally-produced and coherent oscillating fields versus the inverse of the temperature. First we point out that, as expected in the DFSZ case, the thermal saxion and axino yields increase as the temperature is reduced, reaching their maximal value just before their decay. This example clearly shows how both processes (thermal production and decay) happen simultaneously, as previously discussed. The axion thermal production follows a similar behavior, but since the axion is (effectively) stable, its yield remains constant after the thermal production becomes suppressed at T µ. Gravitinos are also produced through thermal scatterings. However, as seen in Fig. 1, their production cross-section peaks at T ∼ T R , much like the saxion/axino production in the KSVZ case. The small increase in the yields around T = 1 TeV is due to the reduction in the number of relativistic SUSY degrees of freedom in the thermal bath, which reduces the entropy density. We also show as dashed lines the respective yields without the inclusion of the inverse decay process. As seen in Fig. 1, the inclusion of inverse decays delays the decay of saxions and axinos, with the effect being larger for saxions, since they tend to decay earlier. Nonetheless, the neutralino and axion relic densities are unchanged, as expected from the discussion above. For the current point chosen, the final neutralino relic density equals its MSSM value and is well above the experimental limits: Ω Z 1 h 2 = 6.8. In the following sections, we will apply the Boltzmann equations presented here to numerically compute the neutralino and axion relic abundances. We will also compute the thermal axion abundance (including the contributions from saxion decays) in order to evaluate its contribution to the number of effective neutrinos (dark radiation), as discussed in Sec. 1.

Benchmark points
In order to compute the dark matter relic abundance in the SUSY DFSZ model we must specify both the PQ and the MSSM parameters. Since the axion supermultiplet interactions are proportional to µ, we consider in our numerical analysis two benchmark MSSM points: one with a small and one with a large value of µ. In the first case (SUA), the neutralino LSP is mostly a higgsino, resulting in a standard underabundance of neutralino cold dark matter (CDM). The second benchmark, which we label SOA, has a standard thermal overabundance of neutralino dark matter, since the neutralino is mostly a bino.
For the SOA case, we adopt the mSUGRA/CMSSM model with parameters (m 0 , m 1/2 , A 0 , tan β, sign(µ)) = (3500 GeV, 500 GeV, −7000 GeV, 10, +). (3. 2) The SOA point has mg = 1.3 TeV and mq 3.6 TeV, so it is just beyond current LHC8 sparticle search constraints. It is also consistent with the LHC Higgs discovery since m h = 125 GeV. The lightest neutralino is mainly bino-like with m Z 1 = 224.1 GeV, and the standard neutralino thermal abundance is found to be Ω MSSM Z 1 h 2 = 6.8, a factor of ∼ 57 above the measured value. Due to its large µ parameter, this point has very high electroweak finetuning [46].
In Fig. 2, we show the solution of the Boltzmann equations for the SUA point with T R = 10 7 GeV, f a = 10 11 GeV, m G = 10 TeV, mã = m s = 5 TeV, θ s = 1, ξ = 1 and θ i = 3.113. We present the evolution of the energy densities of axions and saxions (both CO-and thermally produced), axinos, neutralinos and gravitinos as a function of the scale factor of the universe R/R 0 , where R 0 is the scale factor at T = T R . For this parameter set, the final neutralino abundance is Ω Z 1 h 2 = 0.006 whilst the axion abundance is Ω a h 2 = 0.114, resulting in a total dark matter relic abundance within the measured value. 6 We see that at T = T R (where R/R 0 ≡ 1) the universe is radiation-dominated with smaller abundances of neutralinos, TP axions, axinos and saxions, and even smaller abundances of CO-produced saxions and TP gravitinos. The CO-produced saxions evolve as a non-relativistic matter fluid and so their density diverges from the relativistic gravitino abundance as R increases. Both TP-and CO-populations of saxions begin to decay around R/R 0 ∼ 10 5 , at temperatures (T ∼ 10 2 GeV) well below their masses. Furthermore, since this decay happens before neutralino freeze-out, the TP-neutralino population is unaffected. Somewhat later, but still before neutralino freeze-out, the axino population decays. As a result, the neutralino relic abundance is hardly increased by the injection of neutralinos from axino decays. The axion mass turns on around T ∼ 1 GeV so that the axion field begins to oscillate around R/R 0 2 × 10 7 . The CO-produced axion field evolves as CDM and ultimately dominates the universe at a value of R/R 0 somewhat off the plot. The behavior of the DFSZ axinos and saxions-in that they tend to decay before neutralino freeze-out-is typical of this model for the lower range of f a 10 12 GeV with TeV-scale values of mã and m s [33,32].
Finally, gravitinos are long-lived and decay well after the neutralino freeze-out, at T ∼ O(100) keV. However, for T R = 10 7 GeV, gravitinos typically have a small number density and contribute marginally to the final neutralino relic abundance. Also-due to their small energy density-the gravitino decays do not have any significant impact on big bang nucleosynthesis.
In the following subsections, we compute the neutralino and axion relic abundances for below the IsaReD output due to our fit of the IsaReD σv (T ) function.
the two benchmark points through the numerical integration of the Boltzmann equations presented in Sec. 2. In order to be as general as possible, we will scan over the following SUSY DFSZ parameters: For simplicity, we will fix the initial saxion field strength at s i = f a (θ s ≡ s i /f a = 1) with m G = 10 TeV. Unlike the SUSY KSVZ model, the bulk of our results do not strongly depend on the re-heat temperature (T R ) since the axion, axino and saxion TP rates are independent of this quantity. Nonetheless, the gravitino thermal abundance is proportional to T R and since gravitinos are long-lived they may affect BBN if T R is sufficiently large. In order to avoid the BBN constraints on gravitinos, we choose T R = 10 7 GeV, which results in a sufficiently small (would-be) gravitino abundance. As a result, gravitinos typically do not contribute significantly to the neutralino abundance, as discussed above. For each of the SUA and SOA benchmark points, we consider two different cases: ξ = 0 and ξ = 1. As we can conclude from Eq. 1.7, saxion decays into axions and axinos are turned off if ξ = 0 whereas s → aa and s →ãã decays are allowed for ξ = 1.

Mixed axion/higgsino dark matter: SUA with ξ = 0
In this section, we will examine the SUA SUSY benchmark assuming no direct coupling between saxions and axions/axinos (see Eq. 1.7), which corresponds to ξ = 0. For each parameter set which yields an allowable value of Ω Z 1 h 2 < 0.12, we will adjust the initial axion misalignment angle θ i such that Ω Z 1 h 2 + Ω a h 2 = 0.12, i.e. the summed CDM abundance saturates the measured value by adjusting the initial axion field strength parameter θ i .
Our first results are shown in Fig. 3a where we plot Ω Z 1 h 2 vs. f a for a scan over the parameter space defined in Eq. 3.4. Since for large f a values, saxions and axinos may decay during BBN, we apply the BBN constraints using the bounds from Jedamzik [47] with extrapolations for intermediate values of m X other than those shown in his plots. These constraints depend on the lifetime of the decaying state, its energy density before decaying and the fraction of energy injected as hadrons or color-charged states (R h ). In the DFSZ scenario the dominant decays of saxions are into neutralinos, charginos, Higgs states or gauge bosons. Also, axinos decay into neutralinos or charginos plus gauge bosons or Higgs states. Thus the branching ratio for s → hadrons must be similar to Br(W/Z → quarks) or Br(Higgs → quarks), resulting in R h ∼ 1. So we conservatively take R h = 1 for saxion and axino decays. In Fig. 3a the red points violate BBN bounds on late-decaying neutral relics, while the blue points are BBN safe. The points below the solid gray line at 0.12 are DM-allowed, whilst those above the line overproduce neutralinos and so would be ruled out. The dashed gray line denotes the level of equal axion-neutralino DM densities: each at 50% of the measured abundance. Since, as previously discussed, the thermal production of axions gives a negligible contribution to ∆N eff and, for ξ = 0, there is no axion injection from saxion decays, dark radiation constraints are always satisfied in this case.
For low values of f a ∼ 10 9 − 10 10 GeV, we see that Ω Z 1 h 2 takes on its standard thermal value listed in Table 1. This is because with such a small value of f a , the axino and saxion couplings to matter are sufficiently strong that they always decay before neutralino freeze-out. This behavior was also shown in Ref's [32,33] using semi-analytic calculations. In this region, we expect mainly axion CDM with ∼ 5 − 10% contribution of higgsino-like WIMPs [32]. As f a increases, then saxions and axinos decay more slowly, and often after neutralino freeze-out. The late decays of saxions and axinos increases the neutralino density. If the injection of neutralinos from saxion/axino decays is sufficiently large, the 'supersaturated' decay-produced neutralinos re-annihilate, reducing their density. Although re-annihilation can reduce the neutralino density by orders of magnitude, its final value is always larger than the freeze-out density in the standard MSSM cosmology [48].
As f a increases, the thermal production of axinos and saxions decreases, while the density of CO-produced saxions increases (since we take θ s = s 0 /f a = 1). For f a 10 12 GeV, axinos and saxions are mostly thermally produced and Ω Z 1 h 2 rises steadily with f a mainly due to the increase of axino and saxion lifetimes, resulting in a late injection of neutralinos well after their freeze-out. On the other hand, for f a 5 × 10 12 GeV, the thermal production of axions and axinos becomes suppressed and the main contribution to the neutralino abundance comes from CO-produced saxions and their decay. As seen in Fig. 3, once axinos and saxions start to decay after the neutralino freeze-out (f a 5 × 10 10 GeV), Ω Z 1 h 2 always increases with f a : this is due to the increase in saxion and axino lifetimes and also due to the increase in rate of CO-produced saxions. By the time f a exceeds 10 13 GeV, then always too much neutralino CDM is produced and the models are excluded. BBN constraints do not kick in until f a exceeds ∼ 10 14 GeV. For a given f a value, the minimum value of Ω Z 1 h 2 seen in Fig. 3 happens for the largest saxion/axino masses considered in our scan (20 TeV). This is simply due to the fact that the lifetime decreases with the saxion/axino mass, resulting in earlier decays. As a result, neutralinos are injected earlier on and can re-annihilate more efficiently, since their annihilation rate increases with temperature. Hence, an increase in the axino/saxion mass usually implies a decrease in the neutralino relic abundance (for a fixed f a value).
In Fig. 3b, we show the value of the axion misalignment angle θ i which is needed to obtain Ω Z 1 h 2 + Ω a h 2 = 0.12. For low f a values (∼ 10 9 − 10 10 GeV), rather large values of θ i ∼ π are required to bolster the axion abundance into the range of the measured CDM density. For values of f a ∼ 10 11 − 10 12 GeV, then (perhaps more natural) values of θ i ∼ 1 are required. For f a 10 12 GeV, axions tend to get overproduced by CO-production and so a small value of θ i 0.5 is required for suppression. For even higher f a values, too many neutralinos are produced, so the models are all excluded.

Mixed axion/higgsino dark matter: SUA with ξ = 1
We now discuss the main changes in the results of Fig. 3 if we consider a non-vanishing saxion-axion/axino coupling. For simplicity, we take ξ = 1 where ξ is defined in Eq. 1.7. In this case saxions can directly decay to axions and axinos (if m s > 2mã). The s → aa decay usually dominates over the other decays [33], suppressing BR(s → . . . → Z 1 Z 1 ) and significantly reducing the neutralino injection from saxion decays. As a result, the Figure 3: In a) we plot the neutralino relic density from a scan over SUSY DFSZ parameter space for the SUA benchmark case with ξ = 0. The grey dashed line shows the points where DM consists of 50% axions and 50% neutralinos. In b), we plot the misalignment angle θ i needed to saturate the dark matter relic density Ω Z1a h 2 = 0.12. neutralino relic abundance is usually smaller (for the same choice of PQ parameters) than the ξ = 0 case. Furthermore, the saxion lifetime is reduced (due to the large s → aa width) and saxions tend to decay earlier when compared to the ξ = 1 case.
In Fig. 4a, we once again show Ω Z 1 h 2 vs. f a for the SUA SUSY benchmark but now for ξ = 1. As just discussed, in this case the saxion lifetime is reduced, so the region of f a where saxions/axinos always decay before freeze-out is extended beyond the values generated for the ξ = 0 case. Since BR(s → . . . → Z 1 Z 1 ) is suppressed in the ξ = 1 case, saxions do not significantly contribute to Ω Z 1 h 2 except when f a 10 14 GeV where CO-produced saxions have such large densities that-even though their branching ratio to neutralinos is at the 0.1% level-their decay still enhances the neutralino relic density. For 10 11 GeV f a 10 14 GeV however, Ω Z 1 h 2 is dominated by the thermal axino contribution and the neutralino relic density increases with f a , as in the ξ = 0 case. Once f a 10 13 GeV, the thermal production of axinos becomes strongly suppressed and despite decaying well after neutralino freeze-out, their contribution to Ω Z 1 h 2 starts to decrease as f a increases. This is seen by the turn over of Ω Z 1 h 2 around f a ∼ 10 13 GeV. As f a increases past 10 14 GeV, CO saxions start to contribute to the neutralino relic density, which once again rises with f a .
Another important difference in the ξ = 1 case is the large injection of relativistic axions from saxion decays. For large values of f a , where the density of CO saxions is enhanced, the injected axions have a non-negligible contribution to ∆N eff . In particular, for f a 10 14 GeV, CO saxion decays produce too much dark radiation, so this region (shown by brown points in Fig. 4a) is excluded by the CMB constraints on dark radiation (∆N eff < 1.6). 7 These points are also excluded by overproduction of neutralinos and violation of BBN bounds. We also show as green points the cases where ∆N eff ∼ 0.4 − 1.6 which could explain a possible excess of dark radiation suggested by the combined WMAP9 result. However these points are already excluded by overproduction of dark matter.
Finally, in Fig. 4b, we again plot the value of θ i which is needed by axions so that one matches the measured abundance of CDM, as described in the previous section. Once again, at very low f a , |θ i | ∼ π is required, while for very high f a values ( 10 13 GeV), low |θ i | is required in order to suppress axion CO-production. The range of f a ∼ 10 11 − 10 12 GeV tends to give θ i ∼ 1. Furthermore, since Ω Z 1 h 2 is usually smaller in the ξ = 1 case (when compared to ξ = 0), the CO axion contribution to DM can be larger and higher values of θ i are usually allowed, as seen in Fig. 4b. 3.4 Mixed axion/bino dark matter: SOA with ξ = 0 In this Section, we turn to the SUSY benchmark SOA, which features a bino-like LSP with a standard thermal overabundance Ω Z 1 h 2 = 6.8, i.e. too much dark matter by a factor 57! The SUSY µ parameter has a value of µ = 2598 GeV so this model would be considered fine-tuned in the electroweak sector. However, the large µ-parameter also bolsters the saxion and axino decay rates which are proportional to some power of µ (µ 2 or µ 4 ) in the SUSY DFSZ model [33].
In Fig. 5, we show the coupled Boltzmann calculation of Ω Z 1 h 2 as a function of f a for the SOA benchmark with ξ = 0. At low f a ∼ 10 9 − 10 10 GeV, axinos and saxions decay before neutralino freeze-out, so the model remains excluded due to overproduction of dark matter. As f a increases, neutralinos are only produced at higher and higher rates as their population is bolstered by late time axino and saxion decay, as already observed for the SUA case. In this region, the highest Ω Z 1 h 2 values for a fixed f a are obtained for the smallest m s , mã values, since these correspond to the longest lifetimes. However, once f a 10 14 GeV, a subset of points present the opposite behavior and the neutralino relic abundance actually decreases with f a . This region of parameter space corresponds to Figure 4: In a) we plot the neutralino relic density from a scan over SUSY DFSZ parameter space for the SUA benchmark case with ξ = 1. The grey dashed line shows the points where DM consists of 50% axions and 50% neutralinos. In b), we plot the misalignment angle θ i needed to saturate the dark matter relic density Ω Z1a h 2 = 0.12.
small saxion masses, m s 2m Z 1 , so the decay to neutralinos is kinematically forbidden. As a result (since ξ = 0, saxions do not decay to axions) the only effect of saxion decays is to inject entropy in the early universe. For f a 10 14 GeV, there is a huge rate for saxion production via coherent oscillations and the entropy injection from saxion decays can reduce the neutralino density, resulting in DM-allowed models with Ω Z 1 h 2 < 0.12. We discuss these cases in detail in Sec. 3.6. We also point out that in the SUA model or in the SUSY KSVZ case, such large f a values imply very long-lived saxions, with lifetimes of the order of O(10 s) or greater. As a result, all the solutions with large entropy injection in the SUA case are excluded by BBN contraints. 8 However, for the SOA case, the large µ value enhances the saxion decay rate to Higgs pairs and vector bosons and even at such high f a values, saxions can still decay before BBN starts. Very few points do succumb to BBN constraints (denoted by red points) but these are also excluded due to an overabundance of neutralinos. In Fig. 5 we also see that in the large f a region there is a visible gap (for a fixed f a value) between the branch with a suppression of Ω Z 1 h 2 and the one with an enhanced value of Ω Z 1 h 2 . The lower branch (with Ω Z 1 h 2 20) corresponds to points with low saxion masses, where BR(s → . . . Z 1 Z 1 ) 1, so saxion decays mostly dilute the neutralino relic density. Once m s > 2mt 1 , the s →t 1t1 channel becomes kinematically allowed and there is a sudden increase in Ω Z 1 h 2 , resulting in the gap seen in Fig. 5. Finally, since ξ = 0, axions are only thermally produced resulting in a negligible contribution to ∆N eff , so dark radiation constraints are inapplicable in this case.
3.5 Mixed axion/bino dark matter: SOA with ξ = 1 In Fig. 6 we plot Ω Z 1 h 2 vs. f a for the SOA SUSY benchmark but with ξ = 1. Unlike the SUA case, decays to axions are not always dominant, since Γ(s → aa) ∼ m 3 . Hence saxions dominantly decay to gauge bosons/higgses, except for m s µ. The low f a behavior of Ω Z 1 h 2 is much the same as in the ξ = 0 case: the neutralino abundance is only bolstered to even higher values and thus remains excluded by overproduction of WIMPs. As in the SOA ξ = 0 case, there again exists a set of points with f a 10 15 GeV and with m s 2m Z 1 which is allowed by all constraints. This is possible in the ξ = 1 case, since, for m s µ, saxions mainly decay to higgses and gauge bosons, thus injecting enough entropy to dilute Ω Z 1 h 2 . Points with m s µ, by the PQ breaking scale (θs = s0/fa = 1). As shown in Ref. [36], in the KSVZ case the neutralino relic abundance can be suppressed if one takes s0 fa or θs 1. however, have BR(s → aa) 1, resulting in a large injection of relativistic axions and a suppression of entropy injection. In this case many models start to become excluded by overproduction of dark radiation (brown points) while some also have ∆N eff ∼ 0.4 − 1.6: these points could explain a possible excess of dark radiation except that they also always overproduce neutralino dark matter. Thus, we see that the SUSY DFSZ model with large µ and either small or large ξ along with small m s is able to reconcile the expected value of Peccei-Quinn scale [49] from string theory [50,51] (where f a is expected ∼ m GUT ) with dark matter abundance, dark radiation and BBN constraints.

Mixed axion/bino dark matter with a light saxion
As discussed in the previous sections, the neutralino relic abundance can only be suppressed with respect to its MSSM value if m s 2m Z 1 and f a 10 14 GeV. Here we discuss in detail this case, since it represents the only possibility to reconcile the SOA dark matter scenario with the measured DM abundance. A specific example is shown in Fig. 7, where the evolution of the energy density of various species as a function of the universe scale factor is presented for f a = 4.88 × 10 15 GeV, m s = 456 GeV and mã = 3.3 TeV. For this choice of parameters, the neutralino relic abundance is highly suppressed (Ω Z 1 h 2 = 0.051) but does comprise ∼ 50% of the total DM abundance. The remainding 50% is composed of axions although these require a somewhat small value of the axion mis-alignment angle (θ i = 0.03) in order to suppress the CO axion production. From Fig. 7 we see that the CO-produced saxion energy density dominates over the radiation energy density at R/R 0 ∼ 10 6 and decays at R/R 0 ∼ 10 10 , so that the universe is saxion-dominated during this period. In this case, saxions dominantly decay into SM particles, since the rate for saxion → neutralinos is highly suppressed by the kinematic phase factor (BR(s → Z 1 Z 1 ) ∼ 10 −8 at this point). Therefore, a huge amount of entropy is produced as we can see from the radiation curve (grey), while the neutralino density (blue) is almost unaffected by the saxion decay. As a result, the final neutralino density is given by Ω Z 1 = 0.051 and this can be a viable model, even though the PQ scale is very large. In Fig. 8, we once again scan over the parameter space from Eq. 3.4, but now we focus on the light saxion region, 400 GeV< m s < 500 GeV, where the saxion decay to neutralinos is kinematically suppressed or forbidden. As already seen in Figs. 5 and 6, in this case the large f a region (f a 10 15 GeV) can suppress Ω Z 1 h 2 to values below the observed DM abundance. Fig. 8b shows how the entropy dilution factor (r ≡ S f /S 0 ) increases with f a , reaching values as high as 10 4 , for f a 10 15 GeV.
We note here that one might wonder if the large f a ∼ m GUT region of the SUA model might be DM-allowed if we consider m s < 200 GeV so that saxion decay to SUSY particles is dis-allowed and saxion decay leads only to entropy dilution. Aside from the fact that such light values of m s leads to a large disparity between scalar soft breaking terms, in the SUSY DFSZ model these points should all be BBN dis-allowed as shown by the red points in Fig. 3.

Conclusion
In this paper, we have discussed the evaluation of the relevant Boltzmann equations in the supersymmetrized DFSZ axion model. This is a highly motivated scenario, since it Figure 8: In a), we plot the neutralino relic density vs f a for the scan over the SUSY DFSZ parameter space for the SOA benchmark case with ξ = 0 and 1, but with m s : 400 − 500 GeV. The grey dashed line shows the points where DM consists of 50% axions and 50% neutralinos. In b) we plot the dilution factor r vs. f a . provides solutions to the gauge hierarchy and strong CP problems as well as a solution to the SUSY µ problem while allowing for the Little Hierarchy µ m 3/2 which is expected from combining naturalness considerations with LHC bounds on sparticle masses and the 125 GeV Higgs boson mass. In SUSY DFSZ, axinos and saxions tend to decay to vector bosons, Higgs states and higgsinos. Saxions may also decay into aa orãã, depending on the value of the saxion-axion/axino model dependent coupling ξ. The first of these leads to dark radiation while the second may enhance the neutralino relic density.
In the SUSY DFSZ scenario, the decay widths of saxions and axinos are enhanced for large µ values and their decays may happen at temperatures of the order of their masses. Hence it is crucial to include the inverse decay processes in the Boltzmann equations. Furthermore, since in the SUSY DFSZ case the thermal production of saxions, axions and axinos happen through the freeze-in mechanism, the production and decay processes may happen at similar time scales. In these cases, a precise calculation of the saxion and axino evolution is only possible through the numerical integration of the Boltzmann equations.
Since most of the axion supermultiplet couplings in the SUSY DFSZ model are proportional to µ, we have presented results for two SUSY benchmark points: 1. a natural SUSY model labelled SUA with µ = 110 GeV and a higgsino LSP, and 2. a mSUGRA/CMSSM point (SOA) with µ = 2.6 TeV and a bino-like LSP, resulting in a standard thermal neutralino overabundance. We found that, for the SUA benchmark with ξ = 0, low f a ∼ 10 9 − 10 11 GeV tends to give mainly axion CDM with 5-10% higgsino-like WIMPs. For higher f a (∼ 10 11 − 10 12 GeV), the WIMP density increases and might even dominate the DM abundance. For f a > 10 13 GeV, the model becomes excluded due to overproduction of WIMPs. For SUA with ξ = 1 the contribution of s → aa hastens the saxion decay rate so that saxion decay occurs before neutralino freeze-out over an even larger range of f a . In this case, for sufficiently heavy saxions and axinos, f a ∼ 10 9 − 10 14 GeV is allowed by all constraints. For even higher f a values (f a > 10 14 GeV), the model becomes excluded by overproduction of WIMPs, overproduction of dark radiation and violation of BBN constraints.
For the SOA model, the presence of axions, saxions and axinos typically leads to an enhancement of the neutralino relic abundance for almost the entire f a range, so such models typically remain excluded. The exception comes at very large f a values (∼ 10 15 − 10 16 GeV) with small saxion masses, m s 2m Z 1 . In this case, enormous entropy injection from CO-produced saxions along with their decays to SM particles leads to entropy dilution of the WIMP relic density whilst avoiding BBN and dark radiation constraints.
For all allowed cases, we would ultimately expect both WIMP and axion dark matter detection to occur.