Gravitino or Axino Dark Matter with Reheat Temperature as high as $10^{16}$ GeV

A new scheme for lightest supersymmetric particle (LSP) dark matter is introduced and studied in theories of TeV supersymmetry with a QCD axion, $a$, and a high reheat temperature after inflation, $T_R$. A large overproduction of axinos ($\tilde{a}$) and gravitinos ($\tilde{G}$) from scattering at $T_R$, and from freeze-in at the TeV scale, is diluted by the late decay of a saxion condensate that arises from inflation. The two lightest superpartners are $\tilde{a}$, with mass of order the TeV scale, and $\tilde{G}$ with mass $m_{3/2}$ anywhere between the keV and TeV scales, depending on the mediation scale of supersymmetry breaking. Dark matter contains both warm and cold components: for $\tilde{G}$ LSP the warm component arises from $\tilde{a} \rightarrow \tilde{G}a$, while for $\tilde{a}$ LSP the warm component arises from $\tilde{G} \rightarrow \tilde{a}a$. The free-streaming scale for the warm component is predicted to be of order 1 Mpc (and independent of $m_{3/2}$ in the case of $\tilde{G}$ LSP). $T_R$ can be as high as $10^{16}$ GeV, for any value of $m_{3/2}$, solving the gravitino problem. The PQ symmetry breaking scale $V_{PQ}$ depends on $T_R$ and $m_{3/2}$ and can be anywhere in the range $(10^{10} - 10^{16})$ GeV. Detailed predictions are made for the lifetime of the neutralino LOSP decaying to $\tilde{a}+ h/Z$ and $\tilde{G}+h/Z/\gamma$, which is in the range of $(10^{-1}-10^6)$m over much of parameter space. For an axion misalignment angle of order unity, the axion contribution to dark matter is sub-dominant, except when $V_{PQ}$ approaches $10^{16}$ GeV.


JHEP03(2017)005
Contents 1 Introduction 2 2 Saxion cosmology 6 2.1 High reheat temperature after inflation: T R 10 10 GeV 6 2.2 Low reheat temperature after inflation: T R 10 10 GeV 9 2.3 Field equations 10 3 Axino and gravitino production 12 3.1 Freeze-in production of axinos 12 3.2 UV production of axinos 13 3.3 UV production of gravitinos 15 3.4 Freeze-in production of gravitinos 16 4 Axino and gravitino as the lightest superpartners 16  Perturbative theories with supersymmetry broken at the TeV scale are well-motivated by the hierarchy problem, even if they do not completely solve it, and lead to Higgs boson masses in the region discovered at the LHC. Dark matter could be the lightest superpartner (LSP), cosmologically produced by the freeze-out mechanism. On the other hand, the strong CP problem is elegantly solved by introducing a Peccei-Quinn (PQ) symmetry [1,2] broken at scale V PQ , 1 leading to a light axion degree of freedom a [3,4] that relaxes the CP -violating phaseθ to zero. In this case dark matter could be axions produced by the misalignment mechanism, with V PQ of order 10 12 GeV a motivated possibility. However, the cosmology of these two leading candidates for dark matter, LSPs and axions, is changed enormously in theories that have both weak scale supersymmetry and axions. The axion, a, must be promoted to a superfield and the saxion, s, and the axino,ã, both play central roles in cosmology.
In this work, for reasons discussed below, we focus on DFSZ theories [5,6], where the PQ symmetry forbids the µ term of the minimal supersymmetric standard model (MSSM). At the scale V PQ , PQ breaking induces the µ term as well as a coupling of the axion supermultiplet with the MSSM Higgs superfields with q µ a model dependent parameter defined in appendix A. The superpotential cubic coupling is responsible for axino production in the early universe [7,8], either through decays or inverse decays of charginos and neutralinos,χ →ã. This axino production by the freeze-in (FI) mechanism is IR dominated [9,10], namely most of the axinos are produced at temperatures around the TeV scale. Depending on the fermion content of the PQ breaking sector, a large abundance of axinos can also be produced in the UV, at the temperature T R at the end of inflationary reheating [11,12], analogous to UV production of gravitinos [13][14][15].
In order to make this distinction sharper, we define two different types of theories: DFSZ 0 : the heaviest colored fermion carrying PQ charge is the top quark, so the only source for axino production is the IR dominated freeze-in; DFSZ + : there is at least one heavy (with mass of order V PQ ) colored fermion carrying PQ charge, and thus we also have UV dominated production at T R from gluino scattering off quarks and gluons.
In DFSZ 0 , an IR contribution to axino production via scattering also arises from the supersymmetrized aGG operator generated when the top quark is integrated out [16], but JHEP03(2017)005  For UV production in DFSZ + we fix N DW = 6. Vertical lines correspond to axino freeze-in via decaysχ →ã, followed bỹ a →G a, whereas horizontal lines correspond to UV gravitino production at T R . The thick (thin) portions of the lines refer to a freeze-in contribution smaller (larger) than 50%.
it is suppressed compared to the one from decays and we neglect it. In theories with a low gravitino mass, the decay of neutralinos can also populate gravitinos by the FI mechanism but we find this contribution sub-dominant to the ones mentioned above.
In the absence of a saxion condensate, UV production of both axinos and gravitinos puts a very powerful bound on T R [17]. This is illustrated in figure 1 both for DFSZ 0 (left panel) and DFSZ + (right panel) for gravitino LSP and other superpartners at the TeV scale. Contours that yield the DM relic density are shown in the (V PQ , T R ) plane for four values of the gravitino masses. Even for a gravitino with a weak scale mass, the reheat temperature after inflation is strongly bounded, T R 10 8 GeV. Thus LSP dark matter is typically overproduced in the absence of a saxion condensate, unless V PQ is very large and T R is very low. However, if V PQ is very large so that axino freeze-in is significantly suppressed, the universe is typically overclosed by axions, unless a low value of the axion misalignment angle is selected by an anthropic requirement.
Saxion cosmology can greatly change this conclusion. Axion theories typically have a domain wall problem, which we assume is solved by breaking the PQ symmetry before inflation, and not restoring it afterwards. We define the saxion field so that today the saxion vev is zero. During inflation supersymmetry breaking yields a potential for the saxion, displacing the vev away from today's value. Depending on the sign of the quadratic term, the vacuum value s I is either V PQ or of order the cutoff of the field theory, M * JHEP03(2017)005 V PQ [18,19]. 2 If s I 10 13 GeV or s I ∼ M * , this saxion condensate comes to dominates the energy density of the universe, producing an early matter-dominated (MD) era. When the saxion condensate decays, large entropy is created that has a key effect on both LSP and axion contributions to dark matter. The conventional picture of dark matter survives only for the restricted case of s I ∼ V PQ 10 13 GeV, when the saxion condensate never dominates, allowing the usual favorite cases of LSP freeze-out or axion misalignment with V PQ ∼ 10 12 GeV. However, as we have already mentioned, this case of low V PQ has a serious cosmological problem. Even if we choose T R low enough to suppress UV production, there is still the IR contribution from freeze-in, which cannot be suppressed unless we select T R below the superpartner mass scale. Hence we consider larger s I , giving a saxion condensate MD era that has important consequences for dark matter abundance.
In KSVZ [20,21] axion theories the MSSM µ term is not forbidden by PQ symmetry and the axion multiplet does not couple to the Higgs superfields. Thus the saxion typically has a large decay branching ratio to axions and the decay of the saxion condensate gives an axion contribution to dark radiation that is excluded [22]. Hence we focus on DFSZ theories.
In our scheme the dominant saxion decay is to pairs of Higgs, W or Z bosons, giving a very low reheat temperature of the saxion 3 T Rs ; for example, T Rs ∼ GeV-MeV for V PQ ∼ (10 13 −10 16 ) GeV, respectively. The decay of the saxion condensate has crucial implications for dark matter: • The LSP abundance from freeze-out is greatly diluted [23,24]. The freeze-out mechanism can only give the observed dark matter abundance if it first overproduces LSPs that are then diluted. This could happen by raising superpartners well above the TeV scale.
• For V PQ > 10 13 GeV, the entropy is released after the axion condensate starts to oscillate, diluting misalignment axions. For a misalignment angle of order unity, the value of f a needed for axion dark matter is raised from ∼ 10 12 GeV to ∼ 10 15 GeV [25], suggesting the possibility that PQ and grand unified gauge symmetries are broken together [26].
• The large abundances of gravitinos and axinos produced at T R , and of axinos produced by freeze-in at the TeV scale, are diluted by saxion decays, greatly weakening the constraints of figure 1 and allowing much higher T R and lower V PQ .
• If kinematically allowed, the saxion condensate leads to a late production of superpartners, and therefore LSPs, via such decays s →ãã,ãG,χχ.
Given these crucial effects of the saxion condensate, what are the leading candidates and production mechanisms for dark matter? Misalignment axions with V PQ the scale of grand unification becomes one attractive option [26]. In this paper we consider LSP possibility by assuming that the axion abundance after dilution is sub-dominant. A misalignment angle 2 The possibility sI = 0 requires a special symmetry structure that we do not consider in this paper. 3 We define TR and TRs as the reheat temperatures after inflation and saxion decays, respectively.

