Collective Excitations of Self-Gravitating Bose–Einstein Condensates: Breathing Mode and Appearance of Anisotropy under Self-Gravity

Abstract

We investigate the collective mode of a self-gravitating Bose–Einstein condensate (BEC) described by the Gross–Pitaevskii–Poisson (GPP) equations. The self-gravitating BEC has garnered considerable attention in cosmology and astrophysics, being proposed as a plausible candidate for dark matter. Our inquiry delves into the breathing and anisotropic collective modes by numerically solving the GPP equations and using the variational method. The breathing mode demonstrates a reduction in period with increasing total mass due to the density dependence of the self-gravitating BEC, attributed to the density-dependent nature of self-gravitating BECs, aligning quantitatively with our analytical findings. Additionally, we investigate an anisotropic collective mode in which the quadrupole mode intertwines with the breathing mode. The period of the quadrupole mode exhibits similar total mass dependence to that of the breathing mode. The characteristics of these periods differ from those of a conventional BEC confined by an external potential. Despite the differences in density dependence, the ratio of their periods equals that of the BEC confined by an isotropic harmonic potential. Furthermore, an extension of the variational method to a spheroidal configuration enables the isolation of solely the quadrupole mode from the anisotropic collective mode.

1. Introduction

Dark matter (DM) stands as a pivotal topic in the realms of cosmology, astrophysics, and modern physics. The nature and identity of DM have captivated our attention for some decades. The characteristics of DM entail its nonemissive, nonreflective, and nonabsorptive nature towards light while exerting gravitational influence on matter. Although we cannot directly observe DM in any optical way, evidence of its existence can be inferred from the rotational curve of galaxies [1] and gravitational lensing [2,3].

The Λ cold dark matter (ΛCDM) model, based on the assumption of nonrelativistic DM particles, represents a canonical cosmological framework. It elucidates the large-scale structures of the universe and has been used for various cosmological and astrophysical studies [4]. However, it is widely acknowledged that the ΛCDM model has substantial discrepancies at small scales between observational and theoretical results [5,6]. To address this issue, alternate models have been proposed [7], assuming that, e.g. ultralight scalar particles behave as DM [8,9] or nongravitational self-interaction among DM particles exists [10,11].

An intriguing proposition currently under discussion is the possibility that a DM halo composed of ultralight scalar particles may give rise to a phenomenon known as Bose–Einstein condensation [12–14]. Bose–Einstein condensation is a phenomenon in which a macroscopic quantity of Bose particles in a system condenses into a common quantum state. As a result, the physical state of the Bose–Einstein condensate (BEC) is described by the macroscopic wavefunction. This phenomenon occurs when the temperature of the Bose gas falls below the critical threshold necessary for Bose–Einstein condensation, indicating that the de Broglie wavelength of Bose particles exceeds the average interparticle spacing [15].

The challenges inherent in the ΛCDM model could be alleviated through the occurrence of Bose–Einstein condensation. When the mass of a scalar particle ranges from m ∼ 10⁻²³–10⁻²² eV, it is capable of reproducing large-scale structures similar to the ΛCDM model [16,17]. Furthermore, the Bose gas shows a high critical temperature for Bose–Einstein condensation, typically on the order of 10⁹ K [18], owing to the galactic-scale de Broglie wavelength of approximately 1 kpc [13,19]. Consequently, the DM halo comprising such ultralight scalar particles has the potential to become a BEC in the universe. Due to their ultralight mass, scalar particles manifest their quantum nature on a macroscopic scale within this BEC DM halo, and the uncertainty principle serves to suppress the overpopulation of DM at small scales. A promising candidate for such ultralight DM is, e.g. an axion-like particle proposed in the string theory [20].

A BEC of DM is intriguing within the realm of low-temperature physics, as it posits the existence of a galaxy-scale quantum fluid beyond the limitations of laboratory settings. Quantum fluids such as superfluids ⁴He, ³He, and atomic BECs are defined as fluids wherein quantum effects macroscopically manifest owing to extremely low temperatures. Research on these fluids, including collective modes [21], quantized vortices [22], and quantum turbulence [23,24], has been variously conducted in the field of low-temperature physics [15,25,26]. Particularly, a BEC near zero temperature is quantitatively described by the Gross–Pitaevskii (GP) equation, which governs the time evolution of a macroscopic wavefunction [27,28]. Considering a BEC of ultralight DM, it should be bound by the gravitational potential of the BEC itself, and such a self-gravitating BEC follows the Gross–Pitaevskii–Poisson (GPP) equations [17,29,30].

The GPP equations enable DM particles to incorporate a short-range contact interaction. As mentioned above, such self-interaction of DM particles is proposed to resolve the discrepancies of the ΛCDM model. Observational data from colliding clusters of galaxies [31,32] suggest the possibility that DM exhibits nonnegligible nongravitational self-interaction.

A self-gravitating BEC inherently differs from an atomic BEC due to its trapping potential. In the realm of low-temperature physics, a BEC is generally constrained by external potentials, which are predetermined regardless of the configuration or movement of the BEC. Conversely, a self-gravitating BEC is confined by its gravitational potential, even without any external potential. Then, the deformation of the self-gravitating BEC can change the depth, range, and anisotropy of the gravitational potential. These differences between a self-gravitating BEC and an atomic BEC are quite significant when considering the former as a candidate for astrophysical objects, as it causes quantum hydrodynamical phenomena under the influence of self-gravity.

One of the phenomena affected by the anisotropy of the trapping potential is the collective mode. It is defined as a low-frequency oscillation of a trapped BEC in response to small fluctuations. Despite being a many-body system, the collective mode can be described by representative variables as a one-body motion. For example, within a spherical harmonic potential, a BEC can show the collective mode characterized by the radius, known as the monopole (breathing) mode, and another characterized by two semiaxes or their ratio, termed the quadrupole mode. An important characteristic is the dependency of collective modes on the shape of the trapping potential; in an axisymmetric harmonic potential, the breathing and quadrupole modes are coupled. Given the alteration of the gravitational potential with the oscillations of the BEC in self-gravitating systems, the collective mode of such a BEC is likely to exhibit distinct behavior from conventional scenarios.

Various properties of DM halos have been theoretically studied using the GPP equations in recent studies in astrophysics.

The Thomas–Fermi (TF) approximation of the GP equation has proven valuable in describing the equilibrium state of a BEC and carrying out analytical studies. As for a self-gravitating BEC, this approximation was applied for the first time to compare its consequences with observational data on the rotational curves of galaxies [33]. Subsequently, numerous studies employing this approximation have emerged, including the divergence of the central density in the ΛCDM model called the core-cusp problem [34], gravitational lensing effects [35], rotational deformation [36], and the effects of quantized vortices [37,38].

A self-gravitating BEC has also been investigated through the variational method. In low-temperature physics, the calculation of the collective mode of a BEC using this method is renowned and conventional [15,25]. Using the variational method in the study of self-gravitating BECs, the relationship between the total mass and radius of such systems [29], gravitational collapse [39,40], and the phase transition between dilute and dense phases [41] have been investigated.

Numerical investigations into self-gravitating BECs have recently been reported. The analytical solution of the GPP equations for general scenarios is intricate. Hence, numerical calculation of the GPP equations is required to study self-gravitating BECs. Thus far, the stability of the equilibrium state [42], the process of stabilization [43], the rotational velocity of a test particle within a self-gravitating BEC [44,45], and the collisions among self-gravitating BECs [46] were numerically investigated by leveraging spatial symmetries to mitigate computational costs. However, the latest numerical studies have used a 3D system without spatial symmetry, focusing on quantized vortices in self-gravitating BECs, i.e. their stability [47,48] and collisions between two self-gravitating BECs with quantized vortices [49].

In this work, we study the 3D dynamics of the collective modes of a self-gravitating BEC. We assume a sufficiently large total mass to employ the TF approximation. We prepare the initial state by multiplying the factor on the initial velocity with the equilibrium solution of the GPP equations and implement time evolution by numerically solving the GPP equations. Firstly, the breathing mode of the self-gravitating BEC is induced by radially adding the initial velocity. Our numerical results agree with the analytical results obtained using the variational method. Next, by applying the initial velocity axisymmetrically, we numerically obtain the anisotropic collective mode, in which the quadrupole mode is coupled with the breathing mode.

This paper is organized as follows. Our model and numerical setup are described in Sect. 2. Based on this, the breathing mode of a self-gravitating BEC is studied in Sect. 3. The anisotropic collective mode under self-gravity of the BEC is investigated in Sect. 4. Finally, we conclude this paper in Sect. 5.

2. Model and numerical setup

2.1. GPP model and the equilibrium state

We consider a self-gravitating BEC composed of scalar bosons with mass m and an s-wave scattering length a at zero temperature. In the GPP model, the physical state of the BEC is described by the macroscopic wavefunction |$\psi (\boldsymbol {r},t)$|⁠, and the time evolution is given by the following GPP equations defined as

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}i\hbar \partial _t{\psi } = -\frac{\hbar ^2}{2m}\nabla ^2\psi +[mV+\frac{4\pi \hbar ^2a}{m} |\psi |^2]{\psi },\\ {\nabla }^2V= 4\pi mG\vert {\psi }\vert ^2, \end{array}\right. \end{eqnarray}$$

