Self-Interacting Gravitational Exciton Condensates from Extra-Dimensional Stabilization

1. Introduction

The discovery of the accelerated expansion of the Universe [1,2] and the overwhelming evidence for dark matter from observations such as the cosmic microwave background [3] and galaxy rotation curves [4] have profoundly reshaped modern cosmology. These observations indicate that over 95% of the energy content of the Universe is composed of dark energy and dark matter, phenomena that remain elusive within the Standard Model of particle physics.

Higher-dimensional theories—motivated by string theory, M–theory, and brane–world scenarios [5,6,7]—offer a compelling framework for unifying gravity with the other fundamental forces. In these models, the geometry and dynamics of compact extra dimensions are crucial in determining the effective four–dimensional physics. Dimensional reduction leads naturally to moduli fields, whose fluctuations yield massive scalar excitations (often referred to as gravitational excitons or radions) [8]. Their effective potential,

$$U_{\mathrm{eff}}\left(\phi \right)={\Lambda }_{\mathrm{eff}}+{\displaystyle \frac{1}{2}}{m}^2{\phi }^2+{\displaystyle \frac{\lambda }{4}}{\phi }^4+\cdots ,$$

(1)

emerges from a combination of classical geometric effects, flux stabilization mechanisms [9], and quantum corrections [10,11,12]. Due to their naturally suppressed couplings to Standard Model particles [13,14], these excitons are promising dark matter candidates.

In this work, we extend the conventional effective action by incorporating quartic self–interaction terms as well as higher–derivative corrections, thereby capturing the full nonlinear dynamics of gravitational excitons. This extension is particularly important since it permits the formation of a Bose–Einstein condensate (BEC) of gravitational excitons—a state that could radically alter our understanding of dark matter. The non–relativistic dynamics of such a self–gravitating condensate is effectively described by a Gross–Pitaevskii equation coupled to a Poisson equation [15,16,17,18]. Furthermore, renormalization group (RG) improvement techniques [19,20,21] are employed to incorporate quantum corrections that stabilize the moduli and generate non–trivial potential landscapes.

Our analysis spans several key aspects of the theory. First, by deriving the extended effective action through dimensional reduction and conformal rescaling, we explain the connection between extra–dimensional dynamics and four–dimensional phenomenology. Next, we incorporate quantum corrections and study the RG evolution of the effective couplings, thereby demonstrating how nonlinearities emerge beyond the quadratic approximation. The non–relativistic limit of the RG–improved theory is then taken to derive a Gross–Pitaevskii description of the gravitational exciton condensate. Stationary solutions are analyzed using both variational methods and linear stability techniques (including the Bogoliubov–de Gennes formalism), and the collective excitation spectrum is derived via Bogoliubov theory [22,23].

Notably, our findings indicate that the condensate exhibits a diffuse density profile and a unique excitation spectrum, with a crossover from linear (phonon–like) behavior at low momenta to quadratic free–particle dispersion at high momenta. These properties may lead to observable astrophysical signatures, such as modifications in cosmic microwave background anisotropies and distinctive dynamics in dwarf galaxies [24,25].

In the following sections, we detail the derivation of the effective action, the incorporation of quantum corrections, the formulation of the condensate dynamics, and the stability analysis of stationary solutions. Numerical simulations of the condensate, together with computations of its Bogoliubov excitation spectrum, confirm the theoretical predictions and highlight promising avenues for future observational tests.

2. Extended Effective Action with Higher-Derivative and Self-Interacting Terms

In this section, we develop an extended formulation of the effective action obtained from a higher-dimensional gravitational theory upon compactification. Our approach incorporates quantum corrections, non-trivial topological effects, flux stabilization, and higher-derivative interactions that naturally arise in string-inspired and extra-dimensional models. Our treatment builds upon and extends the methods described in [5,8,26] and further refines the framework discussed in [27,28].

1. Higher-Dimensional Action and Geometrical Setup

We begin with the D-dimensional Einstein–Hilbert action including a bulk cosmological constant and matter contributions,

$$S_D={\displaystyle \frac{1}{2{\kappa }_D^2}}{\int }_M{d}^{D}x\sqrt{\left|g\right|}\left\{R\left[g\right]-2\Lambda \right\}+S_{\mathrm{matter}},$$

(2)

where ${\kappa }_D$ is the gravitational coupling in D dimensions and $\Lambda$ is the cosmological constant. We assume that the underlying manifold factorizes as

$$M=M_0\times M_1\times \cdots \times M_n,$$

(3)

with the external spacetime $M_0$ being $(D_0\equiv d_0+1)$-dimensional (with coordinates $x^{\mu }$, $\mu =0,1,\cdots ,d_0$), and each internal space $M_i$ having dimension $d_i$. Furthermore, every internal metric $g_{mn}^{\left(i\right)}$ is assumed to satisfy the Einstein condition:

$$R_{mn}\left[g^{\left(i\right)}\right]={\lambda }^i{g}_{mn}^{\left(i\right)},m,n=1,\cdots ,d_i,$$

(4)

with constant eigenvalues ${\lambda }^i$. This condition ensures the consistency of the compactification procedure and has been widely employed (see, e.g., [29,30]).

2. Warped-Product Metric Ansatz and Moduli Fields

To capture the dynamics of the extra dimensions, we adopt a warped-product metric ansatz:

$$d{s}^2=g_{\mu \nu }^{\left(0\right)}\left(x\right)d{x}^{\mu }d{x}^{\nu }+\sum _{i=1}^n{e}^{2{\beta }^i\left(x\right)}{g}_{mn}^{\left(i\right)}\left(y_i\right)d{y}_i^{m}d{y}_i^n.$$

(5)

The warp factors $e^{{\beta }^i\left(x\right)}$ represent the dynamical sizes of the internal spaces and are interpreted as moduli fields. In order to study small fluctuations about a stabilized configuration, we define shifted moduli fields:

$${\tilde{\beta }}^i\left(x\right)={\beta }^i\left(x\right)-{\beta }_0^i,$$

(6)

where ${\beta }_0^i$ denote the stabilized (or present-day) values. These moduli will play a central role in the effective external theory and influence the spectral properties of the fluctuation operator (cf. [31]).

3. Conformal Rescaling and Dimensional Reduction

To obtain a $\left(D_0\right)$-dimensional effective action in the Einstein frame, we perform a conformal rescaling of the external metric:

