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

Abstract

We consider F-term hybrid inflation (FHI) and SUSY breaking in the context of a $B-L$ extension of the MSSM that largely respects a global $U\left(1\right)$ R symmetry. The hidden sector Kaehler manifold enjoys an enhanced $SU(1,1)/U\left(1\right)$ symmetry, with the scalar curvature determined by the achievement of a SUSY-breaking de Sitter vacuum without undesirable tuning. FHI turns out to be consistent with the data, provided that the magnitude of the emergent soft tadpole term is confined to the range $(1.2$–$100)$ TeV, and it is accompanied by the production of $B-L$ cosmic strings. If these are metastable, they are consistent with the present observations from PTA experiments on the stochastic background of gravitational waves with dimensionless tension $G{\mu }_{\mathrm{cs}}\simeq (1-9.2)·{10}^{-8}$. The $\mu$ parameter of the 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.

1. Introduction

Supersymmetric (SUSY) hybrid inflation [1] based on F terms, henceforth called F-term hybrid inflation (FHI) for short [2,3], is undoubtedly a well-motivated inflationary model: for reviews, see Refs. [4,5]. The most notable reasons that support our statement above are the following:

• FHI is tied to a renormalizable superpotential uniquely determined by a gauge and global $U\left(1\right)$ R symmetries.
• FHI does not require fine-tuned parameters and trans-Planckian inflaton values.
• FHI can be reconciled with the Planck data [6]—fitted to the standard power-law cosmological model with Cold Dark Matter (CDM) and a cosmological constant ($\mathsf{\Lambda }$CDM)—if we properly take into account not only radiative corrections (RCs) but also corrections originating from supergravity (SUGRA) [7,8,9,10,11,12,13,14], as well as soft SUSY-breaking terms [15,16,17,18,19,20,21,22].
• 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 [23] 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. [24]), it has recently appeared as an interesting ingredient in cosmological model building. This is because the announced data from several pulsar timing array (PTA) experiments [25,26,27]—most notably the NANOGrav 15-year results (NANOGrav15) [28]—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 [29], 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 [30,31,32,33,34,35,36,37,38,39]. In particular, the observations can be interpreted if the CSs are meta- [40,41] or quasi-stable [42,43]. Both types of topologically unstable CSs arise from symmetry breaking $\mathbb{G}\to \mathbb{H}\times U\left(1\right)$, which produces monopoles. The subsequent $U\left(1\right)$ breaking yields CSs that connect monopoles with antimonopoles. Here, we 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\left(1\right)$ symmetry is embedded in a gauge group with a higher rank, such as Pati–Salam [35,36], flipped $SU\left(5\right)$ [34] or $SO\left(10\right)$ [31,32,33]. For this reason, in this work, we focus on FHI realized in a $B-L$ extension of the Minimal SUSY Standard Model (MSSM)—cf. Refs. [35,36,44]—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 monopole production or the metastability of CSs, as was the case, e.g., in Refs. [30,38,39].

We adopt the minimal possible framework [19,20] that supports observationally acceptable FHI. It employs minimal Kähler potential for the inflaton field, RCs and soft SUSY-breaking terms. Among the last ones, the tadpole term plays a crucial role in establishing compatibility with data. Its magnitude can be explained by intertwining [45] the inflationary sector (IS) with a hidden sector (HS) introduced in Refs. [46,47]. Contrary to earlier attempts [48,49,50,51,52,53,54], this HS respects a mildly violated R symmetry that 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/\nu$ of the goldstino superfield—which is related to the geometry of the HS [47]—is constrained to values with $0<\nu <1$.
• SUSY breaking is achieved not only in a Minkowski vacuum—as in the cases of Refs. [48,49,50,51,52,53,54]—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 the order ${10}^{-12}$.
• The sgoldstino is stabilized [49,50,51,52,53,54] to low values during FHI. This fact, together with the selection of a minimal Kähler potential for the inflaton, assists us in resolving the $\eta$-problem. Note that Kähler potentials inspired by string theory were mainly employed in earlier works [48,49,50,51,52,53,54].
• The solution to the $\mu$ problem [55] of the MSSM is achieved by suitably applying [46,47] the Giudice–Masiero mechanism [56,57]. Contrary to similar attempts [48,58,59,60], the $\mu$ 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 the sgoldstino condensate, which decays [61,62,63,64,65] before the onset of Big Bang Nucleosynthesis (BBN) at cosmic temperatures of $(2$–$4)\mathrm{MeV}$ [66] thanks to the aforementioned $\mu$ term. Therefore, our scenario naturally results in prolonged matter domination, which causes a reduction [67,68,69,70] in the spectra of GWs at high frequencies ($f>0.1\mathrm{Hz}$). This fact is welcome since it assists us in avoiding any conflicts with the third run of advanced LIGO-VIRGO-KAGRA (LVK) data [71].
• The SUSY mass scale $\tilde{m}$ is predicted to be close to the $\mathrm{PeV}$ scale [72]. It fits well with the Higgs boson mass, discovered at the LHC, as it is estimated [73] within high-scale SUSY if we assume a relatively low $tan\beta$ and stop mixing. Note that the connection of inflation with SUSY breaking has been extensively discussed in the literature [74,75,76,77,78,79,80,81] over the last several years.

In this feature paper, we further review the model introduced in Ref. [45], focusing exclusively on its implementation within a version of the MSSM endowed with a $B-L$ Higgs sector. We present the complete particle content of the model, paying special attention to the computation of the GWs emerging from CS decay, under the assumption that those are metastable. We also explain the generation of neutrino masses, taking into account SUGRA contributions from Ref. [82].

The remainder of the manuscript is organized as follows: In Section 2, we introduce our framework. Then, we revise the salient features of our model as regards its vacuum in Section 3 and the inflationary era in Section 4. Then, we study the reheating process in Section 5 and the production of GWs from CSs in Section 6. Our predictions for the SUSY mass scale are revealed in Section 7. Our conclusions are summarized in Section 8. Details on the derivation of neutrino masses are given in Appendix A.

2. Model Set-Up

We focus on an extension of the MSSM invariant under the gauge group ${\mathbb{G}}_{B-L}={\mathbb{G}}_{\mathrm{SM}}\times U{\left(1\right)}_{B-L}$, where ${\mathbb{G}}_{\mathrm{SM}}$ is the Standard Model gauge group. The charge assignments under these symmetries of various matter and Higgs superfields are listed in

Table 1. The representations under ${\mathbb{G}}_{B-L}$ and the extra global charges of the superfields of our model.

Superfields	Representations	Global Symmetries
	under ${\mathbb{G}}_{\mathit{B}-\mathit{L}}$	$\mathit{R}$	$\mathit{B}$	$\mathit{L}$
Matter Superfields
$e_i^c$	$(\mathbf{1},\mathbf{1},1,1)$	0	0	$-1$
$N_i^c$	$(\mathbf{1},\mathbf{1},0,1)$	0	0	$-1$
$L_i$	$(\mathbf{1},\mathbf{2},-1/2,-1)$	0	0	1
$u_i^c$	$(\mathbf{3},\mathbf{1},-2/3,-1/3)$	0	$-1/3$	0
$d_i^c$	$(\mathbf{3},\mathbf{1},1/3,-1/3)$	0	$-1/3$	0
$Q_i$	$(\overline{\mathbf{3}},\mathbf{2},1/6,1/3)$	0	$1/3$	0
Higgs Superfields
$H_d$	$(\mathbf{1},\mathbf{2},-1/2,0)$	2	0	0
$H_u$	$(\mathbf{1},\mathbf{2},1/2,0)$	2	0	0
ine S	$(\mathbf{1},\mathbf{1},0,0)$	2	0	0
$\mathsf{\Phi }$	$(\mathbf{1},\mathbf{1},0,2)$	0	0	$-2$
$\overline{\mathsf{\Phi }}$	$(\mathbf{1},\mathbf{1},0,-2)$	0	0	2
Goldstino Superfield
Z	$(\mathbf{1},\mathbf{1},0,0)$	$2/\nu$	0	0

. In particular, the ith generation $SU{\left(2\right)}_{\mathrm{L}}$ doublet left-handed quark and lepton superfields are denoted by $Q_i$ and $L_i$, respectively, whereas the $SU{\left(2\right)}_{\mathrm{L}}$ singlet antiquark [antilepton] superfields are denoted by $u_i^c$ and ${d_i}^c$ [$e_i^c$ and $N_i^c$], respectively. The electroweak Higgs superfields that 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_i^c$s, a pair of Higgs superfields, $\mathsf{\Phi }$ and $\overline{\mathsf{\Phi }}$, which break $U{\left(1\right)}_{B-L}$, and the goldstino superfield Z. In addition to local symmetry, the model also possesses baryon and lepton number symmetries and an R symmetry $U{\left(1\right)}_R$. The latter plays a crucial role in the construction of the superpotential (see Section 2.1) and the Kähler potential (see Section 2.2).

2.1. Superpotential

The superpotential of our model fully respects the symmetries in Table 1. Most notably, it carries R charge 2 and is linear with respect to S and $Z^{\nu }$. It naturally splits into five parts:

$$W=W_{\mathrm{I}}+W_{\mathrm{H}}+W_{\mathrm{GH}}+W_{\mathrm{MSSM}}+W_{\mathrm{MD}},$$

(1)

where the subscripts “I” and “H” stand for the IS and HS, respectively, and the content of each term is specified as follows:

(a) $W_{\mathrm{I}}$ is the IS part of W and reads [1]

$$W_{\mathrm{I}}=\kappa S\left(\overline{\mathsf{\Phi }}\mathsf{\Phi }-M^2\right),$$

(2a)

where $\kappa$ and M are free parameters that can be made positive by field redefinitions.

(b) $W_{\mathrm{H}}$ is the HS part of W and is written as [47]

$$W_{\mathrm{H}}=m{m}_{\mathrm{P}}^2{(Z/m_{\mathrm{P}})}^{\nu }.$$

(2b)

Here, $m_{\mathrm{P}}=2.4\mathrm{ReV}$ is the reduced Planck mass—with $\mathrm{ReV}={10}^{18}\mathrm{GeV}$—m is a positive free parameter with mass dimensions, and $\nu$ is an exponent that can, in principle, acquire any real value if $W_{\mathrm{H}}$ is considered an effective superpotential that is valid close to the non-zero vacuum expectation value (v.e.v) of Z, $\langle Z\rangle$. 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^{{\nu }_s}{Z}^{{\nu }_z}$ can be forbidden in W since the exponent of such a term has to obey the relation

$$\begin{array}{c}\hfill {\nu }_s+{\nu }_z/\nu =1\Rightarrow {\nu }_z=(1-{\nu }_s)\nu ,\end{array}$$

