Moduli stabilization with bulk scalar in nested doubly warped braneworld model

We examine the modulus stabilization mechanism of a warped geometry model with nested warping. Such a model with multiple moduli is known to offer a possible resolution of the fermion mass hierarchy problem in the standard model. A six dimensional doubly warped braneworld model under consideration admits two distinct moduli, with the associated warp factors dynamically generating different physical mass scales on four 3-branes. In order to address the hierarchy problem related to the Higgs mass, both moduli need to be stabilized around their desired values without any extreme fine tuning of parameters. We show that it is possible to stabilize them simultaneously due to the appearence of an effective 4D moduli potential, which is generated by a single bulk scalar field having non-zero VEVs frozen on the 3-branes. We also discuss how the entire mechanism can possibly be understood from a purely gravitational point of view, with higher curvature f(R) contributions in the bulk automatically providing a scalar degree of freedom that can serve as the stabilizing field in the Einstein frame.


Introduction
The genesis of extra dimensional theories in physics dates back to the early 1920s, marked by Kaluza and Klein's original attempt to unify classical electrodynamics and Einstein gravity. In recent decades, higher-dimensional braneworld models have attracted fresh interest [1][2][3][4][5][6][7][8][9][10][11][12][13][14], arising as potential candidates for addressing the gauge hierarchy problem which is inherently linked to the problem of the Higgs mass. As the only fundamental scalar in the Standard Model (SM), the Higgs should experience large radiative corrections to its mass due to its self-coupling and couplings to other matter and gauge fields. These corrections are expected to flow up to the ultraviolet cut-off scale of the underlying quantum field theory. Protecting the Higgs mass from such Planck-scale corrections and keeping it safely within the TeV scale requires an extremely precise fine tuning of parameters, which leads to a naturalness problem within the SM. In this sense, the discovery of a Higgs boson as light as 125 GeV/c 2 [15][16][17] has simultaneously validated the last major prediction of the SM and exposed its central conundrum.
The warped braneworld model proposed by Randall and Sundrum [9] offers an elegant explanation of the large mass hierarchy, based on the premise of a non-factorizable warped fivedimensional spacetime geometry. The extra dimension is taken to be an orbifold with topology S 1 /Z 2 . Two 3-branes are located at the opposite boundaries of the bulk, which is a slice of AdS 5 . One of these branes (the "TeV brane") is identified with our visible universe, on which all the SM fields are assumed to be confined to. Gravity alone is assumed to permeate through the bulk. This places the apparent scale of gravity (M P l ) on the visible brane very close to the fundamental scale (M 5 ). In contrast, the energy density in the bulk induces an exponential warping along the extra dimension, which causes the physical mass m For kr c ∼ 11.54, the electroweak scale is dynamically generated from the fundamental scale. But the original model contains no mechanism to stabilize the magnitude of kr c around its desired value. This problem was addressed by Goldberger and Wise [18], who showed that kr c can be stabilized appropriately by introducing a massive scalar field φ in the bulk with quartic interactions localized on the branes. Integrating out the extra dimensions results in an effective potential V ef f (r c ) which stabilizes the modulus without any excessive fine tuning. This mechanism further allows one to interpret the modulus r c as the vacuum expectation value (VEV) of a dynamical radion field ρ(x µ ), i.e., r c = ρ , with the effective potential V ef f (ρ) driving ρ to settle at its minimum [19]. The mass of the radion and its coupling to matter fields on the visible brane are suggested to be O(TeV), rendering it detectable, in principle, at present generation colliders. Notably, the radion is expected to be lighter than the lowest-lying KK excitations of generic bulk fields, which are typically above O(TeV). This should make the radion the lightest detectable signature of the warped higher dimensional world.
There exists a wealth of work on the dynamics of various bulk fields and KK gravitons in warped spacetime, alongside a variety of generalizations with their own phenomenological and cosmological features (for a small body of examples, see [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]). The original Goldberger-Wise analysis has also been extended along several avenues, eg. by assuming finite values of the brane coupling constants [37], incorporating back reaction of the bulk scalar on the metric [38][39][40], and dynamical stabilization schemes in cosmological backdrops [41][42][43]. However, in spite of its theoretical appeal, the original 5D RS model has been increasingly challenged by experimental data in recent years. The mass of the first excitation of the graviton KK tower in the 5D setting [23] is suggested to be O(TeV). Moreover, the graviton KK modes are expected to couple to visible brane matter with amplified strength, as the lowest-lying excitations cluster close to the TeV brane. These considerations subsequently inspired several unsuccessful attempts to detect the signatures of the graviton modes through various channels at the LHC [44][45][46][47][48][49]. The absence of any such evidence till date has imposed serious constraints on the parameter space of the model. In particular, a small hierarchy m H /M 5 ∼ 10 −2 is required to explain the null results, indicating the potential appearance of new physics at least two orders of magnitude below the fundamental Planck scale. This is particularly problematic since this hierarchy requires r −1 c (which plays the role of the cut-off above which new physics is expected to appear) to be lowered by two orders of magnitude as well. But the sensitive dependence of the exponential warp factor on r c severely restricts this possibility, thereby limiting the efficacy of the five-dimensional setup in addressing the gauge hierarchy problem.

