1 Introduction

The success of general relativity (GR) is widely acknowledged, as the theory has been confirmed by many observations [1, 2] and has been able to predict astrophysical objects like black holes [3]. However, despite its predictive success and elegance, GR does not come flawless. In fact, in order to match observations on galactic and cosmological scales, the theory requires two unknown components, namely dark matter (to explain galactic rotation curves and a dynamical mass on galaxy clusters) and dark energy (to explain the late-time accelerated expansion of the Universe). The lack of detection of these exotic components is the most important challenge faced by GR [4,5,6,7]. Besides, the theory also faces other limitations such as issues with unification with high energy physics [8, 9] and the existence of space-time singularities [10]. These drawbacks are the main motivations for modified theories of gravity.

There are many possibilities for such modified theories of gravity; some proposals consist in a generalization of the Einstein–Hilbert gravitational action [11], other approaches include minimally or non-minimally coupled scalar fields [12, 13] or additional geometric ingredients [14], while other proposals treat the physical constants as dynamical quantities [15, 16]. Among the simplest modified theories of gravity are the so-called f(R) theories which replace the Ricci scalar by a generic function thereof in the action functional [17, 18]. One specific proposal of such theories was initially introduced to deal with the initial conditions problem of the standard Hot Big Bang model, namely the so-called Starobinsky model [19]. This proposal has many interesting virtues as it is in excellent agreement with the most recent data from Planck mission [20]. Besides, it has been shown that requiring f(R) models to be regular at \(R = 0\) leads to a behavior compatible with an effective cosmological constant in a sufficiently curved spacetime which disappears in a flat spacetime [21].

Another promising avenue to solve the aforementioned problematic aspects of GR relies on the extension of f(R) theories with a non-minimal coupling between matter and curvature [22] (hereafter NMC model for short). This approach has a plethora of interesting theoretical and observational implications: It allows for a mimicking effect of dark matter effects at the scale of galaxies [23] and galaxy clusters [24], it accounts for the accelerated expansion of the Universe and provides (through mimicking) a viable unification of dark energy and dark matter [25], and it is compatible with Planck’s inflation data [26], gravitational waves measurements [27], and the modified virial theorem from the spherically relaxed Abell 586 cluster [28]. Besides, it has been shown [29] that this model yields a correction to the Newtonian potential, offering therefore interesting possibilities to test this proposal in the weak field regime.

As most interesting phenomena allowing to test modified theories of gravity, and to constraint their parameters, occur in N-body systems, one needs to implement statistical methods in these modified theories. The kinetic theory has been recently formulated [30] in the context of the NMC alternative gravity model. It has been shown that an extra force term arises in the Boltzmann equation, due to the geodesic deviation force. Later on, this kinetic treatment has been employed, by one of the present authors, to study the so-called Jeans instability [31], i.e., the phenomenon responsible for the collapse of a gravitationally bound system, e.g. the interstellar gas, ultimately leading to star formation. The effects of the NMC alternative gravity model on the Jeans criterion for gravitational instability have been studied, offering new possibilities to test the model and to constraint its parameters.

In this paper, we generalize this kinetic treatment to the quantum regime. For that purpose, we formulate a quantum kinetic approach to the NMC generalization of the Schrödinger–Poisson (SP) model [sometimes called Schrödinger–Newton model] (see also [32,33,34] for numerical simulations of NMC extensions of the SP model). The SP model relies on the combination of the Schrödinger equation and the Poisson equation describing the self-potential. Historically, the SP model was first advocated by Diósi [35] and Penrose [36] as a simple quasi-classical model to introduce quantum effects in gravitational problems. It describes quantum matter confined by gravitational fields and, as such, finds many applications in astrophysics; it describes (hypothesized) boson stars [37,38,39] (which could be a source of ‘exotic’ Laser Interferometer Gravitational-Wave Observatory (LIGO) detections [40]) and it is a central ingredient of scalar field dark matter models [41,42,43,44]. Interestingly, the same model applies equally well to a variety of other systems, such as quantum plasmas [45,46,47] and atomic gases in magneto-optical traps [48,49,50]; a feature that may inspire future laboratory experiments to test alternative gravity theories using condensed-matter gravity analogs [51,52,53].

By employing the Wigner-Moyal procedure [54, 55], allowing to represent quantum states in a classical phase space, we formulate here a quantum kinetic approach generalizing the classical approach to matter confined by gravitational interactions in the context of the NMC alternative gravity model [31], and study the NMC effect on the Jeans instability criterion. This allows us to study the rich interplay between NMC effects, quantum effects, and kinetic effects. The rest of the paper progresses in the following fashion: For completeness, we lay out in Sect. 2 the theoretical background of the NMC gravity model; we present in Sect. 3 the corresponding Wigner kinetic formulation. In Sect. 4, we establish the dispersion relation for Jeans instability and discuss its limits for negligible (zero-temperature) and dominant kinetic effects. In Sect. 5, we analyze in more details specific models and confront the model with data of Bok globules in Sect. 6. We conclude in Sect. 7.

2 Non-minimal matter-curvature coupling model

The NMC gravity model is defined through its action functional [22]

$$\begin{aligned} S=\int d^{4} x \sqrt{-g}\left[ \kappa f_{1}(R)+f_{2}(R) \mathcal {L}\right] , \end{aligned}$$
(1)

where \(\kappa {:}{=} c^4/ 16 \pi G\), \(\mathcal {L}\) is the Lagrangian density of matter fields, and \(f_1(R)\) and \(f_2(R)\) are arbitrary functions of the curvature scalar R. By varying the action (1) with respect to the metric \(g_{\mu \nu }\), one obtains the metric field equations

$$\begin{aligned} \varTheta G_{\mu \nu }=\frac{1}{2 \kappa } f_{2}(R) T_{\mu \nu }+\varDelta _{\mu \nu } \varTheta +\frac{1}{2} g_{\mu \nu }\left[ f_{1}(R)-R \varTheta \right] , \end{aligned}$$
(2)

