1 Introduction

The Jeans instability of sound modes in self-gravitating interstellar clouds is a well known phenomenon which plays a strategic role in understanding not only stars and galaxies formation, but also as a possible justification for the existence of large amount of dark energy and dark matter (DM) in the Universe.

Thereby such a topic represents a long standing problem whose newness is strictly related to the development of several new mathematical models to better fit the interstellar matter [1,2,3,4,5,6,7,8,9]. In particular, as already suggested in earlier papers [4, 10,11,12,13], an interesting research issue may be just given by a different structure for the standard Poisson-type gravity, leading to modified gravitational (MOG) frameworks accounting for the Universe Expansion, via the introduction of Einstein’s cosmological constant \({\varLambda }\).

In order to overcome the well-known Jeans swindle and predict dark energy properties, but also thinking of the recently detected gravitational waves [14], we here propose a non-stationary and non-relativistic Einstein’s equation, which includes \({\varLambda }\) and a (small) gravitational relaxation time. To this aim we introduce a reaction–diffusion equation for the gravitational potential and firstly test its stabilizing role on self-gravitating interstellar matter (ISM) clouds, described as polytropic Euler-Darcy fluids. It follows that this MOG theory accounts for a lower critical wavenumber, here called the Jeans-Einstein wavenumber, \(k^2_{J {\varLambda }} = k^2_J - {\varLambda }\), \(k^2_j\) being the standard Jeans wavenumber; consequently, larger critical values for wavelength, mass and a longer free fall time ensue straightaway to stabilize the medium against collapse formation.

Next we revisit the previous Jeans-Einstein instability criteria for the compressible resisto-elastic MHD (REMHD) scenario proposed in [15], through a non-ideal Ohm’s Law accounting for anomalous resistive features, due to the presence of a magnetic relaxation time.

These magnetic aspects lead to the introduction of a new significant parameter, representing the dimensionless magnetic relaxation time. Thus, in analogy with the Cattaneo number introduced to highlight analogous anomalous thermal effects rather than the diffusive Fourier law (see e.g. [16,17,18,19] and references therein), we here define the Ohm number denoted by \(O_{RF} ={\tau _j}/{t_c}\), where \(\tau _j\) is the magnetic relaxation time and \(t_c\) is the characteristic time equal to the ratio between the characteristic length L and the sound velocity \(c_s\).

The Ohm number plays a key role in the theory: it leads to a modified Reynolds number \(R_{ME} = R_M O_{RF}\) (\(R_M\) being the standard magnetic Reynolds number) for the valuation of the interlacements between the resistive and elastic behavior towards the micro-physics of cosmic plasmas.

In this way, the cosmological constant \({\varLambda }\), interpreted as a new critical wavenumber ascribable to DM, together with the classical plasma-beta (\(\beta \)), which is the ratio between the thermal and magnetic pressures, allow us to highlight an interesting relationship among the most significant parameters of our model (\({\varLambda } \), \(O_{RF}\), \(\beta \), \(R_{ME}\) and the gravitational constant G), against collapse formation. Thanks to such relationship in its limit cases, we may provide both a revisitation of the classical Jeans-Chandrasekhar instability criteria for the ideal frame, as in [3] and a generalization of more recent instability results, within the standard (non-ideal) resistive context [20].

Beyond the star formation process, this new scheme enhances the properties of (damped) torsional Alfvén waves, propagating with a larger effective wavespeed, due to the presence of the magnetic relaxation time, in triggering the intricate network of filamentary structures in every interstellar clouds [21], which in turn share common properties with the cosmic web in the Universe (we are thinking of clusters, filaments, sprouts, spicules, etc.). In fact the framework of the formation of such outstanding skeleton is such that of the gravitational instability.

Interestingly, a larger Alfvén wavespeed yields a smaller Alfvénic time scale compared with the resistive one. So, the ensuing time scalings seem to be more suitable to justify the rapid magnetic energy release on those processes likely to occur in regions of fast MR mechanisms (solar flares and spots, particle acceleration and coronal matter ejection) [22,23,24,25,26,27]. Again, the strategic players are the torsional oscillations in thin flux tubes, interpreted as the recently detected Alfvén waves [28], which, as conjectured in [15], are not affected either by compressional and gravitational aspects.

Also we might argue on the third sound effects. Obviously, the first sound concerns the acoustic waves and one may ascribe a second sound either to thermal waves, due to anomalous heat conduction theories, [29, 30] or to the superfluids sonic modes (as expected in neutron stars [8]). Indeed the presence of a third sound, may be conjectured just thinking of the recently detected noisy gravitational waves, not yet theoretically justified.

The paper is organized as follows. In Sect. 2 we introduce the MOG theory and test it to explain the stability and instability of an infinite homogeneous self-gravitating interstellar gas cloud modeled as a Euler–Darcy gas versus the standard Jeans criteria. In Sect. 3 we generalize the MOG theory to the presence of the magnetic field and propose a new non-ideal MHD (NIMHD) setting, characterized by anomalous resistive features accounting for the dimensionless magnetic relaxation time \(O_{RF}\). Again the general dispersion equation is derived and discussed to gain a special overview of the stability properties towards modified torsional Alfvén waves and longitudinal gravito-magneto-sonic waves against instability ones towards collapse formation. Also, we furnish a comparative study with other standard ideal and resistive MHD settings, enhancing some remarks on possible stability interlacements among the most significant parameters herein involved in Sect. 4. Then in Sect. 5 we read the response of the MOG effects on delayed magnetic properties, described by a Kelvin–Voigt type Ohm’s law. The last section is devoted to final comments.

2 The mathematical formulation of the modified gravitational instability problem: the Jeans–Einstein criteria vs collapse formation

Following the strategy for chemotactic mechanisms, see e.g. [31], but focusing overall on the accelerated expansion of the universe based on particles creation, we here address an evolutive reaction–diffusion equation for the gravitational potential \(\phi \) within the general balance law scheme [32], namely

$$\begin{aligned} \phi _t = - \nabla \cdot {\mathbf {J}}_\phi + r_\phi , \end{aligned}$$
(1)

where the gravitational inflow [33] is the pair of the gravity influx vector \({\mathbf {J}}_\phi \) and the gravity supply \(r_\phi \), defined by

$$\begin{aligned} \begin{aligned} {\mathbf {J}}_\phi = - D_\phi \nabla \phi \,, \qquad r_\phi&= - {\hat{{\varLambda }}} \phi + {\hat{h}} \rho \, . \end{aligned} \end{aligned}$$
(2)

The diffusive mobility \(D_\phi \), likewise the parameters \({\hat{{\varLambda }}}\) and \({\hat{h}}\), are here supposed to be positive constants (for easiness); in this way we are proposing a flux-gradient type constitutive setting for the gravitational flux with reaction terms. By substituting (2) into (1), we obtain the parabolic type gravitational equation

$$\begin{aligned} \tau _\phi \phi _t = {\varDelta } \phi - {\varLambda } \phi + h \rho , \end{aligned}$$
(3)

where \(\tau _\phi = {1}/{D_\phi }\), \({\varLambda } ={{\hat{{\varLambda }}}}/{D_\phi }\) and \( h = {{\hat{h}}}/{D_\phi }\).

It is worth to observe that \(\tau _\phi \) represents a gravitational relaxation time, while \({\varLambda }\) plays the role of the cosmological constant, which means that we are dealing within an expanding universe approach with \(h = 4 \pi G\), G being the gravitational constant. Equation (3) may be viewed as the non-stationary and non-relativistic Einstein’s equation with the cosmological constant.

Indeed, (3) in the limit \(\tau _\phi = 0\) and \({\varLambda } = 0\) reduces to the standard Poisson-type elliptic equation for \(\phi \), which in turn is not satisfied within the zeroth order approximation unless \(\rho = 0\) (Jeans swindle).

In order to test how our MOG theory affects the gravitational instability problem, we start with the simplest reasonable model of the Jeans problem, even so capable to capture the physics for a thorough qualitative understanding. Therefore we describe the infinite and homogeneous self-gravitating interstellar gas cloud as an adiabatic barotropic fluid, in the presence of the (non constant) self-gravitational force, within a Darcy type porous theory. Thus, the governing system reads as follows

$$\begin{aligned} \begin{aligned}&\rho _t + \nabla \rho \cdot {\mathbf {v}} + \rho \nabla \cdot {\mathbf {v}}= 0 \\&\rho \left[ {\mathbf {v}}_t+( \nabla {{\mathbf {v}}}){{\mathbf {v}}}\right] = - p^{\prime }(\rho ) \nabla \rho + \rho \nabla \phi - \rho \nu _D R{\mathbf {v}} \\&\tau _\phi \phi _t = {\varDelta } \phi - {\varLambda } \phi + h \rho \end{aligned} \end{aligned}$$
(4)

where \(\rho \) is the (positive) mass-density, \({\mathbf {v}}\) is the velocity and, neglecting thermal effects, we consider the general density dependent pressure term \(p = p(\rho )\), so that

$$\begin{aligned} \nabla p = p^{\prime }(\rho ) \nabla \rho = c_s^2(\rho ) \nabla \rho \, , \end{aligned}$$