(1)

where |$V(\boldsymbol {r},t)$| denotes the gravitational potential [14,17,29]. Using the Madelung representation written by

$$\begin{eqnarray} \psi (\boldsymbol {r},t)=\sqrt{\frac{\rho (\boldsymbol {r},t)}{m}}e^{i\theta (\boldsymbol {r},t)}, \end{eqnarray}$$

(2)

the GP equation yields the continuity equation and the Euler-like equation, where ρ = m|ψ|² denotes the density and |$\boldsymbol {v}=\hbar \nabla \theta /m$| denotes the velocity field [15]. Then, the total mass is

$$\begin{eqnarray} M=\int d\boldsymbol {r}\rho (\boldsymbol {r})=m\int d\boldsymbol {r}|\psi |^2. \end{eqnarray}$$

(3)

The total energy E is the sum of the kinetic energy K, the potential energy W, and the self-interaction energy I given by [37]:

$$\begin{eqnarray} E&=K+W+I, \end{eqnarray}$$

(4)

$$\begin{eqnarray} K&=\displaystyle\frac{\hbar ^2}{2m}\int d\boldsymbol {r}|\nabla \psi |^2, \end{eqnarray}$$

(5)

$$\begin{eqnarray} W&=\displaystyle\frac{m}{2}\int d\boldsymbol {r}|\psi |^2V, \end{eqnarray}$$

(6)

$$\begin{eqnarray} I&=\displaystyle\frac{2\pi \hbar ^2a}{m}\int d\boldsymbol {r}|\psi |^4. \end{eqnarray}$$

(7)

We can obtain the equilibrium solution of the GPP Eqs. (1) using the TF approximation [33]. When the kinetic energy is negligible, denoted as K ≪ I, the configuration of the equilibrium state is formed by the competition between gravity and self-interaction [15]. The kinetic energy becomes insignificant when the total mass is sufficiently large [50]. The equilibrium configuration is expected to exhibit spherical symmetry under self-gravity. Assuming that the equilibrium state satisfies the ansatz |$\psi (\boldsymbol {r},t)=\phi (r)e^{-i\mu t/\hbar }$|⁠, where μ denotes the chemical potential, the steady solution of Eq. (1) is approximately given by

$$\begin{eqnarray} \rho (r) \simeq \rho _cj_0\left(\frac{\pi r}{R_{\rm TF}}\right)= \rho _c\frac{\sin \left(\pi r/R_{\rm TF}\right)}{\left(\pi r/R_{\rm TF}\right)}~~(r\lt R_{\rm TF}), \end{eqnarray}$$

(8)

where |$r=|\boldsymbol {r}|$| represents the radial coordinate, ρ_c stands for the central density, and j₀ is the 0-th spherical Bessel function. Subsequently, the TF radius R_TF is defined as the minimum radius at which the density becomes zero:

$$\begin{eqnarray} R_{\rm TF}=\pi \sqrt{\frac{\hbar ^2a}{Gm^3}}, \end{eqnarray}$$

(9)

indicating that the size of a massive self-gravitating BEC remains independent of its total mass. There are no densities in the range of r > R_TF. Thus, the cumulative mass profile is expressed as

$$\begin{eqnarray} M\left(r\right)=\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{M_{\rm TF}}{\pi }\left\lbrace \sin \left(\frac{\pi r}{R_{\rm TF}}\right)-\frac{\pi r}{R_{\rm TF}}\cos \left(\frac{\pi r}{R_{\rm TF}}\right)\right\rbrace & (r\lt R_{\rm TF}), \\ M_{\rm TF} & (r\gt R_{\rm TF}), \end{array}\right. \end{eqnarray}$$

(10)

where M_TF is defined as

$$\begin{eqnarray} M_{\rm TF}=\frac{4\rho _cR_{\rm TF}^3}{\pi }. \end{eqnarray}$$

(11)

Given that the TF radius remains unaffected by the total mass, Eq. (11) reveals that only the central density escalates, even with an increase in the total mass. By solving the Poisson equation, the gravitational potential becomes

$$\begin{eqnarray} V(r) =G\int ^r_\infty ds\frac{M(s)}{s^2} =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}-\displaystyle\frac{GM_{\rm TF}}{R_{\rm TF}}\left\lbrace 1+j_0\left(\frac{\pi r}{R_{\rm TF}}\right)\right\rbrace & (r\lt R_{\rm TF}),\\ -\displaystyle\frac{GM_{\rm TF}}{r} & (r\gt R_{\rm TF}). \end{array}\right. \end{eqnarray}$$

(12)

Therefore, using Eqs. (5), (6), (7), (8), and (12), the energy of each equilibrium state can be expressed as

$$\begin{eqnarray} &K\simeq \displaystyle\frac{\pi }{8}\frac{\hbar ^2M_{\rm TF}}{m^2R_{\rm TF}^2}\left\lbrace \mathrm{Si}(\pi )-\pi +\displaystyle\int _0^{\Lambda \pi }dx\displaystyle\frac{x}{\sin x}\right\rbrace , \end{eqnarray}$$

(13)

$$\begin{eqnarray} W=-\frac{3}{4}\displaystyle\frac{GM_{\rm TF}^2}{R_{\rm TF}}, \end{eqnarray}$$

(14)

$$\begin{eqnarray} &I=\displaystyle\frac{\pi ^2}{4}\displaystyle\frac{\hbar ^2aM_{\rm TF}^2}{m^3R_{\rm TF}^3}, \end{eqnarray}$$

(15)

where Si(x) is a sine integral and Λπ(< π) serves as the cutoff to avoid the divergence of the integral for Λ → 1. The divergence is derived from the breakdown of the TF approximation, and Eq. (13) shows the kinetic energy within a spherical region whose radius is obtained by subtracting the depth of the BEC surface from the TF radius. The value of Λ can be estimated using the ratio between the coherence length |$\xi =\sqrt{m/(8\pi a\rho _c)}$| and the TF radius, such that Λ ∼ 1 − ξ/R_TF.

2.2. Numerical setup

The GPP Eqs. (1) are expressed in the following dimensionless formulation:

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}i\tilde{\partial }_t\tilde{\psi }=-\frac{1}{2}\tilde{\nabla }^2\tilde{\psi }+[\tilde{V}-i\tilde{V}_{\rm s}+\tilde{a}|\tilde{\psi }|^2]\tilde{\psi }, \\ \tilde{\nabla }^2\tilde{V}=|\tilde{\psi }|^2, \end{array}\right. \end{eqnarray}$$

(16)

where the tilde symbols represent the dimensionless variables, and as explained later, the term |$-i\tilde{V}_{\rm s}$| in the dimensionless GP equation is artificially added owing to a computational reason. The dimensionless variables are defined as follows: |$\tilde{\boldsymbol {r}}=\boldsymbol {r}(\tilde{\lambda } mc/\hbar ), \tilde{t}=t(\tilde{\lambda }^2mc^2/\hbar ), \tilde{\psi }=\psi \lbrace \sqrt{4\pi G}\hbar /(\sqrt{m}c^2\tilde{\lambda }^2)\rbrace$|⁠, |$\tilde{V}=V/(\tilde{\lambda }^2c^2), \tilde{a}=a\lbrace \tilde{\lambda }^2c^2/(mG)\rbrace$| [44]. The scaling factor |$\tilde{\lambda }$| can adjust the system size using the invariance of the equations for any value.¹ Then, the dimensionless total mass and energy are given by |$\tilde{M}=M\lbrace 4\pi Gm/(\tilde{\lambda }\hbar c)\rbrace$| and |$\tilde{E}=E\lbrace 8\pi Gm/(\hbar c^3\tilde{\lambda }^3)\rbrace$|⁠.

The boundary conditions should be taken into account when dealing with Eq. (16). The GP equation is typically solved under periodic boundary conditions. However, the Poisson equation is inconsistent with the periodic boundary condition owing to the long-range gravitational interaction. To overcome this inconsistency, we propose a novel method.

The GP equation in Eq. (16) is computed using the fourth-order Runge–Kutta and pseudo-spectrum methods under periodic boundary conditions. To solve the Poisson equation, we use the finite difference and Jacobi methods at each time step [51]. Using the Jacobi method, the solution of the Poisson equation in Eq. (16) manifests as the equilibrium solution of the diffusion equation for the function |$\tilde{\mathrm{ V}}^{\prime }(\tilde{\boldsymbol {r}},\tilde{t};\tilde{s})$|⁠:

$$\begin{eqnarray} \frac{d\tilde{\mathrm{ V}}^{\prime }}{d\tilde{s}}(\tilde{\boldsymbol {r}},\tilde{t};\tilde{s}) = \tilde{\nabla }^2\tilde{\mathrm{ V}}^{\prime }(\tilde{\boldsymbol {r}},\tilde{t};\tilde{s}) - |\tilde{\psi }(\tilde{\boldsymbol {r}},\tilde{t})|^2, \end{eqnarray}$$

