Pulsar Timing Array Stochastic Background from light Kaluza-Klein resonances

,


I. INTRODUCTION
A stochastic gravitational wave background (SGWB) signal can originate from the superposition of numerous independent gravitational wave events, often of astrophysical nature, that are too weak to be individually resolved.Alternatively, the SGWB can arise from intrinsically non-local and stochastic sources which may happen in the primordial universe.The detection of a SGWB in the nHz frequency band is a primary goal of pulsar timing array (PTA) experiments.In this frequency range, astrophysics predicts a SGWB signal primarily generated by inspiraling supermassive back hole binaries (SMBHBs) [1][2][3].According to current understanding of galaxy formation and evolution, the SMBHB population produces a SGWB with a strain following a power law characterized by a power index γ ≃ 13/3 [4] and an amplitude A ∼ 10 −15 at frequency of f = 1 year −1 [5,6].Large deviations from this prediction in the detected SGWB would indicate the presence of misunderstood astrophysics, or the breaking of the standard model of cosmology and particle physics.
A few years ago, the PTA collaborations released their analyses of the data collected over the past decade or so [7][8][9][10].These analyses independently revealed the presence of a common red noise contamination in the pulsar arrival timings.Common red noise is a clear signature of a SGWB signal if accompanied by evidence in Hellings and Downs correlation observable [11].The latter was absent in all sets of data, maybe due to insufficient statistics, or overlooked systematics, in the analyses [12].Furthermore, the SGWB that could have explained the observed common red noise would have poorly matched the power-law-like signal produced by SMBHBs for both, the power index [7][8][9][10], and possibly the amplitude [13].
With the additional PTA data collected in the past couple of years, significant progress has been made.For the first time, there is strong evidence of a detected (isotropic) SGWB [14][15][16][17].The common red noise has spectral shape compatible to the one reported in the previous data release, and Hellings and Downs correlation is now present [14][15][16][17].Intriguingly, the observed SGWB signal can be explained in terms of SMBHB population models only by stretching the population parameters slightly beyond what was previously considered reasonable [5,18], and the hypothesis of a SGWB sourced by first-order phase transitions (FOPTs), on the top or without the expected SMBHB SGWB, is more compatible with the NANOGrav data than the SMBHB-only hypothesis [19].
More specifically, a FOPT with gravitational waves primarily triggered by bubble collisions, strength α * ≳ 1, reheating temperature 10 −2 ≲ T * /GeV ≲ 10 1 , and inverse time duration β/H * ≲ 50 is one of the source candidates yielding the highest Bayes factor in the NANOGrav analysis [19].However, in qualitatively similar ranges of values of α * , β/H * and T * , also the FOPT SGWB generated by the sound-waves and turbulence contributions fits well the PTA observations [18,19].
Shortly after the PTA data release of a few years ago, various studies explored the possibility of explaining the observed red noise as a result of a cosmological FOPT [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] (for a recent review see e.g.[35]) 1 .The conclusions drawn from these studies are expected to remain valid even with these data update, as a few extra years of data unlikely could significantly alter the common red noise observed over more than a decade.However, with the evidence of the SGWB reaching 4-σ, it is timely to take a step forward.Indeed, despite the anticipation that PTA experiments were likely on the cusp of strengthening the previous results favoring new physics, explicit particle physics frameworks featuring the required FOPT have been scarce. 2This scarcity highlights the challenges faced by extensions of the Standard Model with FOPTs at the MeV-GeV scale in overcoming existing phenomenological constraints.
In this article, we investigate the feasibility of a BSM setup that we find particularly intriguing.We introduce a simple, yet versatile, warped extradimensional model adapted to the GeV scale.This model predicts unavoidable FOPTs and offers a phenomenology that, at present, remains comfortably consistent with, albeit not invisible to, current or forthcoming experimental searches.While more complex variations can be developed without compromising its fundamental characteristics, our proposed model serves as a promising foundation for investigating the implications of warped extradimensions in the GeV regime relevant for the PTA SGWB discovery.
Our analysis is structured as follows.In Sec.II we introduce the setup and summarize its signatures at some observables.In Sec.III we prove that the model can easily induce a FOPT compatible with the SGWB discovery announced by the PTA Collaborations and, in the parameter region fitting the PTA data, its main phenomenological signatures are broadly compatible with current constraints.In Sec.IV we study possible embedding options of the SM fields in the model, as well as their possible phenomenological signatures due to the interactions of these fields with the graviton and radion sectors of the model.We finally devote Sec.V to our conclusions.