leading to negative values of ${\nu }_z$. This conclusion contradicts our assumptions above.

(c) $W_{\mathrm{GH}}$ is a term that mixes the HS and the $B-L$ gauge fields of the IS. It has the form

$$W_{\mathrm{GH}}=-\lambda m_{\mathrm{P}}{(Z/m_{\mathrm{P}})}^{\nu }\overline{\mathsf{\Phi }}\mathsf{\Phi }$$

(2c)

with $\lambda$ being a real coupling constant. The magnitude of $\lambda$ can be restricted by the DE requirement, as we see below.

(d) $W_{\mathrm{MSSM}}$ is the part of W that contains the usual trilinear terms of the MSSM, i.e.,

$$W_{\mathrm{MSSM}}=h_{ijD}{d}_i^c{Q}_j{H}_d+h_{ijU}{u}_i^c{Q}_j{H}_u+h_{ijE}{e}_i^c{L}_j{H}_d.$$

(2d)

The selected R assignments in Table 1 prohibit the bilinear $\mu H_u{H}_d$ term of the MSSM and other mixing terms—e.g., ${\lambda }_{\mu }S{H}_u{H}_d$ [60], which is frequently employed to generate this $\mu$ term—from being present in $W_{\mathrm{MSSM}}$. This term is generated here via $K_{\mu }$; see Section 2.2 below.

(e) $W_{\mathrm{MD}}$ is the part of W that provides masses to neutrinos:

$$W_{\mathrm{MD}}=h_{ijN}{N}_i^c{L}_j{H}_u+{\lambda }_{{iN}^c}\left(S+{(Z/m_{\mathrm{P}})}^{\nu }\right)\overline{\mathsf{\Phi }}{N}_i^{c2}.$$

(2e)

The first term on the right-hand side of Equation (2e) is responsible for Dirac neutrino masses, whereas for Majorana masses, cf. Refs. [19,24]. The scale of the latter masses is intermediate since $\langle \overline{\mathsf{\Phi }}\rangle \sim 1\mathrm{YeV}$ and $\langle Z\rangle \sim m_{\mathrm{P}}$. The cooperation of both terms leads to the light neutrino masses via the well-known (type I) seesaw mechanism; see also Appendix A.

2.2. Kähler Potential

The Kähler potential respects the ${\mathbb{G}}_{B-L}$, B and L symmetries in Table 1. It has the following contributions:

$$K=K_{\mathrm{I}}+K_{\mathrm{H}}+K_{\mu }+K_{\mathrm{D}}+{\left|Y_{\alpha }\right|}^2,$$

(3)

in which the left-handed chiral superfields of the MSSM are denoted by $Y_{\alpha }$, with $\alpha =1,\dots ,7$, i.e.,

$$\begin{array}{c}\hfill Y_{\alpha }=Q,L,d^c,u^c,e^c,N^c,H_d,\mathrm{and}{H}_u,\end{array}$$

where the generation indices are suppressed.

(a) $K_{\mathrm{I}}$ is the part of K that depends on the fields involved in FHI—cf. Equation (2a). We adopt the simplest possible choice—cf. Refs. [2,19]—which has the form

$$K_{\mathrm{I}}={\left|S\right|}^2+{\left|\mathsf{\Phi }\right|}^2+{\left|\overline{\mathsf{\Phi }}\right|}^2.$$

(4a)

Higher-order terms of the form ${\left|S\right|}^{2{\nu }_S}/m_{\mathrm{P}}^{2{\nu }_S-2}$ with ${\nu }_S>1$ cannot be excluded by the imposed symmetries but may become harmless if $S\ll m_{\mathrm{P}}$ and assuming low enough coefficients.

(b) $K_{\mathrm{H}}$ is the part of K devoted to the HS. We adopt the form introduced in Ref. [47], where

$$K_{\mathrm{H}}=N{m}_{\mathrm{P}}^{2}ln\left(1+\frac{{\left|Z\right|}^2-k^2{Z}_{-}^4/m_{\mathrm{P}}^2}{N{m}_{\mathrm{P}}^2}\right)\mathrm{with}{Z}_{\pm }=Z\pm Z^{*}.$$

(4b)

Here, $k>0$ mildly violates R symmetry, endowing the R axion with phenomenologically acceptable mass. The selected $K_{\mathrm{H}}$ is not motivated by string theory, but it can be considered 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 Section 3—a dS vacuum of the whole-field system with a tunable cosmological constant for

$$N=\frac{4{\nu }^2}{3-4\nu }\mathrm{with}\frac{3}{4}<\nu <\frac{3}{2}\mathrm{for}N<0.$$

(5)

Our favored $\nu$ range, finally, is $3/4<\nu <1$. Since $N<0$, $K_{\mathrm{H}}$ parameterizes the $SU(1,1)/U\left(1\right)$ hyperbolic Kähler manifold for $k\sim 0$.

(c) $K_{\mu }$ includes higher-order terms that generate the needed mixing term between $H_u$ and $H_d$ in the Lagrangian of the MSSM—cf. Refs. [46,47,56]—and has the form

$$K_{\mu }={\lambda }_{\mu }\left({Z^{*}}^{2\nu }/m_{\mathrm{P}}^{2\nu }\right)H_u{H}_d+\mathrm{h}.\mathrm{c}.,$$

(6a)

where the dimensionless constant ${\lambda }_{\mu }$ is taken as real for simplicity.

(d) $K_{\mathrm{D}}$ is an unavoidable term that mixes the observable sector with the HS as follows:

$$K_{\mathrm{D}}={\lambda }_{ij\mathrm{D}}\left(Z^{*\nu }/m_{\mathrm{P}}^{\nu +1}\right)N_i^c{L}_j{H}_u+\mathrm{h}.\mathrm{c}.$$

(6b)

It provides (subdominant) Dirac masses for ${\nu }_i$ [82], as shown in Appendix A.

The total K in Equation (3) enjoys an enhanced symmetry for the $Y^{\alpha },S$ and Z fields, namely,

$$\prod _{\alpha }U{\left(1\right)}_{Y^{\alpha }}\times U{\left(1\right)}_S\times {\left(SU(1,1)/U\left(1\right)\right)}_Z,$$

(7)

where the indices indicate the moduli that parameterize the corresponding manifolds. Thanks to this symmetry, mixing terms of the form $S^{{\tilde{\nu }}_s}{Z}^{*{\tilde{\nu }}_z}$ can be ignored, although they may be allowed by the R symmetry for ${\tilde{\nu }}_z=\nu {\tilde{\nu }}_s$. Most notably, $U{\left(1\right)}_S$ protects $K_{\mathrm{I}}$ from S-dependent terms, which violates the R symmetry, thereby spoiling the inflationary set-up.

3. SUSY and ${\mathbb{G}}_{\mathbf{B}-\mathbf{L}}$ Breaking—Dark Energy

The vacuum of our model is determined by minimizing the F-term (tree-level) SUGRA scalar potential $V_{\mathrm{F}}$ derived [45] from W in Equation (1) and K in Equation (3). Note that D-term contributions to the total SUGRA scalar potential vanish if we confine ourselves to the period during FHI and to the vacuum along the D-flat direction

$$|\overline{\mathsf{\Phi }}|=|\mathsf{\Phi }|\mathrm{which} \mathrm{ensures}{V}_{\mathrm{D}}=\frac{g^2}{2}{\left({\left|\mathsf{\Phi }\right|}^2-{\left|\overline{\mathsf{\Phi }}\right|}^2\right)}^2=0.$$

(8)

Here, g is the unique (considered unified) gauge coupling constant of ${\mathbb{G}}_{B-L}$, and we employ the same symbol for the various superfields $X^{\alpha }=S,Z,\mathsf{\Phi },\overline{\mathsf{\Phi }}$ and their complex scalar components.

As we can verify numerically, $V_{\mathrm{F}}$ is minimized at the ${\mathbb{G}}_{B-L}$-breaking vacuum

$$\left|\langle \mathsf{\Phi }\rangle \right|=\left|\langle \overline{\mathsf{\Phi }}\rangle \right|=M.$$

(9)

It also has a stable valley along $\langle \theta \rangle =0$ and $\langle {\theta }_S/m_{\mathrm{P}}\rangle =\pi$, where these fields are defined by

$$Z=(z+i\theta )/\sqrt{2}\mathrm{and}S=\sigma e^{i{\theta }_S/m_{\mathrm{P}}}/\sqrt{2}.$$

(10)

Substituting Equation (10) into $V_{\mathrm{F}}$ and minimizing it with respect to the various directions, we arrive at the following results:

$$\sigma =-2^{(1-\nu )/2}\left(\lambda (M^2+m_{\mathrm{P}}^2)-m{m}_{\mathrm{P}}\right)z^{\nu }/m_{\mathrm{P}}^{(\nu +1)}\mathrm{and}\langle z\rangle =2\sqrt{2/3}\left|\nu \right|m_{\mathrm{P}},$$

(11)

which yield the constant potential energy density,

$$\langle V_{\mathrm{F}}\rangle ={\left(\frac{16{\nu }^4}{9}\right)}^{\nu }{\left(\frac{\lambda M^2-m{m}_{\mathrm{P}}}{\kappa m_{\mathrm{P}}^2}\right)}^2{\omega }^N{\left(\lambda (M^2+m_{\mathrm{P}}^2)-m{m}_{\mathrm{P}}\right)}^2\mathrm{with}\omega =e^{\frac{\langle K_{\mathrm{H}}\rangle }{N{m}_{\mathrm{P}}^2}}\simeq \frac{2}{3}(3-2\nu ).$$

(12)

Tuning $\lambda$ to a value $\lambda \sim m/m_{\mathrm{P}}\simeq {10}^{-12}$, we may wish to identify $\langle V_{\mathrm{F}}\rangle$ with the DE energy density, i.e.,

$$\langle V_{\mathrm{F}}\rangle ={\mathsf{\Omega }}_{\mathsf{\Lambda }}{\rho }_{\mathrm{c}}=7.3·{10}^{-121}{m}_{\mathrm{P}}^4,$$

(13)

where the density parameter of DE and the current critical energy density of the universe are, respectively, given by [83]

$${\mathsf{\Omega }}_{\mathsf{\Lambda }}=0.6889\mathrm{and}{\rho }_{\mathrm{c}}=2.31·{10}^{-120}{h}^2{m}_{\mathrm{P}}^4\mathrm{with}h=0.6766.$$

(14)

Therefore, we obtain a post-inflationary dS vacuum, which explains the notorious DE problem. Moreover, Equation (11) yields $\langle \sigma \rangle \simeq 0$.

The particle spectrum of the theory at the vacuum in Equations (9) and (11) includes the gravitino ($\tilde{G}$), which acquires mass as follows [47]:

$$m_{3/2}=\langle e^{K_{\mathrm{H}}/2m_{\mathrm{P}}^2}{W}_{\mathrm{H}}/m_{\mathrm{P}}^2\rangle \simeq 2^{\nu }{3}^{-\nu /2}{\left|\nu \right|}^{\nu }m{\omega }^{N/2}.$$

(15a)

Diagonalizing the mass-squared matrix of the field system $S-\mathsf{\Phi }-\overline{\mathsf{\Phi }}-Z$, we also find that the IS acquires a common mass:

$$m_{\mathrm{I}}=e^{K_{\mathrm{H}}/2m_{\mathrm{P}}^2}\sqrt{2}{\left({\kappa }^2{M}^2+{(4{\nu }^2/3)}^{\nu }(1+4M^2/m_{\mathrm{P}}^2)m^2\right)}^{1/2},$$

(15b)

where the second term arises due to the coexistence of the IS with the HS—cf. Ref. [19]. Similar mixing does not appear in the mass spectrum of the HS, which contains the (canonically normalized) sgoldstino (or R saxion) and the pseudo-sgoldstino (or R axion) with the following respective masses:

$$m_z\simeq \frac{3\omega }{2\nu }{m}_{3/2}\mathrm{and}{m}_{\theta }\simeq 12k{\omega }^{3/2}{m}_{3/2}.$$

(15c)

Finally, applying the relevant formulas of Refs. [47,57], we find that $K_{\mu }$ induces a non-vanishing $\mu$ term in the superpotential of the MSSM, whereas $W_{\mathrm{MSSM}}$ and $K_{\mu }+{\left|Y^{\alpha }\right|}^2$ in Equations (2d) and (3) lead to a common soft SUSY-breaking mass parameter $\tilde{m}$ at the vacuum in Equation (11), which indicatively represents the mass level of the SUSY partners. Specifically, we obtain

$$W\ni \mu H_u{H}_d\mathrm{with}\left|\mu \right|={\lambda }_{\mu }{\left(\frac{4{\nu }^2}{3}\right)}^{\nu }(5-4\nu )m_{3/2}\mathrm{and}\tilde{m}=m_{3/2}.$$

(16)

The variant forms of the terms $|Y^{\alpha }|$ in K (see Equation (3)) do not significantly alter our results [45,47].

4. Inflation Analysis

It is well known [1,4] that, in global SUSY, FHI takes place for sufficiently large $\left|S\right|$ values along the F- and D-flat directions of the SUSY potential:

$$\overline{\mathsf{\Phi }}=\mathsf{\Phi }=0,\mathrm{where}{V}_{\mathrm{SUSY}}\left(\mathsf{\Phi }=0\right)\equiv V_{\mathrm{I}0}={\kappa }^2{M}^4\mathrm{are}{H}_{\mathrm{I}}=\sqrt{V_{\mathrm{I}0}/3m_{\mathrm{P}}^2}$$

(17)

which are the constant potential energy density and corresponding Hubble parameter that drive FHI; the subscript 0 means that this is the tree-level value. In a SUGRA context, though, we first check—in Section 4.1—the conditions under which such a scheme can be achieved, and then, in Section 4.2, we give the final form of the inflationary potential. Lastly, we present our results in Section 4.4, imposing a number of constraints listed in Section 4.3.

4.1. Hidden Sector’s Stabilization

The attainment of FHI is possible if Z is well stabilized during it. The relative mechanism is pretty well known [52,53,54]. Due to $V_{\mathrm{I}0}$, Z is transported from its value in Equation (9) to a value well below $m_{\mathrm{P}}$. To determine this, we construct the complete expression for $V_{\mathrm{F}}$ along the inflationary trajectory in Equation (17) and then expand the resulting expression for low $S/m_{\mathrm{P}}$ values, assuming that the $\theta =0$ direction is stable, as in Equation (9). Under these conditions, $V_{\mathrm{F}}$ is minimized for the value

$${\langle z\rangle }_{\mathrm{I}}\simeq {\left(\sqrt{3}·2^{\nu /2-1}{H}_{\mathrm{I}}/m\nu \sqrt{1-\nu }\right)}^{1/(\nu -2)}{m}_{\mathrm{P}}.$$

(18)

This result is in good agreement with its precise value derived numerically. Note that $\nu <1$ ensures a real value of ${\langle z\rangle }_{\mathrm{I}}$ with ${\langle z\rangle }_{\mathrm{I}}\ll m_{\mathrm{P}}$ since $H_{\mathrm{I}}/m\ll 1$.

The (canonically normalized) components of the sgoldstino acquire masses squared, respectively,

$$m_{\mathrm{I}z}^2\simeq 6(2-\nu )H_{\mathrm{I}}^2\mathrm{and}{m}_{\mathrm{I}\theta }^2\simeq 3H_{\mathrm{I}}^2-m^2\left(8{\nu }^2{m}_{\mathrm{P}}^2-3{\langle z\rangle }_{\mathrm{I}}^2\right)\frac{4\nu (1-\nu )m_{\mathrm{P}}^2+(1-96k^2\nu ){\langle z\rangle }_{\mathrm{I}}^2}{2^{3+\nu }\nu m_{\mathrm{P}}^{2\nu }{\langle z\rangle }_{\mathrm{I}}^{2(2-\nu )}},$$

(19a)

whereas the mass of $\tilde{G}$ turns out to be

$$m_{\mathrm{I}3/2}\simeq {\left(\nu {(1-\nu )}^{1/2}{m}^{2/\nu }/\sqrt{3}{H}_{\mathrm{I}}\right)}^{\nu /(2-\nu )}.$$

(19b)

