PeV-Scale SUSY and Cosmic Strings from F-term Hybrid Inflation

We consider F-term hybrid inflation (FHI) and SUSY breaking in the context of a B-L extension of MSSM which largely respects a global U(1) R symmetry. The hidden sector Kaehler manifold enjoys an enhanced SU(1,1)/U(1) symmetry with the scalar curvature determined by the achievement of a SUSY-breaking de Sitter vacuum without ugly tuning. FHI turns out to be consistent with data, provided that the magnitude of the emergent soft tadpole term is confined in the range (1.2-100) TeV, and it is accompanied with the production of B-L cosmic strings. If these are metastable they interpret the present observations from PTA experiments on the stochastic background of gravitational waves with dimensionless tension Gmu~(1-9.2)x10^-8. The mu parameter of MSSM arises by appropriately adapting the Giudice-Masiero mechanism and facilitates the out-of-equilibrium decay of the R saxion at a reheat temperature lower than about 71 GeV. Due to the prolonged matter dominated era the gravitational wave signal is suppressed at high frequencies. The SUSY mass scale turns out to lie in the PeV region.


INTRODUCTION
Supersymmetric (SUSY) hybrid inflation [1] based on F terms, called for short henceforth F-term hybrid inflation (FHI) [2] is undoubtedly a well-motivated inflationary model -for reviews see Ref. [3].The most notable reasons which support our statement above are the following: • FHI is tied to a renormalizable superpotential uniquely determined by a gauge and a global U (1) R symmetries; • FHI does not require fine-tuned parameters and transplanckian inflaton values; • FHI can be reconciled with the Planck data [4] -fitted to the standard power-law cosmological model with Cold Dark Matter (CDM) and a cosmological constant (ΛCDM) -if we take properly into account not only radiative corrections (RCs) but also corrections originated by supergravity (SUGRA) [5][6][7][8] as well as soft SUSY-breaking terms [9][10][11][12][13].
• FHI can be naturally followed by a Grand Unified Theory (GUT) phase transition which may lead to the production of cosmological defects, if these are predicted by the symmetry breaking scheme.In the large majority of GUT breaking chains the formation [14] of cosmic strings (CSs) cannot be avoided.
Although the last feature above is often used to criticize the powerfulness of the embedding of FHI in several GUTs -see e.g.Ref. [15] -, it recently appears as an interesting ingredient in the cosmological model building.This is, because the announced data from several pulsar timing array (PTA) experiments [16] -most notably the NANOGrav 15-years results (NANOGrav15) [17] -provide strong support for the discovery of a gravitational wave (GW) background around the nanohertz frequencies.
Given that the interpretation of this signal in terms of the merging of supermassive black hole binaries is somewhat disfavored [18], its attribution to gravitational radiation emitted by topologically unstable superheavy CSs -which may arise after the end of FHI -attracts a fair amount of attention [19][20][21][22][23][24].In particular, the observations can be interpreted if the CSs are meta- [25] or quasi-stable [26].Both types of topologically unstable CSs arise from the symmetry breaking G → H × U (1), which produces monopoles.The subsequent U (1) breaking yields CSs which connect monopoles with antimonopoles.We here assume that these last ones are inflated away, but can appear on metastable CSs via quantum pair creation.Therefore, metastable CSs can be easily achieved if a U (1) symmetry is embedded in a gauge group with higher rank such as the Pati-Salam [22], the flipped SU (5) [21] or SO (10) [20].For this reason, in this work we focus on FHI realized in a B − L extension of Minimal SUSY Standard Model (MSSM) -cf.Ref. [22,27] -which dilutes possibly preexisting monopoles and is naturally accompanied by the production of a network of CSs.On the other hand, we do not specify the mechanism of the monopole production and the metastability of CSs, as done, e.g., in Ref. [19,24].
We adopt the minimal possible framework [10,11] which supports observationally acceptable FHI.It employs minimal Kähler potential for the inflaton field, RCs and soft SUSY-breaking terms.From the last ones, the tadpole term plays a crucial role in establishing compatibility with data.Its magnitude can be motivated by intertwining [28] the inflationary sector (IS) with a hidden sector (HS) introduced in Ref. [29].Contrary to earlier attempts [30][31][32][33][34] this HS respects a mildly violated R symmetry which is compatible with that adopted for FHI [1].The consequences of the interconnection of the two sectors above are the following: • The R charge 2/ν of the goldstino superfield -which is related to the geometry of the HS [29] is constrained to values with 0 < ν < 1.
• The SUSY breaking is achieved not only in a Minkowski vacuum -as in the cases of Ref. [30][31][32][33][34] -but also in a de Sitter (dS) one which allows us to control the notorious Dark Energy (DE) problem by mildly tuning a single superpotential parameter to a value of order 10 −12 .
• The sgoldstino is stabilized [31][32][33][34] to low values during FHI.This fact together with the selection of a minimal Kähler potential for the inflaton assists us to resolve the η-problem.Note that Kähler potentials inspired by the string theory are mainly employed in earlier works [30][31][32][33][34].
• The solution to the µ problem [35] of MSSM is achieved by suitably applying [29] the Giudice-Masiero mechanism [36,37].Contrary to similar attempts [30,38,39] the µ term here plays no role during FHI but crucially controls the timely decay of the sgoldstino (or R saxion).
• The energy density of the universe is dominated by the energy density of sgoldstino condensate which decays [40][41][42][43][44] before the onset of the Big Bang Nucleosynthesis (BBN) at cosmic temperature (2 − 4) MeV [45] thanks the aforementioned µ term.Therefore, our scenario naturally motivates a prolonged matter domination which causes a reduction [46,47] of the spectrum of GWs at high frequencies (f > 0.1 Hz).This fact is welcome since it assists us to avoid any conflict with the third run advanced LIGO-VIRGO-KAGRA (LVK) data [48].
• The SUSY mass scale m is predicted to be close to the PeV scale [49].It fits well with the Higgs boson mass, discovered at LHC, as it is estimated [50] within high-scale SUSY if we assume a relatively low tan β and stop mixing.Note that the connection of inflation with SUSY breaking has been extensively discussed in literature [51] over the last years.
In this feature paper we review further the model introduced in Ref. [28] focusing exclusively on its implementation within a version of MSSM endowed by a B − L Higgs sector.We present the complete particle content of the model, paying special attention on the computation of the GWs emerging from the CS decay, under the assumption that those are metastable.We also explain the generation of neutrino masses taking into account SUGRA contributions from Ref. [52].
The remaining manuscript is built up as follows: In Sec. 2 we introduce our framework.Following, we revise the salient features of our model as regards its vacuum in Sec. 3 and the inflationary era in Sec. 4.Then, we study the reheating process in Sec. 5 and the production of the GWs from the CSs in Sec. 6.Our predictions for the SUSY-mass scale are exposed in Sec. 7. Our conclusions are summarized in Sec. 8. Details for the derivation of neutrino masses are given in Appendix A. Lastly, the acronyms adopted in the text are arranged in Appendix B.

