Multi-brane cosmology

5D warped extra dimension models with multiple 3-branes can naturally realize multiple hierarchical mass scales which are ubiquitous in physics beyond the Standard Model. We discuss cosmological consequences of such multi-brane models with stabilized radions. It is confirmed that for temperatures below the scale of the IR brane at the end of the extra dimension, we recover the ordinary expansion of the Universe, with the Hubble expansion rate determined by sum of the physical energy densities on all 3-branes where they are localized. In addition, we explore the cosmology for temperatures above the scales of the intermediate and IR branes where the Universe is described by a spacetime with the 3-branes replaced by an event horizon. As the temperature of the Universe cools down, phase transitions are expected to take place, and the intermediate and IR branes come out from behind the event horizon. The Goldberger-Wise mechanism for radion stabilization has a well-known problem of having a supercooled phase transition, which typically does not get completed in time. This problem is even more severe when an intermediate brane is introduced, whose scale is well above TeV, as the corresponding Hubble rate is much larger. We circumvent the problem by employing an alternative mechanism for radion stabilization with dark Yang-Mills fields, which prevents a long supercooling epoch, but still allows the strong first order phase transitions. As a result, the phase transitions in our multi-brane Universe predict a stochastic gravitational wave background with a unique multi-peak signature, which is within the sensitivity reach of future space-based gravitational wave observers. We also show that there are $N-1$ radions for an $N$ 3-brane set-up, unlike a recent claim that there exists only one radion.


Introduction
Theories of physics beyond the Standard Model (SM) often contain new energy scales hierarchically different from the electroweak and Planck scales.Such multiple hierarchical energy scales are naturally realized by introducing new 3-branes with positive tensions into the 5D Randall-Sundrum (RS) spacetime [1] bounded by two 3-branes with positive and negative tensions called UV and IR branes.3-branes have their typical energy scales exponentially different from each other.Warped extra dimension models with multiple 3-branes and their phenomenological applications have been discussed in refs. .According to the AdS/CFT correspondence [27][28][29], the RS model is a holographically dual description of a nearly-conformal strongly-coupled 4D field theory [30,31].Then, the introduction of a new 3-brane corresponds to an extra spontaneous breaking of the conformal symmetry via confinement in the dual 4D picture, i.e. after the first spontaneous breaking of conformal symmetry, the theory flows into a new conformal fixed point [18].The conformal symmetry is spontaneously broken again at the IR scale, and the theory presents two distinct phase transitions.A multi-brane model contains multiple radions because each distance between two 3-branes corresponds to a modulus field.The stabilization of all the radions has been recently established in ref. [32] by a simple extension of the Goldberger-Wise (GW) mechanism [33] introducing a bulk scalar field with brane-localized potentials.1If a multi-brane model is realized in nature, it must predict a consistent cosmological history of the Universe.In the case of the RS model with two 3-branes, there was a question of whether the late cosmology from the Big Bang Nucleosynthesis (BBN) to the present deviates from the ordinary Friedmann-Lemaître-Robertson-Walker (FLRW) Universe.The question has been settled in refs.[36,37].If the radion is stabilized correctly, the conventional FLRW equations can be recovered in the 4D effective theory of the RS model for temperatures below the scale of the IR (TeV) brane.It is then natural to ask if a multi-brane model with stabilized radions can reach the same conclusion, which is discussed in the first part of the present paper.We will confirm that for temperatures below the scale of the IR brane the cosmology of a multi-brane model can reproduce the ordinary FLRW Universe, with the Hubble expansion rate given by the sum of physical energy densities on all the 3-branes.
For temperatures above the typical scale of the IR brane in a multi-brane model, how does the Universe evolve?It has been known in the RS model that the system at high temperatures is described by an AdS-Schwarzschild (AdS-S) spacetime where the IR brane is replaced by an event horizon [38].As the temperature of the Universe cools down, a phase transition from the AdS-S spacetime to the RS spacetime is expected to take place.This phenomenon corresponds to a transition from the deconfined phase to the confined phase in the dual 4D CFT picture.Due to the scale invariant nature, the phase transition generally proceeds through a supercooling phase.Ref. [38] has pointed out that the RS model with the radion stabilized by the GW mechanism suffers from a long supercooling phase which prevents the transition from being completed when a classical treatment of gravity is reliable and also a backreaction effect on the gravitational action from brane-localized potentials of the GW scalar field is negligible.This problem has been analyzed in detail recently [39], and there are several interesting discussions and possible solutions [40][41][42][43][44][45][46][47][48][49][50][51][52].The resulting phase transition predicts the production of a stochastic gravitational wave background probed by future space-based gravitational wave observers such as LISA [53], DECIGO [54], BBO [55], TianQin [56,57] and Taiji [58].
In a multi-brane Universe, it is natural to expect that phase transitions take place and intermediate 3-branes and the IR (TeV) brane come out from behind the event horizon as the temperature cools down (see Fig. 1).For each phase transition to proceed, the bubble nucleation rate has to compete with the Hubble rate at the temperature during the phase transition.Since an intermediate 3-brane has a characteristic energy scale higher than the TeV scale and hence the Universe has a larger Hubble rate during the phase transition, the long supercooling problem becomes more severe.To circumvent the problem, in the present paper, we generalize the radion stabilization mechanism proposed in ref. [51] to a multi-brane model by introducing a dark Yang-Mills gauge field in each subregion of the bulk between two 3-branes.Such a new radion stabilization mechanism makes it possible to avoid a long supercooling epoch, and still allows the strong first order phase transitions.Since each phase transition associated with a 3-brane coming out from behind the event horizon produces a stochastic gravitational wave background with some peak frequency, which is not diluted away by subsequent phase transitions, our multi-brane model predicts a unique multi-peak gravitational wave signature detectable at future space-based observers.
The rest of the paper is organized as follows.In section 2, we demonstrate that one recovers the standard Hubble expansion law of the Universe for a multi-brane model with stabilized radions.Section 3 discusses phase transitions in a multi-brane Universe where 3-branes come out from behind the event horizon.In this section, we focus on the GW mechanism to stabilize radions and show that the bubble nucleation rate for a phase transition associated with an intermediate brane is too small for the transition to get completed.Then, in section 4, we consider a new radion stabilization mechanism with dark gauge fields to address the issue.Section 5 then illustrates a unique gravitational wave signature produced during the phase transitions in our multi-brane model.Finally, we give conclusions and discussions in section 6. Appendix A emphasizes the existence of two radions in the three 3-brane setup while some author claims that there is only one radion.