where \(c_s^2(\rho )\) stands for the (squared) sound speed. In particular we address the polytropic equation of state \( p = A \rho ^\gamma \), \(A > 0\) and \(\displaystyle \gamma = {5}/{3}\) being the adiabatic index: in this frame \( c_s^2 = (\gamma {p})/{\rho }\), corresponds to the adiabatic sound speed.

The positive constant \(\nu _D\) represents the so called kinematical dynamic viscosity and, taking into account the Dupuit–Forchheimer relation and following the Nield’s strategy (see [34, 35]), we have introduced the retardability R, defined as \(R = {\epsilon }/{\kappa }\), \(\kappa \) (\(\kappa > 0\)) and \(\epsilon \) (\(0< \epsilon < 1\)) being the permeability and the porosity, respectively. The presence of a Darcy-type porous term in the momentum Eq. (4)\(_2\) may be easily compared with the friction force introduced in [31] within the different context of chemically driven hydrodynamic motions.

To begin we prerform the standard perturbation approach, as in most of the works on the Jeans instability. Therefore, as a starting point, we consider the stationary and homogeneous quiescent state, described by the constant variables

$$\begin{aligned} {\mathbf {v}}_0 = {\mathbf {0}} \,, \quad \rho = \rho _0 \,, \quad p = p_0 = p(\rho _0) \,, \quad \phi = \phi _0 = \frac{h \rho _0}{{\varLambda }} \, , \end{aligned}$$

so that, once we take into account the expansion of the Universe as currently experienced, we no longer need to invoke the so called Jeans swindle ([13, 36]), due to the stationary approximation of (4)\(_3\) in the limit \({\varLambda } \rightarrow 0^+\).

Then we linearize (4) around this equilibrium state, by introducing the perturbation state \(\delta ^*s = (\delta ^*{\mathbf {v}}, \delta ^*\rho , \delta ^*p, \delta ^*\phi )\). Taking into account that, for small perturbations, the pressure and density perturbations are related through the adiabatic sound speed, i.e. \(\nabla \delta ^*p = c_s^2 \nabla \delta ^*\rho \) with \(c_s^2 =( \gamma {p_0})/{\rho _0}\), the linearized perturbation equations are

$$\begin{aligned} \begin{aligned}&\delta ^*\rho _t + \rho _0 \nabla \cdot \delta ^*{\mathbf {v}}=0 \\&\rho _0\delta ^*{\mathbf {v}}_t= - c_s^2 \nabla \delta ^*\rho + \rho _0\nabla \delta ^*\phi - \rho _0 \nu _D R \delta ^*{\mathbf {v}} \\&\tau _\phi \delta ^* \phi _t = {\varDelta } \delta ^*\phi - {\varLambda } \delta ^*\phi + h \delta ^*\rho \end{aligned} \end{aligned}$$
(5)

To carry out the linear stability analysis, we now introduce a spectral tool in order to investigate the temporal behavior of these perturbations. So, in spatially unbounded domains, we employ the Fourier modes technique and look for perturbations of the dispersive wave form

$$\begin{aligned} \delta ^* s = s_1 e^{i({\mathbf {k}}\cdot {\mathbf {x}}-\omega t)} \,, \end{aligned}$$

where \(s_1\ne 0\) is the constant amplitude, \({\mathbf {k}}\) is the (real) wavenumber vector, with \(|{\mathbf {k}}|^2>0\), and \(\omega \) denotes the frequency of the disturbance; for future convenience we also introduce the growth rate parameter \(\sigma = - i \omega \), which is generally a complex number. Therefore stable modes are those with \(Re(\sigma ) \le 0\) and unstable modes occur when \(Re(\sigma ) > 0\). Among the unstable modes, we may expect monotonically growing unstable ones for \(Im(\sigma ) = 0\), whereas, when \(Im(\sigma ) \ne 0\), overstable modes arise, so that the perturbations oscillate with growing amplitudes. Likewise for stable modes, the model can experience either monotonically decreasing modes (\(Im(\sigma ) = 0\)) or exponentially decreasing oscillatory modes (\(Im(\sigma ) \ne 0\)), accounting for damped perturbations, due to the dissipative mechanisms within the model.

Hence, by using classical identities, we find the following dispersion system for the perturbation amplitudes

$$\begin{aligned} \begin{aligned}&\omega \rho _1 - \rho _0 k v_{1 n}= 0 \\&(\omega + i \nu _D R) v_{1 n} = c_s^2 \rho _1 \frac{k}{\rho _0} - \phi _1 k \\&(\omega + i \nu _D R) v_{1 1} = 0 \\&(\omega + i \nu _D R) v_{1 2}= 0 \\&\left[ \tau _\phi \omega + (k^2 + {\varLambda })i \right] \phi _1 - hi \rho _1 = 0 \end{aligned} \end{aligned}$$
(6)

where we have used the standard decomposition form

$$\begin{aligned} \begin{aligned} {{\mathbf {v}}}_1 = ({{\mathbf {v}}}_1 \cdot {{\mathbf {n}}}) {{\mathbf {n}}} + ({{\mathbf {v}}}_1 \cdot {{\mathbf {t}}}_1) {{\mathbf {t}}} _1 + ({{\mathbf {v}}}_1 \cdot {{\mathbf {t}}}_2) {{\mathbf {t}}} _2 = v_{1 n}{{\mathbf {n}}} + v_{1 1} {{\mathbf {t}}} _1 + v_{1 2} {{\mathbf {t}}} _2\,, \end{aligned} \end{aligned}$$

with \( {\mathbf {n}} = {\mathbf {k}}/{k}\), \(k = |{\mathbf {k}}|\), and \( {\mathbf {t}}_i\) (\(i = 1, 2\)) being two orthonormal transverse vectors, such that \({{\mathbf {t}}_i }\cdot {{\mathbf {n}}}= 0\) and \({{\mathbf {t}}}_i \cdot {{\mathbf {t}}}_j = \delta _{ij}\).

From now on we will suppose

$$\begin{aligned} \omega \notin \left\{ 0, - i \nu _DR, - \displaystyle \frac{(k^2 + {\varLambda })i}{\tau _\phi }\right\} \,, \end{aligned}$$

thus avoiding stationary and stable modes, uninteresting for our purposes. The discussion can be simplified as follows.

From (6)\(_3\)–(6)\(_4\) we have \(v_{1 i}= 0, \; i = 1, 2\) and hence only longitudinal compressible dispersive modes are allowable; Eq. (6)\(_5\) gives us the relation between the amplitudes \(\phi _1\) and \(\rho _1\), whereas from (6)\(_1\) we find the link between \(\rho _1\) and \(v_{1 n}\). Finally, by gathering all the previous results into (6)\(_2\), we derive the dispersion equation

$$\begin{aligned} \omega ^2 - c_s^2 k^2 + \nu _D R i \omega + \frac{h \rho _0 i k^2}{\tau _\phi \omega + \left( k^2 + {\varLambda }\right) i} = 0, \end{aligned}$$

which in terms of \(\sigma \) reads

$$\begin{aligned} \begin{aligned}&\tau _\phi \sigma ^3 + \left( k^2 + {\varLambda } + \nu _DR \tau _\phi \right) \sigma ^2 + \left[ \nu _D R(k^2 + {\varLambda }) + c_s^2 \tau _\phi k^2 \right] \sigma \\&\quad + c_s^2 k^2 \left( k^2 + {\varLambda } - k_j^2\right) = 0 \, , \end{aligned} \end{aligned}$$
(7)

where \(k_j^2 = (4 \pi G {\rho _0})/{c_s^2}\) is the standard Jeans critical wavenumber.

By applying the Routh-Hurwitz criterion, the three solutions of (7) have negative real parts, namely the joint presence of \(\tau _\phi , {\varLambda }\) and \(\nu _D R\) allows to recover three (asymptotically) stable modes if and only if the following condition hold

$$\begin{aligned} k^2 + {\varLambda } - k_j^2 > 0 \, . \end{aligned}$$

It is noteworthy that the same condition should be required in the limit case \(\tau _\phi \rightarrow 0\) and \(\nu _D \rightarrow 0\), corresponding to the (Helmholtz-type) stationary Einstein gravitational equation, and neglecting the Darcy term.

The constant \({\varLambda }\) is therefore strategic for recovering new stability properties of the medium: in particular if \({\varLambda } > k_j^2\), namely \(c_s^2 {\varLambda } > 4 \pi G \rho _0\) (which seems however difficult to be experimentally justified), we would always find (asymptotically stable) damped gravito-sonic modes independently of \(k^2\) vs collapse formation. However, when \({\varLambda } < k_j^2\) we may define a lower critical threshold,

$$\begin{aligned} \displaystyle k_{J{\varLambda }}^2 = k_j^2 - {\varLambda } < k_j^2 \end{aligned}$$

such that, for longitudinal compressible modes, the modified Jeans instability criterion reads as follows \(k^2 < \displaystyle k_{J{\varLambda }}^2\). In this way, with respect to the classical Jeans criterion, we have gained stability also in the range \(\displaystyle k_{J{\varLambda }}^2 \le k^2 \le k_j^2\).

