# Pulsational mode stability in complex EiBI-gravitating polarized astroclouds
with $(r,q)$-distributed electrons

Dipankar Ray and Pralay Kumar Karmakar 11footnotetext: Corresponding author.

###### Abstract

The pulsational mode of gravitational collapse (PMGC) originating from the
combined gravito-electrostatic interaction in complex dust molecular clouds
(DMCs) is a canonical mechanism leading to the onset of astronomical structure
formation dynamics. A generalized semi-analytic model is formulated to explore
the effects of the Eddington-inspired Born-Infeld (EiBI) gravity, non-thermal
$(r,q)$-distributed electrons, and dust-polarization force on the PMGC
stability concurrently. The thermal ions are treated thermo-statistically with
the Maxwellian distribution law, and the non-thermal electrons with the
$(r,q)$-distribution law. The constitutive partially ionized dust grains are
modeled in the fluid fabric. A spherical normal mode analysis yields a
generalized linear PMGC dispersion relation. Its oscillatory and propagation
characteristics are investigated in a judicious numerical platform. It is
found that an increase in the polarization force and positive EiBI parameter
significantly enhances the instability, causing the DMC collapse and vice
versa. The electron non-thermality spectral parameters play as vital
stabilizing factors, and so on. Its reliability and applicability are finally
outlined in light of astronomical predictions previously reported in the
literature.

##  1 Introduction

Interstellar molecular clouds, vast swirling heterogeneous dust and gas
mixtures, serve as the birthplaces for diverse bounded astrophysical structure
evolution, such as stars, planets, and so forth [1, 2, 3, 4]. The fundamental
physics of bounded structure formation has remained as a long-standing obscure
research area yet to be well illuminated. In 1902, Sir James Hopwood Jeans, a
renowned British astronomer and theoretical physicist, has introduced the
Jeans (gravitational) instability criterion to elucidate the gravitational
collapse dynamics in nebular gaseous media. According to his gravitational
collapse theory [2], the stellar structure formation process gets naturally
initiated with the onset of gravitational instability. This instability gets
triggered against the usual hydrostatic equilibrium configuration
conjugationally established by the inward gravitational pressure force and the
outward thermal pressure force. When the mass or length of a molecular cloud
exceeds the critical Jeans mass-length scale limits, the gravitational
pressure force dominates the thermal pressure force. It results in the
collapse of giant molecular clouds fragmented into local cloudlets,
consequently condensed in the form of stellarsimals, planetesimals,
protostructures, and so forth [5, 6, 7, 8].

In self-gravitating complex (dusty) plasmas composed of electrons, ions,
neutral, and charged (positive or negative) dust, the presence of charged dust
grains generates a repulsive electrostatic force that opposes the inward
gravitational force. Additionally, because of the greater mass of the dust
particles against the lighter plasma compositions (electrons and ions), the
self-gravitational effects of the dust species become significant in extensive
astrophysical systems, like dust molecular clouds (DMCs). The DMCs are the
dense sites of the interstellar medium (ISM), consisting of gas and micron to
sub-micron-sized dust (silicate, graphite, magnetite, etc.) [1, 9]. The dust
grains present in the partially ionized DMCs modify the traditional Jeans
instability. Thus, in order to adapt the Jeans instability for structure
formation in complex self-gravitating DMCs, a Jeans-like hybrid mode known as
the pulsational mode of gravitational collapse (PMGC) has been introduced in
1999 by Dwivedi et al. [10]. This hybrid instability emerges in a partially
ionized DMC due to the interaction between gravitational and electrostatic
forces, termed as gravito-electrostatic coupling [1, 10].

Since the initial work on the PMGC instability, several researchers have shown
a keen interest in this mode of gravitational collapse. Consequently, various
scenarios in both linear and non-linear regimes have been explored. Extensive
analyses of the PMGC dynamics in the presence of dust charge fluctuations have
been conducted, demonstrating that dust charge fluctuations play a crucial
role in the PMGC dynamics excitable in the DMCs [1, 11, 12]. Subsequently, the
area of PMGC dynamics have been extended to the non-linear regime. Various
non-linear model formalisms have accordingly been developed in various
conditions [13, 11]. Apart from the above, the effects of various parameters,
such as viscoelasticity [14], turbulence [12], ion drag [15, 16], and non-
thermality [17, 18, 19] on the PMGC dynamics have been investigated. Recently,
this field has further been expanded to modified gravity frameworks, such as
the EiBI gravity [20], zeroth-order gradient effects [21], and so forth.

It is seen in the literature that the impact of polarization force on the PMGC
dynamics, which is a critical aspect of inhomogeneous plasmas, has yet to be
extensively addressed. In non-uniform plasmas with uneven distributions of
electrons and ions, characterized by the absence of dust charge fluctuations,
the deformation of the plasma Debye sheaths gives rise to a new force known as
the dust polarization force ($\vec{F}_{p}$). It is in addition to the
conventional Coulomb electrostatic force ($\vec{F}_{e}=q_{d}\vec{E}$; where,
$q_{d}$ and $\vec{E}$ stand for the dust charge and electric field,
respectively) [22]. This sheath-polarization force acting on the charged dust
particles vanishes entirely in uniform plasma environments. In complex plasmas
comprising of negatively charged dust grains, the Debye sheath is formed by
ions around the dust grains. The distortion of such plasma sheaths results in
the development of polarization force [23].

The expression for the polarization force is given as
$\vec{F}_{p}=-q_{d}^{2}\vec{\nabla}\lambda_{D}/(2\lambda_{D}^{2})$; where,
$\lambda_{D}=\lambda_{Di}/\sqrt{1+(\lambda_{Di}/\lambda_{De})^{2}}$ denotes
the effective Debye length [24, 25]. Here,
$\lambda_{De(i)}=\sqrt{T_{e(i)}/4\pi e^{2}n_{e(i)}}$ represents the electron
(ion) Debye length, with $T_{e(i)}$ being the electron (ion) temperature, $e$
is the elementary charge (electronic), and $n_{e(i)}$ corresponds to the
number density of electrons (ions). By utilizing the approximation
$T_{i}n_{e}\ll T_{e}n_{i}$, which is suitable for negatively charged dust, and
$T_{i}\vec{\nabla}n_{i}=-n_{i}e\vec{\nabla}\phi$ [25], a simplified form for
the polarization force can be derived as
$\vec{F}_{p}=-q_{d}R(n_{i}/n_{i0})^{1/2}\vec{\nabla}\phi$ [26] (refer to
Appendix for details). Here, $R=q_{d}e/(4T_{i}\lambda_{Di0})$ serves as the
polarization interaction parameter that determines the influence of the
polarization force, with $n_{i0}$ and $\phi$ representing the equilibrium ion
number density and electrostatic potential, respectively. It is evident from
here that a higher dust charge (or lower ion temperature) results in a
stronger polarization force and vice versa.

Numerous studies in the literature have relied on the Maxwellian velocity
distribution of the lighter constitutive species to elucidate plasma dynamics
[10, 27, 16, 20]. However, advancements in observational capabilities have
unveiled that the particle distribution of space plasmas significantly
deviates from the Maxwellian velocity distribution. In 1983, flat-top
characteristics in the electron velocity distributions at the boundary of the
Earth’s bow shock have been reported [28]. Multiple spacecraft missions in the
downstream region of the bow shock have subsequently observed this behaviour
of low-energy electrons. The ARTEMIS spacecraft has provided further evidence
for this phenomenon in the geomagnetic tail region [29]. Additionally, similar
flat-top distribution features near the magnetic reconnection regions have
been reported by analyzing data collected from the four Cluster satellites
[30]. Observational studies of such space plasmas have established the
widespread presence of particle populations characterized by high or super-
thermal energy tails that do not adhere to the standard Maxwellian velocity
distribution. These non-Maxwellian (non-equilibrium) distributions exhibit an
approximate power-law dependence in velocity space, implying a substantial
fraction of particles with energies significantly greater than the average
thermal energy value. Examples in this context are the interstellar medium
[31], planetary nebulae [32], planetary magnetospheres [33, 34], solar wind
[35], and so forth.

To model the non-Maxwellian behaviours, various non-thermal velocity
distribution laws have been proposed, such as the kappa distribution [34],
$q$-nonextensive distribution [36], Cairns distribution [37], and so on. In
this context, the kappa distribution law has been demonstrated to be effective
in characterizing the high-energy particle tails. Likewise, the kappa
distribution is also characterized by a single non-thermal spectral index
$(\kappa)$, which typically controls the high-energy tails. To explain both
the flat-top and high-energy tail features simultaneously, Qureshi et al. [38]
have introduced a new bi-spectral non-Maxwellian electron velocity
distribution function called the generalized $(r,q)$-distribution function in
2004. The $(r,q)$-distribution law specifies two non-thermal spectral indices
($r$ and $q$) that regulate the electron population in different parts of the
energy spectrum. Specifically, $r$ controls the electron population at the
flat-tops, while $q$ governs the high-energy tails [39, 38]. The
$(r,q)$-distribution is a generalized form of the kappa distribution, which
can be reduced to the kappa distribution (for $r=0$ and $q\to\kappa+1$) as
well as the Maxwellian distribution (for $r=0$ and $q\to\infty$). Therefore,
the $(r,q)$-distribution thermo-statistically gives more flexibility in
stability analyses as compared with other non-thermal velocity distribution
laws.

Einstein’s general theory of relativity (GTR) has been proven to be one of the
most successful theories explaining gravity. After the introduction of GTR in
1915, Eddington was the very first person to experimentally verify GTR in 1919
by measuring the bending of light during a solar eclipse [40]. Since then, GTR
has passed numerous tests (in weak-field regime) in diverse astrophysical
environments, such as the solar system [41, 42], with binary pulsars [43],
black hole binaries [44, 45], and in proximity to the Galactic Centre [46].
Despite the remarkable accomplishments of GTR, it falls short in tackling
certain inherent challenges. For instance, it fails to explain the rapid
cosmic expansion of the universe, the formation of singularities in the space-
time fabric during the gravitational collapse of stars and in the early
universe [47, 48], as well as the enigma surrounding dark energy and dark
matter [49]. To overcome these critical issues, many extensions of GTR have
been proposed, such as Scalar-Tensor-Vector gravity (STVG) [50], 4D Einstein-
Gauss-bonnet gravity (4DEGB) [51, 52], $f(R)$ gravity [53], Eddington inspired
Born-Infeld (EiBI) gravity [54], etc. After the seminal work on the Jeans
instability of self-gravitating objects in non-relativistic Newtonian gravity
[2], it has been further extended to GTR [55]. Recently, the impact of
different alternative gravity formalisms on the Jeans collapse dynamics has
been thoroughly investigated [56, 57, 58, 59, 60].

