Saxion Cosmology for Thermalized Gravitino Dark Matter

In all supersymmetric theories, gravitinos, with mass suppressed by the Planck scale, are an obvious candidate for dark matter; but if gravitinos ever reached thermal equilibrium, such dark matter is apparently either too abundant or too hot, and is excluded. However, in theories with an axion, a saxion condensate is generated during an early era of cosmological history and its late decay dilutes dark matter. We show that such dilution allows previously thermalized gravitinos to account for the observed dark matter over very wide ranges of gravitino mass, keV<$m_{3/2}$<TeV, axion decay constant, $10^9$ GeV<$f_a$<$10^{16}$ GeV, and saxion mass, 10 MeV<$m_s$<100 TeV. Constraints on this parameter space are studied from BBN, supersymmetry breaking, gravitino and axino production from freeze-in and saxion decay, and from axion production from both misalignment and parametric resonance mechanisms. Large allowed regions of $(m_{3/2}, f_a, m_s)$ remain, but differ for DFSZ and KSVZ theories. Superpartner production at colliders may lead to events with displaced vertices and kinks, and may contain saxions decaying to $(WW,ZZ,hh), gg, \gamma \gamma$ or a pair of Standard Model fermions. Freeze-in may lead to a sub-dominant warm component of gravitino dark matter, and saxion decay to axions may lead to dark radiation.


Introduction
If supersymmetry is relevant for the hierarchy problem, gravitinos, with a mass suppressed by the Planck mass, become an interesting candidate for dark matter, as pointed out by Witten [1]. However, the cosmology of gravitinos has long been viewed as problematic. In 1981 Pagels and Primack found that light gravitinos would overclose the universe if they were heavier than the keV scale [2]. To obtain the dark matter abundance revealed by recent measurements, the gravitino mass must be around 100 eV, which is excluded due to the warmness of the gravitino [3]. These pioneering works assumed that gravitinos, like photons and neutrinos, would be in thermal equilibrium in the very early universe with a high temperature, so that their number density would be given by thermodynamics. Since then it has typically been assumed that, in theories with weak scale supersymmetry, the reheat temperature of the universe after inflation T R is severely restricted [4], to strongly limit the gravitino abundance.
However, in supersymmetric theories with a Peccei-Quinn (PQ) symmetry [5,6] broken at scale V PQ to solve the strong CP problem, the gravitino abundance can be diluted by the late decay of a saxion condensate [7,8] which is generated by supersymmetry breaking during an early era, for example during inflation [9]. (See [10,11] for dilution by thermally produced saxions.) Hence, in this paper we return to the original assumption of Witten, Pagels and Primack that gravitinos were in thermal equilibrium in the very early -1 -JHEP07(2017)125 universe. We take the gravitino to be the Lightest Supersymmetric Partner (LSP) and study constraints on such gravitino dark matter, exploring which regions of the (m 3/2 , V PQ ) plane are preferred and hence relating the scales of supersymmetry and PQ breaking. Most constraints are independent of the very early cosmological history of the saxion oscillations, depending only on the saxion evolution at temperatures less than or of order of the masses of superpartners, especially those of the higgsino and lightest observable supersymmetric particle (LOSP).
Dilution from the saxion condensate allows a much wider set of cosmologies, in particular allowing T R to be arbitrarily high. For comparison, without dilution a light gravitino, m 3/2 < MeV, requires T R below the TeV scale of superpartners. Thus saxion dilution allows high T R scenaria for the interesting case of displaced vertex signals at LHC, and its MATHUSLA extension [15], arising from decays to gravitinos.
Cosmological axion production from the misalignment mechanism [12][13][14] is frequently taken to limit V PQ = f a N DW / √ 2 < ∼ N DW × 10 12 GeV, where f a is the axion decay constant and N DW is the domain wall number. However, with dilution from a saxion condensate this limit is weakened by 3-4 orders of magnitude [16][17][18]. Hence we also explore the abundance of axion dark matter, finding regions of parameter space in both DFSZ [19,20] and KSVZ [21,22] models where it can be a significant component of dark matter.

The cosmological history
In this section we provide an overview of the cosmological evolution of the saxion condensate and the thermal bath, and we give results for the axion abundance.
In the absence of supersymmetry breaking, the saxion field s has no potential. In the early universe, at any era the non-zero energy density of the universe breaks supersymmetry, and hence the form of the saxion potential is a highly model-dependent question. As the universe evolves through inflation, post-inflation, reheating and subsequent eras the saxion potential and its minimum changes leading, in general, to a highly complicated evolution of the saxion condensate. Rather than studying a particular model, we show that the physics relevant for gravitino dark matter depends on the saxion evolution only at temperatures less than or of orderm, the masses of the SM superpartners, specifically the masses of the higgsino and the LOSP, which we take to be O(TeV). Thus, to obtain the main results of this paper the assumption on the cosmological history of the universe at temperatures above the TeV scale is extremely mild: • Before reaching the TeV scale, the saxion field acquired a large displacement from its present value and there was an era where gravitinos were in thermal equilibrium.
Furthermore, for later evolution of the saxion field we assume • From T ∼m ∼ TeV until it decays, the saxion condensate oscillates about the present vacuum in a quadratic potential where m s is the soft supersymmetry breaking mass of the saxion.