Dimensionally speaking, \({\varLambda }\) may be interpreted as a critical wavenumber due to DM, i.e. \({\varLambda } = k^2_{DM}\). Therefore from now on we will refer to \(\displaystyle k_{J{\varLambda }}^2\) as the Jeans-Einstein critical wavenumber, which is typical of our theory and leads to a delayed collapse formation vs damped gravito-sonic waves (for \(k^2 \ge k_{J{\varLambda }}^2\)). Consequently, the new critical wavelength and mass above which the gravitational field of the fluctuations is capable to create collapse formation are

$$\begin{aligned}\begin{aligned} \lambda _{J {\varLambda }}&= \frac{2 \pi }{k_{J {\varLambda }}} \, , \\ M_{J {\varLambda }}&= \frac{4\pi \rho _0}{3} \left( \frac{\lambda _{J {\varLambda }}}{2} \right) ^3 = \frac{4\pi \rho _0}{3} \left( \frac{\pi c_s}{\sqrt{4 \pi G \rho _0}} \right) ^3 \left( 1 - \displaystyle \frac{{\varLambda } c_s^2}{4 \pi G \rho _0} \right) ^{ - \frac{3}{2}} \end{aligned} \end{aligned}$$

respectively. Obviously, these larger thresholds coincide with the classical ones within the traditional Jeans treatment, whenever \({\varLambda } \rightarrow 0^+\).

To sum up, \({\varLambda }\) represents an important stabilizing driver against the onset of unstable gravito-sonic modes, but the instability results are not affected by the gravitational relaxation time \(\tau _\phi \).

3 The Oldroyd-B type Ohm’s Law for resisto-elastic non-ideal MHD plasmas within the MOG theory

Driven by the unquestionable importance of magnetic effects in gravitational instability analysis, in this section we present a special non-ideal MHD background for a self-gravitating interstellar plasma, also preserving the Darcy porous term, here interpreted as a damping friction force.

We are interested in capturing micro- to macro- scales magnetic properties through a generalized form of the resistive Ohm’s law. To this aim we write down the resistive Ohm’s law within the one fluid approach [37]:

$$\begin{aligned} {\mathbf {E}}+ {\mathbf {v}} \times {\mathbf {B}}={\mathbf {R}}\, , \end{aligned}$$

where \({\mathbf {E}}\), \({\mathbf {B}}\) and \({\mathbf {R}}\) denote the electric field, the magnetic field and the resistivity vector respectively. In analogy with other thermal and visco-elastic settings [35], the new resisto-elastic network model is constructed by the additional split of the overall resistive response, due to non-idealness:

$$\begin{aligned} {\mathbf {R}}={\mathbf {R}}_r+ {\mathbf {R}}_e\,, \end{aligned}$$

where \({\mathbf {R}}_r= \eta _r \nabla \times {\mathbf {B}}\) (\(\eta _r > 0\) and constant) represents the classical (parabolic) resistive part, whereas the extra elastic resistivity vector \({\mathbf {R}}_e\) is supposed to satisfy the rate-type equation

$$\begin{aligned} \tau _j\left[ {\displaystyle {\mathbf {R}}_{e_t}}+( \nabla {{\mathbf {R}}_e}){{\mathbf {v}}} - \frac{1}{2}(\nabla \times {{\mathbf {v}}})\times {{\mathbf {R}}_e} \right] +{\mathbf {R}}_e=\eta _e\nabla \times {\mathbf {B}} \, . \end{aligned}$$
(8)

Here \(\tau _j\) denotes the (positive and constant) magnetic relaxation time and \(\eta _e\) is a positive constant.

The total magnetic resistivity, such that the electric conductivity is just the inverse of \(\eta \mu _0\), is now given by \(\eta = \eta _r + \eta _e\).

Note that the \((\nabla \times {{\mathbf {v}}})\)-term is due to the introduction of the objective Jaumann time derivative in place of the convected one. It is worth to stress that this non-linear invariant term does not affect the linear stability analysis we are going to perform, but it may have an important role within the non linear wave propagation tool.

Two remarks are at order:

  • We now have \(\nabla \cdot {\mathbf {R}}\ne 0\) and beyond the diffusion time, typical of the classical resistive scheme, \(\tau _d = {\mu _0 L^2}/{\eta }\) (L being a characteristic length scale), this theory exhibits the new small relaxation time \(\tau _j\) accounting for the lagging response of the current density vector on the magnetic field;

  • Looking at (8), one captures another salient constitutive feature: the ratio \(\eta _e /\tau _j\) will represent a further (transverse) wavespeed to be added to the classical Alfvén wavespeed and hence new smaller Alfvén crossing times to be compared with \(\tau _d\) and \(\tau _j\) ensue.

The Ohm’s law describing these anomalous resistive features may be referred to as the Oldroyd-B type Ohm’s law in its split form.

Interestingly, when \(\eta _r=0\), the theory reduces to a novel hyperbolic non-ideal MHD model, leading to elastic memory-type behaviors through the presence of the collisional relaxation time \(\tau _j\). Anyway the ideal frozen-in condition is broken. The non-ideal behavior becomes very important when the magnetic field gradient is very large with a steepening effect, that is to say in the presence of a strong current density. This exceptional behavior just happens in micro-scale structures, such as current sheets arising from zero magnetic points, where the magnetic field lines exhibit a 3D X-type topology [37, 38]. Moreover. when small-scale fluctuations are incorporated, macroscopic zero magnetic points in the (terrestrial) magnetosphere become a network of microscopic nulls in micro-scale diffusion regimes, thus enhancing the intertwining between micro- to macro- complex structures

With greater reason, when both (parabolic and hyperbolic) effects act jointly, such a mathematical model seems to be more adequate to justify the fast MR mechanisms in the solar atmosphere.

The relevant equations governing these MHD processes are the continuity equation, the momentum equation including the gravitational and Lorentz forces together with the Darcy term and the induction equation incorporating the anomalous resistivity of Oldroyd-B type [15], that is

$$\begin{aligned} \begin{aligned}&\rho _t + \nabla \rho \cdot {\mathbf {v}} + \rho \nabla \cdot {\mathbf {v}}= 0 \\&\rho \left[ {{\mathbf {v}}_t}+( \nabla {{\mathbf {v}}}){{\mathbf {v}}}\right] = - p^\prime (\rho )\nabla \rho - \rho \nu _D R {{\mathbf {v}}} +\rho \nabla \phi +\frac{1}{\mu _0} (\nabla \times {\mathbf {B}}) \times {\mathbf {B}} \\&{{\mathbf {B}}}_t+(\nabla {\mathbf {B}}) {\mathbf {v}}-(\nabla {\mathbf {v}} ) {\mathbf {B}}+(\nabla \cdot {\mathbf {v}}) {\mathbf {B}} =-\eta _r \nabla \times \nabla \times {\mathbf {B}}- \nabla \times {\mathbf {R}}_e \\&\tau _j\left[ {\displaystyle {\mathbf {R}}_{e_t}}+( \nabla {{\mathbf {R}}_e}){{\mathbf {v}}} - \frac{1}{2}(\nabla \times {{\mathbf {v}}})\times {{\mathbf {R}}_e} \right] +{\mathbf {R}}_e=\eta _e\nabla \times {\mathbf {B}} \end{aligned} \end{aligned}$$
(9)

together with the reaction-diffusion gravitational Eq. (3) and supplemented by the usual constraint \(\nabla \cdot {{\mathbf {B}}}= 0\).

3.1 Linear stability analysis

Some aspects of the first order linear stability analysis of an infinite homogeneous self-gravitating plasma, permeated by a uniform magnetic field \({{\mathbf {B}}}_0\), are now to be reexamined in the light of the new constitutive modeling.

Thus we address a small perturbation around the following homogeneous and stationary (quiescent) state

$$\begin{aligned} \rho _0 > 0, \quad {{\mathbf {v}}}_0 = {{\mathbf {0}}}, \quad \displaystyle \phi _0 = \frac{h \rho _0}{{\varLambda }}, \quad {{\mathbf {B}}}_0, \quad {{\mathbf {R}}}_{e_0} = {{\mathbf {0}}} \end{aligned}$$

Following step by step the previous perturbation technique, we arrive at the linearized perturbation system

$$\begin{aligned} \begin{aligned}&\delta ^* \rho _t + \rho _0 \nabla \cdot \delta ^* {{\mathbf {v}}} = 0 \\&\rho _0 \delta ^*{{\mathbf {v}}_t} = - c_s^2 \nabla \delta ^* \rho + \rho _0 \nabla \delta ^* \phi - \rho _0 \nu _D R \delta ^*{{\mathbf {v}}} +\frac{1}{\mu _0} (\nabla \times \delta ^*{{\mathbf {B}}} )\times {{\mathbf {B}}}_0 \\&\nabla \cdot \delta ^*{{\mathbf {B}}}= 0 \\&\delta ^*{{\mathbf {B}}}_t - ( \nabla \delta ^* {{\mathbf {v}}}) {{\mathbf {B}}}_0+ {{\mathbf {B}}}_0 \nabla \cdot \delta ^*{{\mathbf {v}}} = - \eta _r \nabla \times \nabla \times \delta ^* {{\mathbf {B}}} - \nabla \times \delta ^* {{\mathbf {R}}}_e \\&\tau _j \delta ^*{\displaystyle {\mathbf {R}}_{e_t}}+ \delta ^*{{\mathbf {R}}_e} = \eta _e\nabla \times \delta ^*{\mathbf {B}} \\&\tau _\phi \delta ^*\phi _t = {\varDelta } \delta ^* \phi - {\varLambda } \delta ^* \phi + h \delta ^* \phi \end{aligned} \end{aligned}$$

