Probing PeV scale SUSY-breaking with Satellite Galaxies and Primordial Gravitational Waves

We study an inevitable cosmological consequence in PeV scale SUSY-breaking scenarios. We focus on the SUSY-breaking scale corresponding to the gravitino mass $m_{3/2}=100{\rm eV}-1{\rm keV}$. We argue that the presence of an early matter-dominated era and the resulting entropy production are requisite for the Universe with this gravitino mass. We infer the model-independent minimum amount of the entropy production $\Delta$ by requiring that the number of dwarf satellite galaxies $N_{\rm sat}$ in the Milky Way exceed the currently observed value, i.e. $N_{\rm sat}\gtrsim63$. This entropy production is inevitably imprinted on the primordial gravitational waves (pGWs) produced during the inflationary era. We study how the information on the value of $\Delta$ and the time of entropy production are encoded in the pGW spectrum $\Omega_{\rm GW}$. If the future GW surveys observe a suppression feature in the pGW spectrum for the frequency range $\mathcal{O}(10^{-10}){\rm Hz}\lesssim f_{\rm GW}\lesssim\mathcal{O}(10^{-5}){\rm Hz}$, it works as a smoking gun for PeV SUSY-breaking scenarios. Even if they do not, our study can be used to rule out all such scenarios.


I. INTRODUCTION
One of the advantages of studying gravitino cosmology lies in the fact that the gravitino mass m 3/2 is directly related to a supersymmetry (SUSY)-breaking scale √ F , i.e. m 3/2 F/( √ 3M P ) where M P 2.4 × 10 18 GeV. For that reason, understanding effects that gravitinos can potentially have on experimental data can be of great help in studying the SUSY-breaking scale. Then, is there any inevitable physical effect induced by gravitinos so that its absence implies exclusion of a certain SUSY-breaking scale?
Inspired by this question, in this work, we give our special attention to PeV (10 6 GeV) scale SUSY-breaking scenarios where the gravitino mass is so light that gravitinos serve as the lightest SUSY particle (LSP). Particularly, we focus on the gravitino mass range 100eV m 3/2 1keV. 1 In this case, as far as a reheating temperature is greater than a sparticle mass, it is certain that gravitinos were once produced in the thermal bath and exist today in the form of warm dark matter (WDM) with a free-streaming length amounting to O(0.1)Mpc.
Interestingly, for the thermal sub-keV gravitino, the relic abundance Ω 3/2 h 2 is insensitive to the reheating temperature T RH [2], but sensitive to only m 3/2 and the decoupling temperature T 3/2,dec . This means that once the sparticle mass spectrum is fixed, there is a definite prediction for the relic abundance of the thermal sub-keV gravitino WDM. Given the null observation of any sparticle thus far, it is fair to state that the gravitino decoupling temperature is at least O(1)TeV. Therefore, for 100eV m 3/2 1keV, there exists a solid lower bound on Ω 3/2 h 2 ≡ ω 3/2 for each m 3/2 as far as T RH is larger than sparticle masses and no entropy production occurs.
On the other hand, Lyman-α forest observations and the redshifted 21cm signals in EDGES observation [3] give stringent lower bounds on the WDM mass in the case where the whole of DM population consists only of a WDM component. Constraints from each experiment read m wdm > 5.3 keV [4] and m wdm > 6.1 keV [5,6], respectively (See [7] for more conservative lower bound, 1.9keV). Hence, as is well-known, sub-keV gravitino WDM cannot be responsible for 100% DM population today if the Universe went through PeV scale SUSYbreaking. However, since these constraints assume that the whole DM is made of a single species, there is still room for the sub-keV gravitino to be part of DM if it is a subcomponent. Thus, an immediate and natural question is "what is an upper bound on ω 3/2 for a given m 3/2 ∈ [100eV, 1keV] to be consistent with cosmological data at small scales?" In this work, we address this question by making the estimate of the expected number of dwarf satellite galaxies N sat in the Milky Way resulting from the Universe with the mixture of cold dark matter (CDM) and gravitino WDM. Defining the fraction of DM population contributed by the gravitino WDM to be f 3/2 ≡ ω 3/2 /ω DM , we compute N sat based on the matter power spectrum arising from 100(1 − f 3/2 )% CDM and 100f 3/2 % grav-itino WDM. Then by requiring N sat 63, we obtain the maximally allowed f 3/2,max for each m 3/2 . Note that we demand that the rest of the main component DM is of a cold type to avoid too much suppression of the matter growth at small scales. Furthermore, we have the CDM unspecified in our work because it depends on details of a model.
In accordance with our observation that f 3/2,max is smaller than f 3/2 expected in the universe without any mechanism to dilute the gravitino relic abundance, we argue that the entropy production is requisite for m 3/2 ∈ [100eV, 1keV]. We derive the necessary amount of entropy production for each m 3/2 , and then argue that this entropy production is inevitably imprinted on the primordial (inflationary) gravitational wave (pGW) spectrum [8][9][10][11][12]. This means that the feature in the pGW spectrum works as a smoking gun for PeV scale SUSYbreaking scenarios. Since the entropy production must occur after gravitino decoupling, the corresponding frequency of the pGW falls in the best frequency band explored by future experiments. Even if they detect only the pGW spectrum without any suppression feature in the frequency range O(10 −10 )Hz f pGW O(10 −5 )Hz, our study contributes to ruling out all the PeV scale SUSY-breaking scenarios.
The outline of this paper is as follows. In Sec. II, we raise questions we address in this work by making a brief review for m 3/2 = O(100)eV gravitino cosmology. Then in Sec. III, we discuss how to obtain the maximally allowed f 3/2,max for each m 3/2 and discuss the result of the analysis. In Sec. IV, we discuss inevitable signatures on the pGW spectrum with PeV scale SUSY-breaking. We study the deviation of the pGW spectrum from its usual, almost scale-invariant form induced by the presence of an early matter dominated (EMD) era followed by the entropy production via the decay of a heavy degree of freedom. Finally our conclusion is made in Sec. V.
This gravitino mass range is of particular interest in that the relic abundance of the gravitinos is given almost in a model-independent way. On the spontaneous SUSY-breaking, gravitinos become massive by absorbing the Goldstino field. As far as sparticles are present in the MSSM thermal bath, gravitinos remain in equilibrium with the thermal bath. Afterwards, the thermal bath continues to cool down until its temperature (T ) reaches the following gravitino decoupling temperature where g * ρ (T ) is the effective number of relativistic degrees of freedom for the energy density in the MSSM thermal bath at temperature T , and mg is the gluino mass.
For F 10 11 − 10 12 GeV 2 , the low scale gaugemediated SUSY-breaking (GMSB) can explain soft masses of sfermions and gauginos. The gluino mass is dominantly generated by the loop correction contributed from colored messengers, which reads where N mess is the number of messengers, g c is the gauge coupling of the MSSM SU (3) c color gauge group, y is the coupling constant for the interaction between messenger fields and SUSY-breaking field, and M mess is the mass of the messenger. For a perturbative gauge mediation model, in order for the SUSY-breaking vacuum to be stable to date, the condition M 2 mess > > yF needs to be satisfied [13]. This implies that mg can be at most O(10)TeV for m 3/2 ∈ [100eV, 1keV]. Thus, we see that T 3/2,dec mg holds true in Eq. (2). For a given value of m 3/2 , T 3/2,dec thus obtained leads to the following estimate of the relic abundance of gravitinos where Ω 3/2 is the ratio of the present gravitino energy density to the critical energy, h parametrizes the Hubble expansion rate via H 0 = 100hkm/Mpc/sec, g * s (T ) is the effective number of relativistic degrees of freedom for the entropy density in the MSSM thermal bath at temperature T , and T 3/2,0 and T ν,0 are the temperature of gravitinos and neutrinos at present, respectively. Since g * s (T 3/2,dec ) depends on a model-dependent mass spectrum, we cannot obtain a one-to-one correspondence between ω 3/2 and m 3/2 from Eq. (4). However, because g * s in the MSSM is at most g * s 230, 3 it is fair to state that for each m 3/2 the minimum inevitable value of Ω 3/2 h 2 is given by ω 3/2 ω 3/2,min = 10.75 230 Given ω 3/2,min in Eq. (5), we notice that for the Universe with m 3/2 ∈ [100eV, 1keV], as the lightest SUSY particle (LSP), the gravitinos must exist today in the form of DM. Depending on its mass, the gravitino's relic abundance can explain a fraction of DM today or even exceed the DM relic abundance in the absence of any dilution mechanism. Making the estimate of the free-streaming length (λ FS ) of gravitinos via [14], Mpc , (6) it is realized that sub-keV gravitinos with g * s 230 serve as WDM owing to Therefore, if our Universe went through the spontaneous SUSY-breaking at PeV scale, it becomes inevitable today for the DM population to be composed of a CDM candidate and gravitino WDM (mixed DM scenario) or fully of gravitino WDM with a certain mechanism to dilute the relic abundance. If so, one important question to be addressed is whether the presence of gravitino WDM is consistent with various cosmological and astrophysical observables or not.
Regarding this question, for the later possibility, already the mass constraints on WDM from Lyman-α forest observation m thermal wdm > 5.3 keV [4] and from the redshifted 21cm signals in EDGES observations m thermal wdm > 6.1 keV [5,6] require the dilution of the energy density. This implies that for m 3/2 ∈ [100eV, 1keV], gravitinos can reside in the present Universe only as a fraction of the DM population.
Along this line of reasoning, we ask two key questions that have a direct connection to cosmological consequences of PeV SUSY-breaking scenarios and provide us with a powerful cosmological probe of PeV SUSYbreaking scenarios. These questions are 1. For each m 3/2 , what is the minimum amount of the entropy production to make the scenario consistent with observations?
2. How to probe and confirm such an entropy production?
In Sec. III and Sec. IV, we answer these questions by invoking the estimate of the number of satellite galaxies in the Milky Way (MW) and pGWs. In Sec. III, we quantify the necessary amount of the entropy production by ∆ ≡ S/S where S and S are the total entropy before and after the decay of a heavy particle X, respectively, which we assume to take place prior to Big Bang Nucleosynthesis (BBN) era. 4 To this end, we use the observed number of the satellite galaxies N sat in the MW to obtain the maximal allowed fraction of gravitino for each m 3/2 . This eventually provides us with the minimally required dilution factor ∆ min . In Sec. IV, we study the imprint of the EMD era and the entropy production on the pGW spectrum. Regarding testability, we discuss how future GW detection experiments can be employed to either confirm or rule-out PeV-scale SUSY breaking scenarios.