JHEP03(2017)005
3) allows V PQ as high as 10 16 GeV. We work in a generic framework of TeV scale supersymmetry, taking the saxion and axino masses of order the TeV scale, as well as all other superpartners except, perhaps, the gravitino which could be much lighter. Freeze-out gives too low an abundance for a typical TeV-scale superpartner spectrum, although a bino-like LSP is a possibility providing the condensate is not too large. On the other hand, if s decays to superpartners the late decay of the condensate typically overproduces LSP dark matter. Hence we assume these decays are kinematically forbidden and focus on two mechanisms for LSP production: IR freeze-in of axinos from neutralinos and charginos,χ →ã, and UV scattering at T R producing axinos and gravitinos, gg →ã,G. We consider two possibilities for the LSP:ã andG, and perform a systematic analysis for the dark matter abundance over a very wide range of T R and V PQ , identifying regions of both cold and warm dark matter. A key feature of our work is to connect the dark matter cosmology with displaced vertex signals, arising from superpartner production and decay to gravitinos [27] and axinos [28], at the LHC and future colliders.
Another aspect of the gravitino problem is the powerful constraints from Big Bang nucleosynthesis (BBN) arising from decays between the lightest observable sector superpartner (LOSP) and the gravitino [29]. We avoid this by making the axino and gravitino the two lightest superpartners. The LOSP then decays before BBN and if the gravitino is the NLSP it decays toãa, which has no effect on BBN.
There is a considerable literature on the effects of a saxion condensate on LSP abundances. The gravitino and axino problems were solved by the decay of a saxion condensate in ref. [30], which also considered the possibility of gravitino dark matter from inflaton decay. This was further developed by identifying PQ-breaking fields as the waterfall fields of a hybrid inflation model with vevs of order 10 15 GeV [31]. An alternative solution to the gravitino problem involved a very light, keV axino [32,33]. In other work a saxion condensate was used to obtain a PQ breaking scale as large as the unification scale in theories with axino or neutralino LSP [34], and further work considered mixed axion/neutralino dark matter [35]. A KSVZ scheme with a saxion condensate relevant for gravitino dark matter, but at lower values of V PQ than considered in this paper, is given in ref. [36].
We describe the saxion condensate MD era in section 2, by identifying the characteristic temperatures associated to this epoch and giving analytical expressions for them. The production mechanisms for axino and gravitino are discussed in detail in section 3, where we account for the dilution for the saxion condensate and give numerical results for the yields Yã and Y 3/2 as a function of V PQ . In section 4, we discuss key consequences of our scheme with the axino and gravitino as the two lightest superpartners. The NLSP decay leads to a component of dark matter that is warm, with free streaming length as illustrated in figure 7. It is intriguing that this result can be consistent with Lyman-α forest observations while at the same time solving issues with cosmological structures at small scales. Displaced signals at colliders resulting from the decay of the LOSP to axinos and gravitinos are discussed, as well as the axion dark radiation resulting from decay of the saxion condensate.
In sections 5 and 6 we present our results for two classes of SUSY spectra: high scale (i.e. "gravity") and low scale (i.e. "gauge") mediation, respectively. The LSP relic density JHEP03(2017)005 for high scale mediation is shown in figures 9 and 11, whereas the analogous results for low scale are in figures 13 and 15. A remarkable signal of our framework, which holds regardless of the mediation scale, is displaced events at colliders. We study how these lifetimes vary with the supersymmetry breaking parameters in figure 8, and make lifetime predictions in figures 10, 12, 14, and 16. We supplement our work by appendices with useful results employed in our analysis.