(17)

where |$\tilde{s}$| denotes a real parameter. The gravitational potential is obtained from |$\tilde{\mathrm{ V}}(\tilde{\boldsymbol {r}},\tilde{t})=\lim _{\tilde{s}\rightarrow \infty }\tilde{\mathrm{ V}^{\prime }}(\tilde{\boldsymbol {r}},\tilde{t};\tilde{s})$|⁠. Furthermore, the efficacy of the Jacobi method can be notably enhanced through Nesterov’s accelerated gradient method [52–54]. This method is common in the field of machine learning and enables us to rapidly converge the function |$\tilde{\mathrm{ V}^{\prime }}(\tilde{\boldsymbol {r}},\tilde{t};\tilde{s})$| towards the gravitational potential |$\tilde{\mathrm{ V}}(\tilde{\boldsymbol {r}},\tilde{t})$|⁠. The function |$\tilde{\mathrm{ V}^{\prime }}(\tilde{\boldsymbol {r}},\tilde{t};\tilde{s})$| satisfies the Dirichlet boundary condition given by

$$\begin{eqnarray} \tilde{\mathrm{ V}}^{\prime }(\tilde{\boldsymbol {r}}=\tilde{\boldsymbol {R}}_b,\tilde{t};\tilde{s}) = -\frac{\tilde{M}(\tilde{t})}{4\pi |\tilde{\boldsymbol {R}}_b-\tilde{\boldsymbol {R}_c}|} , \end{eqnarray}$$

(18)

where |$\tilde{\boldsymbol {R}}_b$| denotes the position on the boundary of the numerical domain and |$\tilde{\boldsymbol {R}_c}$| denotes the center of the numerical domain. We apply a time-dependent total mass |$\tilde{M}(\tilde{t})$| because, as mentioned later, we consider the probability of emission of particles. The boundary condition (18) is appropriate when the BEC occupies a region smaller than the numerical domain and is located near its center.

In our numerical model, the imaginary potential |$-i\tilde{\mathrm{ V}}_{\rm s}$|⁠, called “sponge” potential, written by

$$\begin{eqnarray} \tilde{\mathrm{ V}}_{\rm s}(\boldsymbol {r}) = \frac{\tilde{\mathrm{ V}}_o}{2}\left\lbrace 2+\tanh \left(\frac{\tilde{r}-\tilde{r}_c}{\tilde{\delta }}\right)-\tanh \left(\frac{\tilde{r}_c}{\tilde{\delta }}\right)\right\rbrace , \end{eqnarray}$$

(19)

is introduced into Eq. (16) [44]. Here, |$\tilde{\mathrm{ V}}_o$| denotes the amplitude. This sponge potential can reduce the number of particles within the range of |$\tilde{r}\gt \tilde{r}_c$|⁠, where |$\tilde{\delta }$| denotes the step width. Functioning as a sink, the sponge potential aids in removing particles emitted from the BEC. As the BEC undergoes motion, some particles acquire large kinetic energy, leading to their escape from the gravitational potential. However, in our system, they return to the BEC due to periodic boundary conditions on the GP equation. To avoid the unphysical situation, the sponge potential becomes imperative.

The initial state is prepared by multiplying the initial phase factor |$\exp [i\tilde{\phi }(\tilde{\boldsymbol {r}})]$| and the equilibrium state |$\tilde{\psi }_{\rm eq}(\tilde{\boldsymbol {r}})$|⁠:

$$\begin{eqnarray} \tilde{\psi }(\tilde{\boldsymbol {r}},\tilde{t}=0) = \tilde{\psi }_{\rm eq}(\tilde{\boldsymbol {r}})e^{i\tilde{\phi }(\tilde{\boldsymbol {r}})}. \end{eqnarray}$$

(20)

The BEC is located at the center of the numerical domain. We obtain the equilibrium state |$\tilde{\psi }_{\rm eq}$| using the imaginary time propagation method of the GPP equations. The initial phase is given by

$$\begin{eqnarray} \tilde{\phi }(\tilde{r}_\perp ,\tilde{z})=\frac{1}{2}\left(\tilde{\alpha } \tilde{r}_\perp ^2+\tilde{\beta } \tilde{z}^2\right), \end{eqnarray}$$

(21)

which gives rise to the initial velocity field |$\tilde{\boldsymbol {v}}=\tilde{\nabla }\tilde{\phi }=(\tilde{\alpha } \tilde{r}_\perp )\boldsymbol {e}_\perp +(\tilde{\beta } \tilde{z})\boldsymbol {e}_z$|⁠. Here, |$\tilde{r}_\perp =\sqrt{\tilde{x}^2+\tilde{y}^2}$| and |$\tilde{z}$| denote the 3D cylindrical coordinates, and the unit vectors along each direction are |$\boldsymbol {e}_\perp$| and |$\boldsymbol {e}_z$|⁠. We determine the initial velocity by |$\tilde{\alpha }$| and |$\tilde{\beta }$| (e.g. it becomes isotropic when |$\tilde{\alpha }=\tilde{\beta }$|⁠).

Finally, the physical and numerical parameters are specified. We consider the case in which each boson has a mass m = 3 × 10⁻²⁴ eV and an s-wave scattering length a = 5.62 × 10⁻⁹⁸λ_c m, where λ_c denotes the Compton wavelength [47]. For the numerical analyses, we use a 3D cube of length |$\tilde{L}_x=\tilde{L}_y=\tilde{L}_z=40$|⁠. The scaling factor is |$\tilde{\lambda }=3\times 10^{-3}$|⁠, ensuring that the numerical domain vastly exceeds the size of the BEC. The spatial grid is configured as |$\tilde{N}_x=\tilde{N}_y=\tilde{N}_z=128$|⁠. The time resolution of the GP equation is |$\Delta \tilde{t}=10^{-3}$|⁠, and the resolution of the parameter |$\tilde{s}$| is |$0.05\times (\tilde{L}_x/\tilde{N}_x)^2\approx 4.88\times 10^{-3}$|⁠. The parameters of the sponge potential are designated as |$\tilde{V}_o=1000, \tilde{r}_c=0.8\times (\tilde{L}_x/2)=16$|⁠, and |$\tilde{\delta }=2\times (\tilde{L}_x/\tilde{N}_x)=0.625$|⁠. Then, |$\tilde{r}_c$| exceeds the dimensionless TF radius |$\tilde{R}_{\rm TF}\approx 9.09$|⁠, which represents the typical size of a BEC. This is reasonable because the sponge potential works outside the BEC. In all our simulations, the change rate of the total mass from its initial value is less than 1%, and particles are hardly emitted from the BEC.

3. Breathing mode

In this section, we investigate the breathing modes of a self-gravitating BEC. For a conventional BEC confined by an external potential, the initial radial velocity field induces the breathing mode [15,25,55]. Likewise, we analytically and numerically investigate the breathing mode of a self-gravitating BEC by setting a spherically symmetric initial phase. We validate our numerical results through comparison with our analytical predictions.

3.1. Analytical calculation by variational method

We apply the variational method in the GPP model to reveal the breathing modes of a self-gravitating BEC. We examine a self-gravitating BEC with a fixed total mass M. Previous studies have used Gaussian [29] and quadratic functions [39] as trial functions of the variational method. We adopt the 0-th spherical Bessel function as the trial function because we posit that the self-gravitating BEC possesses a large total mass, allowing us to use the TF approximation. This approach is the most suitable for such a large total mass.

We prepare a trial function for the variational method. When the self-gravitating BEC is so large that the TF approximation is available, the macroscopic wavefunction of the equilibrium state is derived as

$$\begin{eqnarray} \psi _{\rm eq}(r)=\sqrt{\frac{\pi M}{4mR_{\rm TF}^3}j_0\left(\frac{\pi r}{R_{\rm TF}}\right)} \end{eqnarray}$$

(22)

using Eqs. (8) and (11). Assuming that the density profile maintains its form as in Eq. (22) during the motion of the cloud, the trial function can be obtained as

$$\begin{eqnarray} \psi (\boldsymbol {r},t) = \sqrt{ \frac{\pi M}{4mR(t)^3} j_0\left(\frac{\pi r}{R(t)}\right) } \exp \left[i\frac{mr^2}{2\hbar }H(t)\right] \end{eqnarray}$$

(23)

in the range of r < R(t) [15,39,40], and the radius R(t) and H(t) are treated as variables governing the breathing mode.

We derive the Euler–Lagrange equations for R(t) and H(t). In the GPP model, the Lagrangian can be described by

$$\begin{eqnarray} L=\frac{i\hbar }{2}\int d\boldsymbol {r}\left(\psi ^{*}\partial _t\psi -\psi \partial _t\psi ^{*}\right)-E, \end{eqnarray}$$

(24)

using Eq. (4) [39]. Upon substituting the trial function (23), the Lagrangian transforms into