JHEP07(2017)125
This large saxion condensate plays a crucial role in determining the dark matter abundance. It eventually dominates the energy density of the universe and releases most of its entropy, when Γ s ∼ H and the saxion decays become significant, at temperature The decay rate of the saxion, Γ s , is dependent on the saxion mass and on whether the Higgs doublets carry PQ charge.
In DFSZ models with m s > 2m W , the saxion mainly decays into a pair of Higgs, W or Z bosons, with the rate where we have summed over the final states and assumed the decoupling and large tanβ limits. The PQ charge q µ of the Higgs mass parameter µ is normalized such that all charges of the PQ breaking fields are integers with absolute values as small as possible. We fix q µ = 2 in this paper, as in the minimal supersymmetric DFSZ model. For a lighter DFSZ saxion, m s < 2m W , the main decay channel is into a pair of standard model fermions via mixing with the Higgs, with the rate where N f is the multiplicity of the fermion f (3 and 1 for quarks and leptons, respectively).
Here we have assumed the decoupling and large tanβ limits as well as m s m h . In KSVZ models, the saxion mainly decays into a pair of gluons with a rate (2.5) Here the axion decay constant f a is defined by the axion coupling with the gluon field as 6) and is given by the PQ breaking scale V PQ through the relation f a = √ 2V PQ /N DW . In either case, T Rs is low enough that the gravitinos thermalized at early times are diluted. For most values of V PQ , the decay of the saxion condensate also dilutes axinos and gravitinos from freeze-in (FI) and superpartners from freeze-out. The FI mechanism refers to the process where dark matter (DM) is not in the thermal bath but produced from the decay or scattering of the thermalized particles. For higher values of V PQ , even misalignment axions generated near the QCD phase transition are diluted.
In the next sub-section we discuss aspects of the saxion oscillation matter-dominated era. In the following sub-section we provide a very simple illustration of a possible cosmology for the early evolution of the saxion condensate -the "Decoupled Saxion". In this paper the processes relevant for computing the dark matter abundance are dilution of previously thermalized gravitinos, freeze-in gravitino production and misalignment axions. These processes all occur at low temperatures, T < ∼m . Furthermore, constraints on the theory from overproduction of freeze-in axinos and superpartner freezeout similarly occur at T < ∼m . Hence the results of this paper only depend on the low temperature aspects of the cosmology, not the high temperature aspects. Thus the evolution of the saxion at T >m could be arbitrarily complicated, for example from interactions during inflation or with the thermal bath. Nevertheless, for the late evolution at T < ∼m in the potential (2.1) we need to parameterize the size of the condensate.
In particular well before it decays, the saxion condensate must dominate the energy density producing a matter dominated universe at T T Rs . From T Rs up to some temperature T NA this MD era is non-adiabatic (MD NA ): the radiation density is dominated by the products of recently decayed saxions rather than from pre-existing red-shifted radiation, giving T ∝ 1/a 3/8 in ref. [23]. On the other hand at temperatures above T NA there are so few saxion decays that the MD era is adiabatic (MD A ), with T ∝ 1/a. At T NA the saxion condensate has a size s NA T 4 NA /T 2 Rs m s and we find it convenient to use T NA to describe the strength of the saxion condensate, as it appears directly in the gravitino dilution factor. For decoupled relic particles produced at temperatures above T NA , such as the previously thermalized gravitinos, saxion decays yield a dilution factor 1 as derived in ref. [23], while for relic particles produced at some temperature T between T NA and T Rs the dilution factor is less: (2.8) The condition that fixes T NA follows from requiring that the dilution of previously thermalized gravitinos yields the observed temperature of matter-radiation equality We study gravitino dilution over a very wide range of parameters: V PQ is varied over its entire range from its lower astrophysical bound (see [24] for a review) of 10 9 GeV to M Pl , and the saxion mass is varied over the range of 10 MeV < m s < 10 TeV. Throughout this parameter space, the observed dark matter results from diluting previously thermalized gravitinos and/or gravitinos produced by freeze-in processes, and T NA is constrained to be in the range 10 MeV < T NA < 100 TeV.
In DFSZ theories with large µ and small V PQ , the formulae (2.2) and (2.9) give T NA , T Rs > m s . For temperatures above m s , however, the decay/scattering of the saxion is affected by thermal effects [25][26][27][28], which determine the temperatures T Rs and T NA .

JHEP07(2017)125
For example, when the saxion has a Yukawa interaction ysff with a fermion f , the decay (dissipation) rate of the saxion is given by Γ ∼ y 2 T for T m s . The resultant reheating temperature is given by T Rs ∼ y 2 M Pl , and the dilution factor is D ∼ (T NA /T Rs ) 3 . In the next section we show the constraint on the parameter space in figures 2-4. In the lower part of figure 2, thermal effects determine T Rs and/or T NA . For simplicity we do not show the contours of required T NA if T NA > m s . In figure 4 this thermal effect is irrelevant.

