Ripples in Spacetime from Broken Supersymmetry

We initiate the study of gravitational wave (GW) signals from first-order phase transitions in supersymmetry-breaking hidden sectors. Such phase transitions often occur along a pseudo-flat direction universally related to supersymmetry (SUSY) breaking in hidden sectors that spontaneously break $R$-symmetry. The potential along this pseudo-flat direction imbues the phase transition with a number of novel properties, including a nucleation temperature well below the scale of heavy states (such that the temperature dependence is captured by the low-temperature expansion) and significant friction induced by the same heavy states as they pass through bubble walls. In low-energy SUSY-breaking hidden sectors, the frequency of the GW signal arising from such a phase transition is guaranteed to lie within the reach of future interferometers given existing cosmological constraints on the gravitino abundance. Once a mediation scheme is specified, the frequency of the GW peak correlates with the superpartner spectrum. Current bounds on supersymmetry are compatible with GW signals at future interferometers, while the observation of a GW signal from a SUSY-breaking hidden sector would imply superpartners within reach of future colliders.


Introduction
If supersymmetry is a property of our universe, how will it be discovered? Conventionally, searches for evidence of supersymmetry (SUSY) have focused on the Standard Model, looking for supersymmetric partners of Standard Model particles in direct production at colliders, scattering in dark matter experiments, and virtual effects in precision measurements. Thus far, no evidence has emerged of supersymmetry as it relates to the Standard Model, raising the prospect that it may lie outside the reach of the existing experimental program. Although this would pose a challenge to supersymmetry as a fully natural explanation for the scale of electroweak symmetry breaking, the abundance of remaining motivation (e.g. gauge coupling unification, dark matter, and straightforward string-theoretic embedding) favors continuing the search to shorter and shorter distances. While the LHC and proposed future colliders are promising tools in this search, the immense technical challenges of exploring energies far above the TeV scale in terrestrial experiments suggests casting a broader net. It invites identifying both new ways of accessing shorter distances and new sectors in which supersymmetry may be manifest.
A compelling avenue to shorter distances is to make use of the incredible energies of the Big Bang, searching for the imprint of supersymmetric phenomena on the early universe. In some sense, this is already the path taken by dark matter searches looking for the population of stable superpartners produced in the early universe, but it is not the only cosmological avenue for discovering SUSY. For example, spontaneous breaking of supersymmetry during inflation raises the prospect of observing signals in the threepoint function of primordial curvature perturbations [1,2], although the size of the signal depends on the strength of couplings between SUSY multiplets and the inflaton.
As for new sectors, at least one is guaranteed to exist in a supersymmetric universe: the sector responsible for breaking supersymmetry. Although there are many dynamical mechanisms for breaking supersymmetry, they typically possess a number of generic or universal features which can provide new ways of searching for supersymmetry even when superpartners of the Standard Model are decoupled. These include the goldstino, a goldstone fermion of spontaneous supersymmetry breaking (which becomes the longitudinal mode of the gravitino, the supersymmetric partner of the graviton, once gravity is accounted for), as well as a novel abelian global symmetry called the R-symmetry. The R-symmetry is generically spontaneously broken by the same dynamics that breaks supersymmetry, giving rise to a goldstone boson (the R-axion) and its scalar partner, a pseudo-modulus whose flat potential is protected by supersymmetry. In theories with low-energy supersymmetry breaking (LESB), in which the effects of SUSY breaking are communicated to the Standard Model by forces stronger than gravitation, these states may be accessible on their own. For example, the goldstino couples directly to Standard Model particles and may be produced at colliders, although the current reach of the LHC makes these searches less promising than continuing to look for Standard Model superpartners.
In this paper, we explore a new avenue for discovering supersymmetry in the physics of the early universe: using the stochastic gravitational wave background (SGWB) produced by a first-order phase transition to directly probe the sector responsible for breaking super-symmetry [3]. This makes use of the extraordinary opportunities afforded by the detection of gravitational waves (GW) by the LIGO-Virgo collaboration [4], which has opened a new era in the exploration of the early universe. Sensitivity of current and proposed GW interferometers to stochastic gravitational wave backgrounds broadly motivates identifying beyond-the-Standard Model (BSM) scenarios whose first-order phase transitions may generate such a signal and exploring the complementarity of GW interferometry with other probes of new physics such as present and future colliders. 1 Among the most compelling scenarios for SGWB are those in which a first order phase transition (FOPT) is associated with the breaking of a global or gauge symmetry in the early universe [7][8][9][10][11][12]. As we will show, the supersymmetry-breaking sector is a natural candidate for such a phase transition because it generically possesses at least one pseudoflat complex scalar direction, the pseudomodulus. In our constructions, the phase of this complex scalar direction is associated to the R-symmetry, which is in turn tied to the SUSYbreaking dynamics by many known theorems about SUSY quantum field theories [13][14][15]. This complex scalar direction is lifted by quantum corrections and the resulting potential is likely to possess a metastable minimum at the origin (where R-symmetry is preserved), which will then decay to the true minimum through a FOPT in the early universe. At the true minimum the R-symmetry is broken, consistently with a realistic SUSY spectrum featuring Majorana masses for the fermionic partners of Standard Model gauge bosons.
In our framework, the SUSY-breaking scale √ F correlates directly with the frequency range of the SGWB, such that theories of low-energy SUSY breaking feature a peak frequency accessible at LIGO-Virgo or proposed GW interferometers such as the Laser Interferometer Space Antenna (LISA), Einstein Telescope (ET), Cosmic Explorer (CE), DECi-hertz Interferometer Gravitational wave Observatory (DECIGO), and the Big Bang Observatory (BBO). In fact, a consistent cosmological history (in which the production of gravitinos in the early universe is consistent with the present dark matter abundance and small-scale structure constraints) guarantees that low-energy supersymmetry-breaking phase transitions produce a peak frequency in the range accessible to current and future interferometers [16][17][18][19].
Once a mechanism is specified to mediate supersymmetry breaking to the Standard Model, the scale √ F is also correlated with the spectrum of Standard Model superpartners, allowing the possibility of cross-correlating GW and collider signals. As we will see, the non-observation of SUSY particles at the LHC leaves open the opportunity for seeing SGWB signatures from low-energy SUSY breaking, making this a leading avenue for the discovery of supersymmetry. In return, the observation of a SGWB signal from a lowenergy supersymmetry breaking phase transition would imply the SUSY spectrum to be within the reach of future colliders such as FCC-hh, SPPC, or a high-energy muon collider, highlighting the strong complementarity between such SGWB signals and proposed colliders.
Along the way, we identify a qualitatively new class of natural potentials capable of generating large GW signals from a first-order phase transition, corresponding to the scalar potential along the pseudo-flat direction associated with SUSY breaking. The size of the vacuum energy gap between the metastable and the true vacuum is set by the SUSYbreaking scale √ F , which is necessarily smaller than the SUSY masses of other heavy fields in the sector in order to avoid tachyonic directions. As a consequence, the thermal corrections to the potential along the pseudomodulus direction are well described by a low-T expansion. Another distinctive feature of the pseudomodulus potential is the flatness at large field values, which is ensured by SUSY cancellations independently of the nature of SUSY-breaking deformations around the origin. These features combine to give strong first-order phase transitions, with the strongest transitions arising most naturally in models with two distinct SUSY-breaking scales. As we will discuss, a large amount of tuning would be necessary to realize a similar situation in non-SUSY scenarios, which explains why this possibility has not been explored so far in the literature (see for instance [20] for a collection of potentials giving raise to FOPT for the SM Higgs).
The organization of our paper is intended to highlight the qualitative connections between low-energy supersymmetry breaking, first-order phase transitions, and stochastic gravitational wave signals before progressing into explicit examples, and does not presume deep familiarity with spontaneous supersymmetry breaking. We begin in Sec. 2 by giving a broad overview of the phenomenology of low-energy SUSY breaking, the parametrics of gravitational wave signals from first-order phase transitions, and the relationship between the frequency of the SGWB signal and spectrum of SUSY particles. In Sec. 3 we discuss the general features of the pseudomodulus potential and the properties of a first-order phase transition along this direction, highlighting their novelty compared to commonly-studied potentials for FOPT. We present a simple toy model that captures the main features of concrete SUSY potentials, showing how a promising GW signal from FOPT requires multiple SUSY-breaking scales.
We then proceed to develop a series of increasingly realistic SUSY-breaking hidden sectors featuring FOPT in Sec. 4. In Sec. 4.1, we derive the phase diagram of the simple O'Raifeartaigh model. In Sec. 4.2 we present the simplest single scale model featuring a FOPT, which is a simple deformation of the O'Raifeartaigh model with explicit Rsymmetry breaking. Here, the FOPT takes place between the SUSY-breaking vacuum at the origin (which enjoys an unbroken R-symmetry) and the R-symmetry breaking minimum far away from the origin, where SUSY is restored unless coupled to an additional source of SUSY-breaking. In Sec. 4.3, we develop a fully realistic model by introducing gauge interactions to the O'Raifeartaigh model, such that SUSY is broken in both the metastable vacuum at the origin and true vacuum. The presence of both F -term and D-term SUSYbreaking naturally gives rise to strong GW signals.
In Sec. 5 we further comment on the phenomenology of our setup and the complementarity between GW observatories and colliders, highlighting the sense in which the observation of a SGWB signal in our models would ensure further evidence for SUSY at future colliders. We summarize our qualitative conclusions and future directions in Sec. 6. Technical details are reserved for a series of appendices, including a review of the one-loop thermal effective potential in App. A, approaches to the calculation of the bounce action 2) for definitions). The GW reach is computed by requiring the strength of the SGBW signal at the peak frequency to intersect with the PLI curve of a given GW interferometer (see Appendix C.1 for details). The colored regions show the reach for a signal generated from plasma waves which is generically the dominant one in our scenarios (see Sec. 2.2). The cyan region with β H < 10 imples a large fine-tuning in our setups (see discussion around Eq. (2.18), the parametric discussion in Sec. 3 and the explicit evaluation in the models of Sec. 4). In the red shaded region gravitino pair production is excluded by a γ + MET search at LEP with L = 0.24 fb −1 [21], the gray region is excluded by ATLAS bounds j + MET at √ s = 8 TeV and L = 10.5 fb −1 [22,23]. The dotted gray and dotted red lines are the projection of the γ + MET reach at FCC-hh with √ s = 100 TeV and L = 30 ab −1 and a future high energy lepton collider (HELC) with √ s = 30 TeV and L = 100 ab −1 . The dark green shaded region with dark green arrows indicates the bound on the SUSY-breaking scale derived from the LHC bound on gluinos mg > 2 TeV, requiring the messenger sector to be perturbative. The two dark green and light green bands show the impact of the present LHC bounds [24][25][26][27] and the future FCC-hh reach on gluinos [28] for perturbative messenger sectors with g M ∈ (0.01, 0.1) (see Eq. (2.10) for a definition of g M ). The region between these two lines will be naturally populated by the model discussed in Sec. 4.3 and Sec. 5.2. The dark blue arrows on the right hand side shows the ultralight gravitino window where m 3/2 ≤ 16 eV and the gravitino does not poses any cosmological challenge with κ = 1 (see Eq. (2.2) for a definition) and the gravitino dark matter window where κ 1 and the full gravitino mass is heavier than the gravitino mass contribution set by √ F . The dark cyan region marked as inaccessible in LESB is always excluded by a combination of gravitino overabundance [29] and BBN constraints [30] (see Sec. 5 for details).
in App. B and inputs to our projections for GW interferometers in App. C.

Detectable GW signals in low energy SUSY-breaking
In this section we illustrate the correlation between possible signals at present and future GW interferometers and the phenomenology of low energy SUSY-breaking (LESB). The underlying assumption is that the SGWB is produced through a first-order phase transition controlled by the SUSY-breaking hidden sector. As we will show, this connection is relatively insensitive to the details of the dynamics in the hidden sector. In Sec. 2.1 we summarize the structure and the parametric predictions of LESB theories, while in Sec. 2.2 we go through the field theory inputs that are necessary to compute the spectrum of GWs from a FOPT. In Sec. 2.3 we combine the results of the preceding sections to delineate the parameter space of possible gravitational wave signals from low-energy supersymmetry breaking, illustrated in Fig. 1.