Review of doubly warped braneworld model
Several higher dimensional extensions of the original RS model have been proposed [50][51][52][53], with most of them introducing additional orbifolds with topology S 1 /Z 2 besides the first one. Some of these models, generalized to six dimensions, have interesting cosmological features stemming from dynamical stabilization of the extra spacelike dimensions [54]. A particularly interesting extension [55], capable of addressing the aforementioned issues, emerges in the form of a doubly warped braneworld with topology M(1, 3) × S 1 /Z 2 × S 1 /Z 2 . It results in a "brane-box" configuration, with four 4-branes forming the "walls" of the box and four 3-branes located at the "vertices" where adjacent 4-branes intersect. The bulk is a slice of AdS 6 with cosmological constant Λ 6 ∼ −M 6 , whereM is the fundamental Planck scale. The complete bulk-brane action (S), comprised of the bulk Einstein-Hilbert action (S 6 ) and the brane tension terms (S 5 ), is where R 6 is the six dimensional Ricci scalar, and y and z are angular coordinates charting the extra dimensions. For the current purpose, the visible 3-brane has been assumed to be devoid of any matter field, which would have otherwise introduced a further (S 4 ) term as the contribution of the matter Lagrangian. In a multiply warped setting, the brane tensions can, in general, be functions of the bulk coordinates. This is a departure from the 5D case, where the presence of only one extra dimension rules out any such dependence at the very outset. Using an RS-like ansatz to solve the resulting Einstein equations, one obtains the doubly warped metric where R y and r z are the radii of the orbifolds, and c and k are constants interconnected through This peculiar relation makes it clear that in absence of a considerably large hierarchy between R y and r z , both c and k cannot be simultaneously large. Instead, any one of the two distinct regimes c > k or c << k can emerge. It is not entirely possible to get rid of a little hierarchy between R y and r z , but with suitable choice of c and k, it can be ensured that R y /r z does not exceed O (10). The physical mass scale m In the c > k regime, the choice of c ∼ 10 and k ∼ 0.1 yields warp factors which together generate the TeV scale on the two 3-branes at {y = π, z = 0} and {y = π, z = π}, one of which can be identified with the visible brane. The other couple of 3-branes experiences negligible warping and remains at the Planck scale. This clustering of multiple 3-branes around each of the scales is a salient property of generalized RS models with nested warping, one having important phenomenological implications. Once the metric is fully solved, the brane tensions can be derived using appropriate junction conditions at the orbifold fixed points; V 1 and V 2 retain their z-dependence, whereas V 3 and V 4 turn out to be independent of y.
The coordinate dependence of the former pair can be shown to emerge from suitable scalar fields confined to the corresponding 4-branes. In the large k regime, one needs to account for a phantom scalar field on the y = π brane. This field may serve as a potential candidate for dark energy on the TeV branes. But the presence of such a phantom scalar renders the setup an effective field theory as opposed to a fundamental one. No such phantom field is needed in the large c regime. The doubly warped model has numerous phenomenological advantages over its progenitor model. Firstly, owing to the presence of two extra dimensions, the first excited graviton KK mode turns out to be considerably heavier than that of the 5D case. Moreover, coupling between graviton KK modes and SM fields on the visible brane is largely suppressed compared to the 5D model. Taken together, these features can satisfactorily explain the non-detection of KK gravitons at the LHC so far, without recourse to any small hierarchy between m H andM [56]. At the same time, significant portions of the parameter space of the extended model remain accessible to the LHC, allowing them to be explored in future runs [57]. Secondly, the doubly warped model can offer an explanation of the mass hierarchy among the SM fermions [55]. In the c > k regime, the O(TeV) 4-brane at y = π intersects two other 4-branes at z = 0 and z = π. Assuming SM fermions to be described by five dimensional fields confined to the y = π brane implies natural O(TeV) masses for the fermions. In addition to the z-dependent bulk wavefunction on this 4-brane, the fermionic fields can have kinetic terms on the 3-branes at the two intersection points. These boundary terms can modify the fermion-scalar Yukawa coupling on the {y = π, z = 0} and {y = π, z = π} 3-branes, thereby causing a splitting among the effective fermion masses. This splitting is arguably small, as the natural mass scales of these 3-branes are already clustered close to each other around the TeV scale.