$$\begin{eqnarray} L(R,H)=-\frac{1}{2}\left(1-\frac{6}{\pi ^2}\right)MR^2(\dot{H}+H^2)-U(R). \end{eqnarray}$$

(25)

Here, the effective potential U(R) is defined as

$$\begin{eqnarray} U(R) = \frac{\pi }{8}\frac{\hbar ^2M}{m^2R^2}\mathrm{ F}(\Lambda )-\frac{3}{4}\frac{GM^2}{R}+\frac{\pi ^2}{4}\frac{\hbar ^2aM^2}{m^3R^3} \equiv \frac{C_z}{R^2}-\frac{C_p}{R}+\frac{C_i}{R^3}. \end{eqnarray}$$

(26)

Here the coefficients are |$C_z=\frac{\pi }{8}\frac{\hbar ^2M}{m^2}\mathrm{ F}(\Lambda )$|⁠, |$C_p=\frac{3}{4}GM^2$|⁠, and |$C_i=\frac{\pi ^2}{4}\frac{\hbar ^2aM^2}{m^3}$|⁠, and the function |$\mathrm{ F}(\Lambda ) =\mathrm{Si}(\pi )-\pi +\int _0^{\Lambda \pi }dx\frac{x}{\sin x}$|⁠. Consequently, the Euler–Lagrange equation for H(t) is |$H(t)=\dot{R}(t)/R(t)$|⁠, and the radius R(t) follows the equation² given by

$$\begin{eqnarray} \left(1-\frac{6}{\pi ^2}\right)M\frac{d^2R(t)}{dt^2}=-\frac{dU(R)}{dR}. \end{eqnarray}$$

(27)

The equilibrium solution R_eq is acquired by setting dU(R)/dR = 0, yielding

$$\begin{eqnarray} R_{\rm eq} &=&\frac{C_z+\sqrt{C_z^2+3C_pC_i}}{C_p} \\ &=&\frac{\pi ^2\hbar ^2}{6Gm^2}\frac{\mathrm{ F}(\Lambda )}{M}\left\lbrace 1+\sqrt{1+\frac{36Gam}{\pi ^2\hbar ^2}\frac{M^2}{\mathrm{F}(\Lambda )^2}}\right\rbrace . \end{eqnarray}$$

(28)

Figure 1 describes the result of Eq. (28) with Λ = 0.97 as the cutoff parameter. This figure shows that as the total mass M increases, R_eq monotonically decreases due to gravity and converges to the TF radius defined in Eq. (9) as M → ∞. Thus, our results offer an enhanced understanding of the relationship, particularly for situations with large total masses, compared to previous studies [29,39] employing other trial functions.

Fig. 1.

The relationship between the total mass M and the radius R derived from Eq. (9) with the cutoff parameter set to Λ = 0.97. The horizontal axis shows the total mass, and the vertical axis shows the radius. The cyan solid line shows the analytical result of Eq. (9), and the red dotted line shows the TF radius R_TF.

The breathing mode manifests as a small oscillation of radius around R_eq. The radius is expressed as R(t) = R_eq + δR(t), using its fluctuation δR(t) (|δR(t)| ≪ R_eq). Employing Eqs. (26), (27), and (28), we obtain the following equation of motion:

$$\begin{eqnarray} \frac{d^2\delta {R}(t)}{dt^2} \simeq -\frac{\pi ^2}{\pi ^2-6}\frac{1}{M}\frac{d^2U(R_{\rm eq})}{dR^2}\delta R(t) =-\frac{2\pi ^2}{\pi ^2-6}\frac{C_zR_{\rm eq}+3C_i}{MR_{\rm eq}^5}\delta R(t). \end{eqnarray}$$

(29)

In this context, we assume the total mass of the BEC to be sufficiently large to warrant the TF approximation. Given the negligible contribution of C_z, Eq. (29) is simplified to

$$\begin{eqnarray} \frac{d^2\delta {R}(t)}{dt^2}=-\frac{3}{2\pi (\pi ^2-6)}\sqrt{\frac{G^5m^9}{\hbar ^6a^3}}M\delta R(t), \end{eqnarray}$$

(30)

which implies a harmonic oscillation with a period

$$\begin{eqnarray} T_{\rm B}=2\pi \sqrt{\frac{2\pi (\pi ^2-6)}{3}\sqrt{\frac{\hbar ^6a^3}{G^5m^9}}}\frac{1}{\sqrt{M}}. \end{eqnarray}$$

(31)

The dependence on M has also been delineated in previous studies [29,39]. The prefactor and M-dependence can be confirmed by the following numerical calculations.

3.2. Numerical results

To confirm Eq. (31), we conduct numerical simulations for the total mass M/(10¹⁴M_⊙) = 1, 2, 3, 4, where M_⊙ denotes a solar mass. These values significantly exceed the typical mass of a galaxy. For example, the Andromeda galaxy possesses M ∼ 10¹²M_⊙ [37], and the Milky Way is estimated at M ≃ 0.5–2.0 × 10¹²M_⊙ [56]. However, the substantial total mass facilitates the exploration of BEC oscillations because R_eq hardly changes even if the total mass is reduced (see Fig. 1). For a BEC with the total mass of a typical galaxy, namely M ∼ 10¹⁰–10¹²M_⊙, the timescale of the breathing mode would extend to 10–100 times greater than that presented in this paper, based on Eq. (31). We set the parameters of the initial phase to |$\tilde{\alpha }=\tilde{\beta }=\pm 0.03, 0.1$| to give an isotropic velocity field to a BEC in equilibrium.

Initially, we describe the time evolution of the density and gravitational potential profiles. Despite solving the fully 3D GPP equations in the simulations presented in this section, the configurations of the BECs maintain spherical symmetry throughout the time evolution. Then, we consider the average of the density and gravitational potential over the solid angle to describe the time evolution, expressed as

$$\begin{eqnarray} \rho (r,t) = \frac{1}{4\pi r^2\Delta r} \int _{r-\frac{\Delta r}{2}}^{r+\frac{\Delta r}{2}}dr^{\prime }r^{\prime 2} \int d\Omega \rho (\boldsymbol {r^{\prime }},t) \end{eqnarray}$$

(32)

and

$$\begin{eqnarray} \mathrm{ V}(r,t) = \frac{1}{4\pi r^2\Delta r} \int _{r-\frac{\Delta r}{2}}^{r+\frac{\Delta r}{2}}dr^{\prime }r^{\prime 2} \int d\Omega \mathrm{ V}(\boldsymbol {r^{\prime }},t), \end{eqnarray}$$

(33)

where Ω denotes a solid angle and Δr ≈ 0.22 kpc is set so that |$\Delta \tilde{r}$| matches the spatial resolution, namely |$\Delta \tilde{r}=\tilde{L}_x/\tilde{N}_x$|= 0.3125. Figure 2 shows the time evolution of the density and gravitational potential profiles. The BEC periodically undergoes isotropic expansion and shrinkage. Initially, the BEC expands by an outwardly directed initial velocity field, and the central density decreases until t = 0.46 Myr. Subsequently, the BEC shrinks, and the central density increases until t = 1.39 Myr. Finally, the BEC returns to its initial state at t = 1.86 Myr. The gravitational potential also undergoes simultaneous oscillations with the density profile. The gravitational potential becomes shallower as the central density decreases, and vice versa. Note that at the outskirts, namely for r ≳ R_TF ≈ 6.46 kpc, the gravitational potential remains unchanged throughout the time evolution, maintaining V = −GM/r due to the constancy of the total mass.

Fig. 2.

The density profile and the gravitational potential of the self-gravitating BEC. The total mass is M = 2 × 10¹⁴M_⊙, and the initial phase is set by |$\tilde{\alpha }=\tilde{\beta }=0.1$|⁠. The radial coordinate, representing the distance from the center of our numerical domain, is depicted on the horizontal axis. Each column portrays either the density profile or the gravitational potential at each time: t = 0, 0.46, 1.39, 1.86 Myr, where Myr means a megayear. The magenta points show the data at each time, whereas the black dashed line shows the initial data. The upper row displays the density profile, with the density depicted on the vertical axis, whereas the lower row exhibits the gravitational potential, with its value presented on the vertical axis.

To quantitatively observe the spherical oscillation, we define the effective radius R_eff as

$$\begin{eqnarray} R_{\text {eff}}(t) = \sqrt{ \frac{1}{M}\int d\boldsymbol {r} r^2\rho (\boldsymbol {r},t) }. \end{eqnarray}$$

(34)

Using Eq. (23), it is proportional to the radius R(t). Hence, during the occurrence of the breathing mode in the self-gravitating BEC, the effective radius exhibits an oscillation similar to that of the radius. Figure 3(a) shows the deviation of the effective radius from its initial value, ΔR_eff(t) = R_eff(t) − R_eff(0), for |$\tilde{\alpha }=0.03$|⁠. The graph shows sinusoidal curves for each total mass, where the oscillation period |$T_{R_{\text {eff}}}$| decreases as the total mass M is increased. Figure 3(b) suggests that our results of |$T_{R_{\text {eff}}}$| agree quantitatively with Eq. (31), even for three values of |$\tilde{\alpha }$|⁠. Therefore, we conclude that the self-gravitating BEC induces the breathing mode.

