1 Introduction

Gravity has been conjectured to be the weakest force among all other fundamental interactions a long time ago [1], and it is also conjectured to be weaker than the fifth-force mediated by other scalar fields, namely the dubbed scalar weak gravity conjecture [2] and its strong version [3], from which a refined trans-Planckian censorship conjecture [4] is proposed as the bottomline constraint on the effective theories of cosmic inflation. Only until recently, the weak gravity conjecture is elaborated with a general argument from the positivity of black hole entropy shift from higher-dimension operators [5]. Due to the weakestness of gravity, it is generally difficult to constrain deviation beyond the Einstein gravity (or Newtonian gravitational potential in the non-relativistic limit and the weak gravitational field regime) at small scales [6], despite of numerous constraints from the astrophysical/cosmological scales [7,8,9,10,11,12,13,14,15].

To manifest the small-scale effects in experiments, the Bose–Einstein condensate (BEC) is extensively used to collectively exhibit macroscopic quantum properties that are typically displayed at microscopic scales due to the same quantum identity for all the atoms in a BEC system. However, the cold atom BEC experiments in the terrestrial environment are usually limited by its cooling temperature (nK) and observation time (10 ms–1 s) due to the large gravity pull on Earth. Therefore, there is a growing interest recently [16,17,18] to implement the ultra-cold atom experiment in space to achieve a colder temperature (pK) and longer observation time (10 s–100 s) thanks to the microgravity environment (gravity gradient around \(0.0001\,\mathrm {g/m}\)). The ongoing projects have been proposed by NASA, China, ESA, Germany and France, most of which are focusing on the cold atom interferometry to investigate the modified gravity [19], to search for the ultra-light dark matter [20], and to detect the gravitational waves [20] in the sensitivity range between LISA and LIGO.

In this paper, we will mainly consider the ultra-cold dilute atomic-gas forming a BEC in a spherical harmonic trapping potential with shape oscillation. The long lifetime of BEC state allows for a precise measurement on the shape oscillation frequency, and constraints on the modified gravity is expected in principle. It is worth noting that our preliminary results only serve for theoretical interest since the experimental implements are still quite challenging, which will not be our focus in the current paper.

The paper is organized as follows: In Sect. 2, we derive the general form of the Gross–Pitaevskii equation in the presence of a gravitational potential, from which the frequency deviations in the shape oscillation with respect to the gravity-free case are derived for some illustrative modified gravity theories in Sect. 3. In Sect. 4, we give a preliminary perspective on the experimental constraints on these modified gravity theories. We summarize our results in Sect. 5.

2 Gross-Pitaevskii equation with gravitational potential

Consider N weakly interacting cold atoms of mass m constituting a dilute Bose gas trapped in an external potential \(V_\mathrm {ext}\) that would form BEC in the laboratory with aid of laser cooling and evaporative cooling techniques, its mean-field dynamics of the macroscopic wave function \(\Psi (\vec {r}_1,\cdots ,\vec {r}_N)\approx \Pi _{i=1}^N\psi (\vec {r}_i)\) could be essentially captured by a non-linear Schrödinger equation dubbed the Gross–Pitaevskii equation (GPE),