$$g_{\mu \nu }^{\left(0\right)}\left(x\right)={\Omega }^2\left(x\right){\tilde{g}}_{\mu \nu }^{\left(0\right)}\left(x\right),\mathrm{with}{\Omega }^2\left(x\right)=exp\left[-{\displaystyle \frac{2}{D_0-2}}\sum _{i=1}^n{d}_i{\tilde{\beta }}^i\left(x\right)\right].$$

(7)

Integration over the internal coordinates yields the internal volume,

$$V_{D^{\prime }}=\prod _{i=1}^n{\int }_{M_i}{d}^{d_i}y\sqrt{|g^{\left(i\right)}|},$$

(8)

which relates the D-dimensional gravitational coupling ${\kappa }_D$ to its $\left(D_0\right)$-dimensional counterpart,

$${\kappa }_0^2={\displaystyle \frac{{\kappa }_D^2}{V_{D^{\prime }}}}.$$

(9)

Thus, the effective action takes the preliminary form

$$S_{\mathrm{eff}}={\displaystyle \frac{1}{2{\kappa }_0^2}}{\int }_{M_0}{d}^{D_0}x\sqrt{|{\tilde{g}}^{\left(0\right)}|}\left\{\tilde{R}\left[{\tilde{g}}^{\left(0\right)}\right]-{\overline{G}}_{ij}{\tilde{g}}^{\left(0\right)\mu \nu }{\partial }_{\mu }{\tilde{\beta }}^i{\partial }_{\nu }{\tilde{\beta }}^j-2U_{\mathrm{eff}}\left(\tilde{\beta }\right)\right\},$$

(10)

where the moduli space metric is given by

$${\overline{G}}_{ij}=d_i{\delta }_{ij}+{\displaystyle \frac{d_i{d}_j}{D_0-2}}.$$

(11)

4. Effective Potential and Loop Corrections

The effective potential in Eq. (10) arises from both classical curvature contributions and quantum effects. Explicitly, we have

$$U_{\mathrm{eff}}\left(\tilde{\beta }\right)=exp\left[-{\displaystyle \frac{2}{D_0-2}}\sum _{i=1}^n{d}_i{\tilde{\beta }}^i\left(x\right)\right]\left[-{\displaystyle \frac{1}{2}}\sum _{i=1}^n{\tilde{R}}_i{e}^{-2{\tilde{\beta }}^i\left(x\right)}+\Lambda \right],$$

(12)

with ${\tilde{R}}_i=R_i{e}^{-2{\beta }_0^i}$. In addition, one-loop quantum corrections from heavy Kaluza–Klein modes contribute to the potential through the functional determinant

$$V_{1-\mathrm{loop}}\left(\phi \right)={\displaystyle \frac{1}{2}}Trlog\left(-\square +m^2\left(\phi \right)\right),$$

(13)

which, upon employing zeta-function regularization and the heat kernel expansion [32,33,34], can be expressed as

$$V_{1-\mathrm{loop}}\left(\phi \right)={\displaystyle \frac{1}{2{\left(4\pi \right)}^{D_0/2}}}\sum _{k=0}^{\infty }{a}_k\left(\phi \right)\Gamma \left(k-{\displaystyle \frac{D_0}{2}}\right){\left[m^2\left(\phi \right)\right]}^{{\displaystyle \frac{D_0}{2}}-k}.$$

(14)

This expression, which also determines the renormalization group flow of the effective couplings [35], necessitates analytic continuation for odd $D_0$.

5. Canonical Normalization in the Single Modulus Case

To explain the key features of our framework, we now focus on the single modulus case ($n=1$). In this case, we define

$$\tilde{\beta }\equiv {\tilde{\beta }}^1,d\equiv d_1,\tilde{R}\equiv {\tilde{R}}_1.$$

(15)

We introduce a canonically normalized scalar field $\phi \left(x\right)$ via

$$\phi \left(x\right)=Q\tilde{\beta }\left(x\right),$$

(16)

with the normalization constant Q determined by the condition

$$Q^2\overline{G}=1,\mathrm{with}\overline{G}=d+{\displaystyle \frac{d^2}{D_0-2}}.$$

(17)

6. Perturbative Expansion, Self-Interactions, and Higher-Derivative Corrections

Assuming that the effective potential is stabilized at $\tilde{\beta }=0$ (or equivalently $\phi =0$), we perform a Taylor expansion around the minimum:

$$U_{\mathrm{eff}}\left(\phi \right)=U_{\mathrm{eff}}\left(0\right)+{\displaystyle \frac{1}{2}}{U}_{\mathrm{eff}}^{\prime \prime }\left(0\right){\phi }^2+{\displaystyle \frac{1}{3!}}{U}_{\mathrm{eff}}^{\left(3\right)}\left(0\right){\phi }^3+{\displaystyle \frac{1}{4!}}{U}_{\mathrm{eff}}^{\left(4\right)}\left(0\right){\phi }^4+\cdots .$$

(18)

Motivated by quantum corrections and the integration of heavy modes [27,28], we parameterize the potential as

$$U_{\mathrm{eff}}\left(\phi \right)={\Lambda }_{\mathrm{eff}}+{\displaystyle \frac{1}{2}}{m}^2{\phi }^2+{\displaystyle \frac{\lambda }{4}}{\phi }^4+{\displaystyle \frac{{\lambda }_6}{6!}}{\phi }^6+{\displaystyle \frac{{\lambda }_8}{8!}}{\phi }^8+\mathcal{O}\left({\phi }^{10}\right),$$

(19)

where ${\Lambda }_{\mathrm{eff}}=U_{\mathrm{eff}}\left(0\right)$ and $m^2=U_{\mathrm{eff}}^{\prime \prime }\left(0\right)$, while the higher-order couplings encode further stringy and quantum corrections.

In addition to the standard kinetic term, extra-dimensional effective theories naturally include higher-derivative corrections. The leading correction is given by