Late cosmology
We consider a 5D warped geometry with three 3-branes and study its impact on the late cosmology from the BBN to the present.The discussion is a generalization of that of the original RS model with two 3-branes [36,37].We first present our framework and then analyze conditions under which the late cosmology can reproduce the usual FLRW Universe in the presence of a radion stabilization mechanism.In particular, we study conditions under which the standard Hubble law is recovered.

The framework
Our spacetime geometry is given by R 4 × S 1 /Z 2 with a background metric [36,59], where A, B = 0, i, 4 with i = 1, 2, 3 and x 4 ≡ y, δ ij denotes the Kronecker delta, y ∈ [0, 1/2] is the coordinate of the S 1 /Z 2 orbifold, and n(t, y), a(t, y) and b(t, y) are all functions of the time coordinate t and y.They are determined by the Einstein equation, Here, G AB denotes the Einstein tensor, T AB is the energy-momentum tensor, and κ 2 ≡ 1/(2M where the prime (overdot) denotes a derivative with respect to y (t).We introduce UV and IR branes which extend over R 4 and reside on the two orbifold fixed points at y = 0, 1/2, respectively.In our setup, another intermediate 3-brane is placed at y = y I with 0 ≤ y I ≤ 1/2.The UV, intermediate and IR branes have tensions, V * , V I and V 1/2 , respectively.The energy-momentum tensor on the branes is then given by Here, ρ * , p * denote the energy density and pressure of matter localized on the UV brane, and ρ I , p I and ρ 1/2 , p 1/2 are defined similarly.Due to delta-functions in Eq. (2.4), the Einstein equation (2.2) is solved by decomposing it into singular and non-singular parts.
The functions a , n will be singular, and a , n will have finite jumps at y = 0, y I , 1/2.We can relate them to the localized energy densities as jump conditions.We call a bulk subregion surrounded by the UV and intermediate (intermediate and IR) branes as subregion 1 (subregion 2) and assume bulk cosmological constants, T subregion 1,2 AB = Λ 5  1,2 g AB , for the subregions 1, 2.
Let us now look at a static background solution for the current multi-brane setup, taking the limit of ρ * , ρ I , ρ 1/2 , p * , p I , p 1/2 → 0. We begin with the following ansatz for the metric (2.1): (2.5) with curvature scales k 1,2 , and some constant b 0 setting the compactification radius of the extra dimension prior to orbifolding.We can see that a(t, y), or a(y), is a continuous function with kinks at the branes.Note that there are two radion modes associated with two intervals of the subregions 1, 2, which will be considered in more details in section 4. For the purpose of studying the late cosmology, the current parametrization of the background metric suffices.The jump condition for a(y) at each 3-brane is then obtained from the '00' component of the Einstein equation, where we define the jump in f at a point y as By using the ansatz (2.5), Eq. (2.6) leads to Here, we have introduced ∆k 21 ≡ k 2 − k 1 .It was noted in ref. [32] that ∆k 21 > 0 is required for radion stabilization.Solving the Einstein equation for the subregions 1, 2, we obtain (2.9) Eqs. (2.8), (2.9) relate the brane tensions to the bulk cosmological constants.Without any radion stabilization mechanism, they require two additional fine-tuning conditions aside from the usual tuning for a vanishing four-dimensional cosmological constant.Once radions are stabilized, these extra fine-tuning conditions are removed. 2For a detailed discussion on the stabilization of radions based on the GW mechanism, see ref. [32].
Here, we just assume potentials for the radions generated as a result of some radion stabilization mechanism.The potentials stabilize the compactification radius to b = b 0 .

Getting back the ordinary Universe
In section 2.1, we have summarized the static solution for the three 3-brane construction.
In this limit, the background metric was given by Eq. (2.5).Conditions in Eqs.(2.8), (2.9) are required to be satisfied to maintain this static solution.However, one generally expects a deviation from this static solution in the presence of matter, radiation and vacuum energy.A new solution for the background metric becomes time dependent, and it is not guaranteed that one can recover the standard FLRW Universe and preserve the success of the standard cosmology.To answer this question, we will perform a series expansion of the metric, where Eq. (2.5) with Eq. (2.8), (2.9) is considered as the zero-th order solution.
Let us introduce matter for each brane and consider the energy-momentum tensor on the branes given by Eq. (2.4).We denote the energy densities and pressures on the UV, intermediate and IR branes as ρ * ,I,1/2 , and p * ,I,1/2 , respectively.In the presence of these perturbations, the jump conditions in Eq. (2.6) are altered.In particular, the energy densities ρ * ,I,1/2 modify the jump conditions for a(t, y), and the pressures p * ,I,1/2 appear in the jump conditions for n(t, y).Without a potential to stabilize radions, the system of the Einstein equation together with these jump conditions appears to be over-constrained.As discussed in ref. [36] for the two 3-brane setup, this is the result of demanding the static solution of Eq. (2.5) with Eq. (2.8), (2.9), even in the absence of the radion potential.By taking account of the radion background configurations and the solution of the '44' component of the Einstein equation, we can demonstrate that the effect of the matter perturbations is to shift the intervals between the branes to their new equilibrium positions, while we will discuss that the shifts are sufficiently small.With the new intervals, one looks for a solution for the other components of the Einstein equation.Our goal is to establish a perturbative expansion of the metric with the matter perturbations and study its impact on the late cosmology.We proceed with the following ansatz for the metric tensor: where A(t) is a function of t and (2.11) Note that the time dependence is generated as a result of the matter perturbations.At the leading order, the '04' component of the Einstein equation is solved trivially under the ansatz.Furthermore, the effect of δn(y) does not show up at this order, and hence ρ * ,I,1/2 determine the solution of the remaining components of the Einstein equation.
Our ansatz of Eqs.(2.10), (2.11) warrants some further clarification.As discussed above, initially, we need to consider shifts of the intervals between the branes.For the lighter radion, its potential near the minimum scales as U (b) ∼ m 5 radion (δb/b 0 ) 2 where m radion denotes the radion mass and δb is the shift.Hence, δb scales inversely with the radion mass.The cosmological implication of δb has been considered in detail in refs.[36,37].They concluded that the impact of δb is negligibly small.Concretely, analyzing the radion effective action, ref. [36] obtained δb/b 0 ∝ ρ NR /(m 2 radion M 2 IR ) where M IR denotes the typical scale of the IR brane and ρ NR is the energy density of nonrelativistic matter.For T 10 MeV during the BBN and the TeV scale radion, |δb/b 0 | 10 −20 .Therefore, we have assumed that the radion potential essentially fixes b to its stable value b 0 .The similar discussion holds for the heavier radion associated with the location of the intermediate 3-brane.The ansatz includes δa which determines the deviation of the metric from the standard RS type solution without the perturbations.A possible time dependence of δa can change the Hubble rate H → H + δ ȧ H(1 + δa), because the typical time derivative scale is determined by H.If |δa| 1, we can drop this correction term at the leading order and justify our ansatz consistent with the standard expansion of the Universe.
The Einstein equation for the bulk subregion p (= 1, 2) gives where H = Ȧ/A and Λ 5 p = − 6k 2 p κ 2 from Eq. (2.9).Plugging the ansatz of Eqs.(2.10), (2.11) into Eq.(2.12), we obtain the equation for δa at the leading order, for the subregion p.The solution for δa(y) for each subregion is given by where δa p (y) represents the solution for the subregion p and c 1 , c 2 , c 2 are constants to be determined.We have also used δa 1 (0) = 0.The constant c 2 can be expressed in terms of c 1 , c 2 , H by requiring that δa(y) be continuous across the intermediate 3-brane, namely δa 1 (y I ) = δa 2 (y I ).To determine c 1 , c 2 , H, we consider the jump conditions in Eq. (2.6), including the matter perturbations, Using the ansatz with the relation (2.8), we obtain the conditions for δa (y), Inserting Eq. (2.14) into Eq.(2.16), we obtain where we have defined the appropriate red-shift factors, Ω I ≡ e −k 1 b 0 y I for the intermediate brane and )b 0 y I for the IR brane.Solving Eq. (2.17), the Hubble rate is determined as Here, the 4-dimensional Planck mass is given by [32] Eq. (2.18) is the usual Hubble rate equation, where in the right hand side the physical energy density of the system determines the Hubble rate. 3As expected, the energy density at each brane is warped down by the appropriate exponential suppression factor.The constants in Eq. (2.14) are evaluated to be For the consistency of the calculation, we can check that the existing result for the two 3-brane setup [36] is correctly recovered in the limit where the intermediate 3-brane is absent, namely The deviation of the metric from the static solution is parametrized by δa in Eq. (2.10).Assuming no fine-tuning among the different perturbations, |δa(y)| takes its maximum value at y = 1/2, and hence we demand |δa(1/2)| 1 to recover the standard Hubble expansion equation.Using Eqs.(2.14), (2.18)-(2.22),we can write down the expression for the deviation of the metric at the IR brane, which is well-approximated by expanding in powers of the exponential suppression factors Ω 0∆ , Ω I , Ω 0∆ /Ω I , Here, we have assumed that k 1 , k 2 are not hierarchically different, in particular, Without assuming fine-tuning among the different perturbations, we require each term in Eq. (2.23) to be much smaller than unity, where we have defined the typical intermediate mass scale For the stability of the three 3-brane setup [32], V I should be positive, and hence one requires 0 < r 12 ≤ 1. Evidently, the constraint on the perturbation at the UV brane is ρ * (M IR M pl ) 2 .In a similar manner, from Eq. (2.26), we require the perturbation at the IR brane to satisfy ρ 1/2 k 2 1 M 2 pl /r 12 .When all the SM fields live on the IR brane, the condition (2.26) is satisfied for temperatures below the typical scale of the IR brane, T < M IR .Furthermore, the conditions on ρ * , ρ I in Eqs.(2.24), (2.25) are satisfied if the total energy density is dominated by the contribution from the thermal bath of the SM sector, and ρ * , ρ I give some small portions of the total energy density.Therefore, we can conclude that the late cosmology at temperatures lower than M IR follows the standard evolution.As we will discuss in the next section, the system at temperatures higher than M IR is described by a different spacetime geometry so that the current discussion is not applicable.