Again, the formal application of the Fourier normal modes tool leads to the following dispersive system for the perturbation amplitude \(s_1 = ( \rho _1, {{\mathbf {v}}}_1, \phi _1, {\mathbf {B}}_1, {{\mathbf {R}}}_1)\)

$$\begin{aligned} \begin{aligned}&\omega \rho _1 - \rho _0 k v_{1 n} = 0 \\&(\omega + i \nu _D R) v_{1 n} = c_s^2 \rho _1 \frac{k}{\rho _0} - \phi _1 k + \displaystyle \frac{k}{\mu _0 \rho _0}( {{\mathbf {B}}}_0 \cdot {{\mathbf {B}}}_1) \\&(\omega + i \nu _D R) v_{1 1} = - \displaystyle \frac{k}{\mu _0 \rho _0} B_{0 n} B_{1 1} \\&(\omega + i \nu _D R) v_{1 2} = - \displaystyle \frac{k}{\mu _0 \rho _0} B_{0 n} B_{1 2} \\&(\omega + \eta _r k^2 i) B_{1 1} + k B_{0 n} v_{1 1} - k B_{0 1} v_{1 n} = k {{\mathbf {n}}} \times {{\mathbf {R}}}_1 \cdot {{\mathbf {t}}}_1 \\&(\omega + \eta _r k^2 i) B_{1 2} \ + k B_{0 n} v_{1 2} - k B_{0 2} v_{1 n} = k {{\mathbf {n}}} \times {{\mathbf {R}}}_1 \cdot {{\mathbf {t}}}_2 \\&(\tau _j \omega + i) R_{1 n} = 0 \\&(\tau _j \omega + i) R_{1 1} = - \eta _e k {{\mathbf {n}}} \times {{\mathbf {B}}}_1 \cdot {{\mathbf {t}}}_1 \\&(\tau _j \omega + i) R_{1 2} = - \eta _e k {{\mathbf {n}}} \times {{\mathbf {B}}}_1 \cdot {{\mathbf {t}}}_2 \\&\left( \tau _\phi \omega + (k^2 + {\varLambda })i \right) \phi _1 = hi \rho _1 \end{aligned} \end{aligned}$$
(10)

supplemented by the condition that the magnetic amplitude is always orthogonal to \({\mathbf {n}}\), i.e. \( {{\mathbf {B}}}_{1} \cdot {{\mathbf {n}}} = 0\). As before, for easiness in writing, in (10) we have used the following notations

$$\begin{aligned} B_{0 n} = {{\mathbf {B}}}_{0} \cdot {{\mathbf {n}}} \, , \quad B_{0 i} = {{\mathbf {B}}}_{0} \cdot {{\mathbf {t}}}_i \, , \quad B_{1 i } = {{\mathbf {B}}}_{1} \cdot {{\mathbf {t}}}_i \, , \quad i= 1, 2 \, , \end{aligned}$$

for normal and tangential components;

In order to simplify the picture, henceforth we suppose \(\omega \ne 0\) (stationary case) and also

$$\begin{aligned} \omega \notin \left\{ - i \nu _D R, - i \eta _r k^2, - \displaystyle \frac{i}{\tau _j} \right\} \,, \end{aligned}$$

corresponding to unphysical stable cases.

Therefore, from (10)\(_7\), \(R_{1 n} = 0 \), namely the theory allows only for non zero tangential components of \({{\mathbf {R}}}_1\).

Equations (10)\(_3\) and (10)\(_4\) give us the relations between the transverse components of \({{\mathbf {v}}}_1\) and \({{\mathbf {B}}}_1\),

$$\begin{aligned} v_{1 i} = - \displaystyle k \frac{ B_{0 n} B_{1 i}}{\mu _0 \rho _0 (\omega + i \nu _D R)} \, , \quad i = 1,2 \end{aligned}$$
(11)

and we can solve Eqs. (10)\(_8\) and (10)\(_9\) for \(R_{1 i}\) in terms of \(B_{1 i}\),

$$\begin{aligned} R_{1 1} = \frac{\eta _e k }{(\tau _j \omega + i)}B_{1 2} \, , \quad R_{1 2} = - \displaystyle \frac{\eta _e k }{(\tau _j \omega + i)}B_{1 1} \, . \end{aligned}$$

As expected, the new theory does not affect the relation between \(\rho _1\) and \(\phi _1\):

$$\begin{aligned} \phi _1 = \frac{hi \rho _1}{\tau _\phi \omega + (k^2 + {\varLambda })i } = \frac{hi k \rho _0}{\omega \left( \tau _\phi \omega + (k^2 + {\varLambda })i \right) }v_{1 n} \, . \end{aligned}$$

By gathering all these results into (10)\(_2\), (10)\(_5\) and (10)\(_6\), we easily reduce to the following Cramer system in \({\mathbb {R}}^3\) for the essential amplitudes \(v_{1 n}\), \(B_{1 1}\) and \(B_{1 2 }\)

$$\begin{aligned} \begin{aligned}&\beta v_{1 n} - \frac{k \omega }{\mu _0 \rho _0} \left[ \omega + \frac{i}{\tau _\phi }(k^2 + {\varLambda }) \right] (B_{0 1}B_{1 1} + B_{0 2} B_{1 2}) = 0 \\&\alpha B_{1 1} - B_{0 1} k \left[ \omega ^2 + i \omega \left( \frac{1}{\tau _j} + \nu _D R \right) - \frac{\nu _D R}{\tau _j} \right] v_{1 n} = 0 \\&\alpha B_{1 2} - B_{0 2} k \left[ \omega ^2 + i \omega \left( \frac{1}{\tau _j} + \nu _D R \right) - \frac{\nu _D R}{\tau _j} \right] v_{1 n} = 0 \end{aligned} \end{aligned}$$

with

$$\begin{aligned} \begin{aligned} \beta&= \;\omega ^3 + i \omega ^2 \left( \frac{k^2 + {\varLambda }}{\tau _\phi } + \nu _D R \right) - \omega \left[ c_s^2 k^2 + \frac{\nu _D R}{\tau _\phi } \left( k^2 + {\varLambda }\right) \right] \\&\quad -i \frac{k^2 c_s^2}{\tau _\phi }(k^2 + {\varLambda } - k^2_j) , \\ \alpha&= \;\omega ^3 + i \omega ^2 \left( \frac{1}{\tau _j} + \eta _r k^2 + \nu _D R \right) \\&\quad - \omega \left( c_{ARE}^2 k^2 + \nu _D R \eta _r k^2 + \frac{\nu _D R}{\tau _j} \right) \\&\quad - \frac{i}{\tau _j}k^2 (c_A^2 + \eta _e \nu _D R), \end{aligned} \end{aligned}$$

where we have defined the modified larger (squared) Alfvén wavespeed due to the novel elastic effects

$$\begin{aligned} c_{ARE}^2 = \frac{B_{0 n}^2}{\mu _0 \rho _0} + \frac{\eta _r + \eta _e}{\tau _j} \end{aligned}$$

in place of the classical Alfvén wavespeed \(c_A^2 = {B_{0 n}^2}/{\mu _0 \rho _0}\).

By a quick inspection, we obtain the general dispersion equation

$$\begin{aligned} \alpha \left\{ \alpha \beta -\frac{ k^2 c_B^2\omega \left[ \tau _\phi \omega + i(k^2 + {\varLambda }) \right] \left[ \tau _j\omega ^2 +i \omega (1 +{\tau _j} \nu _D R)- \nu _D R \right] }{\tau _\phi \tau _j} \right\} =0, \end{aligned}$$
(12)

where we have used the classical notation

$$\begin{aligned} c_B^2 = \frac{B_{0 1}^2 + B_{0 2}^2}{\mu _0 \rho _0} \, . \end{aligned}$$

Equation (12) can be split in

$$\begin{aligned} \alpha = 0 \end{aligned}$$

and

$$\begin{aligned} \begin{aligned} \alpha \beta - k^2 c_B^2\omega \left[ \omega + \frac{i(k^2 + {\varLambda })}{\tau _\phi } \right] \left[ \omega ^2 + i \omega \left( \frac{1}{\tau _j} + \nu _D R\right) - \frac{\nu _D R}{\tau _j} \right] = 0. \end{aligned} \end{aligned}$$

3.1.1 Damped non-gravitational torsional Alfvén waves

The first cubic equation \(\alpha = 0\), rewritten in terms of the growth rate parameter \(\sigma \) as