Saxion cosmology
Inflationary dynamics sets the initial conditions for the saxion condensate, typically displacing it by an amount s I from the minimum today. After inflation ends, Hubble friction keeps the saxion field fixed until the universe expands and cools sufficiently that 3H ∼ m s , with m s the saxion mass, when the saxion condensate oscillates around its minimum. The energy density stored in such oscillations red-shifts like non-relativistic matter at a rate slower than radiation, and eventually dominates the energy budget.
The specific temperature where saxion oscillations begin is crucial for the evolution of the saxion condensate. In particular, we identify two main regimes according to the size of the reheat temperature after inflation and describe the cosmologies separately. In the case where T R 10 10 GeV, saxion oscillations start during the radiation-dominated era after inflationary reheating has ended, while in the case where T R 10 10 GeV, they start during inflationary reheating. In addition to the saxion condensate that arises from inflationary dynamics, the saxion can also be generated from the thermal processes such as the scattering with gluinos via the dimension-5 operators generated by the heavy colored fermions with PQ charges in the DFSZ + theories. This contribution is sub-dominant when s I > 10 13−14 GeV. For s I < 10 13−14 GeV, we find that the extra thermal saxions still fail to provide the dilution necessary for the overproduced axinos from UV scattering. Therefore, this thermal contribution of saxions is never relevant in the parameter space of interest.
2.1 High reheat temperature after inflation: T R 10 10 GeV Saxion oscillations begin when the radiation bath has a temperature T osc = 10 where we introduce the effective number of relativistic degrees of freedom in the thermal bath g * (T ) and the reduced Planck mass M Pl = 2.4 × 10 18 GeV. The onset of the saxion oscillation happens during an early Radiation Dominated era (RD ). The energy stored in saxion oscillations red-shifts with the scale factor a as a −3 , namely with the same behavior of non-relativistic matter. If the condensate is long-lived enough, it eventually takes over the radiation energy and the universe enters an early matter-dominated (MD) era at a JHEP03(2017)005 temperature T M . This characteristic temperature is found by imposing that the saxion and the radiation bath equally contribute to the energy density, or in other words by solving the equation where we use the conservation of the total entropy to properly red-shift the saxion energy. Only for the purpose of an analytical estimate, we set g * (T osc ) = g * (T M ) = 228.75, corresponding to the full MSSM field content, giving (2.4) This early MD era consists of two distinct phases, as detailed in refs. [8,26]. In the beginning, the radiation energy density is dominated by the red-shifted initial component. We dub this initial part an adiabatic matter-dominated (MD A ) era. However, radiation constantly produced from saxion visible decays red-shifts slower and in the end accounts for most of the radiation present in the universe. Once the contribution from saxion decays dominates, at a temperature T N A , we enter a non-adiabatic matter-dominated (MD N A ) era, where the total entropy is not conserved. At temperatures below T N A , saxion decays reheat the universe and a large amount of entropy is released. Finally, most of the saxions decay when the Hubble parameter is of the order the saxion decay width, H ∼ Γ s . This is the beginning of the last phase of a radiation-dominated universe (RD), which starts when the radiation bath has the reheat temperature T Rs .
We require that saxion decays to R-odd neutralinos and charginos are kinematically forbidden, and that decays to axions are sub-dominant. The saxion width is thus dominated by decays to MSSM Higgs bosons and longitudinal electroweak gauge bosons. We use the expression in eq. (B.8) for the saxion visible width, valid in the decoupling limit and for large tan β. The transition from MD A to MD N A occurs at temperature   there is no significant entropy production at T Rs . We define V A physically meaningful and useful quantity is the amount of entropy released by saxion decays during the reheating process. We quantify this by introducing the dilution factor D(T i ), defined as the ratio of the entropy after saxion decays to that at an initial temperature T i (2.7) We interpret T i as the temperature when dark matter is produced, e.g. the freeze-in temperature. In the MD N A era, we make use of the scaling a 3 ∝ T −8 , whereas in the MD A and RD eras, Rs , we obtain the following dilution factor   The dilution factor is extremely large for the initial condition s I = M P l , ranging from 10 9 to 10 15 as V PQ increases from 10 10 GeV to 10 16 GeV. On the other hand, the s I = V PQ case has a much small dilution, varying from 1 to 10 9 as V PQ increases from 10 13 GeV JHEP03(2017)005 to 10 16 GeV. Each comoving number density frozen-in or -out before the entropy release, namely at temperatures T i > T N A , gets maximally depleted by D.
The quantity D also allows us to quantify the critical value V . (2.10) 2.2 Low reheat temperature after inflation: T R 10 10 GeV So far we have assumed that saxion oscillations begin during a radiation dominated era. However, this need not be the case as we do not know the temperature of reheating after inflation, which is set by the decay width of the inflaton. We now extend our analysis to low T R . During inflationary reheating, the temperature dependence of the Hubble parameter is different than during a radiation dominated era. An approximate solution to the coupled Boltzmann equations describing the evolution of inflaton and radiation energy densities, at early times before the inflaton decays, gives the expression for the Hubble parameter [23,24] Oscillations for the saxion condensate begin at T osc , when 3H m s . Assuming that T osc is less than the maximum temperature reached during this reheating era, so that H is now given by eq. (2.11), we obtain This equation makes sense only if T osc > T R , and this condition is satisfied as long as T R < T osc 10 10 GeV, namely the oscillation temperature given in eq. (2.2).
The saxion oscillations for T < T osc still red-shift as non-relativistic matter, and once the inflaton finally decays the energy density stored in them results in Afterwards, we have the conventional thermal history described above. The only difference is the initial condition for the saxion energy density as in eq. (2.13). This in turn gives a different expression for T M , which is now found by imposing the condition (2.14) -9 -

JHEP03(2017)005
The solution of this equation simply gives where in the last step we have used the expression in eq. (2.12). This expression breaks down when T M becomes larger than T R , which happens if s I ≥ M P l / √ 3. We avoid this situation because the saxion condensate dominates the energy density before it starts to oscillate and consequently the Universe enters inflation.
The expression for T Rs is of course unchanged, and still given by eq. (2.6). However, the temperature for the transition to the non-adiabatic phase is changed 632 GeV T R 10 10 GeV .

(2.16)
Moreover, also the expression for the dilution factor changes We observe the remarkable fact that the dilution factor is proportional to T R when T i ≥ T N A . The three characteristic temperatures T M , T N A and T Rs are shown in figure 3.

Field equations
The analytical expressions for the characteristic temperatures and the dilution factor are very useful for an order of magnitude estimate of the effect. However, they are not suited for precision calculation of relic densities. In this work we solve numerically the Boltzmann equation system describing the evolution of the saxion condensate coupled to the radiation bath.
We begin with the case where the saxion starts to oscillate after inflationary reheating, presented in section 2.1. The energy density of the saxion condensate, after the onset of oscillations at T osc of eq. (2.2) or at T osc of eq. (2.12), evolves according to The redshift due to the Hubble expansion is accompanied by the term proportional to the saxion total decay width. The radiation bath temperature evolves according to  where g * is the effective number of degrees of freedom contributing to the entropy density. In the limiting case where g * is a constant, the equation takes a more familiar form where the radiation energy density results in This approximation is certainly valid at very high temperature, where the full spectrum is relativistic. At lower temperature the approximation breaks down and the error one makes by using eq. (2.21) is at most of few percent. In our work we use eq. (2.19), and g * (T ) is computed using the masses of the SM particles and of the SUSY particles, assumed degenerate at 1 TeV. Finally, the time-temperature relation can be found by solving the Friedmann equation The initial condition for this case is set at some high temperature T 0 by where T M is defined in eq. (2.4) and our numerical studies are completely insensitive to T 0 .

JHEP03(2017)005
For the case where the saxion starts to oscillate during inflationary reheating, the numerical setup needs to be extended as follows. Firstly, we cannot ignore the inflaton anymore and we couple eqs. (2.18) and (2.19) to the one describing the evolution of the inflaton energy density ρ I , which takes the same form as eq. (2.18) with the modifications ρ s → ρ I and Γ s → Γ I . Furthermore, we need to add the inflaton decay contribution Γ I ρ I to the right-hand side of eq. (2.19), and the inflaton energy density on the right-hand side of eq. (2.22). Secondly, we set the initial conditions for the saxion oscillation at the time 3H(t osc ) = m s , with H in this case dominated by the inflaton energy density. Since inflationary reheating is in an MD era, we can identify 3H(t osc ) = 2/t osc , and thus the initial condition ρ s (t osc ) = m 2 s s 2 I is set at the time t osc = 2/m s .

Axino and gravitino production
In this section we quantify axino and gravitino production by accounting for the saxion condensate effects. We consider three different mechanisms: axino production from freezein, axino UV production (present only in DFSZ + theories) and gravitino UV production.
For each case, we show results for the comoving number density as a function of V PQ , and we consider both s I = V PQ and s I = M * . We also comment on gravitino freeze-in production, a sub-dominant source. The results presented here are completely general and independent of where the axino and gravitino sit in the superpartner spectrum. We apply the framework to two specific spectra corresponding to High Scale and Low Scale mediation in sections 5 and 6 respectively.