Among the various extensions, the EiBI gravity has been successful in
demonstrating the elimination of such singularities in both the primordial
universe and subsequent gravitational collapse scenarios [54, 61]. The
development of EiBI gravity has been motivated by the Born–Infeld action for
nonlinear electrodynamics and it has been formulated using Eddington’s theory
of gravity [56]. Introduced by Bañados and Ferreira in 2010, the EiBI gravity
theory has been scrutinized in several astrophysical scenarios, such as the
solar system [62], compact objects [47, 63], and even in cosmological scales
[64, 65]; hence, it can be treated as a reliable extension of GTR. The
extensions of GTR are generally obtained through modifications to the
Einstein-Hilbert action expressed as (in the unit of $c$)
$S_{\text{EH}}=1/(16\pi G)\int d^{4}x\sqrt{-g}(R-2\Lambda)$ [48, 47, 66].
Here, $g$ represents the determinant of the metric $g_{\mu\nu}$, $R$ is the
Ricci scalar defined as $R=g^{\mu\nu}R_{\mu\nu}$, $G$ is the universal
gravitational constant, and $\Lambda$ denotes the cosmological constant. The
famous Einstein’s field equations can be obtained by varying the Einstein-
Hilbert action with respect to the metric $g_{\mu\nu}$ as
$R_{\mu\nu}-(1/2)g_{\mu\nu}R+\Lambda g_{\mu\nu}=8\pi GT_{\mu\nu}$ [54, 66,
67]; where, $R_{\mu\nu}$ is the Ricci Tensor and
$T^{\mu\nu}=(g)^{-1/2}\partial S_{M}/\partial g_{\mu\nu}$ denotes the stress-
energy tensor with $S_{M}$ being the matter action.

It is noteworthy that following the remarkable works by Deser et al. [68] and
Vollick [69, 70], Bañados and Ferreira have considered a Born-Infled type
action of the form $S_{\text{EiBI}}=1/(8\pi\chi G)\int
d^{4}x\left(\sqrt{|g_{\mu\nu}+\chi
R_{\mu\nu}}|-\lambda\sqrt{-g}\right)+S_{M}[g,\psi_{M}]$ [54, 48]. Here, $\chi$
serves as the independent parameter for EiBI gravity, $\lambda(\neq 0)$ is a
dimensionless constant, and $\psi_{M}$ represents the matter field. It is
noteworthy that in the limit where $\chi R\ll 1$ and
$\lambda=(\Lambda\chi+1)$, the EiBI action reduces to the conventional
Einstein-Hilbert action [47]. The field equations of EiBI gravity can be
obtained from the EiBI action $(S_{\text{EiBI}})$ using the Palatini formalism
by introducing an auxiliary metric $q_{\mu\nu}=g_{\mu\nu}+\chi R$, leading to
$R_{\mu\nu}\approx\Lambda g_{\mu\nu}+8\pi
G\left[T_{\mu\nu}-(1/2)g_{\mu\nu}T\right]+8\pi
G\chi\left[S_{\mu\nu}-(1/4)Sg_{\mu\nu}\right]+\mathcal{O}(\chi^{2})$ [71, 47].
Here, $S_{\mu\nu}=T_{\mu\nu}-(1/2)TT_{\mu\nu}$ is the source tensor, while $S$
and $T$ represent the traces of $S_{\mu\nu}$ and $T_{\mu\nu}$, respectively.
The third term on the right-hand side of the EiBI-modified field equations
introduces corrections to the standard Einstein’s field equations. As $\chi$
approaches zero, the familiar Einstein’s field equations are recovered.

It is clear from the above scenarios that the EiBI gravity, non-thermal
$(r,q)$-distributed electrons, and dust-polarization force play significant
roles in the collapse of a molecular cloud. Motivated by this, a generalized
semi-analytic model is proposed in this paper incorporating the aforementioned
factors to study the PMGC stability dynamics of a complex DMC.

After the introduction given in section 1, the rest of the paper is organized
as follows. Section 2, outlines the physical model formalism and the basic
governing equations describing the system. Section 3, explains the
mathematical method for linearization of the governing equations followed by
the spherical Fourier analysis. A generalized linear dispersion relation
(quartic in degree) is obtained using the method of decoupling (elimination).
In section 4, the numerical outcomes, along with figures, are discussed.
Finally, the results are concluded with a brief summary and future scope of
this investigation in section 5.

##  2 Physical model and mathematical formalism

To study the linear PMGC dynamics, a weakly coupled $(\Gamma_{\text{Cou}}\ll
1)$, inhomogeneous, and globally quasi-neutral self-gravitating DMC is
considered, ignoring the effects of rotation, magnetic field, and tidal action
resulting from the gravitational effects of distant astrophysical objects. It
consists of four species, namely, $(r,q)$-distributed (non-thermal) electrons,
Maxwellian (thermal) ions, neutral and negatively charged dust particles of
identical size $(r_{d})$ and mass $(m_{d})$. The Coulomb coupling parameter
$(\Gamma_{\text{Cou}})$ associated with the system refers to the ratio of dust
potential energy to dust thermal energy and is mathematically expressed as
$\Gamma_{\text{Cou}}=(1/4\pi\epsilon_{0})\left[q_{d}^{2}/(a_{d}T_{d})\right]\left(\text{in
cgs, }\Gamma_{\text{Cou}}=q_{d}^{2}/(a_{d}T_{d})\right)$; where,
$a_{d}=\left[3/(4\pi n_{d})\right]^{1/3}$ denotes the Wigner-Seitz radius of
the dust grains and $T_{d}$ represents the dust temperature (in eV). A system
is said to be weakly coupled for $\Gamma_{\text{Cou}}\ll 1$ and strongly
coupled for $\Gamma_{\text{Cou}}\gg 1$ [72]. The considered weakly coupled
DMCs are indeed partially ionized in nature with an estimated degree of
ionization $\sim 10^{-7}$ [73]. Owing to their low ionization level, both
neutral and charged dust particles coexist with electrons and ions in such
DMCs. A further simplification is considered by ignoring the dust charge
fluctuations which is justifiable when the hydrodynamic time scale
$(\tau_{d})$ significantly exceeds the grain charging time scale $(\tau_{c})$,
i.e., $\tau_{d}\gg\tau_{c}$ [74]. Here, $\tau_{d}=\omega_{pd}^{-1}$,
$\omega_{pd}$ being the dust plasma frequency and $\tau_{c}=\nu_{d}^{-1}$,
$\nu_{d}$ denotes the grain charging frequency. In a typical DMC, the
hydrodynamic time scale is approximately $\tau_{d}\sim 10^{-2}$ seconds, while
the dust charging time scale is $\tau_{c}\sim 10^{-8}$ seconds [75].
Therefore, dust charge fluctuations can be safely disregarded for this
particular DMC.

Figure 1: A schematic showing the fragmentation processes of DMCs leading to
stellar structure formation triggered by the PMGC instability.

The DMC fragmentation process due to PMGC instability can be illustrated with
a schematic diagram as in figure 1. Here, the complex DMC (figure 1a) is
subjected to various external perturbations, such as explosive shock waves
propagating from nearby supernova explosions, cloud-cloud collisional
interactions, etc. These perturbations may produce some highly dense regions
compared to the overall density of the DMC in the form of diverse parametric
inhomogeneities (gradient effects). Consequently, these high-density regions
start to contract due to the inward self-gravitational pull among the
constitutive heavier dust particles, leading to a rise in outward thermal
pressure. However, due to the very low initial density of the cloud, the
thermal pressure can escape from the parent cloud in the form of various
deceptive processes, such as thermal conduction, radiation, etc. [76, 3].
Hence, the outward thermal pressure does not balance the inward gravitational
pressure at the initial stage of collapse. In contrast, the effective
electrostatic repulsion (electrostatic force between the identical charged
particles and polarization force) tries to counter the gravitational collapse
of the cloud. Depending upon the critical Jeans mass-length limit, a gravito-
electrostatic equilibrium is established, forming a stable cloud region
(figure 1c). However, if the mass of the cloud exceeds the Jeans critical mass
limit, the self-gravitational force prevails over the electrostatic repulsion
force. It results in a growing instability in the cloud (figure 1f), leading
to the fragmentation of the original cloud into (smaller) cloudlets (figure
1b). These equilibrium and collapse processes continue till the cloud is
fragmented to a smaller stellar mass scale, resulting in a protostar formation
(figure 1e). As the cloud density and temperature increase with each
condensation stage, they eventually become opaque to the thermal radiation,
and the radiation can no longer escape from the stellar body, eventually
establishing a hydrostatic equilibrium [76].

The electrons are assumed to follow the generalized $(r,q)$ thermo-statistical
distribution law given in generic notations [38, 77] as

$F_{r,q}(v)=\frac{\alpha}{\pi
v_{\text{th}}^{3}}\left(1+\frac{1}{q-1}\left[\frac{v^{2}-2e\phi/m_{e}}{\beta\left(2T_{e}/m_{e}\right)}\right]^{(r+1)}\right)^{-q};$
(2.1)
$\alpha=\frac{3\Gamma[q](q-1)^{-\frac{3}{(2+2r)}}}{4\beta^{\frac{3}{2}}\Gamma\left[q-\frac{3}{2+2r}\right]\Gamma\left[1+\frac{3}{2+2r}\right]},$
(2.2)
$\beta=\frac{3(q-1)^{-\frac{1}{(1+r)}}\Gamma\left[q-\frac{3}{2+2r}\right]\Gamma\left[\frac{3}{2+2r}\right]}{2\Gamma\left[q-\frac{5}{2+2r}\right]\Gamma\left[\frac{5}{2+2r}\right]}.$
(2.3)

The number density of non-thermal $(r,q)$-distributed electrons $(n_{e})$ can
be obtained by integrating (2.1) over a velocity space as [77]

$n_{e}=n_{e0}\left[1+A\frac{e\phi}{T_{e}}+B\left(\frac{e\phi}{T_{e}}\right)^{2}\right];$
(2.4)

$A=\frac{(q-1)^{-\frac{1}{(1+r)}}\Gamma\left[q-\frac{1}{2+2r}\right]\Gamma\left[\frac{1}{2+2r}\right]}{2\beta\Gamma\left[q-\frac{3}{2+2r}\right]\Gamma\left[\frac{3}{2+2r}\right]},$
(2.5)

$B=-\frac{(1+4r)(q-1)^{\frac{-2}{(1+r)}}\Gamma\left[q+\frac{1}{(2+2r)}\right]\Gamma\left[-\frac{1}{2+2r}\right]}{8\beta^{2}\Gamma\left[q-\frac{3}{2+2r}\right]\Gamma\left[\frac{3}{2+2r}\right]}.$
(2.6)

Due to the heavier mass of the ions than that of the electrons, the inertia of
the ions would be larger than that of the electrons. Thus, the ions are
considered to follow the Maxwellian velocity distribution on a slow dust
inertial time scale, with its number density given as [27]

$n_{i}=n_{i0}\exp\left(-\frac{e\phi}{T_{i}}\right).$ (2.7)