Fig. 3.

Oscillation of the effective radius R_eff(t) for each total mass of the self-gravitating BEC. (a) The deviation of the effective radius from the initial/equilibrium one ΔR_eff(t) = R_eff(t) − R_eff(0) when the initial phase is |$\tilde{\alpha }=0.03$|⁠. The horizontal axis shows the time, and the vertical axis shows ΔR_eff(t). Each solid curve is the time evolution of ΔR_eff for different total masses (red: 4 × 10¹⁴M_⊙, blue: 3 × 10¹⁴M_⊙, green: 2 × 10¹⁴M_⊙, and black: 1 × 10¹⁴M_⊙). The black dotted line shows the oscillation criterion, namely ΔR_eff = 0. (b) Total mass dependence on the period of the effective radius |$T_{R_{\text {eff}}}$| for various values of |$\tilde{\alpha }$|⁠. The horizontal axis shows the total mass, and the vertical axis shows the negative squared period |$1/T_{R_{\text {eff}}}^2$|⁠. Each point shows a numerical result (black square: |$\tilde{\alpha }=0.1$|⁠, red solid circle: |$\tilde{\alpha }=0.03$|⁠, and green triangle: |$\tilde{\alpha }=-0.03$|⁠). The blue dashed line shows the analytical result obtained from the variational method, namely the negative squared period of breathing mode |$1/T_{\rm B}^2$|⁠.

4. Appearance of anisotropic mode

4.1. Numerical results

We contemplate the axisymmetric collective mode of a self-gravitating BEC. Its total mass is M/(10¹⁴M_⊙) = 1, 2, 3, 4, the same as that in the breathing mode. These situations can be simulated by configuring |$(\tilde{\alpha },\tilde{\beta })=(\pm 0.03,0),~(0,\pm 0.03)$|⁠. Such parameters introduce an axisymmetric initial velocity field to the BEC in equilibrium and can cause anisotropic oscillation.

To describe the time evolution of the BEC, we obtain the density and gravitational potential profiles by averaging |$\rho (\boldsymbol {r},t)$| and |$\mathrm{ V}(\boldsymbol {r},t)$|⁠:

$$\begin{eqnarray} \rho (r_\perp ,z,t) = \frac{1}{2\pi r_\perp \Delta r_\perp } \int _{r_\perp -\frac{\Delta r_\perp }{2}}^{r_\perp +\frac{\Delta r_\perp }{2}}dr^{\prime }_\perp r^{\prime }_\perp \int _0^{2\pi }d\theta \rho (r^{\prime }_\perp ,\theta ,z,t) \end{eqnarray}$$

(35)

and

$$\begin{eqnarray} \mathrm{ V}(r_\perp ,z,t) = \frac{1}{2\pi r_\perp \Delta r_\perp } \int _{r_\perp -\frac{\Delta r_\perp }{2}}^{r_\perp +\frac{\Delta r_\perp }{2}}dr^{\prime }_\perp r^{\prime }_\perp \int _0^{2\pi }d\theta \mathrm{ V}(r^{\prime }_\perp ,\theta ,z,t), \end{eqnarray}$$

(36)

where Δr_⊥ ≈ 0.22 kpc is determined in the same manner as Δr, namely |$\Delta \tilde{r}_\perp =\tilde{L}_x/\tilde{N}_x=0.3125$|⁠, and θ = tan ⁻¹(y/x). Figure 4 shows the density profile and the gravitational potential. They show synchronous oscillations similar to those in the breathing mode. However, ρ(r_⊥, z = 0, t) is different from ρ(r_⊥ = 0, z, t). At t = 0.73 Myr, the former shows BEC expansion, whereas the latter shows shrinkage. Subsequently, the density profile at t = 2.20 Myr exhibits an inverse behavior: a decrease in ρ(r_⊥, z = 0, t) and an increase in ρ(r_⊥ = 0, z, t). By t = 2.94 Myr, the density profile reverts to a spherical shape similar to the initial profile. Hence, the self-gravitating BEC induces anisotropic oscillations due to the nonspherical initial velocity field. We observe this behavior for other combinations of |$(\tilde{\alpha },\tilde{\beta })$| such as |$\tilde{\alpha }\ne \tilde{\beta }$| (e.g. |$(\tilde{\alpha },\tilde{\beta })=(\pm 0.1,\mp 0.03)$|⁠). Note that the gravitational potential remains in the shape of 1/r at the outskirts throughout our simulations although the density profile anisotropically deforms.

Fig. 4.

The time evolution of the density profile ρ(r_⊥, z) and the gravitational potential V(r_⊥, z). We consider the BEC has a total mass M = 2 × 10¹⁴M_⊙ and is endowed with an initial velocity by |$\tilde{\alpha }=0.1,\tilde{\beta }=0$|⁠. The horizontal axis shows r_⊥ and z. Magenta points show the profiles at z = 0, ρ(r_⊥, z = 0), and V(r_⊥, z = 0), whereas cyan triangles show those at r_⊥ = 0, ρ(r_⊥ = 0, z), and V(r_⊥ = 0, z). The black dashed lines show the initial profile. Each column shows either the density profile or the gravitational potential at each time: t = 0, 0.73, 2.20, 2.94 Myr, where Myr means a megayear. The upper row shows the density profile, and the vertical axis shows the density. The lower row shows the gravitational potential, and the vertical axis shows its value.

To extract the effective degrees of freedom from the anisotropic oscillation, we introduce two quantities |$R_{\text {eff}}^\perp$| and Z_eff defined as

$$\begin{eqnarray} R_{\text {eff}}^\perp (t) = \sqrt{ \frac{1}{M}\int d\boldsymbol {r} r_\perp ^2\rho (\boldsymbol {r},t) }, \end{eqnarray}$$

(37)

$$\begin{eqnarray} Z_{\text {eff}}(t) = \sqrt{ \frac{1}{M}\int d\boldsymbol {r} z^2\rho (\boldsymbol {r},t) }. \end{eqnarray}$$

(38)

These are the effective widths of the BEC in each direction; |$R_{\text {eff}}^\perp$| characterizes the width in the direction of the radial coordinate r_⊥ and Z_eff is the width along the z-axis. Figure 5 shows the complicated oscillatory behavior of these quantities, describing the opposing oscillations between |$R_{\text {eff}}^\perp$| and Z_eff. The BEC exhibits shrinkage along the z-axis when a positive velocity is initially added in the direction of r_⊥, and vice versa. Similarly, the BEC shows shrinkage in the direction of r_⊥ in response to a positive initial velocity in the direction of the z-axis, and vice versa.

Fig. 5.

Oscillations of the effective width |$R_{\text {eff}}^\perp ,~Z_{\text {eff}}$|⁠. The self-gravitating BEC possesses a total mass M = 4 × 10¹⁴M_⊙, and the initial velocity is set by |$(\tilde{\alpha },~\tilde{\beta })=(\pm 0.03,0),~(0,\pm 0.03)$|⁠. The horizontal axis shows the time. The left column shows the result of |$(\tilde{\alpha },~\tilde{\beta })=(\pm 0.03,0)$|⁠, which means the initial velocity in the direction of r_⊥. The green solid line shows the oscillation when the BEC initially expands (⁠|$\tilde{\alpha }\gt 0$|⁠), and the black dashed line shows the oscillation when the BEC initially shrinks (⁠|$\tilde{\alpha }\lt 0$|⁠). Conversely, the right column shows the result of |$(\tilde{\alpha },~\tilde{\beta })=(0,\pm 0.03)$|⁠, indicating the initial velocity along the z-axis. The red solid line shows the oscillation when the BEC initially expands (⁠|$\tilde{\beta }\gt 0$|⁠), and the blue dashed line shows the oscillation when the BEC initially shrinks (⁠|$\tilde{\beta }\lt 0$|⁠). Whereas the upper row shows the time evolution of |$R_{\text {eff}}^\perp$|⁠, the lower row shows the time evolution of Z_eff.

Although |$R_{\text {eff}}^\perp$| and Z_eff manifest intricate oscillations, a harmonic oscillation can be obtained from the ratio |$Z_{\text {eff}}/R_{\text {eff}}^\perp$|⁠, as shown in Fig. 6. Among the four cases, the oscillations have a common amplitude and frequency despite the different behaviors of |$R_{\text {eff}}^\perp ,~Z_{\text {eff}}$|⁠. The ratio |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| indicates the extent to which the BEC axisymmetrically deforms from a spherical shape, similar to the aspect ratio of the BEC. It can be inferred that one of the coupled modes is the quadrupole mode because the harmonic mode of the aspect ratio suggests a quadrupole mode of the BEC [57–59].

Fig. 6.