Freeze-in production of axinos
The freeze-in production of axinos is controlled by neutralino and chargino decays and inverse decays,χ →ã. Explicit decay widths relevant for this case are given in appendix B. The evolution of the axino number density nã is governed by the Boltzmann equation Our goal here is to provide the expression for the collision operators C FI . For a light axino, lighter than all the MSSM superpartners, freeze-in comes from neutralinos and charginos decays to axinos The partial widths for these channels are given in eqs. (B.14) and (B.16), respectively. On the contrary, for a heavy axino, heavier than all the neutralinos and charginos, freeze-in production is dominated by the inverse decays

JHEP03(2017)005
We use the detailed-balance principle to write the collision operator in the Boltzmann equation by using the inverse reaction, namely the axino decay, with partial decay widths given in eqs. (B.15) and (B.17). In the intermediate case, with the axino mass within the neutralinos and charginos masses, we have both types of reactions, but only the ones allowed by kinematics. The freeze-in collision operator is a sum of the possible sources The first two terms account for the processes involving the neutralinos where K 1 is the first modified Bessel function of the second kind. The Heaviside step function θ makes sure that only the kinematically allowed channels are accounted for. The correspondent collision operators for the charginos are The results for the axino comoving density are shown in figure 4, where the saxion dilution is computed using the cosmology of section 2.1 with T R 10 10 GeV. In each of the two panels, we fix the MSSM mass parameters as in the caption, and compute the axino comoving density as a function of V PQ .

UV production of axinos
This gravitino problem is greatly exacerbated in PQ theories with heavy matter since then UV production of axinos also generally occurs. The combination of the two UV production mechanisms provides a particularly powerful upper bound on T R [17]. A complication in computing the UV contribution to axino production is that if the heaviest matter carrying both PQ charges and gauge charges, Φ, has a mass M Φ < T R then the UV production is cutoff at M Φ [16]; thus Y UṼ a is model-dependent. In the DFSZ + theory we take M Φ T R , so that the UV axino production is cutoff at T R . The axino mass is expected to be of order the gravitino mass or larger, in which case, in DFSZ + with V PQ ∼ 10 12 GeV and TeV scale supersymmetry, T R 10 6 GeV [17]. Here we show that the limits from UV production of axinos is greatly ameliorated by the decay of the saxion condensate, with resulting limits depending on (V PQ , mã).   With UV production of axinos cut off at T R , scattering leads to an axino undiluted yield [12] Y UṼ a | + ∼ 2 × 10 −6 g 6 3 ln

JHEP03(2017)005
where g 3 is the strong gauge coupling. As elsewhere, we take N DW = 6. The yields of axinos and gravitinos should not exceed the thermal equilibrium value where ζ is the zeta function and the internal degrees of freedom g is 2 (4) for the axino (gravitino).
To include the effect of saxion dilution, we divide the axino abundance by the dilution factor computed numerically. Since axino UV production occurs before the saxion injects the entropy, the dilution factor is simply the ratio of the total entropy before and after the saxion MD era. The analytic formulas for the dilution factors are also given in eqs. (2.9) and (2.17).
The results for the axino abundance are shown in figure 5, for some values of T R 10 10 GeV that have saxion dilution of section 2.1 and others of section 2.2 with T R 10 10 GeV. We take s I = V PQ (M * = 3 × 10 16 GeV) in the left (right) panel. Note that for a sufficiently low V PQ , the axino reaches thermal equilibrium, in which case we will use the equilibrium value of the yield in eq. (3.12). In particular, the sharp kink in each of the curve is due to the transition from the equilibrium value Y eq to the yield from UV scattering. In the left panel, the smooth change of the slope at V PQ ∼ 10 13−14 GeV is due

JHEP03(2017)005
Axino UV Production T R (GeV)  to the emergence of the saxion MD era demonstrated in figures 2 and 3. In the right panel, the dilution effect is present for the entire range of V PQ . In the case where T R 10 10 GeV, one interesting feature is that the diluted abundance is (almost) independent of T R because both dilution in eq. (2.17) and production in eq. (3.11) are proportional to T R (other than the mild dependence on T R in g 3 ).

UV production of gravitinos
In supersymmetric theories UV production of gravitinos generally limits the reheat temperature after inflation, T R [13][14][15]. The undiluted abundance of gravitinos from scattering at T R is [37,38] Y UV 3/2 ∼ 6 × 10 −12 is the gaugino mass of (U(1), SU(2), SU (3)) and the A-term of the top Yukawa coupling. Similar to section 3.2, we calculate the final gravitino abundance by dividing the yield by the dilution factor computed numerically. The analytic estimate of the dilution factors are also provided in eqs. (2.9) and (2.17).
The numerical results for the gravitino abundance are shown in figure 6 for T R ≥ 10 10 GeV (T R ≤ 10 10 GeV) with saxion dilution of section 2.1 (section 2.2). In the left panel, with PQ because the condensate is too small for the saxions to dominate before they decay. The key feature is the rapid dilution for V PQ > V PQ , leading to dilution at low V PQ . For T R < 10 10 GeV, the final diluted yield is nearly independent of T R (other than γ i (T R ) in eq. (3.13)).

Freeze-in production of gravitinos
Gravitinos are also produced via decays of thermal charginos and neutralinos, N i / C i → G [39]. These freeze-in processes are IR dominated and independent of T R . The resulting yield, which is proportional to the decay width to the gravitino, is enhanced for low values of the gravitino mass. This is just a consequence of the production of longitudinal gravitinos (goldstinos), as manifestly shown in eqs. (B.25)-(B.29). For the same reason, gravitino UV production is also enhanced for small m 3/2 , as shown in eq. (3.13). We checked that the gravitino FI contribution is always sub-dominant compared to those from gravitino UV and axino FI in the parameter space of interest, and therefore we do not consider it in this work.

Axino and gravitino as the lightest superpartners
We classify the superpartner spectra according to the size of the mediation scale for supersymmetry breaking, M mess . Superpartner masses arise from the effective operators where the superfield X, defined in eq. (A.9), has a SUSY breaking F-component. We take the model-dependent coefficients c a and c Q to be comparable c Q c a ; the spectrum is not split.

JHEP03(2017)005
In section 5 we consider M mess of order M Pl , which can be broadly identified with "gravity mediation." The gravitino for this case cannot be much lighter than the other superpartners. This has to contrasted with the low mediation scale case elaborated in section 6, M mess M Pl , where the gravitino is much lighter than other superpartners. We focus our attention on the spectra where the axino and the gravitino are both lighter than all the other superpartners. However, for gravity mediation we do not commit to any relative hierarchy between them, and consider both axino and gravitino LSP cases.
For both high and low mediation scales, axinos and gravitinos are produced in the early Universe through the various mechanisms discussed in the previous section. These processes produce both the NLSP and the LSP, and they all contribute to the present dark matter abundance, which today is made of LSP particles only. Highly relativistic axions are produced in NLSP decays, which make a negligible contribution to dark radiation. Furthermore, since the final products are the LSP and the axion, these decays are not subject to BBN limits [40][41][42]. As we will see in this section, the origin of the current dark matter abundance is typically due to either NLSP or LSP production, unless we consider peculiar parameter space regions.
Before we discuss the high and low scale cases in detail, we highlight the main features of our framework.