Cosmology of the "Decoupled Saxion"
As a particular example of a saxion cosmology at T >m we consider the "Decoupled Saxion", defined by the assumption that the saxion potential is given by eq. (2.10) This could happen if the saxion couples to the thermal bath via either very small dimensionless couplings or through suppressed higher dimension operators. We stress that this is just a simple illustrative example, and is not necessary for the results of the next section.
Taking the reheat temperature after inflation, T R , to be larger than T osc , the saxion field starts to oscillate at T osc with some large amplitude s I that it acquired from some previous era, so that during the adiabatic era following T osc the saxion energy density is (2.11) The universe becomes matter dominated by the saxion condensate at and the subsequent matter-dominated era becomes non-adiabatic at T NA , when the radiation bath becomes dominated by saxion decay products rather than by the red-shifted radiation from inflaton decay, with Hence, in this scenario it is best to describe the strength of the saxion condensate by s I and have T NA as a derived quantity given by (2.13). Furthermore, in this cosmology, D of eq. (2.7) becomes D = T M /T Rs . The gravitino dark matter abundance constraint of eq. (2.9) then leads to (2.14) which we will find yields very large values of s I in the range of (10 14 − 10 18 )GeV and is correlated with other parameters according to   for saxions decaying to pairs of (gluons, electroweak and Higgs bosons) in the (top, bottom) row respectively.
The required large field value may seem to be incompatible with the assumption of the quadratic potential. Actually in models where the saxion mass is generated via quantum corrections, the saxion potential in general becomes logarithmic at a large field value, and hence the assumption breaks down. However, in models where the saxion mass is given by a tree-level superpotential [29][30][31], the saxion potential may be quadratic even for a large field value.

Misalignment axion contribution to dark matter
The axion is produced by the usual misalignment mechanism but now we need to consider the effect of saxion dilution [16][17][18]. The axion field value is initially displaced from the minimum today by an amount f a θ i , where θ i is the misalignment angle. Coherent oscillations of the axion field commence at temperature T (a) osc , when the Hubble rate is equal to the axion mass. We assume that the axion oscillation starts during the MD NA era, which is the case for the parameter spaces constrained by the axion abundance. In this case, the calculation is independent of T NA . We evaluate the axion energy density per entropy density at the end of entropy production [32] ρ a s T Rs = 9 8 osc ) takes into account the temperature dependence of the axion mass. We assume a simple power law m 2 a (T ) = m 2 a (0)(Λ/T ) γ above the QCD scale Λ. The mass takes a constant value, m a (0) = 6 eV (10 6 GeV/f a ), at T < Λ. In other words, if T (a) osc ≤ Λ, ξ = 1; otherwise, ξ = (Λ/T (a) osc ) γ/2 is used to compute the axion energy density today. The analytic formula of ξ was derived in ref. [32].
We predict the axion abundance in terms of T Rs and f a = √ 2V PQ /N DW . We use N DW = 1, while T Rs in eq. (2.2) can be calculated for both DFSZ and KSVZ theories. Figure 1 shows the numerical result of the contours of Ω a h 2 = 0.11 for various misalignment angles θ i and the axion mass index γ obtained from lattice calculations. The region above and to the right of the contour is excluded by axion overproduction. The top (bottom) axis refers to the set of parameters necessary to compute T Rs in DFSZ (KSVZ) theories. For DFSZ, we assume m s in a range where the saxion decays dominantly to W + W − , ZZ and hh. The solid (dashed) lines are for the axion mass index γ = 6.8 (2.7) computed in ref. [33]( [34]). In the regions where the two index lines merge, the axion starts oscillating after its mass is already a constant, i.e. T (a) osc < Λ. Recent lattice calculations [35][36][37][38] show that the axion mass index is well described by the dilute instanton gas approximation, γ 8, in high temperature regimes. The red region is excluded by BBN because of the late decays of the saxion.
Interestingly, including the dependence of ξ on V PQ , one finds that Ω a h 2 decreases (increases) for γ > 4 (γ < 4) when V PQ increases. This explains the different overall slopes of the two index lines, while the detailed features arise due to the rapid change of -6 -JHEP07(2017)125 g * (T ) during the QCD phase transition. In the case of γ = 6.8, the dependence on V PQ and effects of g * (T ) compete with each other, resulting in a nearly vertical contour. We compute g * (T ) from the Boltzmann distribution with the full SM spectrum and the MSSM spectrum degenerate at 1 TeV, and we linearly interpolate g * (T ) across the QCD phase transition, i.e. 100 MeV < T < 300 MeV.
The axion abundance gives an upper bound on the higgsino mass in the DFSZ model. The bound is stringent for m s > 2m W . Note first that the saxion mass cannot be larger than 2µ, since otherwise saxions decay into pairs of higgsinos and result in too large gravitino abundance. Then from the upper x-axis of figure 1, we can derive an upper bound on 2 −1/3 µ 0.8µ. Assuming the fine-tuning in the misalignment angle is no more than 10%, the higgsino mass should be smaller than about 1 TeV for V PQ > ∼ 10 13 GeV. If the misalignment angle takes a randomized value, the axion abundance should be evaluated with the averaged angle, θ mis π/ √ 3. Then the higgsino mass should be smaller than about 200 GeV for V PQ > ∼ 10 12 GeV. This includes the cases where the PQ symmetry is unbroken during inflation, restored after inflation, or the axion field obtains large fluctuation due to the parametric resonance effect from saxion oscillations. The last case is discussed in section 3.4. For m s < 2m W , the upper bound is relaxed. Since a fermion with a mass close to m s dominantly contributes to the decay rate in eq. (2.4), the bound is relaxed roughly by a factor of (100 GeV/m s ) 3/4 in comparison to the case with m s > 2m W .