MODEL SET-UP
We focus on an extension of MSSM invariant under the gauge group where G SM is the Standard Model gauge group.The charge assignments under these symmetries of the various matter and Higgs superfields are listed in Table 1.Namely, the ith generation SU (2) L doublet left-handed quark and lepton superfields are denoted by Q i and L i respectively, whereas the SU (2) L singlet antiquark [antilepton] superfields by u c i and d i c [e c i and N c i ] respectively.The electroweak Higgs superfields which couple to the up [down] quark superfields are denoted by H u [H d ].Besides the MSSM particle content, the model is augmented by seven superfields: a gauge singlet S, three N c i 's, a pair of Higgs superfields Φ and Φ which break U (1) B−L and the goldstino superfield Z.In addition to the local symmetry, the model possesses also the baryon and lepton number symmetries and an R symmetry U (1) R .The latter plays a crucial role to the construction of the superpotential (see Sec. 2.1) and the Kähler potential (see Sec. 2.2).

SUPERPOTENTIAL
The superpotential of our model respects totally the symmetries in Table 1.Most notably, it carries R charge 2 and is linear with respect to (w.r.t.) S and Z ν .It naturally splits into five parts: where the subscripts "I" and "H" stand for the IS and HS respectively and the content of each term is specified as follows: where κ and M are free parameters which may be made positive by field redefinitions.
(b) W H is the HS part of W written as [29] W H = mm 2 P (Z/m P ) ν . (2.2b) Here m P = 2.4 ReV is the reduced Planck mass -with ReV = 10 18 GeV -, m is a positive free parameter with mass dimensions, and ν is an exponent which may, in principle, acquire any real value if W H is considered as an effective superpotential valid close to the non-zero vacuum expectation value (v.e.v) of Z, Z .We assume though that the effective superpotential is such that only positive powers of Z appear.If we also assume that W is holomorphic in S then mixed terms of the form S νs Z νz can be forbidden in W since the exponent of a such term has to obey the relation leading to negative values of ν z .This conclusion contradicts with our assumptions above.
(c) W GH is a term which mixes the HS and the B − L gauge fields of the IS.It has the form with λ a real coupling constant.The magnitude of λ can be restricted by the DE requirement as we see below.
(d) W MSSM is the part of W which contains the usual trilinear terms of MSSM, i.e., (2.2d) The selected R assignments in Table 1 prohibit the presence in W MSSM of the bilinear µH u H d term of MSSM and other mixing terms -e.g.λ µ SH u H d [39] which is frequently employed to generate this µ term.This term is generated here via K µ -see Sec.2.2 below.
(e) W MD is the part of W which is provides masses to neutrinos The first term in the right-hand side of Eq. (2.2e) is responsible for Dirac neutrino masses whereas for the Majorana masses -cf.Refs.[10,15].The scale of the latter masses is intermediate since Φ ∼ 1 YeV and Z ∼ m P .The cooperation of both terms lead to the light neutrino masses via the well-known (type I) seesaw mechanism -see also Appendix A.

K ÄHLER POTENTIAL
The Kähler potential respects the G B−L , B and L symmetries in Table 1.It has the following contributions which left-handed chiral superfields of MSSM denoted by Y α with α = 1, ..., 7, i.e., where the generation indices are suppressed.
(a) K I is the part of K which depends on the fields involved in FHI -cf.Eq. (2.2a).We adopt the simplest possible choice -cf.Ref. [2,10] -which has the form (2.4a) Higher order terms of the form |S| 2ν S /m 2ν S −2 P with ν S > 1 can not be excluded by the imposed symmetries but may become harmless if S ≪ m P and assume low enough coefficients.
(b) K H is the part of K devoted to the HS.We adopt the form introduced in Ref. [29] where Here, k > 0 mildly violates R symmetry endowing R axion with phenomenologically acceptable mass.
The selected K H is not motivated by the string theory but it can be considered as an interesting phenomenological option for two reasons: It largely respects the R symmetry, which is a crucial ingredient for FHI, and it ensures -as we see in Sec. 3 -a dS vacuum of the whole field system with tunable cosmological constant for for N < 0. (2.5) Our favored ν range finally is 3/4 < ν < 1.Since N < 0, K H parameterizes the SU (1, 1)/U (1) hyperbolic Kähler manifold for k ∼ 0.
(c) K µ includes higher order terms which generate the needed mixing term between H u and H d in the lagrangian of MSSM -cf.Ref. [29,36] -and has the form where the dimensionless constant λ µ is taken real for simplicity.
(d) K D is an unavoidable term which mixes the observable sector with the HS as follows It provides (subdominant) Dirac masses for ν i [52] as shown in Appendix A.
The total K in Eq. (2.3) enjoys an enhanced symmetry for the Y α , S and Z fields, namely where the indices indicate the moduli which parameterize the corresponding manifolds.Thanks to this symmetry, mixing terms of the form S νs Z * νz can be ignored although they may be allowed by the R symmetry for νz = ν νs .Most notably, U (1) S protects K I from S depended terms which violates the R symmetry, thereby, spoiling the inflationary set-up.