Warm dark matter from NLSP decays
Dark matter from LSP production is always cold. On the contrary, dark matter from NLSP production and decay could be hot, warm or cold, depending on the ratio of NLSP and LSP masses and the NLSP decay lifetime. The LSP is a gauge singlet extremely weakly coupled to the radiation bath, and thus DM particles coming from NLSP decays just lose their momenta by pure free streaming. This has the effect of potentially washing out cosmological perturbations at small scales through free streaming, and this scenario is severely constrained by Lyman-α forest observations [43][44][45], bounding the free streaming length to be less than 1 Mpc. Interestingly, if the free streaming length is consistent with large scale structure and not too small it can address some large scale structure (LSS) issues that are indicated by simulations of collisionless cold dark matter [46]. Baryonic feedback effects can explain some discrepancies [47,48], although there is much debate on this [49,50].
If the axino is the NLSP, it decays to longitudinal gravitinos (i.e. goldstino) via the effective operators given in eq. For gravitino NLSP, the decays to axinos are mediated by the operator in eq. (B.18). We assume that the decay to saxion and axino final state is kinematically forbidden, and thus the decay width is half of the expression in eq. (B.20), to account for decays to axion and axino only. The resulting lifetime is To simplify the discussion we take the axino mass to be of order    The full result is shown in the right panel of figure 7. In the mã m 3/2 region, the free streaming isocontours are along the lines m 3/2 m 2 a = const, consistently with the approximate expression in eq. (4.5). Lyman-α bounds are evaded for sufficiently heavy gravitino and/or axino, and issues with LSS can be addressed again by TeV scale superpartners.

B-like H-like
A related noteworthy example is the case of gravitinos coming from LOSP freeze-out and decays [51,52], which in our case is made irrelevant by the saxion condensate dilution.

Displaced signals at colliders
If the lightest observable-sector supersymmetric partner (LOSP) is a neutralino N 1 , its decay to the axino and Higgs/longitudinal Z bosons will leave missing energy and a displaced vertex. The "lifetime" for this channel can be obtained from the associated decay width given in eq. (B.14) of appendix B, (4.6) Here, we define the mixing factor kã ≡ |s β R 31 | 2 + |c β R 41 | 2 , with R ij the neutralino mixing matrix (for details see eqs. (B.10) and (B.14)). In the pure Higgsino LOSP limit, with decoupled bino and wino, kã becomes 1/2 and is thus independent of β. On the other hand, kã will be suppressed when the LOSP is bino-or wino-like. For Low Scale mediation with m 3/2 1 MeV, the decay channel of neutralinos to gravitinos becomes sufficiently enhanced to dominate over the axino final state, giving a JHEP03(2017)005 The lifetime that can be probed at the LHC and future colliders depends on the total production cross section of supersymmetric particles, which we quote for the case of degenerate squark and gluino masses,m [53]. Form = (1.5, 2.5) TeV, at √ s = 14 TeV this cross section is of order (100,1) fb, so that planned runs of the LHC will allow ATLAS and CMS to reach cτ of order (100, 10)m. Recently, a surface detector called MATHUSLA [54] has been proposed to search for (ultra) long-lived particles at the LHC and future colliders. In sections 5 and 6 we show predictions for cτ i , (i =ã,G), following from the constraint Ωh 2 = 0.11, in theories with High and Low Scale mediation for the particular point in supersymmtric parameter space of (µ, M 1 , tan β) = (1 TeV, 1 TeV, 2). Here we illustrate the variation in the cτ i as the parameter space changes, always keeping the unified gaugino mass relation, M 2 2M 1 . In figure 8, we display curves for cτG in red, labelled by values of V PQ . The curves for cτã are in purple when the LOSP is Higgsino-like and in orange when it is bino-like, labelled by values of m 3/2 . The LOSP lifetime cτÑ 1 is given approximately by the smallest cτ i . In the left (right) panel, we vary µ (M 1 ) while fixing M 1 (µ). The change of behavior in the lifetime curves at µ = 2 TeV (M 1 = 2 TeV) in the left (right) panel reflects the fact that the mixing factors in eq. (4.6) and eq. (4.7) drastically change as µ becomes larger or smaller than M 1 . The neutralino LOSP decay can lead to observable displaced signals at the LHC and future colliders over a remarkably wide range of parameter space.