The oscillation of the ratio |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| when the total mass of the self-gravitating BEC is M = 4 × 10¹⁴M_⊙. The horizontal axis shows the time, and the vertical axis shows the value of the ratio |$Z_{\text {eff}}/R_{\text {eff}}^\perp$|⁠. The solid lines show cases of initial expansion (green: the direction of r_⊥ and red: along the z-axis), whereas the dotted lines exhibit cases of initial shrinkage (black: the direction of r_⊥ and blue: the direction of the z-axis).

The Fourier transformation serves as a valuable tool for decomposing complicated oscillations into eigenmodes. Figure 7 shows the Fourier transformation of |$R_{\text {eff}}^\perp$| and |$Z_{\text {eff}}/R_{\text {eff}}^\perp$|⁠, as shown in Figs. 5 and 6. Within this representation, two peaks of |$R_{\text {eff}}^\perp$| and a peak of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| can be discerned. The high-frequency peak of |$R_{\text {eff}}^\perp$| corresponds to the period of the breathing mode T_B with the corresponding total mass. On the other hand, the low-frequency peak in |$R_{\text {eff}}^\perp$| closely matches that of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$|⁠. This peak suggests the frequency associated with the quadrupole mode of the BEC. Hence, the self-gravitating BEC causes an anisotropic collective mode wherein the quadruple mode superposes the breathing mode due to the axisymmetric initial velocity field.

Fig. 7.

The Fourier transformation of |$R_{\text {eff}}^\perp$| and |$Z_{\text {eff}}/R_{\text {eff}}^\perp$|⁠. The self-gravitating BEC has a total mass M = 4 × 10¹⁴M_⊙, with the initial velocity given by |$\tilde{\alpha }=0.03,~\tilde{\beta }=0$|⁠. The horizontal axis shows the frequency, whereas the vertical axis shows the amplitude. The solid lines show the results of each transformation (red: |$R_{\text {eff}}^\perp$| and blue: |$Z_{\text {eff}}/R_{\text {eff}}^\perp$|⁠), and the green dashed line shows the analytical/variational result of the breathing mode T_B. The red points show the peaks of |$R_{\text {eff}}^\perp$|⁠, and the blue triangle shows the peak of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$|⁠.

We vary the parameters of the initial phase and total mass to investigate the dependence of the period of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| on these quantities. Although we study the dynamics of various combinations of |$\tilde{\alpha }$| and |$\tilde{\beta }$|⁠, the period remains constant for all the simulations. Thus, the period of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| is independent of the initial velocity field. However, the total mass can change the period of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$|⁠, as shown in Fig. 8. Figure 8(a) shows that |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| exhibits harmonic oscillation similar to R_eff of the breathing mode induced by the isotropic initial velocity field; the period becomes short as the total mass increases. The M-dependence of the period of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| and its comparison with T_B are shown in Fig. 8(b). We find that |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| shows a straight line parallel to T_B in the log-log plot, indicating that the period of the quadrupole mode is also inversely proportional to |$\sqrt{M}$|⁠, similarly to the breathing mode. However, the period of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| is approximately 1.57 times larger than T_B; in other words, the quadrupole mode exhibits lower frequency oscillation compared to the breathing mode. This characteristic can also be observed in atomic BECs confined by an isotropic harmonic potential. The frequency of the breathing mode is |$\sqrt{5/2}$| times that of the quadrupole mode [25]. Hence, our numerical results are consistent with those of typical atomic BECs trapped by external potentials. The ratio of the period of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| to T_B closely resembles that of the BEC confined by an isotropic harmonic potential. Moreover, the periods of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| are in good agreement with |$\sqrt{5/2}T_B$| (see Fig. 8). Thus, the numerical results indicate that the quadrupole mode has the period |$\sqrt{5/2}T_B$|⁠.

Fig. 8.

The oscillation of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| for each total mass. The initial phase is determined by |$\tilde{\alpha }=0.03, \tilde{\beta }=0$|⁠. (a) The time evolution of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| for each total mass. The horizontal axis shows time, whereas the vertical axis shows the value of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$|⁠. The black dotted line shows the ratio when the BEC is in equilibrium; the ratio becomes |$Z_{\text {eff}}/R_{\text {eff}}^\perp =1/\sqrt{2}$|⁠. Each solid line shows the result (red: 4 × 10¹⁴M_⊙, blue: 3 × 10¹⁴M_⊙, green: 2 × 10¹⁴M_⊙, black: 1 × 10¹⁴M_⊙). (b) The dependence of the period of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| on the total mass by the log-log plot. The horizontal axis shows the total mass, whereas the vertical axis shows the period. The black points show the numerical results for each total mass, whereas the blue dashed line shows the variational result of the breathing mode T_B. The red solid line shows |$\sqrt{5/2}T_B$|⁠.

The period of the quadrupole mode in a self-gravitating BEC can also be analytically confirmed by a sum-rule approach [60]. This approach, based on the linear response theory, enables us to derive the frequency of the collective mode without explicitly solving the equation of motion [25]. Applying the approach to the self-gravitating BEC, the angular frequency of the quadrupole mode Ω_Q is generally expressed as

$$\begin{eqnarray} \Omega _Q^2 = \frac{4K-\frac{4}{5}W}{M\langle r^2\rangle }, \end{eqnarray}$$

(39)

where |$\langle r^2\rangle =\int d\boldsymbol {r}r^2\rho (\boldsymbol {r})/M$|⁠. Particularly, within the TF region, the period T_Q = 2π/Ω_Q can be approximated as

$$\begin{eqnarray} T_Q \simeq 2\pi \sqrt{ \frac{5}{4}\frac{M\langle r^2\rangle }{-W} } = 2\pi \sqrt{ \frac{5\pi (\pi ^2-6)}{3} \sqrt{\frac{\hbar ^{6}a^{3}}{G^{5}m^{9}}} } \frac{1}{\sqrt{M}}. \end{eqnarray}$$

(40)

Here, we use Eqs. (8) and (14) to obtain Eq. (40). In Fig. 9, the period of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| is compared with Eq. (40), yielding quantitative consistency. This approach also yields the general angular frequency of the breathing mode Ω_B given by

$$\begin{eqnarray} \Omega _B^2 = \frac{4K+W+9I}{M\langle r^2\rangle }. \end{eqnarray}$$

(41)

Using Eqs. (39), (41), and the virial theorem, 2K + W + 3I = 0, the ratio Ω_B/Ω_Q in the TF regime becomes

$$\begin{eqnarray} \frac{\Omega _B}{\Omega _Q} = \sqrt{ \frac{5}{2} \frac{3+(K/I)}{3+7(K/I)} } \simeq \sqrt{\frac{5}{2}}, \end{eqnarray}$$

(42)

which agrees with our result. Hence, our numerical result regarding the periods of the quadrupole modes is quantitatively consistent with the previous study [60] using a sum-rule approach.

Fig. 9.

The comparison of the period of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| to that from a sum-rule approach. The horizontal axis shows the total mass, whereas the vertical axis shows the negative squared period. The black points show the numerical results for each total mass, whereas the green solid line shows the analytical result written by Eq. (40).

To reproduce our numerical results in more detail, we extend the variational method from the spherical self-gravitating BEC in Sect. 3.1 to the axisymmetric system in the subsequent subsection.

4.2. Analysis of axisymmetric collective mode by variational method

Rindler-Daller and Shapiro extended the TF solution to an ellipsoidal self-gravitating BEC [17,38] based on the ellipsoidal approximation [61]. Given that the self-gravitating BEC takes a spheroidal configuration with semiaxes R_⊥ and Z, the extended TF solution is expressed as

$$\begin{eqnarray} \rho (q) =\frac{\pi M}{4R_\perp ^2Z}j_0\left( \frac{\pi q}{R_\perp }\right), \end{eqnarray}$$

(43)

where M denotes the total mass and q ∈ (0, R_⊥] satisfies (q/R_⊥)² = (r_⊥/R_⊥)² + (z/Z)². Employing Eq. (43), the trial function is set as

$$\begin{eqnarray} \psi (\boldsymbol {r},t) = \sqrt{\frac{\pi M}{4mR_\perp (t)^2Z(t)}j_0\left(\frac{\pi q}{R_\perp (t)}\right)} \exp \left[i\frac{m}{2\hbar }\left(r_\perp ^2H_\perp (t)+z^2H_z(t)\right)\right], \end{eqnarray}$$

(44)

where H_⊥(t) and H_z(t) are variables that provide the velocity field such that

$$\begin{eqnarray} \boldsymbol {v}=\left(r_\perp H_\perp (t)\right)\boldsymbol {e}_\perp +\left(zH_z(t)\right)\boldsymbol {e}_z. \end{eqnarray}$$

(45)

The Lagrangian can be derived by substituting Eq. (44) into its definition (24). Using Eq. (44), the energy components K, W, I defined by Eqs. (5), (6), and (7) are expressed as

$$\begin{eqnarray} & \begin{aligned} K \simeq &\frac{\pi }{24}\frac{\hbar ^2M}{m^2}\mathrm{ F}(\Lambda )\left(\frac{2}{R_\perp (t)^2}+\frac{1}{Z(t)^2}\right) \\ &+\frac{M}{6}\frac{\pi ^2-6}{\pi ^2}\left\lbrace 2H_\perp (t)^2R_\perp (t)^2+H_z(t)^2Z(t)^2\right\rbrace , \end{aligned} \end{eqnarray}$$