AdS-S and phase transitions
Let us now consider the system at high temperatures where the time is Euclidean and compactified on a circle.As discussed for the RS model with two 3-branes [38], our current system at high enough temperatures (above M I ) is described by an AdS-S spacetime with an event horizon.The AdS-S Euclidean metric is given by where L = 1/k 1 and ρ ≥ ρ h with the event horizon located at ρ = ρ h .This metric gives a solution to the Einstein equation when the time periodicity where T h is the Hawking temperature of the blackhole, otherwise we have a conical singularity at the event horizon.The dual 4D picture of the system with the real time is described by a thermal CFT with the temperature T h .The system at low enough temperatures (below M IR ) is described by the RS spacetime discussed in the previous section.For the current purpose, we consider the following form of the metric including two radion degrees of freedom T 1,2 (x) [32]: with the flat 4D metric η µν .This 5D metric is equivalent to the background metric in Eq. (2.1) when the radions are stabilized at the minimum.Explicitly, they are related through φ = b 0 y/ T 1 , for 0 ≤ y ≤ y I , and φ = b 0 y/ T 2 , for y I ≤ y ≤ 1/2.The two 4D radion fields are defined as The stabilization of all the 3-branes amounts to generating a 4D effective potential for the radions with a non-trivial minimum at µ 1,min and µ 2,min , which give the typical Note that the first term requires k 2 > k 1 .In the limit of k 2 → k 1 , this term becomes singular, which signals the disappearance of the intermediate brane.Generally, the mass of the radion µ 1 is around the energy scale of the intermediate brane and much heavier than that of µ 2 determined by the scale of the IR brane.
Since the system behaves differently at high and low temperatures, phase transitions take place between different regimes as the Universe cools down.In the current three 3-brane setup, we can consider three regimes (see Figure 1): (i) the AdS-S spacetime with an event horizon, (ii) the AdS-S spacetime with the intermediate 3-brane and (iii) the time-compactified RS spacetime with the intermediate and IR branes.Under the assumption of a large distance separation between the intermediate and IR branes, µ 2,min µ 1,min , we expect two phase transitions, i.e. the first transition from the regime (i) to (ii) and the second transition from the regime (ii) to (iii).In the dual 4D picture, they correspond to two confinement-deconfinement phase transitions.As the Universe cools down, the thermal CFT undergoes a confinement at the first phase transition and flows into another thermal CFT, which is in the end confined at the second transition.The mixing between two dilatons is suppressed by ∼ (µ 2,min /µ 1,min ) 3 with a conformal symmetry breaking parameter [19].Therefore, in the 5D picture, it is justifiable to assume that two radions are decoupled when µ 2,min µ 1,min .Let us first consider the phase transition from the regime (i) to (ii).As in the case of the two 3-brane model [38], we can expect that the phase transition proceeds through the pure AdS spacetime with the periodic time.That is, moving from the outside to the center of a nucleating bubble, we observe that the horizon goes toward ρ = 0 and then the intermediate brane comes in from ρ = 0 to a value at the minimum.To evaluate the critical temperature of the phase transition, we need to find the free energy of the system in each regime.For the regime (i), the free energy is dominated by the black hole solution (3.1), and we assume that the contribution from a radion potential is subdominant because of its negligible back-reaction to the background geometry.Taking account of the conical singularity effect, the difference between the free energy of the AdS-S spacetime and that of the pure AdS is calculated as [38] ∆F where we have defined which roughly give the numbers of colors in the two corresponding dual CFT theories.The free energy of Eq. (3.5), as a function of T h , has the minimum at T h = T .For the regime (ii), we assume that the critical temperature of the phase transition is somewhat smaller than the typical scale of the intermediate brane so that the heavier radion µ 1 is no longer in the thermal bath.The contribution to the free energy is then approximated by a potential stabilizing the heavier radion.Another contribution may come from the second thermal CFT with the number of colors N 2 in the dual 4D picture.This effect is, however, neglected if its free energy is smaller than that of the first thermal CFT with the number of colors N 1 , i.e.N 1 > N 2 .The intermediate brane is assumed not to contain a lot of light degrees of freedom so that their contribution to the free energy can be safely neglected.In the end, the free energy for the regime (ii) subtracted by that of the pure AdS is given by Here, V r,eff,1 (µ 1 ) denotes a potential stabilizing the heavier radion.We will discuss two different types of the potential below.At the critical temperature T c,1 , the free energy at its minimum for the regime (i) is equal to that for the regime (ii).By using Eqs.(3.5), (3.7), we find which is somewhat suppressed by a large N 1 factor and justifies our assumption.
Under the assumption of a sufficiently large distance separation between the intermediate and IR branes, we can discuss the phase transition from the regime (ii) to (iii) similarly.That is, the relevant part of the free energy for the regime (ii) is now given by the minimum of Eq. (3.5) with N 1 replaced by N 2 , and the free energy for the regime (iii) comes from a potential stabilizing the lighter radion V r,eff,2 (µ 2 ).The critical temperature T c,2 is estimated accordingly.