III. ∆min FROM Nsat
The mixed DM model (MDM), in which the present DM population consists of both CDM and WDM, is parametrized by two quantities: the WDM mass and the fraction of the DM relic abundance today attributable to WDM. For our case with gravitinos being WDM, these are denoted by (m 3/2 , f 3/2 ). For a given set of (m 3/2 , f 3/2 ), the matter power spectrum in MDM scenario (P MDM (k)) can be parametrized as where T (k, m 3/2 , f 3/2 ) is a transfer function.
Relating the matter power spectrum in the CDM case (P CDM (k)) to the one in the MDM case (P MDM (k)), the transfer function contains information for a mass and WDM fraction in a MDM scenario. Suppressing P CDM (k) at small scales, the transfer function is characterized by two quantities (m 3/2 , f 3/2 ). Its structure is determined by (1) a comoving wavenumber k sup beyond which the transfer function starts to deviate from the unity and (2) the depth of suppression T plateau at small scales (for k > k sup ) [18]. k sup is parametrized by the temperature ratio in Eq. (4) and m 3/2 while T plateau is done by f 3/2 . This can be better understood via Fig. 1. By comparing the red solid line and the yellow dashed line, one can see that a larger f 3/2 induces a greater suppression. Moreover, comparison between the red solid line and the blue dotted line shows that k sup becomes smaller for the smaller m 3/2 .
We compute P MDM (k) for a given set of (m 3/2 , f 3/2 ) by using the Boltzmann solver CLASS [19]. For our purpose, we rely on variation of the parameters in ncdm sector of CLASS. For a given m 3/2 , by varying Ω CDM h 2 and Ω ncdm h 2 , we set f 3/2 which we aim for. Also for the case with g * s = 230, we set a present gravitino temperature to be the one deduced from f 3/2 and Eq. (4) via the parameter T ncdm in CLASS.
production takes place before the temperature drops down to 10MeV. As explained later, we consider entropy production from the decay of a heavy particle X whose decay rate determines the decay time (temperature) via Γ X H. Possible candidates of X include particles in the messenger sector in GMSB scenarios [15], particles in the SUSY-breaking sector [16], or the lightest righthanded sneutrino [17]. As one of the ways to constrain (m 3/2 , f 3/2 ), we attend to the estimate of the number of dwarf satellite galaxies in the Milky Way. The approach we adopt in this paper is based on the one given in Ref. [20] (see also Refs. [21][22][23][24][25][26][27][28][29]): Fifteen satellite galaxies were observed by SDSS (Sloan Digital Sky Survey) with the sky coverage f sky 0.28. When this limited sky coverage and eleven classically known satellites are taken into account together, N sat 63 is inferred as the total number of satellites of the Milky Way. Considering the possibility that more satellites will be found in the future surveys, we take N sat 63 as a lower bound for the number of the satellite galaxies that any DM model should satisfy. This set-up gives us a upper bound on f 3/2 for each m 3/2 below which the presence of gravitinos is consistent with the number of satellite galaxies of the Milky Way.
Given P MDM (k), one can make the estimate of the expected number of dwarf satellite galaxies residing in a host halo. Our estimate closely follows Refs. [25,26], which are based on the extended Press-Schechter approach [30,31] with the conditional mass function [32]. Adopting a sharp-k filter for the window function, it is calculated by [25,26] where M i , R i , and S i are the mass, the filter scale, and the variance of the satellite galaxy (for i = s) or of the host halo (for i = h), respectively. While the relation between the mass M i and the filter scale R i is in principle unconstrained in the sharp-k modeling, we follow Ref. [25] and adopt M i = (4π/3) × (cR i ) 3 × Ω m × ρ cr,0 with c = 2.5, which matches with observations best. Here ρ cr,0 is the critical density. The overall normalization is taken to be C n = 45 to reproduce the N-body simulation result [33]. In addition, we take M min = 10 8 h −1 M as the minimum mass of the dwarf satellite galaxies [34] and M h = 1.77 × 10 12 h −1 M as the Milky Way mass based on Ref. [35]. The variance S i is the function of R j and given by In the left panel of Fig. 2, using Eq. (5), we show for each m 3/2 (i) the minimum fraction of DM contributed by the gravitino WDM which is unavoidable in PeV SUSYbreaking scenarios in the absence of the energy dilution (red line), and (ii) the fraction of gravitino satisfying N sat = 63 (blue line). Below the blue line, the number of satellite galaxies is larger than 63. As one can see from the gap between the two lines, even the minimum predicted amount of gravitino exceeds the observationally allowed one. This implies that our Milky Way should be left with too few satellites if there is no any history of dilution of the gravitino energy density. We take this point, therefore, as a strong hint for the presence of an era when the entropy production is made via, for example, a heavy particle decay.
Given Fig. 1 and Eq. (8), now one can grasp in a qualitative manner how the parameters (m 3/2,f 3/2 ) have an effect on N sat . As seen in the integral range in Eq. (8), we consider the mass range 10 8 h −1 M − 10 12 h −1 M as the mass range of the satellite galaxies. This mass range corresponds to the range of the filter scale 0.03h −1 Mpc < R s < 0.66h −1 Mpc. For m 3/2 ∈ [100eV, 1keV], as can be seen in Fig. 1, the suppression in P MDM as compared to P CDM starts to occur for O(0. s , now these imply that suppression in P MDM induced by the presence of sub-keV gravitinos reduces N sat by making the integrand in Eq. (8) smaller than CDM case. Because the smaller m 3/2 and the large f 3/2 give rise to the larger suppression in P MDM , we expect that the criterion N sat > 63 yields more stringent constraint on f 3/2 for a smaller m 3/2 .
In the right panel of Fig. 2, we show the minimum necessary amount of entropy production for the scenario to be consistent with the observed number of the satellite galaxies of the Milky Way. For each m 3/2 , ∆ min = S/S is obtained by dividing f 3/2,min (red line) by f 3/2 associated with N sat = 63 (blue line) in the left panel of Fig. 2. Note that when the precise g * s (T 3/2,dec ) is taken into account, a larger ∆ would be required. For g * s (T 3/2,dec ) < 230, the red line moves upward while the blue line does downward, causing the green line to go up.
Before leaving Fig. 2, one may wonder whether the Lyman-α forest observations can offer a stronger constraint on f 3/2 than shown in the blue line in the left panel of Fig. 2. Indeed in Ref. [36], it is pointed out that the MDM model is very difficult to accommodate Lymanα data unless f 3/2 is tiny. Provided f 3/2 is more severely constrained by Lyman-α forest data, 5 apparently more 5 For instance, for a WDM with m wdm 700eV, it was argued in Ref. [37] that f wdm ∼ 0.15(∼ 0.1) can be consistent with the Lyman-α data (and higher-resolution data). S is the total entropy after a heavy particle decay (entropy production) whileS is the existing entropy before the entropy production took place. entropy production is required than the green line in the right panel of Fig. 2 for each m 3/2 . In any case, the green line remains as the most conservative lower bound for the necessary amount of ∆ to be consistent with a variety of cosmological observations. Given this inevitable requirement any PeV SUSYbreaking scenario confronts, the next important question is how one can confirm the presence of the EMD era and the sudden increase in the radiation energy density. In the next section, we address this question by studying the pGW spectrum. We will see that these non-standard histories leave their trace on the tensor modes produced from the quantum fluctuation during inflation and reentering the horizon before or during the entropy production.
We envision a scenario where a heavy particle X is present in the first radiation-dominated (RD1) era. X is assumed to be out-of-equilibrium from the thermal bath and thus its energy density scales as ρ X ∝ a −3 . When ρ X becomes comparable to the energy density of the existing thermal bath, an EMD era starts and continues until X decays to produce the additional entropy. After the decay, the second radiation-dominated (RD2) era gets started and continues until the matter-radiation equality at z eq 3300 is reached.

IV. PRIMORDIAL GRAVITATIONAL WAVES
Gravitational waves h ij (t, x) in the Friedmann-Robertson-Walker (FRW) background can be written as where τ is the conformal time defined via dt = adτ , and h ij satisfies the traceless and transverse conditions h ii = where λ = +, × are the GW polarizations,k is the unit vector along the three momentum k, and the polarization tensors λ ij are traceless and transverse as with h ij (t, x). The polarization tensors are normalized via λ ij ( λ ij ) * = 2δ λλ . The Fourier components h + k and h × k commonly obey the following time evolution equation in the absence of an anisotropic stress where the dot denotes the time derivative. After the relevant modes get sufficiently sub-horizon, the energy density of pGWs is given as where G (≡ (8πM 2 P ) −1 ) is the Newtonian constant and ... denotes the oscillation average. In this paper we assume the homogeneity and isotropy of the pGWs. As long as we average over a sufficiently large number of oscillations, we may add the spatial average to the definition of ... . Using the Fourier transform of h ij in Eq. (11), |∂h λ k (τ )/∂τ | 2 k 2 |h k (τ )| 2 for the sub-horizon k-modes and translating the spatial average into the ensemble average, one obtains the pGW energy density per logarithmic wavenumber where we define k = √ k · k. The label τ is understood as the time coordinate after taking oscillation average around it. The tensor power spectrum P T (k, τ ) is defined via the ensemble average as (16) which can be decomposed into the primordial part and the transfer function encoding the time evolution The transfer function is defined through h λ k (τ ) = h λ,prim k T T (k, τ ). For the primordial part, P prim T (k) is written in terms of the CMB pivot scale k CMB = 0.002Mpc −1 as where P prim R (k CMB ) 2.1 × 10 −9 is the amplitude of the scalar perturbation spectrum, r ≡ P prim is the tensor-to-scalar ratio and n T is the tensor spectral index. Note that the most recent bound on the tensor-to-scalar ratio is r 0.056 [38] which gives P prim T (k CMB ) 1.176 × 10 −10 . While P prim T (k) is determined by the initial conditions set during the inflationary era as in Eq. (18), T 2 T (k, τ ) reflects physical processes experienced by the k-mode after horizon re-entry. Even within the Standard Model, the pGW spectrum carries rich information on the particle content in/out of the thermal bath [39][40][41][42]. In BSM scenarios we expect much richer structures [43][44][45][46][47][48][49][50]. Similarly in the present scenario, the cosmological effects of the PeV scale SUSY-breaking is encoded in T 2 T (k, τ ) through the entropy production [51][52][53][54][55]. The decay of the pGW amplitude inversely proportional to the scale factor after horizon re-entry implies T 2 T (k, τ ) ∝ a −2 . For the modes re-entering the horizon during deep in either RD1 or RD2 era, we can write the transfer function as [44] where the subscript k denotes the time of re-entry k = a k H k , and the factor 1/2 arises from the oscillation average deep inside the horizon. Thus we obtain the following pGW spectrum expression that holds for these modes We further use the Friedmann equation 3M 2 P H 2 k = ρ rad = g * ρ,k (π 2 /30)T 4 k at the time of horizon re-entry together 6 Depending on the equation of state at the time of horizon reentry, an extra factor appears in Eq. (19) due to the non-zero "thickness" of the re-entry [44]. For the modes re-entering during a RD era, this extra factor takes unity and Eq. (19) holds true.
with the entropy relation g * s,k a 3 k T 3 k = ∆ −1 × g * s,0 a 3 0 T 3 0 (for the modes re-entering during RD1) or g * s,0 a 3 0 T 3 0 (RD2) to obtain the present pGW spectrum Ω GW,0 (k) = Ω rad,0 24 g * ρ,k g * ρ,0 g * s,0 g * s,k Here Ω rad,0 ≡ ρ rad,0 /ρ cr,0 = 4.2 × 10 −5 h −2 is the radiation energy fraction today, and g * ρ and g * s account for the relativistic degrees of freedom for the energy and entropy density, respectively. We use g * ρ,0 = 3.383 and g * s,0 = 3.931 for their present values. While Eq. (21) already tells us the rough behavior of the pGW spectrum at low and high wavenumbers, we parametrize it as 7 Ω GW,0 (k) = Ω rad,0 24 in order to account for the wavenumbers re-entering the horizon during the entropy production.
To obtain C ∆ (k), we explicitly calculate the time evolution of the pGW field. The final spectral shape depends on two parameters: the temperature T ∆ at which the entropy production occurs, and the dilution factor ∆ ≡ S/S. Through H(T ∆ ) = Γ X , the first parameter can be exchanged with the decay rate of the heavy particle (Γ X ). We obtain the value of the second parameter once the value of m 3/2 is specified based on Fig. 2. Different combinations of the two parameters (T ∆ , ∆) give rise to different values of C ∆ (k). This difference in C ∆ (k) enables us to probe the early Universe history through Ω GW (f GW ) with f GW = k/(2πa 0 ).
Before the heavy particle decays, the energy density of the Universe is mainly contributed from radiation (ρ rad ) and the heavy particle X (ρ X ). Thus the time evolution of the Hubble expansion rate in Eq. (13) is contributed by ρ rad and ρ X througḣ When solved together with Eq. (23), Eq. (13) yields h + k (τ ) and h × k (τ ) which in turn produce C ∆ (k). In the left panel of Fig. 3, we show C ∆ (k) obtained by numerically solving Eq. (13) for different choices of the dilution factor ∆. Each different color corresponds to the specified dilution factor. The suppression in C ∆ (k) is observed for the mode k 0.1k dec where k dec = a dec H dec is the mode that re-enters the horizon at the time when the decay of the heavy particle X takes place H dec = Γ X . Note that the value of the dilution factor is reflected in the ratio of the two plateaus which is equal to ∆ −4/3 .
Concerning g * ρ and g * s , since the MSSM sparticle mass spectrum depends on the details for how soft masses are generated, i.e. model-dependent, we do not study how GWs reflect the changes in these two for the temperature T T 3/2,dec = O(1)TeV. Instead, we focus on those modes that re-enter the horizon for T MSSM T 3/2,dec = O(1)TeV. Thus, for the temperature of the Universe of our interest, effectively only the SM particles and gravitinos exist. For the time evolution of g * ρ and g * s contributed by the SM particles, we shall refer to the tabulated data given in Ref. [42]. For the gravitino contribution, we evaluate the following: where we defined x 3/2 = m 3/2 /T and we assume the chemical potential of the gravitino to be negligible. Adding Eqs. (24) and (25) to g * ρ and g * s for the SM, we obtain the final g * ρ and g * s .
We show in the right panel of Fig. 3 the pGW spectrum Ω GW (f GW ) computed from Eq. (22). The pGW frequency and the temperature T hc when the relevant mode re-enters the horizon are related via To see how Ω GW depends on the inflation model, we assume the standard single-field slow-roll inflation of n T = −r max /8 = −0.007 with r max the maximum allowed tensor-to-scalar ratio, and a non-minimal inflation model with n T = 0.4 [56][57][58][59][60][61]. 8 As an example, we take T ∆ = 100MeV with ∆ = 0, 5, 10 and 20. We observe that Ω GW with ∆ = 0 is characterized by the suppression at f GW,∆ corresponding to T ∆ (related via Eq. (26)), which can be used to investigate whether the Universe went through the entropy production during the RD era. Larger ∆ induces a greater suppression in Ω GW for f GW > f GW,∆ . The underlying reason for the suppression is the presence of the EMD era where Ω GW ∝ a −1 . Because Ω GW ∝ a 0 holds true during the RD2 era, given the almost scale-invariant initial spectra on the super-horizon scale, dilution of Ω GW during the EMD era causes the suppression in the spectrum for the modes k k dec . Although only the case with T ∆ = 100MeV is shown in the right panel of Fig. 3, for other choices of T ∆ we expect a similar spectral shape with suppression at the relevant f GW,∆ .
It should be noted that T ∆ is limited to the range 10MeV T ∆ O(1)TeV. The lower bound is to ensure the successful BBN [64,65] (see also references therein). On the other hand, as discussed in Sec. II, T 3/2,dec can be at most O(10)TeV since mg O(10)TeV. Since the entropy production should happen after gravitinos decouple from the thermal bath, we conclude that T ∆ O(1)TeV should be the case. Hence, for a PeV SUSY-breaking scenario, we can make a quite solid argument that the suppression of Ω GW due to the entropy production should occur at f GW,∆ = O(10 −10 )Hz−O(10 −5 )Hz, irrespective of the inflation model assumed.
In the right panel of Fig. 3, given the green shaded region showing the parameter space potentially probed by SKA [66][67][68], we notice that SKA has a chance to see the suppression directly for T ∆ = O(100)MeV when n T = 0.4. However, for other cases with much different values of T ∆ or with n T = −r/8, direct observation of the suppression in Ω GW with SKA may not be possible. In this regard, the detection of CMB B-mode polarization in the future can play a critical role, even if the relevant frequency is too low to directly probe f GW,∆ of our interest. A careful comparison between the amount of pGWs observed at the CMB scales and at high frequencies by space-based interferometers such as LISA [69,70], DE-CIGO [71,72], Taiji [73], TianQin [74,75], and BBO [76][77][78], may indirectly reveal the presence of the suppression in Ω GW . Although this way of indirect investigation cannot pin down T ∆ , the entropy production ∆ can at least be determined. Therefore, the synergy between CMB B-mode polarization surveys and future space-based interferometers still renders PeV SUSY-breaking scenarios testable even if future pulsar timing arrays (such as IPTA [79][80][81] and SKA) are not sensitive enough to probe the entropy production at 10MeV T ∆ O(1)TeV. 9 We end this section by commenting on effects that sub-keV gravitinos can have on the pGW spectrum as freestreaming dark radiation (DR) [82][83][84]. 10 Once gravitinos decouple from the MSSM thermal bath at T 3/2,dec = O(1)TeV, they start to behave as DR. Without a significant late time entropy production after decoupling, ∆N eff contributed by gravitinos (say N 3/2 ) is given by 9 Note that at the GW frequency range relevant to DECIGO and BBO, some or all of the MSSM particles are relativistic. Thus, when probing the suppression of Ω GW due to the entropy production using DECIGO and BBO, the additional suppression caused by the larger g * ρ(T in ) and g * s(T in ) in Eq. (22) should be taken into account. 10 See [85] for another interesting effect of free-streaming radiations on GWs. It is an effect on the IR "causal" part of the spectrum in late time GW production.
where T ν,dec is the neutrino decoupling temperature and g * s (T ν,dec ) = 10.75 is the effective number of degrees of freedom for the entropy density evaluated at T ν,dec . If DR makes a significant contribution to ∆N eff , there can be an overall enhancement of the pGW spectrum (due to the change in the expansion rate and that in the matterradiation equality) and a suppression feature from the non-zero anisotropic stress induced by the free-streaming of DR [44,46]. However, since N 3/2 is too small for the scenarios discussed in the present paper, these effects from N 3/2 are too small to be visible in the right panel of Fig. 3.

V. CONCLUSION
In this work, we demonstrated the necessity of the late time entropy production prior to BBN era for the Universe with PeV scale SUSY-breaking. Our argument is based on the observation that the theoretically predicted relic abundance of gravitino WDM with m 3/2 ∈ [100eV, 1keV] is too large to be consistent with the observed number of satellite galaxies N sat in the Milky Way, provided that there is no energy dilution mechanism. By requiring N sat 63 in the current Universe in which sub-keV gravitino WDM is responsible for a fraction of the DM population, we quantified the minimum amount of the dilution factor ∆ min necessary for observational consistency.
Based on this result, we proposed using the pGW spectrum Ω GW to investigate the presence of an early matterdominated era and entropy production prior to the BBN era. The degree and the characteristic frequency of the suppression in the pGW spectrum, once observed by future pulsar timing arrays such as SKA, give us the information on the entropy production ∆ and the temperature T ∆ when it occurs. Interestingly, the suppression is expected at O(10 −10 )Hz f GW O(10 −5 )Hz in PeV scale SUSY-breaking scenarios, making these scenarios testable via pGW searches. Even if the pGW abundance is too small to be probed by SKA, the synergy between future CMB B-mode experiments and space-based interferometers can alternatively probe the suppression in the pGW spectrum induced by ∆ = 0. If the absence of any suppression feature in the pGW spectrum is established by the comparison between the two frequency regimes f GW < O(10 −10 )Hz and f GW > O(10 −5 )Hz, our study immediately rules out PeV scale SUSY-breaking scenarios.
We have concentrated on the case of m 3/2 = O(100)eV in this paper. However, we stress that our method of using pGWs can be extended to test the low scale SUSY with larger mass range of the gravitino such as 1keV − 1GeV which is predicted in gauge mediation models [86,87].