It is clear from the results above that $m_{\mathrm{I}z}\gg H_{\mathrm{I}}$, and therefore, it is well stabilized during FHI, whereas $m_{\mathrm{I}\theta }\simeq H_{\mathrm{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_{\mathrm{I}\theta }\ll H_{\mathrm{I}}$.

4.2. Inflationary Potential

Expanding $V_{\mathrm{F}}$ for low S values, introducing the canonically normalized inflaton $\sigma =\sqrt{2}\left|S\right|$ and taking into account the RCs [1,4], we derive [45] the inflationary potential $V_{\mathrm{I}}$, which can be cast in the form

$$V_{\mathrm{I}}\simeq V_{\mathrm{I}0}\left(1+C_{\mathrm{RC}}+C_{\mathrm{SSB}}+C_{\mathrm{SUGRA}}\right).$$

(20)

The individual contributions are specified as follows:

(a) $C_{\mathrm{RC}}$ represents the RCs to $V_{\mathrm{I}}/V_{\mathrm{I}0}$, which can be written as [2,4]

$$C_{\mathrm{RC}}=\frac{{\kappa }^2}{128{\pi }^2}\left(8ln\frac{{\kappa }^2{M}^2}{Q^2}+8x^2{tanh}^{-1}\left(\frac{2}{x^2}\right)-4(ln4-x^{4}lnx)+(4+x^4)ln(x^4-4)\right),$$

(21a)

with $x=(\sigma -\sqrt{2}\lambda {\langle Z\rangle }_{\mathrm{I}}^{\nu }{m}_{\mathrm{P}}^{1-\nu }/\kappa )/M>\sqrt{2}$, beyond which the expression above ceases to be valid. Here, we take into account that the multiplicity of the waterfall fields is ${\mathsf{N}}_{\mathrm{G}}=1$ since these are $U{\left(1\right)}_{B-L}$ non-singlets.

(b) $C_{\mathrm{SSB}}$ is the contribution to $V_{\mathrm{I}}/V_{\mathrm{I}0}$ from the soft SUSY-breaking effects [15,16,17,18,19] and is parameterized as follows:

$$C_{\mathrm{SSB}}=m_{\mathrm{I}3/2}^2{\sigma }^2/2V_{\mathrm{I}0}-{\mathrm{a}}_S\sigma /\sqrt{2V_{\mathrm{I}0}},$$

(21b)

where the tadpole parameter reads

$${\mathrm{a}}_S=2^{1-\nu /2}m\frac{{\langle z\rangle }_{\mathrm{I}}^{\nu }}{m_{\mathrm{P}}^{\nu }}\left(1+\frac{{\langle z\rangle }_{\mathrm{I}}^2}{2N{m}_{\mathrm{P}}^2}\right)\left(2-\nu -\frac{3{\langle z\rangle }_{\mathrm{I}}^2}{8\nu m_{\mathrm{P}}^2}\right).$$

(21c)

The minus sign results from the stabilization of $\theta$ (see Equation (10)) at zero and the minimization of the factor $(S+S^{*})=\sqrt{2}\sigma cos({\theta }_S/m_{\mathrm{P}})$, which occurs for ${\theta }_S/m_{\mathrm{P}}=\pi \left(\mathsf{mod}2\pi \right)$; the decomposition of S is shown in Equation (10). We further assume that ${\theta }_S$ remains constant during FHI so that the simple one-field slow-roll approximation is valid. Possible variation in ${\theta }_S$ was investigated in Ref. [20], where they found that acceptable solutions with ${\theta }_S/m_{\mathrm{P}}\ne \pi$ require a significant amount of tuning. The first term in Equation (21b) does not play an essential role in our set-up due to low enough $m_{3/2}$s— cf. Ref. [19].

(c) $C_{\mathrm{SUGRA}}$ is the SUGRA correction to $V_{\mathrm{I}}/V_{\mathrm{I}0}$ after subtracting the one in $C_{\mathrm{SSB}}$. It reads [15,16,17,18,19]

$$C_{\mathrm{SUGRA}}=c_{2\nu }\frac{{\sigma }^2}{2m_{\mathrm{P}}^2}+c_{4\nu }\frac{{\sigma }^4}{4m_{\mathrm{P}}^4}\mathrm{with}{c}_{2\nu }=\frac{{\langle z\rangle }_{\mathrm{I}}^2}{2m_{\mathrm{P}}^2}\mathrm{and}{c}_{4\nu }=\frac{1}{2}\left(1+{\langle z\rangle }_{\mathrm{I}}^2{m}_{\mathrm{P}}^2\right).$$

(21d)

Thanks to the minimality of $K_{\mathrm{I}}$ in Equation (4a) and the smallness of ${\langle z\rangle }_{\mathrm{I}}$, the coefficients above are low enough and allow for FHI to be established.

4.3. Observational Requirements

The analysis of FHI can be performed in the slow-roll approximation if we calculate the slow-roll parameters [6]

$$ϵ=m_{\mathrm{P}}^2{\left(\frac{V_{\mathrm{I}}^{\prime }}{\sqrt{2}{V}_{\mathrm{I}}}\right)}^2\simeq \frac{m_{\mathrm{P}}^2}{2}{\left(C_{\mathrm{RC}}^{\prime }+C_{\mathrm{SSB}}^{\prime }\right)}^2\mathrm{and}\eta =m_{\mathrm{P}}^2\frac{V_{\mathrm{I}}^{″}}{V_{\mathrm{I}}}\simeq m_{\mathrm{P}}^2{C}_{\mathrm{RC}}^{″},$$

(22)

where the derivatives of the various contributions read

$$\begin{array}{cc}\hfill C_{\mathrm{SSB}}^{\prime }& \simeq -{\mathrm{a}}_S/\sqrt{2V_{\mathrm{I}0}},\hfill \end{array}$$

(23a)

$$\begin{array}{cc}\hfill C_{\mathrm{RC}}^{\prime }& \simeq \frac{k{p}^{2}x}{32M{\pi }^2}\left(4{tanh}^{-1}\left(\frac{2}{x^2}\right)+x^{2}ln\left(1-\frac{4}{x^4}\right)\right),\hfill \end{array}$$

(23b)

$$\begin{array}{cc}\hfill C_{\mathrm{RC}}^{″}& \simeq \frac{{\kappa }^2}{32M^2{\pi }^2}\left(4{tanh}^{-1}\left(\frac{2}{x^2}\right)+3x^{2}ln\left(1-\frac{4}{x^4}\right)\right).\hfill \end{array}$$

(23c)

The required behavior of $V_{\mathrm{I}}$ in Equation (20) can be obtained thanks to the relation $C_{\mathrm{RC}}^{\prime }\simeq -C_{\mathrm{SSB}}^{\prime }$, which is established for carefully selecting $\kappa$ (or M) and ${\mathrm{a}}_S$. Apparently, we have $C_{\mathrm{SSB}}^{\prime }<0$ and $C_{\mathrm{RC}}^{\prime }>0$ for ${\sigma }_{★}<{\sigma }_{\max }$ since $|4{tanh}^{-1}\left(2/x^2\right)|>|x^{2}ln(1-4/x^4)|$. On the contrary, $C_{\mathrm{RC}}^{″}<0$, since the negative contribution $3x^{2}ln(1-4/x^4)$ dominates over the first positive one, and so we obtain $\eta <0$, giving rise to acceptably low $n_{\mathrm{s}}$ values.

Our model of FHI can be qualified if we test it against a number of observational requirements:

(a) Thenumber of e-foldings elapsed between the horizon crossing of the pivot scale $k_{★}=0.05/\mathrm{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 Section 5) that FHI is followed, in turn, by matter- and radiation-dominated eras, the relevant condition takes the following form [2,6]:

$$N_{\mathrm{I}★}={\int }_{{\sigma }_{\mathrm{f}}}^{{\sigma }_{★}}\frac{d\sigma }{m_{\mathrm{P}}^2}\frac{V_{\mathrm{I}}}{V_{\mathrm{I}}^{\prime }}\simeq 19.4+\frac{2}{3}ln\frac{V_{\mathrm{I}0}^{1/4}}{1\mathrm{GeV}}+\frac{1}{3}ln\frac{T_{\mathrm{rh}}}{1\mathrm{GeV}},$$

(24)

where the prime denotes derivation with respect to $\sigma$, ${\sigma }_{★}$ is the value of $\sigma$ when $k_{★}$ crosses outside the horizon of FHI, and ${\sigma }_{\mathrm{f}}$ is the value of $\sigma$ at the end of FHI. This is normally obtained by the critical point ${\sigma }_{\mathrm{c}}=\sqrt{2}\left|S_{\mathrm{c}}\right|$; i.e., the end of inflation coincides with the onset of the $B-L$ phase transition. Note that $\langle \sigma \rangle \simeq 0$, as mentioned below in Equation (13), and so it does not disturb the inflationary dynamics that govern the $\sigma$ evolution for $\sigma \ge {\sigma }_{\mathrm{c}}$.

(b) The amplitude $A_{\mathrm{s}}$ of the power spectrum of the curvature perturbation generated by $\sigma$ during FHI must be consistent with the data [83] on the cosmic microwave background (CMB), i.e.,

$$A_{\mathrm{s}}=\frac{1}{12{\pi }^2{m}_{\mathrm{P}}^6}\frac{V_{\mathrm{I}}^3\left({\sigma }_{★}\right)}{|V_{\mathrm{I}}^{\prime }\left({\sigma }_{★}\right){|}^2}\simeq 2.105·{10}^{-9}.$$

(25)

The observed curvature perturbation is generated wholly by $\sigma$ since the other scalars are massive enough during FHI; see Section 4.1.

(c) The remaining observables—the scalar spectral index $n_{\mathrm{s}}$, its running ${\alpha }_{\mathrm{s}}$ and the scalar-to-tensor ratio r—which are calculated by the standard formulas

$$n_{\mathrm{s}}=1-6{ϵ}_{★}+2{\eta }_{★},{\alpha }_{\mathrm{s}}=2\left(4{\eta }_{★}^2-{(n_{\mathrm{s}}-1)}^2\right)/3-2{\xi }_{★}\mathrm{and}r=16{ϵ}_{★},$$

(26)

(where $\xi \simeq m_{\mathrm{P}}^4{V}_{\mathrm{I}}^{\prime }{V}_{\mathrm{I}}^{‴}/V_{\mathrm{I}}^2$ and all the variables with the subscript $★$ are evaluated at $\sigma ={\sigma }_{★}$) must be in agreement with the data. We take into account the latest Planck release 4 (PR4) (including TT, TE, EE + lowE power spectra [83]), Baryon Acoustic Oscillations (BAO), CMB-lensing and BICEP/Keck (BK18) data [84]. Fitting them with $\mathsf{\Lambda }$CDM$+r$, we obtain approximately

$$n_{\mathrm{s}}=0.965\pm 0.009\mathrm{and}r\lesssim 0.032,$$

(27)

at a 95% confidence level (c.l.) with negligible $|{\alpha }_{\mathrm{s}}| \ll 0.01$.

4.4. Results

As deduced from Section 4.1, Section 4.2 and Section 4.3, the inflationary part of our model depends on the parameters

$$\begin{array}{c}\hfill \kappa ,M,m,\lambda ,\nu ,k\mathrm{and}{\lambda }_{\mu }.\end{array}$$

Recall that N is related to $\nu$ via Equation (5). Enforcing Equation (13) fixes $\lambda$ at a rather low value that does not influence our remaining results. Moreover, k exclusively affects $m_{\theta }$ and $m_{\mathrm{I}\theta }$ via Equations (15c) and (19a). Throughout, we select the value $k=0.1$, which ensures the avoidance of massless modes. Here, we present the ${\mathrm{a}}_S$ values as a function of $\kappa$ or M, which assist the comparison of FHI with data—cf. Ref. [19]—and postpone their derivation from m and $\nu$ via Equation (21c) in Section 7. As regards $T_{\mathrm{rh}}$, which influences Equation (24), we adopt a value close to those observed in our set-up, $T_{\mathrm{rh}}\simeq 1\mathrm{GeV}$.

Enforcing Equations (24) and (25), we can restrict M and ${\sigma }_{★}$ as functions of our free parameters $\kappa$ and ${\mathrm{a}}_S$. It is apparent that the allowed ranges of parameters are similar to those explored in Ref. [19], where the HS is not specified. Some deviations are only due to improvements in the determination of the $n_{\mathrm{s}}$ values in Equation (27). The correct values of that quantity are attained if FHI becomes of the hilltop type [19,20], i.e., if $V_{\mathrm{I}}$ is non-monotonic and develops a maximum at ${\sigma }_{\max }$ and a minimum at ${\sigma }_{\min }\gg {\sigma }_{\max }$. For $\sigma >{\sigma }_{\max }$, $V_{\mathrm{I}}$ becomes a monotonically increasing function of $\sigma$, and so its boundedness is ensured. FHI takes place for $\sigma <{\sigma }_{\max }$. The position of ${\sigma }_{\max }$ is predominantly dominated by $C_{\mathrm{RC}}^{\prime }$ and $C_{\mathrm{SSB}}^{\prime }$. On the other hand, ${\sigma }_{\min }$ can be roughly found by the interplay of $C_{\mathrm{SUGRA}}^{\prime }$ and $C_{\mathrm{SSB}}^{\prime }$. Specifically, we obtain [45]

$${\sigma }_{\max }\simeq \frac{{\kappa }^3{M}^2}{4\sqrt{2}{\pi }^2{\mathrm{a}}_S}\mathrm{and}{\sigma }_{\min }\simeq {\left(\frac{{\mathrm{a}}_S{m}_{\mathrm{P}}^4}{\sqrt{2}{c}_{4\nu }\kappa M^2}\right)}^{1/3}·$$

(28)

Note that ${\sigma }_{\max }$ is independent of $C_{\mathrm{SUGRA}}$ and ${\sigma }_{\min }$ of $C_{\mathrm{RC}}$. The attainment of successful hilltop FHI requires the establishment of the hierarchy ${\sigma }_{\mathrm{c}}<{\sigma }_{★}<{\sigma }_{\max }$. From our numerical computation, we observe that, for constant $\kappa$ and ${\mathrm{a}}_S$, $n_{\mathrm{s}}$ decreases as ${\sigma }_{★}$ approaches ${\sigma }_{\max }$. To quantify these tuning values, we define the quantities

$${\Delta }_{\mathrm{c}★}={\sigma }_{★}/{\sigma }_{\mathrm{c}}-1\mathrm{and}{\Delta }_{\max \ast }=1-{\sigma }_{★}/{\sigma }_{\max }.$$

(29)

The naturalness of hilltop FHI increases with ${\Delta }_{\mathrm{c}★}$ and ${\Delta }_{\max \ast }$.

To obtain an impression of the required values of the parameters of the model, we first construct a benchmark table (

Table 2. A benchmark table of our scenario. We fix $\nu =7/8$ (resulting in $N=-49/8$) and $k=0.1$. Recall that $1\mathrm{PeV}={10}^6\mathrm{GeV}$, $1\mathrm{EeV}={10}^9\mathrm{GeV}$, $1\mathrm{ZeV}={10}^{12}\mathrm{GeV}$ and $1\mathrm{YeV}={10}^{15}\mathrm{GeV}$.

Benchmark			Benchmark
Point:	A	B	Point:	A	B
Inputs
ine $M\left(\mathrm{YeV}\right)$	$1.4$	$2.1$	$\kappa$	$0.0005$	$0.001$
$m\left(\mathrm{PeV}\right)$	$0.5$	$3.5$	$\lambda \left({10}^{-12}\right)$	$0.2$	$1.4$
Inflationary Parameters
${\mathrm{a}}_S\left(\mathrm{TeV}\right)$	$2.63$	$25.3$	${\Delta }_{\mathrm{c}★}(\%)$	$2.6$	$8.2$
$H_{\mathrm{I}}\left(\mathrm{EeV}\right)$	$0.23$	$1.05$	${\sigma }_{\max }/M$	$1.49$	$1.66$
${\sigma }_{★}/M$	$1.45$	$1.53$	${\Delta }_{\max \ast }(\%)$	$2.7$	$7.7$
$N_{\mathrm{I}★}$	$40.5$	$40.8$	${\sigma }_{\min }/M$	$35.5$	$30.8$
Observables
$n_{\mathrm{s}}$	$0.965$		$-{\alpha }_{\mathrm{s}}\left({10}^{-4}\right)$	$2.4$	$3.1$
r	$9·{10}^{-13}$	$1.8·{10}^{-11}$	$G{\mu }_{\mathrm{cs}}\left({10}^{-7}\right)$	$2.3$	6
z v.e.v and Particle Spectrum
During FHI			at the Vacuum $\langle z\rangle =1.4m_{\mathrm{P}}$
${\langle z\rangle }_{\mathrm{I}}\left({10}^{-3}{m}_{\mathrm{P}}\right)$	$1.3$	2	$m_{3/2}\left(\mathrm{PeV}\right)$	$0.9$	$6.2$
$m_{\mathrm{I}3/2}\left(\mathrm{TeV}\right)$	$1.2$	$11.2$	$m_z\left(\mathrm{PeV}\right)$	$1.2$	$8.8$
$m_{\mathrm{I}z}\left(\mathrm{EeV}\right)$	0.6	$2.7$	$m_{\theta }\left(\mathrm{PeV}\right)$	$0.9$	$5.6$
$m_{\mathrm{I}\theta }\left(\mathrm{EeV}\right)$	$0.08$	$0.5$	$m_{\mathrm{I}}\left(\mathrm{ZeV}\right)$	$1.7$	$5.2$
Reheating Process
$\mu /\tilde{m}$	$T_{\max }\left(\mathrm{PeV}\right)$		$\mu /\tilde{m}$	$-{\mathsf{\tau }}_{\max }$
3	$0.3$	$2.2$	3	$66.6$	67
$1/3$	$0.1$	$0.7$	$1/3$	$67.4$	$67.6$
$\mu /\tilde{m}$	$T_{\mathrm{rh}}\left(\mathrm{GeV}\right)$		$\mu /\tilde{m}$	$-{\mathsf{\tau }}_{\mathrm{rh}}$
3	$0.21$	$3.5$	3	$28.3$	$31.3$
$1/3$	$0.04$	$0.43$	$1/3$	$26.3$	$29.1$

). There, we display the results of our analysis for two benchmark points (BPs) with $\kappa =0.0005$ (BPA) and $\kappa =0.001$ (BPB) and $n_{\mathrm{s}}$ fixed at its central value in Equation (27). We also employ some representative $\nu$ and k values. In all cases, we obtain $N_{\mathrm{I}★}$∼40 from Equation (24). From the observables listed in Table 2, we infer that $|{\alpha }_{\mathrm{s}}|$ turns out to be of the order ${10}^{-4}$, whereas r is extremely tiny, of the order ${10}^{-11}$, and therefore far outside the reach of the forthcoming experiments devoted to detecting primordial GWs. From the entries in Table 2 related to ${\Delta }_{\mathrm{c}★}$ and ${\Delta }_{\max \ast }$, we notice that ${\Delta }_{\max *}>{\Delta }_{\mathrm{c}★}$; their values may be up to $10\%$ and increase as ${\mathrm{a}}_S$ or $\kappa$ (and M) increases. From the mass spectra arranged in Table 2 (see Section 4.1), we see that $m_{\mathrm{I}3/2}$ is similar to ${\mathrm{a}}_S$∼$\mathrm{TeV}$, whereas the other masses are of the order $\mathrm{EeV}$, whereas at the vacuum, $m_{\mathrm{I}}$ turns out to be of the order $1\mathrm{ZeV}={10}^{12}\mathrm{GeV}$, whereas $m_{3/2},m_z$ and $m_{\theta }$ lie in the $\mathrm{PeV}$ range; see Section 3. In the same table, we find it convenient to accumulate the values of some parameters related to the reheating process, specified in Section 5 below, and the formed CSs (see Section 6).

To further explore the parameter space of our model allowed by the inflationary requirements in Equations (24), (25) and (27), we present the gray-shaded regions in the $\kappa -{\mathrm{a}}_S$ and $M-{\mathrm{a}}_S$ planes; see Figure 1a,b. The boundaries of the allowed areas in Figure 1 are determined by the dashed [dot-dashed] lines corresponding to the upper [lower] bound on $n_{\mathrm{s}}$ in Equation (27); note that, in Figure 1a, we obtain a very narrow strip from the $n_{\mathrm{s}}$ variation. We also display the allowed contours for $n_{\mathrm{s}}=0.965$ by solid lines. The maximum r values are encountered in the upper-right end of the dashed lines—corresponding to $n_{\mathrm{s}}=0.974$. On the other hand, the maximum $|{\alpha }_{\mathrm{s}}|$ values are achieved along the dot-dashed lines, and the minimum value is ${\alpha }_{\mathrm{s}}=-3.2·{10}^{-4}$. Summarizing our findings from Figure 1 for a central $n_{\mathrm{s}}$ in Equation (27), we end up with the following ranges:

$$0.07\lesssim M/\mathrm{YeV}\lesssim 2.6\mathrm{and}0.1\lesssim {\mathrm{a}}_S/\mathrm{TeV}\lesssim 100.$$

(30)

Within the margins above, ${\Delta }_{\mathrm{c}★}$ ranges between $0.5\%$ and $9.5\%$ and ${\Delta }_{\max \ast }$ between $0.4\%$ and $8.2\%$m with the relevant values increasing with M and ${\mathrm{a}}_S$, as deduced from Table 2, too. The lower bounds of the inequalities above are expected to be displaced to slightly higher values due to the post-inflationary requirements (see Equation (37) below), which are not considered here for the sake of generality.

5. Reheating Process

Soon after FHI, the IS and z enter an oscillatory phase about their minima in Equations (9) and (11) and eventually decay via their coupling to lighter degrees of freedom. Note that $\theta$ and ${\theta }_S$ remain well stabilized at their values shown below Equation (9) during and after FHI, and so they do not participate in the phase of damped oscillations. The commencement of the z-dominated phase occurs for values of the Hubble rate $H_{z\mathrm{I}}$∼$m_z$. Given that $\langle z\rangle$∼$m_{\mathrm{P}}$ (see Equation (11)), the initial energy density of its oscillations ${\rho }_{z\mathrm{I}}$ is comparable to the energy density of the universe ${\rho }_{z\mathrm{It}}$ at the onset of these oscillations since

$${\rho }_{z\mathrm{I}}\sim m_z^2{\langle z\rangle }^2\mathrm{and}{\rho }_{z\mathrm{It}}=3m_{\mathrm{P}}^2{H}_{z\mathrm{I}}^2\simeq 3m_{\mathrm{P}}^2{m}_z^2.$$

(31)

Therefore, we expect that z will dominate the energy density of the universe until completing its decay through its weak gravitational interactions. Due to the weakness of these interactions, we expect that the reheating temperature $T_{\mathrm{rh}}$ will be rather low. This is the notorious cosmic moduli problem [61,62] that plagues the vast majority of SUGRA settings.

In our model, though, $T_{\mathrm{rh}}$ adequately increases thanks to the sufficiently large $m_z$ and $\mu$ originating from $K_{\mu }$ in Equation (6a). To show this fact, we estimate $T_{\mathrm{rh}}$ by [85,86]

$$T_{\mathrm{rh}}={\left(72/5{\pi }^2{g}_{\mathrm{rh}*}\right)}^{1/4}{\mathsf{\Gamma }}_{\delta z}^{1/2}{m}_{\mathrm{P}}^{1/2},$$

(32)

where $g_{\mathrm{rh}*}\simeq 10.75-100$ serves as the effective number of relativistic degrees of freedom at $T_{\mathrm{rh}}$. This is achieved for

$${\tilde{\mathsf{\tau }}}_{\mathrm{rh}}\simeq -\frac{2}{3}ln\frac{2}{5}\sqrt{3}{m}_{\mathrm{P}}{\mathsf{\Gamma }}_{\delta z}{\rho }_{z\mathrm{I}}^{-1/2},$$

(33)

where $\tilde{\mathsf{\tau }}=ln\left(R/R_{\mathrm{I}}\right)$, with R being the scale factor of the universe and $R_{\mathrm{I}}$ being its value at the onset of the z oscillations. Also, ${\mathsf{\Gamma }}_{\delta z}$ is the total decay width of the (canonically normalized) sgoldstino

$$\widehat{\delta z}=\langle K_{Z{Z}^{*}}^{1/2}\rangle \delta z\mathrm{with}\delta z=z-\langle z\rangle \mathrm{and}\langle K_{Z{Z}^{*}}\rangle ={\langle \omega \rangle }^{-2}$$

(34)

which predominantly includes contributions from its decay into pseudo-sgoldstinos, $\theta$ and electroweak Higgs, $H_u$ and $H_d$ via the kinetic terms $K_{X{X}^{*}}{\partial }_{\mu }X{\partial }^{\mu }{X}^{*}$, with $X=Z,H_u$ and $H_d$ [62,63,64,65] of the Lagrangian. In particular, we have

$${\mathsf{\Gamma }}_{\delta z}\simeq {\mathsf{\Gamma }}_{\theta }+{\mathsf{\Gamma }}_{\tilde{h}},$$

(35)

where the individual decay widths are found to be

$${\mathsf{\Gamma }}_{\theta }\simeq \frac{{\lambda }_{\theta }^2{m}_z^3}{32\pi m_{\mathrm{P}}^2}\sqrt{1-\frac{4m_{\theta }^2}{m_{3/2}^2}}\mathrm{with}{\lambda }_{\theta }=\frac{\langle z\rangle }{N}=\frac{4\nu -3}{\sqrt{6}\nu }\mathrm{and}{\mathsf{\Gamma }}_{\tilde{h}}=\frac{2^{4\nu -1}}{3^{2\nu -1}}{\lambda }_{\mu }^2\frac{{\omega }^2}{4\pi }\frac{m_z^3}{m_{\mathrm{P}}^2}{\nu }^{4\nu }.$$

(36)

From the expressions above, we readily recognize that ${\mathsf{\Gamma }}_{\delta z}$ is roughly proportional to $m_z^3/m_{\mathrm{P}}^2$, as expected for any typical modulus [61,62]. 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 $\widehat{\delta z}$ into $\tilde{G}$s, selecting $\nu >3/4$, which ensures $m_z<2m_{3/2}$; see Equation (15c). We do so in order to protect our setting from the so-called [87,88] moduli-induced $\tilde{G}$ problem, i.e., the possible late decay of the produced $\tilde{G}$, and problems with the abundance of the subsequently produced lightest SUSY particles—cf. Ref. [62].

If we enforce compatibility between the theoretical and observed values of light-element abundances predicted by BBN, we achieve Ref. [66] a lower bound on $T_{\mathrm{rh}}$, which, for large $m_z$∼$0.1\mathrm{PeV}$, entails

$$T_{\mathrm{rh}}\ge 4.1\mathrm{MeV}\mathrm{for}{\mathsf{B}}_{\mathrm{h}}=1\mathrm{and}{T}_{\mathrm{rh}}\ge 2.1\mathrm{MeV}\mathrm{for}{\mathsf{B}}_{\mathrm{h}}={10}^{-3},$$

(37)

where ${\mathsf{B}}_{\mathrm{h}}$ is the hadronic branching ratio. The bound above is a little softened for larger $m_z$ values.

Taking the $\kappa$ and $m_z$ values allowed by the inflationary part of our model in Section 4.4, we evaluate $T_{\mathrm{rh}}$ as a function of $\kappa$ and delineate the regions allowed by the BBN constraints in Equation (37); see Section 4.3 below. The results of such computation are displayed in Figure 2, where we illustrate the allowed contours in the $\kappa -T_{\mathrm{rh}}$ plane for $\nu =7/8$. This is an intermediate value in the margin selected here $(3/4-1)$. The boundary curves of the allowed region correspond to $\mu =\tilde{m}/3$ or ${\lambda }_{\mu }=0.22$ (dot-dashed line) and $\mu =3\tilde{m}$ or ${\lambda }_{\mu }=1.96$ (dashed line), whereas the solid line is obtained for $\mu =\tilde{m}$ or ${\lambda }_{\mu }=0.65$. Note that the relation between ${\lambda }_{\mu }$ and $\mu /\tilde{m}$ is given in Equation (16). We see that there is an ample parameter space consistent with the BBN bounds in Equation (37), depicted with two horizontal lines. Since the inflationary requirements increase the scale m with $\kappa$ and since m heavily influences $m_z$ and $T_{\mathrm{rh}}$ (see Equation (32)), $T_{\mathrm{rh}}$ increases with $\kappa$. The maximum value of $T_{\mathrm{rh}}$ for the selected $\nu$ is obtained for $\mu =3\tilde{m}$ and is estimated to be

$$T_{\mathrm{rh}}^{max}\simeq 19\mathrm{GeV}$$

(38)

Obviously, decreasing $\mu$ below $\tilde{m}/3$ causes ${\lambda }_{\mu }$, ${\mathsf{\Gamma }}_{\delta z}$ and $T_{\mathrm{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 maximum temperature $T_{\max }$ during this period is much larger than $T_{\mathrm{rh}}$, which can be better considered as the highest temperature during radiation domination [89]. To be quantitative, $T_{\max }$ can be calculated as [86]

$$T_{\max }={\left(3/8\right)}^{3/20}{12}^{1/4}{3}^{1/8}{m}_{\mathrm{P}}^{1/4}{\mathsf{\Gamma }}_{\delta z}^{1/4}{\rho }_{z\mathrm{I}}^{1/8}/g_{\mathrm{rh}*}^{1/4}{\pi }^{1/2}.$$

(39)

It is achieved [86] for $\tilde{\mathsf{\tau }}={\tilde{\mathsf{\tau }}}_{\max }=0.39\ll {\tilde{\mathsf{\tau }}}_{\mathrm{rh}}$. The large hierarchy between $T_{\max }$ and $T_{\mathrm{rh}}$ can be appreciated by their numerical values displayed in Table 2 for BPA and BPB. We see that $T_{\max }$∼$1\mathrm{PeV}$, whereas $T_{\mathrm{rh}}$∼$1\mathrm{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 $\mathsf{\tau }$, which is defined as a function of the redshift z (not to be confused with saxion field z) or the scale factor R as [90]

$$\mathsf{\tau }=-ln(1+\mathrm{z})=ln\left(R/R_0\right),$$

(40)

where the subscript 0 hereafter refers to present-day values. $T_{\max }$ is achieved at ${\mathsf{\tau }}_{\max }$, which can be found from

$${\mathsf{\tau }}_{\max }={\tilde{\mathsf{\tau }}}_{\max }-{\tilde{\mathsf{\tau }}}_{\mathrm{rh}}+{\mathsf{\tau }}_{\mathrm{rh}}\mathrm{since}{R}_{\max }/R_0=\left(R_{\max }/R_I\right)·\left(R_I/R_{\mathrm{rh}}\right)·\left(R_{\mathrm{rh}}/R_0\right)$$

(41)

taking into account Equations (33) and (39). For $\mathsf{\tau }\ge {\mathsf{\tau }}_{\mathrm{rh}}$, $\mathsf{\tau }$ can be found using entropy conservation through the relation

$$\mathsf{\tau }=-ln{\left(g_{s*}/g_{0s*}\right)}^{1/3}\left(T/T_0\right)$$

(42)

where $g_{s*}$ is the entropy effective number of degrees of freedom at temperature T, and $T_0=2.35·{10}^{-13}\mathrm{GeV}$. Its precise numerical values are evaluated by using the tables included in public packages [91,92] and assuming the particle spectrum of the MSSM. For example, at BBN ($T=2\mathrm{MeV}$), ${\mathsf{\tau }}_{\mathrm{BBN}}\simeq -23$, and at the time of radiation–matter equidensity, ${\mathsf{\tau }}_{\mathrm{eq}}\simeq -8.1$, corresponding to ${\mathrm{z}}_{\mathrm{eq}}=3387$. The values of ${\mathsf{\tau }}_{\mathrm{rh}}$ and ${\mathsf{\tau }}_{\max }$ are accumulated in Table 2 for some sample $T_{\mathrm{rh}}$ and $\mu /\tilde{m}$ ratios. We remark that ${\mathsf{\tau }}_{\max }$ values are similar for both $T_{\mathrm{rh}}$, but ${\mathsf{\tau }}_{\mathrm{rh}}$ is smaller for larger $T_{\mathrm{rh}}$.

The energy density of the oscillating z condensate ${\rho }_z$, radiation ${\rho }_{\mathrm{r}}$ and matter ${\rho }_{\mathrm{m}}$ can be evaluated as follows [85,90]:

$${\rho }_z={\rho }_{\mathrm{rh}}{e}^{-3(\mathsf{\tau }-{\mathsf{\tau }}_{\mathrm{rh}})},{\rho }_{\mathrm{r}}=\frac{{\pi }^2}{30}{g}_{*}{T}^4\mathrm{and}{\rho }_{\mathrm{m}}={\mathsf{\Omega }}_{\mathrm{m}0}{\rho }_c{e}^{-3\mathsf{\tau }},$$

(43)

where ${\mathsf{\Omega }}_{\mathrm{m}0}=0.311$, and the energy density at reheating is given by ${\rho }_{\mathrm{rh}}={\rho }_{\mathrm{r}}(T=T_{\mathrm{rh}})$. Taking advantage of the formulas above, we illustrate our cosmological scenario in Figure 3. In particular, we present the evolution of various $log{\rho }_i$—with ${\rho }_i$ normalized to ${\rho }_{\mathrm{c}}$ given in Equation (14)—as functions of $\mathsf{\tau }$, keeping, for simplicity, only the dominant component of the universal energy density for each $\mathsf{\tau }$. As regards $log{\rho }_i$ with $i=z$, we employ the parameters of BPB in Table 2 besides $\mu$, which is $\mu =5\tilde{m}/4$ for $T_{\mathrm{rh}}=1.5\mathrm{GeV}$ (dashed line) or $\mu =2\tilde{m}/5$ for $T_{\mathrm{rh}}=0.5\mathrm{GeV}$ (dot-dashed line). We also use solid lines to show $log{\rho }_i$, with $i=$ r, m. From the plot, it is apparent that our scenario is distinguishable from the standard scenario, according to which radiation domination commences after a high $T_{\mathrm{rh}}$∼$1\mathrm{EeV}$, corresponding to ${\mathsf{\tau }}_{\mathrm{rh}}$∼$-50$. We observe that $log{\rho }_z$ and $log{\rho }_{\mathrm{m}}$ have the same slope since both are proportional to $-3\mathsf{\tau }$, whereas $log{\rho }_{\mathrm{r}}$ is more stiff, since it is proportional to $-4\mathsf{\tau }$. In the plot, we also localize the position of ${\mathsf{\tau }}_{\max }$, which is the same for both $T_{\mathrm{rh}}$s, the two ${\mathsf{\tau }}_{\mathrm{rh}}$ values and the equidensity point ${\mathsf{\tau }}_{\mathrm{eq}}$. Also shown is the $\mathsf{\tau }$ value that corresponds to the CS decay ${\mathsf{\tau }}_{\mathrm{dc}}$ for $r_{\mathrm{ms}}^{1/2}=8$; see Section 6.2. It is located within the radiation-dominated era.

6. Metastable CSs and Early Matter Domination

The $U{\left(1\right)}_{B-L}$ breaking that occurs for $\sigma \simeq {\sigma }_{\mathrm{c}}$ produces a network of CSs that may be stable meta- or quasi- stable. This network has the potential to undergo decay via the Schwinger production of monopole–antimonopole pairs, thereby leading to the generation of a stochastic GW background. Below, we compute the tension of these CSs in Section 6.1 and their GW spectra in Section 6.2 under the assumption that they are metastable.

6.1. CS Tension

The dimensionless tension $G{\mu }_{\mathrm{cs}}$ of the $B-L$ CSs produced at the end of FHI can be estimated by [19,44]

$$G{\mu }_{\mathrm{cs}}\simeq \frac{1}{2}{\left(\frac{M}{m_{\mathrm{P}}}\right)}^2{ϵ}_{\mathrm{cs}}\left(r_{\mathrm{cs}}\right)\mathrm{with}{ϵ}_{\mathrm{cs}}\left(r_{\mathrm{cs}}\right)=\frac{2.4}{ln(2/r_{\mathrm{cs}})}\mathrm{and}{r}_{\mathrm{cs}}={\kappa }^2/8g^2\le {10}^{-2},$$

(44)

where we take into account that $(B-L)\left(\mathsf{\Phi }\right)=2$—cf. Ref. [45]. Here, $G=1/8\pi m_{\mathrm{P}}^2$ is the Newton gravitational constant, and $g\simeq 0.7$ is the gauge coupling constant at a scale close to M. For the parameters in Equation (30), we find

$$0.59\lesssim G{\mu }_{\mathrm{cs}}/{10}^{-8}\lesssim 9.2.$$

(45)

To qualify the result above, we single out the following cases:

(a) If the CSs are stable, the range in Equation (45) is acceptable by the level of the CS contribution to the observed anisotropies of the CMB, which is confined by Planck [93] in the range

$$G{\mu }_{\mathrm{cs}}\lesssim 2.4·{10}^{-7}\mathrm{at} 95\% \mathrm{c}.\mathrm{l}.$$

(46)

On the other hand, the results of Equation (45) are completely excluded by the recent PTA bound, which requires [29]

$$G{\mu }_{\mathrm{cs}}\lesssim 2·{10}^{-10}\mathrm{at} 95\% \mathrm{c}.\mathrm{l}.$$

(47)

(b) If the CSs are metastable, the explanation [29,40,41] of the recent NANOGrav15 data [25,26,27,28] on stochastic GWs is possible for

$$M\gtrsim 0.9\mathrm{YeV}\mathrm{and}\kappa \gtrsim 0.0003.$$

(48)

This is because the $G{\mu }_{\mathrm{cs}}$ values, obtained through Equation (44), from the ranges above are confined to a range dictated by the interpretation of recent data [29]:

$${10}^{-8}\lesssim G{\mu }_{\mathrm{cs}}\lesssim 2.4·{10}^{-4}\mathrm{for}8.2\gtrsim \sqrt{r_{\mathrm{ms}}}\gtrsim 7.5\mathrm{at} 95\% \mathrm{c}.\mathrm{l}.$$

(49)

where the metastability factor $r_{\mathrm{ms}}$ is the ratio of the monopole mass squared to ${\mu }_{\mathrm{cs}}$. Since we do not investigate monopole formation in our work, the last restriction does not impact our parameters. Moreover, the GWs obtained from CSs have to be consistent with the upper bound on their abundance ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ originating from the third advanced LVK observing run [71]:

$${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2\lesssim 6.96·{10}^{-9}\mathrm{for}f=32\mathrm{Hz},$$

(50)

where f is the frequency of the observation. At last, although not relevant for our computation, let us mention for completeness that ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ should be smaller than the limit on dark radiation, which is encoded in an upper limit on $\Delta N_{\mathrm{eff}}$ from BBN and CMB observations [94]

$${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2\lesssim \frac{7}{8}{\left(\frac{4}{11}\right)}^{4/3}{\mathsf{\Omega }}_{\gamma 0}{h}^2=5.6·{10}^{-6}\Delta N_{\mathrm{eff}}\mathrm{with}\Delta N_{\mathrm{eff}}\lesssim 0.28\mathrm{at} 95\% \mathrm{c}.\mathrm{l}.,$$

(51)

taking into account the TT, TE, EE + lowE + lensing + BAO datasets [83]. Here, ${\mathsf{\Omega }}_{\gamma 0}{h}^2\simeq 2.5·{10}^{-5}$ is the photon relic abundance at present, and we assume that ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ develops a flat shape in a broad range of f values such that the f dependence in Equation (51) is very weak.

6.2. GWs from Metastable CSs with Low Reheating

We focus here on case (b) of Section 6.1 and compute the spectra of the GWs 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 Section 4) has important ramifications for the shape of the spectra of GWs. This is because a decaying-particle-dominated era can be approximated by matter domination, which leads to spectral suppression at relatively large frequencies [67,68,69,70].

To verify this fact in the case of our model, we compute the emitted GW background at a frequency f following the standard formula of Ref. [70]:

$${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2=\frac{8\pi h^2}{3H_0^2}f{\left(G{\mu }_{\mathrm{cs}}\right)}^2\sum _{\mathrm{k}=1}^{\infty }{C}_{\mathrm{k}}\left(f\right)P_{\mathrm{k}}\mathrm{with}{H}_0=\sqrt{{\rho }_{\mathrm{c}}/3m_{\mathrm{P}}^2}.$$

(52)

Here, ${\rho }_{\mathrm{c}}$ is given in Equation (14), and $P_{\mathrm{k}}$ is the power spectrum of GWs emitted by the $\mathrm{k}\mathrm{th}$ harmonic of a CS loop. Assuming cusps as the main source of GWs, $P_{\mathrm{k}}$ is given by [70]

$$P_{\mathrm{k}}\simeq \mathsf{\Gamma }/\zeta (4/3){\mathrm{k}}^{-4/3}\mathrm{with}\mathsf{\Gamma }=50\mathrm{and}\zeta (4/3)=3.6.$$

(53)

In Equation (52), ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ is expressed as a sum of contributions from an infinite number of normal modes. On technical grounds, however, we take the sum up to ${\mathrm{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

$$\begin{array}{ccc}\hfill C_{\mathrm{k}}\left(f\right)& =& \frac{2\mathrm{k}}{f^2}\left({\int }_{{\mathsf{\tau }}_{\max }}^{{\mathsf{\tau }}_{\mathrm{rh}}}d\mathsf{\tau }{e}^{5\mathsf{\tau }}{H}_z{n}_{\mathrm{m}}\right.+{\int }_{{\mathsf{\tau }}_{\mathrm{rh}}}^{\min \{{\mathsf{\tau }}_{\mathrm{eq}},{\mathsf{\tau }}_{\mathrm{dc}}\}}d\mathsf{\tau }{e}^{5\mathsf{\tau }}{H}_{\mathrm{st}}{n}_{\mathrm{r}}\hfill \\ & +& {\int }_{{\mathsf{\tau }}_{\mathrm{eq}}}^{{\mathsf{\tau }}_{\mathrm{dc}}}d\mathsf{\tau }{e}^{5\mathsf{\tau }}{H}_{\mathrm{st}}{n}_{\mathrm{rm}}\mathsf{\Theta }({\mathsf{\tau }}_{\mathrm{dc}}-{\mathsf{\tau }}_{\mathrm{eq}})+\left.{\int }_{{\mathsf{\tau }}_{\mathrm{eq}}}^{{\mathsf{\tau }}_{\mathrm{dc}}}d\mathsf{\tau }{e}^{5\mathsf{\tau }}{H}_{\mathrm{st}}{n}_{\mathrm{m}}\mathsf{\Theta }({\mathsf{\tau }}_{\mathrm{dc}}-{\mathsf{\tau }}_{\mathrm{eq}})\right),\hfill \end{array}$$

(54)

where we use $\mathsf{\tau }$ as an integration variable—cf. Ref. [70]—taking into account Equation (40) and $d\mathrm{z}=-e^{-\mathsf{\tau }}d\mathsf{\tau }$. Also, ${\mathsf{\tau }}_{\max }$, ${\mathsf{\tau }}_{\mathrm{rh}}$ and ${\mathsf{\tau }}_{\mathrm{eq}}$ are given in Section 5, and ${\mathsf{\tau }}_{\mathrm{dc}}$ corresponds to the decay $\mathsf{\tau }$ of the CS network, which is estimated from

$${\mathsf{\tau }}_{\mathrm{dc}}=-ln\left({\left(70/H_0\right)}^{1/2}{\left(\mathsf{\Gamma }{\mathsf{\Gamma }}_{\mathrm{dc}}G{\mu }_{\mathrm{cs}}\right)}^{1/4}+1\right)\mathrm{with}{\mathsf{\Gamma }}_{\mathrm{dc}}=4G{\mu }_{\mathrm{cs}}{m}_{\mathrm{P}}^2{e}^{-\pi r_{\mathrm{ms}}}$$

(55)

being the rate per unit length of metastable CSs [49]. For the loop number density per unit length $n(\ell ,t)$—with mass dimension 4—we adopt the expressions [70]

$$\begin{array}{ccc}\hfill n_{\mathrm{r}}(\ell ,t)& =& \frac{0.18}{t^{3/2}{{\ell }_{\mathsf{\Gamma }}}^{5/2}}\mathsf{\Theta }(0.1t-\ell ),\hfill \end{array}$$

(56a)

$$\begin{array}{ccc}\hfill n_{\mathrm{m}}(\ell ,t)& =& \frac{0.27-0.45{(\ell /t)}^{0.31}}{t^2{\ell }_{\mathsf{\Gamma }}^2}\mathsf{\Theta }(0.18t-\ell ),\hfill \end{array}$$

(56b)

$$\begin{array}{ccc}\hfill n_{\mathrm{rm}}(\ell ,t)& =& \frac{0.18t_{\mathrm{eq}}^{1/2}}{t^2{{\ell }_{\mathsf{\Gamma }}}^{5/2}}\mathsf{\Theta }(0.09t_{\mathrm{eq}}-{\ell }_{\mathsf{\Gamma }}-\ell ),\hfill \end{array}$$

(56c)

where the subscripts r and m refer to CSs produced and radiating in radiation- and matter-dominated eras, respectively, whereas rm means that the loops were produced during radiation domination but were radiating during matter domination. Also, ${\ell }_{\mathsf{\Gamma }}$ is the initial length of a loop, which has a length ℓ at a cosmic time t and is given by

$${\ell }_{\mathsf{\Gamma }}=\ell +\mathsf{\Gamma }G{\mu }_{\mathrm{cs}}t\mathrm{with}\mathsf{\Gamma }\simeq 50\mathrm{and}\ell =2\mathrm{k}{e}^{\mathsf{\tau }}/f.$$

(57)

Here, $\mathsf{\Gamma }$ is related to the energy emission from CSs [67,68,69] and introduces some uncertainty in the computation. The Hubble parameter during standard cosmological evolution, $H_{\mathrm{st}}$, and during z oscillations, $H_z$, can be found from the formulas [85,86,89,90]

$$H_{\mathrm{st}}=\frac{1}{\sqrt{3}{m}_{\mathrm{P}}}{\left({\mathsf{\Omega }}_{\mathsf{\Lambda }}{\rho }_{\mathrm{c}0}+{\rho }_{\mathrm{m}}+{\rho }_{\mathrm{r}}\right)}^{1/2}\mathrm{and}{H}_z={\rho }_z^{1/2}/\sqrt{3}{m}_{\mathrm{P}},$$

(58)

where the various ${\rho }_i$ values are given in Equation (43). Lastly, cosmic time as a function of $\mathsf{\tau }$ is written as

$$\begin{array}{ccc}\hfill t\left(\mathsf{\tau }\right)& =& {\int }_{{\mathsf{\tau }}_{\max }}^{\mathsf{\tau }}\frac{d{\mathsf{\tau }}^{\prime }}{H_z}\mathsf{\Theta }(\mathsf{\tau }-{\mathsf{\tau }}_{\max })\mathsf{\Theta }({\mathsf{\tau }}_{\mathrm{rh}}-\mathsf{\tau })+{\int }_{{\mathsf{\tau }}_{\max }}^{\mathsf{\tau }}\frac{d{\mathsf{\tau }}^{\prime }}{H_{\mathrm{st}}}\mathsf{\Theta }(\mathsf{\tau }-{\mathsf{\tau }}_{\mathrm{rh}})\hfill \\ & \simeq & \left\{\begin{array}{cc}2/3H_z\hfill & \mathrm{for}{\mathsf{\tau }}_{\max }\le \mathsf{\tau }<{\mathsf{\tau }}_{\mathrm{rh}},\hfill \\ 2/3(1+w)H_{\mathrm{st}}\hfill & \mathrm{for}\mathsf{\tau }>{\mathsf{\tau }}_{\mathrm{rh}}.\hfill \end{array}\right.\hfill \end{array}$$

(59)

In particular, we have $w=1/3$ or $w=0$ for ${\mathsf{\tau }}_{\mathrm{rh}}<\mathsf{\tau }\le {\mathsf{\tau }}_{\mathrm{eq}}$ or $\mathsf{\tau }>{\mathsf{\tau }}_{\mathrm{eq}}$, respectively.

Armed with the formulas above, we compute ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ for the GWs produced from the CSs formed in our setting under the assumption that they are metastable. Employing the inputs of BPB in Table 2—which yield $G{\mu }_{\mathrm{cs}}=6·{10}^{-8}$—we obtain the outputs displayed in Figure 4. In particular, in Figure 4a, we show ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ as a function of f for $r_{\mathrm{ms}}^{1/2}=8$ and $\mu =\tilde{m}/3$, which yields $T_{\mathrm{rh}}=0.4$ (dotted line) or $\mu =3\tilde{m}$, which results in $T_{\mathrm{rh}}=3.5$ (solid line). On the other hand, for the GW spectra depicted in Figure 4b, we employ $\mu =\tilde{m}$, resulting in $T_{\mathrm{rh}}=1.2\mathrm{GeV}$, and fix $r_{\mathrm{ms}}^{1/2}$ at $7.9$ (dotted line) 8 (solid line) and $8.2$ (dashed line); see Equation (49). In both panels of Figure 4, we see that the derived GW spectra can explain NANOGrav15 data shown with gray, almost vertical lines. We see, though, that, as $r_{\mathrm{ms}}^{1/2}$ increases, the increase in ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ becomes sharper and provides a better fit to the observations. Also, in both panels, the shape of the GW signal suffers a diminishment above a turning frequency $f_{\mathrm{rh}}\sim 0.03\mathrm{Hz}$, which enables us to satisfy Equation (50) more comfortably than in the case with high reheating—cf. Refs. [35,36,38,39]. As $T_{\mathrm{rh}}$ decreases, the reduction in ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ becomes more drastic, in accordance with the findings of Refs. [67,68,69,70]. The plots also show examples of sensitivities of possible future observatories [95,96,97,98,99,100,101,102,103,104], which can test the signals at various f values. Needless to say, the bound in Equation (51) is not depicted since it lies above the ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ values displayed in the plots.

It would be preferable to obtain larger M, and thus $G{\mu }_{\mathrm{cs}}$, values (e.g., $M\simeq 10\mathrm{YeV}$ yields $G{\mu }_{\mathrm{cs}}\simeq {10}^{-6}$) such that the resulting ${\mathsf{\Omega }}_{\mathrm{GW}}{h}^2$ enters the dense part of the NANOGrav15 favorable region. This can be achieved (see, e.g., Refs. [21,35,36]) if we use $\theta =\pi$ and a next-to-minimal version [2] for $K_{\mathrm{I}}$. However, the generation of ${\mathrm{a}}_S$ from Z with $\theta =0$ in Equation (21c) fixes the sign of the second term in Equation (21b) and does not allow for the alternative arrangement mentioned above. On the other hand, next-to-next-to-minimal $K_{\mathrm{I}}$ may accommodate [2,3,13,14,21] such an augmentation of the M value without disturbing the z stabilization during FHI: see Section 4.1. In that case, we expect that ${\mathrm{a}}_S$ would not be constrained by the requirements of successful FHI, and it could be selected so that $T_{\mathrm{rh}}$ lies at a convenient level that allows for the comfortable evasion of the exclusion limits from the LVK [71] experiment. This setting opens up an interesting interconnection between the GW experiments and $\tilde{m}$.

7. Predictions for the SUSY Mass Scale

The ${\mathrm{a}}_S$ values displayed in Figure 1 and Equation (30)—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 $\tilde{m}=m_{3/2}$; see Equation (16). This aim can be achieved by solving Equation (21c) with respect to m as follows:

$$m\simeq {\left(\frac{{\mathrm{a}}_S}{2^{1+\nu }(2-\nu )}\right)}^{(2-\nu )/2}{\left(\frac{3H_{\mathrm{I}}^2}{(1-\nu ){\nu }^2}\right)}^{\nu /4},$$

(60)

where we take into account Equation (18) and the fact that ${\langle z\rangle }_{\mathrm{I}}/m_{\mathrm{P}}$∼${10}^{-3}$. As a consequence, the analytic result above well approximates the numerical one, which is obtained by consistently extracting m as a function of $\kappa$ and M, determined by the conditions in Equations (24) and (25). Note that an iterative process has to be realized while introducing a trial m value, which allows us to use the form of $V_{\mathrm{I}}$ in Equation (20) as input. The aforementioned smallness of ${\langle z\rangle }_{\mathrm{I}}$ causes the diminishment of m compared to ${\mathrm{a}}_S$, therefore rendering $\tilde{m}$ via Equations (16) and (15a) on the order of the $\mathrm{PeV}$ scale. Indeed, for fixed $\nu$, Equation (60) yields m, and then $m_{3/2}$, $m_z$ and $m_{\theta }$ can be easily estimated by Equations (15a) and (15c), whereas $\tilde{m}$ and $T_{\mathrm{rh}}$ can be obtained by Equations (16) and (32). Their numerical values are presented in BPs A and B in Table 2. The magnitude of the derived $\tilde{m}$ values and the necessity for $\mu$∼$\tilde{m}$, established in Section 5, hint toward the high-scale MSSM.

To numerically highlight our expectations above, we display the allowed (gray-shaded) region in the $\kappa -\tilde{m}$ plane, fixing $n_{\mathrm{s}}=0.965$ and varying $\nu$ and $\mu$ within their possible respective margins, $(0.75-1)$ and $(1/3-3)\tilde{m}$; see Figure 5. Along the solid line, we set $\nu =7/8$. From Equation (60), we can convince ourselves that the lower boundary curve of the displayed region is obtained for $\nu \simeq 0.751$, whereas the upper one corresponds to $\nu \simeq 0.99$. Assuming $\mu =\tilde{m}$ as well, we can determine the slice of the area that can be excluded due to the BBN bound in Equation (37). Overall, we find that $\tilde{m}$ turns out to be confined to the range

$$1.2\lesssim {\mathrm{a}}_S/\mathrm{TeV}\lesssim 100\mathrm{and}0.34\lesssim \tilde{m}/\mathrm{PeV}\lesssim 13.6,$$

(61)

whereas $T_{\mathrm{rh}}^{max}\simeq 71\mathrm{GeV}$, attained for $\mu =3\tilde{m}$ and $\nu \simeq 0.99$. The derived allowed margin of $\tilde{m}$ and the employed $\mu$ values render our proposal compatible with the mass of the Higgs boson discovered at the LHC if we adopt the high-scale version of the MSSM as a low-energy effective theory [73]. Indeed, within high-scale SUSY, the updated analysis requires [73]

$$3\mathrm{TeV}\lesssim \tilde{m}\lesssim 0.3\mathrm{ZeV},$$

(62)

for a degenerate sparticle spectrum, $\tilde{m}/3\le \mu \le 3\tilde{m}$, $1\le tan\beta \le 50$ and varying the stop mixing. On the other hand, our setting does not fit well with natural [62] or split [73] SUSY, which assumes $\mu \ll \tilde{m}$.

8. Conclusions

We analyzed the implementation of FHI together with SUSY breaking and CS formation in a B-L extension of the MSSM. We adopted the superpotential and Kähler potential in Equations (1) and (3), applying an approximate R symmetry. The model offers the following interesting achievements:

• Observationally acceptable FHI of the hilltop type, adjusting the tadpole parameter ${\mathrm{a}}_S$ and the $B-L$-breaking scale M in the ranges of Equation (30).
• A prediction of the SUSY mass scale $\tilde{m}$, which turns out to be on the order of PeV.
• An interpretation of the DE problem without extensive tuning. We obtain $\lambda$∼${10}^{-12}$ in Equation (2c).
• An explanation of the $\mu$ term of the MSSM with $\left|\mu \right|$∼$m_{3/2}$ (see Equation (6a)) by invoking an appropriate modification of the Giudice–Masiero mechanism.
• Reheating generated due to the domination of the sgoldstino condensate after the end of FHI. Since $\mu$ is on the order of PeV, the resulting $T_{\mathrm{rh}}$ is higher than its lower bound from BBN.
• $B-L$ CSs that, if metastable, explain the recent observations on GWs rather well. A characteristic of the obtained spectra is their suppression at relatively large frequencies ($f>0.1\mathrm{Hz}$), which is due to the long-lasting matter domination caused by sgoldstino oscillations.
• The generation of neutrino masses via the type I seesaw mechanism supported by $W_{\mathrm{MD}}$ and $K_{\mathrm{D}}$ in Equations (2e) and (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\mathrm{PeV}$) is lighter than $N_i^c$, which acquires mass on the order of $1\mathrm{ZeV}$—cf. Refs. [19,24]—and so the direct decay of z into $N_i^c$ is forbidden. However, there are extensions of the MSSM [105] where the late decay of the sgoldstino may non-thermally generate the baryon asymmetry of the universe. Alternatively, this problem may be overcome by applying improved approaches [106] based on the idea of cold electroweak baryogenesis [107,108]. As regards CDM, the candidacy of the lightest neutralino has to be investigated thoroughly by precisely solving the relevant Boltzmann equations, as in Refs. [85,86,109,110]. If the abundance is low, we can slightly open up the decay channel of z into $\tilde{G}$ so that we obtain a controllable production of neutralinos. Another aspect is the fate of the R axion, which remains stable throughout our setting [111,112].

As regards naturalness, it is a puzzle why SUSY should appear at such a scale, which is higher than the electroweak scale, making the fine-tuning severe, while it is much smaller than the fundamental energy scales, such as the GUT or Planck scales. Our scenario provides a possible solution to this issue: this may be due to the inflationary selection. In particular, the apparent fine-tuning could be a result of the combination of FHI and a bias toward high-scale SUSY in the landscape.