Axino dark matter with low reheating temperature

We examine axino dark matter in the regime of a low reheating temperature, TR, after inflation and taking into account that reheating is a non-instantaneous process. This can have a significant effect on the dark matter abundance, mainly due to entropy production in inflaton decays. We study both thermal and non-thermal production of axinos in the framework of the MSSM with ten free parameters. We identify the ranges of the axino mass and the reheating temperature allowed by the LHC and other particle physics data in different models of axino interactions. We confront these limits with cosmological constraints coming the observed dark matter density, large structures formation and big bang nucleosynthesis. We find a number of differences in the phenomenologically acceptable values of the axino mass mã and the reheating temperature relative to previous studies. In particular, an upper bound on mã becomes dependent on TR, reaching a maximum value at TR ≃ 102 GeV. If the lightest ordinary supersymmetric particle is a wino or a higgsino, we obtain a lower limit of approximately 10 GeV for the reheating temperature. We demonstrate also that entropy production during reheating affects the maximum allowed axino mass and lowest values of the reheating temperature.


Introduction
A variety of cosmological data points to the existence of a non-luminous component of the matter in the Universe, referred to as dark matter (DM). In spite of decades-long efforts aiming at detecting DM particles directly, they have so far remained elusive. Indirect probes, such as the impact of DM on the formation and evolution of large structures in the Universe or on the spectrum of the temperature fluctuations in the cosmic microwave background provide powerful tools to study DM. However, such probes involve only some properties of DM particles, such as their (time-dependent) energy density and free-streaming length, hence the masses of proposed DM candidates span thirty orders of magnitude while their cross-sections for scattering with known particles -forty orders of magnitude (see, e.g., [1] for a recent review).
There are many examples of DM candidates which were not introduced ad hoc, but possess a sound theoretical motivation. They include bosonic fields in coherent motion, such as axions, weakly interacting massive particles (WIMPs), such as the lightest neutralino of the Minimal Supersymmetric Standard Model (MSSM), or extremely weakly interacting particles (EWIMPs), such as gravitinos and axinos.

JHEP11(2015)139
The axion was originally introduced as a solution to the so-called strong CP problem. The smallness of the electric dipole moment of the neutron can be understood if there is a U(1) PQ symmetry, referred to as Peccei-Quinn (PQ) symmetry [2,3], which is spontaneously broken at a scale f a . The very light pseudo-Goldstone boson a associated with this symmetry, called an axion, couples to the gluon anomaly [4] (see also, e.g., [5] for a review), where α s = g 2 s /4π and g s is the strong coupling constant. Various cosmological considerations suggest that f a lies between approximately 10 9 GeV and 10 12 GeV [6].
Two broad frameworks for field theory models of axions have been considered in the literature. In the Kim-Shifman-Vainshtein-Zakharov (KSVZ) model [7,8], the gluon anomaly term (1.1) is induced through a loop of very heavy singlet quarks carrying PQ charges, while in the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) model [9,10] the PQ charges are assigned to the Standard Model fields. These assumptions are sufficient to determine unambiguously axion interactions with the Standard Model fields that are relevant for cosmology, though specific implementations of these models can be rather complicated [5,11]. The Lagrangian describing axion interactions with Standard Model particles is given in appendix A.
In supersymmetric models [12] the axion resides in a chiral multiplet with a fermionic superpartner -the axinoã -whose interactions with the Standard Model particles are related to those of the axion through supersymmetry. See [1,13] for a review, relevant formulae are shown in appendix A. The properties of the axino in models with softly broken supersymmetry can be relevant for cosmology [14][15][16], especially when it is the lightest supersymmetric particle (LSP) and constitutes cold DM [17,18]. Particular scenarios of supersymmetry breaking give predictions for the axino mass [19][20][21]; however, due to a strong model dependence and the absence of an universally accepted scheme of supersymmetry breaking, we adopt a phenomenological approach and treat the axino mass as a free parameter.
The mechanisms for axino generation in the post-inflationary Universe (see [22] for a review) are analogous to those for the gravitino: thermal production (TP) from scatterings and decays of other particles in thermal equilibrium and non-thermal production (NTP) from out-of-equilibrium decays of heavier particles. Although detailed predictions for the axino abundance from TP, Y TP a , strongly depend on the model of axion interactions, it is inversely proportional to f 2 a and does not depend on the axino mass. In this respect, the axino differs significantly from the gravitino: as the latter contains a goldstino component related to broken supersymmetry, the gravitino abundance from TP is inversely proportional to M 2 P and to the square of the gravitino mass. A notable difference between the axino interaction models is that in KSVZ models Y TP a is proportional to the reheating temperature T R defined in terms of the inflaton decay rate Γ φ = π g * (T R )/90 (T 2 R /M Pl ), while in DFSZ models Y TP a does not depend on T R . This can be attributed to the fact that in DFSZ models axinos are mainly produced from decays of thermal particles, while in KSVZ models the main source of axinos are scatterings of strongly interacting particles. Existing numerical analyses [23,24] suggest that, in order not to overclose the Universe