$$\begin{aligned} \begin{aligned}&\sigma ^3 +\left( \frac{1}{\tau _j} + \eta _r + \nu _D R \right) \sigma ^2 + \left[ \left( c_{ARE}^2 + \nu _D R \eta _r\right) k^2 + \frac{\nu _D R}{ \tau _j} \right] \sigma \\&\quad + \frac{k^2}{\tau _j} \left( c_A^2 + \eta _e \nu _D R\right) = 0 \, , \end{aligned} \end{aligned}$$
(13)

provides three novel non-gravitational torsional shear-like Alfvén waves, which are damped due to the simultaneous presence of \(\eta _r\), \(\nu _D R\) and \(\tau _j\).

To sum up, \(\alpha = 0\) allows to show how the non-ideal resisto-elastic Ohm’s law may affect the behavior of the compressible transverse magnetic (dispersive) waves. The application of the Routh–Hurwitz test always guarantees the asymptotic stability properties of the ensuing resisto-elastic magnetic modes, unaffected by (modified or not) gravity. As expected, if we neglect the \(\nu _D R\)-term, the equation \(\alpha = 0\) coincides with the Eq. (12) in [15] governing incompressible transverse modes. By all means the pure torsional behavior of the ideal Alfvén waves is lost.

3.1.2 Longitudinal gravito-magneto-sonic waves versus the delayed collapse formation

Obviously, the gravitational instability analysis is quite concentrated on the sixth grade dispersion equation, which may be easily rewritten in terms of \(\sigma \) as follows

$$\begin{aligned} \sigma ^6 + a_1 \sigma ^5 + a_2 \sigma ^4 + a_3 \sigma ^3 + a_4 \sigma ^2 + a_5 \sigma + a_6 = 0 \, , \end{aligned}$$
(14)

with

$$\begin{aligned} \begin{aligned} a_1&= \frac{k^2 + {\varLambda }}{\tau _\phi } + \frac{1}{\tau _j}+ \eta _r k^2 + 2 \nu _D R \\ a_2&= \left( c_s^2 + c_{ABRE}^2 + \frac{\nu _D R}{\tau _\phi } + \nu _D R \eta _r \right) k^2 + \frac{{\varLambda } \nu _D R}{\tau _\phi } + \frac{\nu _D R}{\tau _j} \\&\quad + \left( \frac{k^2 + {\varLambda }}{\tau _\phi } + \nu _D R \right) \left( \frac{1}{\tau _j} + \eta _r k^2 + \nu _D R \right) \\ a_3&= \left( c_s^2 + c_{AB}^2 + \frac{\nu _D R}{\tau _\phi } + \nu _D R \eta _e \right) \frac{k^2}{\tau _j} +\frac{(\nu _D R )^2}{\tau _j} \\&\quad + \frac{k^2}{\tau _\phi } \left( c_s^2 + c_{ABRE}^2 + \nu _D R \eta _r \right) (k^2 + {\varLambda } - k_{JRED}^2) \\&\quad + \eta _r k^4 \left( c_s^2 + \frac{\nu _D R}{\tau _\phi } \right) + \nu _D R k^2 \left( c_s^2 + c_{ABRE}^2 + \nu _D R \eta _r + \frac{\nu _D R}{\tau _\phi } \right) \\&\quad + \frac{{\varLambda } \nu _D R}{\tau _\phi \tau _j} + \frac{\nu _D R}{\tau _\phi } \left( \frac{k^2 + {\varLambda }}{\tau _j} + {\varLambda } \nu _D R \right) \\ a_4&= \frac{k^2 c_s^2}{\tau _\phi }(k^2 + {\varLambda } - k_j^2) \left( \frac{1}{\tau _j} + \eta _r k^2 + \nu _D R \right) \\&\quad + \frac{k^2}{\tau _j} (c_A^2 + \eta _e \nu _D R)\left( \frac{k^2 + {\varLambda }}{\tau _\phi } + \nu _DR \right) \\&\quad + c_B^2 k^2 \left( \frac{\nu _D R}{\tau _j} + \frac{k^2 + {\varLambda }}{\tau _\phi } \left( \frac{1}{\tau _j} + \nu _D R \right) \right) \\&\quad + \frac{\nu _DR}{\tau _j} \left( c_s^2k^2+ \nu _DR \frac{k^2 + {\varLambda }}{\tau _\phi }\right) \\&\quad + (c_{ARE}^2 + \nu _D R \eta _r) k^2 \left( \left( c_s^2 + \frac{\nu _D R}{\tau _\phi }\right) k^2 + \frac{ {\varLambda } \nu _D R}{\tau _\phi } \right) \\ a_5&= c_B^2 k^2 \frac{\nu _D R}{\tau _j \tau _\phi } (k^2 + {\varLambda })\\&\quad + \frac{k^2}{\tau _\phi }\left[ \left( c_{ARE}^2 + \nu _D R \eta _r\right) k^2 + \frac{\nu _D R}{\tau _j} \right] c_s^2 (k^2 + {\varLambda } - k_j^2) \\&\quad + \frac{k^2}{\tau _j}\left( c_A^2 + \eta _e \nu _D R\right) \left[ \left( c_s^2 + \frac{\nu _D R}{\tau _\phi } \right) k^2+ \frac{{\varLambda } \nu _D R}{\tau _\phi } \right] \\ a_6&= \frac{c_s^2 k^4}{\tau _\phi \tau _j}\left( c_A^2 + \eta _e \nu _D R \right) \left( k^2 + {\varLambda } - k_j^2\right) \end{aligned} \end{aligned}$$

where we have used the notations

$$\begin{aligned} c_{AB}^2 = c_A^2 + c_B^2 \, , \quad c_{ABRE}^2 = c_{ARE}^2 + c_B^2 \, , \end{aligned}$$

and

$$\begin{aligned} k^2_{JRED} = \frac{4 \pi G \rho _0}{c_s^2 + c_{ABRE}^2 + \eta _r \nu _D R} \end{aligned}$$
(15)

to denote the modified lower (with respect to \(k_j^2\)) critical wavenumber, due to the inclusion of all the new effects herein.

In the following we will focus on the analysis of the solutions of the Eqs. (13)–(14) and discuss in detail two significant constitutive subsettings.

3.2 The purely parabolic NIMHD setting within the stationary MOG theory

Within the limit case \(\tau _\phi \rightarrow 0\), \(\tau _j \rightarrow 0\) and \(\nu _D \rightarrow 0\), corresponding to a self-gravitating resistive homogeneous plasma on the expanding universe approach, the dispersion Eq. (14) reduces to the following fourth grade equation

$$\begin{aligned} \begin{aligned}&\sigma ^4 + \eta k^2 \sigma ^3 + \frac{k^2}{k^2 + {\varLambda }}(c_{AB}^2 + c_s^2)(k^2 - k_{JAB{\varLambda }}^2)\sigma ^2 \\&\quad + \eta c_s^2\frac{ k^4}{k^2 + {\varLambda }}(k^2 - k_{j{\varLambda }}^2) \sigma + c_A^2 c_s^2 \frac{k^4}{k^2 + {\varLambda }}(k^2 - k_{j{\varLambda }}^2) = 0 \, , \end{aligned} \end{aligned}$$
(16)

where

$$\begin{aligned} \begin{aligned} k^2_{JAB{\varLambda }} = \frac{4 \pi G \rho _0}{c_s^2 + c_{AB}^2 } - {\varLambda } < k_{j{\varLambda }}^2. \end{aligned} \end{aligned}$$

The stability analysis of (16) generalizes to the expanding universe approach the analogous discussion of Eq. (18) in [20], where in particular the zero-resistivity case (corresponding to the ideal MHD limit) has been thoroughly investigated, both analytically and numerically.

Our primary aim is to stress how the parameters \(\eta \) and \({\varLambda }\) act jointly against collapse formation. It is well known that the Jeans criterion remains unaffected by ideal magnetic effects, except for those modes propagating perpendicular to the direction of the uniform magnetic field, namely \(c_A = 0\) (see e.g. [3]). Now, in this case, besides the marginally stable mode \(\sigma = 0\), we find three longitudinal gravitational modes satisfying

$$\begin{aligned} \begin{aligned}&\sigma ^3 + \eta k^2 \sigma ^2 + \frac{k^2}{k^2 + {\varLambda }}(c_B^2 + c_s^2)(k^2 - k^2_{JB{\varLambda }}) \sigma \\&\quad + \frac{\eta c_s^2 k^4}{k^2 + {\varLambda }}(k^2 - k^2_{J {\varLambda }}) = 0 \end{aligned} \end{aligned}$$

where the notation

$$\begin{aligned} k^2_{JB{\varLambda }} = k^2_{JB} - {\varLambda } < k^2_{J {\varLambda }} \end{aligned}$$

stands for a new lower critical threshold, due to the effective wavespeed, now reduced to \((c_s^2 + c_B^2)\).

Again, the system under study will be (asymptotically) stable if and only if

$$\begin{aligned} k^2 \ge k^2_{J{\varLambda }} \end{aligned}$$

since the coefficients are all positive and also the third condition of the Routh-Hurwitz test is trivially satisfied because it requires the positiveness of

$$\begin{aligned} c_B^2 (k^2 - k^2_{JB{\varLambda }}) + c_s^2(k^2_{J{\varLambda }} - k^2_{JB{\varLambda }}) \, . \end{aligned}$$