Goldberger-Wise potential for radion stabilization
To discuss the phase transitions more concretely, the radion potential needs to be specified.In this section, we consider a single 5D scalar field φ with a bulk mass m φ and localized brane potentials [32] as in the case of the original GW mechanism for the two 3-brane setup [33].The lighter radion potential is similar to that of the two 3-brane case, so we mainly focus on the potential for the heavier radion µ 1 .The low-energy 4D effective potential is obtained after integrating out the extra dimension (see Eq. (A.8) in ref. [32]), Here, the boundary potentials for φ(y) are 2 with dimensionless v UV , v I , a constant term is subtracted from the radion potential and O( 2 ) pieces are dropped.In addition, we have defined , which is a small parameter and ≥ −2 [60].We will take it to be positive for concreteness.Simplifying Eq. (3.9), we arrive at where we have introduced the ratio R v ≡ v I /v UV .The last term ∝ − δµ4 1 v 2 I represents a deviation from the static RS solution, for instance, due to loop corrections. 4Given that δ lies in the interval, V GW has a global minimum (and a maximum) at Note that X min > X max , so that the global minimum of the potential is separated from the origin by a potential barrier situated at µ 1,max .The potential V GW (µ 1 ) has a form µ 4 1 P (µ 1 ) with a polynomial function P , and hence represents an almost marginal operator of conformal dimension 4 + in the dual 4D picture.Therefore, for a small , the conformal symmetry is broken through the slow renormalization group (RG) evolution, so that a hierarchically small mass scale µ 1,min is generated.When the vacuum expectation value (VEV) of φ is large, its backreaction can deform the geometry and the current analysis is not reliable.To avoid this, we impose [44] v To fix µ 1,min in Eq. (3.12) for a given , the ratio R v has to be kept constant, which therefore implies Moreover, in order for the perturbative description of our 5D theory to be reliable, we assume N 1,2 4.4, which ensures that higher powers of Ricci scalars from the quantum gravity effect can be safely neglected [61]. 1. × 10 -29 Figure 2: The effect of δ on the GW potential for the heavier radion.We take δ = 0, 0.5 in the left and right panels, respectively.The bottom panels magnify the regions abound the origin of µ 1 in the top panels.For δ ∈ I δ and δ > , the minimum of the potential gets deeper and the potential barrier is significantly diminished.
The parameter δ makes the minimum of the potential deeper and reduces the potential barrier significantly.Concretely, we perform a series expansion in terms of and δ, and the radion potential (3.10) can be written as For δ , at the leading order, we obtain whereas, for δ = 0, at the leading order, Hence, as long as δ > 2 √ , the minimum of the potential gets deeper, as we show in figure 2. Note that the potential remains shallow because V GW (µ 1 ) is still nearly conformal invariant as a scale only enters as (µ 1 /k 1 ) .