Thermalized gravitino dark matter
In this section we show that dilution by the late decay of a saxion condensate allows the observed dark matter abundance to arise from the thermalized gravitinos of an early epoch over a very wide range of m 3/2 and V PQ . However, gravitinos can be overproduced by reactions occurring at the TeV scale or below: gravitino freeze-in, axino freeze-in and decay to gravitinos, and saxion decays toã +G. We illustrate how these constrain the region where dark matter arises from the primordially thermalized gravitinos. Similarly, we indicate where axions are overproduced, by either early misalignment or parametric resonance during saxion oscillations.
There are several relevant parameters. We show results for essentially complete ranges of (m 3/2 , V PQ ), but choose a few illustrative values for the key parameters (m s , mã, µ) and for other supersymmetry breaking parameters. Our aim is not to provide an exhaustive study of the (m s , mã, µ) space, but to illustrate the wide range that allows thermalized gravitino dark matter and its corresponding rich signals. We examine constraints on the parameter space from other processes creating gravitinos and axions in sub-section B (C) for dominant saxion decays to gluons (Higgs/electroweak bosons).
We comment on the lower bound on the saxion mass. In KSVZ theories, in order for the saxion to decay before the BBN, is required. In DFSZ theories the saxion mass may be smaller due to its effective couplings with standard model fermions through its mixing with the Higgs. However, with such a mixing, the saxion takes away energy from supernovae and changes the duration of the neutrino emission [39][40][41][42]. To prevent this process requires

The maximal parameter space
The unshaded regions of figure 2 show the maximum ranges of (m 3/2 , V PQ ) that allow for thermalized gravitino dark matter from saxion condensate dilution. The left/right panel is for m s = 30 GeV/ 3 TeV and each panel applies to both KSVZ and DFSZ theories. Values of m 3/2 below ≈ keV are excluded by warm dark matter constraints, while values above 100 GeV are possible as long as the gravitino remains the LSP. The gray shaded region of figures 2, 3 and 4 is excluded because the contribution to the vacuum energy from the saxion potential, of order m 2 s V 2 PQ , exceeds that allowed by total supersymmetry breaking F 2 tot = 3 m 2 3/2 M 2 Pl . This bound on V PQ scales as m 3/2 /m s . The bound is saturated in models where supersymmetry and PQ symmetry are simultaneously broken, and the saxion obtains its mass at tree level [29][30][31]. In models where the saxion mass is given by quantum corrections, the bound is stronger by coupling constants and associated loop factors.   The red shaded region of figures 2, 3 and 4 is excluded because the reheat temperature from saxion decays, T Rs , is below 3 MeV destroying the success of BBN [43]. In KSVZ theories, where the dominant saxion decay is s → gg, this bound on V PQ scales as m 3/2 s . In both panels of figure 2, the red shading applies only to KSVZ theories. For the right panel in DFSZ theories, the dominant saxion decay is s → W W, ZZ, hh giving a bound on V PQ that scales as µ 2 /m 1/2 s . We have taken µ = 2 TeV so that this bound on V PQ is larger than 10 16 GeV and does not appear in the figure. In the left panel for DFSZ theories, the dominant saxion decay is s →bb giving a bound on V PQ that scales as µ 2 m 1/2 s and again is larger than 10 16 GeV and does not appear.
The saxion decays into a pair of gravitinos through its mixing with the sgoldstino field or the mixing of the axino with the gravitino. The decay rate of the saxion into a pair of gravitinos is Here κ is an O(1) parameter which depends on the couplings between the PQ breaking field and the supersymmetry breaking field, and may be suppressed if there is an (approximate) Z 2 symmetry in the couplings. Even if κ = O(1), we found that this decay mode does not give additional constraints beyond the gray and red shaded regions. These bounds, leading to the gray and red excluded regions, are inherent to the saxion condensate dilution mechanism and cannot be evaded. Other processes producing gravitinos or axions are frequently important, and may lead to further constraints in the -9 -

JHEP07(2017)125
(m 3/2 , V PQ ) plane, but they depend on other parameters, such as the axino mass, or on details of cosmological evolution at temperatures far above the TeV scale. Hence we omit them from figure 2, which shows the maximal allowed region, and consider them at length in the next two sub-sections and in figure 3 and figure 4. Here we provide a brief qualitative illustration of how these other constraints can be avoided.
Contributions to gravitino dark matter from s →Gã,ãã are avoided by taking mã > m s . A sufficiently large mã also removes constraints from axino freeze-in. Effects on BBN arising from the lightest observable supersymmetric particle (LOSP) freezeout and decay can be made sufficiently mild by having a sneutrino LOSP. With a sneutrino mass of 300 GeV the freezeout abundance is quite small; and neutrinos from decay to νG have only mild effects on BBN [44]. A stau LOSP also makes the BBN constraint mild; a gravitino mass below 10 GeV is allowed. The gravitino freeze-in abundance is controlled by the LOSP mass, and 300 GeV is already large enough to provide a sub-dominant contribution to dark matter. A crucial feature of our scheme is that freeze-in of both axinos and gravitinos are highly suppressed as they occur during a matter dominated era and are subsequently diluted by saxion decays. We also note that in KSVZ theories the decay s → aa must be mildly suppressed for s → gg to dominate. Finally, there is the possibility that during the oscillation of the saxion field inhomogeneities in the axion field are exponentially enhanced by parametric resonance. However, the importance of this effect depends on the very early cosmological evolution of the saxion field, and is model dependent. While all these constraints can be avoided, they are frequently important and we discuss them quantitatively below in sections 3.2 and 3.3.
We find that previously thermalized gravitinos decouple from the thermal bath at a temperature higher than T NA . This means a large amount of entropy is injected only after these gravitinos stop interacting with the bath. As a result, the gravitino abundance is diluted by the factor D of eq. (2.7). Requiring dilution to yield the observed dark matter abundance via eq. (2.9), gives an analytic estimate for T NA . In KSVZ theories saxions decay dominantly to gluons, with a rate given in eq.   Using the full Boltzmann equations for the cosmological evolution and dark matter dilution, we show the numerical results for T NA in figure 2 in the (m 3/2 , V PQ ) plane as solid/dashed contours for saxion decays to gg/(W W, ZZ, hh). Throughout the entire allowed region T NA < O (TeV) so that dilution of previously thermalized gravitinos by decay of the saxion condensate involves cosmology of the TeV era or later. Hence we are able to discuss this scenario in a very general framework, without the need to specify a particular UV theory, by making the two key assumptions listed at the beginning of section II.