In equations (2.1)-(2.7), $\Gamma[z]\left(\text{defined simply as
}\Gamma[z]=\int_{0}^{\infty}e^{-x}x^{z-1}dx\right)$ denotes the so-called
characteristic gamma function [38]. The symbols $r$ and $q$ represent the non-
thermal spectral indices specifying the electronic thermo-statistics. The
variable $\phi$ represents the electrostatic potential,
$v_{\text{th}}=(2T_{e}/m_{e})^{1/2}$ denotes the thermal velocity of the
electrons, $T_{e(i)}$ corresponds to the electron (ion) temperature (in eV),
$m_{e}$ is the mass of an electron, and $n_{e(i)}$ is the number density of
electrons (ions), with $n_{e0(i0)}$ being their corresponding equilibrium
number densities. It is crucial to emphasize that the spectral indices $r$ and
$q$ in the $(r,q)$-distribution law must follow the conditions $q>1$ and
$q(r+1)>5/2$ [38, 77]. In the special case when $r=0$ and $q\to\infty$, the
$(r,q)$-distribution transitions to the Maxwellian distribution; where, the
multi-order non-thermality coefficients A and B become $A=1$ and $B=1/2$.

The dust particles, when emerged in a plasma environment, become electrically
charged (through contact electrification process) due to the random collisions
of lighter constituent species with the dust grains. As the electrons possess
a greater thermal velocity compared to the ions, they reach the dust grain
surface significantly faster than the ions. This rapid accumulation of the
electrons on the grain surfaces makes them negatively charged [72, 78].
However, due to the partially ionized nature of the cloud, not all dust grains
are electrically charged; hence, both neutral and negative dust grains are
present in the system. Thus, the neutral and negatively charged grains are
considered as two different fluids. Furthermore, in the context of studying
the collapse of a dense DMC, it is reasonable to consider the exchange of
momentum between charged and neutral dust grains. The rates of binary
collisions for momentum transfer from the charged to the neutral dust
particles and vice versa are given by $\nu_{cn}\approx\pi
r_{d}^{2}n_{dn0}v_{td}$ and $\nu_{nc}\approx\pi r_{d}^{2}n_{dc0}v_{td}$,
respectively. Here, $r_{d}$ is the radius of the dust grain,
$n_{dco}(n_{dn0})$ is the equilibrium number density of charged (neutral)
dust, and $v_{td}=(2T_{d}/m_{d})^{1/2}$ is their common thermal velocity. The
frictional force term can be neglected for neutral dust grains under the
condition that the Jeans mode frequency $\omega\gg\nu_{cn}$ with
$\nu_{cn}/\nu_{nc}=n_{dn0}/n_{dc0}\gg 1$ [10].

The neutral dust dynamics can be described using the continuity and momentum
equations with the negligible frictional effect respectively as

$\partial_{t}n_{dn}+\vec{\nabla}.(n_{dn}\vec{v}_{dn})=0,$ (2.8)
$\partial_{t}v_{dn}+(\vec{v}_{dn}.\vec{\nabla})v_{dn}=-\vec{\nabla}\psi.$
(2.9)

Again, the charged dust dynamics can be described using continuity and
momentum equations with significant frictional effect respectively as

$\partial_{t}n_{dc}+\vec{\nabla}.(n_{dc}\vec{v}_{dc})=0,$ (2.10)
$\partial_{t}\vec{v}_{dc}+(\vec{v}_{dc}.\vec{\nabla})\vec{v}_{dc}=-\frac{q_{d}}{m_{d}}\vec{\nabla}\phi+\frac{q_{d}}{m_{d}}R\left(\frac{n_{i}}{n_{i0}}\right)^{1/2}\vec{\nabla}\phi-\vec{\nabla}\psi-\nu_{cn}(\vec{v}_{dc}-\vec{v}_{dn}).$
(2.11)

In the above equations (2.8)-(2.11), $n_{dn(dc)}$ and $v_{dn(dc)}$ stand for
the neutral (charged) dust number density and neutral (charged) dust fluid
velocity, respectively. Here, $q_{d}=-Z_{d}e$ denotes the dust charge with
$Z_{d}=|q_{d}/e|$ being the dust charge number and $e$ is the elementary
charge (electronic). Moreover, $R$ represents the polarization interaction
parameter, and $\nu_{cn}$ is the binary collisional rate of momentum transfer
from the charged to the neutral dust grains. Furthermore, $\phi$ and $\psi$
respectively denote the electrostatic and gravitational potential.

The modified Poisson equation serves as the starting point for
phenomenological studies using EiBI gravity. The required Poisson equation can
be obtained from the modified Einstein’s field equation (mentioned in section
1) by taking the weak field limit as [47, 65]

$\nabla^{2}\psi=4\pi
G\rho_{mj}+\left(\frac{\chi}{4}\right)\nabla^{2}\rho_{mj},$ (2.12)

where, $\psi$ represents the gravitational potential, $\rho_{mj}=n_{j}m_{j}\
(j=e,i,dc,dn)$ denotes the matter density, and $\chi$ serves as the EiBI
parameter. Numerous constraints have been imposed on the EiBI parameter for
different astrophysical objects, such as white dwarfs $(\chi\leq 4.86\times
10^{2}$ kg^-1.m^5.s^-2) [47], neutron stars $(\chi<10^{-2}$ kg^-1.m^5.s^-2)
[61, 63], Sun $(\chi\leq 3\times 10^{5}$ kg^-1.m^5.s^-2)[62], and so forth.
The tightest nuclear constraint on the EiBI parameter $(\chi)$ has been
imposed by Avelino in 2012 [65], given as $\chi\leq 10^{-3}$ kg^-1.m^5.s^-2).
It is important to note that the usual gravitational Poisson equation can be
recovered from (2.12) for $\chi=0$.

Finally, the DMC model is closed by utilizing the electro-gravitational
Poisson equations, which provide the potential distributions generated by the
respective density fields, cast as

$\nabla^{2}\phi=4\pi e(n_{e}-n_{i}-\frac{q_{d}n_{dc}}{e}),$ (2.13)
$\nabla^{2}\psi=4\pi
Gm_{d}(n_{dc}+n_{dn}-n_{d0})+\left(\frac{\chi}{4}\right)\nabla^{2}[m_{d}(n_{dc}+n_{dn}-n_{d0})].$
(2.14)

Here, the first term on the right side of (2.14) signifies the effect of
Newtonian gravity, whereas the second term accounts for the modification made
by the EiBI gravity, and $n_{d0}(=n_{dc0}+n_{dn0})$ is the net equilibrium
dust number density, serving as the Jeans swindle.

As the astrophysical self-gravitating plasmas are usually confined in
spherical geometric volumes, a spherically symmetric approximation (radial 1-D
case) is employed rather than a full spherical geometry (spherical 3-D case).
This simplification is justifiable under the assumption of radial symmetry,
which holds when the fluctuation wavelength significantly exceeds the grain-
to-grain distance [27]. The fundamental governing equations (2.8)-(2.14)
defining the DMC dynamics can be written in the spherically symmetric
coordinate $(\rho,t)$ with all usual notations as

$\partial_{t}n_{dn}+\rho^{-2}\partial_{\rho}\left(\rho^{2}n_{dn}v_{dn}\right)=0,$
(2.15)
$\partial_{t}v_{dn}+(v_{dn}\partial_{\rho})v_{dn}=-\partial_{\rho}\psi,$
(2.16)
$\partial_{t}n_{dc}+\rho^{-2}\partial_{\rho}\left(\rho^{2}n_{dc}v_{dc}\right)=0,$
(2.17)
$\partial_{t}v_{dc}+\left(v_{dc}\partial_{\rho}\right)v_{dc}=-\frac{q_{d}}{m_{d}}\partial_{\rho}\phi+\frac{q_{d}}{m_{d}}R\left(\frac{n_{i}}{n_{i0}}\right)^{\frac{1}{2}}\partial_{\rho}\phi-\partial_{\rho}\psi-\nu_{cn}\left(v_{dc}-v_{dn}\right),$
(2.18) $\rho^{-2}\partial_{\rho}\left(r^{2}\partial_{\rho}\phi\right)=4\pi
e\left(n_{e}-n_{i}-\frac{q_{d}n_{dc}}{e}\right),$ (2.19)
$\displaystyle\rho^{-2}\partial_{\rho}\left(\rho^{2}\partial_{\rho}\psi\right)$
$\displaystyle=4\pi Gm_{d}\left(n_{dc}+n_{dn}-n_{d0}\right)$ (2.20)
$\displaystyle\quad+\left(\frac{\chi}{4}\right)\rho^{-2}\partial_{\rho}\left(\rho^{2}\partial_{\rho}[m_{d}(n_{dc}+n_{dn}-n_{d0})]\right).$

Equations (2.15)-(2.20) are the governing equations required to describe the
dynamics of the spherical DMC model. Here,
$\partial_{\rho}\equiv\partial/\partial\rho$ and
$\partial_{t}\equiv\partial/\partial t$ signify the space and time gradient
operators, respectively. It is noteworthy here that the geometrical curvature
effects are responsible for the formation of $1/\rho$ terms, and the planar
equations can be retraced under the geometric approximation $\rho\to\infty$,
$\rho$ being the radial distance.

##  3 Linearization and Fourier analysis

Employing a standard approach of spherical wave analysis [79, 80, 60], a
linear perturbation $(f_{1})$ is introduced to the relevant physical fluid
parameters around their respective homogeneous equilibrium values $(f_{0})$.
The multi-parametric perturbations with amplitude factor $f_{10}$ assume the
mathematical shape of a symmetric spherical wave in a usual symbolism without
any dimensional disparity in the absence of polar and azimuthal contributions
as

$f(\rho,t)=(f_{0}+f_{1})=f_{0}+\rho^{-1}f_{10}\exp\left[-i(\omega
t-k\rho)\right];$ (3.1)
$f=[n_{dn}\;n_{dc}\;n_{e}\;n_{i}\;v_{dn}\;v_{dc}\;\phi\;\psi]^{T},$ (3.2)
$f_{0}=[n_{dn0}\;n_{dc0}\;n_{e0}\;n_{i0}\;0\;0\;0\;0]^{T},$ (3.3)
$f_{1}=[n_{dn1}\;n_{dc1}\;n_{e1}\;n_{i1}\;v_{dn1}\;v_{dc1}\;\phi_{1}\;\psi_{1}]^{T}.$
(3.4)

As already mentioned above, it is repeated broadly that the symbol, $f_{1}$,
represents a small (linear) perturbation of $f$ about $f_{0}$ such that
$f_{1}\approx\rho^{-1}f_{10}\exp[-i(\omega t-k\rho)]$. This sinusoidal
spherical wave has angular frequency $\omega$ and angular wavenumber $k$. It
is noteworthy that $f_{0}$ and $\rho^{-1}f_{10}$ have the same dimensions. $T$
represents the transpose operation thereof. The exponential term embodies the
inherent wave-like nature of the perturbations. The inherent nature of the
adopted linearization approach dictates the exclusion of all higher-order
terms, as their contributions become insignificant in the current model
analysis. The linear differential operators accordingly transform as
$\partial_{\rho}\rightarrow(ik-\rho^{-1})$,
$\partial_{t}\rightarrow(-i\omega)$, and
$\partial_{\rho}^{2}\rightarrow[(-k^{2}+2\rho^{-2})-i(2k\rho^{-1})]$ [80]. The
Fourier analysis is used to transform the fluctuation dynamics from the real
coordination (direct) space $(\rho,t)$ to the wave (reciprocal) space
$(k,\omega)$ in a physically judicious manner. Accordingly, the applied
spherical Fourier analysis reduces the equations (2.1),(2.2), and (2.9)-(2.14)
in their linearized forms given respectively as