(46)

$$\begin{eqnarray} &W = -\displaystyle\frac{3GM^2}{4}\displaystyle\frac{1}{\lbrace R_\perp (t)^2Z(t)\rbrace ^{1/3}}, \end{eqnarray}$$

(47)

$$\begin{eqnarray} &I = \displaystyle\frac{\pi ^2}{4}\displaystyle\frac{\hbar ^2aM^2}{m^3}\displaystyle\frac{1}{R_\perp (t)^2Z(t)}, \end{eqnarray}$$

(48)

where Eq. (47) can be derived using the formula for the gravitational energy of the spheroidal density profile in Ref. [4]. We assume that the deformation of the BEC is sufficiently small.³ Hence, the Lagrangian can be written as

$$\begin{eqnarray} L(R_\perp ,Z,H_\perp ,H_z) =&-\frac{\pi ^2-6}{6\pi ^2}M \left\lbrace 2R_\perp ^2(\dot{H}_\perp +H_\perp ^2) +Z^2(\dot{H}_z+H_z^2) \right\rbrace \\ &-U(R_\perp ,Z). \end{eqnarray}$$

(49)

Here, the effective potential U(R_⊥, Z) is defined as

$$\begin{eqnarray} U(R_\perp ,Z) = \frac{C_z}{3}\left(\frac{2}{R_\perp ^2}+\frac{1}{Z^2}\right)-\frac{C_p}{(R_\perp ^2Z)^{1/3}}+\frac{C_i}{R_\perp ^2Z}, \end{eqnarray}$$

(50)

where C_z, C_p, and C_i are identical to those in Eq. (26). The Euler–Lagrange equations for H_⊥(t) and H_z(t) are |$H_\perp (t)=\dot{R}_\perp (t)/R_\perp (t)$| and |$H_z(t)=\dot{Z}(t)/Z(t)$|⁠. Therefore, the Euler–Lagrange equations for R_⊥(t) and Z(t) are

$$\begin{eqnarray} \left\lbrace \begin{array}{l}\displaystyle\frac{2}{3}\displaystyle\frac{\pi ^2-6}{\pi ^2}M\displaystyle\frac{d^2R_\perp (t)}{dt^2}=-\displaystyle\frac{\partial U(R_\perp ,Z)}{\partial R_\perp },\\ \displaystyle\frac{1}{3}\displaystyle\frac{\pi ^2-6}{\pi ^2}M\displaystyle\frac{d^2Z(t)}{dt^2}=-\displaystyle\frac{\partial U(R_\perp ,Z)}{\partial Z}. \end{array}\right.\\ \end{eqnarray}$$

(51),(52)

These Euler–Lagrange equations show the spherical equilibrium state with R_⊥ = Z = R_eq, which is consistent with the results in Sect. 3.1.

Similarly to the spherically symmetric case, we consider a small oscillation of the semiaxes near R_eq. The fluctuations in the semiaxes are given by R_⊥(t) = R_eq + δR_⊥(t) and Z(t) = R_eq + δZ(t) (|δR_⊥(t)|,  |δZ(t)| ≪ R_eq). The equations of motion for the fluctuations δR_⊥(t) and δZ(t) are given by

$$\begin{eqnarray} \frac{3(\pi ^2-6)}{\pi ^2}MR_{\rm eq}^5\frac{d^2}{dt^2} {\begin{pmatrix}2\delta R_\perp (t) \\ \delta Z(t) \end{pmatrix}} = - {\begin{pmatrix}4C_pR_{\rm eq}^2 & 2(-C_pR_{\rm eq}^2+9C_i) \\ -C_pR_{\rm eq}^2+9C_i & 5C_pR_{\rm eq}^2-9C_i \end{pmatrix}} {\begin{pmatrix}2\delta R_\perp (t) \\ \delta Z(t) \end{pmatrix}} . \end{eqnarray}$$

(53)

When we write the solutions as 2δR_⊥(t) = Aexp [iωt] and δZ(t) = Bexp [iωt], the angular frequency ω satisfies

$$\begin{eqnarray} \omega = \pm \sqrt{ \frac{2\pi ^2}{\pi ^2-6}\frac{C_zR_{\rm eq}+3C_i}{MR_{\rm eq}^5} } \equiv \pm \omega _B \end{eqnarray}$$

(54)

or

$$\begin{eqnarray} \omega = \pm 2 \sqrt{ \frac{\pi ^2}{\pi ^2-6}\frac{C_z}{MR_{\rm eq}^4} } \equiv \pm \omega _Q . \end{eqnarray}$$

(55)

Equation (54) provides the eigenfrequency of the breathing mode compared to the coefficient of Eq. (29). Indeed, when ω = ±ω_B, the eigenmodes satisfy A = 2B, i.e.

$$\begin{eqnarray} {\begin{pmatrix}\delta R_\perp (t) \\ \delta Z(t) \end{pmatrix}} = A {\begin{pmatrix}1 \\ 1 \end{pmatrix}} \exp \left[\pm i\omega _B t\right]. \end{eqnarray}$$

(56)

This indicates that the BEC undergoes either spherical expansion or shrinkage. Conversely, when ω = ±ω_Q, the eigenmodes satisfy A + B = 0, i.e.

$$\begin{eqnarray} {\begin{pmatrix}\delta R_\perp (t) \\ \delta Z(t) \end{pmatrix}} = A {\begin{pmatrix}1 \\ -2 \end{pmatrix}} \exp \left[\pm i\omega _Q t\right]. \end{eqnarray}$$

(57)

We can consider Eq. (57) as the quadrupole mode because it shows that the BEC elongates along the z-axis or expands perpendicular to the z-axis.

We compare Eqs. (56) and (57) with the numerical results. In our numerical simulations, we establish two initial conditions: δR_⊥(t = 0) = δZ(t = 0) = 0 and |$\boldsymbol {\tilde{v}}(\tilde{t}=0,\tilde{r}_\perp ,\tilde{z})=\tilde{\alpha }\tilde{r}_\perp \boldsymbol {e}_\perp +\tilde{\beta }\tilde{z}\boldsymbol {e}_z$|⁠. As a result, the fluctuations δR_⊥(t) and δZ(t) become Consequently, the perturbations δR_⊥(t) and δZ(t) are given by

$$\begin{eqnarray} \delta R_\perp (t) = \frac{\tilde{\lambda }^2c^2mR_{\text {eff}}}{3\hbar } \left\lbrace \frac{2\tilde{\alpha }+\tilde{\beta }}{\omega _B}\sin \left(\omega _B t\right) + \frac{\tilde{\alpha }-\tilde{\beta }}{\omega _Q}\sin \left(\omega _Q t\right) \right\rbrace \end{eqnarray}$$

(58)

and

$$\begin{eqnarray} \delta Z(t) = \frac{\tilde{\lambda }^2c^2mR_{\text {eff}}}{3\hbar } \left\lbrace \frac{2\tilde{\alpha }+\tilde{\beta }}{\omega _B}\sin \left(\omega _B t\right) - 2\frac{\tilde{\alpha }-\tilde{\beta }}{\omega _Q}\sin \left(\omega _Q t\right) \right\rbrace . \end{eqnarray}$$

(59)

The effective widths |$R_{\text {eff}}^\perp$| and Z_eff are, respectively, given by |$R_{\text {eff}}^\perp =\sqrt{2/3}\sqrt{(\pi ^2-6)/\pi ^2}R_\perp$| and |$Z_{\text {eff}}=\sqrt{(\pi ^2-6)/(3\pi ^2)}Z$| using Eq. (44). Thus, we obtain

$$\begin{eqnarray} \begin{split} \frac{Z_{\text {eff}}}{R_{\text {eff}}^\perp } = \frac{Z}{\sqrt{2}R_\perp } \simeq \frac{1}{\sqrt{2}}-\frac{\tilde{\lambda }^2c^2m(\tilde{\alpha }-\tilde{\beta })}{\sqrt{2}\hbar \omega _Q}\sin \left(\omega _Q t\right), \end{split} \end{eqnarray}$$

(60)

which shows that the quadrupole mode induces a harmonic oscillation of |$Z_{\text {eff}}/R_{\text {eff}}^\perp$| near |$1/\sqrt{2}$|⁠. This confirms that the quadrupole mode appears in the numerical simulations (Fig. 6 or Fig. 8(a) in Sect. 4.1).

Our variational calculation of the axisymmetric self-gravitating BEC suggests that solely the quadrupole mode can be extracted through the selection of an appropriate initial phase. As the first term on the right-hand side of Eqs. (58) and (59) vanishes, the exclusive extraction of the quadrupole mode is feasible by setting the initial phase as |$2\tilde{\alpha }+\tilde{\beta }=0$|⁠. Consequently, the quadrupole mode shows