Low-energy SUSY-breaking
Here we briefly review the structure of LESB and its broad parametric predictions, remaining agnostic as to the particular model realization. This general discussion is buttressed by Sec. 5, where we will present the predictions of a simple, explicit model which gives rise to GW signals. For the purposes of this paper, we adopt a phenomenological definition of low-energy SUSY breaking: LESB models are those in which the gravitino is the lightest supersymmetric particle (LSP). This requirement has deep implications for collider searches, precision observables, and cosmology, which we summarize in turn.
The first ingredient in low energy SUSY-breaking scenarios is a hidden sector at a high scale m * , which breaks supersymmetry (and R-symmetry) spontaneously. At energies much below m * , the spontaneous breaking of both supersymmetry and the R-symmetry can be encoded in a model-independent manner through the F -and scalar components of a single chiral superfield where f a is the order parameter for the R-symmetry breaking while √ F is the SUSY breaking scale, corresponding to an R-charge R X = 2 for the superfield X. The parameter θ is a constant, complex anti-commuting two-component spinor enabling component fields of different spin to be united into a single superfield. The Majorana fermion in the multiplet is the GoldstinoG, the goldstone fermion associated with spontaneous SUSY-breaking, while the compact scalar field a is the R-axion, the goldstone boson associated with spontaneous R-symmetry breaking.
Switching on gravity, the Goldstino becomes the longitudinal component of the gravitino via the super-Higgs mechanism [31] while the R-axion is lifted by an unavoidable explicit symmetry-breaking contribution arising from the fine-tuning of the cosmological constant [32,33]. The gravitino and R-axion masses can be written as where M Pl = 2.4 · 10 18 GeV is the reduced Planck scale and we have defined This reflects the fact that the gravitino mass is set by the sum of supersymmetry-breaking contributions from all sectors, corresponding to a total SUSY-breaking scale √ F 0 that may be larger than the scale √ F in the LESB sector under consideration (i.e. κ 1). Similarly, the R-symmetry breaking scale f a may exceed the scale of supersymmetry breaking √ F (i.e. R 1), as is often the case in calculable hidden sectors. In writing the gravity contribution to the R-axion mass in Eq. (2.3) we saturated the upper bound on the superpotential vacuum expectation value (VEV) [34]. In the presence of a possible explicit R-symmetry breaking term in the hidden sector / R , the R-axion mass will receive an extra contribution making the R-axion heavier than the superpartners of Standard Model fields and hence phenomenologically irrelevant. Most of the universal phenomenological predictions of low-energy SUSY breaking follow from the gravitino's role as the LSP [35]. First, the gravitino is the endpoint of every superpartner decay. In particular, the lifetime of the next-to-lightest supersymmetric particle (NLSP) is determined by its decay into the gravitino plus a Standard Model state, where c NLSP is an O(1) coefficient which depends on the particulars of the NLSP. Second, the gravitino may be directly produced in pairs with a rate controlled by dimension-eight contact operators suppressed by 1/F 2 in the m soft m 3/2 limit [21,36]. These operators lead to a total cross section at lepton colliders for pair production in association with a photon of the form where √ s is the beam energy, E min is the minimal photon energy, and θ min is the minimal photon angle with respect to the beam direction. Here we have expanded in E min √ s (see Ref. [21] for the full formula). A similar formula can be derived for σ(pp →GGj) as shown in Ref. [22]. Using these formulas and rescaling the Standard Model backgrounds to higher energies and luminosities, we may determine the sensitivity of missing energy searches at future colliders gravitino pair production; see Sec. 5 for details. These searches lead to projected direct constraints on the SUSY-breaking scale √ F as shown in Fig. 1. Finally, the stability of the gravitino LSP typically results in a cosmological hazard. This is the well known "gravitino problem" of LESB theories [37][38][39]. For sufficiently high reheating temperature (i.e. T r.h. > 45m 2 3/2 M Pl /M 2 3 , where M 3 is the soft mass of the gluino), the gravitino is in thermal equilibrium with the Standard Model bath. At freeze-out, the gravitino is still relativistic and its abundance is bounded from above by small-scale cosmological observables [40,41]. The latter imply m 3/2 16 eV, which corresponds to √ F < 260 TeV. Alternately, if the reheating temperature is low enough, the gravitino is never in equilibrium with the Standard Model bath but is typically overproduced by a combination of UV scattering contributions [42][43][44][45], freeze-in from the decays of superpartners [46], and decay of the NLSP relic abundance after freeze-out [47,48].
In order for SUSY-breaking sectors to generate sizable GW signals, the hidden sector needs to be reheated after inflation. Fixing the reheating temperature T r.h. = √ F and requiring the gravitino to not overclose the universe implies where C UV = 45 √ 5f 3 /(8π 13/2 g 3/2 * ) 4 · 10 −5 , g * 230, and f 3 18 encodes the thermal corrections as computed in Ref. [45]; α eff 0.01 is chosen to match the correct dark matter relic abundance in the pure Higgsino case [49]. Eq. (2.8) reflects a similar expression in Ref. [29], although here we have dropped the freeze-in contribution from superpartner decays because it is always subdominant compared to UV scattering contributions. For κ = 1, the only region where the gravitino is not overabundant for T r.h. = √ F corresponds to √ F < 260 TeV, while for κ 1 one can decouple the gravitino mass and push the SUSY-breaking scale to be as high as √ F 5 · 10 7 GeV. In this case, the upper bound is obtained by combining the overclosure bound, LHC bounds on Standard Model superpartner masses, and the BBN bounds on NLSP decays into the gravitino through the universal two-body decay in Eq. (2.6) as derived in Ref. [30]. This bound could slightly vary depending on the NLSP type and the detailed features of the spectrum, but this does not alter the primary message: requiring a reheating temperature T r.h. = √ F to obtain sufficiently strong gravitational wave signals implies a quite stringent upper bound on √ F as long as the gravitino is required to be the LSP.
Thus far, our discussion has not correlated the scale √ F of supersymmetry breaking with the mass spectrum of Standard Model superpartners. Supersymmetry breaking in the hidden sector is transmitted to the visible sector (which we will take to be the minimal supersymmetric Standard Model, or MSSM, in this paper) through a mediation mechanism. The simplest possibility is to assume that a certain number of messengers N mess in a given representation of the SM gauge group are coupled to the SUSY-breaking field X via the superpotential W mess = y mess XΦΦ. Given this coupling, the R-symmetry breaking scale f a controls the masses of the fermionic messengers, while the SUSY-breaking scale √ F gives an off-diagonal mass to the scalar messengers. The non-supersymmetric splitting between scalar and fermionic messengers is then transmitted to MSSM superfields via Standard Model gauge interactions. The resulting gaugino and squark masses are those of standard gauge mediation [35], where Cf (I) is the quadratic Casimir of the representation of the MSSM sfermionf under the Ith Standard Model gauge group and for simplicity we have considered messengers in the 5 +5 representation of SU (5). The additional coefficient s M 1 appearing in the gaugino masses accounts for the phenomenon of "gaugino screening" [15,50,51]. In the simple scenarios discussed here, the ratio between the gluino and squark soft masses M 3 /mq √ N mess s M 1, so that the most relevant collider bounds at current and future colliders can be framed purely in terms of the gluino mass, assuming the squarks to be decoupled and the lightest gaugino to be the next-to-lightest SUSY particle (NLSP).
Writing the R-symmetry breaking VEV as in Eq. (2.4), the final gluino mass can be simply written in terms of underlying parameters as where we have collected various coefficients into a model-dependent prefactor g M which encodes i) the suppression of the gaugino masses due to f a √ F (i.e. √ R 1), ii) the enhancement for N mess 1, iii) the gaugino screening controlled by s M 1, and iv) the relation of the gluino soft mass to its pole mass, correctly accounting for the one loop running of the gluino soft mass at low energies in the limit of heavy squarks [52,53].
Broadly speaking, Eq. (2.10) establishes an interesting relation between the SUSYbreaking scale and the present and future collider bounds on the gluino. Given the current LHC bound on gluino masses, which ranges between mg 2 − 2.5 TeV [24][25][26][27], Eq. (2.10) indicates the lowest values of the SUSY-breaking scale √ F consistent with data. Depending on the model, g M can span many orders of magnitude, but there are three parametric regimes of interest: • 1 g M 160, which is realized in strongly-coupled messenger sectors that are at the boundary of perturbativity. The upper bound on g M is indeed obtained by requiring the SM gauge couplings and y mess to be perturbative at the scale of the hidden sector.
• g M 1, which is realized in weakly-coupled messenger sectors if M mess √ F and the gaugino masses are not screened. The latter requirement requires non-trivial dynamics in the hidden sector, as shown in Ref. [15].
• g M 1, which is typical of models where the soft masses are suppressed compared to the SUSY-breaking scale because f a √ F and the gaugino masses may be further screened compared to the squark masses. As we will show in Sec. 5, this is the typical situation in simple, explicit setups featuring SGWB signals.
One of the most appealing features of LESB mediated via gauge interactions is the flavor-preserving nature of the MSSM superpartner spectrum. This is because the flavorblind contributions to superpartner masses transmitted by gauge interactions vastly exceeds the omnipresent, flavor-violating "gravity-mediated" contributions. However, these gravity-mediated contributions reflect all contributions to SUSY breaking, while the gaugemediated contributions reflect only the SUSY-breaking in the sector of interest. Thus when the gravitino mass is enhanced by κ 1, the contribution from gravity mediation increases relative to the contribution from gauge mediation, and may eventually run afoul of bounds on flavor violation. In particular, this implies that κ is bounded from below by bounds from flavor-changing neutral currents (FCNCs). For instance, considering the slepton contributions to µ → eγ [54,55] and the squark contributions to ∆m K [56,57] leads to the bounds where κ and R are defined in Eq. (2.4) and we set BR(µ → eγ) < 4.2 · 10 −13 and ∆m K = (3.479 ± 0.001) · 10 −12 MeV, asking for the squark contribution to be less then present experimental uncertainty. These constraints give a robust upper bound on the gravitino mass in our framework. Finally, even in the absence of flavor violating effects, the electric dipole moments arising from the relative phase between gaugino and higgsino masses can be probed in precision experiments such as ACME [58,59]. The current limit are already challenging a CP-violating phase of order 10 −2 with gauginos below the TeV scale and the future experimental program will sensibly improve this reach making it one of the most interesting indirect probes of LESB [60,61].

First order phase transitions and SGWB
Phase transitions in field theory are triggered by the nucleation of vacuum bubbles and their subsequent percolation in the space-time volume. The vacuum bubbles can be found in Euclidean signature as the stationary minimum-energy bounce solutions interpolating between the false and true vacuum [62,63]. In all cases we consider here, the thermal fluctuations will dominate so that the total decay rate per unit volume can be approximated as where S 3 is the 3 dimensional Euclidean action for the O(3)-symmetric bounce [64,65]. The decay rate encodes the probability of true vacuum bubbles to be nucleated in a spacetime region where the false vacuum dominates. The time evolution of the phase transition can be described in terms of different temperatures. First of all, a necessary condition for nucleation is that the universe reaches temperatures below the critical temperature T c , where the false and true vacuum are degenerate. At the nucleation temperature T n < T c , one bubble will nucleate per Hubble volume, corresponding to 2 where we assumed that the phase transition happens during radiation domination so that H 2 (T ) = π 2 g * T 4 90M 2 Pl and normalized g * to the number of the degrees of freedom in the MSSM.
In all the cases we will consider, the last term in Eq. (2.13) can be neglected together with the constant term, so that solving is always a good approximation. Since C(T ) is a slowly-varying function of T , we can further simplify this equation assuming C(T n ) C(T c ); this number is always going to be O(10 2 ) in the temperature range of interest. After one bubble per volume has been nucleated at T n , the bubbles expand to fill the space-time volume. The phase transition is considered to be completed at the percolation temperature T p , when a small fraction of the total volume remains in the false vacuum. For fast phase transitions like the ones discussed here, one can show that T p T n so that we can neglect this difference and take T n as the temperature at which the phase transition completes. This sets the relevant dimensionful scale controlling the frequency range of the SGWB spectrum.
The shape and amplitude of the SGWB spectrum strongly depends on the amount of energy released into GWs during the FOPT, the duration of the phase transition, and the behavior of the bubbles in the cosmic fluid. The two first ingredients can be easily quantified in terms of field theory data via the quantities where ∆V (T n ) is the potential energy difference between the true and the false vacuum at T n . The amount of energy released into GWs is quantified by α, the latent heat relative to the radiation energy density ρ R = π 2 g * T 4 30 [66]. The duration of the phase transition is quantified by β H , the inverse of the typical timescale of the transition normalized with respect to Hubble; it is defined under the assumption that the nucleation rate rises exponentially [64,65] as S(t) = e β H H(t)(t−tn) . Using the approximate nucleation condition in Eq. (2.14) we can write where S (T n ) C in order for the nucleation rate to rise as a function of time, and β H C ∼ 100 unless there is some measure of fine-tuning between the first and the second 2 The nucleation temperature is formally defined by the integral 1 = Tc Tn dT T

Γ(T )
H(T ) 4 , which is well approximated by Eq. (2.13) since Γ(T ) depends exponentially on the temperature. terms of the above expression. To evaluate the fine-tuning associated with β H in explicit models, we define where ∆ p i β H are the individual tunings with respect to the underlying parameters of the theory p i . As we will discuss in Sec. 3, the parametric dependence of the fine-tuning can be derived for our general class of models and then computed explicitly in the models of Sec. 4. As a result, obtaining β H < 10 would imply a large amount of fine tuning (in the sense of being a non-generic prediction of a given model). This is illustrated in Fig. 1, where it provides a meaningful bound on the parameter space of GW signals in LESB .
The dominant production mechanism of gravity waves during the first-order phase transition depends on the dynamics of the bubbles in the cosmic fluid. If the mean free path of the particles is much longer than the width of the bubble wall, the velocity of the wall v w can be determined by equilibrating the pressure on the bubble wall induced by the difference in potential energy ∆V with the friction forces exerted by the surrounding plasma [67,68]. The latter are induced by states whose mass changes in passing from the false to the true vacuum. For v w → 1, the total pressure can be derived in a quasi-classical approximation [68][69][70] and reads where the Lorentz gamma factor is γ = 1/ 1 − v 2 w and the leading-order plasma friction P LO depends on the change in the masses-squared ∆m 2 of all the states in the thermal bath [69]. Since ∆m 2 = m 2 true − m 2 false , the approximate expression in Eq. (2.19) is only valid when both γT m true and T m false . The first condition ensures that particles in the false vacuum have enough energy to pass through the wall, while the second forestalls Boltzmann suppression of the pressure [68]. The next-to-leading order radiation pressure P NLO is instead induced by the change in mass of the vector bosons and it is γ-enhanced for v w → 1, as first derived in [70].
The pressure in Eq. (2.19) determines both how much the bubble wall accelerates as a function of the bubble radius [71,72], and the fraction of the FOPT energy which is in the bubble wall at the time of collision T * (traditionally called k coll ). Since we will be dealing with fast phase transitions, we take T * T p T n .
In the absence of friction, the acceleration of the bubble wall grows linearly with the bubble radius until the gamma factor reaches a terminal value where we took the initial radius to be R 0 R c = 3 2π S 3 (Tn) ∆V (Tn) 1/3 estimated in the thin wall approximation [62], estimated R * as in Ref. [73], and assumed radiation domination. The last term is O(1) in phase transitions which do not have a supercooling phase since ∆V (T n )/T 4 n < 75.6(g * /230).
In the FOPTs discussed here, the bubble growth is generically stopped by plasma effects from heavy states. This is due to a novel effect in which the relevant energy scale for particles interacting with the wall reaches values ∼ γT n much larger than the intrinsic scales associated with the bubble. This is particularly relevant for SUSY-breaking hidden sectors, where there is a large separation of scales that can be spanned by these ultrarelativistic effects. In particular, the bubbles expand and accelerate linearly with the radius until the boost factor is large enough to allow heavy states of mass m true T n to cross the bubble wall. The significant mass change of these states induces a new source of LO friction, This effect is very similar to the pressure term from mixing discussed in Ref. [74], but here we typically pay the Boltzmann suppression of m false . The resulting fraction of the energy in the bubble wall at the time of collisions is generically very suppressed, As such, most of the energy released in the FOPT goes into the plasma, giving rise to sound waves propagating through the cosmic fluid. These sound waves source gravitational waves from the motion of the plasma with an efficiency determined by where we have expanded the general formula of Ref. [75,76] for k coll 1. The resulting GW spectral density is