$n_{e1}=n_{e0}\left(\frac{Ae\phi_{1}}{Te}\right),$ (3.5)
$n_{i1}=-n_{i0}\left(\frac{e\phi_{1}}{T_{i}}\right),$ (3.6)
$n_{dn1}=\frac{n_{dn0}(k^{2}+\rho^{-2})}{\omega^{2}}\psi_{1},$ (3.7)
$v_{dn1}=\left(\frac{ik-\rho^{-1}}{i\omega}\right)\psi_{1},$ (3.8)
$n_{dc1}=\frac{n_{dc0}(k^{2}+\rho^{-2})}{\omega^{2}}\left[\psi_{1}+\left(\frac{q_{d}}{m_{d}}\right)(1-R)\left(1+\frac{i\nu_{cn}}{\omega}\right)^{-1}\phi_{1}\right],$
(3.9)
$v_{dc1}=\frac{(q_{d}/m_{d})(1-R)(ik-\rho^{-1})\phi_{1}+(ik+\rho^{-1})\psi_{1}-\nu_{cn}(v_{dn1})}{i\omega-\nu_{cn}},$
(3.10) $\phi_{1}=\frac{4\pi
e}{k^{2}}\left[n_{i1}-n_{e1}+\left(\frac{q_{d}}{e}\right)n_{dc1}\right],$
(3.11) $\psi_{1}=-\left(\frac{4\pi
e}{k^{2}}-\frac{\chi}{4}\right)m_{d}(n_{dc1}+n_{dn1}).$ (3.12)

Substituting the value of $n_{dn1}$ from (3.7) in (3.12), the perturbed
gravitational potential $(\psi_{1})$ can be obtained as follows

$\psi_{1}=-\left(\frac{4\pi
G}{k^{2}}-\frac{\chi}{4}\right)m_{d}n_{dc1}\left[\frac{4k^{2}\omega^{2}}{4k^{2}\omega^{2}+(k^{2}+r^{-2})\left(4\omega_{J}^{2}-4\omega_{Jc}^{2}-\chi
m_{d}n_{dn0}k^{2}\right)}\right].$ (3.13)

Here, $\omega_{J}=(\omega_{Jc}+\omega_{Jn})=\sqrt{4\pi Gm_{d}n_{d0}}$ and
$\omega_{Jc}=\sqrt{4\pi Gm_{d}n_{dc0}}$ are the effective critical Jeans
frequency for the entire dust species and the critical Jeans frequency for the
charged dust particles [2, 18]. Utilizing the value of perturbed gravitational
potential $(\psi_{1})$ from (3.13) in (3.9), the perturbed charged dust number
density $(n_{dc1}$) can be rewritten as

$\displaystyle n_{dc1}$
$\displaystyle=\frac{n_{dc0}(k^{2}+r^{-2})}{\omega^{2}}\left[\left(\frac{q_{d}}{m_{d}}\right)(1-R)\left(1+\frac{i\nu_{cn}}{\omega}\right)^{-1}\right]$
(3.14)
$\displaystyle\quad\times\left[\frac{4k^{2}\omega^{2}+(k^{2}+r^{-2})\left(4\omega_{J}^{2}-4\omega_{Jc}^{2}-\chi
m_{d}n_{dn0}k^{2}\right)}{4k^{2}\omega^{2}+(k^{2}+r^{-2})\left(4\omega_{J}^{2}-\chi
m_{d}n_{d0}k^{2}\right)}\right]\phi_{1}.$

Finally, substituting the values of $n_{e1}$, $n_{i1}$, and $n_{dc1}$ from
(3.5),(3.6), and (3.14) in (3.11), an unnormalized quartic dispersion relation
is achieved as follows

$\displaystyle\left(\frac{\omega}{\omega_{J}}\right)^{4}+i\frac{\nu_{cn}}{\omega_{J}}\left(\frac{\omega}{\omega_{J}}\right)^{3}+\left[\left(1-\frac{\chi
m_{d}n_{d0}}{4\omega_{J}^{2}}k^{2}-\frac{\omega_{pdc}^{2}(1-R)}{\omega_{J}^{2}(1+k^{-2}\lambda_{D_{r,q}}^{-2})}\right)(k^{2}+r^{-2})k^{-2}\right]\left(\frac{\omega}{\omega_{J}}\right)^{2}$
(3.15) $\displaystyle+i\frac{\nu_{cn}}{\omega_{J}}\left[\left(1-\frac{\chi
m_{d}n_{d0}}{4\omega_{J}^{2}}k^{2}\right)(k^{2}+r^{-2})k^{-2}\right]\left(\frac{\omega}{\omega_{J}}\right)-(1-R)\lambda_{D_{r,q}}^{2}\left(\frac{\omega_{pdc}^{2}}{\omega_{J}^{2}}\right)$
$\displaystyle\times\left(\frac{4\omega_{J}^{2}-4\omega_{Jc}^{2}-\chi
m_{d}n_{dn0}k^{2}}{4\omega_{J}^{2}}\right)\left(\frac{k^{2}+r^{-2}}{1+k^{2}\lambda_{D_{r,q}}^{2}}\right)k^{-2}=0,$

where, $\omega_{pd}=\sqrt{4\pi q_{d}^{2}n_{d0}/m_{d}}$ and
$\omega_{pdc}=\sqrt{4\pi q_{d}^{2}n_{dc0}/m_{dc}}$ are the effective dust-
plasma oscillation frequency (contributed by both the neutral and charged dust
grains) and the charged dust-plasma oscillation frequency (contributed by the
charged dust grains only). The term, $\lambda_{D_{r,q}}$, is the effective
dust-plasma Debye length for $(r,q)$-velocity distribution, cast as

$\lambda_{D_{r,q}}=\left[\frac{T_{e}T_{i}}{4\pi
e^{2}(An_{e0}T_{i}+n_{i0}T_{e})}\right]^{\frac{1}{2}},$ (3.16)

which becomes the standard effective plasma Debye length $(\lambda_{D})$ for
$A=1(r=0,q\to\infty)$.

To check the reliability and validation of (3.15) in comparison with previous
predictions found in the literature, all the newly added modifications are
turned off, i.e., $R=0,\ \chi=0,\ \rho\to\infty$, and $A=1$. Now, using the
approximation $(\lambda_{De}/\lambda_{Di})^{2}\gg 1$, (3.15) exactly reduces
to the same linear dispersion relation (quartic in degree) in the original
form for the PMGC dynamics excited in DMCs obtained by Dwivedi et al. [10], as

$\displaystyle\left(\frac{\omega}{\omega_{J}}\right)^{4}+i\frac{\nu_{cn}}{\omega_{J}}\left(\frac{\omega}{\omega_{J}}\right)^{3}+\left(1-\frac{k^{2}C_{scam}^{2}}{\omega_{J}^{2}}\right)$
$\displaystyle\left(\frac{\omega}{\omega_{J}}\right)^{2}+i\frac{\nu_{cn}}{\omega_{J}}\left(\frac{\omega}{\omega_{J}}\right)$
(3.17)
$\displaystyle-\frac{k^{2}C_{scam}^{2}}{\omega_{J}^{2}}\left(\frac{1}{1+\eta}\right)=0,$

where, $\eta=n_{dc0}/n_{dn0}$ represents the ratio between the equilibrium
charged and neutral dust population densities.
$C_{scam}=\left[(q_{d}/e)^{2}(n_{dc0}T_{i})/(n_{dn0}m_{d})\right]^{1/2}$
refers to the usual phase speed of the supported dust acoustic wave (DAW) or
the so-called acoustic mode (SCAM)[10, 81].

To facilitate the analysis of the intricate PMGC fluctuations, a standard
astrophysical normalization procedure [80, 79] is employed to (3.18), which
converts the dispersion relation from previous wave space $(k,\omega)$ to a
new wave space $(K,\Omega)$ as follows

$\Omega^{4}+a_{3}\Omega^{3}+a_{2}\Omega^{2}+a_{1}\Omega+a_{0}=0;$ (3.18)
$a_{3}=i\left(\frac{\nu_{cn}}{\omega_{J}}\right),$ (3.19)
$a_{2}=\left(1-\frac{\chi
m_{d}n_{dn0}}{4C_{DA}^{2}}K^{2}-\frac{(1-R)\omega_{J}^{2}\lambda_{D_{r,q}}^{2}\Omega_{pdc}^{2}}{K^{2}\omega_{J}^{2}\lambda_{D_{r,q}}^{2}+C_{DA}^{2}}K^{2}\right)\left(K^{2}+\frac{1}{\xi^{2}}\right)K^{-2},$
(3.20)
$a_{1}=i\left(\frac{\nu_{cn}}{\omega_{J}}\right)\left[\left(1-\frac{\chi
m_{d}n_{d0}}{4C_{DA}^{2}}K^{2}\right)\left(K^{2}+\frac{1}{\xi^{2}}\right)K^{-2}\right],$
(3.21) $\displaystyle
a_{0}=-(1-R)\Omega_{pdc}^{2}\left(\frac{\lambda_{D_{r,q}}}{\lambda_{J}}\right)^{2}\frac{C_{DA}^{2}}{(K^{2}\omega_{J}^{2}\lambda_{D_{r,q}}^{2}+C_{DA}^{2})}$
$\displaystyle\left(1-\Omega_{Jc}^{2}-\frac{\chi
K^{2}}{4C_{DA}^{2}}m_{d}n_{dn0}\right)$ (3.22)
$\displaystyle\times\left(K^{2}+\frac{1}{\xi^{2}}\right)^{2}K^{-2}.$

Here, $\xi=\rho/\lambda_{J}$ and $K=k(C_{DA}/\omega_{J})$ are, respectively,
the dimensionless radial distance and angular wavenumber normalized with the
Jeans length $(\lambda_{J})$ and angular Jeans wavenumber
($k_{J}\equiv\omega_{J}/C_{DA}$, $C_{DA}$ being the dust acoustic phase
velocity). Similarly, $\Omega=\omega/\omega_{J}$,
$\Omega_{pdc}=\omega_{pdc}/\omega_{J}$, and
$\Omega_{Jc}=\omega_{Jc}/\omega_{J}$ are, respectively, the dimensionless
angular frequency, charged dust-plasma oscillation frequency, and critical
Jeans frequency for charged dust normalized with the Jeans angular frequency
$(\omega_{J})$.