II. THE MODEL
The candidate model that we propose as an explanation of the FOPT possibly at the origin of the PTA SGWB background, is a slight variation of the Randall-Sundrum (RS) model.The RS model has been the starting point for numerous studies, providing a natural solution of the hierarchy problem, a natural strong FOPT at the EW scale, and a rich phenomenology at the LHC.Here we review the main ideas of the model, and while skipping the intermediate steps of calculations known in the literature, we summarize some key results and adapt them to the scenario we are interested in.In this brief overview we follow Refs.[41][42][43] and references therein.

A. The model at zero temperature
The Randall-Sundrum (RS) model is a 5D setup built on Anti-de Sitter (AdS) spacetime with line element3 with Latin (Greek) indices running over the five (four) spacetime dimensions.Along the r direction, n Minkowskian 4D spaces B a , dubbed branes, are located at r = r a , with a = 0, I 1 , . . ., I m , 1, where B 0 is the UV (Planck) brane, B 1 is the IR brane, and we allow in principle a number m of hypothetical intermediate branes B Ii (i = 1, . . ., m), at intermediate scales between the Planck scale and the scale of the B 1 brane, which will be considered to be at the MeV-GeV scale.The total number of branes would then be n = 2 + m.Due to the warped factor A(r), the energy scale involved in the brane B a2 is exponentially suppressed with respect to the one of the brane B a1 , if a 1 < a 2 , with the suppression depending on the distance r a2 − r a1 .The RS setup thus provides an elegant solution for scenarios featuring hierarchical energy scales.In this paper we will consider two possible scenarios: a model with m = 0 (no intermediate brane), just like the RS scenario, and a model with m = 1, i.e. with an intermediate brane at the TeV scale, but of course more general scenarios can be considered for other purposes.The AdS metric exhibits conformal invariance along r.The localization of the branes requires a stabilization mechanism that breaks this symmetry.The Goldberg-Wise mechanism is a well-known realization of this breaking [44].It involves a (extremely heavy) scalar field, ϕ, that can propagate in the bulk (i.e. the space between the branes) and has bulk potential V (ϕ) and brane potentials Λ a (ϕ) = Λ a + 1 2 γ a (ϕ − v a ) 2 , with Λ a , v a and γ a (≫ 1 in the stiff limit) some positive constants (see Refs. [41,44] for details).In this case, the five-dimensional action of the model reads as where S GHY is the Gibbons-Hawking-York term canceling out the boundary terms of the 5D geometry variations, the induced metric on the branes being ḡµν = e −2A(r) η µν , and finally κ 2 ≡ 1/(2M 3 5 ) where M 5 is the 5D Planck scale, a parameter of order the 4D Planck scale M Pl .At this point, one can see that the potentials of ϕ induce some potentials for r a .In the AdS/CFT correspondence picture, this corresponds to the condensation of the SU (N a ) (conformal) gauge fields of the B a brane.
The values of r a of the potentials can be calculated by solving the equations of motion (EoM) for ϕ and integrating out the heavy degrees of freedom [44].Hereafter, we can limit ourselves to the phenomenology of the two most IR branes with a coordinate difference ∆r 4 .They indeed involve, respectively, the next-to-lowest and the lowest energy scales of the whole picture, and the study of their phenomenology suffices for the purpose of our analysis, namely the FOPT interpretation of the PTA SGWB. 5 As only the difference ∆r matters, we fix the lowest coordinate by convention. 6To determine the potentials of their positions, we solve the EoMs with the "superpotential" method [45].This method enables us to reliably explore the parameter region leading to sizeable backreaction on the metric [46].(In the region with small backreaction, the cosmology of the model tends to get stuck in the wrong phase [43,47]; see discussion later on).By introducing an expansion parameter u ≪ 1 [42], we obtain that r 1 develops the potential [42,44] Ū1,eff (r 1 ) = 2u 2 v2 1 (r 0 1 − r1 ) e 4A0(r 0 1 )−4A0(r1) − 1 e −4A0(r1) + 6λ 1 e −4A0(r1) , where r0 1 is defined by the condition v 1 ≡ v 0 e ur 0 1 , A 0 (r) is the leading order approximation of A(r) at r > 0, and finally va ≡ κv a and r1 ≡ r 1 /ℓ.Notice that in Eq. ( 3), λ 1 is a free negative parameter whose absolute value is in the range O(0.1 − 10).The position of the minimum of the potential Ū 0 1,eff sets the value r min 1 at which the brane B 1 is stabilized.
Alternatively, one can reparametrize the degree of freedom associated to the r 1 coordinate in terms of the radion scalar field χ 1 defined as For consistency, its potential V r1 will have minimum at ρ ≡ e −A0(r m 1 ) /ℓ.Moreover, after some algebraic manipulation of Eq. ( 3), one finds Ū1,eff (r 1 ( χ1 )) , where the function r1 ( χ1 ) is the inverse of Eq. ( 4) and N = (2ℓ) 3/2 π/κ is the number of colors of the SU (N ) symmetry of the AdS/CFT correspondence for the next-to-last brane.The described procedure is not specific of the behaviour of the brane at r 1 .Should we have been interested in energies processes not below the energy of the next-to-last brane, we could have repeated the procedure for the positions r 0 , r I1 , . . .as well, introduced radion fields for each position, and obtained expressions conceptually similar to Eqs. ( 3)-( 5) also for them.