Pl
(1 + α) to account for the reheating of the plasma and β * H is normalized accordingly following Eq. (2.16). The sound wave spectrum is a broken power law which drops like Ω * sw ∼ f 3 for f f * sw , as expected from causality in a radiation dominated universe, and as Ω * sw ∼ f −4 for f f * sw . The high frequency behavior of the spectrum is likely to be affected by the turbulence contribution, whose size is still subject to large theoretical uncertainties [11,77]. Here, we include for simplicity only the sound waves contribution to the GW spectrum in Eq. (2.25), which will mainly determine the  Figure 2. The reach of future GW interferometers in the (α, β H ) plane for two different scales of FOPTs, assuming the signal is dominated by sound waves given by Eq. 2.25. The shaded regions are obtain by requiring the signal at the peak in Eq. (2.28) to be inside the PLI curve of a given experiment. Left: T n = 10 5 GeV corresponds to LESB scenarios with an ultralight gravitino LSP and κ = 1. Right: T n = 10 7 GeV corresponds to LESB scenarios with gravitino DM and κ 1. We will exhibit calculable scenarios of this type in Sec. 4. detectability of a given GW signal. After redshift is taken into account, assuming that the entropy per comoving volume remains constant [66], the GW spectrum today reads where the peak frequency and the power at the peak frequency scale as . (2.28) Here we have taken T * T n and normalized the scalings for α = 0.3, β H = 100 and T n = 10 7 GeV, which will be the typical values for FOPTs related to fully calculable SUSY-breaking hidden sectors explored in the following sections.
Having derived the expected GW spectrum, we can determine the region in the (α, β H ) plane where we expect the SBGW to be detectable at future interferometers. Given the fraction of energy density in GWs today in Eq. (2.26), the sensitivity of a given interferometer is controlled by the time integrated signal-to-noise ratio where Ω noise (f ) is the effective noise of the interferometer within a given frequency band (f min , f max ) and t obs is the observation time. A detectable stochastic GW background is defined to have ρ > 10. The Power Law Integrated (PLI) curves are generated by considering a power law function of the frequency f for the GW signal shape in Eq. (2.29). The PLI curves for each GW interferometer considered here are given in Appendix C for completeness.
In Fig. 2 we show the regions in the (α, β H ) plane where the power at the peak frequency in Eq. (2.28) lies within the reach of future interferometers for two different nucleation temperatures. Low nucleation temperatures such as T n = 10 5 GeV can be probed over a wide frequency range depending on β H (i.e. the duration of the FOPT) while high nucleation temperatures such as T n = 10 7 GeV will be accessible only at future high frequency interferometers such as Advanced LIGO (A-LIGO) [16], the Einstein Telescope (ET) [17] and the Cosmic Explorer (CE) [18,19]. In the next section we show that T n ∼ √ F in our LESB scenarios, so that high nucleation temperatures in fully calculable SUSYbreaking scenarios correspond to superpartners lying out of the reach of the LHC.

LESB in the future: GW interferometers vs. colliders
We are now ready to establish a connection between the SGWB signals and SUSY-breaking phenomenology described in the previous two sections. The first step is to relate the nucleation temperature relevant for the SGWB signal to the scales in a SUSY-breaking hidden sector. As we will see, the nucleation temperature is essentially set by the SUSYbreaking scale √ F in our scenarios. Focusing on FOPT where the barrier between the false and the true vacuum is present at T = 0, S 3 is bounded from below by a constant and we can define T min as the temperature where is monotonic for T > T min , the solution of the equation above is unique. The nucleation temperature is then bounded from above by T c , where β H → ∞ and α in Eq. (2.15) is suppressed and dominated by d∆V (Tn) dT . It is further bounded from below by T min where β H → 0 and α is dominated by ∆V (T n ).
More importantly, T n can be directly related to the SUSY-breaking scale √ F which sets the size of the O(3)-symmetric bounce action. A simple way of seeing this is to note that the bounce action at T n is itself set by the scale of relevant features in the potential, where c 3 is a model-dependent function of the parameters controlling the shape of the potential which we assume to be temperature independent for simplicity (an approximation that is certainly justified if T n is close enough to T min ). Here we assume that c 3 /C ∼ O(1), an assumption that will turn out to be justified analytically in Sec. 3 and numerically in the explicit models of Sec. 4.
Having established a relation between T n and the SUSY-breaking scale √ F , we identify two different viable regions of the LESB parameter space satisfying the following simple requirements: • the gravitino is the lightest supersymmetric particle (LSP) as required by LESB, and • the reheating temperature T r.h. is as high as √ F to generate GW signals from the hidden sector. We take T r.h. = √ F in Fig. 2.11 to maximize the allowed parameter space.
The two viable regions satisfying the above requirements are Gravitino Dark Matter window: where 260 TeV < √ F 50 PeV and κ 1 so that the gravitino mass is larger than the nominal value set by the F -term in Eq (2.1). The upper bound on the SUSY-breaking scale is obtained by combining the constraints on gravitino overabundance in Eq. (2.8), BBN constraints on NLSP decays, and the LHC bound on the gluino mass mg > 2 TeV. The precise upper bound is potentially dependent on further model-building epicycles; the value here is meant to be indicative. In this window, the gravitino abundance can match the observed dark matter relic abundance today, while the soft masses are still dominated by the gauge mediation contributions in Eq. (2.9) so that flavor constraints are under control when Eq. (2.11) is satisfied. Perturbative gauge mediation models with f a √ F and gaugino screening will naturally live in the upper end of this window for √ F 1 − 50 PeV. As shown in Fig. 1, the future reach on gluinos at FCC-hh [28] could provide a further direct test of these models. Future interferometers in the LIGO frequency band such as A-LIGO [16], ET [17] and CE [18,19] have the unique opportunity to probe these scenarios as long as the thermal transition to the SUSY-breaking vacuum is associated with a sufficiently strong FOPT (see Fig. 2). In the rest of the paper, we discuss explicit scenarios of this type.
Ultralight gravitino window: where m 3/2 < 16 eV and √ F < 260 TeV. This region has no cosmological issues for κ = 1, but it requires g M 1 to satisfy the LHC bound on the gluino mass given the low SUSY-breaking scale (see Eq. (2.10) for definition and comments). A lower bound on the gravitino mass can be derived from direct searches for gravitino pair production at LEP in γ + MET and at the LHC in j + MET. As shown in Fig. 1, present direct bounds on the gravitino are not competitive with the bound on √ F obtained by requiring mg > 2 TeV and perturbativity in the messenger sector. The HL-LHC will not improve much on that. Future colliders -in particular, high energy lepton colliders (HELCs) -can drastically improve the reach on gravitino pair production and meaningfully probe this window even if MSSM superpartners remain inaccessible. As shown in Fig. 2, these scenarios can be probed across a wide frequency range by future GW interferometers depending on the strength and the duration of the FOPT. Building explicit calculable models in this window presents challenges [78,79], and we leave a study of possible GW signals for a future work.

Anatomy of the SUSY-breaking phase transition
In this section we describe the generic features of FOPT occurring in calculable SUSYbreaking hidden sectors. First, we discuss how a large class of perturbative hidden sectors can be encoded in the effective field theory of the universal pseudomodulus, which is the scalar component x of the chiral superfield X in Eq. (2.1), universally related to the spontaneous breaking of supersymmetry [13][14][15].
Second, we show how the flatness of the pseudomodulus potential gives rise to a new class of FOPTs with a very distinctive feature: the nucleation temperature is generically small compared to the SUSY mass scale, T n ≤ m * , so that the thermal potential is well approximated in the low-T expansion. As we will discuss, non-supersymmetric realizations of this class of FOPT typically entail a large amount of fine-tuning.
Finally, we derive parametric estimates for T n , α and β H for this new class of FOPTs using the triangular barrier approximation [80,81] and comment on a universal feature of bubble dynamics in our FOPTs. The observations of this section will find a concrete realization in the working examples of Section 4.

The SUSY-breaking pseudomodulus
The existence of flat directions is a trademark of hidden sectors with spontaneous SUSY breaking. Here we focus on a large class of SUSY breaking sectors where the dynamics of both SUSY and R-symmetry breaking can be embedded in a single chiral superfield X parametrized as in Eq. (2.1) where the R-charges of the components are respectively The scalar component x (the universal pseudomodulus) tracks the breaking of the R-symmetry, while F sets the SUSY breaking scale. 3 The phase transition occurs along x from a local minimum at the origin x = 0 (where R-symmetry is preserved) to the T = 0 vacuum of the theory where x = f a and R-symmetry is broken. Hence x = f a is the order parameter of the phase transitions of interest here, parameterizing the spontaneous breaking of the R-symmetry.
In hidden sectors which admit a weakly coupled description, the phase transition can be fully described by studying the effective potential of the pseudomodulus x, whose mass is typically well below the mass m * of the heavy SUSY states in the hidden sector. As we will see, the unique features of the pseudomodulus potential leave a strong imprint on the properties of the phase transition. The full effective potential for the pseudomodulus can be written as Figure 3. Qualitative features of the pseudomodulus potential relevant to the FOPT in SUSYbreaking hidden sectors. Left: Sketch of the zero-temperature potential as described in Sec. 3.1, exhibiting the following features: i) the distance between the two minima is larger than their potential difference, f 4 a ∆V , and ii) the height of the peak between the two minima is loopsuppressed compared to the potential difference, V P

2)
∆V . An explicit realization of this potential is presented in Sec. 3.4. The tree level potential (dashed blue) generated by explicit R-symmetry breaking destabilizes the origin, giving rise to a minimum at x = f a where the R-symmetry is further spontaneously broken by the VEV of x. Quantum corrections (dashed red) generate a local minimum at origin. Right: Behavior of the temperature corrections described in Eq. (3.5) at T = 0, T = T c , and T = T n . The thermal corrections give a contribution to the potential at the origin which at T n is typically much smaller than F 2 . The barrier and the true vacuum are essentially unchanged. The approximations in Sec. 3.2 are then justified.
where V 0 (x) encodes the zero-temperature quantum corrections and V T (x) the thermal ones.
The zero-temperature part of the effective potential V 0 (x) is flat at tree level, up to explicit R-symmetry breaking effects. Along this so-called F -flat direction, the size of the potential energy is set by supersymmetry breaking, V ∼ F 2 . Interactions that explicitly violate the R-symmetry typically destabilize the origin and give a slope to the pseudomodulus potential at tree level, but these features are usually small compared to the scale √ F . At one loop, quantum corrections lift the pseudomodulus potential; these corrections are present even in the absence of explicit R-symmetry breaking. The combination of tree-level explicit R-symmetry breaking and one-loop quantum corrections give rise to the schematic zero-temperature potential shown in Fig. 3. Assuming the quantum corrections exceed the R-symmetry breaking effects, at zero temperature this creates a metastable vacuum at the origin that is separated by a barrier from the true vacuum at x true = f a . The energy difference between the two vacua ∆V is proportional to the SUSY-breaking scale. The barrier is located at a distance x P from the origin; at this point, the barrier height is V P . The essential features characterizing the zero temperature potential are: • The potential is flat. This means that the distance f a in field space between the false vacuum and the true vacuum is larger than the size of the potential energy difference ∆V : where this hierarchy assumes that R-symmetry breaking effects are parametrically smaller than the loop corrections. This will be manifest in the toy model of Sec. 3.4. Under this assumption, the flatness of the potential is a direct consequence of the fact that SUSY loop corrections asymptote to a logarithm at large field values (see for instance [83]). Obtaining similar quantum corrections in non-supersymmetric theories with a field-independent mass gap is notoriously difficult without fine-tuning. 4 • The barrier between the two vacua is small. Given that the potential is generated by loop effects (and subleading R-symmetry breaking effects), the size of the barrier V P is one-loop suppressed with respect to the energy difference between the true and the false vacuum ∆V : where λ eff ∼ O(1) should be thought of as the effective coupling determining the height of the barrier. The position of the barrier x P is model-dependent, but will not play a critical role in the determination of the bounce action as long as Eq. (3.3) is satisfied.
We now turn to the finite temperature corrections. First, as it is well known, finite temperature effects break SUSY and thus significantly modify the pseudo-modulus potential. The thermal effects are dominated by the loops of heavy fields in the SUSY hidden sector coupled to x, whose mass is of order m * . Since m * is by construction larger than the SUSY-breaking scale √ F setting the zero-temperature potential, the relevant temperatures for the phase transition are smaller than the mass scale m * of the particles running in thermal loops. This implies that the correct approximation of the thermal potential is the low temperature expansion (see Appendix A for explicit formulas). This makes the finite-temperature potential of the pseudomodulus qualitatively different from ordinary non-SUSY models, where the high temperature approximation applies since the typical scalar potential curvature is of the same order of the highest mass scale in the theory.
The effect of thermal corrections on the pseudomodulus potential takes the schematic form where we assumed the presence of N degrees of freedom with masses-squared ∼ λ 2 x 2 + m 2 * . Notice that N counts all the heavy degrees of freedom, both bosonic and fermionic, which contribute with the same sign to the thermal potential. This enhances the importance of thermal effects compared to zero-temperature loops, where cancellations occur between states of different statistics.
The thermal correction constitutes a negative contribution to the potential which is maximal (in absolute value) at the origin of the pseudomodulus, when x ∼ 0 and the Boltzmann suppression factor is minimized. As a consequence, thermal corrections in our scenarios have an exponentially larger impact at the origin relative to the true vacuum or the barrier. This behavior is explicitly shown in the right panel of Fig. 3. As we will show in Sec. 4 small deviations from this generic feature can be induced by heavy states becoming lighter at large field values of the pseudomodulus.