Hence, differently from the ideal MHD setting, the presence of resistive magnetic effects will not decrease the instability in the direction orthogonal to \({\mathbf {B}}_0\), except for the modes \(k^2 = k^2_{J {\varLambda }}\). In fact in this case, besides the marginally stable mode \(\sigma = 0\), we get a second grade dispersion equation whose solutions always have negative real parts.

On the other hand, when \(c_A^2 > 0\) (namely for modes parallel to \({\mathbf {B}}_0\)), the discussion of (16) accounts for a larger effective speed, just due to the presence of \(c_{AB}^2\). We know that the Routh-Hurwitz test asserts that the four solutions to (16) have negative real parts, corresponding to four damped oscillatory modes, if and only if the following requirements

$$\begin{aligned} \begin{aligned}&\eta k^4 \left[ c_{AB}^2 (k^2 - k_{JAB{\varLambda }}^2) + c_s^2 (k^2_{J {\varLambda }} - k_{JAB{\varLambda }}^2) \right]> 0 \, , \\&\frac{\eta ^2 k^8}{k^2 + {\varLambda }} c_s^2c_B^2(k^2 - k^2_{J {\varLambda }}) > 0 \end{aligned} \end{aligned}$$

hold, together with the positiveness of all the coefficients.

Indeed, the condition for stable gravito-magneto-sonic modes reads as \(k^2 > k^2_{J {\varLambda }}\), since \((k^2_{J {\varLambda }} - k_{JAB{\varLambda }}^2) > 0\) is always true.

Also, as before, it is a simple matter to show that in the limit case \(k^2 = k^2_{J {\varLambda }}\) we find \(\sigma = 0\) as a (double) marginally stable mode. By all means, the sufficient condition for gravitational collapse formation becomes

$$\begin{aligned} k^2 < k^2_{JAB{\varLambda }} \, , \end{aligned}$$

thus leading to a new larger critical wavelength

$$\begin{aligned} \lambda _{JAB{\varLambda }} = \frac{2\pi }{ k_{JAB{\varLambda }} }\,. \end{aligned}$$

As a consequence, larger critical mass and free fall time follow. Obviously the purely oscillatory modes of the standard Jeans model are lost; also, the instability arises in the form of overstable oscillations, versus the classical exponential growth.

As regards the behavior of the non-gravitational and non-compressible transverse Alfvèn modes, characterized by [see e.g. (13)]

$$\begin{aligned} \sigma ^2 + \eta k^2 \sigma + c_A^2 k^2 = 0 \end{aligned}$$

we point out that the time dependence of the related perturbations always contains decreasing exponential parts, so that only stable modes are present, even if the purely oscillating perturbations are lost, as before. In particular we find either two exponentially decreasing modes or two damped oscillatory ones, according to the sign of the discriminant. We believe that it is just the case of a negative discriminant (corresponding to damped oscillating modes), namely

$$\begin{aligned} c_A^2 > \frac{1}{4} \eta ^2 k^2 \, , \end{aligned}$$

to be the most interesting from a physical point of view. Such a condition is generally fulfilled due to the strength of the magnetic field against the smallness of the resistivity.

In conclusion the effect of dissipative non-ideal processes is to damp the evolution of perturbations, preserving unchanged the key values of the (modified or not) Jeans wavenumbers. Meanwhile the same role is played by the presence of dissipative mechanisms due to the bulk viscosity, as in [4]; thus, besides the damping effects on the density perturbations, one could also study top-down schemes for (very) small and (very) strong resistive values from a numerical point of view.

3.3 The hyperbolic NIMHD setting versus the IMHD one, within the stationary MOG theory

We now limit our attention to the case \(\eta _r = 0\), \(\nu _D = 0\) and \(\tau _\phi \rightarrow 0\), but \(\eta _e = \eta > 0\) and \(\tau _j > 0\). The dispersion Eqs. (13) and (14) reduce to the forms

$$\begin{aligned}&\sigma ^3 + \frac{1}{\tau _j} \sigma ^2 + c_{ARE}^2 k^2 \sigma + \frac{c_A^2}{\tau _j} k^2 = 0 \, , \end{aligned}$$
(17)
$$\begin{aligned}&\sigma ^5 +\frac{1}{\tau _j} \sigma ^4 + \frac{k^2}{k^2 + {\varLambda }} (c_{ABE}^2 + c_s^2)(k^2 - k^2_{JE {\varLambda }}) \sigma ^3 \\&\quad + \frac{k^2}{k^2 + {\varLambda }} (c_{AB}^2 + c_s^2)(k^2 - k^2_{JAB{\varLambda }}) \frac{\sigma ^2}{ \tau _j} \\&\quad +\frac{k^4}{k^2 + {\varLambda }} c_s^2 c_{ARE}^2(k^2 - k^2_{J{\varLambda }}) \sigma \\&\quad + \frac{1}{\tau _j} \frac{k^4}{k^2 + {\varLambda }} c_s^2 c_A^2(k^2 - k^2_{J{\varLambda }}) = 0 \, , \end{aligned}$$
(18)

where \(c_{ABE}^2 = c_{ARE}^2 + c_B^2 \) and

$$\begin{aligned} k^2_{JE {\varLambda }} = k^2_{JE} - {\varLambda } = \displaystyle \frac{4 \pi G \rho _0}{c_s^2 + c_{ABE}^2} - {\varLambda }\,. \end{aligned}$$

Equation (17) describes three modified non-gravitational transverse Alfvén modes, which are always stable, even if their purely torsional behavior is quite lost, due to the damping effects related to the current relaxation time \(\tau _j\). Moreover, the larger Alfvén wavespeed suggests the presence of smaller time scales, useful to justify fast MR mechanisms. As expected, these modes are not affected by the self-gravitational potential and hence, focusing on collapse formation, we converge our attention to Eq. (18), which describes five gravito-magneto-elasto-sonic modes vs collapse formation.

For a rapid understanding, we first address the special cases of parallel and orthogonal propagation to the direction of \({{\mathbf {B}}_0}\). It is just an easy matter to show that the modified Jeans-Einstein criteria again remain unchanged either for \({{\mathbf {B}}}_0 \cdot {{\mathbf {n}}} \ne 0\) or \({{\mathbf {B}}}_0 \cdot {{\mathbf {n}}} = 0\). In fact, when \({{\mathbf {B}}}_0 \cdot {{\mathbf {n}}} = 0\), Eq. (18) splits in \(\sigma = 0\) (marginally stable mode) and a reduced fourth grade equation whose known term still include the factor \((k^2 - k^2_{J{\varLambda }})\), due to the elastic correction to the standard Alfvén speed.

The stability analysis may be analytically provided via the Routh-Hurwitz test, as done in the previous sections but, to best enhance the reliability of our hyperbolic NIMHD setting within the MOG theory, a numerical approach could deepen the discussion both mathematically and physically.

However to better compare the behavior of our NIMHD model with respect to the classical IMHD one (Lundquist model), we may carry out a further crucial simplification, working within weakly (\(\sigma \tau _j<< 1\)) and strongly (\(\sigma \tau _j>> 1\)) coupling limits [5]. Likewise to viscoelastic effects [7], in the strongly coupling limit the dispersion Eqs. (17) and (18) become

$$\begin{aligned}&\sigma ^2 + c^2_{AE} k^2 = 0 \\&\sigma ^4 + \frac{k^2}{k^2 + {\varLambda }} (c_{ABE}^2 + c_s^2)(k^2 - k^2_{JE {\varLambda }}) \sigma ^2 \\&\quad +\frac{k^4}{k^2 + {\varLambda }} c_s^2 c_{ARE}^2(k^2 - k^2_{J{\varLambda }}) = 0 \end{aligned}$$

whereas, in the weakly coupling limit, the equations remain the same with \(c_A^2\) in place of \(c^2_{AE}\). In this way we recover the torsional behavior of the (classical and shear-like) transverse Alfvén waves and read the gravitational response on the four coupled fast and slow (magneto-sonic and magneto-elasto-sonic) longitudinal modes.

Beyond the additional presence of the cosmological parameter \({\varLambda }\), again and differently from the classical conclusions as in [3], the stability analysis is not affected by the direction of \({\mathbf {B}}_0\). In particular, the asymptotic stability is recovered for all \(k^2\) such that \(k^2 \ge k^2_{J{\varLambda }}\), whereas gravitational collapse arises for \(k^2 < k^2_{JE {\varLambda }}\), thus leading to larger critical values \( \lambda _{JE{\varLambda }} = ({2 \pi })/k_{JE{\varLambda }}\) and \(M_{JE {\varLambda }}= ({4 \pi } \rho _0 \lambda ^3_{JE {\varLambda }})/{3}\) for the wavelength and mass, respectively.

4 Stabilizing interlacements among the most significant (dimensionless) numbers of the new theory