SUSY AND G B−L BREAKING -DARK ENERGY
The vacuum of our model is determined by minimizing the F-term (tree level) SUGRA scalar potential V F derived [28] from W in Eq. (2.1) and K in Eq. (2.3).Note that D-term contributions to the total SUGRA scalar potential vanish if we confine ourselves during FHI and at the vacuum along the D-flat direction Here g is the unique (considered unified) gauge coupling constant of G B−L and we employ the same symbol for the various superfields X α = S, Z, Φ, Φ and their complex scalar components.
As we can verify numerically, V F is minimized at the G B−L -breaking vacuum It has also a stable valley along θ = 0 and θ S /m P = π, where these fields are defined by Substituting Eq. (3.3) in V F and minimizing it w.r.t the various directions we arrive at the results which yield the constant potential energy density Tuning λ to a value λ ∼ m/m P ≃ 10 −12 we may wish identify V F with the DE energy density, i.e., where the density parameter of DE and the current critical energy density of the universe are respectively given by [58] Ω Λ = 0.6889 and ρ c = 2.31 • 10 −120 h 2 m 4 P with h = 0.6766.
The particle spectrum of the theory at the vacuum in Eqs.(3.2) and (3.4) includes the gravitino ( G) which acquires mass [29] Diagonalizing the the mass-squared matrix of the field system S − Φ − Φ − Z we also find out that the IS acquire a common mass where the second term arises due to the coexistence of the IS with the HS -cf.Ref. [10].Similar mixing does not appear in the mass spectrum of HS which contains the (canonically normalized) sgoldstino (or R saxion) and the pseudo-sgoldstino (or R axion) with respective masses Applying finally the relevant formulas of Refs.[29,37], we find that K µ induces a non-vanishing µ term in the superpotential of MSSM whereas W MSSM and K µ + |Y α | 2 in Eqs.(2.2d) and (2.3) lead to a common soft SUSY-breaking mass parameter m, at the vacuum of Eq. (3.4), which indicatively represents the mass level of the SUSY partners.Namely, we obtain Variant form of the terms |Y α | in K -see Eq. ( 2.3) -do not alter essentially our results [28,29].

INFLATION ANALYSIS
It is well known [1,3] that in global SUSY FHI takes place for sufficiently large |S| values along a F-and D-flat direction of the SUSY potential are the constant potential energy density and correspoding Hubble parameter which drive FHI -the subscript 0 means that this is the tree level value.In a SUGRA context, though, we first check -in Sec.4.1 -the conditions under which such a scheme can be achieved and then in Sec.4.2 we give the final form of the inflationary potential.Lastly, we present our results in Sec.4.4 imposing a number of contraints listed in Sec.4.3.

HIDDEN SECTOR'S STABILIZATION
The attainment of FHI is possible if Z is well stabilized during it.The relative mechanism is pretty well known [34].Z is transported due to V I0 form its value in Eq. (3.2) to a value well below m P .To determine this, we construct the complete expression for V F along the inflationary trajectory in Eq. (4.1) and then expand the resulting expression for low S/m P values, assuming that the θ = 0 direction is stable as in Eq. (3.2).Under these conditions V F is minimized for the value This result is in good agreement with its precise value derived numerically.Note that ν < 1 assures a real value of z I with z I ≪ m P since H I /m ≪ 1.
The (canonically normalized) components of sgoldstino, acquire masses squared, respectively, ) whereas the mass of G turns out to be It is clear from the results above that m Iz ≫ H I and therefore it is well stabilized during FHI whereas m Iθ ≃ H I and gets slightly increased as k increases.However, the isocurvature perturbation is expected to be quite suppressed since it becomes observationally dangerous only for m Iθ ≪ H I .

INFLATIONARY POTENTIAL
Expanding V F for low S values, introducing the canonically normalized inflaton σ = √ 2|S| and taking into account the RCs [1,3] we derive [28] the inflationary potential V I which can be cast in the form The individual contributions are specified as follows: (a) C RC represents the RCs to V I /V I0 which may be written as [1-3] , beyond which the expression above ceases to be valid.Here we take into account that the multiplicity of the waterfall fields is N G = 1 since these are U (1) B−L non-singlets.
(b) C SSB is the contribution to V I /V I0 from the soft SUSY-breaking effects [9] parameterized as follows: where the tadpole parameter reads The minus sign results from the stabilization of θ -see Eq. (3.3) -at zero and the minimization of the factor (S + S * ) = √ 2σ cos(θ S /m P ) which occurs for θ S /m P = π (mod 2π) -the decomposition of S is shown in Eq. (3.3).We further assume that θ S remains constant during FHI so that the simple onefield slow-roll approximation is valid.Possible variation of θ S is investigated in Ref. [11] where they found that acceptable solutions with θ S /m P = π require a significant amount of tuning.The first term in Eq. (4.5b) does not play any essential role in our set-up due to low enough m 3/2 's -cf.Ref. [10].
(c) C SUGRA is the SUGRA correction to V I /V I0 , after subtracting the one in C SSB , which is Thanks to the minimality of K I in Eq. (2.4a) and the smallness of z I the coefficients above are low enough and allow for FHI to be established -cf.Ref. [9,10].