where \(G_{\mu \nu }\) is the Einstein tensor and \(\varTheta {:}{=}\left( f_{1}^{\prime }(R)+\frac{f_{2}^{\prime }(R) \mathcal {L}}{\kappa }\right) \), with the primes denoting derivatives with respect to the curvature scalar \(f_{i}^{\prime }(R) \equiv d f_{i}(R) / d R\) and \(\varDelta _{\mu v} \equiv \nabla _{\mu } \nabla _{v}-g_{\mu v} \square \). Note that for the particular choice \(f_1(R)=R\) and \(f_2(R)=1\), one recovers GR. The trace of the metric field equations reads as

$$\begin{aligned} \varTheta R-2 f_{1}(R)=\frac{1}{2 \kappa } f_{2}(R) T-3 \square \varTheta . \end{aligned}$$
(3)

Hence, the energy-momentum tensor is not covariantly conserved

$$\begin{aligned} \nabla _{\mu } T^{\mu v}=\left( g^{\mu v} \mathcal {L}-T^{\mu v}\right) \nabla _{\mu } \ln f_{2}(R). \end{aligned}$$
(4)

This is in fact one of the sticking features of the NMC model. It implies that for a perfect fluid, with \(T_{\mu \nu }=( \rho + P) u^{\mu } u^{v}-P g_{\mu v}\), a test particle will not follow geodesic lines, due to the presence of an extra force term in the geodesics equation. That is [22],

$$\begin{aligned} \frac{d u^{\alpha }}{d s}+\varGamma _{\mu v}^{\alpha } u^{\mu } u^{v}=f^{\alpha }, \end{aligned}$$
(5)

where \(u^{\mu }\) is the 4-velocity of the particle and the extra-force per unit mass \(f^{\alpha }\) reads as

$$\begin{aligned} f^{\alpha }=\frac{1}{\rho +P}\left[ \frac{f_{2}^{\prime }(R)}{f_{2}(R)}\left( \mathcal {L}_{m}-P\right) \nabla _{v} R-\nabla _{v} P\right] V^{\alpha v}, \end{aligned}$$
(6)

\(V^{\alpha \nu }=g^{\alpha \nu }+u^{\alpha } u^{\nu }\) being the projection operator. In the case of dust \(P=0\) where the Lagrangian has a clear choice \(\mathcal {L}=- \rho c^2\), in the Newtonian regime, the geodesic equation reads as [29]

$$\begin{aligned} \frac{d^{2} x^{i}}{d t^{2}}=\partial _{i}\left[ \frac{g_{t t}+1}{2}-\ln \left| f_{2}\right| \right] -\frac{\partial _{i} P}{\rho +P }, \end{aligned}$$
(7)

from which one may define a NMC potential [29] \(\varPhi _c {:}{=} \ln |f_2 |\). We note that we shall use units where \(c=1\) throughout the rest of the paper, until Sect. 6, where we shall recover the International System of Units.

As shown in [31], for the obvious choice for the metric field

$$\begin{aligned} g_{\mu \nu }={\text {diag}}(-1-2 \varPhi , 1-2 \varPsi , 1-2 \varPsi , 1-2 \varPsi ) \end{aligned}$$
(8)

where both \(|\varPhi |\) and \(|\varPsi |\ll 1\) correspond to the Bardeen gauge invariant potentials, the 00-components of the Ricci tensor are \(\delta R_{00}=\nabla ^{2} \varPhi \) and the scalar curvature \(\delta R = 2 \nabla ^{2}(2 \varPsi -\varPhi )\). In this case, the two potentials can be decoupled into two Poisson-like equations [31]