$$\Delta {\mathcal{L}}_{\mathrm{deriv}}=\alpha {\left({\tilde{g}}^{\left(0\right)\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi \right)}^2,$$

(20)

with $\alpha$ a coupling constant of appropriate mass dimension. Such corrections affect the ultraviolet behavior of the theory and may modify scattering amplitudes and vacuum stability [20,36].

Spectral Analysis of the Fluctuation Operator.

To study the impact of higher-derivative terms on the spectrum of moduli fluctuations, consider the eigenvalue equation for the modified operator:

$$\mathcal{O}\delta \phi =\lambda \delta \phi ,\mathrm{with}\mathcal{O}=-\square +m^2+\xi \tilde{R}\left[{\tilde{g}}^{\left(0\right)}\right]+\alpha {\square }^2.$$

(21)

Employing a plane-wave ansatz $\delta \phi \sim exp\left(i{k}_{\mu }{x}^{\mu }\right)$, the dispersion relation becomes

$$\alpha {\left(k^2\right)}^2+k^2+m^2+\xi \tilde{R}=0.$$

(22)

The quartic nature of this equation in $k^2$ requires the use of methods akin to Cardano’s approach for depressed quartic equations [37], revealing branch cuts in the propagator and indicating potential instabilities unless $\alpha$ is sufficiently suppressed. This analysis mirrors the nonlocal operator techniques discussed in [31].

7. Non-Minimal Couplings and Interactions with Matter

The interplay between curvature and moduli dynamics is incorporated via a non-minimal coupling between the scalar field and the Ricci scalar:

$$\Delta {\mathcal{L}}_{\mathrm{NM}}=\xi \tilde{R}\left[{\tilde{g}}^{\left(0\right)}\right]{\phi }^2,$$

(23)

with $\xi$ a dimensionless parameter. Such couplings naturally emerge in theories with conformal symmetry and have important implications in both early- and late-time cosmology [38,39]. Additional interactions with matter fields are also included. For example, a Yukawa interaction with fermions is given by

$$\Delta {\mathcal{L}}_{\mathrm{Yukawa}}=g_Y\overline{\psi }\psi \phi ,$$

(24)

and a coupling to gauge fields is described by

$$\Delta {\mathcal{L}}_{\mathrm{gauge}}=\zeta \phi F_{\mu \nu }{F}^{\mu \nu }.$$

(25)

These interactions not only open novel decay channels but also influence the renormalization group evolution of the effective couplings [40,41].

8. Flux Stabilization, Topology, and Non-Perturbative Contributions

If the internal space exhibits non-trivial topology (e.g., orbifolds or Calabi–Yau manifolds) or supports background fluxes [42,43], additional contributions to the effective potential arise. Flux-induced superpotentials may introduce exponential or other non-polynomial corrections that aid in moduli stabilization. Non-perturbative effects, such as instanton contributions and vacuum tunneling, further enrich the vacuum structure, potentially leading to metastable vacua of phenomenological interest [44].

9. Final Extended Effective Action

Collecting all contributions, the complete effective action for the canonically normalized scalar field $\phi$ in $D_0$ dimensions is given by

$$\begin{array}{cc}\hfill S_{\mathrm{eff}}& ={\displaystyle \frac{1}{2{\kappa }_0^2}}{\int }_{M_0}{d}^{D_0}x\sqrt{|{\tilde{g}}^{\left(0\right)}|}\left\{\tilde{R}\left[{\tilde{g}}^{\left(0\right)}\right]-2{\Lambda }_{\mathrm{eff}}-\xi \tilde{R}\left[{\tilde{g}}^{\left(0\right)}\right]{\phi }^2\right\}\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\int }_{M_0}{d}^{D_0}x\sqrt{|{\tilde{g}}^{\left(0\right)}|}\{{\tilde{g}}^{\left(0\right)\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi +m^2{\phi }^2+\lambda {\phi }^4\hfill \\ \hfill & +{\lambda }_6{\phi }^6+{\lambda }_8{\phi }^8+\alpha {\left({\tilde{g}}^{\left(0\right)\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi \right)}^2+\cdots \}\hfill \\ \hfill & +S_{\mathrm{int}},\hfill \end{array}$$

(26)

where $S_{\mathrm{int}}$ collects additional interaction terms (such as the Yukawa and gauge couplings given in Eqs. (24) and (25)), as well as further corrections from fluxes, topological effects, and higher-loop contributions.

This extended effective action forms the foundation for exploring a broad spectrum of phenomena, including gravitational exciton dynamics, moduli stabilization, inflationary scenarios, and late-time cosmic acceleration. Notably, the inclusion of higher-order self-interactions and derivative corrections permits the possibility of Bose–Einstein condensates (BECs) of gravitational excitations, with potential observational signatures [24].

Figure 1 displays a representative effective potential,

$$U_{\mathrm{eff}}\left(\phi \right)={\Lambda }_{\mathrm{eff}}+{\displaystyle \frac{1}{2}}{m}^2{\phi }^2+{\displaystyle \frac{\lambda }{4}}{\phi }^4+{\displaystyle \frac{{\lambda }_6}{6!}}{\phi }^6+{\displaystyle \frac{{\lambda }_8}{8!}}{\phi }^8,$$

for the parameter choices

$${\Lambda }_{\mathrm{eff}}=0,m^2=1,\lambda =1,{\lambda }_6=0.1,{\lambda }_8=0.01.$$

The presence of sextic and octic terms yields a rich vacuum structure that may facilitate nontrivial phase transitions and metastable states.

3. Quantum Corrections and Renormalization Group Effects

Having established the classical extended effective action, we now incorporate quantum corrections within the framework of effective field theory. Our analysis employs the Coleman–Weinberg mechanism alongside gravitational loop effects to derive one-loop corrections and study the ensuing renormalization group (RG) evolution of both self-couplings and the moduli parameters controlling the compactification scale [19,20,21].

1. One-Loop Effective Potential: Matter and Graviton Contributions

Starting from the classical effective potential,

$$U_{\mathrm{eff}}^{\mathrm{cl}}\left(\phi \right)={\Lambda }_{\mathrm{eff}}+{\displaystyle \frac{1}{2}}{m}^2{\phi }^2+{\displaystyle \frac{\lambda }{4}}{\phi }^4,$$

(27)

quantum corrections are obtained by evaluating the one-loop functional determinant. For a scalar field fluctuating about the background $\phi$, the matter loop contribution is given by the Coleman–Weinberg formula:

$$V_{1-\mathrm{loop}}^{\mathrm{matter}}\left(\phi \right)={\displaystyle \frac{1}{64{\pi }^2}}{\left[M^2\left(\phi \right)\right]}^2\left[ln\left({\displaystyle \frac{M^2\left(\phi \right)}{{\mu }^2}}\right)-{\displaystyle \frac{3}{2}}\right],$$

(28)

with the field-dependent mass-squared defined as

$$M^2\left(\phi \right)=m^2+3\lambda {\phi }^2,$$

(29)

and $\mu$ representing the renormalization scale. A more detailed derivation follows from the background field method. In dimensional regularization the functional determinant

$$\Delta U_{\mathrm{eff}}^{\mathrm{matter}}\left(\phi \right)={\displaystyle \frac{1}{2}}\int {\displaystyle \frac{d^{D_0}k}{{\left(2\pi \right)}^{D_0}}}ln\left[k^2+M^2\left(\phi \right)\right],$$

reduces to Eq. (28) after proper subtraction of divergences [45,46].

Graviton loops, computed within effective quantum gravity using e.g. the Vilkovisky–DeWitt formalism, contribute additional logarithmic and non-local corrections suppressed by inverse powers of the Planck mass $M_{Pl}$. Schematically, we express the gravitational contribution as

$$V_{1-\mathrm{loop}}^{\mathrm{grav}}\left(\phi \right)={\displaystyle \frac{1}{64{\pi }^2}}{f}_{\mathrm{grav}}(\phi ,\mu ,M_{Pl}),$$

(30)

so that the full one-loop effective potential becomes

$$\begin{array}{cc}\hfill U_{\mathrm{eff}}\left(\phi \right)& =U_{\mathrm{eff}}^{\mathrm{cl}}\left(\phi \right)+\Delta U_{\mathrm{eff}}^{\mathrm{matter}}\left(\phi \right)+\Delta U_{\mathrm{eff}}^{\mathrm{grav}}\left(\phi \right)\hfill \\ \hfill & ={\Lambda }_{\mathrm{eff}}+{\displaystyle \frac{1}{2}}{m}^2{\phi }^2+{\displaystyle \frac{\lambda }{4}}{\phi }^4\hfill \\ \hfill & +{\displaystyle \frac{1}{64{\pi }^2}}{\left[m^2+3\lambda {\phi }^2\right]}^2\left[ln\left({\displaystyle \frac{m^2+3\lambda {\phi }^2}{{\mu }^2}}\right)-{\displaystyle \frac{3}{2}}\right]\hfill \\ \hfill & +{\displaystyle \frac{1}{64{\pi }^2}}{f}_{\mathrm{grav}}(\phi ,\mu ,M_{Pl}).\hfill \end{array}$$

(31)

Here, the function $f_{\mathrm{grav}}(\phi ,\mu ,M_{Pl})$ encapsulates the graviton loop corrections computed via methods described in [47,48,49].

2. Emergent Scale Invariance and the Dilaton Mode

An intriguing feature of the quantum-corrected potential is that, over a suitable range of parameters, the logarithmic terms induce an approximate scale invariance. In this regime the potential becomes nearly flat over an extended region, thereby facilitating spontaneous symmetry breaking and the emergence of a light dilaton mode [10,46,50,51]. Such a dilaton may address the hierarchy problem and have implications for dark energy phenomenology.

3. Renormalization Group Flow in a Multidimensional Moduli Space

Quantum corrections render the couplings scale dependent. The renormalization group evolution is governed by the Callan–Symanzik equation,

$$\left(\mu {\displaystyle \frac{\partial }{\partial \mu }}+{\beta }_{m^2}{\displaystyle \frac{\partial }{\partial m^2}}+{\beta }_{\lambda }{\displaystyle \frac{\partial }{\partial \lambda }}-\gamma \phi {\displaystyle \frac{\partial }{\partial \phi }}\right)U_{\mathrm{eff}}\left(\phi \right)=0,$$

(32)

with the beta functions for the scalar mass and self-coupling defined as

$${\beta }_{m^2}\equiv \mu {\displaystyle \frac{d{m}^2}{d\mu }}={\displaystyle \frac{3\lambda }{16{\pi }^2}}{m}^2+{\beta }_{m^2}^{\mathrm{grav}},{\beta }_{\lambda }\equiv \mu {\displaystyle \frac{d\lambda }{d\mu }}={\displaystyle \frac{9{\lambda }^2}{16{\pi }^2}}+{\beta }_{\lambda }^{\mathrm{grav}}.$$

(33)

Additional gravitational corrections ${\beta }_{m^2}^{\mathrm{grav}}$ and ${\beta }_{\lambda }^{\mathrm{grav}}$ are computed within the effective theory of gravity [20,21,45,52]. In a multidimensional moduli space, the RG flow becomes even richer, with couplings mixing as

$$\mu {\displaystyle \frac{d\overrightarrow{g}}{d\mu }}=\mathbf{B}\left(\overrightarrow{g}\right),$$

(34)

where $\overrightarrow{g}=\{m^2,\lambda ,\cdots \}$ and the matrix $\mathbf{B}$ encodes contributions from both matter and graviton loops. Fixed points and critical behavior in this flow may lead to cascades of phase transitions and are crucial for extra-dimensional stabilization.

4. Non-Perturbative Phenomena and Vacuum Tunneling

Beyond perturbation theory, non-perturbative effects such as instanton-induced vacuum tunneling between distinct compactification vacua can alter the effective potential dramatically. In the dilute gas approximation, instanton contributions modify the potential as

$$\Delta U_{\mathrm{eff}}^{\mathrm{inst}}\left(\phi \right)\sim A\left(\phi \right)exp\left[-{\displaystyle \frac{S_{\mathrm{inst}}\left(\phi \right)}{\hslash }}\right],$$

(35)

where $S_{\mathrm{inst}}\left(\phi \right)$ is the instanton action and $A\left(\phi \right)$ is the fluctuation determinant. These non-analytic terms can trigger first-order phase transitions, as originally discussed in [53,54,55].

5. Quantum-Corrected Parameters and Renormalization Conditions

To quantify the impact of quantum corrections, we expand the effective potential about its minimum at $\phi =0$ and define the renormalized parameters via

$$m_{\mathrm{eff}}^2=U_{\mathrm{eff}}^{\prime \prime }\left(0\right)=m^2+\Delta m^2,{\lambda }_{\mathrm{eff}}={\displaystyle \frac{1}{6}}{U}_{\mathrm{eff}}^{\left(4\right)}\left(0\right)=\lambda +\Delta \lambda .$$

(36)

A straightforward calculation yields, for the mass correction,

$$\Delta m^2={\displaystyle \frac{3\lambda }{32{\pi }^2}}{m}^2\left[ln\left({\displaystyle \frac{m^2}{{\mu }^2}}\right)-1\right]+\Delta m_{\mathrm{grav}}^2,$$

(37)

with analogous expressions for $\Delta \lambda$, including both matter and gravitational contributions [20,21]. The renormalization conditions guarantee that physical observables remain independent of the arbitrary scale $\mu$ when an RG-improved effective potential is employed.

6. Numerical Analysis and RG Scale Dependence

To illustrate the interplay between quantum corrections and RG evolution, we perform a numerical minimization of the effective potential for several renormalization scales. Table 1 lists the computed location of the minimum ${\phi }_{\min }$, the value of the potential at the minimum $U_{\min }$, and the quantum-corrected mass squared $m_{\mathrm{eff}}^2$ for $\mu =0.5$, 1.0, and 2.0.

Figure 2 and Figure 3 illustrate our findings. In Figure 2, we compare the classical potential (Eq. (27)) with the one-loop corrected potential (Eq. (31)) at $\mu =1.0$. Figure 3 displays the evolution of the effective potential as the renormalization scale is varied, thus reflecting the running of the couplings.

In summary, the interplay between matter and gravitational loop corrections, the emergence of approximate scale invariance, and non-perturbative phenomena enrich the structure of the effective potential. The resulting RG flow, with its fixed points and possible phase transitions, offers valuable insights into exciton dynamics, extra-dimensional stabilization, and potential cosmological implications.

4. Non-Relativistic Limit and Condensate Dynamics

Having derived the RG–improved extended effective action in the relativistic regime, we now focus on its non-relativistic limit and the corresponding condensate dynamics. This limit is particularly relevant for astrophysical applications (e.g., dark matter halos, boson stars) where gravitational exciton modes are expected to condense. For simplicity, we restrict the full theory to four spacetime dimensions ($D_0=4$) and assume a fixed Minkowski background, ${\tilde{g}}_{\mu \nu }^{\left(0\right)}={\eta }_{\mu \nu }$, thereby neglecting gravitational backreaction on the condensate dynamics.

1. Reduction to the Exciton Sector

In the non-relativistic limit the exciton sector is well described by the action

$$S_{\phi }=-{\displaystyle \frac{1}{2}}\int d^{4}x\left\{{\eta }^{\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi +m_{\mathrm{eff}}^2{\phi }^2+{\lambda }_{\mathrm{eff}}{\phi }^4\right\},$$

(38)

where $m_{\mathrm{eff}}$ and ${\lambda }_{\mathrm{eff}}$ are the quantum–corrected parameters determined from the RG analysis in previous sections. Variation of Eq. (38) with respect to $\phi$ yields a modified Klein–Gordon equation. Expanding the field in Fourier modes and performing a systematic expansion in powers of $\left|\mathbf{k}\right|/m_{\mathrm{eff}}$ (cf. [15,16]) sets the stage for a non-relativistic reduction.

2. Non-Relativistic Ansatz and Derivation of the Gross–Pitaevskii Equation

In the low-energy regime the field $\phi$ is dominated by its positive-frequency components. We therefore adopt the standard ansatz [23]:

$$\phi (\mathbf{x},t)={\displaystyle \frac{1}{\sqrt{2m_{\mathrm{eff}}}}}\left[e^{-i{m}_{\mathrm{eff}}t/\hslash }\psi (\mathbf{x},t)+e^{i{m}_{\mathrm{eff}}t/\hslash }{\psi }^{*}(\mathbf{x},t)\right],$$

(39)

where $\psi (\mathbf{x},t)$ is a slowly varying complex amplitude relative to the rapid oscillations encoded in $e^{\pm i{m}_{\mathrm{eff}}t/\hslash }$. Substituting Eq. (39) into the modified Klein–Gordon equation and neglecting second-order time derivatives of $\psi$ (i.e. assuming ${\partial }_t\psi \ll m_{\mathrm{eff}}\psi$) leads, after averaging over the fast oscillations, to the Gross–Pitaevskii (GP) equation:

$$i\hslash {\partial }_t\psi (\mathbf{x},t)=\left[-{\displaystyle \frac{{\hslash }^2}{2m_{\mathrm{eff}}}}{\nabla }^2+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{\left|\psi (\mathbf{x},t)\right|}^2\right]\psi (\mathbf{x},t),$$

(40)

with the effective nonlinearity parameter given by

$${\lambda }_{\mathrm{eff}}^{\mathrm{NR}}={\displaystyle \frac{3{\lambda }_{\mathrm{eff}}}{2m_{\mathrm{eff}}}}.$$

(41)

The numerical prefactor in Eq. (41) arises from the particular normalization in Eq. (39) and the projection of the quartic term onto the slowly varying mode [56].

3. Inclusion of Self-Gravity: The Coupled GP–Poisson System

In many astrophysical and cosmological settings the condensate is self-gravitating. To account for this, we couple the GP equation to the Poisson equation for the gravitational potential $\Phi (\mathbf{x},t)$. The resulting coupled system is

$$\begin{array}{cc}\hfill i\hslash {\partial }_t\psi (\mathbf{x},t)& =\left[-{\displaystyle \frac{{\hslash }^2}{2m_{\mathrm{eff}}}}{\nabla }^2+m_{\mathrm{eff}}\Phi (\mathbf{x},t)+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{\left|\psi (\mathbf{x},t)\right|}^2\right]\psi (\mathbf{x},t),\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\nabla }^2\Phi (\mathbf{x},t)& =4\pi G{m}_{\mathrm{eff}}{\left|\psi (\mathbf{x},t)\right|}^2.\hfill \end{array}$$

(43)

Here, the condensate density is defined as $\rho (\mathbf{x},t)=m_{\mathrm{eff}}{\left|\psi (\mathbf{x},t)\right|}^2$. The GP–Poisson system has been extensively studied in the context of Bose–Einstein condensate (BEC) dark matter and boson star formation (see, e.g., [57,58,59]).

4. Finite-Temperature Effects and Phase Transition Dynamics

Realistic condensates are subject to thermal fluctuations. Finite-temperature field theory techniques [60,61] enable the derivation of a temperature–dependent GP equation. Near a critical temperature $T_c$, the condensate vanishes and the system undergoes a second–order phase transition characterized by universal critical exponents [56]. In our framework, finite-temperature corrections can be incorporated either by adding a stochastic noise term to Eq. (40) (as in [62]) or by promoting the effective coupling ${\lambda }_{\mathrm{eff}}^{\mathrm{NR}}$ to a temperature–dependent function. These thermal effects are particularly significant in the early Universe, where they may influence structure formation and the stabilization of extra dimensions [63].

5. Rotating Condensates and Vortex Nucleation

Angular momentum plays a crucial role in many astrophysical scenarios. In a rotating frame the GP equation acquires an additional term corresponding to the Coriolis force. The modified equation is

$$i\hslash {\partial }_t\psi (\mathbf{x},t)=\left[-{\displaystyle \frac{{\hslash }^2}{2m_{\mathrm{eff}}}}{\nabla }^2+m_{\mathrm{eff}}\Phi (\mathbf{x},t)-\Omega L_z+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{\left|\psi (\mathbf{x},t)\right|}^2\right]\psi (\mathbf{x},t),$$

(44)

where $\Omega$ is the angular velocity and $L_z=-i\hslash (x{\partial }_y-y{\partial }_x)$ denotes the z–component of the angular momentum operator. Rotation can induce the nucleation of quantized vortices and the formation of vortex lattices, phenomena well documented in laboratory BECs [64,65] and proposed as mechanisms affecting galactic rotation curves and boson star dynamics [57,58].

6. Numerical Simulations and Coupling with External Fields

To investigate the condensate dynamics under realistic conditions, we performed numerical simulations of the coupled Gross-Pitaevskii (GP)–Poisson system, as described by Eqs. (42) and (43). These simulations take into account several factors: finite-temperature effects, modeled by introducing a stochastic noise term in the GP equation; self-gravity, which is implemented through a self-consistent solution of the Poisson equation; and potential rotational effects, incorporated by including the Coriolis term in Eq. (44), with the rotation parameter $\Omega$ set to zero unless stated otherwise.

Table 2 summarizes key observables from our simulations, including the maximum density, mean density, and the extremal values of the gravitational potential at representative simulation times $t=10$ and $t=20$.

Figure 4 and Figure 5 display snapshots of the condensate density and gravitational potential at times $t=10$ and $t=20$, respectively. At $t=10$, pronounced density fluctuations (localized peaks) correlate with deep gravitational potential wells, reflecting the competition between self–gravity and quantum pressure. By $t=20$, these fluctuations merge and dissipate, yielding a smoother density profile and a shallower potential indicative of equilibration.

In summary, by taking the non-relativistic limit of the RG–improved effective action, we arrive at a Gross–Pitaevskii description of gravitational exciton condensates. The coupled GP–Poisson system, enriched by finite-temperature corrections and potential rotational effects, offers a robust framework for exploring astrophysical phenomena such as dark matter halo formation, boson star evolution, and vortex dynamics. Future work will incorporate external fields and investigate nontrivial topological defects within this condensate framework.

5. Stationary Solutions and Stability Analysis

To understand the long–term behavior and astrophysical implications of gravitational exciton condensates, it is crucial to study their stationary configurations and analyze their stability. In this section, we investigate both spherically symmetric and more complex solutions (including rotating, anisotropic, and solitonic states), which can model realistic structures such as dark matter halos or boson stars. We further perform a linear stability analysis via the Bogoliubov–de Gennes formalism and discuss the Vakhitov–Kolokolov criterion for stability [66,67,68].

1. Stationary Ansatz and Time–Independent Equations

We begin by seeking stationary solutions of the coupled Gross–Pitaevskii–Poisson system. Assuming a harmonic time dependence,

$$\psi (\mathbf{x},t)=e^{-i\mu t/\hslash }{\psi }_0\left(\mathbf{x}\right),$$

(45)

where $\mu$ is the chemical potential and ${\psi }_0\left(\mathbf{x}\right)$ is time–independent. Substituting Eq. (45) into the time–dependent Gross–Pitaevskii equation (cf. Eq. (42)) yields the time–independent GP equation:

$$\mu {\psi }_0\left(\mathbf{x}\right)=\left[-{\displaystyle \frac{{\hslash }^2}{2m_{\mathrm{eff}}}}{\nabla }^2+m_{\mathrm{eff}}\Phi \left(\mathbf{x}\right)+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{\left|{\psi }_0\left(\mathbf{x}\right)\right|}^2\right]{\psi }_0\left(\mathbf{x}\right),$$

(46)

with the gravitational potential determined self–consistently by

$${\nabla }^2\Phi \left(\mathbf{x}\right)=4\pi G{m}_{\mathrm{eff}}{\left|{\psi }_0\left(\mathbf{x}\right)\right|}^2.$$

(47)

The condensate wavefunction is normalized via

$$\int d^{3}x{\left|{\psi }_0\left(\mathbf{x}\right)\right|}^2=N,$$

(48)

where N is the total number of excitons. Equations (46) and (47) thus form a nonlinear eigenvalue problem for ${\psi }_0\left(\mathbf{x}\right)$ and $\mu$.

2. Spherically Symmetric Configurations

For spherically symmetric solutions, we set ${\psi }_0\left(\mathbf{x}\right)={\psi }_0\left(r\right)$ and $\Phi \left(\mathbf{x}\right)=\Phi \left(r\right)$, with $r=\left|\mathbf{x}\right|$. In spherical coordinates, the Laplacian simplifies to

$${\nabla }^{2}f\left(r\right)={\displaystyle \frac{d^{2}f}{d{r}^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{df}{dr}}.$$

(49)

Thus, the time–independent GP equation reduces to

$$\mu {\psi }_0\left(r\right)=\left[-{\displaystyle \frac{{\hslash }^2}{2m_{\mathrm{eff}}}}\left({\displaystyle \frac{d^2}{d{r}^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{d}{dr}}\right)+m_{\mathrm{eff}}\Phi \left(r\right)+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{\left|{\psi }_0\left(r\right)\right|}^2\right]{\psi }_0\left(r\right),$$

(50)

and the corresponding Poisson equation becomes

$${\displaystyle \frac{1}{r^2}}{\displaystyle \frac{d}{dr}}\left(r^2{\displaystyle \frac{d\Phi \left(r\right)}{dr}}\right)=4\pi G{m}_{\mathrm{eff}}{\left|{\psi }_0\left(r\right)\right|}^2.$$

(51)

Multiplying Eq. (50) by ${\psi }_0^{*}\left(r\right)$ and integrating yields the energy functional of the system. A virial theorem can be derived by scaling arguments [57], providing further insights into the balance between kinetic, self–interaction, and gravitational energy contributions.

3. Extensions to Rotating and Anisotropic States

Many astrophysical systems exhibit non–spherical morphologies. To model such systems, one may generalize the stationary ansatz by incorporating angular momentum:

$$\psi (\mathbf{x},t)=e^{-i\mu t/\hslash }{e}^{i\ell \theta }{\psi }_0(r,\theta ,\varphi ),$$

(52)

where ℓ denotes the quantized angular momentum. This ansatz modifies the kinetic energy by introducing a centrifugal barrier and, for anisotropic deformations, one may adopt a variational ansatz with anisotropic Gaussian profiles [64]. Localized solitonic states (such as Q–balls) can also be obtained by minimizing the energy functional

$$E\left[{\psi }_0\right]=\int d^{3}x\left[{\displaystyle \frac{{\hslash }^2}{2m_{\mathrm{eff}}}}|\nabla {\psi }_0{|}^2+{\displaystyle \frac{m_{\mathrm{eff}}}{2}}\Phi |{\psi }_0{|}^2+{\displaystyle \frac{{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}}{2}}{\left|{\psi }_0\right|}^4\right],$$

(53)

subject to the normalization constraint (48). Variational methods and the Bogoliubov–de Gennes formalism [23] can then be employed to map the stability region in parameter space, taking into account the scale dependence induced by RG corrections [20,21].

4. Linear Stability Analysis

To assess the dynamical stability of a stationary solution ${\psi }_0\left(\mathbf{x}\right)$, we consider small perturbations:

$$\psi (\mathbf{x},t)=e^{-i\mu t/\hslash }\left[{\psi }_0\left(\mathbf{x}\right)+u\left(\mathbf{x}\right)e^{-i\omega t}+v^{*}\left(\mathbf{x}\right)e^{i\omega t}\right].$$

(54)

Inserting this into the time–dependent GP equation and linearizing in u and v, we obtain the Bogoliubov–de Gennes equations. The resulting eigenvalue problem for the excitation frequencies $\omega$ determines the stability of the stationary state. In particular, the Vakhitov–Kolokolov (VK) criterion,

$${\displaystyle \frac{dN}{d\mu }}<0,$$

(55)

provides a necessary condition for stability [67]. A violation of this condition typically signals the onset of collapse or modulational instability [66].

5. Numerical Stationary Solutions and Convergence

To obtain explicit stationary solutions, we numerically solve Eqs. (50) and (51) using an imaginary time evolution method on a finite–difference grid in spherical coordinates. In a typical simulation, convergence is achieved after a simulation time of approximately 225.80 (in dimensionless units), yielding a ground state characterized by a chemical potential $\mu \approx 0.05057$ and a total energy $E\approx 0.00587$. The stationary density profile $|{\psi }_0{\left(r\right)|}^2$ exhibits a maximum of about 0.00336 with a mean density near 0.000564, while the gravitational potential $\Phi \left(r\right)$ attains a minimum value of roughly $-0.00588$ at $r=0$ and smoothly approaches zero at large r. These results indicate that the condensate relaxes into a diffuse ground state with a central density peak and a shallow potential well.

Figure 6 (left panel) shows the converged stationary density (blue solid line) and gravitational potential (red dashed line) as functions of r, and the right panel illustrates the convergence of the energy functional during the imaginary time evolution, confirming that a stable state is reached.

In summary, our stability analysis demonstrates that gravitational exciton condensates support a rich variety of stationary configurations. The combination of numerical solutions in the spherically symmetric case, extensions to rotating and anisotropic states, and a linear stability analysis based on the Bogoliubov–de Gennes equations provide a robust framework for exploring how RG–improved quantum corrections affect the stabilization of self–gravitating condensates. These results are essential for modeling dark matter structures and boson stars in astrophysical contexts.

6. Bogoliubov Excitations and Dispersion Relation

Collective excitations play a key role in determining the dynamical response and stability of gravitational exciton condensates. In this section, we analyze small fluctuations about a homogeneous condensate state and derive the corresponding Bogoliubov dispersion relation. This analysis not only confirms the presence of phonon–like modes at low momenta but also delineates the crossover to free–particle behavior at higher momenta, thereby providing insight into both the low–energy collective dynamics and high–energy single–particle excitations [22,23].

1. Linearization about a Homogeneous Condensate

Assuming that the condensate is spatially homogeneous, we set

$${\psi }_0\left(\mathbf{x}\right)=\sqrt{n_0},$$

(56)

with constant density $n_0$, and neglect any gravitational potential variations over the relevant length scales. To study small perturbations, we write

$$\psi (\mathbf{x},t)=e^{-i\mu t/\hslash }\left[\sqrt{n_0}+\delta \psi (\mathbf{x},t)\right],$$

(57)

where $\mu$ is the chemical potential. Decomposing the fluctuation into plane–wave modes,

$$\delta \psi (\mathbf{x},t)=u_{\mathbf{k}}{e}^{i(\mathbf{k}·\mathbf{x}-\omega t)}+v_{\mathbf{k}}^{*}{e}^{-i(\mathbf{k}·\mathbf{x}-\omega t)},$$

(58)

and substituting into the time–dependent Gross–Pitaevskii equation [22,23], we retain only linear terms in the fluctuations. This procedure yields the Bogoliubov–de Gennes equations:

$$\begin{array}{cc}\hfill \hslash \omega u_{\mathbf{k}}& =\left[{ϵ}_k+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0\right]u_{\mathbf{k}}+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0{v}_{\mathbf{k}},\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill -\hslash \omega v_{\mathbf{k}}& =\left[{ϵ}_k+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0\right]v_{\mathbf{k}}+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0{u}_{\mathbf{k}},\hfill \end{array}$$

(60)

where the free–particle kinetic energy is defined as

$${ϵ}_k={\displaystyle \frac{{\hslash }^2{k}^2}{2m_{\mathrm{eff}}}}.$$

(61)

2. Derivation of the Bogoliubov Dispersion Relation

Writing Eqs. (59)–(60) in matrix form,

$$\left(\begin{array}{cc}{ϵ}_k+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0& {\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0\\ -{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0& -\left({ϵ}_k+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0\right)\end{array}\right)\left(\begin{array}{c}{u}_{\mathbf{k}}\\ v_{\mathbf{k}}\end{array}\right)=\hslash \omega \left(\begin{array}{c}{u}_{\mathbf{k}}\\ v_{\mathbf{k}}\end{array}\right),$$

(62)

non–trivial solutions exist if the determinant of the coefficient matrix vanishes. This condition leads to

$${\hslash }^2{\omega }^2={\left[{ϵ}_k+{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0\right]}^2-{\left[{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0\right]}^2={ϵ}_k\left({ϵ}_k+2{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0\right).$$

(63)

This is the celebrated Bogoliubov dispersion relation originally derived in the context of dilute Bose gases [22].

3. Long–Wavelength Limit: Phonon Modes

In the long–wavelength limit where ${ϵ}_k\ll {\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0$, Eq. (63) simplifies to

$$\hslash \omega \approx \hslash k\sqrt{{\displaystyle \frac{{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0}{m_{\mathrm{eff}}}}}.$$

(64)

This linear dispersion defines the speed of sound in the condensate,

$$c_s=\sqrt{{\displaystyle \frac{{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0}{m_{\mathrm{eff}}}}},$$

(65)

so that the low–momentum excitation spectrum is given by

$$\omega \approx c_{s}k.$$

(66)

These phonon–like modes are essential for understanding superfluid behavior and the low–energy response of the condensate [56].

4. Finite-Temperature and Nonlinear Corrections

While the derivation above is valid at zero temperature and within the linear regime, finite-temperature effects can lead to damping (via mechanisms such as Landau damping) and modify the dispersion relation by establishing a Landau critical velocity below which excitations remain undamped [56]. Furthermore, nonlinear corrections—arising from higher–order terms in the fluctuation expansion—can result in the formation of solitons, shock waves, and other coherent structures. Recent studies have also proposed that Bogoliubov excitations might couple to gravitational waves, providing an intriguing observational window into the properties of gravitational exciton condensates [24,25].

5. Numerical Results and Discussion

Figure 7 displays the Bogoliubov dispersion relation $\hslash \omega \left(k\right)$ as computed from Eq. (63) for representative parameters: ${\lambda }_{\mathrm{eff}}^{\mathrm{NR}}=1.5000$ and $n_0=2.5\times {10}^{-3}$, yielding a speed of sound $c_s\approx 6.1237\times {10}^{-2}$. For example, at a low wavevector $k\approx 0.10$, the free–particle energy is ${ϵ}_k\approx 5.0335\times {10}^{-3}$ and the frequency is $\omega \approx 7.9428\times {10}^{-3}$; at a higher wavevector $k\approx 2.01$, we find ${ϵ}_k\approx 2.0134$ and $\omega \approx 2.0171$. These results clearly illustrate the crossover from linear (phonon–like) behavior at low momenta to the quadratic dispersion characteristic of free particles at higher momenta.

In conclusion, the derivation of the Bogoliubov dispersion relation confirms that gravitational exciton condensates support well–defined collective excitations. The emergence of a linear, sound–like regime at low momenta and the crossover to quadratic (free–particle) behavior at higher momenta provide a robust framework for further investigations into finite–temperature effects, nonlinear dynamics, and possible couplings to gravitational phenomena with potential astrophysical signatures.

7. Conclusion and Discussion

In this work, we have developed a comprehensive theoretical framework for gravitational exciton condensates arising from extra-dimensional stabilization. Starting from a higher-dimensional gravitational action, we derived an effective four-dimensional model in which the conformal moduli of the internal spaces manifest as massive scalar fields (gravitational excitons). By extending the effective action to include self-interacting terms—incorporating higher-order and derivative couplings as well as non-minimal curvature interactions—we provided a natural mechanism for both the stabilization of the extra dimensions and the generation of an effective cosmological constant.

Quantum corrections were introduced via the Coleman–Weinberg mechanism, and gravitational loop effects were incorporated to yield a renormalization-group (RG) improved description of the model. Our analysis revealed that the scale-dependent effective mass and coupling parameters can shift the phase boundaries between stable and unstable regimes, potentially preventing collapse in parameter spaces where the classical theory might otherwise predict instability.

Focusing on the non-relativistic limit, we derived the coupled Gross–Pitaevskii–Poisson system that governs the dynamics of a gravitational exciton Bose–Einstein condensate (BEC). Numerical simulations of this system demonstrated that the condensate relaxes into a stable, low-energy state characterized by a diffuse density profile and a shallow gravitational potential well. The conservation of the mean density during the evolution confirmed particle-number conservation, while the gradual smoothing of density fluctuations and the flattening of the gravitational potential corroborated the expected long-term behavior of a self-gravitating BEC.

Our stability analysis proceeded by seeking stationary solutions via imaginary time evolution. In the spherically symmetric case, we obtained a ground state with a chemical potential $\mu \approx 0.05057$ and an energy $E\approx 0.00587$, with the density peaking at approximately 0.00336 and a corresponding shallow gravitational potential minimum of about $-0.00588$. These results underscore that the condensate, while exhibiting a central density enhancement, is overall diffuse—a feature that is consistent with the expected behavior of dark matter candidates.

Furthermore, our Bogoliubov analysis of small fluctuations about a homogeneous condensate revealed that the excitation spectrum follows the well-known dispersion relation

$${\hslash }^2{\omega }^2={ϵ}_k\left({ϵ}_k+2{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0\right),$$

with the free-particle energy ${ϵ}_k={\hslash }^2{k}^2/\left(2m_{\mathrm{eff}}\right)$. In the long-wavelength limit, the dispersion relation becomes linear, $\omega \approx c_{s}k$, where the speed of sound is given by $c_s=\sqrt{{\lambda }_{\mathrm{eff}}^{\mathrm{NR}}{n}_0/m_{\mathrm{eff}}}$. The numerical results, which yield $c_s\approx 6.1237\times {10}^{-2}$, validate the emergence of phonon-like excitations in the condensate. We further discussed how finite-temperature effects, nonlinear excitations (such as solitons and shock waves), and coupling to gravitational waves might modify the Bogoliubov spectrum, potentially offering novel observational signatures in cosmic microwave background anisotropies and the dynamics of dwarf galaxies.

Overall, our study demonstrates that gravitational exciton condensates possess a rich and robust phenomenology, supported by both analytical derivations and numerical simulations. The interplay between extra-dimensional stabilization, quantum corrections, and nonlinear self-interactions not only provides a viable mechanism for stabilizing the internal dimensions and generating an effective cosmological constant but also suggests that these condensates may serve as promising dark matter candidates. The diffuse density profiles, long lifetimes, and weak interactions of gravitational excitons align well with astrophysical observations, while the predicted Bogoliubov excitations and their modified dispersion relations open up new avenues for indirect detection via gravitational wave observations.

Future work will extend this analysis to non-spherically symmetric and rotating configurations, investigate finite-temperature and damping effects in greater detail, and explore the formation of topological solitons. Such studies will further explain the role of gravitational exciton condensates in the evolution of cosmic structures and contribute to our understanding of the dark sector of the Universe.