Does the phase transition get completed?
Now, let us discuss the phase transition from the regime (i) to (ii), assuming the GW potential V GW (µ 1 ).In order for the phase transition to be completed, during the phase transition, the bubble nucleation rate should be greater than the Hubble rate per unit horizon time and volume, where the action S is dominated by O(4) symmetric bubbles for temperatures smaller than µ 1,min , and Γ 0 denotes the relevant energy scale of the potential raised to the fourth power.For our purpose of getting an estimate, the zero temperature action for the O(4) symmetric bubbles suffices, Here, µ 1,r is the field value which the field tunnels to and determined by minimizing S 4 .This can be done numerically, but for our rough estimation purpose, we can take µ 1,r ∼ µ 1,min as in ref. [44].The Hubble rate during the phase transition is H µ 2 1,min /M pl , which is valid for temperatures near or below the critical point.Note that the larger the intermediate brane mass scale is, the larger the Hubble rate is, and the bubble nucleation rate has to compete with this larger expansion rate for the phase transition to proceed.Therefore, in order for the phase transition to be completed, we require where we assume Γ 0 µ 4 1,min .It is evident that there is a tug-of-war between the back reaction condition in Eq. (3.13), and Eq.(3.20), namely that in order to satisfy the back reaction condition, we prefer a large N 1 , but the bubble nucleation action scales as N 4  1 , which should be small in order for the phase transition to proceed.In addition, −V GW (µ 1,min )/µ 4 1,min is preferred to be large, by taking , δ as large as possible, to make the potential deep.
Although δ helps deepening the potential, it is bound to live in the interval I δ .Moreover, one can not make arbitrarily large, keeping µ 1,min fixed, as v UV increases rapidly and eventually goes beyond the v UV N 1 regime.Now, looking at the right hand side of the inequality (3.20), we see that increasing µ 1,min makes it difficult for the phase transition to be completed.Therefore, without the help of other sources of conformal breaking, the transition from the regime (i) to (ii) does not get completed.The numerical calculation has also confirmed this conclusion.

Radion stabilization with dark gauge fields
In this section, we consider an alternative mechanism to stabilize the radions in the three 3-brane model to make the phase transitions completed correctly.Ref. [51] has presented a radion stabilization mechanism for the two 3-brane model with a bulk confining gauge field, where the deviation of the IR brane tension from the static condition contributes to a quartic potential for the radion, which is balanced by a potential generated from the gauge field confinement.We here generalize this mechanism for our setup and discuss the phase transitions with the new radion potentials.