Within the realm of ultra low-frequency (ULF) fluctuations, (3.19) transforms
to an equation characterized by a vanishing propagatory component
$(\Omega_{r}=0)$, signifying the wave condensation and a non-vanishing decay
(growth) component $(\Omega_{i}\neq 0)$, indicating wave collapse [18]. This
implies that the wave no longer exhibits propagatory characteristics in the
sense of classical wave mechanics. However, the wave undergoes either
exponential decay or growth over time, even in the absence of propagation.
Hence, the final dispersion relation after separating the real $(\Omega_{r})$
and imaginary $(\Omega_{i})$ parts can be written as

$\Omega_{i}=\frac{(1-R)(1+K^{2}\xi^{2})(-4C_{DA}^{2}\omega_{J}^{2}+\chi
m_{d}n_{d0}\omega_{J}^{2}K^{2}+4C_{DA}^{2}\omega_{Jc}^{2})C_{DA}^{2}\lambda_{D}^{2}\Omega_{pdc}^{2}}{\nu_{cn}\xi^{2}\omega_{J}(4C_{DA}^{2}-\chi
m_{d}n_{d0}K^{2})(C_{DA}^{2}+K^{2}\lambda_{D}^{2}\omega_{J}^{2})}.$ (3.23)

Equation (3.23) is the required linear dispersion relation for the PMGC
instability in a self-gravitating complex DMC due to the combined influence of
the polarization force and $(r,q)$-distributed electrons within the EiBI
gravity framework.

##  4 Results and discussions

A semi-analytic model is developed to assess the PMGC instability by
incorporating the effects of the EiBI gravity, $(r,q)$-distributed electrons,
and dust-polarization force. The derived quartic polynomial dispersion
relation (3.15) serves as an investigative tool for the oscillatory and
propagatory dynamics associated with the PMGC instability. As detailed in
section 3, the derived quartic dispersion relation (3.15) is reduced to a
linear dispersion relation using the approximation of extremely low-frequency
fluctuations. Finally, (3.23) is utilized for further numerical analysis, and
the resulting graphical representations (figures 2-10) are respectively
interpreted. Various input values of significant physical parameters are
adopted from reliable sources in the literature. It includes $m_{e}=9.1\times
10^{-28}$ g, $m_{i}=1.6\times 10^{-24}$ g, $m_{d}=4\times 10^{-12}$ g,
$e=4.80\times 10^{-10}$ esu, $q_{d}=-50e$, $r_{d}=1.5\times 10^{-4}$ cm,
$n_{e0}=1.2\times 10^{6}$ cm^-3, $n_{i0}=4.95\times 10^{6}$ cm^-3, $n_{dc0}=2$
cm^-3, $n_{dn0}=4$ cm^-3, $T_{e}=1$ eV, $T_{i}=0.8$ eV, and $T_{d}=0.1$ eV
[78, 24, 82, 83]. Accordingly, the Coulomb coupling parameter
$\Gamma_{\text{Cou}}$ is estimated as $\sim 0.011(\ll 1)$, which indicates a
weakly coupled DMC. Additionally, the Jeans angular frequency $(\omega_{J})$
is calculated as $4.49\times 10^{-9}$ rad . s^-1, the dust thermal velocity
$(V_{Td})$ is found as $8.9\times 10^{-2}$ cm . s^-1, and the dust-plasma
oscillation frequency $(\omega_{pd})$ is determined as $1.04$ rad . s^-1.
Furthermore, the effective dust-plasma Debye length for $(r,q)$-distribution
$(\lambda_{D})$ is obtained as $2.65\times 10^{-1}$ cm.

As in figure 2, we illustrate the impact of the polarization force on the PMGC
instability in the presence of $(r,q)$-distributed electrons within the
considered EiBI gravity framework. Here, the Jeans-scaled growth rate
$(\Omega_{i})$ of the PMGC instability is plotted against the Jeans-scaled
angular wavenumber $(K)$ for different values of polarization parameter
$(R=0,0.2,0.4,1.2,\text{ and }1.6)$. The EiBI parameter is kept fixed at
$\chi=2\times 10^{6}$ g^-1.cm^5. s^-2, the Jeans-scaled radial distance at
$\xi=30$, and the multi-order non-thermality coefficient at $A=1.4$ (for $r=0$
and $q=5)$. It demonstrates that for $R<1$, the attenuation of the PMGC
instability at shorter wavelengths (acoustic-like) is significantly more
prominent compared to the longer wavelengths (gravitational-like), whereas,
for $R>1$, the growth of the instability is more significant for shorter
wavelengths (acoustic-like), i.e., larger K-values. It is clear from figure 2
that an increase in the value of $R$ reduces the damping nature of the PMGC
instability, leading to a less stable system (most stable for $R=0$ and least
stable for $R=1.6$). This indicates that within the EiBI gravity framework,
the polarization force serves as a destabilizing factor for negatively charged
dust. Previous studies have also reported a similar destabilizing nature of
the polarization force on gravitational collapse due to the Jeans instability
in simplified model configurations [24, 78, 84]. It is important to note here
that, for the considered weakly coupled DMC, with the above-mentioned values
of the physical parameters, the polarization interaction parameter is found as
$R\sim 1.52\times 10^{-5}$, and hence the polarization force $(F_{p})$ should
not affect the stability of the system significantly.

Figure 2: Variation of the growth rate $(\Omega_{i})$ with the angular
wavenumber $(K)$ for different values of the polarization parameter $(R)$ on
the Jeans scale. The various curves link to: (a) $R=0$ (black solid line), (b)
$R=0.2$ (cyan dashed line), (c) $R=0.4$ (red dashed-dotted line), (d) $R=1.2$
(blue dotted line), and $R=1.6$ (magenta dotted pentagram line), respectively.

Figure 3 is same as figure 2, but portrayed in a wave space defined by Jeans-
scaled radial distance $(\xi)$ and Jeans-scaled angular wavenumber $(K)$. In
this scale of $\xi$, no significant effect over the PMGC fluctuation is
observed; however, changes in $R$-values distinctly affect the
$\Omega_{i}$-value. At $\xi=30$, the Jeans-scaled imaginary angular frequency
has a minimum value $(\Omega_{i}=-7.6\times 10^{-2})$ for $R=0$ and a maximum
value $(\Omega_{i}=4.07\times 10^{-2})$ for $R=1.6$, suggesting the growth of
the instability escalates with $R$. This supports the findings from figure 2.
The destabilizing effect of the polarization force and its underlying physics
can be deduced from the effective electrostatic force $(F_{\text{eff}})$
acting on the dust grains. The effective electrostatic force acting on the
dust grains subjected to an external electric field E is given by
$F_{\text{eff}}=(F_{e}+F_{p})=qE-
q_{d}^{2}\nabla\lambda_{D}/(2\lambda_{D}^{2})$ [22, 25]. Typically, the
polarization force has a minimal impact on system stability. However, it
becomes significant as the dust grain size increases [85, 25, 24]. Large-sized
dust grains tend to have higher charges $(q_{d})$, leading to an increased
$R$-value (since $R\propto q_{d}$). Thus, with the rise of the sheath-
polarization force, the counteracting electrostatic force that opposes the
inward gravitational pull weakens. Consequently, the gravitational force
becomes predominant over the electrostatic force, leading to system
instability. When $R\geq 1$, the dust polarization force exceeds the
corresponding electrostatic force. It causes the effective electrostatic force
$(F_{\text{eff}})$ to fail as a stabilizing force against the gravitational
force [86], ultimately triggering the self-gravitational collapse of the
interstellar cloud.

Figure 3: Same as figure 2, but in a wave space defined by Jeans-scaled radial
distance $(\xi)$ and Jeans-scaled angular wavenumber $(K)$. The different
panels correspond to (a) $R=0$, (b) $R=0.1$, (c) $R=0.2$, and (d) $R=1.6$,
respectively.

In a similar way, figure 4 depicts a color-spectral top-view projection of
figure 3. This graphic projection is utilized for a more detailed examination
of the conjugated relationship of $\Omega_{i}$ with $\xi$. Interestingly, it
is found that in the vicinity of the central region of the DMC $(\xi\sim 0)$,
the variation of $\Omega_{i}$ suddenly appears to disappear. This modal
behaviour indicates the possible signatures of a microphysical correlation
depictable with $\Omega_{i}=f(\xi)$.

Figure 4: Same as figure 3, but in a color-spectral top-view projection.

Figure 5 displays similar features as figure 3, but shown in a magnified view
of $\xi$ and an extended range of $K$ (encompassing the intervals $\xi=0-3$
and $K=0-50$). Apart from the polarization effects noticed in figures 2 \- 4,
a distinct spatial variation in $\Omega_{i}$ is observed in the low-$\xi$
region $(\sim 0-0.2)$. It indicates that in the proximity of the spherical DMC
core, the PMGC instability shows a rapidly escalating tendency for both $R<1$
and $R\geq 1$. However, beyond $\xi\sim 0.2$, the instability reaches a
plateau and spreads uniformly across the system. This plateau region can be
interpreted as the result of the interaction between gravity and effective
electrostatic force. As $\xi$ continues to increase, the effective
electrostatic force within the cloud becomes sufficient enough to oppose the
gravitational collapse of the DMC.

Figure 5: Same as figure 3, but in a zoomed-in form $(\xi=0-3$, and $K=0-50)$.

Figure 6 illustrates the EiBI gravity effect on the PMGC instability in the
presence of polarization force and $(r,q)$-distributed electrons. For any
self-gravitationally bounded astrophysical object of size $L$, the value of
$\chi$ should be less than $GL^{2}$, $G$ being the universal gravitational
constant [65]. For our system having modified critical Jeans length,
$\lambda_{J}\sim 1.37\times 10^{8}$ cm, the constraint is found to be
$\chi<1.33\times 10^{9}$ g^-1.cm^5.s^-2. The $\chi$-value is varied in the
range from $-4\times 10^{6}\ \unit{g^{-1}.cm^{5}.s^{-2}}\text{ to }4\times
10^{6}$ g^-1.cm^5.s^-2. It is noteworthy here that, $\chi=0$ serves as the
Newtonian reference in the non-local gravity, which is in fact, a crucial
feature for assessing the observed trends. All the other parameters are kept
the same as mentioned at the beginning of this section, except $R=0.2$ and
$A=1.4\ (r=0,q=5)$.

Figure 6: Same as figure 2, but for different values of the EiBI parameter
$(\chi)$. The various curves respectively link to: (a) $\chi=0$ (black solid
line, taken as the Newtonian reference line), (b) $\chi=2\times 10^{6}$
g^-1.cm^5.s^-2 (cyan dashed line), (c) $\chi=4\times 10^{6}$ g^-1.cm^5.s^-2
(red dashed-dotted line), (d) $\chi=-2\times 10^{6}$ g^-1.cm^5.s^-2 (blue
dotted line), and (e) $\chi=-4\times 10^{6}$ g^-1.cm^5.s^-2 (magenta dotted
pentagram line).