JHEP11(2015)139
in the KSVZ scheme, T R should be at most a few orders of magnitude larger than the electroweak scale.
The contribution of the NTP to the axino abundance comes from decays of the lightest ordinary supersymmetric particles (LOSP) which underwent freeze-out, hence Y NTP a = Y LOSP , or in terms of cosmological parameters, Depending on the masses of the axino and the LOSP, the LOSP lifetime may be so long (above seconds) that the highly energetic LOSP decay products can affect successful predictions of the big bang nucleosynthesis (BBN).
Although there has been remarkable progress in the understanding of axino physics and cosmology in the last two decades, the latest LHC data and the increasing precision of reconstructing the history of the Universe motivate extending the existing analyses in two ways. Previous studies of axino DM focused on axino interactions, while some specific assumption about the MSSM mass spectra were made for simplicity. Also, the low T R regime was studied assuming instantaneous reheating and some fixed typical values of the abundance of axino DM originating from NTP. Given many constraints that the LHC data put on the MSSM, it is worthwhile to ask what are the allowed and excluded ranges of the parameters of the general MSSM with axino DM. It is also known that the altered expansion rate of the Universe during reheating may result in a significant change of DM abundance [25][26][27][28][29][30][31][32][33][34].
In this paper, we address the two issues mentioned above. We identify the phenomenologically viable parameter ranges of the MSSM with axino DM and calculate accurately the axino abundance [35] not only during the radiation dominated (RD) period, but also during the phase of reheating after inflation [33]. We extend the analysis presented in [31] firstly by taking into account NTP of axinos for various candidates for the NLSP. In addition, we include in our study possible U(1) Y contributions to axino TP that can play an important role for light neutralino NLSP in the regime of low T R . We also present results for a wide range of SUSY spectra satysfying current experimental constraints.
The paper is organized as follows. In section 2, we examine how non-instantaneous reheating affects the predictions for axino DM. In section 3, we show the results of our numerical study of a 10-parameter version of phenomenological MSSM (p10MSSM) and identify the ranges of values of the axino DM mass and the reheating temperature that are consistent with data including large scale structure formation and big bang nucleosynthesis (BBN) constraints. We present the conclusions in section 4. A number of technical issues are relegated to the appendices. Appendix A summarizes the interaction Lagrangians of the Standard Model or MSSM fields with the axion and the axino. Appendix B presents some details of our calculation of TP of axinos and appendix C describes some phase space integrals necessary to address the free-streaming length of axinos from NTP. In appendix D, we describe the details of our numerical procedure.