Modulus stabilization in doubly warped model
Like the 5D RS model, the action of the 6D scenario contains no dynamics which can stabilize the extra dimensional moduli around their desired values. In absence of such an underlying stabilizing mechanism, the braneworld model alone cannot be considered adequate. Motivated by the success of the 5D Goldberger-Wise mechanism, it becomes natural to seek a similar approach for the 6D model that might stabilize both c and k (or equivalently, R y and r z ) simultaneously. While other phenomenological aspects of multiply warped spacetimes are well-studied [58][59][60][61], there has been little work in this direction so far, with the notable exception of [62]. The latter study proposes disjoint stabilization mechanisms for c and k with the help of two separate bulk and brane-localized fields in the c > k regime. In the other regime, taking c << k makes the metric almost conformally flat, wherefore only r z needs to be stabilized satisfactorily, with R y either left unstabilized (which is justified because of the negligibly small value of c) or stabilized with the help of another brane-localized field. Whether stabilization of both moduli can be achieved with the help of a single bulk scalar field, has, however, been an open question so far. In this paper, we attempt to address this very question.

Dynamics of the bulk field
For our current purpose, we choose to work in the c > k regime, which, as noted already, obviates the need to introduce a phantom scalar (which raises phenomenological difficulties) on the y = π brane to explain the coordinate dependence of its tension. Additionally, the relative smallness of k allows us to view the entire setup as a not-too-large departure from the 5D model, helping identify the key points of deviation from the latter clearly. A study of the complementary c << k domain using tools from supersymmetric quantum mechanics appears in [62].
Analogous to the original Goldberger-Wise scenario, we consider a bulk field propagating freely through the extra dimensional y − z bulk, and interacting only at the locations of the four 3-branes through quartic self-interaction terms. The action (S GW ) of the bulk field is given by an immediate generalization of the original 5D Goldberger-Wise action. Here,g 4 is the induced metric on the i th 3-brane, with u i being the VEV of the bulk field (with mass dimension [u i ] = +2) and λ i the corresponding coupling constant ([λ i ] = −4). With its dynamics confined to the higher dimensional space, the bulk field is essentially "frozen" on the corner branes. Defining a(y) = exp(−c|y|) and b(z) = sech(kπ)cosh(kz), we extremize this action with respect to φ to obtain the following equation of motion.
Away from the boundaries, the contributions of the interaction terms vanish. Assuming a separable solution of the form φ(y, z) = φ 1 (y)φ 2 (z), the bulk equation of motion can be reduced to the following pair of uncoupled ODEs, where −4α 2 is the separation constant.
An interesting physical interpretation of α emerges immediately. From the equation for φ 1 , it is apparent that 2α plays the role of the mass of the φ 1 field. Also, as the k → 0 (and simultaneously r z → 0) limit implies b(z) → 1 and mathematically reduces the 6D metric to the familiar singly warped form, it must also enforce 2α → m in order to reduce this equation to the 5D Goldberger-Wise bulk equation of motion. As for the equation for φ 2 , it is satisfied trivially in this limit, as φ 2 (z) → 1 leads merely to m 2 = 4α 2 . However, from a physical standpoint, one cannot of course allow k to be arbitrarily small, as any effective field theory in the semiclassical approach strictly remains valid only for r z ≥ M −1 , with quantum gravity effects dominating for r z < M −1 . In this sense, the k → 0 limit is not a physically tenable one, but comes with a lower cut-off regulated by Λ 6 . The bottom line is simply that for sufficiently small k, one can expect α to be reasonably close to m, which is a fact that proves useful in due course. Eq. (9) can be solved exactly for the component fields φ 1 (y) and φ 2 (z), leading to the following general solution for the bulk field under the condition of Z 2 orbifold symmetry.
φ(y, z) = e 2c|y| Ae νc|y| + Be −νc|y| DP l n (tanh(k|z|)) + EQ l n (tanh(k|z|)) sech where A, B, D and E are four arbitrary constants (arising on account of each equation being a second order ODE), and P l n and Q l n are associated Legendre functions of the first and second kind respectively. The quantities ν, n and l are defined as Assuming the magnitudes of the brane coupling constants to be large (i.e. λ i → ∞) allows the identification of simple energetically favourable configurations to serve as boundary conditions. The very structure of the interaction terms, given by as such a configuration on the corner brane at (y 0 , z 0 ). This approach leads to the following four equations, where the shorthand τ z = tanh(k|z|) has been introduced for convenience.
Eq. (14)−Eq. (17) are not all linearly independent, and as such, only allow three of the constants to be solved in terms of the remaining fourth. This poses no problem though, as the forms of the solutions still allow φ(y, z) to be determined uniquely. In the large c regime, the normalized solutions of the individual component fields are where c 2 = 1, and the other three constants are given by It is interesting to note that, apart from the overall normalization factor of u 1 , the dynamics of φ(y, z) is controlled not by the absolute values of the VEVs but only by their ratios. This feature is reminiscent of the 5D mechanism. Further, Eq. (18)−Eq. (22) dictate u 4 /u 2 = u 3 /u 1 . So the two ratios appearing explicitly in the solution are sufficient to completely specify all the six possible ratios among the four VEVs, thus ruling out any ambiguity.