Further constraints in KSVZ theories with s → gg
We explore further constraints on KSVZ theories from freeze-in of gravitinos and freezeout of the LOSP. First we perform a similar calculation of T NA contours as in figure 2, with numerical results given in figure 3 for m s = (1, 10) TeV in the (left, right) panel. The allowed parameter space is the white (unshaded) region. The red and gray regions excluded by BBN and the consistency of supersymmetry breaking are discussed in section 3.1. The dot-dashed (dashed) lines of figure 3 and figure 4 give contours of s I = V PQ (10 V PQ ), using eq. (2.13). s I is the initial saxion field value in the "Decoupled Saxion" cosmology of section 2.2 and the importance of these contours for axion parametric resonance is discussed in section 3.4.
In addition to gravitinos thermalized during an early epoch, the decays of supersymmetric particles to gravitinos also contributes to the final abundance via the FI mechanism. The freeze-in process is IR dominated and terminates when the abundance of supersymmetric partners becomes exponentially suppressed as the temperature falls below their masses. Since gravitinos interact with all multiplets via the goldstino interaction, this FI -11 -

JHEP07(2017)125
contribution is determined mainly by the LOSP and hence the LOSP mass. For illustration purposes, we assume a higgsino LOSP.
The FI abundance is proportional to the decay rate, which is enhanced for low gravitino masses and high parent particle masses. As shown in figure 3, the FI contribution is absent for µ = 1 TeV because the higgsino LOSP decay is inefficient. However, with µ = 6 TeV in the right panel, freeze-in gravitinos become the dominant source in the brown shaded region. Furthermore, since the freeze-in of gravitinos occurs below T NA in this region, the dilution factor is given by eq. (2.8) and depends only on µ and T Rs but not T NA . As a result, the brown region is excluded because the saxion decays too early to provide sufficient dilution for FI gravitinos, regardless of T NA .
The regions shaded in magenta in figure 3 are excluded by BBN due to late decays of higgsinosH produced via freeze-out. The subsequent decays of higgsinos place constraints on high gravitino masses, where the decay is late, and on low V PQ , where the saxions decay early. The higgsino mass, µ, is increased from the left to the right panel, increases both the higgsino decay rate and the higgsino FO abundance. The BBN constraint shifts to the right and upward.

Further constraints in DFSZ theories with s → hh/ZZ/W W
The analysis for KSVZ and DFSZ theories has both similarities and differences, as seen by comparing figures 3 and 4. In figure 4 we fix m s = 300 GeV, while µ and mã take on different values in the four panels, such that in the upper (lower) panels, mã > m s (mã < m s ). Numerical solutions for contours of T NA are given in the (m 3/2 , V PQ ) plane. The white (unshaded) regions give the allowed parameter space, which is significantly limited by several constraints as discussed below.
While the gray shaded regions requiring m 3/2 V PQ < √ 3F tot are the same as in KSVZ theories, the red shaded region from BBN limits on saxion decay are much milder, especially for m s > 2m W . The saxion decay rate to gluons, eq. (2.5), is loop-suppressed relative to the decay rate to electroweak bosons, eq. (2.3). In DFSZ theories, this faster saxion decay rate increases T Rs so that the region with T Rs < 3 MeV excludes only the largest values of V PQ .
As in the KSVZ case, DFSZ theories can have small regions at low V PQ excluded by over-closure from gravitino freeze-in (brown region at low m 3/2 ) and from BBN limits from the decays of freeze-out higgsinos,H → hG, (magenta region at large m 3/2 ). The former disappears in the left panels because the higgsino decay rate to gravitinos depends strongly on µ. The latter disappears in the lower panels because the higgsinos from FO decay before BBN to axinos which then decay harmlessly to aG.
A key feature of DFSZ theories is additional production processes involving axinos, if mã is not too large, giving the large green, purple and orange regions of figure 4.
Axino FI and decay to gravitinos yields too much dark matter in the green regions at low V PQ . The freeze-in of axinos (and gravitinos) from decays of the higgsino LOSP occurs during the MD NA era or after saxion reheating for low enough V PQ and m 3/2 . This implies that the relevant dilution factor is (2.8) and is insensitive to T NA and solely determined by the freeze-in temperature and T Rs .     Axino FI and decay is excluded by BBN in the purple regions of figure 4. In the top left panel the relevant decay chain isã →Hh followed byH → hG. In the top right panel the higgsino is heavier than the axino so the relevant decay chain isã → sG followed by s → hh/ZZ/W W . In the lower two panels the axino is the NLSP and the only decay mode isã → aG, which is harmless. Note that, in both purple and magenta regions, the exclusion from BBN is coming from the long lifetime of the decay to gravitinos at large m 3/2 .
In the lower two panels, the axino mass is sufficiently small that a new saxion decay channel opens up: s →ãG. This leads to excessive dark matter in the large orange regions, -13 -