Thermal production of axinos
Our first aim is to understand the effects of low T R and of non-instantaneous reheating on axino abundance. Therefore, we will first discuss TP of axinos, focusing mainly on the KSVZ model. For concreteness, we will present the results obtained for the gluino and squark masses of mg = mq = 1 TeV, and we shall assume that the C aW W = 0 (see appendix A for details regarding parameters and notations), unless indicated otherwise. We shall relax this assumption about the MSSM mass spectrum later.
In the high-T R regime, the scatterings associated with the SU(3) c group dominate axino TP. 1 Were the reheating process instantaneous, we could identify the reheating temperature T R with the temperature T RD at which the radiation dominated epoch began. However, with the RD epoch having been preceded by a reheating epoch, during which the energy density of the oscillating and decaying inflaton dominated the Universe, there is an additional contribution to Y TP a originating from a modified relation a(t) ∼ T −8/3 between the scale factor of the Universe a(t) and the temperature T . A straightforward, but technical and rather involved calculation presented in appendix B leads to the conclusion that this additional contribution is about 1/6 of the standard high T R result. Loosely speaking, the existence of a reheating phase before the RD epoch effectively extends the period of relic production and available temperature range -and therefore the axino abundance increases.
One comment is in order here. In ref. [32], which studied the same situation as described above, a constant reduction of Y TP a by a factor of about 1/4 was reported. This apparent difference results from a different choice of a reference point. In ref. [32] the results obtained for instantaneous and non-instantaneous reheating were compared for the same value of T R , while we believe that it is more appropriate to compare both abundances at the same value of T RD , which (in contrast to T R being a convenient shorthand notation for the inflaton decay rate) has a clear physical meaning of the temperature that marks the beginning of the standard radiation dominated epoch. As shown in appendix B, we find an approximate relation T RD ∼ T R /2, so we can use these two temperatures interchangeably when referring to orders of magnitude.
For intermediate values 10 2 GeV < ∼ T RD < ∼ 10 4 GeV, there is a phase space suppression of the scattering terms associated with TeV-scale superparticles. For a given value of T RD this effect is smaller in the case of non-instantaneous reheating, as the additional contribution to Y TP a stabilizes the result and the ratio of the abundances calculated for non-instantaneous and instantaneous reheating becomes slightly larger than 7/6 ≈ 1.17. This can be seen in the left panel of  For T RD smaller than the masses of strongly interacting particles (herein 1 TeV) the ratio falls, as the contribution to axino TP from scatterings becomes small with respect to the contribution from squark and gluino decays. This is because, unlike scatterings, the decays depend rather weakly on the details of reheating as shown in the right panel of figure 1. In principle, larger temperature values attainable with non-instantaneous reheating may lead to a larger equilibrium number density of decaying squarks or gluinos and therefore to an increase of Y TP a , but this effect is less important than the phase space suppression in the case of scatterings. Hence, as long as Y TP a is dominated by decays, the ratio in the left panel of figure 1 falls below 1.17. Figure 2 shows examples of the dependence of axino TP on the gluino and the squark masses.
The increase of the ratio of the abundances for the instantaneous and noninstantaneous reheating scenarios can be sizable only for T RD < ∼ 100 GeV, if the parameter C aY Y parametrizing the coupling between the axino and the U(1) Y gauge bosons and gauginos is equal to zero or if bino-like neutralino is heavier than about 500 GeV. Otherwise, axino TP is at low temperatures dominated by the decays of the light bino which is again rather insensitive to the details of reheating if T RD is comparable to the bino mass.
However, for such low values of T RD , even for vanishing C aY Y and/or heavy bino, TP often gives an abundance which is a very small fraction of that required for DM, so details of reheating are irrelevant in that case. One should also remember that for very low values of T RD axinos are decoupled from kinetic but not from chemical equilibrium and their distribution function can differ from that for the equilibrium case [36], but for axinos this only happens for the DM abundance dominated by NTP, so this effect has no consequences for our study.
We can see that the main effect of non-instantaneous reheating on the axino TP originates from a modified contribution from scattering. Therefore, the details of reheating  play a minor role in DFSZ models, as in this case TP is typically dominated by thermal higgsino decays [23,24].