The effective potential
Having obtained the solution of the bulk field, one needs to substitute it back in equation Eq. (7) and integrate over y and z in order to obtain the stabilizing potential V ef f (c, k), or equivalently, V ef f (R y , r z ). From the resultant effective 4D action S ef f , the definition S ef f = − d 4 xV ef f (c, k) allows the potential to be read off directly. In the large λ i limit, this amounts to evaluating only the bulk contribution. First, let us define the following dimensionless parameters.
These quantities roughly estimate the ratios between the mass parameters of the bulk field and the fundamental Planck scale. From the semiclassical standpoint, the admissible range of each ratio should be 0 < µ i < 1. In terms of these ratios, the parameters from Eq. (11)−Eq. (13) can be re-expressed as ν = 2 1 + µ 2 2 cosh 2 (kπ), n = ν − 0.5, and l = 5/2 1 + 4µ 2 1 /25. Upon substituting the solution from Eq. (10) in the action, the integral over y can be readily evaluated, but the presence of the special functions makes the z-integral analytically intractable. Defining the dimensionless potentialṼ ef f = V ef f /u 2 1 , and making use of Eq. (3) to eliminate R y and r z in favour of c and k, the following form emerges for the potential where the functions F 1 (k) and F 2 (k), arising out of the z-integral, are given explicitly by There is little choice but to evaluate F 1 (k) and F 2 (k) numerically for different values of µ 1 , µ 2 , and u 4 /u 2 , which, by this point, serve as three of the fundamental parameters of the model. The fourth parameter u 2 /u 1 , contained in the coefficient c 1 , contributes chiefly to the c-dependence of the potential. Equipped with the potential given by Eq. (24)−Eq. (26), we are in a position to demonstrate the existence of a simultaneous minimum ofṼ ef f in c and k over some region of the parameter space that doesn't require excessive fine tuning.