It is clear from figure 6 that an elevation in the negative EiBI parameter
$(\chi)$ results in the deceleration of the PMGC instability, thereby
enhancing the stability of the system in comparison to those in the Newtonian
scenarios. Conversely, an increase in the positive EiBI parameter $(\chi)$
drives the system towards a state of higher instability. Therefore, it can be
deduced that the positive EiBI parameter functions as a destabilizing element,
while the negative EiBI parameter acts as a stabilizing factor for the system
being analyzed. Previous literature has also reported a similar impact of the
EiBI gravity on the Jeans instability, aligning with the findings of this
study [57].

Figure 7: Same as figure 6, but in a wave space defined by Jeans-scaled radial
distance $(\xi)$ and Jeans-scaled angular wavenumber $(K)$. The different
panels correspond to (a) $\chi=0$ (taken as the Newtonian reference), (b)
$\chi=2\times 10^{6}$ g^-1.cm^5.s^-2, and (c) $\chi=-2\times 10^{6}$
g^-1.cm^5.s^-2, respectively.

Figure 7 depicts similar characteristics to those in figure 6. In this
representation, the spatial aspects of the PMGC instability are illustrated by
jointly varying the Jeans-scaled radial distance $(\xi)$ and Jeans-scaled
wavenumber $(K)$. Here, only two values of the EiBI parameter ($\chi=\pm\
2\times 10^{6}$ g^-1.cm^5.s^-2) are examined alongside the Newtonian case
$(\chi=0)$. All other pertinent physical parameters remain unchanged from
those utilized in figure 6. Similar to the findings in figure 6, it is noted
that a positive (negative) EiBI parameter $(\chi)$ acts as a destabilizing
(stabilizing) factor in the presence of non-thermal $(r,q)$-distributed
electrons and sheath-polarization of the Debye sheaths around the dust grains
present in the DMC. Additionally, a color-spectral top-view projection of the
same is shown in figure 8.

Figure 8: Same as figure 7, but in a colour-spectral top-view projection.

Figure 9 elucidates the influence of $(r,q)$-distributed electrons on the PMGC
instability within the EiBI gravity framework, taking the polarization force
into account. The variation of Jeans-scaled imaginary angular frequency
$(\Omega_{i})$ with the Jeans-scaled angular wavenumber $(K)$ is depicted for
different values of $q\ (=3,5,10)$ with (a) $r=0$ and (b) $r=1$. In the
limiting condition $r=0$ and $q\to\infty$, the value of $A$ converges to $1$.
Hence, the line corresponding to $r=0$ and $q\to\infty$ serves as the
Maxwellian reference. This numerical analysis maintains consistent fixed
equilibrium values as used before except $R=0.2$ and $\chi=2\times 10^{6}$
g^-1.cm^5.s^-2.

It is observed in figure 9a that for $r=0$, the different curves exhibit a
compelling propensity to converge towards the Maxwellian reference as the
value of $q$ increases. This outcome is completely inline with the previous
prediction founded on the $(r,q)$-distribution laws. Interestingly, the
results diverge significantly from the previous case when $r=1$. It is seen in
figure 9b the $q$-value increments deviate the curves away from the Maxwellian
reference, indicating the emergence of non-Maxwellian behavior. Moreover, it
is evident from figure 9 that $q$-value increments dampen the instability of
the system for both the cases $(r=0,1)$. Nevertheless, for $r=0$, the
Maxwellian system is observed to be more stable than the non-Maxwellian
system, whereas, for $r=1$, the non-Maxwellian system is more stable.

Figure 9: Same as figure 2, but for different values of $q$ with (a) $r=0$ and
(b) $r=1$. The various curves respectively link to : (a) $r=0,\ q\to\infty$
(black solid line, taken as the Maxwellian reference line), (b) $q=3$ (cyan
dashed line), (c) $q=5$ (red dashed-dotted line), and (d) $q=10$ (blue dotted
line).

Figure 10: Same as figure 9, but for different values of $r$ with (a) $q=2$
and $q=6$. The various curves respectively link to: (a) $r=0,\ q\to\infty$
(black solid line, taken as the Maxwellian reference line), (b) $r=2$ (cyan
dashed line), (c) $r=4$ (red dashed-dotted line), and (d) $r=6$ (blue dotted
line).

In a similar way, in figure 10, we portray the same as figure 9, but for
varying values of $r\ (=2,4,6)$ with (a) $q=2$ and (b) $q=6$. As figure 9, it
is evident that an increase in the $r$ spectral parameter enhances the system
stability accompanied by a rise in the non-Maxwellian (non-thermal)
characteristics. A comparative analysis of figures 9-10 reveals a consistent
trend. Both the non-thermal spectral parameters ($r$ and $q$) play a
significant role in influencing the system stability and non-Maxwellian
characteristics. While elevating the $q$ value consistently improves stability
across different $r$ values and vice versa, it also steers the system away
from the Maxwellian velocity distribution. Additionally, it is notable that at
longer wavelengths (gravitational-like), the system maintains closer proximity
to the Maxwellian characteristics, whereas, at shorter wavelengths (acoustic-
like), the deviation (from thermality) becomes more pronounced. The stability
of the non-Maxwellian (non-thermal) systems can be comprehended by noting that
as the electron energy increases, their mobility also rises and vice versa.
Consequently, there is a rapid accumulation of electrons at the dust grain
surfaces, leading to an escalation in the dust charge. The enhanced dust
charge at the cost of thermal loss of electrons results in an increased
electrostatic repulsion among the dust grains, counterbalancing the
gravitational forces acting inward.

##  5 Conclusions

The proposed study explores the collective impact of the EiBI gravity,
$(r,q)$-distributed electrons, and dust-polarization force on the pulsational
mode of gravitational collapse (PMGC) in complex spherical dust molecular
cloud (DMC) fluids. Application of a standard spherical normal mode analysis
in the EiBI-modified DMC results in a unique form of a generalized linear
quartic dispersion relation. The derived linear dispersion relation, if all
the newly added astronomical sophistication are switched off, fairly
corroborates with the previous results reported in the literature [10]. After
its necessary validation checkup, the derived PMGC dispersion relation is
numerically illustrated with judicious parametric inputs to explore diverse
PMGC stability features in different realistic astronomical circumstances.

The PMGC instability is found to be significantly more pronounced for longer
wavelengths (gravitational-like) than the shorter wavelengths (acoustic-like).
Additionally, the instability exhibits a rapidly growing tendency in the
vicinity of the DMC core. It shows a saturating propensity from the cloud
center outwards. It is seen that the polarization force and the positive EiBI
parameter act as destabilizing factors, while the negative EiBI parameter acts
in favour of the DMC stability. The non-thermal $(r,q)$-distributed electrons
produce a DMC stabilizing influence. It is important to note that the
deviation from the Maxwellian (thermal) characteristics is significant at
shorter wavelengths (acoustic-like), whereas, at longer wavelengths
(gravitation-like), there is no substantial distinction between the Maxwellian
(thermal) and non-Maxwellian (non-thermal) features in the thermo-statistical
perspective. However, the non-Maxwellian systems are found to be more stable
than the Maxwellian ones.

The outcomes of the PMGC instability analyses presented in this work could be
utilized to enhance the understanding the diverse processes of astrophysical
structure formation in the H ii regions (ionized hydrogen regions) of the
molecular clouds, such as Sh2-87 [87], Sh2-88B [88], Sh2-235 [89], etc. The
Jeans-like instability in the EiBI gravity framework discussed in this study
may lead to the fragmentation of extremely dense interstellar DMCs into
smaller substructures that cannot withstand their own gravitational forces.
These results may be helpful in further probing the EiBI gravity role in
various other similar structures. It may include various star forming
molecular cloud zones in our galaxy, such as the Central Molecular Zone [90],
the Orion molecular cloud [91], etc. While there have not been any substantial
observational studies of the PMGC model reported yet, the recently launched
JWST enables us to observe star-forming regions beyond the Milky Way and its
satellite galaxies [92]. Already, the JWST has uncovered galaxies with $z>10$
($z$ being the redshift). Additionally, the Atacama Large Millimeter Array
(ALMA) facilitates the observation and analysis of gas and star formation
process in galaxies with $z>4$, on scales less than a kiloparsec [93]. These
developments could open up new avenues for future investigations in light of
the real-time insightful observations in the relevant PMGC research areas.

It is important to note here that the effects, like viscoelasticity (in a
strongly coupled system), dust charge fluctuation, ion drag, turbulence, and
magnetic field are not considered in this work. Additionally, there are
numerous other extensions of general relativity as well as several other
velocity distributions, the effects of which on the behaviour of instabilities
are yet to be explored. Therefore, it is conclusively admitted that there is
potential scope for future improvements to the presented model by
incorporating such important plasma-fluidic and thermo-statistical features of
realistic astronomical significance.

## Acknowledgments

The Department of Physics at Tezpur University is gratefully acknowledged for
the invaluable cooperation. Special appreciations are due to all the members
of the Astrophysical Plasma and Nonlinear Dynamics Research Laboratory
(APNDRL) for their collaborative efforts and insightful discussions, which
have enriched the intellectual environment of our research pursuits. Special
acknowledgment is given to Mr. Souvik Das of Department of Physics, Tezpur
University for his meticulous attention to detail and unwavering support in
manuscript formatting. Finally, Mr. Siddharth Saikia, affiliated to University
of New South Wales, is also hereby acknowledged for his help in writing
elementary parts of the manuscript.

## Appendix: Polarization force for negatively charged dust

The expression for the polarization force is given by [25, 22]

$F_{p}=-\frac{q_{d}^{2}\nabla\lambda_{D}}{2\lambda_{D}^{2}};$ (5.1)

where,

$\lambda_{D}=\lambda_{Di}\times\left(1+\frac{\lambda_{Di}^{2}}{\lambda_{Di}^{2}}\right)^{-1/2}=\left(\frac{k_{B}T_{i}}{4\pi
e^{2}n_{i}}\right)^{1/2}\times\left(\frac{n_{i}T_{e}}{n_{i}T_{e}+n_{e}T_{i}}\right)^{1/2}.$
(5.2)

For a complex (dusty) plasma comprising of negatively charged dust grains,
$n_{i}\gg n_{e}$ and $T_{e}\gg T_{i}$, hence we can judiciously use the
approximation $n_{e}T_{i}\ll n_{i}T_{e}$ in (5.2)

$\lambda_{D}=\left(\frac{k_{B}T_{i}}{4\pi
e^{2}n_{i}}\right)^{1/2}\approx\lambda_{Di}\hskip 2.84526pt.$ (5.3)

Now,

$\displaystyle\nabla\lambda_{D}$
$\displaystyle=\nabla\lambda_{Di}=\nabla\left(\frac{k_{B}T_{i}}{4\pi
e^{2}n_{i}}\right)^{1/2}$ (5.4)
$\displaystyle=-\frac{1}{2}\left(\frac{k_{B}T_{i}}{4\pi
e^{2}}\right)^{1/2}n_{i}^{-3/2}\nabla n_{i}\hskip 2.84526pt.$