Axion dark radiation
Saxions can decay to axions with a rate given by eq. (B.2) if κ does not vanish due to symmetry. Using the branching ratio of the saxion to the visible sector and to axions, we predict the amount of dark radiation to be [26]  (4.8) The Planck experimental bound [57] is ∆N eff < 0.6. The proposed CMB Stage-IV experiment [58] can be sensitive to ∆N eff = 0.03.
5 Results for high scale or "gravity" mediation  occurs via decays of neutralinos and charginos with rates proportional to 1/V 2 PQ given by eqs. (3.7) and (3.9). On the other hand, gravitinos are populated by the UV scattering of quarks, gluons, and their superpartners with an abundance given by eq. (3.13), which is proportional to T R . Both production sources are heavily diluted by the decay of the saxion condensate, which also makes the LOSP freeze-out and decay contribution negligible.
We compute the total DM abundance from these LSP and NLSP production sources in terms of V PQ and T R and draw contours of Ωh 2 = 0.11 in the (V PQ , T R ) plane in figure 9. After decaying to the LSP, the NLSP number density is transferred to that of the LSP, and therefore we only need the LSP mass to compute the final DM abundance. In the left two panels, we take s I = V PQ , and in the right two panels we take s I = M * and give contours for M * = (10 15 , 3 × 10 16 , 10 18 ) GeV. The top two panels have T R ≥ 10 10 GeV and are therefore described by the cosmology in section 2.1, while the bottom two panels have T R ≤ 10 10 GeV and are described by the cosmology in section 2.2.
In the upper two panels, for each value of s I two contours are shown. The blue one is an example of an axino LSP, while the red one is an example of gravitino LSP. The vertical parts of the contours have axino FI as the dominant production mechanism and hence are independent of T R , while the parts of the contour with positive constant slope have UV gravitino production dominate. When freeze-in occurs above T N A , the dilution factor from the saxion condensate is proportional to s 2 I V PQ , which is much larger in the right panel than in the left panel. This means that the production needs to be much larger in the right panel than the left, and hence the contours in the right panel are at much lower V PQ . This also explains why in the right panel, as s I is increased from 10 15 GeV to 3 × 10 16 GeV, the contours move to lower V PQ . However, for s I > 10 16 GeV a new regime is entered where freeze-in occurs during the MD N A era and dilution becomes independent of T M and therefore of s I . Hence, for s I = 3 × 10 16 GeV and 10 18 GeV, the blue and red vertical contours differ slightly only because of a different m LSP .
When UV gravitino production dominates, the difference in the slopes of the contours results from the dilution factor, scaling as V 3 PQ in the left panel and V PQ in the right panel. At large T R where UV gravitino production dominates, the blue contours are well above JHEP03(2017)005 the red ones; this is because Y 3/2 ∝ 1/m 2 3/2 and the blue contours have larger m 3/2 and hence need larger T R to compensate. On the other hand, at lower T R where axino FI dominates, Yã is independent of mã and m 3/2 , so that the red and blue Ωh 2 contours differ only because of m LSP .
Each contour is divided into thick and thin parts. The thick parts indicate cold dark matter (CDM) where NLSP production is sub-dominant. The thin parts label the warm dark matter (WDM) case where the component of dark matter from the NLSP decay constitutes more than 50% of the total abundance. The thin lines may be relevant for understanding possible difficulties with pure CDM, as in core-cusp, too big to fail and missing satellite problems.
A key point emerges from comparing the contours in the two upper panels of figure 9 with the contour for High Scale mediation (m 3/2 = 100 GeV) in figure 1, where the saxion condensate is absent. The saxion condensate increases the maximum value of T R from 10 8 GeV to 10 16 GeV. This allows very high reheat temperatures after inflation, 4 so that baryogenesis may occur at very high temperatures. The upper bound on T R now arises because inflation only gives a saxion condensate if PQ symmetry is broken before inflation. The precise constraint on T R is dependent on the model for the PQ phase transition and on the model for reheating after inflation. For instantaneous reheating we expect the condition to typically be T R < V PQ corresponding to the unshaded region of figure 9. However, certain theories may have PQ breaking before inflation even for T R somewhat higher than V PQ , so that for these theories the lightly shaded region also becomes physical. We expect the dark shaded region to be unphysical in all models. On the other hand, if the inflaton decay rate is slow, T R could be much below the energy scale of inflation, lowering the shaded bands and reducing the maximal allowed value of T R . A reduction by a few orders of magnitude would still allow T R to be sufficiently large for leptogenesis.
Another key point emerges from comparing the contours in the two upper panels of figure 9 with the contour for High Scale mediation (m 3/2 = 100 GeV) in figure 1. For large s I = M * , the saxion condensate lowers the minimal value of V PQ by several orders of magnitude. An important part of the axino problem is that the minimal value of V PQ is so large -4 × 10 14 GeV for m 3/2 = 100 GeV -that misalignment axions overclose the universe unless the misalignment angle is very small. This difficulty is removed with a saxion condensate because the decay of the condensate also dilutes misalignment axions, so that they typically give sub-dominant contributions to dark matter for f a 10 15 GeV and misalignment angles order unity [19,26]. Since V PQ is larger than f a by N DW / √ 2, which is often an order of magnitude, axions are typically sub-dominant on all the contours of figure 9, although they could give a comparable contribution when V PQ is of order (10 15 −10 16 ) GeV.
In the lower two panels, for s I = 10 15 GeV and 3 × 10 16 GeV, freeze-in occurs above T N A , so that the dilution factor is proportional to T R , as seen in eq. (2.15), and the axino FI abundance depends on T R . In particular, a decrease in T R is compensated by an increase in JHEP03(2017)005 V PQ at a rate that depends on the cosmological era during FI. However, for s I = 10 18 GeV much of the contour is still vertical as freeze-in occurs during MD N A . For very high s I there is a very robust prediction for V PQ from axino freeze-in dark matter.
If the LOSP is a neutralino, its lifetime, eq. (4.6), can be predicted by determining V PQ from the dark matter abundance. The prediction shown in figure 10 is obtained by inverting the two axes in the right panels of figure 9 and converting the V PQ axis to the lifetime using eq. (4.6). The left (right) panel corresponds to the cosmology of T R less (greater) than 10 10 GeV discussed in section 2.2 (section 2.1). As explained for figure 9, for axino FI a larger saxion condensate (higher s I ) leads to a lower V PQ , which in turn gives a shorter LOSP lifetime. However, this behavior does not continue for an arbitrarily high s I : once s I is sufficiently high, the axino freeze-in occurs during the MD N A era and the abundance becomes insensitive to s I [8]. Consequently, for T R > 3 × 10 9 GeV and any s I > 10 16 GeV, there is a very robust prediction of cτ LOSP 10m. In the left panel, as T R drops FI transitions to occurring in MD A , so that the dilution factor also drops and V PQ and cτ LOSP increase. As in figure 9, the thin parts of contours give warm dark matter from NLSP decay.
We do not show the LOSP lifetime for s I = V PQ . The large values of V PQ indicated by the left panels of figure 9 lead to cτ ∼ 10 (5−7) m.

The DFSZ + theory
In addition to the axino FI and gravitino UV contributions existing in DFSZ 0 theories, there is also axino production from UV scattering discussed in section 3.2 for DFSZ + theories. In a setting similar to figure 9, we show the results for the total abundance in figure 11 as a function of T R and V PQ . Gravitino production is everywhere sub-dominant. In fact UV axino production dominates everywhere, except for T R less than about 10 6 GeV where axino freeze-in becomes important. The blue and red curves, corresponding to axino and gravitino LSP, nearly coincide; they differ only because the LSP mass differs by a factor of two between the curves.
The parametric dependence of the contours can be understood from Ωh 2 ∝ Yãm LSP /D, where, for UV production, the dilution factor D ∝ s 2 I V PQ (1, T R ) for (T R 10 10 GeV, T R 10 10 GeV). For UV axino production, Yã ∝ (T R /V 2 PQ , 1), with the constant value applying only if the equilibrium abundance is reached, which happens above and to the left of the dashed line. For T R > 10 11 GeV, the contour for s I = 10 18 GeV is vertical, corresponding to an equilibrium axino abundance, while the contours with s I = 10 15 , 3 × 10 16 GeV have T R ∝ V 3 PQ , and those with s I = V PQ are steeper with T R ∝ V 5 PQ . For 10 6 GeV < T R < 10 10 GeV, both Yã and D are linear in T R , so that the contours are vertical. At very low T R the contours become sloped as the axino yield becomes dominated by freeze-in; in this low T R region they are identical to the curves in figure 9.
In figure 12, we predict the neutralino LOSP lifetime, eq. (4.6), from the values of V PQ fixed by the observed dark matter abundance, in a way completely analogous to figure 10, with the left (right) panel for T R less (greater) than 10 10 GeV. The key feature in DFSZ + is that, for all T R > 10 6 GeV, UV axino production dominates dark matter production and hence, compared to DFSZ 0 , more dilution (i.e. higher s I ) is required for a given V PQ . Figure 12 shows that for s I of order, or larger than, the scale of supersymmetric grand 6 Results for low scale or "gauge" mediation neutralinos decay to gravitinos with an inverse partial width of order (m 3/2 / GeV) 2 sec as given in appendix B.4, so that for m 3/2 1 GeV decays to gravitinos occur before BBN and this part of the gravitino problem disappears. For such light gravitinos, our results apply even if the axino is not the NLSP, increasing the generality of the predictions for gravitino dark matter. For m 3/2 > 1 GeV, a non-NLSP axino has a sizable branching ratio tohh, even for the highest values of V PQ of order 10 15−16 GeV, so that in this case the axinos must be the NLSP.