Stabilizing k
Owing to the explicit form of c 1 from Eq. (20), both the linear and quadratic terms in c 1 appearing in the coefficients of F 1 (k) and F 2 (k) in Eq. (24) are suppressed at least by O(e −4cπ ). For c ∼ 10, this suppression factor is nearly of order 10 −68 . So the leading order terms within both coefficients are the c 2 2 terms which suffer no such suppression. In order to study the k-dependence of the potential, it suffices to approximateṼ ef f by taking only the dominant contributions offered by these two terms. Recalling c 2 = 1, the potential effectively reduces to a single-variable function of k given bỹ Establishing the stabilization of k is thus tantamount to locating a suitable minimum of F (k). The plots in figure 1 show the behaviour of F (k) for various combinations of the parameters. For all the choices of parameters, F (k) admits a broad minimum at some k < 0.5. The location of this minimum is physically important, as the range 0.1 ≤ k ≤ 0.6 is of particular interest to scenarios attempting to explain the non-detectability of graviton KK modes at the LHC [56] [57], or exploring the phenomenology of off-brane SM fields that extend into the bulk [60] [61]. As special cases, it can be checked that µ 2 → 0 produces an asymptotically decaying F (k) with no finite minimum, whereas u 4 = u 2 admits only k = 0 as the global minimum. As expected, the most pronounced dependence of the minimum is on the VEV ratio u 4 /u 2 , which governs the potential at the leading order.

Stabilizing c
The previous plots reveal the existence of a fairly large parameter space which can stabilize k around its desired value. ButṼ ef f (c, k) also needs to stabilize c, which plays the leading role in determining the degree of warping. To this end, we need to reinstate the previously neglected c-dependence inṼ ef f (c, k) in order to locate a suitable minimum along c. Owing to the form of c 1 in Eq. (20), one might expect such a minimum to occur at c 1 = 0, which would be close in spirit to the 5D Goldberger-Wise result. Upon closer inspection though, it becomes clear that the existence of such a minimum depends crucially on the value of µ 2 . For considerably small values of µ 2 which render ν ≈ 2 + k , where k = µ 2 2 cosh 2 (kπ) is small enough so that O( 2 k ) and higher are negligible compared to unity, the deviation δṼ ef f for any arbitrarily small deviation c = c 0 − δ from the alleged minimum c 0 can be estimated as For excessively small k and vanishingly small δ, this quantity is negative. This clearly rules out a minimum at c 0 . The trouble can be traced to the existence of the O(c 1 ) terms iñ V ef f (c, k). While the 5D stabilization mechanism ensured automatic cancellation of the pair of O(c 1 ) terms in the potential, the current model offers no such way out by default. However, if some particular combination of parameters can result in such cancellation (without excessive fine tuning of course), the issue would be resolved. Once k has been stabilized at some specific k min , the condition for the O(c 1 ) terms to cancel each other translates to As F 1 (k) and F 2 (k) themselves depend on all three parameters in a rather complicated manner, this equation is transcendental and cannot be solved analytically. However, it can be checked numerically that for every choice of µ 1 , the allowed parameter space contains approximate solutions to Eq. (29), some of which are presented in Table 1.  Each of the combinations ensures that the non-negative O(c 2 1 ) terms dominate over the O(c 1 ) terms by nearly two orders of magnitude, leading to a proper minimum c min (satisfying c 1 = 0) with the following approximate form One can immediately estimate the corresponding u 1 /u 2 ratio which produces the desired value of c min ∼ 11.54. As given in Table 1, its magnitude depends strongly on µ 1 . In order to stabilize k between 0.1 and 0.6, the minimum admissible magnitude of u 1 /u 2 is O(10), which is slightly larger than the requirement in 5D. For larger values of µ 1 , a minor hierarchy between u 1 and u 2 becomes increasingly prominent. This hierarchy has its physical origin in the fact that we have essentially made the component fields φ 1 (y) and φ 2 (z) yield two distinct warping scales k ∼ 0.1 and c ∼ 10 in spite of having chosen their mass parameters (m and α) to be of comparable magnitudes! The somewhat large u 1 /u 2 ratio is precisely the price we need to pay for that. As a reality check, it is instructive to check that choosing u 1 ∼ u 2 in Eq. (30) would have resulted in c min ∼ 0.1, just as this argument suggests. The mutual cancellation of the O(c 1 ) terms involves a certain degree of fine tuning among the three parameters which control F 1 (k) and F 2 (k). While the situation is clearly more delicate than the 5D case, this tuning need not be very extreme. The O(c 2 1 ) terms can dominate and provide c with a stable minimum if the cancellation of the O(c 1 ) terms in Eq. (29) is accurate roughly up to O(10 −2 ), which is sufficient to suppress the latters' contribution by O(10 −2 ). For c ∼ 10, the critical values of the parameters need to be accurate at most up to O(10 −3 ), while larger values of c are somewhat more likely to ameliorate the situation due to increased suppression. This is no worse than the fine tuning associated with the choice of m/M for a given VEV ratio, that is required to generate the TeV scale accurately in the 5D case. Moreover, for some fixed ν(k), the tuning associated with u 2 /u 1 can be even more lenient due to the logarithmic dependence of c min on u 2 /u 1 . So simultaneous stabilization of c and k depends on a small degree of fine tuning among the parameters µ 1 , µ 2 , and u 4 /u 2 , which constitutes a crucial aspect of the extended Goldberger-Wise mechanism in the doubly warped scenario. This feature is not surprising since one should physically expect a two-level tuning for the stabilization of two distinct moduli in a spacetime with nested warping. As the stabilization of k alone requires negligible tuning (as evident from Fig. 1), the stabilization of c justifiably involves both levels. It might be interesting to investigate if incorporating back reaction or quartic self-interaction terms within the bulk can lead to improvements for this tuning requirement.