Therefore,

$\displaystyle\frac{\nabla\lambda}{\lambda_{D}^{2}}$
$\displaystyle=-\frac{1}{2}\left(\frac{k_{B}T_{i}}{4\pi
e^{2}}\right)^{1/2}n_{i}^{-3/2}\nabla
n_{i}\times\left(\frac{1}{\lambda_{D}^{2}}\right)$ (5.5)
$\displaystyle=-\frac{1}{2}\left(\frac{k_{B}T_{i}}{4\pi
e^{2}n_{i}}\right)^{1/2}n_{i}^{-1}\nabla
n_{i}\times\left(\frac{1}{\lambda_{D}^{2}}\right)$
$\displaystyle=-\frac{1}{2}\lambda_{D}n_{i}^{-1}\nabla
n_{i}\times\left(\frac{1}{\lambda_{D}^{2}}\right)$
$\displaystyle=-\frac{1}{2\lambda_{D}}\left(\frac{\nabla
n_{i}}{n_{i}}\right).$

Considering the ions to be Maxwellian, we can use $T_{i}\nabla
n_{i}=-n_{i}e\nabla\phi$ in (5.5) as

$\frac{\nabla\lambda}{\lambda_{D}^{2}}=-\frac{1}{2\lambda_{D}}\left(\frac{-n_{i}e\nabla\phi}{n_{i}T_{i}}\right)=\frac{1}{2\lambda_{D}}\left(\frac{e\nabla\phi}{T_{i}}\right).$
(5.6)

Finally, substituting $\nabla\lambda_{D}/\lambda_{D}^{2}$ in (5.1), we get

$\displaystyle F_{p}$
$\displaystyle=-\frac{q_{d}^{2}}{2}\times\frac{1}{2\lambda_{D}}\left(\frac{e\nabla\phi}{T_{i}}\right)$
(5.7) $\displaystyle=-\frac{q_{d}^{2}}{2}\times\frac{1}{2}\left(\frac{4\pi
n_{i}e^{2}}{k_{B}T_{i}}\right)^{\frac{1}{2}}\times\left(\frac{e\nabla\phi}{T_{i}}\right)$
$\displaystyle=-\frac{q_{d}^{2}}{4}\times\left(\frac{4\pi
n_{i0}e^{2}}{k_{B}T_{i}}\right)^{\frac{1}{2}}\times\left(\frac{n_{i}}{n_{i0}}\right)^{\frac{1}{2}}\times\left(\frac{e\nabla\phi}{T_{i}}\right)$
$\displaystyle=-q_{d}\times\left(\frac{q_{d}e}{4T_{i}\lambda_{Di0}}\right)\times\left(\frac{n_{i}}{n_{i0}}\right)^{\frac{1}{2}}\times\nabla\phi=-q_{d}R\left(\frac{n_{i}}{n_{i0}}\right)^{\frac{1}{2}}\nabla\phi,$

where, $R=\left(q_{d}e/4T_{i}\lambda_{Di0}\right)$ is known as the
polarization interaction parameter.