B. The model at finite temperature
As mentioned earlier, our interest lies in the phenomenology and the FOPT associated with the lowest B 1 and nextto-lowest B N branes of the complete RS picture7 .Correspondingly, we limit ourselves to consider temperatures T below the scale of the brane B N .Only the dynamics of the B N and B 1 is then relevant, and we now focus on it (see Ref. [43] for the finite-temperature picture in the case with several branes).
In quantum field theory at finite temperature, the temperature T replaces the time coordinate and the x 0 component of the spacetime gets compactified in a circle.In this regime, the 5D gravity action leading to the AdS solution in Eq. ( 1) now admits two geometric/gravitational solutions [47,48]: the previous AdS solution of Eq. ( 1) where the x 0 direction is compactified, say ds 2 RS,T its line element, and the so called Anti de Sitter-Schwarzschild (AdS-S) solution.The latter corresponds to the metric of black hole with line element with the event-horizon singularity at r = r h > 0. Notice that both solutions are equivalent at r → ∞.The existence of the brane of next-to-lowest energy scale, B N (with r N stabilized), is then allowed in both cases.Instead, the presence of brane B 1 at r = r 1 is permitted in both solutions only when r h > r 1 .Indeed, roughly speaking, for r h < r 1 , in the physical space containing the brane B N , the event horizon masks the brane B 1 .Then, when the zero-temperature RS action is heated up to the AdS-S background, the relevant 5D action is the one in Eq. ( 2) with the B 1 terms omitted.
The degrees of freedom in the AdS-S and RS T phases are then different: in the former, the light fields localized in L 1 disappear while the SU (N ) symmetry is restored; in the latter, the SU (N ) fields condensate and the fields of the lightest energy scale emerge.Due to this feature the AdS-S and RS T phases are also respectively dubbed as "deconfined" and "confined" phases.
Similarly to what happens to r 1 in the zero-temperature case, the gravity action itself does not specify the value of r h ; it is the backreaction of the non-gravity content that stabilizes it.The metric solution with r h at its minimum r min h is the "on-shell" (or classical) gravitational solution.This implies that, if the scalar field T h takes the role of the degree of freedom along the direction r h (and hence the potential of r h is translated into a potential of T h ), at the minimum of the T h potential we must find ⟨T h ⟩ = T , with T being the Hawking temperature of the black hole.In analogy with the radion, we define 8 To determine the position at which r h stabilizes, we follow Ref. [42].As previously done for the zero temperature case, we introduce a small expansion parameter u and use the "superpotential" method [45] to solve the EoM.In the solution we take into account that the blackening factor h(r) must satisfy the boundary and regularity conditions h(0) = 1 and h(r h ) = 0.This leads to and dA 0 (r) from which, after multiple algebraic manipulations, one can derive the finite-temperature effective potential, or equivalently, the free energy of T h .Here we report on its value at the minimum, which is located at T h = T as expected (see Ref. [42] for the full expression): The quantity F BH min (T ) is the "geometric contribution" (i.e.due to the singularity) to the total free energy of the model in the AdS-S phase [47].No similar contribution arises in the free energy of the AdS phase.However, a further contribution comes from the plasma made of the field content of the two phases.Accounting for them leads to [49] where g eff d/c is the number of relativistic degrees of freedom in the deconfined/confined phase, and is the potential gap between the two phases in the limit T → 0. The main quantities controlling the free energies are then: the value N of the SU (N ) conformal group entering in F AdS−S min (T ), the radion potential parameters v 0 , v 1 , λ 1 , ρ modulating E 0 , and the field content of L bulk , L 0 and L 1 which sets g eff d and g eff c .At this point one may already correctly guess that the details of the brane and bulk Lagrangians will play a minor role in our final conclusions because the field content only appears in the number of relativistic species.