In astrophysical plasmas the interaction with the magnetic field may be discussed via the introduction of the so called plasma-beta (\(\beta \)), which measures the ratio of the thermal (plasma) pressure to the magnetic pressure. In the politropic context, we may define \(({\gamma } \beta )/{2}= {c_s^2}/{c_{AB}^2}\). We observe that for a perfectly conducting plasma, within the ideal Ohm’s law, \(\beta \) is the only parameter able to describe the strength of the magnetic field. Taking into account that astrophysical environments are often very different, \(\beta \) can vary in a broad range and hence we can meet both the physics of very high-\(\beta \) (>> 1) and very low-\(\beta \) (<< 1), including also the intermediate case \(\log \beta = 0\), corresponding to the interstellar medium. At a low-enough compression rate, the magnetic pressure becomes higher than the thermal pressure, namely the magnetic pressure is dominant and the β << 1 case is realized. On the contrary, for a high-enough compression rate, it is the gaseous pressure to become higher and the β >> 1 case occurs. The \(\beta \rightarrow \infty \) case corresponds to the simplified scheme discussed in Sect. 2.

The further presence of a finite (even small) plasma resistivity is read via the magnetic Reynolds number \(R_M = {(L c_s)}/{\eta }\), by taking the sound velocity as the characteristic one. The anomalous resistive effects, here incorporated through the \({\mathbf {R}}_e\) term, lead to the definition of the elasto-magnetic Reynolds number \(R_{ME} = R_M O_{RF}\), \(O_{RF}\) being the (new) dimensionless magnetic relaxation time, namely \(O_{RF} ={\tau _j}/{t_c}\), where \(\tau _J\) is the magnetic relaxation time and \(t_c ={L}/{c_s}\) is the characteristic time. As already mentioned, this correction allows for an enhancement of the role of \(\eta \), due to \({\tau _j}\), towards micro-instabilities, which are often excited only in reconnection regions, and this may be the key topic to justify fast MR processes.

Neglecting for simplicity the Darcy term, the previous results show that the largest effective wavespeed involved herein leads to the lowest critical wavenumber (15), which can be suitably rewritten as follows

$$\begin{aligned} k^2_{JRE{\varLambda }} = \displaystyle \frac{4 \pi G \rho _0}{c_s^2 + c_{ABRE}^2 } - {\varLambda } = \frac{k_j^2}{1 + d} \left[ 1- \frac{{\varLambda } c_s^2}{4 \pi G \rho _0} (1 + d) \right] \, , \end{aligned}$$

where

$$\begin{aligned} d = \displaystyle \left( \frac{2}{\gamma \beta } + \frac{1}{R_{ME}} \right) \,. \end{aligned}$$

This form suggests the possible experimental validity of the following interesting lower and upper bound estimates for the cosmological constant

$$\begin{aligned} \frac{4 \pi G \rho _0}{c_s^2(1 + d)}< {\varLambda } < \frac{4 \pi G \rho _0}{c_s^2} \, . \end{aligned}$$

Thanks to the above relation, we may conclude that these anomalous magnetic effects due to the Oldroyd-B type version of the resistive Ohm’s Law have a stabilizing role on the onset of gravitational instabilities but, even for transverse propagation, they do not alter the sufficient condition for a delayed collapse formation, previously recovered in Sect. 2 for hydrodynamic flows. Obviously this conclusion does not change in the presence of the dissipative Darcy term. In particular, within the ideal MHD setting (i.e. \(d = {2}/({\gamma \beta })\) ) for transverse propagation, we would recover stability always. Therefore the presence of ideal magnetic effects combined with the Einsteinian correction could completely stabilize the problem, so leading to novel Jeans-Einstein-Chandrasekhar criteria [3].

Also, by a quick inspection, when \({\mathbf {B}}_0 \cdot {\mathbf {n}} = 0\) and for \({\varLambda } \rightarrow 0^+\), it turns out that the lower Jeans-Chandrasekhar threshold is controlled by \(\beta \) as follows

$$\begin{aligned} k_{JC}^2 = \frac{\gamma \beta k_J^2}{\gamma \beta + 2}\, . \end{aligned}$$

For \(\log \beta = 0\) an upper bound for \( k_{JC}^2\) is \({k_{J}^2}/{2}\) within the adiabatic regime (\(\gamma = {5}/{3}\)), while \({k_{J}^2}/{3}\) is just its value in the isothermal case (\(\gamma = 1\)). Indeed, for very low-\(\beta \), \( k_{JC}^2\) is much lower than \( k_{J}^2\). On the other hand, if \(\beta \) is sufficiently high, for example in the limit \(\beta \rightarrow + \infty \), we find \(k_{JC}^2 \rightarrow k_{J}^2\), as expected.

5 The Kelvin–Voigt type Ohm’s Law within the MOG effects

It is noteworthy that, for constant parameters and within the linear theory, by employing the phase-lag concept together with “backwards” temporal shiftings [39], a Kelvin–Voigt type model is always derivable from a Jeffreys-type one (see e.g. [35]). In our scheme this procedure leads to a non stationary Ohm’s Law in the Kelvin–Voigt form

$$\begin{aligned} {\mathbf {E}}+ {\mathbf {v}}\times {\mathbf {B}}=\eta \nabla \times \left[ {\mathbf {B}}- \tau _j(1-r_j) {\mathbf {B}}_t\right] \end{aligned}$$
(19)

where, as before, \(\eta \mu _0\) is the inverse of the electrical conductivity and \(\tau _j(1-r_j)\) represents the difference delay time parameter, \(r_j\) being the so called (dimensionless) retardation time. As in the Kelvin–Voigt type thermal and viscoelastic contexts, the previous equation describes ”backwards responses” of \(\nabla \times {\mathbf {B}}\) on \({\mathbf {J}}\), so the well-posedness of this new MHD background, in the absence of the self-gravity potential, is always guaranteed by retarded responses, corresponding to the physically reasonable setting \(r_j > 1\).

Interestingly, in this case (19) may be also viewed as the linear version of the generalized Ohm’s Law investigated in [40] for MR solutions and hence it seems a significant proposal from an astrophysical point of view.

In the new frame the induction Eq. (9)\(_3\) reads

$$\begin{aligned} {{\mathbf {B}}}_t+\nabla \times ({\mathbf {B}} \times {\mathbf {v}}) =\eta \left( {\varDelta } {\mathbf {B}}- \tau _j(1-r_j) {\varDelta } {\mathbf {B}}_t \right) \end{aligned}$$

so that the spectral analysis arguments yield the vectorial equation

$$\begin{aligned} \left( \omega (1 - k^2 \eta \tau _j(1-r_j)) + \eta i k^2\right) {\mathbf {B}}_1 = k {\mathbf {n}} \times ({\mathbf {B}}_0 \times {\mathbf {v}}_1) \end{aligned}$$

accounting only for the transverse relations

$$\begin{aligned} \left( \omega (1 - k^2 \eta \tau _j(1-r_j)) + \eta i k^2 \right) B_{11}&= k v_{1n} B_{01} - k B_{0n} v_{11} \\ \left( \omega (1 - k^2 \eta \tau _j(1-r_j)) + \eta i k^2 \right) B_{12}&= k v_{1n} B_{02} - k B_{0n} v_{12} \end{aligned}$$

Henceforth, in order to converge to the point quickly, we work within the stationary MOG theory (\(\tau _\phi = 0\)), neglecting also the dissipative porous effects (\(\nu _D = 0\)). Starting from (11), it is a straightforward matter to recover again a Cramer system in \({\mathbb {R}}^3\) for the essential amplitudes \((v_{1n}, B_{11} , B_{12} )\), so that we read off the general dispersion equation, as follows

$$\begin{aligned} \alpha _{KV} \left[ \alpha _{KV}\left( \omega ^2 - \frac{c_s^2 k^2}{k^2 + {\varLambda }} (k^2 - k^2_{J {\varLambda }}) \right) - c_B^2 k^2 \omega ^2 \right] = 0 \end{aligned}$$
(20)

with