$$\begin{eqnarray} \delta R_\perp (t) = \frac{\tilde{\alpha }\tilde{\lambda }^2c^2mR_{\rm eq}}{\hbar \omega _Q} \sin \left(\omega _Q t\right) \end{eqnarray}$$

(61)

and

$$\begin{eqnarray} \delta Z(t) = -2 \frac{\tilde{\alpha }\tilde{\lambda }^2c^2mR_{\rm eq}}{\hbar \omega _Q} \sin \left(\omega _Q t\right), \end{eqnarray}$$

(62)

resulting in harmonic oscillations of |$R_{\text {eff}}^\perp$| and Z_eff in opposite phase. Indeed, we can numerically extract only the quadrupole mode when the total mass is M = 4 × 10¹⁴M_⊙ and the initial phase is set as |$\tilde{\alpha }=0.03$| and |$\tilde{\beta }=-0.06$|⁠, as shown in Fig. 10. This figure shows that |$\Delta R_{\text {eff}}^\perp \equiv R_{\text {eff}}^\perp (t)-R_{\text {eff}}^\perp (0)$| and ΔZ_eff ≡ Z_eff(t) − Z_eff(0) oscillate monotonically in opposite phases, with periods identical to those shown in Fig. 6. Additionally, the amplitude of ΔZ_eff is approximately 1.40 times larger than that of |$\Delta R_{\text {eff}}^\perp$|⁠, aligning with our analysis, as Eqs. (61) and (62) predict that the ratio of the amplitudes of ΔZ_eff and |$\Delta R_{\text {eff}}^\perp$| is |$\sqrt{2}$|⁠. Therefore, we successfully extract only the quadrupole mode of a self-gravitating BEC.

Fig. 10.

The variation in the effective widths |$\Delta R_{\text {eff}}^\perp =R_{\text {eff}}^\perp (t)-R_{\text {eff}}^\perp (0)$| and ΔZ_eff = Z_eff(t) − Z_eff(0) when a total mass is M = 4 × 10¹⁴M_⊙ and the initial phase is set by |$\tilde{\alpha }=0.03$| and |$\tilde{\beta }=-0.06$|⁠. The horizontal axis shows time and the vertical axis shows |$\Delta R_{\text {eff}}^\perp$| and ΔZ_eff. The blue solid line shows the time evolution of |$\Delta R_{\text {eff}}^\perp$|⁠, and the dashed red line shows the time evolution of ΔZ_eff. The dotted black line shows that |$\Delta R_{\text {eff}}^\perp =\Delta Z_{\text {eff}}=0$|⁠.

Although this variational method qualitatively agrees with the numerical results, the frequency or period of the quadrupole mode diverges from them. The ratio of ω_Q to ω_B is given by

$$\begin{eqnarray} \frac{\omega _Q}{\omega _B} =2\sqrt{\frac{C_zR_{\rm eq}}{C_pR_{\rm eq}^2+3C_i}}, \end{eqnarray}$$

(63)

using Eqs. (28), (54), and (55). However, this ratio almost vanishes because we use |$C_pR_{\rm eq}^2\simeq 3C_i$| and C_zR_eq ≪ C_i in the TF approximation. This is different from the numerical results indicating |$\omega _Q/\omega _B=\sqrt{2/5}$| (see Fig. 8(b)). This inconsistency arises because the trial function in Eq. (44) is based on the TF solution, which neglects kinetic energy and becomes invalid near the surface of the BEC. However, Eqs. (13) and (55) show that the kinetic energy of the equilibrium state predominantly influences the frequency of the quadrupole mode. Notably, the quadrupole mode does not oscillate the central density, and neither does it involve potential nor self-interaction energy contributions. Consequently, it only changes the density near the surface where the TF approximation fails. Hence, a trial function based on the TF approximation prevents the quantitative evaluation of the kinetic energy, which affects the frequency of the quadrupole mode. To accurately compare ω_Q with the numerical results, we need to introduce a first-order correction to the kinetic energy of the equilibrium state of the self-gravitating BEC. This investigation will be conducted in a future study.

5. Conclusions

We study the collective modes of self-gravitating BECs. In particular, we focus on the breathing and anisotropic collective modes. In this study, we show the following three aspects.

Firstly, a self-gravitating BEC can induce a breathing mode by the introduction of an isotropic initial velocity field to its equilibrium state. Due to the time-dependent density profile of the BEC, the gravitational potential near the center of the BEC oscillates in this mode. The oscillation is characterized by a harmonic oscillation of the radius of the BEC, similarly to that of a conventional BEC trapped by an external potential. In the self-gravitating BEC, as the total mass M increases, the amplitude decreases and the period becomes short in proportion to |$1/\sqrt{M}$|⁠. These distinctive properties distinctly reflect the density dependence of self-gravity.

Secondly, the axisymmetric initial velocity field yields an anisotropic collective mode. It is a mixture of the quadrupole mode and the breathing mode. The quadrupole mode can be extracted by focusing on the oscillation of the aspect ratio of the BEC. The amplitude and period in the quadrupole mode decrease as the total mass increases, similarly to the radius of the BEC in the breathing mode. Particularly, the period of the quadrupole mode has the same M-dependence, in proportion to |$1/\sqrt{M}$|⁠, as that of the breathing mode. This property of the self-gravitating BEC differs from that of a conventional BEC trapped by an external potential. Although the M-dependence of the period is different for the self-gravitating BEC and that in an external isotropic harmonic potential, the ratio of the period of the quadrupole mode to that of the breathing mode takes the same value of |$\sqrt{5/2}$| in both cases. These properties of the period of the quadrupole mode are also shown analytically in a previous study employing a sum-rule approach [60].

Thirdly, we can also extract only the quadrupole mode from the anisotropic collective mode by setting an appropriate initial velocity. We extend the variational method, in which the trial function is based on the TF solution, to a spheroidal configuration. With this approach, we succeed in describing the anisotropic collective mode of the self-gravitating BEC and proposing appropriate initial conditions for the pure quadrupole mode. However, this extended variational method does not work for the evaluation of the kinetic energy of the equilibrium state. Consequently, the frequency of the quadrupole mode cannot be accurately estimated because it depends critically on the kinetic energy of the equilibrium state. To quantitatively reproduce the frequency, it is necessary to introduce a first-order correction to the kinetic energy.

In instances where the oscillation exhibits a large amplitude, the density dependency inherent in the self-gravitating BEC probably changes the M-dependence of the periods and other characteristics of the collective mode. In this study, we consider the amplitude to be small. However, the large deformation of a self-gravitating BEC affects the gravitational energy. Accounting for this deformation, a collective mode with a large amplitude can cause various nontrivial phenomena triggered by the difference in whether the spheroidal configuration is oblate or prolate.

If the self-gravitating BEC explains a DM halo, we expect that the collective oscillations of the self-gravitating BEC provide evidence that the DM halo is composed of BEC. This is a similar situation to early studies on atomic BECs, whose collective mode played an important role in determining whether a system was a BEC. For instance, the collective modes investigated in the present paper may be caused by the collision and merging of BECs. The timescales of these oscillations are shorter than the dynamical timescales of the cluster of galaxies, the breathing modes excite acoustic waves in baryonic matter, and the anisotropic modes, including quadrupole moments, excite gravitational waves. The frequencies of these oscillations are determined mainly by the self-interaction of the Bose particles, and these frequencies depend on the total mass of the BEC, the boson mass, and the s-wave scattering length. Hence, these waves provide important information about the Bose particle and the BEC in the universe. Observational study of such waves with long wavelengths is a challenging work.

One of our next interests is the collective modes of a rotating self-gravitating BEC. In the universe, objects are naturally rotating rather than in static states. The equilibrium state of a self-gravitating BEC is affected by rotation owing to its density dependence. For example, the presence of a quantized vortex within a BEC locally pushes the density through rotation, changing the density profile of the equilibrium state. Since the gravitational potential depends on the density distribution, it differs from that of a BEC without a quantized vortex. Hence, the self-gravitating BEC with rotation or a quantized vortex likely provides distinct phenomena of the collective modes from those in the present work. Clarifying this property is an interesting next task.

Acknowledgement

We thank M. Kobayashi and Y. Sano for their fruitful discussions of our numerical calculations.

Funding

This work was supported by Japan Science and Technology Agency (JST) and the establishment of university fellowships towards the creation of science technology innovation, Grant Number JPMJFS2138. M.T. acknowledges the support from Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Number JP23K03305 and JP22H05139. This work was partially supported by Osaka Central Advanced Mathematical Institute: MEXT Joint Usage/Research Center on Mathematics and Theoretical Physics JPMXP0619217849.

Footnotes

1

We introduce |$\tilde{\lambda }$| to keep the values of the dimensionless variables within a suitable range for the numerical calculations.

2

These Euler–Lagrange equations imply that the motion of the spherical self-gravitating BEC resembles that of the Friedmann model because Eq. (27) corresponds to the Friedmann equation. Here, R and H are the scale factor and the Hubble parameter, respectively.

3

In general, the gravitational energy of a spheroid is affected by its deformation (see Ref. [4]). However, we neglect these effects because, in the present study, we address only small oscillations.