$$\begin{aligned} \left\{ \begin{array}{l} \left( 3 \alpha \nabla ^{4}-\nabla ^{2}\right) \varPhi =\left( \alpha \gamma -\frac{\beta }{2}\right) \nabla ^{2} \rho -\frac{\gamma }{4} \rho \\ \left( 3 \alpha \nabla ^{4}-\nabla ^{2}\right) \varPsi =\left( \frac{\alpha \gamma +\beta }{2}\right) \nabla ^{2} \rho -\frac{\gamma }{4} \rho \end{array}\right. \end{aligned}$$
(9)

where \(\rho \) is the mass density (the term \(\nabla ^2 \rho \) originating from the non-minimal coupling) and the model parameters are given by \(\alpha {:}{=}f_1^{''}(0)/f_1'(0)\), \(\beta {:}{=}f_2'(0)/(\kappa f_1'(0))\) and \(\gamma {:}{=}f_2(0)/(\kappa f_1'(0))\). Note that by setting \(f_2(R)=0\), one recovers f(R) theory while by further setting \(f_1(R)=R\) (equivalently, \(\alpha =0\)), one recovers \(\varPsi = \varPhi \) and the standard Poisson equation is recovered.

3 The Wigner equation in the Newtonian limit of the NMC

In a kinetic approach, one considers the phase space spanned by the space and momentum coordinates \((\textbf{r},\textbf{p})\). A state of the given system is characterized by the one-particle distribution function \(f(\textbf{r},\textbf{p},t)\), such that \(f(\textbf{r},\textbf{p},t)d^3\textbf{r}d^3\textbf{p}\) gives the number of particles in the volume element \(d^3\textbf{r}\) about the position \(\textbf{r}\) and with momentum in the range \(d^3\textbf{p}\) about \(\textbf{p}\). The space-time evolution of the one-particle distribution function f in the phase space is given by the Boltzmann equation. For a collisionless medium in the presence of a gravitational potential \(\varPhi \) and the coupling potential \(\varPhi _c {:}{=} \ln |f_2 |\) arising from the NMC coupling, the Boltzmann equation reads as [30, 31]

$$\begin{aligned} \frac{\partial f}{\partial t}+ \frac{\textbf{p}}{m} \cdot \nabla f- m \nabla \left( \varPhi +\varPhi _{c}\right) \cdot \frac{\partial f}{\partial \textbf{p}}=0, \end{aligned}$$
(10)

where m is the mass of the self-gravitating particles. A quantum version of the collisionless Boltzmann equation (10) can be constructed by employing Wigner functions and the so-called Wigner–Moyal procedure [54, 55] which allows representing quantum mechanics into classical phase space. For that purpose, we consider the NMC extension of the Schrödinger equation, that is

$$\begin{aligned} i \hbar \frac{\partial \psi (\textbf{r},t)}{\partial t}=-\frac{\hbar ^{2}}{2 m} \nabla ^2 \psi (\textbf{r},t) + V \psi (\textbf{r},t), \end{aligned}$$
(11)

where we have defined a generic potential \(V \equiv m (\varPhi + \varPhi _c)\) as the sum of the gravitational potential \(\varPhi \) and the nonminimally coupled potential \(\varPhi _c\). To obtain the quantum counterpart of the Boltzmann equation, we follow the so-called Wigner–Moyal procedure [54, 55] (see also [56,57,58]). To do this, let us first define the Wigner distribution functionFootnote 1 associated with the wave-function \(\psi \) as follows

$$\begin{aligned} W(\textbf{r}, \textbf{p}; t)= & {} \int \frac{d \textbf{y}}{(2 \pi \hbar )^3} \exp (i \textbf{p}. \textbf{y} / \hbar ^3)\psi ^{*}(\textbf{r}+\textbf{y} / 2, t) \nonumber \\{} & {} \times \psi (\textbf{r}-\textbf{y} / 2, t), \end{aligned}$$
(12)

where \(\textbf{p} \equiv m \textbf{v}\) is the particle momentum. The Wigner function (12) is simply the Fourier transform of the auto-correlation function corresponding to the wave-function \(\psi \). It is standard practice, in purely gravitational problems, to normalize the Wigner function to the mass density. That is

$$\begin{aligned} \int d \textbf{p} W(\textbf{r}, \textbf{p}; t)=\left| \psi (\textbf{r}, t)\right| ^{2} = \rho (\textbf{r},t), \end{aligned}$$
(13)

where \(\rho (\textbf{r},t)\) denotes the mass density. Following the Wigner–Moyal procedure, one may write the Schrödinger equation (11) in the form of a kinetic equation as follows (see for instance [56] for detailed calculations)

$$\begin{aligned} \left( \frac{\partial }{\partial t}+\frac{\textbf{p}}{m} \frac{\partial }{\partial \textbf{r}}\right) W-\frac{2 (\varPhi + \varPhi _c)}{\hbar } \sin \left( \frac{\hbar }{2} \frac{\overleftarrow{\partial }}{\partial \textbf{r}} \frac{\overrightarrow{\partial }}{\partial \textbf{p}}\right) W=0, \end{aligned}$$
(14)

where the sine operator is defined in terms of its Taylor expansion and the arrows indicate the sense according to which the operators act; derivatives with respect to the momentum act forward on the Wigner function while derivatives with respect to the position act backward on the potential. Note that in the limit \(\hbar \rightarrow 0\), Eq. (14) reduces to the NMC Boltzmann equation (10). This formal limit corresponds to the case of a slowly varying potential, changing significantly only over a length-scale much larger than the de Broglie wavelength, such that \(\sin \Lambda \simeq \Lambda \).

In what follows, we shall use the Wigner kinetic equation (14) together with the NMC Poisson-like equations (9) to address the collective behavior of a self-gravitating medium and the associated Jeans instability. This quantum kinetic treatment allows to investigate the interplay between NMC effects, quantum effects, and kinetic effects, covering therefore a wide range of situations.

4 Jeans instability

Using the Wigner equation (14), we follow here the standard procedure yielding the dispersion relation (see for instance [31]). We restrict ourselves to the case of linear perturbations and write the relevant quantities as a small perturbation, in the form of plane waves, around homogeneous and stationary quantities (given by \(W_0 (\textbf{p})\), \(\varPhi _0=0\), and \(\varPhi _{c ~0}=0\)). That is

$$\begin{aligned} W({\textbf {r}}, {\textbf {p}} ; t)&\equiv W_{0}({\textbf {p}})+\tilde{W} \exp [i({\textbf {k}}. {\textbf {r}}-\omega t)] \nonumber \\ \varPhi ({\textbf {r}}, t)&\equiv \varPhi _0 + \tilde{\varPhi } \exp [i({\textbf {k}}.{\textbf {r}} -\omega t)], \nonumber \\ \varPhi _c({\textbf {r}}, t)&\equiv \varPhi _{c0}+ \tilde{\varPhi } \exp [i({\textbf {k}}.{\textbf {r}} -\omega t)]. \end{aligned}$$
(15)

Upon linearizing the Wigner equation (14) and making use of the so-called Jeans swindle, i.e., by considering that the potentials are sourced only by the density perturbations and not by the density background,Footnote 2 we obtain the following dispersion relation

$$\begin{aligned}{} & {} 1+\frac{(\alpha \gamma -\beta / 2) k^{2}+\gamma / 4}{k^{2}+3 \alpha k^{4}} \frac{m}{\hbar } \nonumber \\{} & {} \times \int _{- \infty }^{\infty } \frac{ f_0(p+\hbar k/2)-f_0(p-\hbar k/2)}{p k/m-\omega } d {p}=0, \end{aligned}$$
(16)

where we have defined \(f_0(p)\) as the one-dimensional projected (marginal) distribution along the axis parallel to the wave-vector \(\textbf{k}\). That is,

$$\begin{aligned} f_0(p)=\iint W_{0}\left( p, \mathbf {p_{\perp }} \right) {d \mathbf {p_{\perp }}}, \end{aligned}$$
(17)

where p and \(\mathbf {p_{\perp }}\) stand for the parallel and perpendicular components of the momentum respectively, i.e.,

$$\begin{aligned} \textbf{p}=p \frac{\textbf{k}}{k}+\textbf{p}_{\perp }. \end{aligned}$$
(18)

It may be noted that, in the limit \(\hbar \rightarrow 0\), the dispersion relation (16) reduces to the classical one [31]

$$\begin{aligned} 1+\frac{(\alpha \gamma -\beta / 2) k^{2}+\gamma / 4}{k^{2}+3 \alpha k^{4}} \frac{1}{m} \int _{-\infty }^{\infty } \frac{ {\partial f_0}/{\partial {p}}}{{p}- m\omega /k} d {p}=0. \end{aligned}$$
(19)

The dispersion relation (16) is very general and contains both quantum and kinetic effects. It may be interesting to study particular limits where kinetic effects are negligible or dominant.

4.1 The zero temperature limit

It may be interesting to analyze the zero-temperature limit (negligible kinetic effects). In this limiting case, there is no dispersion in the momentum (velocity) space, and the equilibrium distribution \(f_0\) shrinks to a Dirac delta, i.e.,

$$\begin{aligned} f \equiv \rho _0 \delta (p). \end{aligned}$$
(20)

In this case, the integral dispersion relation (16) reads as

$$\begin{aligned}{} & {} 1- \frac{(\alpha \gamma -\beta / 2) k^{2}+\gamma / 4}{k^{2}+3 \alpha k^{4}} \frac{m \rho _0}{\hbar } \nonumber \\{} & {} \quad \left[ \frac{1}{\omega +\hbar k^2 /2m}-\frac{1}{\omega -\hbar k^2 /2m} \right] =0. \end{aligned}$$
(21)

After some simple algebraic manipulations, Eq. (21) can be written as

$$\begin{aligned} \omega ^2=-\frac{(\alpha \gamma -\beta / 2) k^{2}+\gamma / 4}{1+3 \alpha k^{2}} \rho _0 + \frac{\hbar ^2 k^4}{4m^2}. \end{aligned}$$
(22)

The latter equation shows that quantum effects (i.e., quantum pressure) act against gravitational instability, preventing gravitational collapse to occur. This may be made more transparent by observing that, in the limit corresponding to GR, i.e., \(\alpha =\beta =0\) and \(\gamma =1/\kappa \), Eq. (22) reduces to the well-known dispersion relation (see for instance [65, 66])

$$\begin{aligned} \omega ^2=- \varOmega _J^2 + \frac{\hbar ^2 k^4}{4m^2}, \end{aligned}$$
(23)

where \(\varOmega _J=\sqrt{4 \pi G \rho _0}\) is the so-called Jeans frequency.

Note that for \(\omega ^2>0\), the frequency is real and the perturbation behaves as \({\text {e}}^{-i \omega t}\), i.e., harmonic waves, while for \(\omega ^2<0\), the frequency is imaginary (\(\omega = i \gamma \)) and the perturbation behaves as \({\text {e}}^{\gamma t}\), i.e., it evolves exponentially with time with a rate \(\gamma \) (with \(\gamma = \pm \sqrt{- \omega ^2}\)). There is a growing mode (\(\gamma >0\)) and a decaying mode (\(\gamma <0\)). The growing mode is responsible for Jeans instability. The sign of \(\omega ^2\) in (22) determines the threshold value of k separating between an oscillatory regime (\(\omega \) real) and exponential growth or instability (\(\omega \) imaginary). Setting \(\omega =0\), Eq. (22) can be solved for \(k^2\). By defining \(\xi =m^2 \rho _0/\hbar ^2\), the general solution, provided that \(\alpha \ne 0\), can be written as

$$\begin{aligned} k^2= -\frac{1}{9\alpha }\left[ 1+\varDelta +\frac{1-18 \alpha \xi (\beta -2 \alpha \gamma )}{\varDelta }\right] ~, \end{aligned}$$
(24)

where

$$\begin{aligned} \varDelta ^3= & {} 1-\alpha \xi \left( 27\beta +\frac{135}{2}\alpha \gamma \right) \\{} & {} +\sqrt{\left( 1-\alpha \xi \left( 27\beta +\frac{135}{2}\alpha \gamma \right) \right) ^2-\left( 1-18\alpha \xi (\beta -2\alpha \gamma ) \right) ^3}. \end{aligned}$$

This is the critical wave number separating between oscillatory and unstable modes. Perturbations characterized by a wave number smaller than the critical wave number (24) are unstable.

Before closing this subsection, it is instructive to note that, although using the language of the kinetic theory, kinetic effects have not been considered so far, since the equilibrium distribution \(f_0(p)\) has been identified with a Dirac delta (20). This limit can be studied in a more straightforward way by relying on the quantum hydrodynamic model (QHM). To show that, let us go back to the Schrödinger equation (11) and write the wave-function in polar form, that is

$$\begin{aligned} \psi (\textbf{r},t)= \sqrt{\rho (\textbf{r},t)} \exp \left( i S(\textbf{r},t) / \hbar \right) , \end{aligned}$$
(25)

where \(\rho = \left| \psi \right| ^{2}\) is the density and

$$\begin{aligned} S(\textbf{r},t)=-i \frac{\hbar }{2} \ln \left( \frac{\psi }{\psi ^{*}}\right) \end{aligned}$$
(26)

is the phase of the wave-function. Following the Madelung-Bohm procedure [67, 68], one may define a fluid velocity field such that

$$\begin{aligned} \quad \textbf{u} =\frac{1}{m} \nabla S, \end{aligned}$$
(27)

which ensures that the current is given by

$$\begin{aligned} \textbf{j} =\frac{\hbar }{2 m i}\left[ \psi ^{*}(\nabla \psi )-\psi \left( \nabla \psi ^{*}\right) \right] = \rho \textbf{u} . \end{aligned}$$
(28)

Substituting the wave-function (25) into the Schrödinger equation (11) and splitting apart the real and imaginary parts, one obtains

$$\begin{aligned} \frac{\partial \rho }{\partial t}+\nabla \left( \rho \textbf{u}\right)&=0, \nonumber \\ m \left( \frac{\partial \textbf{u}}{\partial t}+\textbf{u} \nabla \textbf{u} \right)&=-{\nabla \varPhi } -{\nabla \varPhi _c} - \nabla Q, \end{aligned}$$
(29)

where

$$\begin{aligned} Q \equiv -\frac{\hbar ^{2}}{2 m} \left( \frac{\nabla ^{2}\sqrt{\rho }}{\sqrt{\rho }}\right) =-\frac{\hbar ^{2}}{4 m}\left[ \frac{\nabla ^2 \rho }{\rho }-\frac{1}{2} \frac{(\nabla \rho )^{2}}{\rho ^{2}}\right] \end{aligned}$$
(30)

is the so-called quantum Bohm potential. The first equation in Eq. (29) is the continuity equation while the second one is the quantum Euler equation. The set of equations (29) may be employed following standard lines (linearization and decomposition in Fourier modes) to obtain the dispersion relation (22). Through these lines, it is more explicit that the term \(\hbar ^2 k^4/4m^2\) in Eq. (22) arises from the quantum pressure force \(-\nabla Q\) which acts against gravity, stabilizing the self-gravitating instability for small wavelengths.

4.2 Kinetic effects

To analyze the effect of the NMC gravity model in astrophysical scenarios, it may be necessary to account for purely kinetic effects, i.e., finite-temperature corrections to the distribution function. For that purpose, we go back to the general dispersion relation (16), by considering a general distribution \(f_0(p)\) other than the Dirac delta (20). We assume an even distribution \(f_0(p)\), that is

$$\begin{aligned} \int _{- \infty }^{\infty } p f_{0}(p) \textrm{d} p=0, \end{aligned}$$
(31)

which is characteristic for equilibrium and nearly equilibrium situations, and consider small quantum effects (\(\hbar k/2 \ll p\)). In this case, one may Taylor expand the functions \(f_0({p+ \hbar k/2})\) and \(f_0({p- \hbar k/2})\) and write the integral in the dispersion relation (16) as

$$\begin{aligned}{} & {} \int _{- \infty }^{\infty } \frac{f_0({p+ \hbar k/2})-f_0({p- \hbar k/2})}{p k/m- \omega } d {p} \nonumber \\{} & {} \approx \hbar k \int _{- \infty }^{\infty } \frac{f_0'(p) dp}{p k/m- \omega } + \frac{\hbar ^3 k^3}{24} \int _{- \infty }^{\infty } \frac{f_0''(p) dp}{p k/m- \omega }. \end{aligned}$$
(32)

We consider the limit of long wavelengths (\(p \gg m \omega / k\)), so that the singularity in the denominator is avoided. By performing a Taylor expansion of \((p-m \omega / k)^{-1}\) and keeping only terms with small values of k (long wavelength), one obtains the following dispersion relation

$$\begin{aligned} \omega ^{2} \simeq -\frac{(\alpha \gamma -\beta / 2) k^{2}+\gamma / 4}{1+3 \alpha k^{2}} \rho _0 +3 \frac{k_B T}{m} k^{2}+\frac{\hbar ^{2} k^{4}}{4 m^{2}}, \end{aligned}$$
(33)

where we have defined the temperature in the kinetic sense, i.e.,

$$\begin{aligned} {k_B T} \equiv \frac{\int p^{2}/m f_{0}(p) \textrm{d} p}{\int f_{0}(p) \textrm{d} p}, \end{aligned}$$
(34)

\(k_B\) being the Boltzmann constant. The dispersion relation (33) generalizes Eq. (22) by incorporating the effect of a non-vanishing velocity dispersion (thermal corrections). It remains valid as long as quantum corrections are small. It shows that, in addition to the effect of quantum pressure, thermal corrections also act against the instability, preventing the gravitational collapse from occurring. In this case, the solutions for Eq. (33) when \(\omega =0\) (separating between oscillatory and unstable modes) are more elaborated, and can be cast as

$$\begin{aligned} k^2= & {} -\frac{1}{9\alpha }\left[ (1+36\alpha \xi \epsilon )-\frac{A}{\left( B+\sqrt{B^2+A^3}\right) ^{1/3}}\right. \nonumber \\{} & {} +\left. \left( B+\sqrt{B^2+A^3}\right) ^{1/3}\right] ~, \end{aligned}$$
(35)

where we have defined \(\xi {:}{=}\frac{m^2 \rho _0}{\hbar ^2}\), as before, and \(\epsilon {:}{=}\frac{k_B T}{m\rho _0}\), \(A{:}{=}18\alpha \xi \left( (\beta -2\alpha \gamma )+6\epsilon \right) -1-72\xi \alpha \epsilon \xi -1296\xi ^2\alpha ^2\epsilon ^2\) and \(B{:}{=}1-\alpha \xi (27\beta +\frac{135}{2}\alpha \gamma )-54\xi \epsilon \alpha (1+18\xi \alpha (\beta -2\alpha \gamma ))-1944\xi ^2\alpha ^2\epsilon ^2+46656\xi ^3\alpha ^3\epsilon ^3\). Note that when the nonzero temperature corrections are turned off, we get a smooth transition to the previous scenario (Sect. 4.1), namely \((B+\sqrt{B^2+A^3})=\varDelta ^3\), and \(A\rightarrow -1+18 \alpha \xi (\beta -2 \alpha \gamma )\). Note also that in the limit \(\hbar \rightarrow 0\), one recovers the results of [31] (the classical kinetic regime).

5 Specific gravity models

In this section, we shall analyze some relevant particular situations, such as f(R), pure non-minimal coupling, and \(2\alpha \gamma = \beta \) models, likewise Ref. [31] in the classical context of Jeans instabilities for this alternative gravity model. These models allow for understanding the differences between quantum corrections with and without nonzero temperature kinetic effects.

5.1 f(R) theories

In this case, we have \(\beta =0\) and \(\gamma =1/\kappa \), which for \(\alpha \ne 0\) in the zero temperature limit yields:

$$\begin{aligned} k^2=-\frac{1}{9 \alpha }\left[ 1+\varDelta _1+\frac{\left( 1+36\xi \alpha ^2/\kappa \right) }{\varDelta _1}\right] ~, \end{aligned}$$
(36)

where \(\varDelta _1^3=1-\frac{135\xi \alpha ^2/\kappa }{2}+\sqrt{\left( 1-\frac{135\xi \alpha ^2/\kappa }{2}\right) ^2-\left( 1+36\xi \alpha ^2/\kappa \right) ^3}\). On the other hand, for \(\alpha =0\) which is the case of GR, we have solved Eq. (22) and we get two constant solutions \(k^2_{\pm }=\pm \sqrt{\frac{\xi }{\kappa }}\), from which only the positive one is physical. If \(\alpha >0\), then \(k^2<0\) and imaginary solutions for the modified Jeans mass are found. On the other hand, \(\alpha <0\) yields positive solutions, however, we should recall that the condition \(f_1''>0\) is demanded for avoiding Dolgov-Kawasacki instabilities in f(R) theories, which implies that the denominator of the definition of the \(\alpha \) parameter should be negative, \(f_1'(0)<0\). In particular, if we take the limit \(\alpha \rightarrow -\infty \), we find the behavior of \(\tilde{k}_J^2\) as a function of \(\xi /\kappa \), which signals the characteristics of the system under analysis, is as shown in Fig. 1.

Fig. 1
figure 1

Behavior of the squared modified Jeans wavenumber in f(R) theories as a function of the system parameter \(\xi \kappa \) in the limit \(\alpha \rightarrow -\infty \) for different values of \(\epsilon \)

Full size image

When kinetic corrections are included, for \(\alpha =0\), we have:

$$\begin{aligned} k^2_{\pm }=-6\epsilon \xi \pm \sqrt{(6\epsilon \xi )^2+\xi /\kappa }~, \end{aligned}$$
(37)

where the physical solutions correspond to \(k^2_+{:}{=}-6\epsilon + \sqrt{(6\epsilon )^2+\xi /\kappa }\), provided that \(\xi /\kappa > 0\). For \(\alpha \ne 0\), we have:

$$\begin{aligned} k^2= & {} -\frac{1}{9\alpha }\left[ 1+36\alpha \xi \epsilon -\frac{\bar{A}}{\left( \bar{B}+\sqrt{\bar{B}^2+\bar{A}^3}\right) ^{1/3}}\right. \nonumber \\{} & {} \left. + \left( \bar{B}+\sqrt{\bar{B}^2+\bar{A}^3}\right) ^{1/3}\right] ~, \end{aligned}$$
(38)

where we defined \(\bar{A}{:}{=}-1-36\frac{\alpha ^2}{\kappa }\xi +36\alpha \xi \epsilon -1296\xi ^2\alpha ^2\epsilon ^2\) and \(\bar{B}=1-\frac{135}{2}\frac{\alpha ^2}{\kappa }\xi -54\alpha \xi \epsilon +1944\frac{\alpha ^3}{\kappa }\xi ^2\epsilon -1944\alpha ^2\xi ^2 \epsilon ^2+46656\alpha ^3\xi ^3\epsilon ^3\). In limit of \(\alpha \rightarrow - \infty \), we get:

$$\begin{aligned} k^2= & {} \frac{2}{3} \left( \root 3 \of {-216 {(\xi \kappa )}^3 \epsilon ^3+\sqrt{-{(\xi \kappa )}^3 \left( 27 \xi \kappa \epsilon ^2+1\right) }-9 {(\xi \kappa )}^2 \epsilon }\right. \nonumber \\{} & {} +\left. \frac{{(\xi \kappa )} \left( 36 {\xi \kappa } \epsilon ^2+1\right) }{\root 3 \of {-216 {(\xi \kappa )}^3 \epsilon ^3+\sqrt{-{(\xi \kappa )}^3 \left( 27 \xi \kappa \epsilon ^2+1\right) }-9 {(\xi \kappa )}^2 \epsilon }}-\right. \nonumber \\{} & {} \left. -6 \xi \kappa \epsilon \right) \end{aligned}$$
(39)

5.2 Pure non-minimal matter-curvature coupling

This case is characterized by setting \(\alpha =0\). In the absence of kinetic corrections, the solution for the Jeans wavenumber is given by:

$$\begin{aligned} k^2_{\pm }=-\beta \xi \pm \sqrt{\gamma \xi +\beta ^2\xi ^2 }~, \end{aligned}$$
(40)

however, only \(k^2_+\) yields a physical solution as it provides real solutions for k provided \(\gamma >0\), as we want to preserve a positive gravitational coupling during a gravitational collapse scenario (Fig. 2).

When kinetic effects are added, the solution becomes:

$$\begin{aligned} k^2_{\pm }=-\beta \xi -6\epsilon \xi \pm \sqrt{\gamma \xi +\left( 6\epsilon \xi +\beta \xi \right) ^2} \end{aligned}$$
(41)

The positive branch, \(k^2_+\), yields positive solutions provided that \(\gamma >0\).

Fig. 2
figure 2

Behavior of the squared modified Jeans’ wavelength with only quantum correction, or with quantum and kinetic effects for a pure non-minimal matter-curvature coupling model

Full size image

5.2.1 \(2\alpha \gamma =\beta \)

This is a special case where a given combination of parameters allows for a simpler analysis. The solutions for the Jeans’ wavenumber are given by:

$$\begin{aligned} k^2=-\frac{1}{9 \alpha }\left[ 1+\varDelta _2+\frac{1}{ \varDelta _2}\right] \end{aligned}$$
(42)

where \(\varDelta _2^3=1-\frac{243 \alpha ^2 \xi \gamma }{2}+\sqrt{-243 \xi \alpha ^2 \gamma +\left( \frac{243}{2}\right) ^2 \xi ^2\alpha ^4\gamma ^2}\).

To ensure real solutions for \(k^2\), a further condition can be found: \(\alpha ^2\gamma =\frac{\alpha \beta }{2}>\frac{4h^2}{243m \rho _0}\). Moreover, \(\varDelta _2\) needs to be negative to give positive solutions for \(k^2\).

When kinetic effects are added, the solution becomes:

$$\begin{aligned} k^2= & {} -\frac{1}{9\alpha }\left[ 1+36\alpha \epsilon \xi -\frac{A_2}{\left( B_2+\sqrt{B_2^2+A_2^3}\right) ^{1/3}}\right. \nonumber \\{} & {} +\left. \left( B_2+\sqrt{B_2^2+A_2^3}\right) ^{1/3}\right] \end{aligned}$$
(43)

where \(A_2{:}{=}108\alpha \epsilon \xi -(1+36\alpha \epsilon \xi )^2\) and \(B_2{:}{=}1-\frac{243}{2}\alpha ^2\gamma \xi -54\alpha \epsilon \xi -1944\alpha ^2\epsilon ^2\xi ^2+46656\alpha ^3\epsilon ^3\xi ^3\). However, as it can be numerically found, there are no real solutions for k, either with quantum corrections alone, or with kinetic effects, hence this scenario does not provide a good physical description of a gravitationally collapsing system.

6 Astrophysical systems

One way to assess the physical viability of the previous solutions is to test with regions in the Universe that can experiment star formation. One of such examples is in Bok globules, which are nearby isolated clouds of interstellar gas and dust with simple shapes. They have characteristic temperatures of the order of 10K, and masses of \(10 M_{\odot }\) which are close to their corresponding Jeans’ masses.

The Jeans’ mass \(M_J\) is defined as the mass inside a sphere of diameter \(2\pi /k_J\), where the Jeans’ wavenumber reads \(k_J^2=\frac{4\pi G \rho _0}{c_s^2}\), being \(c_s\) the sound speed. Thus, these astrophysical systems are perfect candidates to test alternative theories of gravity as well as infer whether kinetic and quantum corrections are observable.

Moreover, for the NMC alternative theories of gravity, we have a modified Jeans’ mass:

$$\begin{aligned} \tilde{M}_J=\left( \frac{k_J^2}{\tilde{k}^2_J}\right) ^{3/2}M_J~, \end{aligned}$$
(44)

where \(\tilde{k}_J\) is the modified Jeans’ mass and corresponds to the solutions we have obtained in the previous section, \(\tilde{k}_J{:}{=}k(\omega =0)\).

Noting that \(\rho _0=n_{H_2}\mu m_p\), being \(n_{H_2}\) the particle number density and \(\mu =2.33\) the mean molecular weight for Bok globules [69], and \(m_p\) the proton mass, we can find that the typical parameters for this system are \(\epsilon =\mathcal {O}(10^{20})\) and \(\xi =\mathcal {O}(10^{-1})\). This means that the kinetic corrections, \(\epsilon \), are several orders of magnitude higher than the quantum ones, \(\xi ^{-1}\), for Bok globules, hence we can neglect the latter ones. Therefore, these considerations together with the phenomenological factor \(\tilde{M}_J=\left( \frac{2}{5}\right) M_J\) found in Ref. [31] as a sufficient condition for matching with the Bok globules stability observations [69], we can find saturating bounds for the parameters of the gravity models parameters, namely \(\alpha , ~\beta , ~\gamma \). That is, for alternative theories of gravity, the following equation, provided the observational saturation bound, has to be solved:

$$\begin{aligned} \left( \frac{k_J^2}{\tilde{k}^2_J}\right) ^{3/2}=\frac{2}{5}~. \end{aligned}$$
(45)

Let us now consider the standard scenario of General Relativity and assess the modifications that appear in the NMC gravity model for the three subclasses analyzed in Sect. 5. To this end, we shall resort to power-law functions for both \(f_1(R)\) and \(f_2(R)\) as it has been shown to correctly address cosmological and astrophysical problems such as dark matter, dark energy, inflation, gravitational waves or the cosmic virial theorem [23,24,25,26,27,28].

6.1 General relativity

GR together with kinetic corrections yields a modified Jeans wavenumber of the form:

$$\begin{aligned} k^2=\frac{1}{12\epsilon \kappa }~, \end{aligned}$$
(46)

and as expected these corrections counteract the gravitational instability, hence we get modified Jeans’ masses, \(M_J^{(T)}\), which are higher than the observed masses, as depicted in Table 1. This results in expecting that kinetic corrections counteract the gravitational collapse, hence providing stability. On the other hand, there are regions which are observed being unstable, hence opening an avenue for allowing alternative theories of gravity, whose additional terms can ensure a suitable gravitational collapse in some scenarios.

We further note that the difference between the computed Jeans mass in General Relativity reported by Ref. [70] and our values lies in the fact that we performed a long wavelength expansion, which resulted in a multiplying factor 3 in the term \(k_B T\) of the dispersion relation. This occurs either in a quantum kinetic regime as in our problem [52], or in a pure classical kinetic approach by performing such expansion. Nevertheless, we shall use the values of Ref. [70] as a direct comparison to ours in the scenarios of General Relativity and f(R) theories.

Table 1 Kinetic temperature, particle number density, mass, Jeans’ mass and observed stability from several Bok globules [69, 70], with the General Relativity predictions when kinetic effects are included
Full size table

6.2 f(R) theories

For f(R) theories, the solution reads:

$$\begin{aligned} k^2_{\pm }=\frac{\alpha /\kappa -3 \epsilon \pm \sqrt{(\alpha /\kappa ) ^2 +3 \alpha /\kappa \epsilon +9 \epsilon ^2}}{18 \alpha \epsilon } ~, \end{aligned}$$
(47)

being \(k^2_+\) the only viable physical solution.

By solving Eq. (45), we get two different solutions for \(\alpha \) for each Bok globule, as shown in Table 2.

Table 2 Kinetic temperature, particle number density, mass, Jeans’ mass and observed stability from several Bok globules [69, 70], and the saturation bounds found for \(\alpha \) in f(R) theories
Full size table

We can explore a particular model of this scenario, namely the Starobinsky model \(f_1(R)=R+aR^2\) [19], for which \(\alpha =2a\), and the results obtained follows from the above Table. In particular, we can note that in the Starobinky model \(2a\approx 3,84\times 10^{13}\) which is remarkably similar to the values found for \(\alpha _2\).

Furthermore, for \(f_1(R)=R+aR^n\), for \(n>1\), we have \(\alpha =0\), and we expect to have the same results as in the General Relativity case. We note that this occurs since we have developed the Jeans analysis for a low field expansion, as it is standard. However, if the background spacetime had more curvature/nonlinear effects, the corresponding, and much more elaborated, form of \(\alpha \) would differ from its GR counterpart.

6.3 Pure NMC

The pure non-minimal matter-curvature coupling gravity model yields the following solution for the Jeans wavenumber:

$$\begin{aligned} k^2=\frac{\gamma }{2\beta + 12 \epsilon }~, \end{aligned}$$
(48)

which in its turn, by solving equation Eq. (45), we get the following relation:

$$\begin{aligned} \gamma = 50^{1/3} (\beta +6 \epsilon ) k_J^2~. \end{aligned}$$
(49)

Let us consider a linear pure non-minimal coupling function, namely \(f_2(R)=1+b R\), which has been recently constrained to have a coupling \(|b| < 2\times 10^{-12}m^2\) [71]. This model has \(\gamma =1/\kappa \) and \(\beta =b/\kappa \), hence we can estimate the modified Jeans’ mass as \(\tilde{M}_J=5.2M_J\). Therefore, we can conclude that this linear functions is not suitable for Bok globules.

We can further consider a general function \(f_2(R)=1+b R^n\), for \(n>1\), for which \(\gamma =1/\kappa \) and \(\beta =0\), which also fails for explaining Bok globules. Therefore, the non-minimal coupling alone, given by a power-law function, is not sufficient to explain the observed (in)stability of Bok globules.

6.4 \(2\alpha \gamma =\beta \)

In this particular case of gravity theories, we get:

$$\begin{aligned} k^2_{\pm }=\frac{-\epsilon \pm \sqrt{\alpha \gamma \epsilon +\epsilon ^2}}{6\alpha \epsilon }~, \end{aligned}$$
(50)

where only the positive sign solution is physically meaningful, thus yielding a relation between \(\alpha \) and \(\gamma \), by solving Eq. (45), which is given by:

$$\begin{aligned} \gamma = \frac{90 \root 3 \of {20} \alpha \epsilon k_J^4+3\root 3 \of {50}\epsilon k_J^2\left( 4+675 \alpha ^3 k_J^2+\left| 4-675 \alpha ^3 k_J^6\right| \right) {}^{\frac{2}{3}}}{\left( 4+675 \alpha ^3 k_J^2+\left| 4-675 \alpha ^3 k_J^6\right| \right) {}^{\frac{1}{3}}}~. \end{aligned}$$
(51)

Let us consider the case of \(f_1(R)=R+aR^2\) and \(f_2(R)=1+bR\), for which a and b are related to each other according to the relation \(2\alpha \gamma =\beta \). For this choice of functions, we have \(\alpha =2a\), \(\beta =b/\kappa \) and \(\gamma =1/\kappa \), hence \(4a=b\). Thus:

$$\begin{aligned} b=\frac{10^{2/3}-60 \root 3 \of {5} \epsilon \kappa k_J^2}{225 \root 3 \of {2} \epsilon \kappa k_J^4}~, \end{aligned}$$
(52)

which for the Bok globules reported in Table 1, we get an estimate of \(b=\mathcal {O}(-10^{29})~m^2\), which is a too strong coupling constant, hence this specific combination of the functions \(f_1\) and \(f_2\) seems to be not viable.

7 Conclusions

In this paper, we have studied the phenomenon of Jeans gravitational instability in the context of non-minimal matter-curvature coupling alternative gravity model, in the quantum regime. In the weak-field limit, the model reduces to a modified Schrödinger–Poisson (or Schrödinger–Newton) model, that we have treated kinetically by relying on the use of Wigner functions, allowing to represent quantum states in the classical phase space. Our results generalize those of [31] and enable to study the interplay between non-minimal matter-curvature effects, finite temperature effects, and quantum effects. We have discussed special cases of the model (general relativity, f(R) theories, pure non-minimal coupling, and \(2\alpha \gamma = \beta \)) and have used the data of Ref. [69] for Bok globules to constraint the parameters of the model.

Bok globules constitute, in fact, an excellent laboratory to test modified gravity models because their mass is of the same order as their Jeans mass; hence a small deviation from the Jeans criterion leads to a different prediction for their stability. Our approach remains nevertheless valid for compact objects where quantum effects are dominant, opening up novel avenues to test this alternative gravity model. Besides, given the universality of the Schrödinger–Poisson model, which applies equally well to other media, such as plasmas and cold atomic clouds, our results may guide future laboratory experiments to emulate alternative theories of gravity in condensed-matter gravity analogs (see for instance [72, 73]).