C. The phase transition picture
As just explained, the model in Eq. ( 2) has two competitive phases at temperatures below the next-to-last energy scale of the model, or equivalently for r h > r N .Both phases have the brane B N at r = r N , but the brane B 1 only exists in the AdS phase with the value of r 1 already stabilized.It is like that the event horizon at r h somewhat hides the brane B 1 at r 1 when r h < r 1 .It is however worth stressing that the two geometries match at r → ∞ as well as at r 1 = r h → +∞.In the coordinate plane {r h , r 1 }, the deconfined and confined phases can then be described as 9 {r min h , +∞} and {+∞, r min 1 }, and one can move from the former to the latter via the path {r min h , +∞} → {+∞, +∞} → {+∞, r min 1 } along which the metric is kept regular.In terms of the fields χ and T h , this path is equivalent to To determine the value of the free energy at the intermediate step, one can notice that the configuration {T h = 0, χ 1 = 0} is an on-shell solution when T = 0, which implies that the free energy at the intermediate step is equal to F d (0) and is hence larger than the free energies of the final and initial phase.The phase transition between these two phases is then of first order.In particular, the direction of the phase transition occurring along the history of the universe, is the one above described, since |F c | ≪ |F d | at high enough T (but still fulfilling the condition r h > r N ) 10 . 8We remind that T is assumed to be below the energy scale of Ba controlled by the distance ra − r 0 (where ra is the position of any other brane whose energy scale is known).This is consistent with the described picture taking r h > r 0 as T is also the position of the potential minimum of T h and, due to Eq. ( 7), sets the scale of r h . 9Here we denote by r min 1 the value of r 1 at the minimum of the potential, while r min h is given by the value of r h at the minimum of the potential as a function of T h in the AdS-S phase where T min h = T and thus r min h /ℓ = − log(πℓT ). 10 If one considered generically high temperatures, the additional phases with r I i < r h < r I i+1 , would enrich the picture.Multiple FOPTs would then be possible [28], but if there is the hierarchy r I i ≪ r I i+1 , only one would dominate the SGWB at PTA frequencies.This is why we focus only on one of them in the present analysis.
A peculiarity of this FOPT is that, potentially, a short epoch of inflation can precede the onset of the bubble nucleation, and the phase transition can end up with a sizable reheating.Indeed, the considered FOPT is often very strong, leading to a sizable supercooling and thus to T n ≪ T c , with T c implicitly defined by F c (T c ) = F d (T c ).In the deconfined phase, inflation starts at the temperature T i , obtained by imposing that the energy density in the deconfined phase, be dominated by the vacuum energy E 0 .This yields which leads to N e = log(T i /T n ) e-folds of inflation provided T i > T n .This entropy injection would dilute all abundances, in particular the baryon-over-entropy ratio which is fixed at the BBN epoch to η B = n b /s ∼ 10 −10 .Assuming that baryogenesis takes place at the electroweak phase transition (electroweak baryogenesis), or earlier (e.g.leptogenesis), this means that the entropy injection triggered by the FOPT has to be taken into account in the baryogenesis mechanism, and the required value of the generated η B has to compensate for this late entropy production.
Moreover, the energy balance requires that at the end of the FOPT, the energy density in the deconfined phase, ρ d (T n ), is converted into radiation density in the confined phase, ρ c (T R ), which implies11 where T R is the final reheating temperature.
Overall, these expressions show that T R , T c and T n are mildly sensitive to variations of g eff c and g eff d due to the fourth power dependence.This is an extra hint to the fact that our numerical results of the FOPT are almost independent of the particular model setup that specifies L bulk , L 0 and L 1 .