The DFSZ 0 theory
Dark matter production in DFSZ 0 theories is dominantly from axino freeze-in and Gravitino UV production as explained in section 5.1. Here we are concerned with Low Scale mediation of supersymmetry breaking, as in gauge mediation, and hence take the gravitino to be the LSP. Axinos produced from freeze-in subsequently decay and produce a warm component of gravitino dark matter.
In figure 13, we show contours of Ω 3/2 h 2 = 0.11. Following the same setup as in figures 9 and 11, the thick (thin) parts of the curves indicate that dark matter arises dominantly from LSP (NLSP) production. The gray region is excluded because PQ breaks after inflation. The upper (lower) panels are for the cosmology of T R greater (lower) than 10 10 GeV discussed in section 2.1 (section 2.2). The left (right) panels assume s I = V PQ (M * = 3 × 10 16 GeV). It is important to remember that the dilution factor is larger for larger V PQ .  Figure 13. Contours of Ωh 2 = 0.11 from axino freeze-in and gravitino UV production. We fix q µ = 2, D = 4, tan β = 2, M 2 /2 = M 1 = µ = 1 TeV, and m s = 600 GeV. The top (bottom) row is for the cosmology with T R ( )10 10 GeV discussed in section 2.1 (section 2.2).
Gravitino UV production given by eq. (3.13) is much enhanced for low m 3/2 . For low enough m 3/2 and high enough T R , gravitino UV production is so efficient that gravitinos thermalize, with an undiluted yield independent of m 3/2 and T R given in eq. (3.12). With dilution in eq. (2.9), the final abundance reads In the upper panels, the vertical thick lines correspond to this equilibrium gravitino production. As T R decreases, the gravitino may become non-equilibrium with an abundance decreasing with T R , as in eq. (3.13), corresponding to the sloped parts of the thick curves in the upper panels. The diluted abundance is given by which is valid when the quantity inside the square brackets is less than unity. Lastly, the thin vertical curves label the axino FI domination. Since axino FI production in eqs. (3.7) and (3.9) is independent of mã, the result generically applies for any mã between 650 GeV− 1 TeV as required by the free-streaming length for cold dark matter and neutralino decay kinematics respectively. In the lower panels, T R < 10 10 GeV and the dilution factor in the cosmology described in section 2.2 scales linearly with T R . Therefore, the curves for equilibrium gravitino production (with a constant undiluted yield) are now at a slope. On the other hand, nonequilibrium gravitino UV production is close to being proportional to T R as well, and with dilution, the final yield becomes independent of T R , corresponding to the vertical thick parts. Lastly, the axino FI production dominates at the vertical thin parts of the curves.
In figure 14, the LOSP lifetimes can again be predicted from eq. (4.6) with the values of V PQ that give the observed dark matter abundance. The left (right) panel corresponds to the cosmology of T R less (greater) than 10 10 GeV discussed in section 2.2 (section 2.1). It is worth noting that we are choosing different values of the gravitino mass in these JHEP03(2017)005 two panels to demonstrate interesting LOSP decay signals for the relevant regions. In particular, we open up a new decay channel when the gravitino mass is less than O(1) MeV -the neutralino LOSP decays to the gravitino and γ/h/Z. This decay channel, with the rate given in eq. (4.7), dominates in the purple portions of the curves and gives the very well known displaced vertex signal of gauge mediation. Our framework yields a cosmology for this scenario with T R far above the TeV scale. As in figure 13, the thin parts of the lines are dominated by the axino FI contribution and lead to warm dark matter because of the late axino decay to axions and gravitinos. For clarity, only one value of s I is shown but there exists a set of correlated m 3/2 and s I that can lead to LOSP displaced signals. In particular, a low m 3/2 enhances the production and thus requires a larger s I for dilution. However, once m 3/2 is sufficiently low for the gravitino to be thermalized, any lower m 3/2 results in a decrease in its energy density and a smaller s I is needed. It is remarkable that, due to the interplay between the two decay channels, the LOSP lifetimes for low m 3/2 are always within the reach of current and future colliders.

The DFSZ + theory
The dark matter abundance is shown in figure 15 as contours of Ωh 2 = 0.11. Since the thick parts of the curves indicate gravitino UV domination, it can be seen that axino production never dominates when T R 10 5−6 GeV. Thus, the features of these parts are identical to those in figure 13. To see gravitino UV domination, we first notice from figure 11 that axino UV production generally dominates that of axino FI. Above black dashed lines, UV axinos thermalize and Y eq a /Y eq G = 1/2 based on eq. (3.12). In comparing the non-equilibrium axino UV production to that of the gravitino In the lower panels, since the dilution factor becomes unity in regions where T R 10 5−6 GeV and V PQ 10 14−15 GeV, the gravitino UV abundance decreases with T R and the curve becomes thin and vertical, corresponding to axino FI. As axino FI and UV production in eqs. (3.7) and (3.9) and eq. (3.11) are independent of mã, the result generically applies for any mã between 650 GeV − 1 TeV as required by the free-streaming length for cold dark matter and neutralino decay kinematics, respectively.
Finally, we also make predictions of neutralino LOSP lifetimes in figure 16. Similar to figure 14, the relevant decay modes are N 1 →ã (red) and N 1 →G (purple), with the latter dominating for low m 3/2 . The prediction of the total decay rate, when dominated by Γ N 1 →G , depends on m 3/2 and is insensitive to dark matter production and dilution.
On the other hand, once m 3/2 is sufficiently large, N 1 →ã dominates and the lifetime prediction is set by the value of V PQ that gives Ωh 2 = 0.11 in figure 15. The thin (thick) lines label warm (cold) dark matter. The left (right) panel corresponds to the low (high) T R cosmology studied in section 2.

A Axion supermultiplet interactions
In this appendix we develop a general framework to describe effective interactions of the axion supermultiplet. We consider both self-interactions and couplings to MSSM fields. In the next appendix, we use these results to compute the decay widths used in this work. At energies below the PQ breaking scale the theory is still approximately supersymmetric. We assume the PQ symmetry to be broken by a set of chiral superfields Φ i , which we expand around their vacuum expectation values (vev) The axion a fills the supermultiplet A as explicitly given in eq. (1.1), together with its superpartners, the saxion s and the axinoã. We normalize the PQ charges q i of the PQ breaking fields in such a way that they are all integers and |q i | as small as possible. Within

JHEP03(2017)005
this convention, the effective PQ breaking scale is defined as follows The effect of a PQ rotation with an angle α on the axion superfield is the following The interactions of the axion supermultiplet are significantly constrained by this shift invariance.

A.1 Color anomaly
The PQ symmetry is broken by a color anomaly, which generates the low-energy interaction between the axion superfield A and the supersymmetric QCD field strength W α . The domain wall number N DW appearing in the above expression is the QCD anomaly coefficient of the PQ symmetry. Finally, we define the axion decay constant Upon expanding in component fields the supersymmetric expression in eq. (A.4), we identify the effective interaction between the axion a and the QCD field strength G

A.2 Supersymmetric interactions
We assume the fields in eq. (A.1) to be canonically normalized, and the Kähler potential reads This function only depends on the PQ invariant combination A + A † , consistently with the invariance under the shift in eq. (A.3). The axion superfield A is canonically normalized. Holomorphy and PQ invariance forbid superpotential self-interactions for A. However, the axion superfield A can appear in the superpotential of DFSZ theories, where the µ term is PQ charged. We are allowed to write the PQ invariant operator Here, we denote q µ as the model-dependent PQ charge of the µ term.