Non-thermal production of axinos
For low enough values of T RD axino TP is highly suppressed and it is the non-thermal contribution to the axino relic density that dominates. However, NTP is also affected by an additional entropy production during the reheating period provided that T RD is low enough so that the LOSPs freeze out before the RD epoch (see, e.g., a recent discussion about similar issue in the case of gravitino DM [33]). 2 This results in a suppression of Ω LOSP h 2 below the value Ω LOSP h 2 (high T RD ) obtained in the standard cosmological scenario where the LOSP freeze-out occurs in the RD epoch. Consequently, for a fixed mã/m LOSP ratio in eq. (1.2), Ω NTP a h 2 also decreases. The lower the reheating temperature, the longer the period between the LOSP freeze-out and the beginning of the RD epoch, so the LOSPs are effectively diluted. As a result Ω NTP a h 2 decreases with T RD . This is illustrated in figure 3 where we show both TP and NTP contributions to Ωãh 2 for two selected SUSY spectra with bino LOSP mass equal to 100 GeV and ∼ 1 TeV, while squark and gluino masses are set to 1 TeV. The rest of the SUSY spectrum is in both cases chosen such as to obtain LOSP yield at freeze-out corresponding to ΩBh 2 (high T R ) = 0.1 or 10 4 .

BBN and WDM constraints
The BBN constraints for the axino can be analyzed similarly to those for the gravitino (see e.g. [39][40][41][42][43][44][45][46]); they are typically mild. This is because the lifetime of the LOSP decaying to the axino usually hardly exceeds 0.1 sec., unless one considers very light LOSP, allows for a very strong mass degeneracy, mã ≈ m LOSP or considers a LOSP whose 2-body decays to  axinos are forbidden (i.e. bino with C aY Y = 0). However, if the neutralinos are too light, they have a very small relic abundance, which is further suppressed by a small value of T R , and there are simply too few decaying LOSPs around to put primordial nucleosynthesis in danger (in this case, the correct axino DM abundance must originate from TP). The only exception is a scenario with a very light axino and a very light bino LOSP with extremely small annihilation rate, suppressed by large slepton masses. In this case, the relic abundance of bino LOSPs can be very large and the correct axino DM relic density is obtained from NTP by a suppression with a small axino mass. A more dramatic change appears when we set C aY Y = 0, thereby assuming that the axino interacts only with the SU(3) c gauge sector. For bino LOSP, the 2-body bino decay into axino and a photon or Z boson are then disallowed and the dominant bino decay channel B → qqã involves an effective axino-quark-squark interaction vertex discussed in [47]. Applying the formulae obtained therein, we can estimate the bino lifetime as (see also appendix C) (2.1) With a longer LOSP lifetime, the BBN constraints are significantly more severe. For small enough values of mã, axinos may become warm dark matter (WDM) with a free streaming length which is too large to account for observed structures in the Universe. 3 Thermally and non-thermally produced axinos have different rms velocities and hence these constraints from structure formation depends on the fraction of the WDM axinos in DM (for a discussion about non-thermally produced WDM see, e.g., [49,50]). We use the 95% CL JHEP11(2015)139 exclusion limits from ref. [51]. For TP domination these limits translate to mã < ∼ 5 keV, while the case of NTP domination and the mixed case require a more careful analysis In particular, for C aY Y = 0 the calculation of the rms axino velocity requires an integration over 3-and 4-particle phase space (bino and stau LOSPs, respectively); we provide technical details of this computation in appendix C.

Axino with low reheating temperature in the MSSM
Having discussed the predictions for TP and NTP of axinos for different values of the reheating temperature T R and for different patterns of the MSSM mass spectra, we are now ready to present and discuss the results of our extensive numerical scan. There are shown in figure 4 for f a = 5 × 10 9 GeV and 10 11 GeV. In comparison with the results of previous studies, we find a number of differences in the shape of the allowed region of the axino mass and the reheating temperature. We discuss them separately for the bino LOSP and for the wino and higgsino LOSP.