The model
Let us introduce two SU (N H,i ) gauge fields residing respectively in the subregions i = 1, 2, as described in figure 3. Their 5D action with the metric defined in Eq. (3.2) is given by where F M N,i and g 5,i denote the 5D gauge field strength and coupling for the SU (N H,i ).
In addition, we consider deviations of the brane tensions from their static tuned values in Eq. (2.8), namely Their resulting contributions to the action are Here, g * ,I,1/2 are the induced metrics on the UV, intermediate and IR branes, respectively.We now perform the Kaluza-Klein (KK) decomposition and integrate out the extra dimension to obtain the 4D effective action for the zero-mode gauge fields.From the subregion 1, we find where µν,1 denotes the field strength of the zero-mode of the SU (N H,1 ) gauge field.Taking account of the 1-loop RG flow from the scale k 1 to a scale Q µ 1 , the 4D gauge coupling of the zero-mode gauge field for the subregion 1 becomes Here, b YM,1 = 11 3 N H,1 for the pure Yang-Mills theory.In the dual 4D picture, the first term in the right hand side can be understood as the contribution from the CFT degrees of freedom, which vanishes below the scale µ 1 .Thus, one can define ; A natural expectation from the CFT is b CFT,1 −ξ 1 N 1 with a positive constant ξ 1 .The confinement scale is given by the scale where the 4D gauge coupling blows up, where Λ H,1 denotes the confinement scale for µ 1 = µ 1,min .This equation is valid only when the confinement scale Λ H,1 (µ 1 ) is somewhat smaller than the radion field value µ 1 .
We parametrize our ignorance of the threshold between the confinement and deconfinement phases with γ c,1 : Here, the scale µ c,1 is defined such that Eq. (4.6) is valid for µ 1 ≥ µ c,1 .As the 4D effective description breaks down for the confinement scale larger than the mass of the lightest KK mode, one expects γ c,1 π.For µ 1 < µ c,1 , on the other hand, we expect the confinement scale to become independent of the radion field value, Equating Eq. (4.6) with Eq. ( 4.8) at the scale µ c,1 , we find . (4.9) The dark gluon condensates contribute to the heavier radion potential, where the non-vanishing trace of the stress-energy tensor T µ 1,µ appears as a result of the conformal anomaly.Therefore, the total effective potential for the heavier radion is given by ; for µ 1 > µ c,1 ; for µ 1 < µ c,1 . (4.11) Here, λ 1 ≡ 4δV I /k 4 1 , and this quartic coupling of µ 1 can be balanced against the term generated from the condensation, leading to a non-trivial minimum for the heavier radion.For n 1 < 1, the radion potential has a global minimum, 12) The constant δV * in Eq. (4.11) can be chosen to reproduce the observed vanishingly-small cosmological constant.
As the heavier radion µ 1 has been stabilized at µ 1 = µ 1,min , we look at the evolution of the theory from there on.At the second stage of the evolution, the SU (N H,2 ) gauge field in the subregion 2 confines, and similarly, we obtain the 4D effective action for the zero-mode gauge field, From the definition in Eq. (3.3), we have the following relation: Hence, the 4D gauge coupling of the zero-mode of the SU (N H,2 ) gauge field at a scale with b YM,2 = (11/3)N H,2 .Note that the heavier radion µ 1 has been already stabilized, and g 4,2 does not depend on µ 1 dynamically but just on µ 1,min .The confinement scale for the SU (N H,2 ) is then given by where b CFT,2 ≡ − 8π 2 k 2 g 2 5,2 ; As in the case of the SU (N H,1 ), the dark gluon condensates lead to the lighter radion potential.In addition, the deviation of the IR brane tension δV 1/2 generates an effective potential of the form -λ 2 4 µ 4 2 with λ 2 = 4δV 1/2 /k 4 2 for the lighter radion.Therefore, the similar discussion for the stabilization of µ 2 follows as long as n 2 < 1.

Phase transitions
As discussed in section 3, our three 3-brane setup at finite temperature can have three regimes: (i) the AdS-S spacetime with an event horizon, (ii) the AdS-S spacetime with the intermediate 3-brane and (iii) the time-compactified RS spacetime with the intermediate and IR branes.Then, we expect the phase transition from the regime (i) to (ii) and the transition from the regime (ii) to (iii), as described in figure 1, and need to investigate whether those phase transitions are correctly completed in the current model of radion stabilization with two dark gauge fields.
Let us first consider the phase transition from the regime (i) to (ii).For T T c,1 , the dominant contribution to the O(4)-symmetric bounce action comes from the gradient energy and is reliably estimated by the thick-wall approximation [62], where the critical temperature T c,1 in Eq. (3.8) is given by The tunneling point µ t,1 is found by minimizing the action, ∂S ∂µ t,1 = 0, and for a low T , it is obtained as If the phase transition is correctly completed, the nucleation temperature, T n,1 , is given by the solution of the equation, with the Hubble rate H µ 2 1,min /M pl .The left panel of figure 4 shows the parameter region that the bounce action is small enough to compete with the Hubble rate and hence the phase transition from the regime (i) to (ii) is safely completed.We can see that there exists a large viable region in the current model, unlike the case of the GW mechanism.The underlying reason for this result can be understood as follows.The radion potential generated by the GW mechanism is shallow as it is suppressed by in Eqs.(3.16), (3.17), and 1 in order to produce a hierarchically smaller mass scale than the Planck scale.Furthermore, one cannot make the VEVs of the GW scalar on the branes arbitrarily large without violating the backreaction condition in Eq. (3.13).In contrast, from Eq. (4.12), it is evident that the radion potential generated via the dark gauge field confinement can be much deeper, and hence can allow for a smaller bubble nucleation action.In terms of the dual CFT picture, the breaking of the scale invariance by a small parameter results into a much shallower dilaton potential for the GW-type stabilization, whereas the dark gauge field confinement breaks the scale invariance strongly at low energies.The backreaction problem is circumvented as well, thanks to the asymptotically-free nature of the dark gauge field.In each panel, the gray shaded region shows the excluded region that the bubble nucleation action is too large to compete with the Hubble rate.A cautionary remark is that the boundary of the gray shaded region is only an analytical estimate of the zero-temperature action.The blue shaded regions correspond to N Hi > N i , where the finite temperature effects of dark gauge fields become important (this is an approximate estimation, and can be modified by O(1) correction coefficient).The orange regions, corresponding to n i > 0.8, have µ c,i smaller than the QCD scale, and its effect has to be included.The left most violet exclusion regions correspond to the regime where quantum gravity effects become important.The red circle and violet triangle in each panel denote benchmark points for the discussion of the gravitational wave production (see figure 5).In the left and right panels, the black dashed contours depict the heavier and lighter radion masses, respectively.The remaining parameters are chosen as λ 1,2 = 1, ξ 1,2 = 0.5, and γ c1,2 = π.
After the first phase transition is completed, the intermediate brane appears and the Universe gets reheated.Assuming that the vacuum energy in the false vacuum is transformed into the radiation in the true vacuum and using the definition of the critical temperature in Eq. (3.8), the reheating temperature is estimated as where g * (T RH,1 ) denotes the number of degrees of freedom at T RH,1 .As the reheating temperature is above the energy scale of the IR brane, we enter into the phase where the IR brane is replaced by a horizon.The minimum of the free energy in this phase increases as the Universe is cooled, and eventually, becomes comparable to the potential minimum for the lighter radion, which defines as before the critical temperature T c,2 for the second phase transition from the regime (ii) to (iii).The similar argument as in the case of the phase transition from the regime (i) to (ii) can be applied to the second transition.In the right panel of figure 4, we show the parameter region that the bounce action is small enough to compete with the Hubble rate and the phase transition from the regime (ii) to (iii) is completed.Before closing this section, let us mention about the possibility for the production of dark glueballs, which are formed below the confinement scales of the dark gauge fields.However, the produced glueballs are unstable and immediately decay into lighter particles in the thermal bath.In more detail, let us first consider dark glueballs from the subregion 1.They, for example, couple to (zero-mode) dark gluons in the subregion 2 by where M I ∼ k 1 e −k 1 y I , Φ 1 denotes the dark glueballs from the subregion 1 and the coupling constant has been omitted.With M I Λ H,1 , the glueballs Φ 1 decay into lighter particles immediately and do not remain in the Universe.Similarly, dark glueballs from the subregion 2 couple to the SM particles, e.g. through interactions between the dark glueballs and the photons, and decay quickly.