First order phase transitions in the low-T expansion
Given the shape of our potential as shown in Fig. 3, we can approximate the bounce action in the triangular barrier approximation [80,81]. Within this approximation, we will be able to capture the parametric behavior of the FOPTs analytically and in Sec. 4 we will show how the full analytical solution in explicit models reflects the general features explored here. A more in-depth discussion about the computation of the bounce action in the various cases can be found in Appendix B.
Taking the false vacuum to be at the origin of field space, we can write the bounce action in the triangular barrier approximaton as where we have defined the variables and the functions The expansion for r λ → 0 is justified as long as both Eq. (3.3) and Eq. (3.4) are satisfied and the true vacuum VEV f a sets the largest scale in the pseudomodulus potential. The triangular barrier approximation can be extended beyond the region set by f a /x P > g(r λ ), but the range of validity of Eq. (3.6) is sufficient to capture the parametrics of the phase transitions of interest. We give the full expression of the triangular barrier approximation in Appendix B.
The triangular barrier approximation depends in general on only five parameters characterizing the potential : the three values of the potential at the critical points, V ± , V P , and the position of the two critical points, f a , x P . For a given theory we can compute these temperature-dependent quantities explicitly, and find that the bounce action in Eq. (3.6) is an excellent match to the full numerical result. 5 At the leading order in the r λ → 0 expansion, the bounce action is independent of x P ; our analytical estimates will assume that this holds. To further simplify our analytical treatment, we approximate the thermal potential in Eq. (3.5) by only including thermal corrections at x = 0, where the exponential suppression is minimized, and neglecting the temperature dependence of V − and V P . Within this approximation we obtain where we set V + to be exactly zero at zero temperature so that its (strictly negative) value is purely controlled by the thermal corrections at the origin. The value of the potential at the true vacuum is −∆V , and independent of temperature in this approximation. With these approximations, the bounce action becomes simply and we are now ready to describe the shape of S 3 /T as a function of T . First we define the critical temperature T c , where the thermal corrections at the origin balance the zerotemperature potential difference between the two minima: where W(x) is the Lambert function, defined as the solution to the equation At large x the function W(x) behaves approximately like 3/4 log(1 + x), and this simple approximation can be used for all practical purposes here (see Appendix A for a short summary of the properties of the Lambert function). Using this, the low-T expansion will apply in regions of parameter space where where we have normalized the number of degrees of freedom in the thermal loops to the typical order of magnitude we will find in the explicit examples of Sec. 4. The low-T approximation is then valid whenever Eq. (3.13) is satisfied, making it a generic feature of the pseudomodulus potential where the vacuum energy is protected from quantum corrections induced by heavy SUSY states. From the definition of T c in Eq. (3.12), we can immediately see that the triangular approximation in Eq. (3.11) breaks down in this regime and should be extended (see Appendix B). However, the nucleation temperature in our setup is generically very far from T c , so that Eq. (3.11) is always a good approximation at the temperatures relevant for the FOPT. As the temperature decreases below T < T c , S 3 /T decreases as long as |V 0 T | > V P , since |V 0 T | decreases exponentially with the temperature. When the temperature approaches T min defined in Eq. (2.30), then |V 0 T | V P and S 3 /T attains a minimum value. As the temperature decreases further below T min , S 3 /T grows as 1/T .
Plugging the simplified bounce action in Eq. (3.11) into the T min definition in Eq. (2.30), we can easily obtain an analytic expression for T min which reads where to obtain the first expression we assumed T min 0.48m * and the second inequality follows from approximating the Lambert function W (x) 3/4 log(x + 1), assuming N ∼ O(10) and using the Eq. (3.4) for the scaling of V P with λ eff ∼ O(1) and ∆V ∼ F 2 . Higher values of λ eff or a suppressed value of ∆V will lead to a reduction of the hierarchy between T min and T c . The latter cases are less interesting from the point of view of the expected GW signal.
We are now ready to verify that there exists a nucleation temperature T n where S 3 /T satisfies the nucleation condition Eq. (2.14). As discussed in Eq. (2.30), the nucleation temperature is always within the interval (T min , T c ). Scenarios where T n is closer to T min have a larger α (see Eq. (2.15)) and a smaller β H (see Eq. (2.16)), favorable for generating an observable GW signal. Understanding the scaling of T n with respect to T min and T c thus provides valuable information about the strength of the FOPT.
Even approximating the nucleation condition in Eq. (2.14) with a constant C, solving the equation analytically with respect to T using S 3 /T given by Eq. (3.11) is not possible. We may, however, expand in V P /|V 0 T | 1 and solve for T n order by order in this expansion. This is always a good approximation as long as T n does not approach T min too closely. At first order, writing T n = T 0 n (1 + δT 1 n ) we find and we have again assumed T 0 n < 0.48m * . Given that the argument of the Lambert function is much larger than one, T 0 n depends only logarithmically on the parameters N, C, f a , ∆V , and can be taken proportional to m * for simplicity.
The leading scaling of T n with respect to the parameters shaping the potential is captured by the leading corrections proportional to V P in (3.15). Indeed, we observe that by increasing V P (i.e. the height of the barrier), or by increasing f a , the nucleation temperature decreases, approaching the region of parameter space where nucleation does not occur. The border between the nucleation and the non-nucleation areas is the portion of parameter space which is optimal for gravitational waves, since it is where β H is minimal. This behavior is in good agreement with the numerical results of Sec. 4, and one can verify that Eq. (3.15) reproduces the behavior of the full numerical result when properly matched to the models in Sec. 4 up to an overall scaling of the bounce action.

α, β H and fine-tuning
Now we can use our prediction for T n to compute the parameters characterizing the FOPT: • Within our analytical approximation, the temperature corrections only affect the potential at the origin of field space and are exponentially suppressed for T < m * .
Therefore, we approximate α as where the scaling of T n can obtained by using (3.15). Within this approximation, ∆V is temperature-independent and the largest values of α correspond to T n closer to T min .
• The inverse time scale of the phase transition can be computed explicitly from (3.11), giving One can easily verify that if T n = T min , then β H 1 within the small V P expansion. Moreover, the exponential dependence on T n makes β H very sensitive to the underlying parameters.
We now use the approximate β H formula in Eq. (3.18) to estimate the β H -tuning defined in Eq. (2.18). We compute first the tuning with respect to V P , which is encoded in Eq. (3.18) through the dependence of T n on V P . At leading order in V P /m 4 * 1 we obtain where in the last step we used the fact that since T 0 n < 10 21 m * and β H C in the interesting region of parameter space. The tuning associated with the barrier height is the dominant one, given that the tuning with respect to vacuum distance f a is suppressed by an extra factor of T 0 n /m * . As in Eq. (2.17), we see that the natural value of β H is β H C(T n ) 100 for the scales of interest in this study. Smaller values of β H can be obtained at the price of fine-tuning the barrier height at the percent level. This might imply an even larger tuning with respect to the fundamental parameters of a given model, as we will show in a concrete example in Section 4.2.

A toy example: fine-tuning vs. single SUSY-breaking scale
We now present a simple toy model which captures most of the features of the pseudomodulus potential in the explicit SUSY-breaking hidden sectors we will encounter in Sec. 4. We take the zero-temperature potential to be which reproduces the shape of the potential sketched in Fig. 3. The first term captures tree-level effects, while the second term captures one-loop quantum corrections. The x potential is flat at tree-level up to R-symmetry breaking operators parametrized by / R . 6 SUSY-breaking corrections induced by heavy fields lift the x potential around the origin, giving a mass to the pseudomodulus, but ultimately become subdominant for x √ F where SUSY is restored in the direction associated to the F -term. This large-field behavior is a unique characteristic of SUSY models.
As long as the explicit R-symmetry breaking is parametrically small, the position of the true vacuum and the zero temperature difference energy between the true vacuum and false vacuum are where we have introduced the parameter κ D to allow the scale controlling the difference in vacuum energy to vary relative to the scale controlling the loop corrections along the pseudomodulus potential. We will exhibit a concrete realization of such a model in Sec. 4.3.
Requiring the potential to be flat as in Eq.
Following the triangular barrier prescription, we need to find the position of the barrier and the value of the potential at the barrier; for the toy model these take the form From the last equation we see that the loop suppression of the zero-temperature barrier V P , as assumed in Eq (3.4), is here an automatic consequence of the fact that the pseudomodulus direction is lifted by quantum corrections. For a single-scale model (i.e. κ D = 1) the position of the peak x P is fixed in terms of the one of the true vacuum f a , while for a two-scale model, κ D 1 can enhance the hierarchy between x P and f a . We are now ready to use the triangular barrier approximation in Eq. (3.6) to compute the bounce action and the features of the FOPT between the origin and the true vacuum. is justified. If the general features of the bounce action characterize the FOPT in the low-T expansion discussed above, this simple toy model allows us to say something more precise about the scaling of the energy released during the FOPT. From Eq. (3.17) we have where we normalized the number of relativistic degrees of freedom at T n to be close to the MSSM value and we substituted T n ∼ T 0 n ∼ 0.5m * , which is the natural value of the nucleation temperature unless either V P or f a are tuned to suppress it (see Eq. (3.15)). In a 6 As shown in Sec. 4.2, the potential controlled by / R can be obtained from a marginal operator breaking R-symmetry in the superpotential. Similarly, one could study explicit R-breaking operators of arbitrary dimension in the superpotential W / R = / R X n nΛ n−3 which correspond to tree level potentials of the form . These types of operators would naturally be generated by UV dynamics as in Ref. [87]. single-scale model where κ D = 1, tuning T n √ F is the only way to enhance the strength of the FOPT. The same tuning will allow β H to be small. Conversely, in a two-scale model of SUSY-breaking, having κ D 1 can compensate the suppression in Eq. (3.25) without any tuning. We will show an explicit example of this class of hidden sectors in Sec. 4.3. These are clearly the best candidates to be probed by future GW interferometers.

Explicit Models
In this section we provide two working examples of the general idea described in the previous sections. Both models are straightforward deformations of the minimal O'Raifeartaigh model, which is the simplest theory of chiral superfields that breaks SUSY spontaneously [88].
The O'Raifeartaigh model involves three chiral superfields, namely the SUSY-breaking field X containing the pseudomodulus and two messenger fields Φ 1,2 . The dynamics are determined by three parameters: the SUSY-breaking scale √ F , the SUSY-preserving mass m of the messengers, and the coupling λ between the three fields. To set the stage for our analysis, we begin in Sec. 4.1 by determining the phase diagram of the minimal O'Raifeartaigh model which can be described as a function of the dimensionless parameter The model exhibits a rich phase structure as a function of temperature and the underlying parameters; for y F ∼ 1 the origin of the pseudomodulus is the global minimum at all T and no interesting phase transitions occur, while for y F 1 a second minimum develops away from the origin that may become the global minimum at intermediate temperatures, leading to a variety of phase transitions. Unfortunately, as we will see, none of these phase transitions are sufficiently strongly first-order to generate an observable GW signal. However, this minimal O'Raifeartaigh model serves as the foundation for SUSY-breaking hidden sectors that do generate observable GW signals.
In Sec. 4.2, we present the simplest SUSY-breaking hidden sector featuring a strong FOPT like the ones describe in Sec. 3. This hidden sector involves a marginal deformation in the superpotential of the minimal O'Raifeartaigh model, breaking the R-symmetry explicitly and obtaining a pseudomodulus potential very similar to the one described in the toy model in Sec. 3.4. We show that in such a simple single-scale model, α will be generically suppressed as predicted in Eq. (3.25), and discuss quantitatively the fine-tuning of β H defined in Eq. (2.18). Phenomenologically, this model is unsatisfactory since the global minimum restores SUSY, although this may be remedied by the introduction of external SUSY-breaking effects.
In Sec. 4.3, we show how both the shortcomings of the simple model of Sec. 4.2 are resolved in hidden sectors with two SUSY-breaking scales, in keeping with our expectations from Sec. 3.4. We make this concrete by gauging a U (1) flavor symmetry of the messengers in the minimal O'Raifeartaigh model, which admits an additional source of SUSY breaking via the Fayet-Iliopoulos term. This additional "D-term" supersymmetry breaking provides a second SUSY-breaking scale, which both ensures that supersymmetry is broken everywhere on the pseudomoduli space and increases α, leading to observable GW signals.