III. ANALYSIS OF THE PHASE TRANSITION AND COMPATIBILITY WITH THE SGWB PTA DATA
A FOPT is characterized by the following thermodynamic quantities: the inverse time duration β/H * , the strength α * , the nucleation temperature T n .From these, other quantities follow, e.g.T R in Eq. (13).We compute the thermodynamic quantities by solving numerically the O(3) and O(4) bounce equations along the aforementioned path The procedure is standard in the FOPT literature and we omit its description (see Ref. [42] for the details of the procedure we adopt).
Once nucleated, the FOPT bubbles expand and eventually collide.During this process, they perturb the plasma from its equilibrium.Both the bubble collisions and the plasma motion generate a SGWB.As the plasma motion can occur via turbulence and sound waves, in total three dynamical mechanisms contribute to the final SGWB production: the "bubble collision", the "sound waves" and the "hydrodynamic turbulence" [50,51].For each of them, the produced SGWB frequency shape at the time of production has its own dependence on α * , β/H * and v w , with v w being the bubble expansion velocity.The temperature T R enters to take into account the redshift due to the expansion of the Universe between the times of the signal production and detection [51].
In their recent analyses, the NANOGrav, EPTA and Indian PTA experimental collaborations also provide the interpretation of the observed SGWB signal in view of a FOPT source.They consider each mechanism once at the time.In principle, these mechanisms are not mutually exclusive so that combinations of them, weighted by some efficiency factors, are possible.This assumption is however reasonable for the PTA analyses.Indeed, the PTA measurements are sensitivity to a rather narrow GW frequency band and, inside it, a single contribution does not dominate over the others only for peculiar FOPT combinations.Likewise, we analyze one SGWB contribution at the time too.
Specifically, the NANOGrav Collaboration considers the "bubble collision" and "sound wave" contributions with their SGWB predictions parametrized in terms of the thermodynamic parameters α * , β/H * , T R , H * R * = 3 √ 8π/(β/H * ) with the assumption v w = 1 [19]. 12The analysis considers both the presence or absence of a SMBHB SGWB FIG. 1. Scatter plots of the inputs for the RS framework parameters (λ1, N , ρ) and the resulting FOPT quantities (TR, H * R * , Tn, Ti).Blue points fall in the NANOGrav 95% favorite region of the "bubble collision only" hypothesis, the region covered by blue + red points corresponds to the "bubble collision + SMBHB" hypothesis, and the gray points belong to none of the previous regions.In the top panel, the thin (outer) and thick (inner) lines enclose the NANOGrav 95% C.L. favorite regions for the "bubble collision + SMBHB" hypothesis and "bubble collision only" hypothesis, respectively.component accompanying the FOPT SGWB contribution.We report their results in the top panels of Figs. 1 and 2. In particular, the thin-dark line in the top panel of Fig. 1 encloses the NANOGrav 95% C.L. favorite region for the "bubble collision+SMBHB" hypothesis in the parameter space {T R , H * R * } (after marginalizing over α * ).The thick-dark line does the same for the "bubble collision only" hypothesis, and "sound wave only" hypothesis in Fig. 1 and Fig. 2, respectively.We do not consider the NANOGrav 95% C.L. region for the "sound wave+SMBHB" as it practically provides no constraint (see Fig. 8 of Ref. [19]).For the sake of readability, we do not display the 95% C.L. regions in the planes {α * , T R } and {α * , H * R * }, but we impose their fulfillment.In the regime relevant in our analysis, namely α * ≫ 10, this corresponds to request −1.5 < log 10 (H * R * ) < −0.85 and −2.8 < log 10 (T * /GeV) < −1.6 for the "sound wave only" hypothesis, and no further constraint for the "bubble collision+SMBHB" or "bubble collision only" hypotheses [19].
The analysis in terms of the "hydrodynamic turbulence" contribution is carried out by the European PTA and Indian PTA Collaborations [18].They consider the 'hydrodynamic turbulence" SGWB expressed as a function of the thermodynamic parameters β/H * , T R and an effective parameter Ω * encapsulating the spectrum dependence on To determine whether, and in case in which parameter region, the radion FOPT provides a fundamental explanation of the PTA SGWB signal, we perform a scan on the parameter space of the model.Specifically, we explore the parameter region N ∈ [10,30], log 10 (ρ/GeV) ∈ [−3, 1], −λ 1 ∈ [0.3, 2.7], log 10 (v 0,1 ) ∈ [−2, 2] with flat priors.For concreteness, we fix g eff c (T ) = g eff d (T ) ≈ g eff SM (T ).As previously stressed, our results are stable under (reasonable) variations of g eff c and g eff d (cf.Eqs. ( 12) and ( 13)).Order-of-magnitude variations of v0 and v1 have no large impact of the FOPT thermodynamic parameters either.A strong dependence on v0,1 instead appears in the ratio m rad /ρ, and it is approximately given by m rad /ρ ≈ (v 1 /15) log (v 1 /v 0 ).Moreover, the parameter u of the model turns out to have a dependence with v0,1 which is given by u = 1 r0 1 log (v 1 /v 0 ) with r0 1 ≈ 40.Then m rad /ρ ≈ 2.7v 1 u, so that the radion mass vanishes in the limit u → 0. Notice that we have considered v1 > v0 in the parameter region to ensure that m 2 rad > 0 and u > 0. For each point of the scan, we determine the FOPT thermodynamic parameters.Figures 1, 2 and 3 summarize the outcome of the analysis.The gray [red] {blue} points show the results for whole sample [the fraction of the sample in the 95% C.L. region of the FOPT+SMBHB hypothesis] {the fraction of the sample in the 95% C.L. region of the FOPT-only hypothesis}.Every blue point overlaps a red and a gray point, and every red point overlaps a gray point.We do not report on α * as it always turns out to be much larger than 10 (the FOPT phenomenology is insensitive to variations of α * at α * ≳ 10).We do not calculate v w but, as argued in Ref. [42], it is reasonable to expect v w ≃ 1 in a radion FOPT with α * ≫ 10.Neither we determine the efficiency factor suppressing the GW production from turbulence, as its determination is unknown in much simpler scenarios than the radion FOPT.We postulate it to be maximal for simplicity, and leave for future studies the inclusion of this extra parameter in our analysis.
Figures 1, 2 and 3 carry a wealth of information.The energy scale of the brane B 1 , represented by ρ, must be in the MeV-GeV range.Such a scale is even smaller if the sound-wave regime applies.This may constitute a challenge for some models built upon the generic RS setup here investigated (we remind that our FOPT results are practically independent of L bulk and L a ).In fact, it typically turns out ρ ∼ T R ∼ (10 2 − 10 3 )T n and, on top of this, some e-folds of inflation before the FOPT are unavoidable.This requires some caution to guarantee the SM sector to reach to the standard-cosmology conditions before BBN.We also highlight that there seems to be slight preference (i.e. higher density) for low N and large λ 1 , but different priors should be tested before clarifying this aspect.