Gravitational wave generation
The first order phase transitions in our multi-brane model produce stochastic gravitational wave backgrounds.We now assume the radion stabilization mechanism with dark gauge fields, discussed in section 4, and estimate their energy densities.Following ref. [63], a strong first-order phase transition can be described in terms of two quantities: the latent heat released α, which is defined as the ratio of the vacuum energy density released ρ vacuum to the radiation density ρ rad at the bubble nucleation temperature T n , namely Here, T RH is the reheating temperature, and t * denotes the time when gravitational waves are produced.We assume that t * is approximated by the time corresponding to the nucleation temperature.One feature of the current radion stabilization mechanism is that the nucleation temperature is not much smaller than the critical temperature.This means that a supercooling phase between the critical and nucleation temperatures is short compared to the corresponding one from the Goldberger-Wise stabilized model, so that gravitational waves generated from the first transition from the regime (i) to (ii) are not completely diluted away during the second transition from the regime (ii) to (iii) 5 .Since the second phase transition also produces gravitational waves, the total gravitational wave density is estimated as where g * ,s denotes the effective number of degrees of freedom in entropy, the two terms in the right hand side come from the two consecutive phase transitions, and we have assumed instantaneous reheating.Clearly, if we had T n,2 T RH,2 , the gravitational wave density generated from the first transition would be significantly suppressed.The gravitational wave density produced by each of the two phase transitions can be decomposed into contributions from the bubble collision, sound wave and turbulence of the thermal plasma, namely Ω GW,i h 2 = Ω col,i h 2 + Ω sw,i h 2 + Ω turb,i h 2 ; i = 1, 2 . (5. 3) The background of the strongly interacting SU (N H,i ) gauge field provides the thermal plasma that results in a terminal velocity for the expanding bubble wall [51,63].Therefore, the dominant contribution comes from the bulk motion of the fluid, in the form of the sound wave and turbulence.We use the numerical approximation and modeling presented in ref. [63].The contribution of the sound wave is then given by where i = 1, 2. For our purposes, we set the efficiency factor κ sw 1 and the bubble wall velocity v w 1.Note that as α i 1, its dependence drops out in Eq. (5.4).f sw denotes the peak frequency, redshifted accordingly.The contribution of the turbulence is (5.5) We take κ turb 0.05κ sw , conservatively.For more details on these analytical approximations, we point the reader to ref. [63] and references therein.
Figure 5 shows the gravitational wave spectrum with unique two peaks predicted in our model.The peak towards the higher characteristic frequency comes from the first phase transition from the regime (i) to (ii).In table 1 we show the relevant temperature scales and gravitational wave parameters for the chosen benchmark points.Note that this peak is damped down due to the supercooling during the second phase transition, while this is not washed away hopelessly beyond the sensitivity limit of future gravitational wave probes such as LISA, DECIGO, and BBO.If the intermediate and IR brane mass scales lie at relatively lower scales, then the associated reheating temperatures become smaller, and with appropriate choices of the parameters in the viable region, the two peaks can lie in the sensitivity range covered by LISA, as shown in figure 6.