Warm up: The O'Raifeartaigh model at finite temperature
In the minimal O'Raifeartaigh model, the pseudo-modulus is stabilized at the origin by quantum corrections. Since the R-symmetry is unbroken in the global minimum at T = 0, one would expect that including finite temperature corrections will not induce any phase transitions. Instead, the dynamics of the O'Raifeartaigh model at finite temperature presents rich features that we discuss here in detail (see Refs. [89,90] for earlier works on related issues).
Having in mind applications to gauge mediated SUSY breaking, we consider the vectorlike version of the minimal O'Raifeartaigh model, which is described by the superpotential encoding the interactions of the SUSY-breaking chiral superfield X and two vector-like sets of messenger superfields Φ i ,Φ i (i = 1, 2). The first term is a tadpole ensuring that supersymmetry is broken at the scale √ F , while the second term encodes interactions among the fields with strength λ. We take the masses of the two pairs of messengers to be equal for simplicity. The superpotential above enjoys an unbroken R-symmetry under which X carries R[X] = +2, as well as a U (1) D flavor symmetry under which the messengers Φ andΦ have opposite charges (see Fig. 4 right for a summary table with the full charge assignment).
The potential for the scalar components of the chiral superfields is where X = x √ 2 denotes the scalar component of the pseudomodulus in the notation of Eq. (3.1). For λF ≤ m 2 , the tree level vacuum of the theory is at φ i =φ i = 0 with x undetermined, and SUSY is broken at a scale √ F . Radiative corrections from loops of the messenger fields φ andφ generate a potential for x that stabilizes it at the origin, and thus the global vacuum at zero temperature lies at φ i =φ i = 0 and x = 0. Note that the one-loop corrections have the shape described in Sec. 3, being polynomial close to the origin of the pseudomodulus potential and logarithmic for large field values. Expanding for y F ≡ λF m 2 ∼ 1 we obtain where we have fixed the renormalization scale to the messenger mass m. The thermal corrections to the x potential can be added with standard formulas that we review in the Appendix A. The shape of the thermal corrections is set by the x dependence of the mass eigenvalues for the scalar and fermionic components of the messengers. From (4.2) we can distinguish two classes of mass-squared eigenvalues: i) the ones growing quadratically with x, and ii) the ones decreasing as 1/x 2 and asymptotically going to zero in the large-x region. Specifically, the fermionic eigenvalues scale as and the bosonic eigenvalues are split in pairs around the fermionic ones, e.g. at the origin the bosonic eigenvalues are {m 2 , m 2 , m 2 + λF, m 2 − λF }. The behavior of the full spectrum as a function of x is shown in Figure 4 (right). We also observe that at large x, the spectrum asymptotes to a supersymmetric one. For low temperatures (i.e. T < m), the induced thermal corrections are a decreasing function of x, since they are mainly controlled by the lightest eigenstates. These corrections are mildly Boltzmann suppressed at large x and modify the pseudo-modulus potential as soon as T 4 ∼ λ 2 F 2 16π 2 . For larger temperatures, the contribution from the other mass eigenstates and in particular from the ones growing with x become relevant, and the thermal potential is a growing function of x. Hence at temperatures T ∼ m we expect the global minimum to be at the origin of the field space. However, for intermediate temperatures the thermal corrections can make the origin of the field space unstable, leading to a very rich evolution of the potential with temperature.  The thermal corrections compete with the loop corrections in the large x region (see Eq. (4.5)), eventually leading to a minimum of the potential at where x is obtained using the high-T expansion for the thermal potential up to T 2 , assuming 2 bosons and 2 fermions with masses-squared 2m 4 λ 2 x 2 , and T is an estimate of the temperature where the new minimum can be the global one. The latter is estimated by requiring the temperature corrections at x to be comparable to the height of the one loop potential. If T is close to m, then the neglected contributions from the states whose masses grow with x 2 lifts again the minimum at x , which will then never be the global minimum at any temperature. In conclusion, we expect that depending on the hierarchy between λF and m 2 , the minimum at x could become the global minimum in a certain temperature range around T . This complicated phase diagram is well summarized in Fig. 5, where we have fixed the ratio F m 2 to a representative value and explore the dynamics of the model as a function of the temperature and coupling λ.
For large λ, corresponding to y F ∼ 1, the minimum at X is never the global minimum of the scalar potential (green region in the plot). For small λ, i.e. y F 1 two phase transitions occur while lowering the temperature. Specifically, at very high temperature the global minimum is at the origin, as explained above. At intermediate temperatures the global minimum is at X , and finally at zero temperature the global vacuum is again at the origin. The corresponding two phase transitions can be first or second order. We have explored the parameter space of the model for different values of y F and λ, and found that these phase transitions are never strongly first order (i.e. small β H and large α) in the regime of perturbative λ.
Although the minimal O'Raifeartaigh model is itself not a good candidate for a strong FOPT, it nonetheless provides the foundation for simple variations that are. We explore these variations in the following subsections, restricting our attention to the region of parameter space in which the minimal O'Raifeartaigh model exhibits a simple thermal history corresponding to the green region in Figure 5. Deformations of the minimal O'Raifeartaigh model will endow this region with phase transitions as a function of temperature, while avoiding the complications of new minima arising from the interplay of thermal and loop corrections shown in the red and blue regions.

O'Raifeartaigh model with explicit R-symmetry breaking
Now we turn to a simple, concrete realization of a SUSY-breaking hidden sector whose pseudomodulus potential exhibits the properties exlpored in Sec. 3. This model simply amounts to deforming the minimal O'Raifeartaigh model studied in the previous section with the following marginal, R-symmetry-breaking term in the superpotential: The complete tree-level scalar potential of the model is and assuming y F ≤ 1, the global minimum sits at x true = 2F and φ i =φ i = 0.
In contrast to the minimal O'Raifeartaigh model, the R-symmetry-breaking deformation destabilizes the origin at tree level and restores supersymmetry in the true vacuum. The radiative corrections are identical to the ones in the O'Raifeartaigh at zeroth order in , and they tend to stabilize the pseudo-modulus at x = 0, competing with the tree-level contributions induced by the deformation. Close to the origin, the effective potential for the pseudomodulus obtained by integrating out the φ i andφ i fields reads (up to quartic order) where again we have approximated the loop corrections in the leading order in y F = λF m 2 ∼ 1. For < λ 3 16π 2 (log 4 − 1), the radiative corrections are sufficient to create a metastable vacuum at the origin of x. In this regime, the contribution to the quartic is always negligible. Along the pseudomodulus direction, there is now a true vacuum created by the R-symmetry-breaking deformation and a false vacuum created by radiative corrections. The height and location of the barrier between these two vacua may be approximated as This approximation is valid if x P m λ , that is if there is a cancellation between the two terms in m 2 eff such that m 2 eff λ 3 96π 2 F . The global minimum far away from the origin is not modified by the quantum corrections since SUSY is effectively restored there (we will come back to this point in Sec. 4.2.2) and stays at x true = 2F , so that ∆V = F 2 .
In summary, this hidden sector provides a concrete realization of the toy model discussed in Sec. 3.4. A direct consequence of having a single SUSY-breaking scale √ F is that the potential difference ∆V and the quantum corrections determining the barrier are both controlled by the same scale. This corresponds to κ D = 1 in the toy model of Sec. 3.4, and typically leads to suppressed α as we will show below. The formulae above allow straightforward matching of the model parameters onto the variables entering in the triangular barrier bounce action of Sec. 3.2. In the following, we will compare our analytical expectations with the full numerical analysis of the FOPT from the origin to the x true vacuum.

First order phase transition dynamics
We now study the model at finite temperature with an eye towards the dynamics of the phase transition associated to R-symmetry breaking. The thermal corrections to the X potential are equivalent to the ones that we studied in the simplest O'Raifeartaigh model, up to small corrections proportional to . The main difference is that the thermal effects are added on top of a zero-temperature potential described in the previous section, where the global minimum is far away from the origin. If we restrict to the parameter space where λF/m 2 ∼ 1 (the green region of Fig. 5), the role of the thermal corrections is to stabilize the origin at high temperature. Lowering the temperature, the thermal history is very similar to the one described in section 3: the negative thermal contributions at the origin decrease in absolute value until we reach T c , where the minimum at x = 0 is degenerate with the minimum at x true . An analytic estimate of this temperature can be obtained following Eq. (3.12). By further lowering the temperature, the thermal corrections become more and more negligible and one recovers the zero-temperature potential with a local minimum at the origin separated from the true vacuum by a loop-induced barrier.
Bounce action and nucleation temperature The next step in determining the phase transition dynamics is to compute the bounce action and the nucleation temperature. In the left panel of Fig. 6, we show the numerical result for the nucleation temperature T n as a function of the two dimensionless couplings of the model, having fixed y F ≡ λF/m 2 = 3/4 Numerically we see that T n ∼ √ F as assumed in the general discussion around Eq. 2.31. As can be seen from the explicit formula in Eq. 3.15, this feature is a consequence of the fact that F cannot be arbitrarily decoupled from m 2 if we require y F ∼ 1 and perturbativity of λ.
Keeping fixed, we see that the nucleation temperature T n decreases when λ increases, up until reaching the no-nucleation zone. Indeed, increasing λ makes the barrier between the two vacua higher, and the bounce action larger, decreasing the likelihood that the phase transition completes. When the bounce action increases, the nucleation temperature lowers, approaching T min at the border of the no-nucleation region. Conversely, decreasing λ at fixed shifts the nucleation temperature towards T c . In Fig. 6 we plot the quantity which indicates whether T n is closer to T c or T min . As we will see this quantity is strongly correlated with the strength of the signal. The main effect of on the pseudomodulus potential is to set the distance in field space between the origin and the true vacuum. The barrier is also -dependent, but away from the region where the effective mass in Eq. (4.10) changes sign, the effect of varying is negligible. Decreasing makes x true = f a larger. From Eq. (3.11) we see that the bounce action grows, making it more difficult for the phase transition to occur. This explains why lowering at fixed λ causes the nucleation temperature to decrease until the no-nucleation region is reached.
The interesting area of the parameter space is the sliver between the no-nucleation region and the region where there is not a barrier at T = 0. Within this sliver, λ 3 16π 2 ∼ and the effective mass of the pseudomodulus in Eq. (4.10) is small and positive (in units of √ F ∼ ∆V 1/4 ). The whole region shrinks for small because the true vacuum is pushed to large field values, and there is no nucleation unless the effective mass at the origin is tuned to be small.
The scaling of T n / √ F with the Lagrangian parameters can be captured by the analytic approximations presented in Sec. 3. We match the generic parameterization of Sec. 3 using the expressions in (4.10) and (4.11), giving where we have approximated the radiative corrections for y F ∼ 1. This expression qualitatively reproduces the left panel of Fig. 6, up to an overall normalization of the bounce action (corresponding to a shift in C).
For large , the rightmost term in parentheses in (4.13) is always O(1), and hence the variation of T n is largely controlled by the prefactor. The ∼ √ λ scaling of the numerator is balanced by the log λ scaling in the denominator, and the resulting prefactor of T n is essentially flat in λ and only decreases with decreasing . In the small region, the rightmost term in parentheses in (4.13) becomes smaller than 1 and controls the shapes of the T n contours, in agreement with the left panel of Fig. 6. α, β H and fine-tuning The microscopic properties of the FOPT dynamics are encoded in the two parameters α and β H , which correspond to the energy release and the duration of the phase transition. In the central and rightmost panels of Fig. 6, we show the behavior of β H and α in the (λ, ) plane. Here the scaling of α is essentially dictated by the scaling of T n / √ F , since α ∼ 30 g * π 2 ∆V T 4 n and ∆V ∼ F 2 . We see that in the parameter space explored here, we cannot reach large values of α except in the thin sliver towards small and λ where T n ∼ T min and T min is minimized with respect to √ F . This unfortunate feature is a generic prediction of a single-scale SUSY-breaking hidden sector, as discussed in Sec. 3.4.
As shown in Fig. 6, β H is small in the regions of the parameter space at the border of the no-nucleation zone where T n ∼ T min . Getting closer and closer to this boundary, one can achieve β H 100 at the price of a large tuning of the model parameters as discussed in Sec. 2.2. The fine-tuning is dominated by the tuning of the barrier V P between the two minima. Substituting the dependence of V P on the Lagrangian parameters, we can estimate the the tuning of β H with respect to λ as (4.14) The same fine-tuning can be computed numerically using the prescription of Eq. (3.19).
We show the results in Fig. 7, where we see that β H 100 corresponds to ∆ β H ∼ 10 3 , which is larger than the estimate derived above. The numerical results also confirm our expectation that the tuning associated with the λ parameter (setting the height of the barrier) dominates relative to the tuning associated with the parameter (setting the location of the true vacuum). Comparing the left and right panels in Fig. 7, it is apparent that the tuning grows in the region of small , where α is larger.

Phenomenological Challenges
As discussed in the previous section, the model presented here is not optimal for generating a sizable SGWB signal. Indeed, by comparing the resulting values of α and β H in Figure  6 to the values displayed in Figure 2, it is clear that the typical value of α is too small to lead to a detectable signal. Of course, this issue can be resolved by going in a tuned region of the parameter space where a very small and an appropriately fine-tuned λ give T ∼ T min and a suppressed T min compared to √ F . However, it is fair to say that, in general, a perturbative single-scale SUSY-breaking hidden sector cannot lead to strong SGWB signals. For this reason, we do not display the SGWB for this model, although it may easily be inferred from the α, β H , T n plots in Figure 6.
As shown in the simple toy model of Sec. 3.4, the suppression of α is a consequence of the fact that in a single-scale SUSY-breaking hidden sector one cannot significantly separate ∆V from T 4 n . Notice that this conclusion hinges on requiring y F ∼ 1, which is necessary to avoid the region shown in Fig. 5 where thermal corrections at high temperatures induce new minima. A more careful study of the dynamics of single scale models for y F 1 is left for future work.
In addition to the α suppression, the true vacuum in this model restores SUSY, so that the phase transition is not genuinely a "SUSY-breaking phase transition"; this sector, as presented, cannot be responsible for the SUSY breaking transmitted to the MSSM. However, this is not a fatal obstruction, as the vacuum energy in the true vacuum far from the origin can be easily lifted by coupling to another source of SUSY breaking, very much in the spirit of [91]. If we require the new source of SUSY-breaking to not significantly affect the dynamics of the phase transition, the parametrics of the model will not significantly deviate from those presented here and the resulting α will be still suppressed. Interestingly, our analysis seems to point towards SUSY-breaking hidden sectors with multiple dynamical field directions and scales. In the next section, we will exhibit the simplest model of this type, leaving a more thorough exploration of the different possibilities for future study.