$$\begin{aligned} i\hbar \frac{\partial }{\partial t}\psi (t,\vec {r}\,)=&\left[ -\frac{\hbar ^2}{2m}\vec {\nabla }^2+V_\mathrm {ext}(\vec {r}\,)\right. \nonumber \\&\left. +N\int \mathrm {d}^3\vec {r'}\,\,V_\mathrm {int}(\vec {r}-\vec {r'})|\psi (t,\vec {r'})|^2\right] \psi (t,\vec {r}\,), \end{aligned}$$
(1)

where the interacting potential \(V_\mathrm {int}\) consists of the usual s-wave scattering potential as well as a gravitational potential,

$$\begin{aligned} V_\mathrm {int}(\vec {r}-\vec {r'})=g \delta ^3(\vec {r}-\vec {r'})+V_G(|\vec {r}-\vec {r'}|) \end{aligned}$$
(2)

with \(g \equiv 4\pi \hbar ^2a/m\) characterizing the strength of the s-wave scattering of length a. Due to the extremely low temperature achieved in the ultra-cold atom experiments, the atoms in BEC are almost collision-less so that all higher-order partial-wave collisions are suppressed at zero collision energy in a short-ranged potential. The external trapping potential \(V_\mathrm {ext}\) is assumed to realize a perfect spherically symmetric form in space,

$$\begin{aligned} V_\mathrm {ext}(\vec {r}\,)=\frac{1}{2}m\omega _0^2r^2, \end{aligned}$$
(3)

defining a characteristic length scale for the BEC ground state by

$$\begin{aligned} a_0\equiv \sqrt{\frac{\hbar }{m\omega _0}}. \end{aligned}$$
(4)

For the gravitational potential probed by the ultra-cold atom experiment in space, we will focus on the following three kinds of the short-range modifications [21] of the gravitational inverse-square law in the non-relativistic limit and the weak gravitational regime:

  1. 1.

    Power-law potential:

    $$\begin{aligned} V_G^\mathrm {power}(r)=-\frac{Gm^2}{r}\left[ 1+\beta _k\left( \frac{1\,\mathrm {mm}}{r}\right) ^{k-1}\right] , \end{aligned}$$
    (5)

    which could be produced by the simultaneous exchange of multiple massless bosons in the higher-order exchange processes. For example, the \(k=2\) case corresponds to the simultaneous exchange of two massless scalar bosons [22]; the \(k=3\) case corresponds to the simultaneous exchange of massless pseudo-scalar particles between two fermions with the \(\gamma _5\)-couplings [23]; the \(k=5\) case corresponds to the simultaneous exchange of two massless pseudoscalars with the \(\gamma _5\gamma ^\mu \partial ^\mu \) couplings [23] such as the axion or other Goldstone bosons; and the fractional k’s are expected from the unparticle exchange [24].

  2. 2.

    Yukawa interaction:

    $$\begin{aligned} V_G^\mathrm {Yukawa}(r)=-\frac{Gm^2}{r}\left[ 1+\alpha \,\mathrm {e}^{-r/d}\right] , \end{aligned}$$
    (6)

    which could be produced by the exchange of natural-parity bosons between unpolarized bodies with the boson mass \(\hbar c/d\).

  3. 3.

    Fat graviton:

    $$\begin{aligned} F_\mathrm {fat}(r)&=-\frac{Gm^2}{r^2}\left[ 1-\exp \left( -\frac{r^3}{\lambda ^3}\right) \right] , \end{aligned}$$
    (7)
    $$\begin{aligned} V_G^\mathrm {fat}(r)&=-\frac{Gm^2}{r}\left[ 1-E_{4/3}\left( \frac{r^3}{\lambda ^3}\right) \right] \end{aligned}$$
    (8)

    with the exponential integral function \(E_\nu (z)=\int _1^\infty \mathrm {e}^{-zt}t^{-\nu }\mathrm {d}t\). Here the graviton is conjectured to be a “fat” object with size \(\lambda \) [25], and the gravitational force falls off rapidly to zero at \(r\rightarrow 0\) limit [21].

2.1 Equation of motion

The GPE could be derived from the variation of the action [26,27,28]

$$\begin{aligned} S=&\,\int \mathrm {d}tL=\int \mathrm {d}t\int \mathrm {d}^3\vec {r}\,\,\mathcal {L}, \end{aligned}$$
(9)
$$\begin{aligned} \mathcal {L}=&\,\frac{i\hbar }{2}\left( \psi \frac{\partial \psi ^*}{\partial t}-\psi ^*\frac{\partial \psi }{\partial t}\right) +\frac{\hbar ^2}{2m}\vec {\nabla }\psi ^*\cdot \vec {\nabla }\psi +V_\mathrm {ext}|\psi |^2\nonumber \\&\,+\frac{g N}{2}|\psi |^4+\frac{N}{2}|\psi |^2\int \mathrm {d}^3\vec {r'}\,\,V_G(|\vec {r}-\vec {r'}|)|\psi (t,\vec {r'})|^2 \end{aligned}$$
(10)

with respect to \(\delta \psi \) and \(\delta \psi ^*\). For a normalized wave function ansatz,

$$\begin{aligned} \psi (t,\vec {r}\,)=\frac{\mathrm {e}^{i\left[ \gamma (t)+B(t)r^2\right] }}{(\sqrt{\pi }\sigma (t))^{3/2}}\mathrm {e}^{-\frac{r^2}{2\sigma (t)^2}}, \end{aligned}$$
(11)

the width parameter \(\sigma (t)\) and phase parameters \(\gamma (t)\) and B(t) could be solved from the corresponding Euler–Lagrange equations of Lagrangian

$$\begin{aligned} L=&\int \mathrm {d}^3\vec {r}\,\,\mathcal {L}=\hbar \dot{\gamma }+L_G+\frac{g N}{4\sqrt{2\pi ^3}\sigma ^3}\nonumber \\&+\frac{3}{2}\sigma ^2\left( \hbar \dot{B}+\frac{2\hbar ^2}{m}B^2+\frac{\hbar ^2}{2m\sigma ^4}+\frac{1}{2}m\omega _0^2\right) , \end{aligned}$$
(12)

where

$$\begin{aligned} L_G&=\!\frac{N}{2}\int \mathrm {d}^3\vec {r}_1\,|\psi (t,\vec {r}_1)|^2\int \mathrm {d}^2\vec {r}_2\,V_G(|\vec {r}_1\!-\!\vec {r}_2|)|\psi (t,\vec {r}_2)|^2\nonumber \\&=\frac{N}{2\pi ^3}\int \mathrm {d}^3\vec {r}_1\mathrm {d}^3\vec {r}_2V_G(|\vec {r}_1-\vec {r}_2|)\frac{1}{\sigma ^6}\mathrm {e}^{-\frac{r_1^2+r_2^2}{\sigma ^2}} \end{aligned}$$
(13)

will be analytically evaluated later with specific form of the gravitational potential.

The total time-derivative term \(\hbar \dot{\gamma }(t)\) in (12) does not admit a dynamical equation, and the rest of Euler-Lagrange equations for \(\sigma (t)\) and B(t) read

$$\begin{aligned}&B(t)=\frac{m\dot{\sigma }(t)}{2\hbar \sigma (t)}, \end{aligned}$$
(14)
$$\begin{aligned}&\hbar \dot{B}+\frac{2\hbar ^2}{m}B^2-\frac{\hbar ^2}{2m\sigma ^4}+\frac{1}{2}m\omega _0^2=\frac{g N}{4\sqrt{2\pi ^3}\sigma ^5}-\frac{1}{3\sigma }\frac{\mathrm {d}L_G}{\mathrm {d}\sigma }, \end{aligned}$$
(15)

respectively, which could be combined into a single equation of motion (EOM)

$$\begin{aligned} m\ddot{\sigma }=-m\omega _0^2\sigma +\frac{\hbar ^2}{m\sigma ^3}+\frac{g N}{2\sqrt{2\pi ^3}\sigma ^4}-\frac{2}{3}\frac{\mathrm {d}L_G}{\mathrm {d}\sigma }. \end{aligned}$$
(16)

After adopting following dimensionless quantities,

$$\begin{aligned} \tau =\omega _0t, \quad \nu =\frac{\sigma }{a_0}, \end{aligned}$$
(17)

the EOM (16) becomes

$$\begin{aligned} \frac{{\text{ d }}^2\nu }{{\text{ d }}\tau ^2}+\frac{{\text{ d }}U}{{\text{ d }}\nu }=0, \end{aligned}$$
(18)

which effectively describes a particle moving in a dimensionless effective potential

$$\begin{aligned} U(\nu )=\frac{1}{2}\left( \nu ^2+\frac{1}{\nu ^2}\right) +\sqrt{\frac{2}{\pi }}\frac{N}{3\nu ^3} \frac{a}{a_0}+\frac{2}{3}\frac{L_G}{\hbar \omega _0}. \end{aligned}$$
(19)

To evaluate \(L_G\) in (19), we first make replacements of variables as

$$\begin{aligned} \vec {r}=\frac{1}{\sqrt{2}}\left( \vec {r}_1-\vec {r}_2\right) , \quad \vec {R}=\frac{1}{\sqrt{2}}\left( \vec {r}_1+\vec {r}_2\right) , \end{aligned}$$
(20)

so that \(r_1^2+r_2^2=r^2+R^2\), and then \(L_G\) becomes

$$\begin{aligned} L_G\!=\!\frac{N}{2\pi ^3\sigma ^6}\int _0^\infty \mathrm {d}R(4\pi R^2)\mathrm {e}^{-\frac{R^2}{\sigma ^2}}\int _0^\infty \mathrm {d}r(4\pi r^2)V_G(\sqrt{2}r)\mathrm {e}^{-\frac{r^2}{\sigma ^2}}. \end{aligned}$$
(21)

After employing the following dimensionless quantities

$$\begin{aligned} x=\frac{r}{a_0}, \quad a_0^G=\frac{Gm^2}{\hbar \omega _0}\equiv \frac{c}{\omega _0}\frac{m^2}{m_\mathrm {Pl}^2}, \end{aligned}$$
(22)

with the Planck mass \(m_\mathrm {Pl}\equiv \sqrt{\hbar c/G}\), one arrives at

$$\begin{aligned} \frac{L_G}{\hbar \omega _0}=\frac{N}{2\pi ^{3/2}}\frac{a_0^G}{a_0}\frac{I(\nu )}{\nu ^3}, \end{aligned}$$
(23)

where the integration \(I(\nu )\) reads

$$\begin{aligned} I(\nu )=\int _0^\infty \mathrm {d}x(4\pi x^2)\left[ \frac{a_0}{Gm^2}V_G(\sqrt{2}a_0x)\right] \mathrm {e}^{-\frac{x^2}{\nu ^2}}. \end{aligned}$$
(24)

Therefore, the EOM could be solved for given effective potential of form

$$\begin{aligned} U(\nu )=\frac{1}{2}\left( \nu ^2+\frac{1}{\nu ^2}\right) +\sqrt{\frac{2}{\pi }} \frac{N}{3\nu ^3}\frac{a}{a_0}+\frac{N}{3\pi ^{3/2}}\frac{a_0^G}{a_0}\frac{I(\nu )}{\nu ^3}. \end{aligned}$$
(25)

If \(U(\nu )\) admits a local minimum \(\nu _\mathrm {min}\), the width of BEC \(\sigma (t)=\nu (t)a_0\) would experience the shape oscillation with frequency \(\omega =\sqrt{U''(\nu _\mathrm {min})}\omega _0\) around \(\nu _\mathrm {min}\). An illustration for the effective potential is shown in Fig. 1 with magnetically manipulated scattering length \(a=0\) (as we will assume later), and the cases without gravity (black solid line), with the Newtonian gravity only (blue solid line), and with extra Yukawa interaction (red dashed line), which will be presented in detail in the following sections in addition to the cases with the extra power-law potential and the fat-graviton potential.

Fig. 1
figure 1

The effective potential for BEC oscillation in the ultra-cold atom experiment with magnetically manipulated scattering length \(a=0\). The gravitational potentials are illustrated for the cases without gravity (black solid line), with Newtonian gravity only (blue solid line), and with Yukawa interaction (red dashed line)

2.2 Shape oscillation frequency

To estimate the contributions of each term in the effective potential (25), one could adopt the typical value of \(a_0\sim 10^{-4}\,\mathrm {cm}\) and \(N\sim 10^6{-}10^8\), which is experimentally achievable currently, as well as \(a_0^G\) taken to be

$$\begin{aligned} a_0^G=\frac{c}{\omega _0}\left( \frac{m}{m_\mathrm {Pl}}\right) ^2=1.747\times 10^{-28}\,\mathrm {cm}\left( \frac{\omega _0}{\mathrm {Hz}}\right) ^{-1}A^2, \end{aligned}$$
(26)

where the atomic mass number A is usually around order of hundred, and the trapping frequency \(\omega _0\) could be adjusted in the range of \(50{-}10,000\) Hz. Therefore, for typical choice of scattering length \(a\approx a_0\gg a_0^G\) due to the \(m/m_\mathrm {Pl}\) suppression in the ratio \(a_0^G/a_0\sim 10^{-24}{-}10^{-22}\), the second term would dominate over the third term in (25), namely, the s-wave scattering would totally erase the trace of gravitational potential. To manifest the gravitational effect, one could tune down the scattering length a magnetically to zero [29] (Hereafter we will assume this idealistic condition for our preliminary studies) so that we can get rid of the second term in (25), and the effect of additional gravitational potential could be extracted from the frequency deviation of shape oscillation [30,31,32,33] in the absence of gravity.

2.2.1 Without gravity

In the case without gravity [27], the effective potential reads

$$\begin{aligned} U(\nu )=\frac{1}{2}\left( \nu ^2+\frac{1}{\nu ^2}\right) , \end{aligned}$$
(27)

and the corresponding EOM

$$\begin{aligned} \frac{{\text{ d }}^2\nu }{{\text{ d }}\tau ^2}+\nu -\frac{1}{\nu ^3}=0 \end{aligned}$$
(28)

is solved by

$$\begin{aligned} \nu _0(\tau )=\frac{1}{\nu _i}\sqrt{\frac{1+\nu _i^4+(\nu _i^4-1)\cos 2\tau }{2}} \end{aligned}$$
(29)

with initial condition \(\nu (\tau =0)=\nu _i\). The corresponding local minimum \(\nu _\mathrm {min}=1\) could be solved from the zeros of the effective potential \(U'(\nu )=0\), around which the oscillation frequency is exactly \(\omega (\nu _\mathrm {min})\equiv \sqrt{U''(\nu _\mathrm {min})}\omega _0=2\omega _0\) to the second order. Higher order terms in the expansion of \(\nu (\tau )\) around \(\nu _\mathrm {min}\), which should be extracted by recording the full time evolution of shape oscillation in the experiment, could be ignored in principle by controlling the initial size \(\nu _i\) near \(\nu _\mathrm {min}\).

2.2.2 With Newtonian gravity only

As for the shape oscillation with only Newtonian gravity included, the effective potential is computed as

$$\begin{aligned} U(\nu )=\frac{1}{2}\left( \nu ^2+\frac{1}{\nu ^2}\right) -\sqrt{\frac{2}{\pi }}\frac{Na_0^G}{3a_0\nu }, \end{aligned}$$
(30)

from which the oscillation frequency is found to be deviated from \(2\omega _0\) by

$$\begin{aligned} \frac{\omega }{2\omega _0}=\sqrt{1+\frac{1}{4}\left( \frac{1}{\nu _\mathrm {min}^4}-1\right) }, \end{aligned}$$
(31)

where the local minimum \(\nu _\mathrm {min}\) is the root of \(U'(\nu )=0\), namely,

$$\begin{aligned} \frac{1}{\nu _\mathrm{{min}}^4}-1=\sqrt{\frac{2}{\pi }}\frac{Na_0^G}{3a_0\nu _\mathrm {min}^3}. \end{aligned}$$
(32)

Due to the suppression factor \(a_0^G/a_0\ll 1\), deviations on both the local minimum \(\nu _\mathrm {min}\) from 1 and oscillation frequency \(\omega \) from \(2\omega _0\), respectively, are extremely small.

Nevertheless, if the oscillation frequency could be measured to a high accuracy, we are able to obtain the Newtonian gravitational constant in this BEC experiment. In fact, the combination of (31) and (32) gives rise to an expression for the Newtonian gravitational constant as

$$\begin{aligned} G=\sqrt{\frac{\pi }{2}}\frac{3\sigma _\mathrm {min}^3}{mN}(\omega ^2-4\omega _0^2), \end{aligned}$$
(33)

from which the derivative with respect to frequency change reads

(34)

To further evaluate the part of \(\mathrm {d}\sigma _\mathrm {min}/\mathrm {d}\omega \), one first rewrites (32) as

$$\begin{aligned} a_0^4=\sigma _\mathrm {min}^4+\sqrt{\frac{2}{\pi }}\frac{mN}{3\omega _0^2}G\sigma _\mathrm {min}, \end{aligned}$$
(35)

which, after taking derivative, becomes

$$\begin{aligned} 0=4\sigma _\mathrm {min}^3\frac{\mathrm {d}\sigma _\mathrm {min}}{\mathrm {d}\omega }+\sqrt{\frac{2}{\pi }}\frac{mN}{3\omega _0^2}\frac{\mathrm {d}(G\sigma _\mathrm {min})}{\mathrm {d}\omega }. \end{aligned}$$
(36)

Here the derivative term \(\mathrm {d}(G\sigma _\mathrm {min})/\mathrm {d}\omega \) could be replaced by taking the derivative of (33) after multiplied with \(\sigma _\mathrm {min}\), namely,

$$\begin{aligned} \frac{\mathrm {d}(G\sigma _\mathrm {min})}{\mathrm {d}\omega }=3\sqrt{2\pi }\frac{a_0^4}{mN}\frac{\omega \omega _0^4}{(\omega ^2-3\omega _0^2)^2}. \end{aligned}$$
(37)

Then \(\mathrm {d}\sigma _\mathrm {min}/\mathrm {d}\omega \) could be obtained from above two equations as

$$\begin{aligned} \frac{\mathrm {d}\sigma _\mathrm {min}}{\mathrm {d}\omega }=-\frac{\omega \sigma _\mathrm {min}^5}{2\omega _0^2a_0^4}. \end{aligned}$$
(38)

Now (34) could be exactly calculated as

$$\begin{aligned} \frac{1}{G}\frac{\mathrm {d}G}{\mathrm {d}\omega }=\frac{\omega ^3}{2(\omega ^2-4\omega _0^2)(\omega ^2-3\omega _0^2)} \end{aligned}$$
(39)

without explicitly finding the root of (32). Finally, by approximating \(\omega \approx 2\omega _0\) in the factors \(\omega +2\omega _0\) and \(\omega ^2-3\omega _0^2\), we arrive at a preliminary estimation for the relative error of G from the relative error of \(\omega \) as

$$\begin{aligned} \frac{\Delta G}{G}\approx \frac{\Delta \omega }{\omega }\bigg /\left( \frac{\omega }{2\omega _0}-1\right) , \end{aligned}$$
(40)

which is suppressed by the frequency deviation \(\omega /(2\omega _0)-1\).

As we will see later at the Newtonian limit of modified gravity theories, the frequency deviation \(\omega /(2\omega _0)-1\) is usually of the size \(10^{-14}\) for \(N=10^9\) just above the current achievable \(N\sim 10^6{-}10^8\), which requires the relative error on the \(\omega \) as small as \(10^{-18}\) if we want to measure \(\Delta G/G\) up to precision of \(10^{-4}\). This corresponds to a resolution of frequency measurement \(\Delta \omega \sim 10^{-18}\omega \sim 10^{-16}{-}10^{-14}\,\mathrm {Hz}\) provided that the value of \(\omega _0\) ranges from \(50{-}10000\,\mathrm {Hz}\). However, the practical resolution of frequency measurement could be relaxed since the oscillation frequency could be measured and calibrated over numerous oscillation periods for sufficiently long lifetime of the BEC state like that in space. To our knowledge, despite of the atom interferometric determination of the Newtonian gravitational constant [34, 35], our proposal for the measurement on Newtonian gravitational constant from the cold atom BEC experiment with shape oscillation has not been investigated in the literature, although it is experimentally more challenging.

On the other hand, with Newtonian gravitational constant known from other experiments, we can further probe the realms of modified gravity theories by measuring the deviation of the oscillation frequency. Note that the impact on the relative error of the oscillation frequency due to the Newtonian gravitational potential is much smaller than the frequency deviation as seen from (40), namely,

$$\begin{aligned} \frac{\Delta \omega }{\omega }\approx \frac{\Delta G}{G}\left( \frac{\omega }{2\omega _0}-1\right) \ll \left( \frac{\omega }{2\omega _0}-1\right) , \end{aligned}$$
(41)

as long as the frequency deviation is small, \(\omega \approx 2\omega _0\), and the Newtonian gravitational constant could be measured from other experiments to relatively high accuracy, \(\Delta G/G\ll 1\).

3 Shape oscillation for modified gravitational potential

3.1 Power-law potential

For an additional power-law potential (5) to the Newtonian gravity, the corresponding effective potential is obtained as

$$\begin{aligned} U(\nu )&=\frac{1}{2}\left( \nu ^2+\frac{1}{\nu ^2}\right) \nonumber \\&\quad -\sqrt{\frac{2}{\pi }}\frac{Na_0^G}{3a_0\nu }\left[ 1+\beta _k\left( \frac{1\,\mathrm {mm}/a_0}{\sqrt{2}\nu }\right) ^{k-1}\Gamma \left( \frac{3-k}{2}\right) \right] \end{aligned}$$
(42)

for \(k<3\). For \(k\ge 3\), the integral (24) is not well-defined. Similar to (31), the oscillation frequency could be found to be deviated from \(2\omega _0\) by

$$\begin{aligned} \frac{\omega }{2\omega _0}&=\sqrt{1+\frac{2-k}{4}\left( \frac{1}{\nu _\mathrm {min}^4}-1\right) +\frac{k-1}{4}\sqrt{\frac{2}{\pi }}\frac{Na_0^G}{3a_0\nu _\mathrm {min}^3}}, \end{aligned}$$
(43)

where the local minimum \(\nu _\mathrm {min}\) is determined by \(U'(\nu )=0\) as

$$\begin{aligned} \frac{1}{\nu _\mathrm {min}^4}-1=\sqrt{\frac{2}{\pi }}\frac{Na_0^G}{3a_0\nu _\mathrm {min}^3}\left[ 1+k\beta _k\left( \frac{1\,\mathrm {mm}/a_0}{\sqrt{2}\nu _\mathrm {min}}\right) ^{k-1}\Gamma \left( \frac{3-k}{2}\right) \right] . \end{aligned}$$
(44)

3.2 Yukawa interaction

For an additional Yukawa interaction (6) to the Newtonian gravity, the corresponding effective potential is obtained as

$$\begin{aligned} U(\nu )=&\frac{1}{2}\left( \nu ^2+\frac{1}{\nu ^2}\right) -\sqrt{\frac{2}{\pi }}\frac{Na_0^G}{3a_0\nu }\nonumber \\&\times \left\{ 1+\alpha \left[ 1-\sqrt{\pi }\left( \frac{\beta \nu }{\sqrt{2}}\right) \mathrm{{e}}^{\frac{\beta ^2\nu ^2}{2}}\mathrm{{Erfc}}\left( \frac{\beta \nu }{\sqrt{2}}\right) \right] \right\} , \end{aligned}$$
(45)

where \(\beta \equiv a_0/d\) and \(\mathrm {Erfc}(z)\equiv 1-\mathrm {Erf}(z)\) is the complementary error function. The oscillation frequency could be found to be deviated from \(2\omega _0\) by

$$\begin{aligned} \frac{\omega }{2\omega _0}=&\left\{ 1+\left( \frac{1}{\nu _\mathrm {min}^4}-1\right) \left( 1+\frac{\beta ^2\nu _\mathrm {min}^2}{4}\right) \right. \nonumber \\&\left. -\sqrt{\frac{2}{\pi }}\frac{Na_0^G}{a_0}\frac{3+3\alpha +\beta ^2v_\mathrm {min}^2}{12v_\mathrm {min}^3}\right\} ^{1/2}, \end{aligned}$$
(46)

where the local minimum \(\nu _\mathrm {min}\) is determined by \(U'(\nu )=0\) as

$$\begin{aligned} \frac{1}{\nu _\mathrm {min}^4}-1=&\sqrt{\frac{2}{\pi }}\frac{Na_0^G}{3a_0\nu _\mathrm {min}^3}\left[ 1+\alpha -\alpha \beta ^2\nu _\mathrm {min}^2\right. \nonumber \\&\left. +\sqrt{\frac{\pi }{2}}\alpha \beta ^3\nu _\mathrm {min}^3\mathrm {Erfc}\left( \frac{\beta \nu _\mathrm {min}}{\sqrt{2}}\right) \right] . \end{aligned}$$
(47)

3.3 Fat graviton

For an additional fat-graviton potential (8) to the Newtonian gravity, the corresponding effective potential is obtained as

$$\begin{aligned} U(\nu )=&\,\frac{1}{2}\left( \nu ^2+\frac{1}{\nu ^2}\right) -\sqrt{\frac{2}{\pi }}\frac{Na_0^G}{3a_0\nu }\nonumber \\&\times \left[ 1-\frac{1}{2\beta ^2\nu ^2}\Gamma \left( \frac{5}{3}\right) \,_2F_2\left( \frac{1}{2},\frac{5}{6};\frac{2}{3},\frac{3}{2};-\frac{1}{54\beta ^6\nu ^6}\right) \right. \nonumber \\&-\frac{1}{45\beta ^4\nu ^4}\Gamma \left( -\frac{2}{3}\right) \,_2F_2\left( \frac{5}{6},\frac{7}{6};\frac{4}{3},\frac{11}{6};-\frac{1}{54\beta ^6\nu ^6}\right) \nonumber \\&\left. -\frac{1}{56\beta ^6\nu ^6}\,_3F_3\left( 1,\frac{7}{6},\frac{3}{2};\frac{4}{3},\frac{5}{3},\frac{13}{6};-\frac{1}{56\beta ^6\nu ^6}\right) \right] , \end{aligned}$$
(48)

where \(\beta \equiv a_0/\lambda \) and \(\,_pF_q(a_1,\ldots ,a_p;b_1,\ldots ,b_q;z)\) is the generalized hypergeometric function. The oscillation frequency could be found to be deviated from \(2\omega _0\) by

(49)

where the local minimum \(\nu _\mathrm {min}\) is determined by \(U'(\nu )=0\) as

(50)

4 Experimental perspective

To maximize the effect from the gravitational potential on the effective potential (25), N should be sufficiently large to balance the suppression factor \(a_0^G/a_0\) (chosen as \(3\times 10^{-22}\) specifically in this section). Therefore, we will take N as \(10^9\), just one order above the experimentally achievable \(N\sim 10^6{-}10^8\) currently, which might be feasible in the far future in space, although it is regarded here just as a preliminary theoretical perspective. We will also take N up to \(10^{12}\), although it is not feasible in practice. The smaller number of N could be compensated by a larger ratio of \(a_0^G/a_0\) with the larger atomic mass m and the smaller trapping frequency \(\omega _0\), although the two conditions are usually difficult to be satisfied simultaneously in the experiment.

Fig. 2
figure 2

The frequency deviation in the shape oscillation with respect to the gravity-free case for the power-law potential (left panel) and fat-graviton potential (right panel) with particles number \(N=10^9\) (red solid lines) and \(N=10^{12}\) (blue solid lines) beyond the current achievable values \(N\sim 10^6{-}10^8\). The shaded regions are excluded by current experiments (left panel and \(\lambda >98\,\upmu \mathrm {m}/0.914\) in the right panel) and naturalness argument (\(\lambda <20\,\upmu \mathrm {m}/0.914\) in the right panel). The dashed lines are for negative values of \(\omega /\omega _0{-}2\)

Fig. 3
figure 3

Experimental perspective for the modified Newtonian gravitational potential with Yukawa interaction. In the first panel, with small deviation in the oscillation frequency with respect to the gravity-free case, the time evolution of the BEC radius deviates with power-law in time. The absolute value of oscillation frequency deviations with respect to the parameter space of Yukawa interaction are shown in the second and forth panels for different values of condensed particle number with the white curves indicated by the current experimental exclusion limits. There is a small parameter region in the upper left corner with flipped sign for the frequency deviation when the Yuakwa interaction range comparable or smaller than the characteristic length scale of the BEC state as also shown in the third panel. The dashed lines are for negative values of \(\omega /\omega _0{-}2\). In the fifth panel, the potential minimum \(\nu _\mathrm {min}\) is depicted with respect to \(\alpha \) for different interaction range d, which experiences cross-over transition for \(d>a_0\) (blue curves) but first-order phase transition for \(d<a_0\) (red curves). In the last panel, the absolute deviation of oscillation frequency is presented with respect to a largely extended parameter space supported by d down to \(10^{-12}\) m and \(\alpha \) up to \(10^{30}\). The white curves are borrowed from [36, 37] for comparison with experimental constraints

4.1 Power-law potential

For the power-law potential (5), we will take \(k=2\) for illustration. The deviation of the oscillation frequency (43) with respect to the gravity-free oscillation frequency \(2\omega _0\) is shown in the left panel of Fig. 2 for the particle numbers \(N=10^9\) (red line) and \(N=10^{12}\) (blue line), where the shaded regions \(\beta _2>1.3\times 10^{-3}\) [38], \(\beta _2>4.5\times 10^{-4}\) [21], and \(\beta _2>3.7\times 10^{-4}\) [39] are excluded by the existing experiments. Unfortunately, despite of the smallness of frequency deviation, it is also insensitive to the size of coefficient \(\beta _k\). Therefore, the ultra-cold BEC experiment in space cannot constrain the modified gravity with power-law gravitational potential.

4.2 Yukawa interaction

The experimental perspective for the gravitational potential with Yukawa interaction (6) are presented in Fig. 3. The deviations of the oscillation frequency (46) with respect to the gravity-free oscillation frequency \(2\omega _0\) are shown for the particle numbers \(N=10^9\) (the second panel) and \(N=10^{12}\) (the fourth panel) with respect to the current exclusion limits [38,39,40,41,42,43,44,45,46,47,48,49,50]. The frequency deviations are insensitive to both d and \(\alpha \) when \(\alpha \lesssim \mathcal {O}(10)\), which are as small as \(10^{-14}\) (the second panel) and \(10^{-11}\) (the fourth panel), respectively. For \(d\gtrsim 10^{-6}\,\mathrm {m}\), the frequency deviation grows with power-law in the range of \(\alpha \gtrsim \mathcal {O}(10)\). For \(d<10^{-6}\,\mathrm {m}\), the frequency deviation flips a sign around \(\alpha \sim \mathcal {O}(10)\), after which the absolute value of frequency deviation still grows with power-law in \(\alpha \). This could also be seen in the third panel for some illustrative values of N and d, where the dashed lines denote for the negative value of the frequency deviation \(\omega /\omega _0-2\). Note that the regime of flipped sign appears when d is smaller than the characteristic length scale \(a_0\simeq 10^{-6}\,\mathrm {m}\) of the BEC state. Except for the region of flipped sign, the frequency deviation would be insensitive to the interaction range d for given \(\alpha \), which could greatly push current exclusion limits.

To detect such a small frequency deviation, for example in the regime of \(\alpha \lesssim \mathcal {O}(10)\) with \(\omega /\omega _0-2\simeq 10^{-14}\) for \(N=10^9\) and \(\omega /\omega _0-2\simeq 10^{-11}\) for \(N=10^{12}\) , we need the resolution of frequency in the measurement as

$$\begin{aligned} N=10^{9}&: \Delta \omega \simeq 10^{-14}\omega \simeq (1-100)\,\mathrm {pHz}, \end{aligned}$$
(51)
$$\begin{aligned} N=10^{12}&: \Delta \omega \simeq 10^{-11}\omega \simeq (1-100)\,\mathrm {nHz}, \end{aligned}$$
(52)

where the trapping frequency is usually in the range of 50–10,000 Hz. The above requirement for the resolution of frequency could be either relaxed in the regime where the absolute value of frequency deviation \(|\omega /\omega _0-2|\) grows with the power-law in \(\alpha \), or enhanced in the regime where the frequency deviation \(\omega /\omega _0-2\) flips a sign.

However, the actual requirement for the resolution of frequency measurement could be further relaxed universally in the whole parameter space of d and \(\alpha \) by virtue of the long lifetime of ultra-cold BEC in space. In fact, the difference between the measured time-evolution \(\nu (\tau )\) of oscillating BEC and the gravity-free solution \(\nu _0(\tau )\) grows with power-law in time \(\tau \) as seen in the first panel of Fig. 3. Furthermore, for trapping frequency in the range of 50–10,000 Hz, the measurement of BEC oscillations up to \(\tau \,{\sim }\,10^4\) corresponds to a condensate lifetime of 1–100 s that is exactly comparable to the free-fall time achieved with the help of microgravity environment. Therefore, even if the current resolution of frequency-measurement is not high enough to determine the oscillation frequency within a single oscillation period, the long-time measurement of oscillating BEC in space allows us to calibrate the measured frequency from previous oscillation period to higher precision over the next numerous oscillation periods. In our specific case, the measurement of oscillations up to \(\tau \,{\sim }\,10^4\) corresponds to \(n\,{\sim }\,10^3\) periods, which could further relax the resolution of frequency measurement down to

$$\begin{aligned} N=10^9&: \Delta \omega _{nT}\sim n\Delta \omega \simeq (1-100)\,\mathrm {nHz}, \end{aligned}$$
(53)
$$\begin{aligned} N=10^{12}&: \Delta \omega _{nT}\sim n\Delta \omega \simeq (1-100)\,\upmu \mathrm {Hz}. \end{aligned}$$
(54)

As matter of a preliminary investigation, we also study the Yukawa interaction in a largely extended parameter space supported by the interaction range d down to \(10^{-12}\) m and strength \(\alpha \) up to \(10^{30}\). For the cases of \(d>a_0\) and \(d<a_0\), the potential minimum \(\nu _\mathrm {min}\) experiences a cross-over and first-order phase transitions with increasing \(\alpha \) as shown with blue and red curves in the fifth panel of Fig. 3, respectively. The absolute deviation of oscillation frequency in this extended parameter space is depicted in the last panel of Fig. 3, which is relatively small for the blue and purple regions so that the impact on the relative error of the oscillation frequency from the Newtonian gravitational potential is negligible compared to the measured deviation of oscillation frequency according to (41). Therefore, this parameter regions could be promising to be ruled out by future ultra-cold atom experiments. The rest of the parameter regions, although largely affected by the Newtonian potential, were already ruled out mostly by the current constraints reviewed in [36, 37] and references therein. Furthermore, the idealistic condition we assume for neglecting the second term in (25) could be relaxed for large \(\alpha \) in the extended parameter space, and the further investigation for the full competitions of all three terms in (25) will be reserved for future work.

4.3 Fat graviton

For the fat-graviton potential (8), the deviation of oscillation frequency (49) with respect to the gravity-free oscillation frequency \(2\omega _0\) is shown in the right panel of Fig. 2 for the particle numbers \(N=10^9\) (red line) and \(N=10^{12}\) (blue line), where the dashed lines denote for the negative value of the frequency deviation \(\omega /\omega _0-2\), and the shaded region \(\lambda >98\,\upmu \mathrm {m}/0.914\) is excluded by the experiments [21, 41]. The other shaded region \(\lambda <20\,\upmu \mathrm {m}/0.914\) is argued by naturalness in [25]. Similar to case with power-law potential, the frequency deviation is also insensitive to \(\lambda \) regardless of the flipped sign around \(\lambda \simeq 2\times 10^{-6}\,\mathrm {m}\), which is also roughly the characteristic length scale \(a_0\simeq 10^{-6}\,\mathrm {m}\) of BEC state. Therefore, the ultra-cold BEC experiment in space also cannot constrain the modified gravity with fat-graviton potential in the regime of interest \(20\,\upmu \mathrm {m}/0.914<\lambda <98\,\upmu \mathrm {m}/0.914\).

5 Conclusion and discussions

Contrary to the native expectation that gravity cannot be further constrained at small scales since it is weakest among all other fundamental/quintessential interactions, the ultra-cold BEC experiments in space with microgravity environment could explore the parameter space in the modified Newtonian potentials because this experiment could significantly extend the observation time so as to precisely measure the microscopic properties of gravity via its macroscopic manifestation on the BEC state with shape oscillation. Newtonian gravity and several modified gravity theories are examined in detail. The results with Yukawa interaction are optimistic in probing the parameter space greatly beyond the current experimental exclusion limits. However, the results with the Newtonian gravity, power-law potential and fat-graviton potential are quite pessimistic. Several comments are given below:

Firstly, our proposal of measurement on Newtonian gravitational constant is new to our knowledge in the sense of the use of cold-atom BEC with shape oscillation, although it is experimentally more challenging than the current use of atom interferometer.

Secondly, the long lifetime of BEC state in space buys us more time to measure and calibrate the oscillation frequency over numerous periods so that the actual resolution of frequency measurement could be relaxed compared to the frequency resolution in a single period.

Thirdly, our estimation for the frequency resolution is quite conservative, which could be further relaxed for \(\alpha \gtrsim \mathcal {O}(10)\) in the regime of \(d\sim 10^{-6}{-}10^{-5}\,\mathrm {m}\) where frequency deviation grows with power-law in \(\alpha \). This greatly improves the current exclusion limits.