Bino LOSP
Firstly, similarly to ref. [35] we find an upper bound on T R and a lower bound on mã, beyond which axino DM becomes too warm. For T R < ∼ 10 6 GeV the region allowed by the BBN or WDM constraints extends to smaller values of mã ∼ 10 −4 GeV than in previous analyses. For T R < ∼ 10 3 GeV these points correspond to dominant NTP from a very large relic density of bino LOSP with a very small annihilation cross-section dominated by a t-channel exchange of multi-TeV sleptons. Such points, are however excluded by both BBN constraints (bino LOSP has a long lifetime and a sizable hadronic branching fraction, cf. [52]) and by WDM constraints (with dominant NTP, one needs mã > 30 MeV). Eventually, the lower bounds for mã at fixed T R are similar to those obtained in [35], but in our analysis they result from incorporating additional cosmological constraints.
Secondly, looking at the rightmost part of each plot in figure 4, we see an upper bound on mã. In this respect our result visually resembles that of ref. [35], but this similarity follows from very different physical assumptions. In ref. [35] three typical choices of fixed axino abundance from NTP, Y NTP a , were considered, which led to upper bounds on axino mass for Y NTP a > 0. When NTP was the dominant source of axino DM, these bounds on mã did not depend on T R . In contrast, in our analysis we obviously do not make such a restrictive assumption and the maximum value of mã that we obtain is related to the maximum value of the LOSP mass. More precisely, the largest value of mã is obtained for T R ∼ O(10 2 GeV), corresponding to the freeze-out temperature of the bino LOSP. This maximum mã corresponds to the largest available values of the bino LOSP mass and the largest available slepton masses (cf. table 1). These two features have the same effect: one expects a larger relic density for heavy particles and there is also a big suppression of the annihilation cross-section due to large masses of the intermediate particles. We show these aspects of our results in the left panel of figure 5, where we show how the allowed region changes with different assumptions about the largest possible bino and slepton masses. Furthermore, we see a decrease of the maximum allowed mã for T R falling below O(10 2 GeV). This can be understood taking into account that during reheating there is entropy production due to inflaton decays. As a consequence, the LOSP relic density becomes suppressed, if T R falls significantly below the freeze-out temperature. For a given value of T R the suppression is stronger for binos with larger masses. This confines the allowed region to smaller values of bino mass and, as mã < m B , to smaller values of axino mass.
We also note that, for T R < ∼ 10 4 GeV the calculations for TP of axinos cannot be fully trusted, as the SU(3) c gauge coupling becomes large (see [13] and references therein). This poses no problem for T R < ∼ 10 2 GeV, as NTP of axinos is then dominant, but in the window 10 2 GeV < ∼ T R < ∼ 10 4 GeV the upper bounds on mã or T R should be treated as approximate. Changing fã mostly leads to a shift of the allowed region in the (mã, T R ) plane, as shown in the right panel of figure 4.
The difference between regions of the (mã, T R ) plane allowed in the KSVZ and DFSZ models is presented in the right panel of figure 5. For large values of T R > ∼ 10 5 GeV it can be traced back to the fact that Ω TP a h 2 scales as mãT R and mã for the KSVZ and DSFZ models, respectively [23,24]. For smaller values of T R the bulk of the allowed region is very similar for both models. For a fixed mã > ∼ 1 MeV, there is still a small difference between the largest allowed T R corresponding to the TP dominance. This results from different sources of thermal contributions to axino DM from decays. In KSVZ models the most important decays are those of colored particles; these are loop-suppressed with respect to decays of neutralinos with a non-negligible higgsino fraction in DFSZ models. Therefore, in DFSZ models TP of axinos from decays is much more efficient which may lead to overproduction of axino DM for values of mã and T R corresponding to the correct axino DM density within the KSVZ scheme.