O'Raifeartaigh model with gauge interactions
In the previous subsection we analyzed a simple model displaying a first order phase transition associated with the breaking of the R-symmetry. However, there were two aspects that were not completely satisfactory: i) SUSY breaking in the global minimum had to be added as a further deformation, and ii) the phase transition was generically not strong enough to generate a sizable signal. Both issues were related to the fact that there was only one SUSY breaking scale in the problem. In this subsection we resolve these issues in a hidden sector where the global minimum breaks both SUSY and R-symmetry spontaneously and the presence of two SUSY breaking scales leads to a strong FOPT from the origin to the true minimum. It is well-known that adding gauge interactions to SUSY breaking models with chiral superfields modifies the potential and typically leads to a new SUSY-and R-symmetrybreaking vacuum at large field values (see e.g. [14]). As a prototype of this class of models we consider the simplest realization, which consists of the vector-like O'Raifeartaigh model of the previous sections where the anomaly-free U (1) D flavor symmetry defined in the right panel of Fig. 4 is gauged. The model we consider has been studied at zero temperature in [92]. The qualitative features that we find here are generic to models where SUSY is broken through the interplay of F -and D-term effects.
The field content and superpotential are the same as those introduced in (4.2). The gauging of the U (1) D symmetry contributes new terms in the scalar potential from the D-term contribution. The F -and D-term contributions to the potential together give where g is the gauge coupling of the U (1) D symmetry and we have also included a UV Fayet-Iliopoulos (FI) term D/g. This FI term contributes a second source of SUSY breaking that will strengthen the GW signal. Note that the model contains, in addition to the O'Raifeartaigh degrees of freedom, a gauge boson and gaugino associated with the U (1) D vector multiplet. We focus on the regime where y F 1 and we do not discuss the origin of the FI term here. 7 The scalar potential at zero temperature As a first step, we analyze the zerotemperature vacuum structure and map it onto the parameterization of Section 3. Neglecting the gauge dynamics, the tree-level potential has a minimum at φ i =φ i = 0 where SUSY is broken everywhere along the F -flat pseudomoduli space parameterized by X. Including the gauge interactions, the minimization of the D-term part of potential favors configurations where theφ i fields acquire a VEV to compensate for the FI term D/g. This results in a tension between the minimization of the F -term and D-term contributions to the potential. While the F -term can never be set to zero, one can find a runaway direction in field space which leads, asymptotically, to the vanishing of the D-term.
First, we can solve for the F -terms of Φ 1 andΦ 2 by taking On this solution the scalar potential simplifies to (4.18) Note that φ 1 ,φ 2 have vanishing R-charge, so the only direction where the R-symmetry is spontaneously broken is along x. In order to visualize the shape of the scalar potential and the approach to the runaway, we show in Fig. 8 the tree level scalar potential as a function of x, as well as the values of φ 1 andφ 2 as a function of x. The scalar potential is flat around the origin and then turns to the runaway direction along which the D-term diminishes. The turning point along x is where the fields φ 1 and φ 2 acquire a non-vanishing VEV. The VEV of φ 1 is different from zero since the potential energy is most efficiently minimized if φ 1 partially cancels the first term in (4.18) as well as minimizing the D-term. The VEV of φ 1 is suppressed by a factor ∼ gF λD with respect to the VEV ofφ 2 . An analytic estimate of the scalar potential can then be captured by working at zeroth order in the VEV of φ 1 . In this approximation, and focusing on the parameter region where gD/m 2 < 1, the effective mass-squared for theφ 2 field is X-dependent, and turns negative at the transition point x trans where the fieldφ 2 develops a VEV. 7 The inclusion of a fundamental FI term is not strictly required to obtain a strong FOPT. A very similar potential for the pseudomodulus can be obtained by considering two different masses for the messengers and working in the regime where λF > m1m2. Models with multiple F -terms would also lead to similar conclusions.
The potential for x is flat for x ≤ x trans , while for x ≥ x trans it can be obtained by integrating outφ 2 , (4.20) Since we work in the small-g regime, radiative corrections from the gauge sector may be neglected, such that the 1-loop corrections are the same as the ones discussed in the previous sections (see Eq.s (4.4) and (4.5)). They have two effects, namely i) they create a local minimum at the origin, and ii) they generate a global minimum at large x values along the D-flat direction.
The barrier between the two vacua is approximately at x x trans where we can estimate the one-loop potential simply by the large field behaviour in (4.5), giving Combining the approximate tree level potential in (4.20) with the loop corrections in (4.5), we find that the true vacuum at large field values lies at where the difference in potential energy between the two minima is dominated by the D-term contribution. This completes the matching of the potential of this model to the general discussion of Section 3. Note that here the SUSY-breaking F -term controls the height of the barrier in Eq. (4.21), while the SUSY-breaking D-term sets the potential energy difference as in Eq. (4.22) . This implies that the phase transition can have sizable values of α, as we will see in the numerical analysis.

First order phase transition dynamics
We now turn to the finite-temperature corrections and compute the parameters associated with the phase transition. Note that the spectrum is similar to the O'Raifeartaigh model with the addition of the gauge boson and the gaugino of the U (1) D symmetry. These additional states are massless in the false vacuum and massive in the true vacuum, so they contribute to making the origin the global minimum at high temperatures. We numerically evaluate the one-loop and the thermal corrections to the scalar potential, and then compute the bounce action for tunneling from the false vacuum to the true vacuum. In Appendix B.1.2 we present the triangular barrier approximation for this model and compare it with the full numerics. Even though the bounce profile in field space involves three different fields, i.e. (x, φ 1 ,φ 2 ), in our numerical scan, we approximate the bounce as one-dimensional, neglecting the contribution from the φ 1 ,φ 2 directions. As detailed in Appendix B.2, we checked the single field approximation against the full 3d bounce action computed numerically with both FindBounce [85] and CosmoTransitions [86]. As a result, the single field approximation gives a good description of the bounce as long as Figure 8. Left: Tree level and one-loop scalar potential as a function of the pseudomodulus direction x minimizing the directions φ 1 andφ 2 . The dashed blue line shows the tree level potential which is flat around the origin and develops a runaway at x P x trans (see Eq, (4.19)). Quantum corrections generate a local minimum at the origin as shown by the black solid line in the small quadrant and a global minimum far away in field space indicated with a green dashed line. The difference in energy density is ∆V  . The GW signal weakens going from red to blue. In the gray shaded region at the bottom R 1d/3d > 0.5 and as described in Eq. (B.29) our 1d approximation is expected to break down. The grey region on the left is excluded by the perturbativity of λ below m. In the white region the nucleation condition in Eq. (2.13) cannot be satisfied. φ 1 ,φ 2 are smaller than X at the bounce release point (defined as the starting point of the tunneling set at r = 0, where the kinetic terms of all the fields are exactly zero). In order to estimate where we expect sizable deviations from the multidimensional contribution, we borrow some intuition from the triangular barrier approximation, where S 3 /T scales as ∼ X 3 , and define where X(0), φ 1 (0) andφ 2 (0) are the field distances from the origin computed at the release point r = 0. In Fig. 9 we show the region where R 1d/3d > 0.5 and we expect deviations of 50% or more from our one-dimensional estimate of the bounce action. As we can see, this region is not phenomenologically relevant since it is quite far from the interesting region for GW signals. In Fig. 9 we show the behavior of T n , α, and β H in the (λ, g) plane, having fixed √ F = 30 PeV and y F = 3/4 as in the previous model and set y D = 1/5. Keeping fixed the ratios in Eq. (4.24), the triangular barrier parameters scale as As a consequence of these scalings, using Eq. (3.11) it is straightforward to see that for fixed λ the boundary of the nucleation region is reached for large g, while for fixed g the boundary lies at small λ. The shape of the nucleation temperature T n can be captured by a simple analytic formula after rewriting Eq. (3.15) in terms of the theory parameters: This expression reproduces the contours in Fig. 9 (left) up to overall normalization. In the middle and right panels of Fig.9 we show the contours for β H and α, respectively. The main difference compared to the model in Sec. 4.2 is that even if T n ∼ √ F , it is possible to obtain sizable values of α because ∆V is here controlled by the D-term. Approaching the boundary of the nucleation zone without fine-tuning the theory parameters by more than O(1) we can reach β H ∼ 100 and α ∼ 0.3 − 0.4, which we use as a benchmark for our summary plot in Fig. 1.
The interplay of the two SUSY-breaking scales √ F and √ D is an essential ingredient for a strong FOPT. This is illustrated in Fig. 10, where we show the behavior of T n , α, and β H in the (λ, D/F ) plane, having again set √ F = 30 PeV and y F = 3/4, and now fixing g = 0.1. In this scaling the FOPT is essentially independent of λ, and one can see clearly that the separation of D from F is the crucial ingredient for a sufficiently strong phase transition. Notice that the required separation is O(1) and therefore not obviously in tension with theoretical bounds on large D-terms [93]. Strictly speaking, these bounds do not apply to our simple model, where a tree level Fayet-Iliopoulos term makes the Ferrara-Zumino multiplet not gauge invariant [94]. However they would have applied if we were to UV complete this model to a full-fledged model of dynamical SUSY-breaking or for instance if we were to explore the second branch of the model with two different messengers masses and λF > m 1 m 2 .

Gravitational Wave spectrum and phenomenology
Having shown that a strong FOPT can be achieved without fine-tuning in a SUSY-breaking hidden sector with at least two SUSY-breaking scales, we now turn to the gravitational wave signal itself, again using the simple model presented in the previous section as a benchmark. As discussed in Sec. 2.2, computing the SGWB signal requires understanding the macroscopic dynamics of the vacuum bubbles expanding in the plasma. This is essentially determined by the balance of the energy ∆V released in the FOPT and the pressure effects from the plasma (see Eq. (2.19)). If pressure effects stop the bubbles before they collide, most of the SGWB signal will be sourced by the energy released in the plasma.
In Sec. 2.2 we described a quite unique friction mechanism at work in our class of models. This mechanism is a direct consequence of two peculiar features of the pseudomodulus potential: i) the nucleation temperature T n ∼ √ F is set by exponentially suppressed temperature corrections to be smaller than the typical scale of the heavy states in the theory m, and ii) the true vacuum VEV is typically the larger scale in the problem and controls the mass variation ∆m 2 /m 2 ∼ λ 2 f 2 a /m 2 1 of the heavy states from the false to the true vacuum. These two properties together imply that when the vacuum bubbles accelerate enough, γT n > m true and the heavy states can cross the bubble wall. Their crossing switches on a new pressure effect which is generically larger than ∆V and immediately stops the bubble runaway. This last statement can be checked explicitly with the parametric dependence of the simple model described here. The heavy state pressure term in Eq. (2.21) scales as where we used the scaling of the true vacuum as a function of the theory parameters in Eq. (4.22) and approximated T n √ F for simplicity (this approximation is numerically correct up to an O(1) factor as shown by the dashed contours in Fig. 9). Inside the exponential, we should take the lightest heavy states in the plasma m false ∼ √ m 2 − λF which are of course less Boltzmann suppressed and dominate the friction. Comparing this quantity with the energy released in the phase transition ∆V = D 2 /2 we can get the range of the gauge coupling g such that this friction prevents the bubble runaway, Plugging in the typical numbers for our phase transition (F/D ∼ 1/5, y F ∼ 3/4 and m false / √ F √ λ 2.5) indicates that the vacuum bubbles are always stopped in the range of interest for the gauge coupling g for perturbative values of λ. The predicted boost factor at equilibrium in this case is where again f a is defined in this model by Eq. (4.22). As a final remark, we notice that the NLO friction induced by gauge degrees of freedom radiated through the wall never dominates over the one from heavy states in the interesting range of the gauge coupling g.
Given that the bubble runaway is always prevented, the dominant SGWB comes from sound waves in the plasma. The predicted energy fraction as a function of frequency at GW interferometers has been discussed in Eq. (2.25) and below. Putting everything together, in Fig. 11 we compare our model predictions with the PLI curves for future GW interferometers derived in Appendix C.1. This clearly demonstrates that SUSY-breaking hidden sectors with multiple SUSY-breaking scales can generate stochastic signals detectable at future GW interferometers. Moreover, it makes explicit the expected correlation between the SUSY-breaking scale and the peak frequency of the resulting SBGW. All that remains is to explore the full range of viable SUSY-breaking scales (and hence signal frequencies), as well as the correlation between signals at GW interferometers and other experiments. In the next section, we will bound the SUSY-breaking scale from above around ∼few tens of PeV by computing the gravitino cosmological abundance. By specifying a mediation mechanism, we will also use the explicit hidden sector presented here to show how the SUSY-breaking scale determines the spectrum of MSSM superpartners, thereby correlating signals at GW interferometers and future colliders.

Phenomenology
Having demonstrated that the first-order phase transition in a SUSY-breaking hidden sector can generate an observable GW signal, we now turn to complementary aspects of hidden sector phenomenology that shape the motivated parameter space and suggest additional experimental tests in the event of a signal at GW interferometers. We begin with universal features that are intrinsic to the hidden sector itself and independent of the mediation mechanism that connects the hidden sector to the MSSM. This includes key aspects of gravitino cosmology, where we will see that the requirement T r.h. = √ F implies an upper bound on √ F even if m 3/2 receives extra contributiosn from other SUSY-breaking sectors as in Eq. (2.2). We also explore the prospects for collider searches for the gravitino (independent of the MSSM spectrum), finding that future high energy lepton colliders could probe almost the entirety of the light gravitino window (i.e. m 3/2 < 16 eV) by directly producing gravitino pairs. We then relate the parameters of the hidden sector to the spectrum of the MSSM, which requires specifying details of the mediation mechanism. Here we consider the prototypical example of gauge mediation via vector-like messengers, where the parameter space for observable GW signals generates a superpartner spectrum within reach of future proton-proton colliders such as FCC-hh.  [24][25][26][27] and the future FCC-hh reach on gluinos [28] for perturbative messenger sectors with g M ∈ (0.01, 0.1) (see Eq. (2.10) for a definition of g M ). The dark blue dashed lines show the values of κ = F/F 0 , the ratio between the total SUSY-breaking scale F 0 and one controlling the soft masses (see Eq. (2.2)). As discussed in Eq. (2.11), we expect constraints on flavor changing neutral currents to exclude κ 10 −8 as indicated by the dark blue arrows. The dark magenta thick line indicates the BBN bound on the higgsino NLSP decaying to gravitino plus hadrons as obtained in [30].