IV. POSSIBLE EMBEDDINGS AND THEIR PHENOMENOLOGICAL SIGNATURES
As previously explained, the thermodynamic characteristics of the radion FOPT are almost independent of the particular embedding, namely L bulk and L a .The phenomenology instead strongly depends on it.Before identifying some promising embedding options, we summarize some key results and ideas useful for our purpose.Given that the scale of the IR brane B 1 is at the MeV-GeV scale, the bulk Lagrangian L bulk cannot contain any SM fields as their KK resonance masses would be too low and should have already been discovered at LHC and through electroweak precision measurements.In other words we will concentrate the SM in some 4D brane with cutoff higher than the TeV scale.Two simple options appear at this level.
1.The simplest option is having a model with two branes (n = 2), the UV brane B 0 (at the Planck scale) and the IR brane B 1 (at the sub-GeV scale).The SM should be localized in the UV brane, with Planckian cutoff, so that the running of SM couplings is logarithmic and unification of gauge couplings can take place at high scales provided the appropriate matter is added to the SM.In this model the warp factor does not solve the hierarchy problem, but the role of the warping is precisely to create the scale of the IR brane for which the SGWB is eventually being detected at PTA Collaborations.The hierarchy problem could be mitigated e.g. if the SM in the UV brane is completed by supersymmetry, in which case we will have the bonus of gauge coupling unification at scales ∼ 10 16 GeV.In the brane B 1 only some hypothetical dark sector can eventually live.Only the graviton (and radion) sectors are propagating in the bulk of the extra dimensions with KK modes at the sub-GeV scale.
2. An alternative option is a model with three branes (n = 3), the UV brane B 0 (at the Planck scale), the most IR brane B 1 (at the sub-GeV scale) and an intermediate brane B T (at the TeV scale) 13 .Models with multiple intermediate branes have been studied e.g. in Refs.[52,53].The three branes are located at r 0 ≪ r T ≪ r 1 , or using conformal coordinates kz ≃ e kr , at , where ρ T ≃ TeV.In this model the SM is located in the B T brane, so that the hierarchy problem is solved by the warp factor.This model is an IR extension of the RS model.The graviton is propagating in the whole bulk providing graviton KK modes at the sub-GeV scale.There are two radions: the radion r T (x) which stabilizes the B T brane at ⟨r T ⟩ = ρ T , at temperatures T ∼ ρ T , with support in the region [r 0 , r T ] and KK masses in the sub-TeV scale; and the radion r(x), which stabilizes the B 1 brane with respect to B T at ⟨r⟩ = ρ at temperatures T ∼ ρ, with 13 For notational convenience we are renaming the intermediate brane B I at the TeV scale in this section as B T .
support in the region [r T , r 1 ] and KK masses in the sub-GeV scale.For hierarchically related phase transitions, i.e. for ρ T ≫ ρ, we can assume that, at temperatures T ≳ ρ, the first phase transition has already settled the heavy radion r T (x) in its vacuum expectation value ⟨r T ⟩ = ρ T , while keeping the light radion in the symmetric phase ⟨r⟩ = 0, just leaving a vacuum energy, so that the second transition may proceed without any influence of the heavy radion, which has been integrated out.We then assume a two-step phase transition.The brane B 1 is then stabilized by the radion r(x) leaving the SGWB detected by the PTA Collaborations.
A. The graviton sector The model has a graviton sector, which propagates in the whole bulk, with a massless zero mode and an infinite number of resonances with masses m n such that m 1 ≃ 3.9ρ . . . .It generates a deviation of the Newtonian potential given by [54,55] V where which is dominated by the first KK mode.The experimental constraints on Yukawa-type deviations from Newtonian gravity yield for KK gravitons 1/m 1 ≲ 30 µm [56][57][58], or This is a very mild constraint as a consequence of the exponential suppression of the correction.
The graviton coupling to matter in the brane B b is given by The equations of motion of the graviton sector can be found e.g. in Ref. [59].There is a massless zero mode, the physical graviton, while the mass spectrum of the KK modes is given by m n = x n ρ, where x n is determined by the condition J 1 (x n ) = 0.As for the wave function of KK modes, after making the expansion in modes The coupling of KK modes to the B b brane is then written as As for the possible three branes we can have in the above models, the coupling is: This coupling is very mildly suppressed for the present models where ρ is at the sub-GeV scale.In RS-type models where the SM is localized at the B 1 brane and ρ is at the TeV scale, this reproduces the usual result that the KK modes are coupled with the strength ∼ 1/TeV.In our models we only assume that a dark sector could be localized on the B 1 brane.
• For z = z T we have where we have normalized the ratio k/M Pl to 0.1, and made use of the property J 2 (x) = x 2 /8 + O(x 4 ) valid for x ≪ 1.The prefactor of ϵ n (z T ) gives the usual 1/TeV coupling of KK gravitons in RS models, while ϵ n (z T ) is an extra suppression factor.For instance, for the first KK mode we have that x 1 ≃ 3.9 and J 2 (x 1 ) ≃ 0.4 so that the suppression factor, for ρ T ≳ 1 TeV and ρ ≲ 1 GeV, is given by ϵ 1 (z T ) ≲ 5 • 10 −10 .This tiny coupling makes it possible to evade present ATLAS bounds of m 1 ≳ 2 TeV [60], valid for ϵ 1 (z T ) = 1.
• For z = z 0 the suppression factor stems from Eq. ( 21) by replacing z T → z 0 .It is so tiny (ϵ n (z 0 ) ∼ 10 −55 ) that in practice graviton KK modes are decoupled from the SM, when the latter is localized in the UV brane B 0 , as in model 1) above.