JHEP07(2017)125
again showing the critical importance of the axino mass in DFSZ theories. Note there is a complementarity between s →ãG (for mã < m s ) leading to excessive dark matter and a → sG (for m s > mã) leading to BBN problems.
We have not displayed a panel for the spectrum m s > mã > µ, where saxions decay into axinos (giving an orange over-closure region), and axinos decay into higgsinos (giving a purple BBN region). In this case, the excluded regions are roughly the unions of these regions from the upper and lower panels of figure 4. Since we require m s < 2µ, this applies to a small range of m s .

A limit on the saxion condensate from axion parametric resonance
The constraints derived in the previous sub-section require only the two mild assumptions declared early in section 2. To discuss the constraint from axion parametric resonance we add an assumption about evolution well above the TeV scale: • The saxion condensate oscillates about the present vacuum in a quadratic potential from an initial amplitude larger than V PQ .
In this case the saxion oscillation, through its coupling with the axion field, enhances the fluctuation of the axion field via the parametric resonance effect [45,46]. Once the fluctuation of the axion field becomes much larger than V PQ , the axion field takes spatially varying random values, leading to the formation of domain walls after the QCD phase transition [47][48][49][50][51][52][53][54]. If N DW = 1, these domain walls are unstable and decay into axions. If N DW > 1, these domain walls are stable, dominate the energy density of the universe, and hence are excluded. For illustration, we consider the "Decoupled Saxion" of section 2.2. Saxion oscillations significantly enhance the fluctuations of the axion field modes with physical wave numbers around k ∼ m s via parametric resonance. When the saxion field has dropped from s I to s and undergone N osc (s) oscillations, the fluctuation of the angular direction is given roughly by where e µ (µ = O(1)) is the growth rate per oscillation. Here the factor of (H inf /2πs I ) is the primordial fluctuation of the angular direction produced during inflation. Assuming the universe is radiation-dominated during these early oscillations, The number of the oscillations grows at small s, as the Hubble scale is smaller. However, this scaling breaks down for s < ∼ V PQ . Additional fluctuations created once s falls below V PQ are small, ∆θ < O(1). This can be seen easily by energy conservation. The energy density of the fluctuation is

JHEP07(2017)125
and cannot be larger than the energy density of the saxion oscillations, m 2 s s 2 . Thus ∆θ < O(1) for fluctuations produced at s < V PQ . The growth in the fluctuation cuts off as s drops below V PQ , and the condition for domain walls not to be produced is δθ(s ∼ V PQ ) < O(1), giving an upper bound on s I , Below the dashed lines in figures 3 and 4 with labels "s I = 10 V PQ ", this condition is violated and we expect axion fluctuations with δθ > O(1). For N DW > 1, the regions below the dashed lines are excluded. For N DW = 1, the axion misalignment angle is randomized and these regions are subject to the constraint given in figure 1, which is shown by pink shadings. These constraints may be avoided in cosmologies that violate the assumption itemized at the beginning of this sub-section.
In the DFSZ theory with N DW > 1, from figure 4 the only allowed parameter region has large V PQ . In four dimensional grand unified theories, symmetries which control the µ term of the Higgs doublets must be broken at the unification scale [55][56][57]. It is illuminating that the constraint from parametric resonance also points towards a large PQ breaking scale.
Note that this constraint is derived by evaluating (3.7) and (3.8) at s ∼ V PQ ; the details of the evolution prior to this is irrelevant. Thus the constraint applies provided the saxion condensate oscillates about the present vacuum in a quadratic potential from an initial amplitude larger than V PQ . However, to phrase the constraint in a way that is independent of the earlier saxion evolution requires a reinterpretation of s I in eq. (3.10). On the left hand side of eq. (3.10), s I is a parametrization of the strength of the saxion condensation, and can be rewritten in terms of T NA via eq. (2.13). In this more general formulation of the bound, the positions of the dashed lines and pink shaded regions are not changed. In the log in the middle of eq. (3.10), s I parametrizes the size of the primordial fluctuation of the angular direction, and should be replaced by the PQ symmetry breaking scale during inflation. It affects the bound only logarithmically.