Phenomenological implications
The minor hierarchy between u 1 and u 2 (or equivalently between u 3 and u 4 ) can be better interpreted if we arrange the magnitudes of the four VEVs in proper order. As the boundary conditions require u 1 /u 2 = u 3 /u 4 , the VEVs obtained in each case satisfy The first two values are very close to each other, as are the last two, with the aforesaid hierarchy pushing the pairs apart. Interestingly, Eq. (31) reflects the order of physical mass scales on the corresponding corner branes (on which these classical values of φ are defined) for large c and small k, as can be checked using Eq. (4) . This can be understood physically as the massive bulk scalar, once frozen on the boundary branes, naturally tends to act against the warping induced by the bulk energy density. Consequently, the resulting mass scales on corner branes associated with larger VEVs can be generally expected to be higher than those with smaller VEVs. This feature is visible in the 5D Goldberger-Wise mechanism as well, albeit in a greatly tempered form. In the present case, it is more pronounced as it also incorporates the clustering effect observed among the physical mass scales. Due to the smallness of k, phenomenological features associated with c should be sufficiently similar to those of the 5D model. At first glance, this might make the smallness of u 2 compared to u 1 for larger values of µ 1 appear alarming, as the mass of the associated radion in the 5D model is proportional to the ratio between the classical value of φ on the TeV brane and the fundamental scale. In particular, one might worry that such small u 2 could bring the predicted radion mass down to the MeV scale or lower, thereby giving rise to phenomenological issues. It must be borne in mind, however, that we have neglected the back reaction of the bulk field. On occasions, this can cause underestimation of the radion mass by a few orders of magnitude [38] [39]. Furthermore, unlike the 5D model, one can exploit the 6D setup immensely by allowing SM fields to extend into the bulk. This possibility, which constitutes one of the most interesting features of the doubly warped model, might imply further corrections to the radion mass due to interactions of φ with these fields. There can also be significant mixing between φ and the scalar degrees of freedom which give rise to coordinate dependent 4-brane tensions. Any realistic study of radion phenomenology needs to address these questions in order to estimate the magnitude of the radion mass accurately.
Besides, due to the presence of two distinct moduli in the model, one should physically expect the appearance of two fundamental radions with different masses and widths. Owing to the smallness of k, the second radion should not suffer any significant suppression from the fundamental scale. This is echoed by the form of the stabilizing potentialṼ ef f (c, k) from Eq. (24), where F (k) is usually O(1) around its minimum. However, the aforementioned subtleties need to be accounted for in this case too before its phenomenology can be deduced conclusively. Due to the evidently complicated nature of such a study, we defer it to a future work.