B. The radion sector
The radion corresponds to scalar perturbations F (r, x) of the metric as ds 2 = e −2(A+F ) η µν −(1+2F ) 2 dr 2 .Its coupling to matter localized on the B b brane is given by In the limit of small backreaction the relation between the canonically normalized 4D radion r(x), which stabilizes the B 1 brane, and the metric perturbation F (z, x) is given by [61] where Λ r = √ 6ρ T provides the radion coupling strength in the RS model, and ϵ(z b ) is an extra suppression factor.The mass of the radion is given by m 2 rad = ρ 2 /Π rad .The function Π rad is a cumbersome expression with dΠ rad /dλ 1 > 0 and Π rad (λ 1 ) ∼ 10 2 −10 4 for λ 1 ∈ [−10, −0.1], but a simple phenomenological good enough approximation is provided in Sec.III.We refer the reader to Eq. (31) in Ref. [42] for the full expression that we adopt in our analysis.
• For z b = z 1 we can see that the suppression gets reduced from the typical RS result, as the coupling of the radion has the strength 1/ρ.This coupling would only apply to a hypothetical dark sector localized on the B 1 brane.
• For z b = z T we can see that the suppression factor is given by This radion decays into the SM matter in model 2) above, for which the SM is localized in the B T brane.In particular its coupling to two gluons gg (leading to the main decay channel) is governed by the interaction Lagrangian given by [61,62] where we have included the contribution from the localized trace anomaly (first term in c), and from triangular diagrams involving quarks Q such that 4m 2 Q ≫ m 2 r (second term in c), considering the second and third generation quarks.The decay width, and lifetime, into gg is given by where, in the last equality, we have used the numerical values Λ r = 2.5 TeV, m r = 0.1ρ, and normalized ρ to 1 GeV.Investigating the phenomenology and cosmology of such light feebly coupled radion is beyond the scope of the present paper and should be done in a separate publication.
• Finally for z b = z 0 , which is relevant for model 1) above, the coupling of the radion to the SM, localized in the B 0 brane, has a huge suppression factor given by ϵ(z 0 ) = z 2 0 /(z 1 z T ), and taking typical values ρ ≃ 1 GeV, ρ T ≃ 1 TeV we get ϵ(z 0 ) ≃ 10 −33 and the previous calculation would lead to a lifetime of the radion for the γγ channel larger than the present age of the universe.This particle is then stable in cosmological times and could play the role of dark matter provided its abundance can be conveniently controlled.However if there is some dark sector in the B 1 brane (containing e.g.dark matter) it could easily decay into it.Again, constructing a complete model with a dark sector is not going to affect appreciably the FOPT and is outside the scope of the present paper.
As we aim at addressing the PTA observation, we do not investigate extensions of this proposal.Adding new light fields in the B 1 brane may turn a viable option to address the Dark Matter puzzle, and many of the 4D extensions of the SM investigated in the literature can be accommodated in the UV brane.Most of such extensions are not expected to qualitatively change our main conclusions obtained in the minimal scenario.