Discussions
Warped extra dimensional models with multiple 3-branes can naturally accommodate hierarchically different mass scales that often appear in physics beyond the SM.We have contemplated the cosmological effects of such multi-brane warped models.In particular, it was confirmed that the usual FLRW Universe is recovered below the IR mass scale, and the Hubble expansion parameter is determined by the sum of brane-localized energy densities.We have also analyzed the first-order confinement-deconfinement phase transitions in the three 3-brane models and pointed out that the problem of the completion of the phase transitions gets worse, in comparison with the case of the two 3-brane model, with the usual GW stabilization mechanism.In hindsight, this is expected, as the bubble nucleation rate has to compete with the even larger Hubble rate corresponding to the intermediate mass scale.To circumvent the problem, we first developed an extension of dark gauge field stabilization of multiple branes.Here we utilized two dark gauge fields, which reside in the bulk subregions 1, 2, respectively.The confinement of these gauge fields provides the strong breaking of the scale invariance, and generates a potential for the radions deeper than that of the GW mechanism, which allows for the phase transitions to get completed.We have then analyzed the spectrum of gravitational waves generated by the first-order phase transitions.With the radion stabilization with two dark gauge fields, the supercooling epoch can be small enough to allow for the possibility of a two peak gravitational wave frequency spectrum whose amplitude is within the projected sensitivity limits of future space-based gravitational wave observers such as LISA, DECIGO and BBO for generic choices of the parameters.
In the present paper, we have assumed that the two phase transitions are decoupled from each other.This assumption is valid as long as N 1 > N 2 , and µ 1,min µ 2,min .In the dual 4D picture, it can be understood by the fact that the mixing between the two corresponding dilatons are proportional to (µ 2,min /µ 1,min ) 3 .It would be interesting to consider a scenario where the two radions are not completely decoupled from each other, and therefore, it is conceivable that the dynamics of the two phase transitions would be phenomenologically richer in terms of gravitational wave signatures.We leave it for a future study.
When the system is in the regime (ii), we considered the AdS-S spacetime with the Euclidean time compactified on a circle, and did not keep track of its time evolution as we have discussed in section 2 where the system is at a temperature below the scale of the IR brane.However, it may be important to understand the time evolution of the background spacetime in the regime (ii), if the system accommodates a topological object, such as a cosmic string, formed by some spontaneous symmetry breaking at the intermediate brane or in the subregion 1.The behavior of such a topological object depends on how the background spacetime evolves in time.Furthermore, it is interesting to note that a network of cosmic strings acts as a durable source of gravitational waves from the time of their production so that the resulting gravitational wave background stretches across a wide range of frequencies, which encodes the evolution of the background spacetime.The second phase transition from the regime (ii) to (iii) must also affect the evolution of cosmic strings and their production of gravitational waves.Since the confinement-deconfinement phase transitions themselves generate gravitational waves, the final spectrum is given by the mixture of those from cosmic strings and the phase transitions.We may discover the early Universe based on our multi-brane model by looking at the frequency spectrum of gravitational waves!refs.[34,35], it is unclear where the condition of E = 0 on the intermediate brane comes from.Ref. [34] simply takes it for granted that E = 0 condition is correct, and the author worries about an ambiguity due to a nonzero E , while ref.[32] does not encounter such an ambiguity.In more detail, around Eq. (6.4) in ref. [32], the authors get proper radion kinetic terms without any ambiguity for E on the intermediate brane.
Ref. [35] argues that if the radion is not stabilized, E = 0 on the intermediate brane can be taken.Following the discussion of ref. [32], however, this implies that another radion degree of freedom appears in the three 3-brane setup, at least, for the case without radion stabilization.This statement cannot be correct because the number of radion degrees of freedom should be determined regardless of their stabilization.In ref. [32], the condition of E = 0 is crucial to see multiple radion degrees of freedom, i.e. if we consider E = 0, we can only see one radon degree of freedom.Ref. [14] has also discussed that E on each brane is gauge independent, and a non-trivial E is obtained on the intermediate brane.
Let us consider what would happen if E = 0 were taken by hand on the intermediate brane.Then, we would get the same forms of the bulk equation and the boundary condition in every subregion as those for the two 3-brane case, i.e. the form of the bulk equation is the same as Eq.(6.3) of ref. [37] and the form of the boundary condition is the same as Eq.(6.4) of ref. [37].If we use the same forms of the bulk equation and the boundary condition in every subregion but take e.g.different bulk cosmological constants in different subregions (i.e.k 1 = k 2 ) different radion mass eigenvalues are obtained in different subregions, which contradicts the claim that there is one radion in the theory.This inconsistency does not appear in the calculation of ref. [32].Note that k 1 = k 2 corresponds to the zero tension case for the intermediate brane, i.e. the intermediate brane is absent, and the inequality k 2 > k 1 is required in order to avoid a tachyonic solution [14].Since E = 0 at the intermediate brane implies k 1 = k 2 , we conclude that it is an inconsistent condition.
Finally, let us describe our intuitive understanding of the existence of two radions in the three 3-brane setup.For the two 3-brane setup, the radion degree of freedom can be considered as the fluctuation of the IR brane compared to the UV brane (see e.g.appendix of ref. [12]).Now, if we introduce another brane (the intermediate brane), there is also the fluctuation of the intermediate brane compared to the UV brane, which

Figure 1 :
Figure 1: A schematic diagram of phase transitions via bubble nucleation in the three 3-brane model.There are three regimes: (i) the AdS-S spacetime with an event horizon, (ii) the AdS-S spacetime with the intermediate 3-brane and (iii) the time-compactified RS spacetime with the intermediate and IR branes.The first phase transition from the regime (i) to (ii) takes place through the nucleation of the intermediate brane bubbles.The second transition from the regime (ii) to (iii) takes place through the nucleation of the IR brane bubbles.

Figure 3 :
Figure 3: The schematic description of our setup to stabilize radions utilizing dark gauge fields.

Γ ≲ H 4 N 1 < 4 . 4 N 1 < N H 1 n 1 4 N 2 < 4 . 4 N 2 < N H 2 n 2 Figure 4 :
Figure4: The allowed (white) parameter regions that the phase transition from the regime (i) to (ii) (left) and the transition from (ii) to (iii) (right) get completed in the model with dark gauge fields.In each panel, the gray shaded region shows the excluded region that the bubble nucleation action is too large to compete with the Hubble rate.A cautionary remark is that the boundary of the gray shaded region is only an analytical estimate of the zero-temperature action.The blue shaded regions correspond to N Hi > N i , where the finite temperature effects of dark gauge fields become important (this is an approximate estimation, and can be modified by O(1) correction coefficient).The orange regions, corresponding to n i > 0.8, have µ c,i smaller than the QCD scale, and its effect has to be included.The left most violet exclusion regions correspond to the regime where quantum gravity effects become important.The red circle and violet triangle in each panel denote benchmark points for the discussion of the gravitational wave production (see figure5).In the left and right panels, the black dashed contours depict the heavier and lighter radion masses, respectively.The remaining parameters are chosen as λ 1,2 = 1, ξ 1,2 = 0.5, and γ c1,2 = π.