Displaced vertices, kinks and saxion resonances at colliders
We discuss the following colliders signals resulting from LOSP decays -displaced vertices and kinks involving gravitinos or axinos, and the saxion resonance. For illustration we assume a neutralino (χ 0 ) or right-handed stau (τ R ) LOSP.
The conventional signals of displaced vertices involving gravitinos result from NLSP decay via interactions suppressed by the mediation scale of supersymmetry breaking, with a decay length given by (4.1) Examples include the neutralino (stau) NLSP decaying to a gravitino and h, γ, Z (τ ). This conventional signal applies whenever the LOSP decays to the axino are kinematically -15 -JHEP07(2017)125 forbidden or sufficiently suppressed. In figure 2 and the upper left panel of figure 4 the axino is heavier than the LOSP, and in figure 3 the axino may be heavier than the LOSP. Hence, in all these cases a displaced signal from (4.1) can occur. In the lower left panel of figure 4, this signal competes with LOSP decays to axinos, described below. We may also observe displaced vertices and/or kink signals if the axino is lighter than the LOSP. Any MSSM particles produced at the colliders will first cascade down to the LOSP. Through the axino interactions with higgsinos (gauginos) for DFSZ (KSVZ) theories, the LOSP will decay into the axino and SM particles.
In DFSZ models with a neutralino LOSP, the neutralino decays to the axino and the Higgs/Z boson -(a) of table 1 -with where Cχ 0H is the higgsino component of the neutralino LOSP [23]. This is applicable to all panels of figure 4, except the upper left since we need µ > mã. Below the gray regions, this decay mode is typically more efficient than that to the gravitino final state, as can be seen from the smaller suppression scale m s V PQ Similar to conventional gauge mediation with a stau NLSP, this decay may leave a kink signal. For large stau masses, the stau instead dominantly decays to ν τ +W +ã or τ +Z/h+ã -(c) of table 1. The latter decay is observed as a kink where Z/h is emitted.
In KSVZ models, a neutralino LOSP decays to (γ/Z)ã -(d) of table 1 -and is observed as a displaced vertex. Since this mode is loop-induced the decay rate is smaller that the one in eq. (4.2), typically by a factor of 10 5-6 . On the other hand, a stau LOSP decays into τ Rã -(e) of table 1 -through axino-bino mixing arising from the nonzero electroweak D term, leaving a pure kink. For large stau masses this mixing becomes quadratically suppressed by the electroweak vev, so that the 3-body final state τ R (γ/Z)ã becomes favored -(f) of table 1. The stau also has 3-body decays linearly suppressed by the electroweak vev: τ L (Z/γ)ã and ν τ Wã -(g) of table 1 -the latter has W appearing at a displaced vertex. All of the above modes are suppressed by loop factors as well as three-body phase space factors or the ratio between the electroweak and SUSY scales. The decay of the LOSP into axinos is sub-dominant near the gray-shaded regions of figures 2, 3 and 4 and is dominant far enough from these regions.
If the saxion is heavier than the axino then the axino decays invisibly to an axion and a gravitino. However, an interesting and unique signal arises when m s < mã < m LOSP , e.g. in the upper right panel of figure 4. Since the lower limit on the saxion mass is of order (10 MeV, 1 GeV) for (DFSZ, KSVZ) theories, this can occur in a wide region of parameter space. In this case, the LOSP is produced at the collision point, travels a distance cτ LOSP→ã and decays to the axino and SM particles, leaving a displaced vertex or a kink if cτ LOSP→ã is in an appropriate range. The axino then travels some other distance cτ NLSP→G before it finally decays to the gravitino and the saxion/axion, with the saxion decaying to hh/W W/ZZ/f f /gg for DFSZ and gg for KSVZ, leaving a displaced vertex for an appropriate cτ NLSP→G . Remarkably, the saxion can be observed as a resonance despite its feeble coupling with the SM. This particular decay mode has a distinctive feature of multi-jets from Higgs/Z boson resonances, taus, missing energy, and a saxion resonance.

A warm component of dark matter
We have considered three sources for gravitinos: decoupling from the thermal bath, freezein by higgsino decays, and decays of freeze-in axinos. Given the observed DM abundance, dilution is always large enough that the thermally decoupled gravitinos satisfy warm DM constraints if m 3/2 > ∼ O(keV) [58]. At larger m 3/2 these gravitinos rapidly become cold. On the other hand, gravitinos produced from the FI decays of higgsinos can be warm even if m 3/2 is larger than a keV. These FI gravitinos give a significant component of DM only at low V P Q and m 3/2 , for example near the boundary of the brown regions of figures 3 and 4. For this range of parameter space, FI gravitinos are produced in the MD NA era and hence are diluted less than thermal gravitinos, leading to the larger free-streaming length of eq. (A.5), as discussed in appendix A. If FI gravitinos dominate DM, the warm DM bound on m 3/2 becomes somewhat more stringent, as shown in eq. (A.4). For m 3/2 ∼ O(10 keV) near the brown regions of figures 3 and 4, we predict a mixture of cold and warm dark matter from thermal and FI gravitinos, respectively.

JHEP07(2017)125
In DFSZ theories, the abundance of gravitinos from FI axinos can be comparable to that of thermally decoupled gravitinos near the green regions of figure 4. The FI axinos become non-relativistic before decaying to gravitinos, giving them momenta larger than the thermal gravitinos. When m 3/2 mã, the free streaming length of gravitinos becomes independent of m 3/2 [59] and can be approximated by where the NLSP is the higgsino (axino) in the left (right) panels of figure 4. The analyses from Lyman-α forest [60][61][62] place an upper bound on the dark matter free streaming length, λ FS < ∼ 1 Mpc. As a result, for m NLSP ≈ 650 GeV there is a sizable component of warm dark matter in the parameter space close to the boundaries of the green regions. These warm gravitinos lead to possible signals in the Lyman-α observations. It has been shown that warm dark matter can solve the small scale structure problems although baryon feedback may also play a role. (See [63] for a review.)

Dark radiation
Axions may contribute to dark radiation in both KSVZ and DFSZ theories because the saxion can decay to a pair of relativistic axions via the trilinear coupling in the Kähler potential [32,64]. This decay rate depends on the model-dependent parameter κ = i q 3 where q i and v i are the PQ charge and vev of each PQ breaking field, leading to an effective number of relativistic neutrinos  In KSVZ theories, we take κ < ∼ O(0.1) to be compatible with the current Planck constraint of ∆N eff = 0.6 [65]. A small κ can arise from an approximate Z 2 symmetry or fine tuning. For DFSZ theories, the constraint on κ is much relaxed because the saxion decay to the visible sector is more efficient; in fact, for m s < 2µ, as required to forbid the decay of the saxion into a pair of higgsinos, the Planck constraint is satisfied even if κ ∼ O(1). In all three cases of (4.6), part of parameter space is accessible to the CMB-S4 experiment.