Gravitino cosmology vs future colliders
The gravitino overabundance is a well known problem of LESB scenarios [37][38][39]. This problem is exacerbated in our setup, because having sizeable GW signals from the SUSYbreaking hidden sector requires the reheating temperature T r.h. to be at least as high as the SUSY-breaking scale, enhancing the gravitino production from scattering as detailed in Eq. (2.8). In light of this tension, here we delve into further detail about the two viable scenarios sketched Sec. 2.3. Since T r.h. ∼ √ F , which is much larger than the scale of the soft masses, the main player in determining the final gravitino abundance is the production from UV scattering computed in [42][43][44][45]. The final yield can be written as where the production through gluon-gluino scattering dominates over the other channels and we have substituted T r.h. √ F , which is the lowest reheating temperature compatible with our scenario. Following Eq. (2.2), we assume that the gravitino mass m 3/2 is set by an independent SUSY-breaking scale F 0 = F/κ, possibly higher than the one setting the soft spectrum (i.e. κ 1).
Ultralight gravitino window vs. pair production at future colliders If the gravitino mass and the soft spectrum are set by the same SUSY-breaking scale √ F , the yield scales as Y UV 3/2 ∼ M Pl / √ F . For sufficiently low SUSY-breaking scales, the yield becomes just the equilibrium one, Y UV 3/2 > Y eq , where Y eq = n eq 3/2 /s = 1.8 × 10 −3 . The gravitino is a thermal relic as long as √ F 1/3 , which corresponds to √ F 8.6 × 10 7 GeV for M 3 = 2 TeV. Moreover, since √ F m 3/2 , the gravitino is relativistic at freeze-out and its abundance today is constrained by measurements of the matter power spectrum at short scales [40,41]. The current bounds imply m 3/2 16 eV , F 260 TeV . (5. 2) The above requirement identifies the ultralight gravitino window. Although it is unquestionably challenging to decouple the soft spectrum from the LHC in this window (see [78,79] for attempts in this direction), it is interesting to ask whether future colliders can test this window in a model-independent fashion through direct pair production of the longitudinal component of the gravitino, the goldstino. This production rate depends directly on √ F even when the MSSM superpartners are decoupled, and so provides a direct experimental test of the SUSY-breaking sector.
The projected sensitivity to gravitino pair production at both hadron and lepton colliders is displayed in Figure 1. For the bound at future lepton colliders, we consider a high energy lepton collider operating at √ s = 30 TeV. Assuming minimal cuts on the photon kinematics (E γ > 50 GeV, |η γ | < 2.4), the signal cross section from Eq. (2.7) is Applying the same minimal cuts on the photon, the SM background (estimated with Mad-Graph5 [95,96]) is σ SM 2pb. In contrast to LEP, at high energy lepton colliders the SM background is dominated by W W fusion while the Drell-Yan process with an ISR photon is negligible. We can then derive a lower bound on the scale of SUSY breaking as displayed in Figure 1 given an assumed integrated luminosity, namely which is still an order of magnitude away from entirely closing the ultralight gravitino window. However, improved analyises and new cosmological data could strengthen the gravitino mass bound by an order of magnitude, potentially closing the ultralight gravitino window completely. For instance, Ref. [97] already claims a bound on the gravitino mass of m 3/2 < 4.7 eV; although the robustness of this bound is subject to interpretation, improved limits from Planck data are likely to be comparable. In order to estimate the reach of future hadron colliders, we perform a rescaling of the limits discussed in [23], based on the mono-photon search of ATLAS [98], which constrain √ F 850GeV with 20.3fb −1 at √ s = 8 TeV. 8 As at lepton colliders, the signal cross section for gravitino pair production in association with a photon σ(pp →GGγ) scales as s 3 /F 4 at hadron colliders [22]. To estimate the limit attainable at √ s = 100 TeV, we first compute the ratio of the signal cross sections at √ s = 100 TeV and √ s = 8 TeV, taking the partonic signal cross section to scale as σ sig ∼ŝ 3 /F 4 and assuming that the p T,γ ≥ 125 GeV cut at √ s = 8 TeV is increased to p T,γ ≥ 1 TeV at √ s = 100 TeV. We additionally compute the ratio of background cross sections, assuming the partonic background cross section scales as σ bkg ∼ 1/ŝ. Using the √ s = 8 TeV signal and background predictions in [98] and the above ratios, we find the expected limit at √ s = 100 TeV to be which is the one displayed in Figure 1. Even with this aggressive estimate, the reach of high energy hadron colliders is limited compared to the reach of high energy lepton colliders because the signal cross section at the former only grows withŝ 3 , while at the latter it grows as s 3 .
Gravitino Dark Matter window If the SUSY-breaking scale √ F 0 setting the gravitino mass exceeds the scale √ F of the hidden sector, we can treat m 3/2 as a free parameter and access an interesting region where the gravitino is never in thermal equilibrium with the SM. This could arise naturally from additional sequestered sectors that break supersymmetry at higher scales. As shown in Fig. 12, we also require this new source of SUSY-breaking to not spoil the defining phenomenological features of LESB, namely i) the gravitino is still the LSP, and ii) the soft masses are dominated by the flavor-diagonal contribution from gauge mediation.
Requiring the gravitino avoid thermalization, Y UV 3/2 < Y eq , we obtain an upper bound on √ F at fixed gravitino mass which can be cast as an upper bound on the ratio between the two SUSY-breaking scales, κ = F/F 0 : Here we have used the expression for the gravtino mass in Eq. (2.2) and the one for the gluinos in Eq. (2.10), where the parameter g M encodes the model-dependence of the latter. If the gravitino is never in thermal equilibrium, we can assume (as usual in freeze-in scenarios) that the gravitino sector is not directly reheated after inflation and the gravitino abundance is frozen-in through scattering of SM states with their superpartners. Setting the gravitino abundance to explain the DM abundance today, we can predict the gluino mass in the (m 3/2 , √ F ) plane, which corresponds to the red lines of Fig. 12 where the gravitino accounts for the total DM abundance today at fixed gluino mass. The current LHC bounds on the gluino mass set a boundary of our parameter space, which is shown in Fig. 12. 9 The second scaling in Eq. (5.7) shows the value of κ required to achieve a given gluino mass. As shown in Fig. 12, the parameter space of interest has κ between (10 −8 , 10 −4 ), where smaller values of κ would not open up more parameter space and in any event would be in tension with FCNC constraints as discussed in Eq. (2.11).
The bound at larger gravitino masses (the gray band on the r.h.s. on Fig. 12) is given by the requirement that the gravitino be the LSP. A stronger bound is derived from BBN constraints on the freeze-out abundance of the NLSP decaying into gravitinos. We have computed the NLSP freeze-out abundance assuming the NLSP is a pure higgsino NLSP and applied the BBN bound of Ref. [30] given the NLSP lifetime in Eq. (2.6). The triangular-shaped region where the gravitino could be DM can be probed by both GW interferometers and future colliders, as shown in Fig. 12. This highlights the potential for future colliders to determine whether a SUSY-breaking phase transition is the source of a SGWB signal observed at GW interferometers.

A complete model of gauge mediation
Finally, we can correlate the GW signals of the SUSY-breaking hidden sector with the superpartner spectrum of the MSSM by specifying a mediation mechanism. In order to embed the model of Sec. 4.3 into a successful model of gauge mediation, we work in terms of a simple generalization of the gauged O'Raifeartaigh model which allows a natural embedding of both the gauged U (1) D symmetry and the SM gauge group into the flavor symmetry of the messengers. This requires M copies of the vector-like messengers Φ andΦ coupled to the singlet X, so that the superpotential is identical to the one of the O'Raifeartaigh model in Eq. (4.2), but where now the fields are intended as vectors with M components. The minimal setup requires M = 6 so that the superpotential enjoys an SU (6) symmetry, where a U (1) subgroup of the SU (6) is the gauged U (1) D with a non vanishing Fayet-Ilipoulos term, while the SM gauge group lies inside the remaining global SU (5) such that the messengers can be taken to transform in the 5 +5 representation of SU (5) as in standard gauge mediation scenarios [35].
The mass matrix of the messenger fields is where f a is the VEV of the pseudomodulus given in Eq. (4.22). Integrating out the messengers, one can compute the soft masses for the MSSM following the general formulas in [99]. The scalar masses follow the standard gauge mediation scaling discussed in Eq. (2.9), while it is worth explicitly writing the parametric dependence of the gluino soft mass in the notation of Eq. (2.9): Here we have expanded in λf a m F and identified the gaugino screening factor s M in this model. Since we typically have y F ∼ 1 in our scenarios (in order to remain in the green region of Fig. 5) the gaugino screening factor does not provide significant suppression, but interestingly it is generic for models like ours where the messengers mass matrix is never singular along the pseudomodulus direction [15,51]. Abandoning this requirement, one could avoid gaugino screening at the price of opening up messenger field directions where the SM gauge group is spontaneously broken in the UV [100].
Substituting the value of f a in Eq. (4.22) and taking as benchmark values a typical point with α ∼ 0.3 and β H ∼ 100 from Fig. 9, the gaugino pole mass is This shows that the band between the the present exclusion at the LHC and the future reach of FCC-hh can be populated with simple, concrete models featuring strong SGWB signals within the reach of future high-frequency interferometers such as A-LIGO, ET and CE.

Conclusions
We began by asking if future gravity wave detectors could provide a new window into supersymmetry by probing SUSY-breaking hidden sectors in a region not yet excluded by LHC searches. The answer to this question is well summarized in Fig. 1, which shows the complementarity of future gravitational wave interferometers and colliders in probing scenarios of low-energy supersymmetry breaking (LESB). Fortuitously, the cosmological history of the gravitino -a key degree of freedom in LESB scenarios -bounds the SUSYbreaking scale from above, so that the viable parameter space lies within reach of both high-frequency GW interferometers and high-energy colliders. The underlying assumption in Fig. 1 is that the SUSY-breaking hidden sector actually undergoes a strong first-order phase transition. The remainder of the paper has been devoted to demonstrating, on general grounds, the circumstances under which strong FOPTs can be produced in SUSY-breaking hidden sectors. We have focused on phase transitions along the pseudomodulus direction, which as a universal feature of spontaneous SUSY-breaking is guaranteed to exist in a vast class of SUSY-breaking hidden sectors. Remarkably, the generic features of the pseudomodulus potential gave rise to a parametrically new way of realizing strong first-order phase transitions in field theory.
The novelty of the pseudomodulus FOPT is a consequence of the flatness of the treelevel potential accompanied by the presence of a mass gap for the heavy states, which makes the theory calculable everywhere in field space. Since the mass gap is supersymmetric, it does not destroy the flatness of the potential at large field values. The two resulting features of this setup are that i) the nucleation temperature is well below the scale of the heavy states, so that the low-T expansion applies, and ii) the pressure from Boltzmann-suppressed heavy states in the plasma is responsible for stopping the vacuum bubble runaway. The dominant GW signal then comes from the energy released in the plasma during the phase transition.
The strength of the GW signal depends on finer details of the hidden sector dynamics. However, we found that multiple SUSY-breaking scales in the hidden sector are a necessary condition for generating strong GW signals without fine tuning of the theory parameters. This result is quite general and can be obtained analytically without reference to a specific model. For the sake of concreteness, we presented an explicit model for a hidden sector generating a strong GW signal, where SUSY is broken by both an F -term and a Dterm. Detailed predictions for the GW signal and superpartner spectrum in this model substantiate the general phenomenological observations of Fig. 1. scheme is given by where F = 1 (0) for fermions (bosons), the number of degrees of freedom associated with the particle i is g i = 1/2/3 for real scalars, fermions and vectors, respectively, and c i = 3 2 ( 5 2 ) for scalars/fermions (vectors). The thermal one-loop potential is given by where the thermal functions for both species are with z i ≡ m i /T . These functions can only be fully evaluated numerically, but admit analytical approximations for large and small |z 2 |. In the high-temperature limit, z 2 1 and the thermal functions are where a b = π 2 exp (3/2 − 2γ E ) and a f = 16π 2 exp (3/2 − 2γ E ). The low temperature limit (i.e. z 2 1) can be approximated in terms of modified Bessel functions of the second kind where m is high enough such that the series converge. For T low enough, we can take only the first term in the series and further expand the modified Bessel function to the leading order K ν (z) z→∞ π 2z e −z to get Within this approximation we can obtain a simple expression for V th (x, T ) which is valid at the leading order in the low-T expansion: where it is important to notice that bosons and fermions contribute with the same (negative) sign to the effective potential. The approximation above is used in Section 3 to derive an analytical scaling of the dynamics of FOPTs.
Lambert function As a consequence of the low-T expansion the equations we will be dealing with have the typical form where W(z) is the Lambert function, which is defined such that W(z)e W(z) = z. Without entering into the details of the interesting properties of this function, we restrict our interest to finding a good approximation for it using simpler functions. First, we will consider W(z) for strictly positive arguments. Second, we note that W → 0 for z → 0 and that W(z) ∼ log z for z → ∞. By inspection one finds that a good approximation of W(z) is given simply by The relative difference between the Lambert function and the approximation with the logarithm in (A.10) is at most ∼ 1/4 (for z → 0 and z → ∞) and smaller (in absolute value) in intermediate regions. For practical purposes in analytic estimations of relevant quantities, we will hence often consider the approximation in (A.10).

B Bounce action computation schemes
The transition of a quantum system from a meta-stable vacuum state to the true vacuum can be driven either by quantum tunneling or by thermal fluctuations. In the FOPTs describe in this paper the latter are always dominant. The probability of thermal tunnelling is described semi-classically by Eq. 2.12 and it is exponentially dependent on the classical O(3)-symmetric bounce solution [64,65]. In this appendix we review both the analytical and the numerical approaches we used to study the behavior of the O(3)-symmetric bounce in our FOPTs. First, in Sec. B.1, we describe in detail the triangular barrier approximation introduced in Ref. [80] and its generalization to the O(3)-symmetric case [81]. We compare this approximation with the behavior of the full bounce action computed numerically with the FindBounce package [85,101] and the CosmoTransitions code [86]. In Sec. B.1.1 we discuss an optimization of the triangular barrier approximation which leads to excellent agreement with the full numerical computation. In Sec. B.1.2 we consider the O'Raifeartaigh model with gauge interactions of Sec. 4.3 as our main case study.
Second, in Sec. B.2 we discuss the single field approximation of the bounce action in the model of Sec. 4.3. We study numerically the behavior of the bounce action in the full three-dimensional field space and compare it with the single field approximation, deriving where we expect the latter to deviate sensibly from the full solution.