A.3 SUSY breaking interactions
Finally, we account for SUSY breaking. We find it convenient to employ the non-linear field where F is the SUSY breaking scale andη is the associated goldstino. SUSY breaking is transmitted to the PQ sector through the higher dimensional operator with M mess the mass scale of the particles coupling the two sectors. We expect this operator to be generated by Planck scale dynamics [17], and therefore we have the upper limit on the mediation scale M mess M Pl . A consequence of this operator is a contribution to the axino mass where we have used the known relation m 3/2 = F/( √ 3M Pl ). As discussed in the main text of this work, it is not natural to have an axino much lighter than the gravitino, unless there is some good reason to suppress the size of the coefficient c AAX (e.g. sequestering).
Furthermore, we also require the presence of an effective Bµ term necessary for a successful electroweak symmetry breaking with c B F = Bµ.

B Decay widths
Saxion decays are responsible for reheating the universe and producing dark radiation. Decays and inverse decays of neutralinos and charginos generate an axino freeze-in abundance. The axino decay to the gravitino creates a warm DM population and vice versa. Finally, neutralino decays to the axino or gravitino lead to displaced collider events. We compute the decay widths for all these processes, using the interactions derived above.

B.1 Saxion decays
The saxion can decay to three possible final states: axions, axinos and Higgs bosons. The first two processes are mediated by the axion multiplet self-interactions, more specifically the cubic Kähler potential term in eq. (A.7). Once we expand in component fields, we find the operators

JHEP03(2017)005
Here, we define the dimensionless parameter κ = i q 3 i v 2 i /V 2 PQ . For models with only a single PQ breaking field, or theories with more than one but all with the same PQ charge, we have κ = 1. In more general cases κ is a free parameter. If κ is non-zero, axinos can be copiously produced from saxion decays. Neglecting the final state masses, the decay widths are More importantly for the reheating process, the saxion can decay to Higgs bosons and longitudinal electroweak gauge bosons. These decays are mediated from the supersymmetric interactions arising from the superpotential in eq. (A.8) as well as the SUSY breaking interactions arising from eq. (A.12). The resulting scalar potential interactions are The Higgs doublets H u and H d contain three Goldstones (G ± and G 0 ) providing electroweak gauge bosons longitudinal modes, two CP-even (h and H) and one CP-odd (A) neutral scalars and one charged scalar (H ± ). We consider the decoupling limit where the non-SM fields are heavy. In such a limit, we find it convenient to introduce the doublets The SM-like Higgs boson h lives within the multiplet H SM , which takes the electroweak symmetry breaking vev. The decoupling limit holds as long as m A m Z , and in such a limit the connection between gauge eigenstates and mass eigenstates reads where we introduce the ratio between the two vevs tan β = v u /v d and with H (c) i = iσ 2 H * i . We report here the decay width in such a decoupling limit and for large tan β (for more general expressions see appendix A of ref. [26]) (B.8) Here, D counts the number of kinematically allowed final state particles. For a SM Higgs sector we have D = 4, whereas D = 8 for the full MSSM Higgs sector.

B.2 Neutralino and chargino decays to axinos
We start by defining the mixing matrices for the MSSM neutralinos and charginos. We collect the neutralinos into the four-dimensional array and the connection with mass eigenstates N i is achieved through the rotation Likewise, we group the charginos into the two different two-dimensional arrays containing positively and negatively charged states, respectively. We separately rotate the two chargino arrays and as expected we find that the mass eigenvalues for the two arrays are the same. The positively and negatively charged mass eigenstates fill a Dirac fermion. For neutralinos decays and inverse decays we have |s β R 3i | 2 + |c β R 4i | 2 , (B.14) where R is the rotation matrix defined in eq. (B.10). Likewise, for charginos with U and V defined in eqs. (B.12) and (B.13).

JHEP03(2017)005
The supercurrent J µ contains both chiral superfields with components (φ, χ) as well as gauge superfields with components (V, λ). Here, F is the field strength of the vector field V . The gravitino decay width to chiral multiplet components reads (B.21) Here, the multiplicity factor is (8, 3, 1) for the final state (gluino, wino, bino).

B.4 Axino and neutralino decays to gravitinos
In the last part of this appendix, we consider decays to final states involving gravitino. We assume the gravitino to be much lighter than the decaying particle, so that we can approximate the process with decays to longitudinal gravitinos, in accordance with the equivalence theorem. The process is correctly described by the effective Lagrangian 5 where we consider again both chiral (φ, χ) and vector (V, λ) supermultiplets. We note that the above interactions can also be derived from the full supergravity result in eq. (B.18), by identifying the longitudinal component of the gravitinõ Upon using the Goldstino equations of motion [62,63] and the relation m 3/2 = F/( √ 3M Pl ), we recover the interaction between the Goldstino and the flat-space global SUSY supercurrent.
The first case we discuss is the axino decay to axion and gravitino. The associated axino decay to saxion and gravitino is assumed to be kinematically forbidden. Upon using the general supercurrent result, and accounting for only the axion final state, we find The last cases we discuss are the ones relevant for displaced collider signatures. We only consider decays of the lightest neutralino N 1 , since all other R-odd particles produced JHEP03(2017)005 at collider will promptly decay to it. The mass eigenstate N 1 is related to the gauge eigenstates through the rotation in eq. (B.10). For decays to photons we are only sensitive to Goldstino interactions with the neutral gauginos, and we find Γ N 1 →G γ = |R 11 c w + R 21 s w | 2 C Free streaming of warm DM component Our framework provides a warm dark matter source through production of NLSP particles and subsequent decay NLSP → LSP a, where a is a (nearly) massless axion. We always assume that if the LSP is the gravitino (axino), then the NLSP is the axino (gravitino).
In this appendix, we provide the calculation of the free streaming length for such a warm component a(t) dt = 2 t eq a 2 eq aeq aτ v(a) da . (C.1) Here, τ NLSP is the NLSP lifetime, whereas t eq and a eq are the time and scale factor at matter-radiation equality, respectively. In the second equality, we defined a τ = a(τ NLSP ) and we changed integration variable by using the relation for a radiation dominated universe a ∝ t 1/2 , justified if NLSP decays happen after BBN (i.e. τ NLSP 1 sec). The LSP velocity v(a) after decays is just a consequence of free streaming. The initial energy and momentum at the decay time τ NLSP follow from four-momentum conservation It is convenient to derive an approximate expression for λ FS by identifying the scale factor value a NR , correspondent to the time when the free streaming LSP enters the nonrelativistic regime. The LSP velocity in eq. (C.3) can be approximated as follows v(a) 1 a < a NR a NR /a a ≥ a NR . (C.6) We find a NR by imposing p LSP (a NR ) m LSP and we find a NR a τ p τ LSP /m LSP . The free streaming length, as defined in eq. (C.1), approximately reads λ FS 2 t eq a eq a τ a eq p τ LSP m LSP 1 + log a eq a τ m LSP p τ LSP . (C.7) We evaluate this expression by using the known values t eq /a eq 93 Mpc, t eq 2.2 × 10 12 sec, and the time dependence of the scale factor a(t) = a eq (t/t eq ) 1/2 . In the m NLSP m LSP limit, the free streaming length reads