Insight from higher curvature gravity
One of the most compelling advantages of having the size of the extra dimensions set by a single bulk field (as opposed to separate bulk and brane localized fields) is that it allows a purely gravitational interpretation of the stabilizing mechanism. Conventional wisdom suggests that the Einstein-Hilbert action, which provides an effective low energy description of gravity, needs to be amended with additional higher curvature terms respecting diffeomorphism invariance at sufficiently high energy scales. The warped geometry model which is being considered here has a large cosmological constant ( ∼ M P ) in the bulk and as a result the inclusion of higher curvature terms is a natural choice. Two broad classes of such higher curvature theories are the quasi-linear Lanczos-Lovelock models and f (R) models. While Lanczos-Lovelock models enjoy the benefit of being naturally ghost-free [63][64][65], the mathematically simpler f (R) models, equipped with specific conditions to ensure freedom from ghosts, pass some of the cosmological tests [66][67][68][69][70][71]. Furthermore, any given f (R) action typically admits a dual scalar-tensor representation [72][73][74][75]. In the so-called Einstein frame (related to the Jordan frame through a conformal transformation), the situation is equivalent to that of a scalar fieldφ coupled minimally to gravity, alongside a potential U (φ) whose form is determined by the functional form of f (R). For singularity-free metrics which are not experiencing rapid evolution, this equivalence holds physically [76][77][78][79][80].
In recent works [81] [82], it has been shown how such a scalar degree of freedom, arising solely from gravity in the 5D RS model, can play the role of the bulk field in the Goldberger-Wise scheme, thus obviating the need to introduce the latter by hand. By choosing f (R) = R + aR 2 − |b|R 4 (where a and b are coupling constants satisfying a > 0 and a > |b| to ensure freedom from ghosts), the potential can be arranged to contain both quadratic and quartic terms, with the latter encapsulating the effects of back reaction. Although the technique can be readily extended to spacetimes of arbitrary dimensionality, the inclusion of back reaction quickly renders multiply warped settings intractable in their full generality. In the following analysis, we discuss an extension of this technique to the doubly warped model, by retaining only the lowest-order correction term in f (R). To that end, we choose f (R) = R + aR 2 , which provides the action where κ 6 is the six-dimensional gravitational constant. In order to arrive at the Einstein frame, one conventionally makes a detour [75] through the intermediate Jordan frame representation where one introduces the auxiliary field χ and defines ψ = f (χ). For f (χ) = 0, the equation of motion for χ from Eq. (33), i.e. the on-shell condition, imposes χ = R. This makes Eq. (33) equivalent to the original f (R) action in Eq. (32). In order to reduce this to the minimally coupled Einstein frame representation, we apply the conformal transformation where the last equality follows from the on-shell condition. The action then transforms to Substituting f (R) = R + aR 2 in Eq. (34) and inverting the relation yields R(φ), which can be plugged immediately in Eq. (35) to obtain U (φ) as follows.
The minimum of U (φ) occurs atφ = 0, as evident from U (0) = 0 and U (0) > 0. Expanding U (φ) about this minimum, the leading order non-vanishing contribution comes from the quadratic term, with all subsequent terms increasingly suppressed by higher powers of κ 6 .
Having identified 1/5a with the mass squared of the scalar modeφ, the action from Eq. (35) reduces precisely to the sum of the Einstein-Hilbert action and the Goldberger-Wise action in the 6D bulk, withφ in the Einstein frame remarkably playing the role of the bulk field. While [82] uses this approach to deal with the full back-reacted problem (by including the R 4 term in the f (R) Lagrangian) for the 5D RS model in the Jordan and Einstein frames separately, such an exact treatment is too complicated to be feasible in the 6D case. Instead, with the results from the gravitational sector at hand, one can use heuristic arguments to bridge the gap with the existing bulk field method from the preceding sections. As the physical origin of the higher curvature correction(s) can be traced to the bulk energy density, it stands to reason that the resulting scalarφ should be an explicit function of the coordinates y and z only. With the metric solution in the Einstein frame given by Eq. (2), the equation of motion forφ is identical to Eq. (8) away from the boundaries. This preventsφ from having any non-trivial dynamics on the 3-branes, which is in agreement with the bulk field prescription. As an immediate corollary, we obtain four constant values ofφ (having mass dimension +2) serving as fixed boundary values on the four branes. The results are completely analogous to Eq. (14)−Eq. (17), with the only difference being that the current approach renders these boundary conditions exact, whereas in the earlier formulation they were valid only in the large coupling regime. With the full solution Eq. (18)−Eq. (22) forφ subsequently at hand, the rest of the analysis can proceed exactly as before, leading to the stabilization of both moduli under appropriate choices of parameters.