## References

  * [1] B.P. Pandey, J. Vranješ, S. Poedts and P.K. Shukla, _The pulsational mode in the presence of dust charge fluctuations_ , _Physica Scripta_ 65 (2002) 513.
  * [2] J.H. Jeans, _The stability of a spherical nebula_ , _Phil. Trans. A. Math. Phys. Eng. Sci._ 199 (1902) 1.
  * [3] D. Prialnik, _An Introduction to the Theory of Stellar Structure and Evolution_ , Cambridge University Press (2009).
  * [4] A.R. Choudhuri, _Astrophysics for physicists_ , Cambridge University Press (2010).
  * [5] F.F. Chen et al., _Introduction to plasma physics and controlled fusion_ , vol. 1, Springer (1984).
  * [6] A.R. Choudhuri, _The physics of fluids and plasmas: an introduction for astrophysicists_ , Cambridge University Press (1998).
  * [7] P.P. Eggleton, _Chandrasekhar’s book an introduction to the study of stellar structure_ , _Pramana_ 77 (2011) 97.
  * [8] F. Hoyle, _On the fragmentation of gas clouds into galaxies and stars._ , _Astrophysical Journal, vol. 118, p. 513_ 118 (1953) 513.
  * [9] K. Avinash, B. Eliasson and P. Shukla, _Dynamics of self-gravitating dust clouds and the formation of planetesimals_ , _Physics Letters A_ 353 (2006) 105.
  * [10] C. Dwivedi, A. Sen and S. Bujarbarua, _Pulsational mode of gravitational collapse and its impact on the star formation_ , _Astronomy and Astrophysics, v. 345, p. 1049-1053 (1999)_ 345 (1999) 1049.
  * [11] K. Devi, A. Nag, J. Paul, U. Sarma and P.K. Karmakar, _Nonlinear pulsational mode dynamics in strongly correlated dust molecular clouds_ , _Journal of Taibah University for Science_ 17 (2023) 2202605.
  * [12] P.K. Karmakar and A. Haloi, _Pulsational instability of complex charge-fluctuating magnetized turbulent astroclouds_ , _Astrophysics and Space Science_ 362 (2017) 1.
  * [13] B. Borah and P. Karmakar, _Pulsational mode fluctuations and their basic conservation laws_ , _Advances in Space Research_ 55 (2015) 416.
  * [14] P. Dutta and P.K. Karmakar, _Viscoelastic pulsational mode_ , _Astrophysics and Space Science_ 362 (2017) 1.
  * [15] A. Haloi and P. Karmakar, _Nonlinear gravito-electrostatic waves in self-gravitating complex plasma in presence of ion-drag effects_ , _Astrophysics and Space Science_ 358 (2015) 1.
  * [16] W.-H. Yang, H. Chen and S.-Q. Liu, _Effect of ion drag on a pulsational mode of gravitational collapse_ , _Communications in Theoretical Physics_ 72 (2020) 075504.
  * [17] B. Bhakta, M. Gohain and P.K. Karmakar, _Pulsational mode behaviour in complex non-extensive non-thermal viscous astrophysical fluids_ , _Europhysics Letters_ 119 (2017) 25001.
  * [18] P.K. Karmakar and P. Dutta, _Evolutionary pulsational mode dynamics in nonthermal turbulent viscous astrofluids_ , _Astrophysics and Space Science_ 362 (2017) 203.
  * [19] P.K. Karmakar and J. Dutta, _Athermal pulsational mode dynamics with negative ions_ , _Chinese Journal of Physics_ (2024) .
  * [20] M. Das and P.K. Karmakar, _Dynamics of pulsational mode in the eibi gravity fabric_ , _Chinese Journal of Physics_ 87 (2024) 232.
  * [21] B. Borah, A. Haloi and P. Karmakar, _Atypical gravito-electrostatic fluctuations in the presence of active ion-inertial dynamics_ , _Journal of Plasma Physics_ 82 (2016) 905820206.
  * [22] S. Hamaguchi and R.T. Farouki, _Polarization force on a charged particulate in a nonuniform plasma_ , _Phys. Rev. E_ 49 (1994) 4430.
  * [23] S. El-Labany, W. El-Taibany, A. El-Tantawy and N. Zedan, _Effects of double spectral electron distribution and polarization force on dust acoustic waves in a negative dusty plasma_ , _contributions to plasma physics_ 60 (2020) e202000049.
  * [24] B. Dolai and R. Prajapati, _Effects of dust-charge gradient and polarization forces on the waves and jeans instability in strongly coupled dusty plasma_ , _Physics Letters A_ 384 (2020) 126462.
  * [25] S.A. Khrapak, A.V. Ivlev, V.V. Yaroshenko and G.E. Morfill, _Influence of a polarization force on dust acoustic waves_ , _Phys. Rev. Lett._ 102 (2009) 245004.
  * [26] K. Ashrafi, A. Mamun and P. Shukla, _Polarization force for different dusty plasma situations_ , _Journal of Plasma Physics_ 80 (2014) 1.
  * [27] P. Karmakar, _Nonlinear stability of pulsational mode of gravitational collapse in self-gravitating hydrostatically bounded dust molecular cloud_ , _Pramana_ 76 (2011) 945.
  * [28] W. Feldman, R. Anderson, S. Bame, S. Gary, J. Gosling, D. McComas et al., _Electron velocity distributions near the earth’s bow shock_ , _Journal of Geophysical Research: Space Physics_ 88 (1983) 96.
  * [29] D. Zhao, S. Fu, G.K. Parks, W. Sun, Q. Zong, D. Pan et al., _Electron flat-top distributions and cross-scale wave modulations observed in the current sheet of geomagnetic tail_ , _Physics of Plasmas_ 24 (2017) .
  * [30] Y. Asano, R. Nakamura, I. Shinohara, M. Fujimoto, T. Takada, W. Baumjohann et al., _Electron flat-top distributions around the magnetic reconnection region_ , _Journal of Geophysical Research: Space Physics_ 113 (2008) .
  * [31] G. Livadiotis, D. McComas, M. Dayeh, H. Funsten and N.A. Schwadron, _First sky map of the inner heliosheath temperature using ibex spectra_ , _The Astrophysical Journal_ 734 (2011) 1.
  * [32] Y. Zhang, X.-W. Liu, R. Wesson, P. Storey, Y. Liu and I. Danziger, _Electron temperatures and densities of planetary nebulae determined from the nebular hydrogen recombination spectrum and temperature and density variations_ , _Monthly Notices of the Royal Astronomical Society_ 351 (2004) 935.
  * [33] S. Krimigis, J. Carbary, E. Keath, T. Armstrong, L. Lanzerotti and G. Gloeckler, _General characteristics of hot plasma and energetic particles in the saturnian magnetosphere: Results from the voyager spacecraft_ , _Journal of Geophysical Research: Space Physics_ 88 (1983) 8871.
  * [34] V.M. Vasyliunas, _A survey of low-energy electrons in the evening sector of the magnetosphere with ogo 1 and ogo 3_ , _Journal of Geophysical Research_ 73 (1968) 2839.
  * [35] J. Gosling, J. Asbridge, S. Bame, W. Feldman, R. Zwickl, G. Paschmann et al., _Interplanetary ions during an energetic storm particle event: The distribution function from solar wind thermal energies to 1.6 mev_ , _Journal of Geophysical Research: Space Physics_ 86 (1981) 547.
  * [36] C. Tsallis, _Possible generalization of boltzmann-gibbs statistics_ , _Journal of statistical physics_ 52 (1988) 479.
  * [37] R. Cairns, A. Mamum, R. Bingham, R. Boström, R. Dendy, C. Nairn et al., _Electrostatic solitary structures in non-thermal plasmas_ , _Geophysical Research Letters_ 22 (1995) 2709.
  * [38] M. Qureshi, H. Shah, G. Murtaza, S. Schwartz and F. Mahmood, _Parallel propagating electromagnetic modes with the generalized (r, q) distribution function_ , _Physics of Plasmas_ 11 (2004) 3819.
  * [39] S. Zaidi, M. Qureshi and S. Khalid, _Effect of nonthermal electron distributions on dust acoustic solitons in cometary plasmas_ , _AIP Advances_ 14 (2024) .
  * [40] F. Dyson, A. Eddington and C. Davidson, _A determination of the deflection of light by the sun’s gravitational field, from observations made at the total eclipse of may 29, 1919_ , _Memoirs of the Royal Astronomical Society, Vol. 62, p. A1_ 62 (1923) A1.
  * [41] I.I. Shapiro, _Solar system tests of general relativity: recent results and present plans_ , _General relativity and gravitation_ (1990) 313.
  * [42] C.M. Will, _The confrontation between general relativity and experiment_ , _Living reviews in relativity_ 17 (2014) 1.
  * [43] I.H. Stairs, _Testing general relativity with pulsar timing_ , _Living Reviews in Relativity_ 6 (2003) 1.
  * [44] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _Binary black hole mergers in the first advanced ligo observing run_ , _Phys. Rev. X_ 6 (2016) 041015.
  * [45] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _Observation of gravitational waves from a binary black hole merger_ , _Phys. Rev. Lett._ 116 (2016) 061102.
  * [46] A. Hees, T. Do, A.M. Ghez, G.D. Martinez, S. Naoz, E.E. Becklin et al., _Testing general relativity with stellar orbits around the supermassive black hole in our galactic center_ , _Phys. Rev. Lett._ 118 (2017) 211101.
  * [47] S. Banerjee, S. Shankar and T.P. Singh, _Constraints on modified gravity models from white dwarfs_ , _Journal of Cosmology and Astroparticle Physics_ 2017 (2017) 004.
  * [48] P. Banerjee, D. Garain, S. Paul, R. Shaikh and T. Sarkar, _A stellar constraint on eddington-inspired born–infeld gravity from cataclysmic variable binaries_ , _The Astrophysical Journal_ 924 (2022) 20.
  * [49] S. Jana, G.K. Chakravarty and S. Mohanty, _Constraints on born-infeld gravity from the speed of gravitational waves after gw170817 and grb 170817a_ , _Phys. Rev. D_ 97 (2018) 084011.
  * [50] J.W. Moffat, _Scalar–tensor–vector gravity theory_ , _Journal of Cosmology and Astroparticle Physics_ 2006 (2006) 004.
  * [51] D.c.v. Glavan and C. Lin, _Einstein-gauss-bonnet gravity in four-dimensional spacetime_ , _Phys. Rev. Lett._ 124 (2020) 081301.
  * [52] C. Zanoletti, B.R. Hull, C.D. Leonard and R.B. Mann, _Cosmological constraints on 4-dimensional einstein-gauss-bonnet gravity_ , _Journal of Cosmology and Astroparticle Physics_ 2024 (2024) 043.
  * [53] A.A. Starobinsky, _A new type of isotropic cosmological models without singularity_ , _Physics Letters B_ 91 (1980) 99.
  * [54] M. Bañados and P.G. Ferreira, _Eddington’s theory of gravity and its progeny_ , _Phys. Rev. Lett._ 105 (2010) 011101.
  * [55] E.M. Lifshitz, _On the gravitational stability of the expanding universe_ , _J. Phys._ 10 (1946) 116.
  * [56] Q. Yang, L. Tan, H. Chen and S. Liu, _Jeans instability analysis of viscoelastic astrofluids in eddington-inspired-born–infeld (eibi) gravity_ , _New Astronomy_ 99 (2023) 101947.
  * [57] W.-H. Yang, Y.-Z. Xiong, H. Chen and S.-Q. Liu, _Jeans gravitational instability with $\kappa$-deformed kaniadakis distribution in eddington-inspired born–infield gravity_, _Chinese Physics B_ 29 (2020) 110401.
  * [58] J. Vainio and I. Vilja, _Jeans analysis of bok globules in f (r) gravity_ , _General Relativity and Gravitation_ 48 (2016) 129.
  * [59] K.-R. He, _Jeans analysis with $\kappa$-deformed kaniadakis distribution in f (r) gravity_, _Physica Scripta_ 97 (2022) 025601.
  * [60] A. Hazarika, M. Das and P.K. Karmakar, _Polytropic spherical astrocosmic fluid stability in the eibi gravity with strong collective correlative effects_ , _Chinese Journal of Physics_ (2024) .
  * [61] P. Pani, T. Delsate and V. Cardoso, _Eddington-inspired born-infeld gravity: Phenomenology of nonlinear gravity-matter coupling_ , _Phys. Rev. D_ 85 (2012) 084020.
  * [62] J. Casanellas, P. Pani, I. Lopes and V. Cardoso, _Testing alternative theories of gravity using the sun_ , _The Astrophysical Journal_ 745 (2011) 15.
  * [63] P. Pani, V. Cardoso and T. Delsate, _Compact stars in eddington inspired gravity_ , _Phys. Rev. Lett._ 107 (2011) 031101.
  * [64] T. Harko, F.S. Lobo and M.K. Mak, _Bianchi type i cosmological models in eddington-inspired born–infeld gravity_ , _Galaxies_ 2 (2014) 496.
  * [65] P.P. Avelino, _Eddington-inspired born-infeld gravity: Astrophysical and cosmological constraints_ , _Phys. Rev. D_ 85 (2012) 104053.
  * [66] T. Clifton, P.G. Ferreira, A. Padilla and C. Skordis, _Modified gravity and cosmology_ , _Physics reports_ 513 (2012) 1.
  * [67] B. Schutz, _A first course in general relativity_ , Cambridge university press (2022).
  * [68] S. Deser and G. Gibbons, _Born-infeld-einstein actions?_ , _Classical and Quantum Gravity_ 15 (1998) L35.
  * [69] D.N. Vollick, _Palatini approach to born-infeld-einstein theory and a geometric description of electrodynamics_ , _Phys. Rev. D_ 69 (2004) 064030.
  * [70] D.N. Vollick, _Born-infeld-einstein theory with matter_ , _Phys. Rev. D_ 72 (2005) 084026.
  * [71] P. Pani and T.P. Sotiriou, _Surface singularities in eddington-inspired born-infeld gravity_ , _Phys. Rev. Lett._ 109 (2012) 251102.
  * [72] P.K. Shukla and A.A. Mamun, _Introduction to dusty plasma physics_ , _Plasma Physics and Controlled Fusion_ 44 (2002) 395.
  * [73] F.H. Shu, F.C. Adams and S. Lizano, _Star formation in molecular clouds: observation and theory_ , _Annual review of astronomy and astrophysics_ 25 (1987) 23.
  * [74] M.R. Gupta, S. Sarkar, S. Ghosh, M. Debnath and M. Khan, _Effect of nonadiabaticity of dust charge variation on dust acoustic waves: Generation of dust acoustic shock waves_ , _Phys. Rev. E_ 63 (2001) 046406.
  * [75] C. Jian-Hong and W. Nan-Xia, _Effects of adiabatic dust charge fluctuation and particles collisions on dust-acoustic solitary waves in three-dimensional magnetized dusty plasmas_ , _Communications in Theoretical Physics_ 51 (2009) 524.
  * [76] D. Ward-Thompson and A.P. Whitworth, _An introduction to star formation_ , Cambridge University Press (2011).
  * [77] S. Ali, W. Masood, H. Rizvi, R. Jahangir and A.M. Mirza, _Contribution of the generalized (r, q) distributed electrons in the formation of nonlinear ion acoustic waves in upper ionospheric plasmas_ , _AIP Advances_ 11 (2021) .
  * [78] P. Dutta, P. Das and P. Karmakar, _Stability analysis of non-thermal complex astrofluids in the presence of polarized dust-charge fluctuations_ , _Astrophysics and Space Science_ 361 (2016) 322.
  * [79] S. Dasgupta and P.K. Karmakar, _The jeans instability in viscoelastic spherical astrophysical fluid media_ , _Astrophysics and Space Science_ 364 (2019) 1.
  * [80] D. Kalita and P.K. Karmakar, _Jeans instability in nonideal mhd plasma clouds with geometric curvature effects_ , _IEEE Transactions on Plasma Science_ 49 (2021) 2042.
  * [81] C. Dwivedi, _Is dust acoustic wave a new plasma acoustic mode?_ , _Physics of Plasmas_ 4 (1997) 3427.
  * [82] P.K. Singha and P.K. Karmakar, _Bimodal instability dynamics in nonthermal complex tridust astroclouds towards structure formation_ , _Europhysics Letters_ 128 (2020) 35001.
  * [83] P.K. Shukla and L. Stenflo, _Jeans instability in a self-gravitating dusty plasma_ , _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 462 (2006) 403.
  * [84] R. Prajapati, _Effect of polarization force on the jeans instability of self-gravitating dusty plasma_ , _Physics Letters A_ 375 (2011) 2624.
  * [85] M. Asaduzzaman and A. Mamun, _Roles of polarization force and nonthermal electron on dust-acoustic waves in an inhomogeneous dusty plasma with positively charged dust_ , _Physics of Plasmas_ 19 (2012) .
  * [86] H. Chen, S. Zhou, R. Luo and S. Liu, _Nonlinear dust acoustic waves with polarization force effects in kappa distribution plasma_ , _Japanese Journal of Applied Physics_ 56 (2016) 016101.
  * [87] S. Biswas, S. Mondal, A. Hoque, A. Panja, T. Baug and R. Das, _Star formation in the h ii region sh2-87: Evidence of global hierarchical collapse_ , _The Astrophysical Journal_ 961 (2024) 176.
  * [88] L. Deharveng, D. Nadeau, A. Zavagno and J. Caplan, _The stellar content of the compact h ii region sh2-88b_ , _Astronomy and Astrophysics-A &A_ 360 (2000) 1107.
  * [89] M. Kirsanova, P. Boley, A. Moiseev, D. Wiebe and R. Uklein, _3d structure of the h ii region sh2-235 from tunable-filter optical observations_ , _Monthly Notices of the Royal Astronomical Society_ 497 (2020) 1050.
  * [90] X. Lu, Y. Cheng, A. Ginsburg, S.N. Longmore, J.D. Kruijssen, C. Battersby et al., _Alma observations of massive clouds in the central molecular zone: Jeans fragmentation and cluster formation_ , _The Astrophysical Journal Letters_ 894 (2020) L14.
  * [91] A. Palau, L.A. Zapata, C.G. Román-Zúñiga, Á. Sánchez-Monge, R. Estalella, G. Busquet et al., _Thermal jeans fragmentation within 1000 au in omc-1s_ , _The Astrophysical Journal_ 855 (2018) 24.
  * [92] J. Peltonen, E. Rosolowsky, T.G. Williams, E.W. Koch, A. Dolphin, J. Chastenet et al., _Jwst reveals star formation across a spiral arm in m33_ , _Monthly Notices of the Royal Astronomical Society_ 527 (2024) 10668.
  * [93] J. Freundlich, _Star formation across cosmic time_ , _Fundamental Plasma Physics_ 11 (2024) 100059.