B.1 Triangular barrier approximation of the bounce action
The d-dimensional Euclidean action for n scalar fields φ i is where Ω d = 2π d/2 /Γ(d/2) and the equations of motion for φ i arë The boundary conditions defining the bounce solution arė where the conditions at r → ∞ ensures that the solution starts with zero kinetic energy from the false vacuum and stops at r = 0 with zero kinetic energy. Here, we are interested in solving this equation for a single field (n = 1) and in d = 3. The idea behind the triangular barrier approximation is to approximate the potential as a piecewise linear function anchored at three points: the false vacuum, the top of the barrier and the true vacuum. In the following, we review in some details the derivation of the 3d bounce action in the triangular barrier approximation, and we refer to the original paper [80] for more details and to [102] for the derivation in arbitrary dimensions.
Following the notation of [80], we define the false vacuum to be at the field position φ + with potential V + ; the peak of the triangular barrier to be at φ P with potential V P ; and the true vacuum to be at φ − with potential V − . It is then convenient to define the magnitudes of the gradients of the potential by so that V (φ) = ±λ ± on either side of the barrier, precisely In order to solve the equation of motions we have to specify the boundary conditions. At a large radius R + the field attains the false vacuum, so we have Then, we have to specify the boundary conditions at the start of the tunneling. There are two possibilities: 1. The field immediately start to roll at r = 0 and hence we impose where the initial field value φ 0 is the undetermined release point. This is the valid regime if one finds that φ 0 ≤ φ − .
2. Otherwise, the field sits in the true vacuum for r < R 0 and then starts rolling. In this second case the boundary conditions are We begin with the analysis of the first case. Imposing the previous boundary conditions one finds the following solutions for the equations of motion in the different radius intervals Then we impose that the two solutions match at R P and also that their first derivatives match at R P , obtaining the following conditions where we have introduced r λ ≡ λ − λ + . Then we insert into the action the solutions (B.10) and we integrate from r = 0 to r = R + with the appropriate potential (see (B.5)). From this computation we have to subtract the action for the case in which the field sits at the false vacuum from r = 0 to r = R + , that is we compute all in all (B.14) Using the equations in (B.13) we can rewrite the result as a function of the parameters of the potential to obtain The condition to select the first case (i.e. φ 0 > φ − ) can also be rewritten by employing again formula (B.13) as

(B.16)
We then analyse the second case. Imposing the boundary conditions we get the solutions By imposing matching of the fields and the derivative at R = R P we get the following equations for the unknown parameters R 0 , R P , R + We then compute the bounce action by inserting the solutions and integrating from r to R + , and after some rearrangements we get where R + and R 0 are functions of parameters of the scalar potential through the implicit equations (B.20). So we conclude that in both possible cases, the computation of the bounce action in the TBA approximation only needs to specify the critical points of the potential, that is the metastable vacuum, the peak of the barrier and the true vacuum. For the scenarios studied in this paper the first case (B.15) is the relevant one, and the validity condition is reported also in (3.6) (note that in the conventions in the main text we always choose the metastable vacuum to be at the origin of the field space). In this Appendix we have nevertheless reviewed both cases for completeness.
In the following we will compute the TBA bounce action by both evaluating the potential numerically and approximating it analytically. These approximations can be compared with the results of the full-fledged numerical bounce action computation.

B.1.1 Optimized triangular bounce
Studying the evolution of the bubble profiles for the actions computed numerically we note that the release point φ 0 is typically closer to the potential barrier than to the true vacuum at φ − in our setups. Therefore, we introduce here a modified version of the TBA that takes into account this feature, leading to a better agreement with the full numerics than the analytic formulas presented above. The idea is to allow the minimum of the potential to be a free parameter rather than to fix it at the true vacuum, allowing the TBA to more closely represent the shape of the potential, where the slope is closer to linear.
The TBA bounce action is computed by replacing φ − → φ 0 and V − (φ − ) → V − (φ 0 ), where φ 0 is now an arbitrary point along the potential in the interval φ eq < φ 0 < φ − (here φ eq > φ P is the point after the potential barrier where V (φ eq ) = V + ). This procedure defines a function S 3 /T (φ 0 ) that we can minimize over φ 0 in the allowed interval. The resulting minimum is the sought bounce action, that we have dubbed the optimal TBA.

B.1.2 Triangular bounce for the O'Raifeartaigh model with gauge interactions
In this section, we specify the discussion to the triangular barrier approximation for the model of Sec. 4.3 and compute the inputs needed for the TBA bounce action. We denote the usual combination of parameters as y F ≡ λF m 2 and y D ≡ gD m 2 . The local minimum at the origin is approximated as In blue we show the TBA computed using the numerical scalar potential and optimized with the procedures explained in the text. The red lines are the standard TBA approximations as described in Appendix B.1. The green line is the standard TBA evaluated on the analytic approximation of the scalar potential (as detailed in the text) and taking only the zeroth order term in the expansion for small r λ (as in Eq. (3.6)). This last approximation is the one used in Section 3 to derive analytic estimates.
The peak of the barrier is located at Finally, the true vacuum location and energy are given by where we set the renormalization scale µ = m. Within our approximation the thermal effects only enter at the origin, and at the top of the barrier, where they act to lower the potential relative to the true vacuum, and the potential difference between the top of the barrier and origin, respectively. In Figure 13 we consider two benchmarks with very different r λ at T n and show the bounce action S 3 /T as a function of the temperature, computed in different approximations. The black line is computed using the fully numerical thermal effective potential and the mathematica package "FindBounce" [85]. The blue line is obtained with the TBA evaluated on the full-numerical scalar potential and optimized with the procedure explained above. The red line is the TBA (as computed in Appendix B.1) evaluated on the full-numerical scalar potential. Finally, the green is the TBA evaluated on the analytical approximation of the critical points of the scalar potential as explained above, and moreover keeping only the leading order term in the small r λ expansion in Eq. (3.6). This is the approximation used to derive the analytic formulas in Section 3.
First, we see that the optimal TBA reproduces almost perfectly the numerical bounce computation. The standard TBA can predict well the location of T min but the overall normalization can be off up by a factor of ∼few, and in particular it does not agree with the numerical result in the vicinity of T c . However, the different trend in the overall shape of the bounce action in the two benchmarks (e.g. very flat around T min in the left one) is also captured in the standard TBA approximation.
Then, we note that the simplest approximation of the standard TBA reproduces very well the TBA when r λ is small. This was not obvious a priori since, besides expanding the TBA at leading order in small r λ , we have: i) assumed that the only temperature dependence is in the height of the potential at the origin; ii) employed the low-T approximation of the scalar potential to estimate it. The agreement between the red and green curve in the left panel of Fig. 13 hence confirms the fact that low-T is the correct approximation to employ, as discussed at length in Section 3. When r λ is larger (right plot) the simplest approximation (green line) clearly deviates from the standard TBA (red line), but nevertheless capture approximately the location of T min and the shape of the numerical results. It is important to observe that even if the normalization of the bounce action and its raising towards T c are not exactly reproduced by the approximations employed, they can still track the changes of the bounce action (shape deformations and overall size) as a function of the fundamental parameters of the model. This elucidates why the analytic estimates obtained in Section 3 can capture the scaling of T n in the different models as discussed in Section 4.2 and 4.3.

B.2 Single field approximation of the multi-field bounce action
As discussed previously, the model presented in Sec. 4.3, contains more than one dynamical degree of freedom that actually enters into the bounce action computation. Namely, the fields X, φ 1 andφ 2 vary along the minimal potential energy trajectory in field space which connects the two minima of the potential. This implies that the full bounce solution is that of a three field problem, which is in general only solvable numerically. In this Appendix we explain why in the model under study the bounce action can be effectively approximated with a one field problem and what is the regime of validity of such approximation.
The bounce action involving the three fields is where T i is the kinetic energy associated with each field, which is an additive quantity. In our analysis we approximate this action by minimizing V (X, φ 1 ,φ 2 ) along all three directions and by solving the bounce equation only for X. This corresponds to neglect the contribution from the kinetic energy of the other two fields φ 1 andφ 2 . Since the kinetic . Left: Bubble profiles for the single field and three field bounce actions as a function of the bubble radii. The release point for the single field (pseudo flat) direction is practically identical in both schemes, while the contribution coming from motion along the other field directions is generally small, but non-zero. Right: The bubble trajectory in field space for the single field and three field bounce action. The x-axis represents motion along the pseudo flat direction while the y-axis represents motion along the larger of the other two field directions (φ 2 ). It is clear from both plots that the motion along the pseudo flat direction is mostly unaffected by the motion along the other directions, and therefore the single field path approximation is viable. energy of these fields is related to the potential energy along the same directions by the equation of motion, we expect that the kinetic energy contributions of φ 1 andφ 2 can be consistently neglected if their VEVs along the bounce trajectory are small compared to the one of X. This is typically what happens in this model as we show in the left panel of Figure 14, where we plot the bubble profiles of the full three field problem computed using FindBounce (red, blue and green line). We also show for comparison (in dashed black) the bubble profile that we obtain for the single-field bounce solution for X, which is essentially identical to the one in the three-field solution. In the right panel of Figure 14 we show a two dimensional slice of the field path along the bounce solution, in the X,φ 2 plane (we showφ 2 since the φ 1 vev is smaller and hence it has even a smaller impact on the value of the bounce action). We see that the X trajectory is almost unchanged by the addition of the second field. We hence conclude that it is typically a robust approximation in this model to neglect the kinetic contribution of the fields φ 1 andφ 2 and to solve the one-field problem.
Nevertheless, we would like to estimate the range of validity of our approximation. As mentioned, we expect the difference to come from the kinetic terms of φ 1 andφ 2 , which will be non negligible if the size of the φ 1 andφ 2 VEV's compared to the one of X is not negligible. In particular, we would like to estimate the impact of this approximation in the overall bounce action. We hence use intuition from the TBA where the size of the bounce action is proportional to the cubic power of the field displacement. We define the following ratio R 1d/3d def = X 3 (r) X 2 (r) + φ 2 1 (r) +φ 2 2 (r) where X(0), φ 1 (0) andφ 2 (0) are the field distances from the origin computed at the release point r = 0. The ratio R 1d/3d provides a measure of the relative error of the single-field bounce action computation against the full three-field one. Indeed, we have also crosschecked numerically in several benchmarks that the difference in the values of the bounce action is negligible when R 1d/3d is small. In the main body of the paper, we will define the region where the 1-field approximation breaks down as when the quantity R 1d/3d > 0.5, corresponding approximately to a relative error of ∼ 50% of the single field approximation compared to the full three field solution.

C Sensitivity of GW interferometers
In this section we briefly discuss the interpretation and generation of the sensetivity curves used to define detection of a GW signal. We follow standard definitions and conclusions obtained in [103][104][105][106], see also [107] and references therein. The detection sensitivity for GW background for a given experimental setup, is given by the integrated signal-to-noise ratio (SNR) over an observation time interval t obs as where S is the mean signal, N 1/2 = S 2 IJ − S IJ 2 is the average noise, I, J = 1, 2 indicate coupled detectors, n det distinguishes between experiments that aim at detecting the signal by means of an auto-correlation of a single detector (n det = 1) or a crosscorrelation of a couple of detectors (n det = 2) measurement. The effective noise strain can be written in terms of the noise strain power spectrum D N (f ) and the frequency dependent detector response function R(f ). The latter quantities are the ones typically reported by the experimental collaborations. The detector response becomes more intricate in the case of correlated detectors [106], where an overlap detection function must be computed. We perform the appropriate procedure when computing the relevant SNR. The signal/noise strains can be rewritten in terms of the signal/noise spectrum density as where in this paper Ω S (f ) will be the energy density of a SGWB produced by the SUSYbreaking FOPT in the early universe redshifted till today (see Eq. (2.26)). The frequency interval (f min , f max ) is determined by the bandwidth of each experiment, typically related to the length scale of each detector. This integrated quantity will imply detection if it surpasses a predefined threshold value ρ 2 ≥ ρ 2 thr set by baysian probabilistic measures [103], usually taken between 3 and 10.
By defining a frequency dependent shape function which models the expected signal one can define the sensitivity of a given experiment even though this procedure does not provide a precise statistical indication of the expected sensitivity as the signal to noise ratio in Eq. (C.3). Depending on how well the signal shape is approximated, sensitivity curve based on the shape function approximation can provide an instructive visual tool for estimating detection in a frequency dependent way. In the following section we discuss the most common method of generating sensitivity curves for 1OPT GW signal detection.

C.1 PLI curves
Power Law Integrated (PLI) curves are generated by considering a power law function of the frequency f for the GW signal shape. The most common assumption, is a power law of the form where b is known as the spectral index. Taking this assumption we obtain the integrated SNR as [108,109] h Finally, we can define the PLI sensitivity curve by maximizing the integrated SNR over the spectral index b, that is which gives the threshold value for the signal at each frequency. A curve which crosses the threshold value at a given frequency will therefore represent a detectable signal, assuming an approximate power law behavior. Given that GW signals from FOPTs have a broken power law shape, the method described here is typically appropriate to visualize the sensitivity of a given experiment. In Fig. 15 we summarize the PLI curves for the different GW interferometers considered here.

C.1.1 Experimental parameters
For completeness we recompute here the PLI curves as described in the previous section for ground based interferometers such as The Laser Interferometer Gravitational Wave Observatory (LIGO) [16], the Einstein Telescope (ET) [17], the Cosmic Explorer (CE) [18,19] and the Atom Interferometer Observatory and Network (AION) [118], as well as for space based detectors such as The Laser Interferometer Space Antenna (LISA) [110], the Big  Figure 15. PLI curves for the different experiments considered here with the threshold value for the signal to noise ratio is conservatively fixed to ρ thr = 10. The required data to derive a PLI curve for every experiment are collected in Table 1 Table 1. Summary of the experimental parameters used in generating PLI curves for GW detection in this work. Auto/cross-correlation measurement is indicated by n det = 1(2), the observation time and bandwidth for each experiment are presented above. The detector response R(f ), multiple detector overlap function Γ IJ (f ) and noise strain spectrum D noise (f ), S noise (f ) are extracted from the references herein.
In order to determine the reach of a given experiment, we need to know the frequency band, the response function of the detector R(f ) within this band, the measured noise at every accessible frequency D I noise (f ), the time of observation t obs and the number of coupled detectors n det . We report the extracted parameters for the various experiments in the Table 1.