$$\begin{aligned} \alpha _{KV} = \omega ^2 \left[ (1 - k^2 \eta \tau _j(1-r_j)\right] + \eta i k^2 \omega - c_A^2 k^2 \end{aligned}$$

Equating to zero the first factor, we have

$$\begin{aligned} (1 - \omega ^2_{\eta j}) \sigma ^2 + \omega _\eta ^2 \sigma + \omega _A^2 = 0 \end{aligned}$$

where, for convenience, we have introduced the notations

$$\begin{aligned} \omega _{\eta j}^2 = \eta \tau _j (1-r_j)k^2, \quad \omega _\eta ^2 = \eta k^2, \quad \omega _A^2 = c_A^2 k^2 \end{aligned}$$

for the characteristic frequencies.

This equation describes the non-gravitational and non-compressional Kelvin–Voigt type Alfvén regimes and it is just the same equation as in [15] for the incompressible Alfvén modes. Meanwhile, in the simplified limit \(\omega _\eta ^2 \rightarrow 0\), we recover the standard purely oscillatory Aflvén modes of the IMHD theory, whereas, for \(\omega _{\eta j}^2 \rightarrow 0\) and \(\omega _{\eta }^2 > 0\), we obtain again the two damped Alfvén modes, driven by the resistive NIMHD frame (see e.g. Sect. 3.2).

More generally, the equation admits the solutions

$$\begin{aligned} \sigma _{\pm } = \frac{1}{2(1 - \omega ^2_{\eta j})} \left[ - \omega _\eta ^2 \pm \sqrt{{\varDelta }} \right] \end{aligned}$$

with

$$\begin{aligned} {\varDelta } = \omega _\eta ^4 - 4(1 - \omega ^2_{\eta j})\omega _A^2 \end{aligned}$$

As we can see, the stability discussion is now affected by the sign of \((1-\omega ^2_{\eta j})\) and so also unstable modes arise. By inspection, when \((1 -\omega ^2_{\eta j}) < 0\) we always have \({\varDelta } >0\), namely two purely real values for \(\sigma \) of different signs, corresponding to exponentially growing and decreasing regimes (i.e. instability). On the other hand when \((1 -\omega ^2_{\eta j} )> 0\) since \(r_j > 1\), the sign of \({\varDelta }\) is strictly related to the limit cases of very ”strong” magnetic field and ”small” anomalous resistivity, due to the presence of \(\omega ^2 _{\eta j}\). Anyway, independently of \(k^2\), we register non-compressible (asymptotically) stable modes. More precisely, the new constitutive theory is responsible of two damped oscillatory modes for \({\varDelta } > 0\), whereas \(for {\varDelta } \le 0\) two exponentially damped modes are allowable. Again, the two purely oscillating modes are lost, due to the presence of the Kelvin–Voigt type resistivity.

In particular, for retardation effects, it is likely that \({\varDelta } < 0\) be the most physically relevant case, even thinking of fast MR mechanisms [40]. Meanwhile the presence of these damped Alfvén oscillatory modes is now dictated by the condition

$$\begin{aligned} \omega _A^2 > \frac{\omega _\eta ^4}{4(1- \omega ^2_{\eta j})} \end{aligned}$$

which weakens the analogous requirement, exhibited in the purely resistive context.

Indeed the same stability results follow even for \(r_j < 1\), when the positiveness of \((1 -\omega ^2_{\eta j}) \) is assured by the condition \(k^2 <{1}/[\eta \tau _j (1-r_j)]\), so that the previous unstable modes are recovered for \(k^2 > {1}/[{\eta \tau _j (1-r_j)}]\).

In analogy with the previous setting, the gravitational response of the Kelvin–Voigt type Ohm’s law rests on the second factor in (20), which may be rewritten as

$$\begin{aligned} \begin{aligned}&(1 - \omega _{\eta j}^2)\sigma ^4 + \omega _\eta ^2 \sigma ^3 + \frac{k^2(1 - \omega _{\eta j}^2) \omega ^2_{JKV {\varLambda }}}{k^2 + {\varLambda }} \sigma ^2 \\&\quad + \frac{k^2\omega ^2_\eta \omega ^2_{J {\varLambda }} }{k^2 + {\varLambda }} \sigma + \frac{k^2 \omega ^2_A \omega ^2_{J {\varLambda }}}{k^2 + {\varLambda }} = 0 \end{aligned} \end{aligned}$$
(21)

Here, focusing on the case \((1 - \omega _{\eta j}^2) > 0\), we have reasonably defined the critical frequencies

$$\begin{aligned} \begin{aligned}&\omega ^2_{J {\varLambda }} = \omega ^2_s + {\varLambda } c_s^2 - 4 \pi G \rho _0 \,, \\&\omega ^2_{J KV {\varLambda }} = \omega ^2_{KV} + \omega ^2_s + {\varLambda }(c^2_{KV} + c_s^2) - 4 \pi G \rho _0 \end{aligned} \end{aligned}$$

with

$$\begin{aligned} c^2_{KV} = \displaystyle \frac{c^2_{AB} }{1 - \omega ^2_{\eta j}} \end{aligned}$$

being the modified (squared) magnetic wavespeed, due to the new constitutive theory.

As expected, the interplay between \({\varLambda }\) and \(r_j\) turns out to be strategic in the discussion. In particular, in the limit case \(\omega ^2_{J {\varLambda }} = 0\) (i.e. \( k^2 = k^2_{J {\varLambda }}\)), besides the double root \(\sigma = 0\) (corresponding to marginally stable modes), Eq. (21) reduces to

$$\begin{aligned} (1 - \omega _{\eta j}^2)\sigma ^2 + \omega _\eta ^2 \sigma + \frac{k^2}{k^2 + {\varLambda }} (1 - \omega _{\eta j}^2) ( \omega ^2_{KV} + {\varLambda } c^2_{KV} )= 0 \end{aligned}$$

so that two (exponentially stable) Kelvin–Voigt type purely magnetic regimes, unaffected by the gravitational effects, always follow.

More generally, Kelvin–Voigt type gravito-magneto-sonic waves are allowable whenever \(\omega ^2_{J {\varLambda }} > 0\), the finite- and zero- resistive frames being recovered for \(\omega ^2_{\eta j} = 0\) but \(\omega ^2_\eta > 0\) and for \(\omega ^2_\eta = 0\), respectively. To sum up, we find four unstable coupled compressible modes for \(\omega ^2_{JKV{\varLambda }} < 0\), whereas four (asymptotically) stable coupled regimes are always possible when \(\omega ^2_{J{\varLambda }} \ge 0\). On the other hand, whenever \((1 - \omega ^2_{\eta j}) < 0\) and in compliance with other single-phase lag theories working with a difference delay parameter in thermo-viscoelasticity, the backwards Ohm’s law yields four unstable coupled modes, independently of (modified or not) gravitational effects.

It is worth to note that \(\omega ^2_\eta \) and \(\omega ^2_{\eta j}\) are the inverse time scales of the novel Ohm’s theory: the electrical conductivity is likely to be very slow in protostellar disks, due to a low temperature and hence, since the magnetic field lines will not be frozen-in during collapse, particle collisions within the plasma are accounted for, as it happens experimentally [41].

In our opinion, this special mathematical setting might be suitable to model microphysical processes in most astrophysical plasmas [42].

6 Conclusions

  • In the present paper we have discussed the effects of a non-relativistic gravity with the cosmological constant \({\varLambda }\) firstly on the sound modes of interstellar gas clouds and then in the presence of special non-ideal magnetic properties, allowing for lagging behaviors of the magnetic field.

  • The presence of \({\varLambda }\) is responsible of a new critical wavenumber ascribable to DM; as a consequence, we may define suitable new critical thresholds for wavenumber, wavelength and mass, herein referred to as the Jeans–Einstein thresholds, which stabilize the medium against collapse formation.

  • By the way, the interaction with the magnetic field may be discussed via the introduction of the so called plasma-beta (\(\beta \)). More specifically, a key role of the new theory is played by the Ohm number \(O_{RF}\), representing the dimensionless magnetic relaxation time, which modifies the Reynolds number \(R_M\) in an elastic Reynolds number \(R_{ME}\), in order to control the anomalous resistive features towards the micro-physics of cosmic plasmas. Also, the larger Alfvén effective wavespeed due to our theory leads to shorter Alfvén transitional/crossing times to be compared with standard diffusion times via a modified Lundquist number S, which seems more favorable to capture the micro-to macro- scales properties of fast MR processes (see [37, 39]).

  • We highlight an interesting relationship among all the significant (dimensionless) numbers of the theory, whose validity, within the ideal limit and for transverse propagation, allows for the full stabilization of the medium and hence leads to a revision of the Jeans–Chandrasekhar instability criteria. We have thoroughly analyzed the stabilizing effects of all the physical parameters herein involved against the onset of gravitational instability, even under both strongly and weakly coupling limits, as in viscoelastic settings. Interestingly, differently from other (ideal or resistive) MHD frameworks, our modified criteria are unaffected by the direction of the equilibrium magnetic field \({\mathbf {B}}_0\), namely they are not influenced by either longitudinal and transverse modes.

  • As a result of a delayed collapse formation, the gas cloud would have fragmentations, thus stimulating star formation, and the presence of further small (gravitational and magnetic) relaxation times could be responsible for the intricate network of filamentary structures of the ISM, where the large amount of DM may be hidden, as recently observed in [43]. Given the importance of magnetic fields within the formation and the evolution of filaments, this cosmic web may be ascribed to the enhanced twisting properties displayed by the non gravitational elastic-type torsional Alfvén waves, here exhibited. In the light of their recent observational evidence, these modified, even damped Alfvén waves propagating along thin current sheets, besides their undisputed involvement in solar corona heating problems, could be an important solar interplanetary driver of the global thermospheric perturbations and may receive much attention also in the exciting field of solar magnetoseismology.

  • A special attention has been devoted to stable longitudinal (fast and slow) elasto-magneto-sonic waves which, differently from Alfvén waves, are affected by (modified or not) gravity.

  • Another dimensionless number, the so called retardation time, arises for the Kelvin–Voigt type Ohm’s Law. Hence the ensuing revisited time scalings approach results more reliable to justify also the rapid magnetic energy release on those non-ideal processes likely to occur in regions of fast MR (such as solar flares and spots, particle acceleration and coronal matter ejection) [40].

  • This rheological scheme for self-gravitating plasmas furnishes not only a more realistic description of gravitational collapse formation, but anomalous resistivity (AR), through more favorable timescales, may shed new light on the theoretical study of MR by bridging the gap between micro- to macro-scale dynamics in narrow diffusion regions around X-type magnetic nulls.