V. CONCLUSIONS
Recent results on the detection of a SGWB in the nHz frequency band by the pulsar timing array collaborations lead, as a possible BSM explanation, to the existence of a strongly-first order (dark) phase transition with scale in the GeV range.We have first proposed in this paper a modelization of such phase transition in a simple warped five-dimensional model with a UV brane, at the Planck scale, and a dark brane at the GeV scale.The phase transition is the confinement/deconfinement phase transition where the radion field fixing the inter brane distance undergoes a transition from its deconfined state (above the nucleation temperature) in the 5D BH configuration, to its confined state (below the nucleation temperature) in the thermal 5D warped space.The model can be considered as a discretization of the RS2 model [63] by means of the introduction of a brane at the GeV scale.This model cannot support SM fields propagating in the bulk of the extra dimension, so that the SM (as in the RS2 model) has to be localized in the UV brane.On the other hand the dark brane can contain a dark sector, which the radion field interact with, and might contain dark matter candidates.The details of the matter content of the dark sector are not relevant for the present paper (as soon as they interact with the SM only through gravitational interactions) so we do not specify them here.As it is clear from its conception, as the SM is localized on the UV brane, the model does not solve the hierarchy problem (whose solution would more conventionally require the introduction of supersymmetry in the UV brane) but just introduces the dark brane to make contact with the recent PTA results.
Furthermore, models with multiple hierarchies [53] have been studied in the literature, and can be equally used to describe the PTA data along with solving the hierarchy problem.For instance between the UV and dark brane, one could introduce a TeV brane, where the SM can be localized, as in the original RS framework.In that case there should be two strongly-first order phase transitions induced by the two radion excitations (describing the fluctuations of the TeV and dark branes, respectively): i) a heavy sub-TeV radion, with support in the region between the UV and TeV scales; and, ii) a light sub-GeV radion, with support in the region between the TeV and dark branes.As gravitons propagate along the whole extra dimension, their KK modes are at the GeV scale.The first phase transition, driven by the heavy radion r T (x), should happen at temperatures of the TeV and reproduces the confinement/deconfinement phase transition studied in RS models [47].This FOPT should generate a SGWB that should be detected at LISA and ET [41].The second confinement/deconfinement phase transition, driven by the light radion r(x), should happen at temperatures of the order of the GeV.The latter phase transition, that has been analyzed throughout this paper, is straightforward as we can integrate out the GeV (or heavier) KK modes and remain with the SM and the light radion as the only degrees of freedom.However for the former one, as the heavy radion is normally heavier than the graviton KK modes, has to happen in the presence of the irreducible background of the graviton KK modes.They can affect the phase transition in a way that has to be thoroughly taken into account.The details of this analysis are postponed for future investigations.

FIG. 2 .
FIG. 2. As in Fig. 1 but for the "sound wave only" hypothesis.The "sound wave + SMBHB" hypothesis is omitted as its 95% C.L. favorite region practically covers the whole, considered parameter region.In the top panel, the grey points enclosed by the black line are not in blue because they fall outside the NANOGrav 95% C.L. region in the plane {α * , TR} and {α * , H * R * }, which reduces to −1.5 < log 10 (H * R * ) < −0.85 and −2.8 < log 10 (T * /GeV) < −1.6 when α ≫ 10.