OBSERVATIONAL REQUIREMENTS
The analysis of FHI can be performed in the slow-roll approximation if we calculate the slow-roll parameters [4] where the derivatives of the various contributions read The required behavior of V I in Eq. (4.4) can be obtained thanks to the relation C ′ RC ≃ −C ′ SSB which is established for carefully selecting κ (or M ) and a S .Apparently, we have On the contrary, C ′′ RC < 0, since the negative contribution 3x 2 ln(1 − 4/x 4 ) dominates over the first positive one, and so we obtain η < 0 giving rise to acceptably low n s values.
Our model of FHI can be qualified if we test it against a number of observational requirements.Namely: (a) The number of e-foldings elapsed between the horizon crossing of the pivot scale k ⋆ = 0.05/Mpc and the end of FHI has to be adequately large for the resolution of the horizon and flatness problems of standard Big Bang cosmology.Taking into account -see Sec. 5 -, that FHI is followed in turn, by a matter and a radiation dominated era, the relevant condition takes the form [2,4]: where the prime denotes derivation w.r.t.σ, σ ⋆ is the value of σ when k ⋆ crosses outside the horizon of FHI and σ f is the value of σ at the end of FHI.This is normally obtained by the critical point , the end of inflation coincides with the onset of the B − L phase transition.Note that σ ≃ 0, as mentioned below Eq. (3.6), and so it does not disturb the inflationary dynamics which governs the σ evolution for σ ≥ σ c .
(b) The amplitude A s of the power spectrum of the curvature perturbation generated by σ during FHI must be consistent with the data [58] on cosmic microwave background (CMB), i.e., (4.9) The observed curvature perturbation is generated wholly by σ since the other scalars are massive enough during FHI -see Sec.4.1.
(c) The remaining observables -the scalar spectral index n s , its running α s , and the scalar-totensor ratio r -which are calculated by the following standard formulas (where ξ ≃ m 4 P V ′ I V ′′′ I /V 2 I and all the variables with the subscript ⋆ are evaluated at σ = σ ⋆ ) must be in agreement with data.We take into account the latest Planck release 4 -including TT,TE,EE+lowE power spectra [58] -, Baryon Acoustic Oscillations (BAO), CMB-lensing and BICEP/Keck data [59].Fitting it with ΛCDM+r we obtain approximately n s = 0.965 ± 0.009 and r 0.032, ( at 95% confidence level (c.l.) with negligible |α s | ≪ 0.01.

RESULTS
As deduced from Sec. 4.1 -4.3, the inflationary part of our model depends on the parameters κ, M, m, λ, ν, k and λ µ .
Recall that N is related to ν via Eq.(2.5).Enforcing Eq. (3.6) fixes λ at a rather low value which does not influence our remaining results.Moreover, k affects exclusively m θ and m Iθ via Eqs.(3.8c) and (4.3a).We select throughout the value k = 0.1 which assures the avoidance of massless modes.We here present the a S values as a function of κ or M which assist the confrontation of FHI with data -cf.
Ref. [10] -and postpone their derivation from m and ν via Eq.(4.5c) in Sec. 7. As regards T rh , which influences Eq. (4.8), we adopt a value close to those met in our set-up, T rh ≃ 1 GeV.Enforcing Eqs.(4.8) and (4.9) we can restrict M and σ ⋆ as functions of our free parameters κ and a S .It is apparent that the allowed ranges of parameters is similar to those explored in Ref. [10], where the HS is not specified.Some deviations are only due to the improvements on the determination of the n s values in Eq. (4.11).The correct values of that quantity are attained if FHI becomes of hilltop type [10,11] .I.e., V I is non-monotonic and develops a maximum at σ max and a minimum at σ min ≫ σ max .For σ > σ max V I becomes a monotonically increasing function of σ and so its boundedness is assured.FHI takes place for σ < σ max .The position of σ max is predominantly dominated by C ′ RC and C ′ SSB .On the other hand, σ min can be roughly found by the interplay of C ′ SUGRA and C ′ SSB .Namely, we obtain [28]   Recall that 1 PeV = 10 6 GeV, 1 EeV = 10 9 GeV, 1 ZeV = 10 12 GeV and 1 YeV = 10 15 GeV.
we observe that, for constant κ and a S , n s decreases as σ ⋆ approaches σ max .To quantify somehow the amount of these tunings, we define the quantities The naturalness of the hilltop FHI increases with ∆ c⋆ and ∆ max * .
To get an impression of the required values of the parameters of the model, we construct first the benchmark table Table 2.We display there the results of our analysis for two benchmark points (BPs)  with κ = 0.0005 (BPA) and κ = 0.001 (BPB) and n s fixed to its central value in Eq. (4.11).We also employ some representative ν and k values.In all cases, we obtain N I⋆ ∼ 40 from Eq. (4.8).
From the observables listed in Table 2 we infer that |α s | turns out to be of order 10 −4 , whereas r is extremely tiny, of order 10 −11 , and therefore far outside the reach of the forthcoming experiments devoted to detect primordial GWs.From the entries of Table 2 related to ∆ c⋆ and ∆ max * , we notice that ∆ max * > ∆ c⋆ , their values may be up to 10%, and increase as a S or κ (and M ) increases.From the mass spectra arranged in Table 2 -see Sec.4.1 -we see that m I3/2 is similar to a S ∼ TeV whereas the other masses are of order EeV whereas at the vacuum m I turns out to be of order 1 ZeV = 10 12 GeV whereas m 3/2 , m z and m θ lie in the PeV range -see Sec. 3. In the same Table we find it convenient to accumulate the values of some parameters related to the reheating process, specified in Sec. 5 below, and the formed CSs -see Sec. 6.
To explore further the parameter space of our model allowed by the inflationary requirements in Eqs.(4.8), (4.9) and (4.11) we present the gray shaded regions in the κ − a S and M − a S plane -see Fig. 1-(a Within the margins above, ∆ c⋆ ranges between 0.5% and 9.5% and ∆ max * between 0.4% and 8.2% with the relevant values increasing with M and a S , as deduced from Table 2 too.The lower bounds of the inequalities above are expected to be displaced to slightly upper values due to the post-inflationary requirements -see Eq. (5.7) below -which are not considered here for the shake of generality.

REHEATING PROCESS
Soon after FHI, the IS and z enter into an oscillatory phase about their minima in Eqs.(3.2) and (3.4) and eventually decay via their coupling to lighter degrees of freedom.Note that θ and θ S remain well stabilized at their values shown below Eq. (3.2) during and after FHI and so they do not participate in the phase of damped oscillations.The commencement of z-dominated phase occurs for values of Hubble rate H zI ∼ m z .Given that z ∼ m P -see Eq. (3.4) -, the initial energy density of its oscillations ρ zI is comparable with the energy density of the universe ρ zIt at the onset of these oscillations since ρ zI ∼ m 2 z z 2 and ρ zIt = 3m 2 P H 2 zI ≃ 3m 2 P m 2 z . (5.1) Therefore, we expect that z will dominate the energy density of the universe until completing its decay through its weak gravitational interactions.Due to weakness of these interactions we expect that the reheating temperature T rh will be rather low.This is the notorious cosmic moduli problem [40,41] which plagues the vast majority of the SUGRA settings.
In our model though T rh adequately increases thanks to the large enough m z and µ originated by K µ in Eq. (2.6a).To show this fact we estimate T rh by [55] where g rh * ≃ 10.75 − 100 counts for the effective number of the relativistic degrees of freedom at T rh .This is achieved for where τ = ln (R/R I ) with R the scale factor of the universe and R I its value at the onset of the z oscillations.Also Γ δz is the total decay width of the (canonically normalized) sgoldstino which predominantly includes the contributions from its decay into pseudo-sgoldstinos, θ and electroweak higgs, H u and H d via the kinetic terms K XX * ∂ µ X∂ µ X * with X = Z, H u and H d [41][42][43][44] of the Lagrangian.In particular, we have where the individual decay widths are found to be From the expressions above we readily recognize that Γ δz is roughly proportional to m 3 z /m 2 P as expected for any typical modulus [40,41].We explicitly checked that z-decay channels into gauge bosons through anomalies and three-body MSSM (s)particles are subdominant.On the other hand, we kinematically block the decay of δz into G's selecting ν > 3/4 which ensures m z < 2m 3/2 -see Eq. (3.8c).We do so in order to protect our setting from the so-called [53] moduli-induced G problem, i.e., the possible late decay of the produced G and the problems with the abundance of the subsequently produced lightest SUSY particles -cf.Ref. [41]. ) n s = 0.965 Allowed strip in the κ − T rh plane compatible with the inflationary requirements in Sec.4.3 for n s = 0.965.We take ν = 7/8 and µ = m (solid line), µ = m/3 (dot-dashed line) or µ = 3 m (dashed line).The BBN lower bounds on T rh for hadronic branching ratios B h = 1 and 0.001 are also depicted by two thin lines.
If we enforce the compatibility between theoretical and observational values of light element abundances predicted by BBN we achieve Ref. [45] a lower bound on T rh which, for large m z ∼ 0.1 PeV, entails where B h is the hadronic branching ratio.The bound above is a little softened for larger m z values.Taking κ and m z values allowed by the inflationary part of our model in Sec.4.4, we evaluate T rh as a function of κ and delineate the regions allowed by the BBN constraints in Eq. (5.7) -see Sec.4.3 below.The results of a such computation are displayed in Fig. 2, where we design allowed contours in κ − T rh plane for ν = 7/8.This is an intermediate value in the selected here margin (3/4 − 1).The boundary curves of the allowed region correspond to µ = m/3 or λ µ = 0.22 (dotdashed line) and µ = 3 m or λ µ = 1.96 (dashed line) whereas the solid line is obtained for µ = m or λ µ = 0.65.Note that the relation between λ µ and µ/ m is given in Eq. (3.9).We see that there is an ample parameter space consistent with the BBN bounds in Eq. (5.7) depicted with two horizontal lines.Since the inflationary requirements increase the scale m with κ and m heavily influences m z and also T rh -see Eq. (5.2) -T rh increases with κ.The maximal value of T rh for the selected ν is obtained for µ = 3 m and is estimated to be Obviously, decreasing µ below m/3, causes λ µ , Γ δz and T rh to decrease too, and the slice cut from the BBN bound increases.
It is worth emphasizing that the reheating stage is not instantaneous.In particular, the maximal temperature T max during this period is much larger than T rh , which can be better considered as the largest temperature of the radiation domination [54].To be quantitative, T max can be calculated as [55] T max = (3/8) 3/20 12 1/4 3 1/8 m (5.9) It is achieved [55] for τ = τmax = 0.39 ≪ τrh .The large hierarchy between T max and T rh can be appreciated by their numerical values displayed in Table 2 for BPA and BPB.We see that T max ∼ 1 PeV whereas T rh ∼ 1 GeV.As a consequence, the electroweak sphalerons are still operative and so baryogenesis via leptogenesis is in principle feasible.To obtain a more complete picture of our post-inflationary cosmological setting we introduce the logarithmic time τ which is defined as a function of the redshift z (not to be confused with saxion field z) or the scale factor R as [57] where the subscript 0 is refereed hereafter to present-day values.T max is achieved at τ max which can be found from taking into account Eqs.(5.3) and (5.9).For τ ≥ τ rh , τ can be found using the entropy conservation through the relation where g s * is the entropy effective number of degrees of freedom at temperature T and T 0 = 2.35 • 10 −13 GeV.Its precise numerical values are evaluated by using the tables included in public packages [56] and assuming the particle spectrum of the MSSM.E.g, at BBN (T = 2 MeV) τ BBN ≃ −23 and at the time of radiation-matter equidensity τ eq ≃ −8.1 corresponding to z eq = 3387.The values of τ rh and τ max are accumulated in Table 2 for some sample T rh and ratios µ/ m.We remark that τ max values are similar for both T rh but τ rh is smaller for larger T rh .The energy density of the oscillating z condensate ρ z , radiation ρ r and matter ρ m may be evaluated as follows [55,57] where Ω m0 = 0.311 and the energy density at reheating is given by ρ rh = ρ r (T = T rh ).Taking advantage from the formulae above we illustrate in Fig. 3 our cosmological scenario.Namely, we present the evolution of the various log ρ i -with ρ i normalized to ρ c given in Eq. (3.7) -as functions of τ keeping, for simplicity, only the dominant component of the universal energy density for each τ.As regards log ρ i with i = z we employ the parameters of BPB in Table 2 besides µ which is µ = 5 m/4 for T rh = 1.5 GeV (dashed line) or µ = 2 m/5 for T rh = 0.5 GeV (dot-dashed line).We also draw with solid lines log ρ i with i = r, m.From the plot it is apparent that our scenario is distinguishable from the standard scenario according to which the radiation domination commences after a high T rh ∼ 1 EeV corresponding to τ rh ∼ −50.We observe that log ρ z and log ρ m have the same slop since both are proportional to −3τ whereas log ρ r is more stiff since it is proportional to −4τ.In the plot we also localize the position of τ max , which is the same for both T rh 's, the two τ rh values and the equidensity point τ eq .Shown is also the τ value which corresponds to the CS decay τ dc for r 1/2 ms = 8 -see Sec.6.2.It is located within the radiation dominated era.

METASTABLE CSS AND AN EARLY MATTER DOMINATION
The U (1) B−L breaking which occurs for σ ≃ σ c produces a network of CSs which may be stable meta-or quasi-stable.This network has the potential to undergo decay via the Schwinger production of monopole-antimonopole pairs leading thereby to the generation of a stochastic GW background.We below compute the tension of these CSs in Sec.6.1 and their GW spectrum in Sec.6.2 under the assumption that these are metastable.

CS TENSION
The dimensionless tension Gµ cs of the B − L CSs produced at the end of FHI can be estimated by [10,27] ǫ cs (r cs ) with ǫ cs (r cs ) = 2.4 ln(2/r cs ) and r cs = κ 2 /8g 2 ≤ 10 −2 , ( where we take into account that (B − L)(Φ) = 2 -cf.Ref. [28].Here G = 1/8πm 2 P is the Newton gravitational constant and g ≃ 0.7 is the gauge coupling constant at a scale close to M .For the parameters in Eq. (4.14) we find 0.59 Gµ cs /10 −8 9.2.( To qualify the result above, we single out the cases: (a) If the CSs are stable, the range in Eq. (6.2) is acceptable by the level of the CS contribution to the observed anisotropies of CMB which is confined by Planck [60] in the range On the other hand, the results of Eq. (6.2) are completely excluded by the recent PTA bound which requires [18] Gµ cs 2 • 10 −10 at 95% c.l. (6.4)(b) If the CSs are metastable, the explanation [18,25] of the recent NANOGrav15 data [16,17] on stochastic GWs is possible for M 0.9 YeV and κ 0.0003.(6.5)This is, because the obtained Gµ cs values, through Eq. (6.1), from the ranges above is confined in the range dictated by the interpretation of the recent data [18] 10 −8 Gµ cs 2.4 • 10 −4 for 8.2 √ r ms 7.5 at 95% c.l. (6.6)where the metastability factor r ms is the ratio of the monopole mass squared to µ cs .Since we do not investigate the monopole formation in our work, the last restriction does not impact on our parameters.Moreover, the GWs obtained from CSs have to be consistent with the upper bound on their abundance Ω GW h 2 originating from the advanced LVK third observing run [48] Ω GW h 2 6.96 • 10 −9 for f = 32 Hz, ( where f is the frequency of the observation.At last, although not relevant for our computation, let us mention for completeness that Ω GW h 2 should be smaller than the limit on dark radiation which is encoded in an upper limit on ∆N eff from BBN and CMB observations [61] Ω GW h 2 7 8 4 11

GWS FROM METASTABLE CSS WITH LOW REHEATING
We focus here on case (b) of Sec.6.1 and compute the spectrum of the GW produced by the CSs.The presence of the long-lasting matter domination obtained in our set up due to the z oscillations after the end of FHI -see Sec. 4 -has important ramifications in the shape of the spectrum of GW.This is, because a decaying-particle-dominated era can be approximated by a matter domination which leads to a spectral suppression at relatively large frequencies [46,47].
To verify this fact in the case of our model we compute the emitted GW background at frequency f following the standard formula of Ref. [47] Here ρ c is given in Eq. (3.7) and P k is the power spectrum of GWs emitted by the k th harmonic of a CS loop.Assuming cusps as the main source of GWs, P k is given by [47] P k ≃ Γ/ζ(4/3)k −4/3 with Γ = 50 and ζ(4/3) = 3.6.(6.10) In Eq. (6.9) Ω GW h 2 is expressed as a sum of the contributions from an infinite number of normal modes.On technical ground, however, we take a sum up to k max ∼ 10 5 to achieve a sufficiently accurate result.Following our cosmological scenario, the number of loops emitting GWs, observed at a given frequency f , can be found from  [17] and LVK [48] -and future -SKA [62], THEIA [63], µAres [64], LISA [65], Taiji [66], TianQin [67], BBO [68], DECIGO [69], ET [70] and CE [71] -experiments.
where we use as integration variable τ -cf.Ref. [47] -taking into account Eq. (5.10) and dz = −e −τ dτ.Also τ max , τ rh and τ eq are given in Sec. 5 and τ dc corresponds to the decay τ of CS network, which is estimated from τ dc = − ln (70/H 0 ) 1/2 (ΓΓ dc Gµ cs ) 1/4 + 1 with Γ dc = 4Gµ cs m 2 P e −πrms (6.12) the rate per unit length of the metastable CSs [31].For the loop number density per unit length n(ℓ, t) -with mass dimension 4 -we adopt the expressions [47] n r (ℓ, t) = 0.18 ) ) where the subscripts r and m are refereed to CSs produced and radiating in a radiation-and a matterdominated era respectively whereas rm means that the loops produced during the radiation domination, but radiating during the matter domination.Also ℓ Γ is the initial length of a loop which has length ℓ at a cosmic time t and is given by ℓ Γ = ℓ + ΓGµ cs t with Γ ≃ 50 and ℓ = 2ke τ /f.(6.14) Here Γ is related to the energy emission from CS [46] and introduces some uncertainty in the computation.The Hubble parameter during the standard cosmological evolution, H st , and during the z oscillation, H z , can be found from the formulas [54,55,57] where the various ρ i are given in Eq. (5.13).Lastly, the cosmic time as a function of τ is written as In particular, we have w = 1/3 or w = 0 for τ rh < τ ≤ τ eq or τ > τ eq respectively.Armed with the formulae above we compute Ω GW h 2 for the GWs produced from the CS formatted in our setting under the assumption that these are metastable.Employing the inputs of BPB in Table 2 -which yield Gµ cs = 6 • 10 −8 -we obtain the outputs displayed in Fig. 4. Namely, in Fig. 4-(a) we show Ω GW h 2 as a function of f for r 1/2 ms = 8 and µ = m/3 which yields T rh = 0.4 (dotted line) or µ = 3 m which results to T rh = 3.5 (solid line).On the other hand, for the GW spectra depicted in Fig. 4-(b) we employ µ = m resulting to T rh = 1.2 GeV and fix r 1/2 ms to 7.9 (dotted line) 8 (solid line) and 8.2 (dashed line) -see Eq. (6.6).In both panels of Fig. 4 we see that the derived GW spectra can explain NANOGrav15 data shown with gray almost vertical lines.We see though that, as r 1/2 ms increases, the increase of Ω GW h 2 becomes sharper and provide better fit to the observations.Also, in both panels the shape of GW signal suffers a diminishment above a turning frequency f rh ∼ 0.03 Hz which enables us to satisfy Eq. (6.7) more comfortably than in the case with high reheating -cf.Ref. [22,24].As T rh decreases, the reduction of Ω GW h 2 becomes more drastic in accordance with the findings of Ref. [46,47].The plots also show examples of sensitivities of possible future observatories [62][63][64][65][66][67][68][69][70][71] which can test the signals at various f values.Needless to say, the bound in Eq. (6.8) is not depicted since it is lies above the Ω GW h 2 values displayed in the plots.
It would be preferable to obtain larger M , and so Gµ cs , values (e.g., M ≃ 10 YeV yields Gµ cs ≃ 10 −6 ) such that the resulting Ω GW h 2 enters the dense part of the NANOGrav15 favorable region.This can be achieved -see e.g.Ref. [12,22] -if we use θ = π and a next-to-minimal version [2] for K I .However, the generation of a S from Z with θ = 0 in Eq. (4.5c) fixes the sign of the second term in Eq. (4.5b) and does not allow for the alternative arrangement mentioned above.On the other hand, next-to-next-to-minimal K I may accommodate [2,8,12] such an augmentation of the M value without disturbing the z stabilization during FHI -see Sec.4.1.In that case we expect that a S would not be constrained by the requirements of successful FHI and it could be selected so that T rh lies at the convenient level which allows for the comfortable evasion of the exclusion limits from the LVK [48] experiment.This setting opens up an interesting interconnection between the GW experiments and m.

PREDICTIONS FOR THE SUSY-MASS SCALE
The a S values displayed in Fig. 1 and Eq.(4.14) -which result from the observational constraints to FHI -give us the opportunity to gain information about the mass scale of SUSY particles through the determination of m = m 3/2 -see Eq. (3.9).This aim can be achieved by solving Eq. (4.5c) w.r.t.m as follows where we take into account Eq. (4.2) and the fact that z I /m P ∼ 10 −3 .As a consequence, the analytic result above approximates well the numerical one which is obtained by extracting consistently m as a function of κ and M determined by the conditions in Eqs.(4.8) and (4.9).Note that an iterative process has to be realized introducing a trial m value which allows us to use as input the form of V I in Eq. (4.4).The aforementioned smallness of z I causes a diminishment of m compared to a S rendering therefore m via Eqs.(3.9) and (3.8a) of the order of PeV scale.Indeed, for fixed ν, Eq. (7.1) yields m and then m 3/2 , m z and m θ can be easily estimated by Eqs.(3.8a) and (3.8c) whereas m and T rh by Eqs.(3.9) and (5.2).Their numerical values for the BPs A and B are presented in Table 2.The magnitude of the derived m values and the necessity for µ ∼ m, established in Sec. 5 hints towards the high-scale MSSM.
To highlight numerically our expectations above, we display the allowed (gray shaded) region in the κ − m plane fixing n s = 0.965 and varying ν and µ within their possible respective margins (0.75 − 1) and (1/3 − 3) m -see Fig. 5. Along the solid line we set ν = 7/8.From Eq. (7.1) we can convince ourselves that the lower boundary curve of the displayed region is obtained for ν ≃ 0.751, whereas the upper one corresponds to ν ≃ 0.99.Assuming also µ = m we can determine the slice of the area which can be excluded due to the BBN bound in Eq. (5.7).In all, we find that m turns out to be confined in the range 1.2 a S /TeV 100 and 0.34 m/PeV 13.6, ( whereas T max rh ≃ 71 GeV attained for µ = 3 m and ν ≃ 0.99.The derived allowed margin of m and the employed µ values render our proposal compatible with the mass of the Higgs boson discovered in LHC, if we adopt as a low energy effective theory the high-scale version of MSSM [50].Indeed, within high-scale SUSY, updated analysis requires [50] 3 TeV m 0.3 ZeV, (7.3) for degenerate sparticle spectrum, m/3 ≤ µ ≤ 3 m, 1 ≤ tan β ≤ 50 and varying the stop mixing.On the other hand, our setting does not fit well with natural [41] or split [50] SUSY which assume µ ≪ m.

CONCLUSIONS
We analyzed the implementation of FHI together with the SUSY breaking and the CS formation in a B-L extension of MSSM.We adopted the super-and Kähler potentials in Eqs.(2.1) and (2.3) applying an approximate R symmetry.The model offers the following interesting achievements: • Observationally acceptable FHI of hilltop-type adjusting the tadpole parameter a S and the B − L breaking scale M in the ranges of Eq. (4.14).
• A prediction for the SUSY-mass scale m which turns out to be of the order of PeV.
• Reheating generated due to the domination of the sgoldstino condensate after the end of FHI; since µ is of the order PeV, the resulting T rh is higher than its lower bound from BBN.
• B − L CSs which, if metastable, explain rather well the recent observations on GWs.A characteristic of the obtained spectra is their suppression at relatively large frequencies (f > 0.1 Hz) which is due to the long-lasting matter domination caused by the sgolstino oscillations.
• Generation of neutrino masses via the type I see-saw mechanism supported by W MD and K D in Eqs.(2.2e) and (2.6b).
A potential shortcoming of our proposal is that baryogenesis is made difficult due to the low reheat temperature.Non-thermal leptogenesis is not operative since the sgoldstino (with mass of order 1 PeV) is lighter than N c i which acquire mass of order 1 ZeV -cf.Ref. [10,15] -and so the direct decay of z into N c i is forbidden.However, there are extensions of MSSM [72] where the late decay of the sgoldstino may generate non-thermally the baryon asymmetry of the universe.Alternatively, this problem may be overcome applying improved attempts [73] based on the idea of cold electroweak baryogenesis [74].As regards the CDM, the candidacy of lightest neutralino has to be investigated thoroughly solving precisely the relevant Bolzmann equations, as in Ref. [55,75].In case that the abundance is low we can open up slightly the decay channel of z into G so that we obtain a controllable production of neutralinos.Another aspect is the fate of R axion which remains stable throughout our setting [76].
As regards naturalness, it is a puzzle why the SUSY should appear at such scale, which is higher than the electroweak scale making the fine-tuning severe, while it is much smaller than the fundamental energy scale such as the GUT or Planck scales.Our scenario provides a possible solution to this issue: this may be due to the inflationary selection.Namely, the apparent fine-tuning could be a result of combination of the FHI and a bias toward high-scale SUSY in the landscape.

APPENDIX B ABBREVIATIONS
For further clarity we list here the abbreviations used in this manuscript.

(
a) W I is the IS part of W which reads[1]

FIGURE 1 :
FIGURE 1: Allowed (shaded) regions as determined by Eqs.(4.8), (4.9), and (4.11) in the κ − a S (a) and M − a S (b) plane.The conventions adopted for the various lines are also shown.
) and Fig.1-(b).The boundaries of the allowed areas in Fig.1are determined by the dashed [dot-dashed] lines corresponding to the upper [lower] bound on n s in Eq. (4.11) -note that in Fig. 1-(a) we obtain a very narrow strip from the n s variation, We also display by solid lines the allowed contours for n s = 0.965.The maximal r's are encountered in the upper right end of the dashed linescorresponding to n s = 0.974.On the other hand, the maximal |α s |'s are achieved along the dot-dashed lines and the minimal value is α s = −3.2• 10 −4 .Summarizing our findings from Fig. 1 for central n s in Eq. (4.11) we end up with the following ranges: 0.07 M/YeV 2.6 and 0.1 a S /TeV 100.(4.14)

56 FIGURE 3 :
FIGURE 3: Evolution of log ρ i with i = z -dashed line for T rh = 1.5 GeV and dot-dashed line for T rh = 0.5 GeV -and i = r, m (solid lines) as a function of τ = − ln(1 + z).Shown are also τ max , τ rh , τ eq and τ dc for κ = 0.001 and a S = 25.3TeV.

TABLE 1 :
The representations under G B−L and the extra global charges of the superfields of our model.