Wino, higgsino and stau LOSP
In figure 6 we show the allowed regions of the (mã, T R ) plane for the wino and higgsino LOSP and in figure 7 the same for the stau LOSP, with constraints imposed in the same way as described in section 3.1. We again show the results for two values of f a = 5 × 10 9 GeV and 10 11 GeV. We restrict our analysis to the KSVZ model, since we have seen that for low T R the predictions of the KSVZ an DFSZ models do not differ very much. The main difference with respect to the bino LOSP case (figure 4) is that T R is now bounded from below. This is because for winos, higgsinos and staus (unlike for binos) the masses of the states mediating annihilations cannot be much larger than the masses of annihilating particles. Hence, the annihilation cross-section of these particles cannot be made very small by increasing the mass of the intermediate states. Nonetheless, by varying MSSM parameters, we find examples of axino DM and the higgsino or wino or stau LOSP in which the observed DM abundance originates mainly from TP (upper parts of the allowed regions) or NTP (lower parts of the allowed regions). This may seem at odds with a recent analysis of axino DM in the DFSZ model with the higgsino LOSP [34] which studied both TP (freeze-in of axinos) and NTP (LOSP freeze-out), and concluded that TP always dominates. This apparent difference can be traced to the fact that in ref. [34] two examples of the MSSM mass spectrum are studied, while our numerical analysis extends to larger values of the masses of the supersymmetric particles and thus allows larger LOSP JHEP11(2015)139  relic abundances. This in turn can give rise to the observed axino DM abundance via NTP even with additional entropy production from inflaton decays.
In the case of the stau LOSP one can obtain significant contribution to the axino relic density from TP also for low values of T R ∼ 10 GeV. It is due to decays of light stau LOSPs being still in thermal equilibrium. The mass of the axino is then typically significantly lower than mτ 1 , which suppresses the NTP contribution. However, one needs to take into account that such a scenario is constrained by the LHC searches for direct production of staus with missing energy [53]. When treating this we calculate the relevant production cross section with MadGraph5 aMC@NLO [54].

Direct and cascade decays of the inflaton field
In our analysis so far, we have made an assumption that there are no direct and cascade decays of the inflaton field to axino DM particles. Here, we would like to study the impact of such decays on the allowed ranges of axino mass and reheating temperature, following the model-independent approach used in [28,29].
Our results are presented in figure 7 and 8 where we plot the allowed regions for the stau, bino and higgsino LOSP for different values of the dimensionless parameter η = b · (100 TeV/m φ ), where b is an average number of axino DM particles per inflaton decay and m φ is the inflaton mass at the minimum of the potential. As inflaton decays provide an additional non-thermal component of axino DM (see a recent discussion in the case of gravitino [33]), the allowed region becomes extended towards smaller values of T R at largest allowed values of mã. This is because the additional NTP from inflaton decays allows for a smaller contribution to axino DM density originating from LOSP decays, hence -for a fixed set of the MSSM parameters -for a smaller T R and a larger suppression of LOSP abundance by entropy production in inflaton decays.

Conclusions and outlook
In this paper we have studied the impact of a low reheating temperature T R on thermal and non-thermal production of axino DM, taking into account the non-instantaneous nature of the reheating process. We also extended previous studies by analyzing wide ranges of phenomenologically acceptable parameters of the 10-parameter version of phenomenological MSSM instead of presenting the results for a single typical parameter choice. Comparing our results with previous works, we found a number of differences in the allowed ranges of the axino mass mã and the reheating temperature. In particular, depending on the choice of the axion model and the choice of the MSSM parameters, we showed that BBN

JHEP11(2015)139
constraints can exclude large portions of the parameter space corresponding mainly to nonthermal production of axino DM relevant for low T R . We also demonstrated how entropy production during reheating affects the upper limits on the axino mass for a given range of the MSSM parameters.
The are a few directions in which the analysis presented herein could be extended. In our work, we relaxed the simplified assumption of instantaneous reheating, we still required that the maximum temperature during the inflaton-dominated period was sufficiently larger than the reheating temperature that the supersymmetric particles could reach thermal equilibrium. Although it is a realistic requirement, one can also envision scenarios with a very small energy density at the end of inflation. Additionally, it has recently been noted that the maximum temperature during reheating may not be as large as previously estimated [55]. Such cases are not included in our study and we leave them for future work.
An issue which requires some care in models with low T R is the origin of the primordial baryon asymmetry. Since we work in a supersymmetric setup, the Affleck-Dine mechanism [56] (see, e.g., [57] for a review) is a feasible options, though a detailed discussion is beyond the scope of our study.