Discussions
The prospect of stabilizing both moduli of a doubly warped Randall-Sundrum braneworld model using a single bulk scalar field has been studied. Such a mechanism is crucial for a complete resolution of the gauge hierarchy problem in a higher dimensional scenario, and needs to supplement the gravitational part of the action giving rise to the warped metric. While the approach taken here is essentially a direct generalization of the Goldberger-Wise mechanism to six dimensions, the presence of nested warping brings out additional subtleties and constraints on the parameter space. As demonstrated, these constraints do not necessarily involve any extreme fine tuning of the fundamental parameters. In the c > k regime, which is phenomenologically preferred as it requires no brane-localized phantom field to explain the coordinate dependence of the 4-brane tensions, the effective potential admits a true minimum in k around k min ∼ 0.1 without any significant tuning. The stabilization of c, on the other hand, requires tuning on two different levels: firstly among the parameters µ 1 , µ 2 , and u 4 /u 2 (for the appearance of a suitable c min ), and secondly for u 1 /u 2 (to ensure c min ∼ 12). This can be interpreted physically as a consequence of the doubly warped structure of the underlying spacetime. A further departure from the singly warped model is that in order to achieve 0.1 ≤ k min ≤ 0.6 and c min ∼ 12, the minimum admissible magnitude of u 1 /u 2 is O(10), which is one order of magnitude higher than the analogous requirement in case of the 5D mechanism.
The bulk scalar approach is especially attractive as it allows room for a purely gravitational interpretation, with higher curvature contributions in the bulk automatically giving rise to the required scalar mode and its potential in the Einstein frame. Since the bulk of such warped geometry model is endowed with large bulk cosmological constant , the contributions from higher curvature terms become natural. This motivates us to include higher curvature terms in the bulk such as in f (R) model. This feature distinguishes it from certain other stabilization schemes, eg. stabilization of the two moduli with two distinct bulk and brane-localized fields. For appropriate choices of f (R), a variety of bulk scalar potentials can be generated, the simplest of which is the quadratic potential considered here. While the singly warped model could be solved exactly in presence of non-negligible back reaction, the doubly warped model becomes unsolvable as the field equations turn out to be non-linear PDEs. Nevertheless, the well-known R 2 correction is sufficient to make contact with the conventional approach on physical grounds. In soothe, at sufficiently high energy scales, gravity alone appears capable of both warping spacetime and determining the degree of warping. This possibility is intriguing, as it may produce observable TeV scale signatures of higher curvature gravity that can be explored in future collider experiments. In principle, foremost among them should be the O(TeV) radion mass associated with the larger modulus, which would be of fundamental importance in estimating/constraining the magnitude(s) of the higher curvature coupling(s).
The current study can be extended along various avenues. As pointed out earlier, a realistic study of radion phenomenology in a multiply warped background needs to include interactions of the bulk field with higher dimensional fermionic and gauge fields into account, alongside mixing with the brane-localized scalars responsible for coordinate dependent brane tensions. These effects, which constitute salient features of geometries with nested warping, may introduce a plethora of non-trivial modifications as far as collider signatures are concerned. Theoretically, it is worth investigating how other viable choices of f (R), or other classes of higher curvature theories (eg. Einstein-Gauss-Bonnet gravity), incorporate higher order phenomena like significant back reaction and affect the stabilization scheme. Such studies would necessarily have to rely on numerical techniques due to the complexity of the model. As plausible alternatives, one can also attempt to explain the origin of the stabilizing potential from quantum and/or thermodynamic perspectives, proceeding along the lines of [83][84][85][86][87]. Finally, it would be interesting to study the role of multiple dynamical radion fields in various cosmological contexts as well, eg. in the inflationary and bouncing settings, alongside their signature on the Cosmic Microwave Background.