Conclusions
The main results of this paper are shown in figures 2, 3 and 4. Gravitinos that were thermalized early in the universe and later diluted from the decay of a saxion condensate -18 -

JHEP07(2017)125
provide an excellent candidate for dark matter; they are subject to several important constraints, but these leave large allowed regions in the (m 3/2 , V PQ ) plane; as large as shown in figure 2. Results are shown for a very wide range of (m 3/2 , V PQ ) and saxion mass m s , and are independent of almost all UV model-dependence.
In KSVZ theories a large part of the (m 3/2 , V PQ ) parameter space is allowed, although there is a strong upper bound on V PQ for the saxion condensate to decay before BBN, as shown in figure 3. For this upper bound on V PQ to be larger than 10 9 GeV, the saxion mass must be larger than O(1) GeV.
In DFSZ theories the efficient interaction between the axion and Higgs multiplets weakens the upper bound on V PQ , but puts strong constraints on (m 3/2 , V PQ ) from a variety of processes, as shown by shading in figure 4. However, these additional constraints can be removed, for example by making the axino mass sufficiently large. The reduced upper bound on V PQ allows for a saxion mass as small as O(10) MeV, where both BBN and astrophysical bounds are included.
If the saxion begins its oscillation with a field value larger than V PQ , parametric resonance may induce large fluctuations of the axion field. In theories with domain wall number N DW > 1, these large fluctuations lead to the formation of disastrous stable domain walls. This is avoided for large V PQ , which typically has a larger region in DFSZ theories, although in KSVZ theories N DW = 1 is more easily obtained.
We have also estimated the axion abundance from the misalignment angle. For a sufficiently small saxion decay rate, the axion abundance is also diluted. For V PQ > 10 12 GeV and θ mis order unity, some dilution is required, placing an upper bound on the decay rate, as shown in figure 1. In DFSZ theories the higgsino mass is bounded from above accordingly.
Here we summarize possible signals of our scenario: • If the LOSP is lighter than the axino, it will decay to light gravitinos at displaced vertices. We have provided a cosmology with high T R for this well-known signal, for example, of low-scale gauge mediation.
• If the LOSP is heavier than the axino, the axion multiplet participates in the decay chain of MSSM particles produced at colliders (table 1), leaving displaced vertices and/or kinks. If the saxion is lighter than the axino it is produced through axino decay and can be observed as a resonance.
• In some parameter regions gravitinos are dominantly produced via freeze-in processes. Such gravitinos may behave as warm dark matter even if m 3/2 > keV.
• The saxion condensate also decays to relativistic axions, leading to a non-zero dark radiation abundance.
A combination of measurements, especially from displaced vertices or kinks at colliders, could constrain the theory and narrow the prediction for V P Q , which may allow an independent probe from axion physics.

JHEP07(2017)125
A Warm dark matter from freeze-in gravitinos If dark matter is initially thermalized and decouples from the bath while relativistic, m DM < T dec , its energy at decoupling is of order the decoupling temperature T dec . A sufficiently light dark matter particle will affect structure formation via its large velocity at matterradiation equality, leading to a lower bound on m DM [58]. Nonetheless, this bound is different when dark matter is produced from freeze-in decays during a MD era instead of thermal decoupling during a RD era. In this section, we investigates how the freeze-in scenario affects the constraint.
In general, the DM abundance today is related to the initial one at production by where D is the dilution factor, s the entropy density, and the yield is parametrized by in units of the thermal equilibrium value Y th of relativistic Weyl fermions. For freeze-in (thermal decoupling), < 1 ( = 1). We are concerned with the free-streaming length of dark matter so we study how dilution affects the momentum red-shift; where D is substituted using eq. (A.1). This gives the momentum p f at any temperature For gravitinos that originate from higgsino decays at the freeze-in temperature T FI = mH /x FI , with x FI ∼ 2 − 5, p i mH /2, whereas p i T dec for those that decouple from the thermal bath at temperature T dec . Specifically, the ratio p i /s 1/3 i in eq. (A.3) becomes a constant for each production mechanism and is larger by a factor of x FI /2 for freezein than thermal decoupling. This implies that the free-streaming length in the case of freeze-in is enhanced in comparison. The constraint from Lyman-α gives a lower bound on p f /m DM and thus a lower bound on m DM ∝ (x FI /2) 3/4 −1/4 . Although the freeze-in gravitino phase space distribution is different from that of thermally decoupled gravitinos, we expect a bound of similar order applies. Therefore, the constraint on m 3/2 for freeze-in can be obtained from rescaling the result of ref. [58] that assumes thermal decoupling thermal decoupling

JHEP07(2017)125
One can estimate the yield for freeze-in production by In the brown regions of figures 3 and 4, FI gravitinos are overproduced. We investigate the parameter space immediately above these brown regions so that the warm FI gravitino abundance is sizable. In this case, gravitinos freeze-in during the MD A era so H FI ≈ π g * (T FI ) 3 √ 10 Note that is necessarily less than unity because the FI yield cannot exceed the thermal value. If the above expression gives > 1, gravitinos stay in thermal equilibrium until T ∼ mH , where they decouple as the higgsino LOSPs becomes non-relativistic and exponentially depleted. Using the numerical results of T NA in figures 3 and 4, one finds that m 3/2 ∼ O(10 keV) can lead to warm freeze-in gravitinos in the parameter space near the edges of the brown regions.