Acknowledgments
We thank K.-Y. Choi, L. Covi and K. Kohri for discussions and comments. ST would also like to thank E. M. Sessolo for a helpful discussion regarding the LHC constraints. This work has been funded in part by the Welcome Programme of the Foundation for Polish Science. LR is also supported in part by a STFC consortium grant of Lancaster, Manchester, and Sheffield Universities. ST is also partly supported by the National Science Centre, Poland, under research grant DEC-2014/13/N/ST2/02555. The use of the CIS computer cluster at the National Centre for Nuclear Research is gratefully acknowledged.

A Axion and axino interactions
The effective axion interaction Lagrangian after integrating out all heavy PQ charged fields can be written, to the lowest order terms in 1/f a , as where (following a partial integration over on-shell quark fields) the c 1 term can be reabsorbed into the c 2 term. The KSVZ case can be identified with c 1 = 0, c 2 = 0, c 3 = 0, while the DFSZ one with c 1 = 0, c 3 = 0, c 2 = 0. General axion models can have both c 2 = 0 and c 3 = 0. C aW W and C aY Y are model-dependent parameters that correspond to axino-gaugino-gauge boson anomaly interactions for the U(1) Y and the SU(2) L groups, respectively. In a supersymmetrized version of an axion model [14][15][16] the real scalar axion field a resides in a chiral supermultiplet since it is a gauge singlet. The other members of the JHEP11(2015)139 axion supermultiplet are the fermionic superpartner axinoã and the real scalar field saxion s that provides a remaining bosonic degree of freedom on-shell.
The interaction Lagrangian for the axion supermultiplet can be obtained by supersymmetrizing eq. (A.1). In particular, the axino-gaugino-gauge boson and the axino-gauginosfermion-sfermion interaction terms are given by [18,35] wheref D andf denote sfermions carrying non-zero T 3 and Y , respectively. A generic form of interactions between the axion and matter supermultiplets was considered in [58]. In particular, it was pointed out that, for v PQ > T M Φ , where M Φ is the mass of the heaviest PQ-charged and gauge-charged supermultiplet Φ, the axino-gauginogauge boson interaction term is suppressed by M 2 Φ /T 2 . This is particularly important for the DFSZ axino, where Φ corresponds to the Higgs supermultiplets and therefore M Φ = µ (the higgsino mass). The dominant contribution to axino TP is then associated with a higgsino decay to the axino and the Higgs boson that is described by [24,58,59] L ef

B Calculation of axino TP with low reheating temperature
In scenarios with non-instantaneous reheating in calculating Y TP a one has to take into account a modified expansion rate of the Universe. Below we briefly describe the methodology that can be used to calculate the axino TP yield.
Axino TP yield with non-instantaneous reheating. It results in a modification of temperature dependence on the scale factor T (a). The Boltzmann equation can then be written as dXã dT where Xã = a 3 nã . . In that case the present-day axino abundance can be written as: where T up corresponds to an effective upper limit in the integration (in practice it is sufficient to use T up (5 − 10)T RD , as for larger temperatures TP of axinos is more efficient, but the fast expansion of the Universe in the reheating period dilutes away all the axinos produced at that early times). We also used A = a/a I = aT R and where R is related to the energy density of radiation ρ R by R = ρ R a 4 .

