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

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.


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 non-emissive, non-reflective, and non-absorptive 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 non-relativistic 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, for example, ultralight scalar particles behave as DM [8,9] or non-gravitational 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][13][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 −23 − 10 −22 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 9 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, for example, 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 4 He, 3 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 non-negligible non-gravitational 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 selfgravitating 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 semi-axes 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 proves valuable in describing the equilibrium state of a BEC and carrying out the analytical studies.As for a self-gravitating BEC, this approximation was applied for the first time to compare the 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 selfgravitating 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 three-dimensional system without spatial symmetry, focusing on quantized vortices in self-gravitating BECs, such as their stability [47,48] and collisions between two self-gravitating BECs with quantized vortices [49].
In this work, we study the three-dimensional dynamics of the collective modes of a selfgravitating 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.First, 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 Sec.II, Based on this, the breathing mode of a self-gravitating BEC is studied in Sec.III.Anisotropic collective mode under self-gravity of the BEC is investigated in Sec.IV.Finally, we conclude this paper in Sec.V.

GROSS-PITAEVSKII-POISSON MODEL AND THE EQUILIBRIUM STATE
We consider a self-gravitating BEC composed of scalar bosons with mass m and an swave scattering length a at zero temperature.In the GPP model, the physical state of the BEC is described by the macroscopic wavefunction ψ(r, t), and the time evolution is given by the following GPP equations defined as where V (r, t) denotes the gravitational potential [14,17,29].Using the Madelung representation written by the GP equation yields the continuity equation and the Euler-like equation, where ρ = m|ψ| 2 denotes the density and v = ℏ∇θ/m denotes the velocity field [15].Then, the total mass is 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] We can obtain the equilibrium solution of the GPP equations (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 ψ(r, t) = ϕ(r)e −iµt/ℏ , where µ denotes the chemical potential, the steady solution of Eq. ( 1) is approximately given by where r = |r| represents the radial coordinate, ρ c stands for the central density, and j 0 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: 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 where M TF is defined as 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 Therefore, using Eqs.( 5), ( 6), ( 7), (8), and ( 12), the energy of each equilibrium state can be expressed as 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 ξ = m/(8πaρ c ) and the TF radius, such that Λ ∼ 1 − ξ/R TF .

NUMERICAL SETUP
The GPP equations (1) are expressed in the following dimensionless formulation: where the tilde symbols represent the dimensionless variables, and as explained later, the term −i Ṽs in the dimensionless GP equation is artificially added due to a computational reason.The dimensionless variables are defined as follows: r = r( λmc/ℏ), t = [44].The scaling factor λ can adjust the system size using the invariance of the equations for any value1 .
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 Ṽ ′ (r, t; s): where s denotes a real parameter.The gravitational potential is obtained from Ṽ (r, t) = lim s→∞ Ṽ ′ (r, t; s).Furthermore, the efficacy of the Jacobi method can be notably enhanced through Nesterov's accelerated gradient method [52][53][54].This method is common in the field of machine learning and enables us to rapidly converge the function Ṽ ′ (r, t; s) toward the gravitational potential Ṽ (r, t).The function Ṽ ′ (r, t; s) satisfies the Dirichlet boundary condition given by where Rb denotes the position on the boundary of the numerical domain and Rc denotes the center of the numerical domain.We apply a time-dependent total mass M ( 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 Ṽs , called "sponge" potential, written by is introduced into Eq.( 16) [44].Here, Ṽo denotes the amplitude.This sponge potential can reduce the number of particles within the range of r > rc , where δ 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 φ(r)] to the equilibrium state ψeq (r): The BEC is located at the center of the numerical domain.We obtain the equilibrium state ψeq using the imaginary time propagation method of the GPP equations.The initial phase is given by φ which gives rise to the initial velocity field ṽ = ∇ φ = (αr ⊥ )e ⊥ + ( β z)e z .Here, r⊥ = x2 + ỹ2 and z denote the 3D cylindrical coordinates, and the unit vectors along each direction are e ⊥ and e z .We determine the initial velocity by α and β (e.g., it becomes isotropic when α = β).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 one percent, and particles are hardly emitted from the BEC.

BREATHING MODE
In this section, we investigate the breathing modes of 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.

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 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 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 using Eq. ( 4) [39].Upon substituting the trial function (23), the Lagrangian transforms into Here, the effective potential U (R) is defined as Here the coefficients are , and the radius R(t) follows the equation2 given by The equilibrium solution R eq is acquired by setting dU (R)/dR = 0, yielding Fig. 1 describes the result of Eq. ( 28) with Λ = 0.97 as the cut-off 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.
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 Fig. 1 The relationship between the total mass M and the radius R derived from Eq. ( 9) with the cut-off 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 .
Eqs. ( 26), (27), and (28), we obtain the following equation of motion: 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 which implies a harmonic oscillation with a period 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.

NUMERICAL RESULTS
To confirm Eq. (31), we conduct numerical simulations for the total mass M/(10 14 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 12 M ⊙ [37], and the Milky Way is estimated at M ≃ 0.5 − 2.0 × 10 12 M ⊙ [56].However, the substantial total mass facilitates the exploration of BEC oscillations because R eq hardly changes even if the total < l a t e x i t s h a 1 _ b a s e 6 4 = " v R E z g 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " v R E z g 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " v R E z g 1 mass is reduced (see Fig 1).For a BEC with the total mass of a typical galaxy, namely M ∼ 10 10 − 10 12 M ⊙ , the time scale 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 α = β = ±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 three-dimensional 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 and where Ω denotes a solid angle and ∆r ≈ 0.22 kpc is set so that ∆r matches the spatial resolution, namely ∆r = Lx / Ñx =0.3125.Fig. 2 shows the time evolution of the density and To quantitatively observe the spherical oscillation, we define the effective radius R eff as 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.Fig. 3(a) shows the deviation of the effective radius from its initial value, ∆R eff (t) = R eff (t) − R eff (0), for α = 0.03.The graph shows sinusoidal curves for each total mass, where the oscillation period T R eff decreases as the total mass M is increased.Fig. 3(b) suggests that our results of T R eff agree quantitatively with Eq. ( 31), even for three values of α.Therefore, we conclude that the self-gravitating BEC induces the breathing mode.
To describe the time evolution of the BEC, we obtain the density and gravitational potential profiles by averaging ρ(r, t) and V (r, t): and where ∆r ⊥ ≈ 0.22 kpc is determined in the same manner as ∆r, namely ∆r ⊥ = Lx / Ñx = 0.3125, and θ = tan −1 (y/x).Fig. 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, while 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 (α, β) such as α ̸ = β (e.g.(α, β) = (±0.1,∓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 14 M ⊙ and is endowed with an initial velocity by α = 0.1, β = 0.The horizontal axis shows r ⊥ and z.Magenta points show the profiles at z = 0, ρ(r ⊥ , z = 0) and V (r ⊥ , z = 0), while 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 ⊥ eff and Z eff defined as These are the effective widths of the BEC in each direction; R ⊥ eff characterizes the width in the direction of the radial coordinate r ⊥ and Z eff is the width along the z-axis.Fig. 5 shows the complicated oscillatory behavior of these quantities, describing the opposing oscillations between R ⊥ eff 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.The horizontal axis shows the time.The left column shows the result of (α, β) = (±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 (α > 0), and the black dashed line shows the oscillation when the BEC initially shrinks (α < 0).Conversely, the right column shows the result of (α, β) = (0, ±0.03), indicating the initial velocity along the z-axis.The red solid line shows the oscillation when the BEC initially expands ( β > 0), and the blue dashed line shows the BEC initially shrinks ( β < 0).While the upper row shows the time evolution of R ⊥ eff , the lower row shows the time evolution of Z eff .
Although R ⊥ eff and Z eff manifest intricate oscillations, a harmonic oscillation can be obtained from the ratio Z eff /R ⊥ eff , as shown in Fig. 6.Among the four cases, the oscillations have a common amplitude and frequency despite the different behaviors of R ⊥ eff , Z eff .The ratio Z eff /R ⊥ eff indicates the extent to which the BEC axisymmetrically deforms from Fig. 6 The oscillation of the ratio Z eff /R ⊥ eff when the total mass of the self-gravitating BEC is M = 4 × 10 14 M ⊙ .The horizontal axis shows the time, and the vertical axis shows the value of the ratio Z eff /R ⊥ eff .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).
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][58][59].
The Fourier transformation serves as a valuable tool for decomposing complicated oscillations into eigenmodes.Fig. 7 shows the Fourier transformation of R ⊥ eff and Z eff /R ⊥ eff , as shown in Fig. 5 and Fig. 6.Within this representation, two peaks of R ⊥ eff and a peak of Z eff /R ⊥ eff can be discerned.The high-frequency peak of R ⊥ eff corresponds to the period of the breathing mode T B with the corresponding total mass.On the other hand, the lowfrequency peak in R ⊥ eff closely matches that of Z eff /R ⊥ eff .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.
We vary the parameters of the initial phase and total mass to investigate the dependence of the period of Z eff /R ⊥ eff on these quantities.Although we study the dynamics of various combinations of α and β, the period remains constant for all the simulations.Thus, the period of Z eff /R ⊥ eff is independent of the initial velocity field.However, the total mass can change the period of Z eff /R ⊥ eff , as shown in Fig. 8. Fig. 8 (a) shows that Z eff /R ⊥ eff exhibits harmonic oscillation similar to R eff of the breathing mode induced by the isotropic initial Fig. 7 The Fourier transformation of R ⊥ eff and Z eff /R ⊥ eff .The self-gravitating BEC has a total mass M = 4 × 10 14 M ⊙ , with an the initial velocity given by α = 0.03, β = 0.The horizontal axis shows the frequency, while the vertical axis shows the amplitude.The solid lines show the results of each transformation (red: R ⊥ eff and blue: Z eff /R ⊥ eff ), and the green dashed line shows the analytical/variational result of the breathing mode T B .The red points show the peaks of R ⊥ eff , and the blue triangle shows the peak of Z eff /R ⊥ eff .
velocity field; the period becomes short as the total mass increases.The M -dependence of the period of Z eff /R ⊥ eff and its comparison with T B are shown in Fig. 8 (b).We find that Z eff /R ⊥ eff 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 √ M , similar to the breathing mode.However, the period of Z eff /R ⊥ eff 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 be also observed in atomic BECs confined by an isotropic harmonic potential.The frequency of the breathing mode is 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 eff /R ⊥ eff to T B closely resembles that of the BEC confined by an isotropic harmonic potential.Moreover, the periods of Z eff /R ⊥ eff are in good agreement with 5/2T B (see Fig. 8).Thus, the numerical results indicate that the quadrupole mode has the period 5/2T 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 The oscillation of Z eff /R ⊥ eff for each total mass.The initial phase is determined by α = 0.03, β = 0. (a) The time evolution of Z eff /R ⊥ eff for each total mass.The horizontal axis shows time, while the vertical axis shows the value of Z eff /R ⊥ eff .The black dotted line shows the ratio when the BEC is in equilibrium; the ratio becomes solid line shows the result (red: 4 × 10 14 M ⊙ , blue: 3 × 10 14 M ⊙ , green: 2 × 10 14 M ⊙ , black: The dependence of the period of Z eff /R ⊥ eff on the total mass by the log-log plot.The horizontal axis shows the total mass, while 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 5/2T B .frequency of the quadrupole mode Ω Q is generally expressed as where ⟨r 2 ⟩ = drr 2 ρ(r)/M .Particularly, within the TF region, the period T Q = 2π/Ω Q can be approximated as Here, we use Eqs.( 8) and ( 14) to obtain Eq. (40).In Fig. 9, the period of Z eff /R ⊥ eff is compared with Eq. ( 40), yielding quantitative consistency.This approach also yields the general angular frequency of the breathing mode Ω B given by Fig. 9 The comparison of the period of Z eff /R ⊥ eff to that from a sum-rule approach.The horizontal axis shows the total mass, while 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).
Using Eqs. ( 39), (41), and the virial theorem, 2K + W + 3I = 0, the ratio Ω B /Ω Q in the TF regime becomes 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 approach.
To reproduce our numerical results in more detail, we extend the variational method from the spherical self-gravitating BEC in Sec.III-A to the axisymmetric system in the subsequent subsection.

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 semi-axes R ⊥ and Z, the extended TF solution is expressed as where M denotes the total mass and q ∈ (0, R ⊥ ] satisfies (q/R ⊥ ) 2 = (r ⊥ /R ⊥ ) 2 + (z/Z) 2 .
Employing Eq. ( 43), the trial function is set as where H ⊥ (t) and H z (t) are variables that provide the velocity field such that 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 where Eq. ( 47) can be derived using the formula for the gravitational energy of the spheroidal density profile in [4].We assume that the deformation of the BEC is sufficiently small3 .Hence, the Lagrangian can be written as Here, the effective potential U (R ⊥ , Z) is defined as 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 ⊥ (t) = Ṙ⊥ (t)/R ⊥ (t) and H z (t) = Ż(t)/Z(t).Therefore, the Euler-Lagrange equations for R ⊥ (t) and These Euler-Lagrange equations show the spherical equilibrium state with R ⊥ = Z = R eq , which is consistent with the results in Sec.III-A.
This indicates that the BEC either undergoes spherical expansion or shrinkage.Conversely, when ω = ±ω Q , the eigenmodes satisfy A + B = 0, i.e.
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.
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) vanish, the exclusive extraction of the quadrupole mode is feasible by setting the initial phase as 2α + β = 0. Consequently, the quadrupole mode shows and resulting in harmonic oscillations of R ⊥ eff and Z eff in opposite phase.Indeed, we can numerically extract only the quadrupole mode when the total mass is M = 4 × 10 14 M ⊙ and the initial phase is set as α = 0.03 and β = −0.06,as shown in Fig. 10.This figure shows that 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 ∆R ⊥ eff , aligning with our analysis, as Eqs.( 61) and (62) predict that the ratio of the amplitudes of ∆Z eff and ∆R ⊥ eff is √ 2. Therefore, we successfully extract only the quadrupole mode of self-gravitating BEC.
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 ω using Eqs.( 28), ( 54) and ( 55).However, this ratio almost vanishes because we use C p R 2 eq ≃ 3C i and C z R eq ≪ C i in the TF approximation.This is different from the numerical results indicating ω Q /ω B = 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.

CONCLUSIONS
We study the collective modes of the self-gravitating BECs.In particular, we focus on the breathing and anisotropic collective modes.In this study, we show the following three aspects.
First, 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 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, similar to 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/ √ 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, similar to the radius of BEC in the breathing mode.Particularly, the period of the quadrupole mode has the same M -dependence, in proportion to 1/ √ 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 the one 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 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 to 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 time scales of these oscillations are shorter than the dynamical time scales 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 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 quantized vortex likely provides distinct phenomena of the collective modes from those in the present work.Clarifying this property is an interesting next task.

Finally, the physical
and numerical parameters are specified.We consider the case in which each boson has a mass m = 3 × 10 −24 eV and an s-wave scattering length a = 5.62 × 10 −98 λ c m, where λ c denotes the Compton wavelength [47].For the numerical analyses, we use a 3D cube of length Lx = Ly = Lz = 40.The scaling factor is λ = 3 × 10 −3 , ensuring that the numerical domain vastly exceeds the size of the BEC.The spatial grid is configured as Ñx = Ñy = Ñz = 128.The time resolution of the GP equation is ∆ t = 10 −3 , and the resolution of the parameter s is 0.05 × ( Lx / Ñx ) 2 ≈ 4.88 × 10 −3 .The parameters of the sponge potential are designated as Ṽo = 1000, rc = 0.8 × ( Lx /2) = 16 and δ = 2 × ( Lx / Ñx ) = 0.625.Then, rc exceeds the dimensionless TF radius RTF ≈ 9.09, which Time [Myr] < l a t e x i t s h a 1 _ b a s e 6 4 = " v R E z g 1 L T P + U o l 8 G A e 4 f w x x I B g 2 o

Fig. 2
Fig.2The density profile and the gravitational potential of the self-gravitating BEC.The total mass is M = 2 × 10 14 M ⊙ , and the initial phase is set by α = β = 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, while the black dashed line shows the initial data.The upper row displays the density profile, with the density depicted on the vertical axis, while the lower row exhibits the gravitational potential, with its value presented on the vertical axis.

Fig. 3
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 α = 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 14 M ⊙ , blue: 3 × 10 14 M ⊙ , green: 2 × 10 14 M ⊙ and black: 1 × 10 14 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 eff for various values of α.The horizontal axis shows the total mass, and the vertical axis shows the negative squared period 1/T 2 R eff .Each point shows the numerical result (black square: α = 0.1, red solid circle: α = 0.03 and green triangle: α = −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 2 B .
Fig. 8The oscillation of Z eff /R ⊥ eff for each total mass.The initial phase is determined by α = 0.03, β = 0. (a) The time evolution of Z eff /R ⊥ eff for each total mass.The horizontal axis shows time, while the vertical axis shows the value of Z eff /R ⊥ eff .The black dotted line shows the ratio when the BEC is in equilibrium; the ratio becomes Z eff /R ⊥ eff = 1/ 0) and ∆Z eff = Z eff (t) − Z eff (0) when a total mass is M = 4 × 10 14 M ⊙ and the initial phase is set by α = 0.03 and β = −0.06.The horizontal axis shows time and the vertical axis shows ∆R ⊥ eff and ∆Z eff .The blue solid line shows the time evolution of ∆R ⊥ eff , and the dashed red line shows the time evolution of ∆Z eff .The dotted black line shows that ∆R ⊥ eff = ∆Z eff = 0.