JHEP11(2015)139
The scattering term. In order to deal with the scattering contribution to Y TP a we express Σ scat in terms of the scattering cross-section σ(s), following [60]), and obtain: We then change the order of integration and decompose the above integral into where I 1 and I 2 are given by expressions very similar to (B.4), but with T integrated from T 0 to (m 1 + m 2 )/x and from (m 1 + m 2 )/x to T up , respectively. Then I 2 ≈ 0 since T 0 ≈ 0. For the remaining integral I 1 one obtains A careful analysis of eq. (B.7) in the reheating period shows that one can make an approximation where c −7 and T RD ∼ 0.5 T R . In practice we find more exact values of a and T RD numerically, but they depend only slightly on the model parameters.
High T R limit of the scattering term. The integral (B.6) can be rewritten as a sum of three integrals, schematically represented by One can verify that in the limit T RD → ∞, we have J 1 → 0, since the range of external integration shrinks to zero, while the integrand does not diverge. The second integral, J 2 , corresponds to the the standard result obtained in the instantaneous reheating approximation. In the limit of high reheating temperature the inner integral can be simplified to

JHEP11(2015)139
(σ depends on t via m eff , see, e.g., [35]) where we noticed that the integral is mainly determined by the values of the integrand in high-s limit in which, to a good approximation, σ(s, t) = σ(t). For the third integral, J 3 , we similarly note that inner integration leads to In the integration range T = √ s/t > T RD and therefore (B.11) becomes where we assumed T up = cT RD with c high enough so that effectively T up can be replaced by ∞ in the integration. The remaining (external) integrals for both J 2 and J 3 are the same. Hence for high T R J 3 J 2 Eq. (B.13) is valid for each contribution to the scattering term.
The decay term. In the case of the decay term we substitute Γ n eq i (see, e.g., [60]) into (B.2) and obtain Once again we change the order of integration and find that one term is negligible, while the other leads to where f is given by eq. (B.8) with √ s replaced by E. The inner integral where m i /T up ≤ t ≤ ∞ can be calculated analytically. Depending on the value of t one obtains for (B.17) and for t ≥ m i /T RD (the temperature can be either larger or smaller than T RD ) One can verify that in the case of instantaneous reheating the standard result [60] is rederived.

C Phase space integrals for non-thermally produced axinos
In this appendix we provide results for both the LOSP lifetime and the present-day rms velocity of axinos in the case of C aY Y = 0. Depending on the nature of the LOSP this requires analysis of 3− of 4−body decays. It is straightforward to verify that eq. (C.1) can be simplified to eq. (2.1) for m 2 B m 2 q . When treating the WDM constraint we calculate the present-day rms velocity [50] by where C R = 1 for the "right" stau and C L 4.45 for the "left" stau, while functions f (x), g(x) and h(x) are shown in figure 10. They all tend to 1 for x 1, i.e., for mã mτ m B , mq. The present-day rms velocity is equal to h(x) Figure 11. Auxiliary functions in the formula for the present-day rms velocity of axinos produced in late-time stau decays if C aY Y = 0.

D Description of the numerical analysis
Here we explain some details of our numerical analysis of the scenario of axino DM with low reheating temperatures of the Universe in the context of the MSSM. A study of a completely general MSSM would be at the same time complicated and unnecessary, hence we select a 10-parameter version of the MSSM (p10MSSM) which has practically all the relevant features of the general model. These adjustable parameters of the model and their ranges are specified in table 1. Our choice is closely related to that of [62] (see discussion therein), except that we keep both the wino mass M 2 and the bino mass M 1 free. We scan the parameter space of p10MSSM following the Bayesian approach. The numerical analysis was performed using the BayesFITS package which utilizes Multinest [63] for sampling the parameter space of the model. Mass spectra were calculated with SOFTSUSY-3.4.0 [64], while B-physics related quantities with SuperIso v3.3 [65].
The constraints imposed in scans are listed in table 2. The LHC limits for supersymmetric particle masses were implemented following the methodology described in [62,73]. The DM relic density for low T R was calculated by solving numerically the set of Boltzmann equations, as outlined in [25,33]. In order to find the point where WIMPs freeze out, we adapted the method described, e.g., in [74] to the scenario with a low reheating temperature, extracting the relevant functions entering the Boltzmann equations ( σv eff and σv eff E eff , cf. ref. [33]) with appropriately modified MicrOMEGAs v3.6.7 [75].
The LUX upper limits [72] have been implemented as a hard cut.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.