© ESO 2018

1. Introduction

One of the most intriguing features of magneto-hydrodynamics (MHD) consists of the generation of unstable behaviors from the coupling between wave propagation and plasma inhomogeneities (Biskamp 1993). A relevant example of such instability mechanisms corresponds to the so-called magneto-rotational instability (MRI), firstly discovered by Velikhov (1959) and applied to the astrophysical context by Chandrasekhar (1960); it results from the coupling of Alfvén waves with the differential rotation of the plasma.

It is easy to understand the interest that MRI has raised over the last five decades toward its astrophysical applications, especially since it was clarified (Balbus & Hawley 1991; see also Balbus & Hawley 1998 and Balbus 2003) that it constitutes the basic instability able to trigger the turbulence in stellar accretion disks. Indeed, the existence of MRI in accreting structures is crucial to justify the assumption of an effective visco-resistive plasma, which forms part of the basis of the original Shakura idea of accretion onto compact objects (Shakura 1973; Shakura & Sunyaev 1973). Most of the studies on MRI have been locally developed, that is, they consider small wave lengths of the perturbations with respect to the typical variation scale of the background configuration. For a satisfactory discussion of MRI in the global approach, see Papaloizou & Szuszkiewicz (1992). The most important results concerning MRI come from the analysis of the background plasma profiles, according to the idea that vertical shear weakly couples to the perturbation dynamics. Nonetheless, for recent developments concerning MRI in stratified configurations, revising the original approach of Balbus (1995), see Montani et al. (2016) and Cianfrani & Montani (2017). Finally, for studies of the morphology of MRI in the presence of dissipative effects, like viscosity and resistivity, we refer to the reviews by Balbus & Hawley (1998), and Shakura & Postnov (2015) and to Carlevaro et al. (2017) for recent investigations.

In this paper, we apply MRI to a different context with respect to stellar accretion disks, to examine whether or not it could play a role in the structure formation across the Universe. Indeed, the plasma nature of the Universe before recombination and the existence of primordial magnetic fields (Barrow et al. 2007) suggest that the cosmological perturbation dynamics, and therefore the structure formation, can be influenced by MHD effects. In particular, in the recombined-matter-dominated Universe, we can expect that differential rotation significantly enters the nonlinear gravitational collapse of overdensities across the Universe. In this respect, it is important to stress that the Universe remains ionized even after recombination for one part in one hundred thousand. Therefore, the emergence of MRI in the linear and nonlinear phases of the overdensity collapse is a reliable phenomenon for the early Universe, especially due to the tight coupling between ions and neutral hydrogen on many relevant cosmological scales in the era between the recombination and structure formation (Lattanzi et al. 2012). However, the cosmological implementation (provided in Sect. 6) of the present results, obtained for a steady linearized MHD system, has to be regarded as a qualitative hint for further investigations which include the detailed role of the Universe expansion and the relativistic perturbation dynamics.

In Sadhukhan et al. (2016), the role played by MRI in the linear stability of a plasma infinite filament is considered in order to investigate possible implications on the morphology of the filament in different regimes. For further studies on the interaction between self-gravitating and magnetic instabilities, see for instance Fromang et al. (2004a,b,c). Here, we consider the same type of problem, but we refer to a nebula of plasma treated as a differentially rotating nonstratified configuration, in order to evaluate how the corresponding Jeans scale (Jeans 1902) is affected by MRI. The main issue of such an analysis is to demonstrate how, starting with an axisymmetric profile (i.e., a configuration of plasma in differential rotation within a weak magnetic field), the linear gravitational collapse is able to induce anisotropic features of the overdensity, up to some percent of the background density and of order unity for the perturbation itself. This linear and local investigation suggests that the extreme nonlinear regime of the collapse could be significantly affected by such initial anisotropies and the emergence of filament-like structures could be justified by the nonlinear evolution of such a linear initial condition. The coupling of the MRI with the Jean mechanism of self-gravitational collapse (i.e., the instability due to the presence of differential rotation in a weakly magnetized plasma) could, in principle, influence the process of fragmentation of the nebula since we show that the critical scale of instability is fixed both by gravitational and magneto-rotational effects. However, we clarify how the anisotropy feature of the observed linear growth rates is essentially due to an intrinsic property of the pure Jeans mechanism. Indeed, the largest growth-rate values correspond to small perturbation wave numbers, where MRI is almost suppressed. On the other hand, those wave numbers corresponding to a leading character of MRI exist is the region where the Jeans growth rates strongly decrease: such two instabilities turn out to coexist within the nebula instability profile. For a discussion of similar questions in the general relativistic sector, see Germani & Tsagas (2006), where the evolution of small linear anisotropies is investigated in a specific case, demonstrating their enhancement by the gravitational collapse.

2. Basic equations

Let us consider a self gravitating plasma axial structure (nebula) described in standard cylindrical coordinates (r, φ, z), whose self-gravity potential is denoted by ϕ. The nebula is taken in differential rotation with angular velocity Ω = Ω(r), and it is embedded into a purely vertical magnetic field ${\mathit{B}}_0=B_0{\widehat{\mathit{e}}}_{\mathrm{z}}$. Since we are considering a weak dependence on the vertical direction of all the background quantities, they are taken as a function of r only and the nebula stratification is, on this level, neglected.

The theory properly describing the physics governing such a system is therefore the ideal MHD, whose fundamental equations are the mass, momentum, and magnetic flux conservation laws$$\begin{array}{cc}& 0={\partial }_t\rho +\mathrm{\nabla }·(\rho \mathit{\nu }),\hfill \end{array}$$(1) $$\begin{array}{cc}& 0=\rho {\partial }_t\mathit{\nu }+\rho (\mathit{\nu }·\mathrm{\nabla })\mathit{\nu }+\mathrm{\nabla }(p+B^2/8\pi )-\rho \mathrm{\nabla }\varphi -(\mathit{B}·\mathrm{\nabla })\mathit{B}/4\pi ,\hfill \end{array}$$(2) $$\begin{array}{cc}& 0={\partial }_t\mathit{B}-\mathrm{\nabla }\wedge (\mathit{\nu }\wedge \mathit{B}),\hfill \end{array}$$(3)

respectively. Here, ρ denotes the mass density of the fluid, p the hydrostatic pressure, while ν and B (B = | B|) are the velocity and magnetic field, respectively. In addition, we consider the Poisson equation describing the self-gravity of the nebula, that is,$$\begin{array}{c}\hfill 0={\mathrm{\nabla }}^2\varphi -4\pi G\rho ,\end{array}$$(4)

with G being the Newton constant and, finally, we specify the equation of state (EoS) as$$\begin{array}{c}\hfill 0=p-{\nu }_{\mathrm{s}}^2\rho ,\end{array}$$(5)

where ν_s is the sound speed and we consider an isothermal relation between pressure and mass density (ν_s ≃ const.).

2.1. Background equilibrium

Let us now describe the fundamental equations governing the time-independent background equilibrium, characterized by a set of variables indicated via the subscript 0, in order to distinguish from the dynamical perturbed quantities of the following section denoted by a subscript 1.

Assuming that the background magnetic field B₀ is sufficiently small to have no effect on the steady equilibrium of the plasma, we can easily fix the equation governing its gravostatic configuration as$$\begin{array}{cc}& 0=r{\mathrm{\Omega }}^2(r)-({\partial }_r{p}_0)/{\rho }_0-{\partial }_r{\varphi }_0,\hfill \end{array}$$(6) $$\begin{array}{cc}& 0=({\partial }_r(r{\partial }_r{\varphi }_0))/r-4\pi G{\rho }_0(r),\hfill \end{array}$$(7) $$\begin{array}{cc}& 0={\partial }_r{p}_0-{\nu }_{\mathrm{s}}^2{\partial }_r{\rho }_0.\hfill \end{array}$$(8)

Here, we neglected the z dependence of the problem, while all the φ derivatives vanish because of axial symmetry. Actually, the z dependence is assumed negligible here, bearing in mind that the perturbations will have a small wavelength, making it difficult to explore the vertical shear of the configuration. Nonetheless, the role of the vertical gradients can be, in principle, important in order to fix the background profile. In this case, the system above must include vertical gradients as well as the vertical force balance.

3. Linear perturbation theory

Let us now face the problem of characterizing the linear stability of the considered plasma configuration. In what follows, we consider a local approximation, that is, the typical perturbation length is considered much smaller than the length scale of the background variation. The background quantities and their radial derivatives are taken at a fiducial radius $r=\overline{r}$, and they are treated as constants in the considered perturbative problem. According to this approximation, the perturbation equations, associated to the system of Eqs. (1)–(5), are written as$$\begin{array}{cc}& 0={\partial }_t{\rho }_1+{\rho }_0(\mathrm{\nabla }·{\mathit{\nu }}_1)+{\mathit{\nu }}_0·\mathrm{\nabla }{\rho }_1,\hfill \end{array}$$(9) $$\begin{array}{cc}& 0={\rho }_0({\partial }_t{\mathit{\nu }}_1+({\mathit{\nu }}_0·\mathrm{\nabla }){\mathit{\nu }}_1)+\mathrm{\nabla }(p_1+2{\mathit{B}}_0·{\mathit{B}}_1/8\pi )\hfill \\ \hfill & -{\rho }_0\mathrm{\nabla }{\varphi }_1-({\mathit{B}}_0·\mathrm{\nabla }){\mathit{B}}_1/4\pi ,\hfill \end{array}$$(10) $$\begin{array}{cc}& 0={\partial }_t{\mathit{B}}_1+{\mathit{B}}_0(\mathrm{\nabla }·{\mathit{\nu }}_1)-({\mathit{B}}_0·\mathrm{\nabla }){\mathit{\nu }}_1+({\mathit{\nu }}_0·\mathrm{\nabla }){\mathit{B}}_1,\hfill \end{array}$$(11) $$\begin{array}{cc}& 0={\mathrm{\nabla }}^2{\varphi }_1-4\pi G{\rho }_1,\hfill \end{array}$$(12) $$\begin{array}{cc}& 0=p_1-{\nu }_{\mathrm{s}}^2{\rho }_1,\hfill \end{array}$$(13)

respectively.

Since the background does not depend on time, we search WKB solutions writing each axisymmetric perturbation quantity A₁ as $A_1(t,r,z)=\overline{a}{\mathrm{e}}^{i(k_{r}r+k_{z}z-\omega t)}$ with $\overline{a}=\mathrm{const}.$ and the wave vector as k = (k_r, 0, k_z). In this way, the following natural replacements take place: ∇A₁ = ikA₁ and ∂_tA₁ = − iωA₁. The system of Eqs. (9)–(13) can be rewritten as$$\begin{array}{cc}& 0=-i\omega {\rho }_1/{\rho }_0+i{k}_r{\nu }_{1r}+i{k}_z{\nu }_{1z},\hfill \end{array}$$(14) $$\begin{array}{cc}& 0=-i\omega {\nu }_{1r}-2\mathrm{\Omega }{\nu }_{1\varphi }+i{k}_r{p}_1/{\rho }_0+i{k}_r{\varphi }_1\hfill \\ \hfill & +i{k}_r{\nu }_A^2{B}_{1z}/B_0-i{k}_z{\nu }_A^2{B}_{1r}/B_0,\hfill \end{array}$$(15) $$\begin{array}{cc}& 0={\kappa }^2{\nu }_{1r}/2\mathrm{\Omega }-i\omega {\nu }_{1\varphi }-i{k}_z{\nu }_A^2{B}_{1\varphi }/B_0,\hfill \end{array}$$(16) $$\begin{array}{cc}& 0=-i\omega {\nu }_{1z}+i{k}_z{p}_1/{\rho }_0+i{k}_z{\varphi }_1,\hfill \end{array}$$(17) $$\begin{array}{cc}& 0=-i\omega B_{1r}-i{k}_z{B}_0{\nu }_{1r},\hfill \end{array}$$(18) $$\begin{array}{cc}& 0=-i\omega B_{1\varphi }-B_{1r}{\partial }_{lnr}\mathrm{\Omega }-i{k}_z{B}_0{\nu }_{1\varphi },\hfill \end{array}$$(19) $$\begin{array}{cc}& 0=-i\omega B_{1z}+i{k}_r{B}_0{\nu }_{1r},\hfill \end{array}$$(20) $$\begin{array}{cc}& 0=-k^2{\varphi }_1-4\pi G{\rho }_1,\hfill \end{array}$$(21) $$\begin{array}{cc}& 0=p_1-{\nu }_{\mathrm{s}}^2{\rho }_1,\hfill \end{array}$$(22)

where we introduced the standard definitions of the epicyclic frequency and Alfén velocity$$\begin{array}{c}\hfill {\kappa }^2=4{\mathrm{\Omega }}^2+{\partial }_{lnr}{\mathrm{\Omega }}^2,{\nu }_{\mathrm{A}}^2=B_0^2/4\pi {\rho }_0,\end{array}$$(23)

respectively.

It is worth reiterating that the Rayleigh criterion (Balbus & Hawley 1998), in the absence of differential rotation, states that the nebula is unstable if κ² < 0. Therefore, in order to select the role played by MRI in the stability of the system, in what follows we take positive epicyclic frequencies κ² > 0. The system above is linear algebraic and homogeneous in the perturbed quantities, and therefore it admits a nontrivial solution only if the corresponding determinant vanishes, therefore leading to the dispersion relation$$\begin{array}{cc}\hfill 0& ={\omega }^6-{\omega }^4({\kappa }^2+{\omega }_0^2+{\omega }_A^2+{\omega }_{\mathit{Az}}^2)\hfill \\ \hfill & +{\omega }^2({\omega }_{\mathit{Az}}^2(2{\omega }_0^2+{\omega }_A^2+{\partial }_{lnr}{\mathrm{\Omega }}^2)+{\kappa }^2{\omega }_0^2{\omega }_{\mathit{Az}}^2/{\omega }_A^2)\hfill \\ \hfill & -{\omega }_0^2{\omega }_{\mathit{Az}}^2({\omega }_{\mathit{Az}}^2+({\partial }_{lnr}{\mathrm{\Omega }}^2){\omega }_{\mathit{Az}}^2/{\omega }_A^2).\hfill \end{array}$$(24)

Here, to simplify the notation, we introduce the following frequencies: ${\omega }_0^2={\nu }_{\mathrm{s}}^2{k}^2-4\pi G{\rho }_0$ is the typical frequency appearing in the Jeans (self-gravitation) instability, ${\omega }_A^2={\nu }_{\mathrm{A}}^2{k}^2$ is the Alfvén frequency and ${\omega }_{\mathit{Az}}^2={\nu }_{\mathrm{A}}^2{k}_z^2$ represents the Alfvén parameter ${(\mathit{k}·({\nu }_{\mathrm{A}}{\widehat{\mathit{e}}}_z))}^2$, which modulates the magnetic tension.

We can now rewrite Eq. (24) as a function of the angle θ between k and B₀. Therefore, we define $\chi ={\omega }_{\mathit{Az}}^2/{\omega }_A^2=k_z^2/k^2={cos}^2\theta$ and, substituting this expression in Eq. (24), we finally get$$\begin{array}{cc}\hfill 0& ={\omega }^6-{\omega }^4({\kappa }^2+{\omega }_0^2+{\omega }_A^2(1+\chi ))\hfill \\ \hfill & +{\omega }^2({\omega }_A^2(2{\omega }_0^2+{\omega }_A^2+{\partial }_{lnr}{\mathrm{\Omega }}^2)+{\kappa }^2{\omega }_0^2)\chi \hfill \\ \hfill & -{\omega }_0^2{\omega }_A^2({\omega }_A^2+{\partial }_{lnr}{\mathrm{\Omega }}^2){\chi }^2,\hfill \end{array}$$(25)

where, clearly, 0 ≤ χ ≤ 1. We stress that, turning off gravity, i.e., G → 0, taking ${\nu }_{\mathrm{s}}^2\to \infty$ (which is the limit toward the Boussinesq approximation Balbus & Hawley 1998) and fixing $\mathit{k}=k_z{\widehat{\mathit{e}}}_z$, the expression above reduces to the dispersion formula for the MRI (Balbus & Hawley 1991).

4. Discussion of the dispersion relation

Before analyzing the solutions of the obtained dispersion relation in detail, we investigate some relevant simplified cases.

4.1. The limit χ = 0

We start by considering a perturbation $\mathit{k}=k_r{\widehat{\mathit{e}}}_r$, implying χ = 0. The dispersion relation takes the following simplified form.$$\begin{array}{c}\hfill {\omega }^2-({\kappa }^2+{\omega }_0^2+{\omega }_A^2)=0.\end{array}$$(26)

The instability condition ω² < 0 gives$$\begin{array}{c}\hfill k_r^2<{\overline{k}}_J^2\equiv \frac{4\pi G{\rho }_0-{\kappa }^2}{{\nu }_{\mathrm{s}}^2+{\nu }_{\mathrm{A}}^2},\end{array}$$(27)

for which the associated Jeans length is$$\begin{array}{c}\hfill {\overline{\lambda }}_J\equiv \frac{2\pi }{{\overline{k}}_J}=2\pi \sqrt{\frac{{\nu }_{\mathrm{s}}^2+{\nu }_{\mathrm{A}}^2}{4\pi G{\rho }_0-{\kappa }^2}}·\end{array}$$(28)

This solution is physically meaningful only if κ² < 4πGρ₀: when this relation is satisfied, Eq. (28) provides the Jeans scale for a rotating, self-gravitating, and magnetized plasma in correspondence to a radial propagation of the perturbations (absence of magnetic tension). Furthermore, the limit κ² → 0 and ${\nu }_{\mathrm{A}}^2\to 0$ reproduces the standard Jeans criterium.

From Eq. (28), we see that the magnetic field contrasts the collapse, since the corresponding Jeans length is larger than the nonmagnetized one. The epicyclic frequency of the fluid contrasts the collapse too and can make the system stable for any perturbation if κ² ≥ 4πGρ₀. This squared frequency is always intended to be positive, in order to satisfy the Rayleigh criterion. However, negative values for κ² can increase the instability of the system, inducing a collapse at smaller wavelengths.

4.2. Perturbative solution of the dispersion relation

Since the aim of our analysis is to verify the role played by MRI in the stability of a self-gravitating medium, we consider Eq. (25) nearby a static gas cloud by turning on a small magnetic field and epicyclic frequency. More precisely, this means ${\omega }_0^2\gg {\omega }_A^2$ and ${\omega }_0^2\gg {\kappa }^2$. Moreover, we take χ ≪ 1 in order to have small magnetic tension. The differential rotation parameter is taken at half magnitude between ${\omega }_0^2$ and the other frequencies to preserve its relevant role in the MRI instability, i.e., ${\omega }_0^2>{\partial }_{lnr}{\mathrm{\Omega }}^2\gg {\omega }_A^2,{\kappa }^2$.

Under these assumptions, keeping only the first-order perturbations, the dispersion formula can be rewritten as$$\begin{array}{c}\hfill {\omega }^4-{\omega }^2[{\kappa }^2+{\omega }_0^2+{\omega }_A^2(1+\chi )]+\frac{{\omega }_A^2}{{\omega }_0^2}{\partial }_{lnr}{\mathrm{\Omega }}^2=0.\end{array}$$(29)

Considering the solution with the + sign, because it represents a modified Jeans frequency, we get$$\begin{array}{c}\hfill {\omega }^2={\omega }_0^2+{\kappa }^2+{\omega }_A^2(1+\chi )+\frac{{\omega }_A^2}{{\omega }_0^2}{\partial }_{lnr}{\mathrm{\Omega }}^2.\end{array}$$(30)

Given ω₀² < 0, the role of MRI naturally emerges. As magnetic tension and differential rotation are turned on, the root tends to be less negative if ∂_{ln r}Ω² < 0, which is the condition for MRI. Therefore, the obtained instability mitigates the pure Jeans one. Instead, a positive differential rotation parameter makes the system even more unstable.

4.3. General case

Let us now analyze the dispersion relation by plotting its numerical solutions in terms of the model parameters. In what follows, we expect to be able to recognize, in the profile of the growth rate γ (the positive imaginary part of ω) as a function of k, both the MRI and the Jeans behavior. Naming the unstable solutions ${\omega }_{\mathrm{J}}^2$ and ${\omega }_{\mathrm{M}}^2$ (since their behaviors match Jeans and MRI modes, respectively), we consider the hybrid root defined as ${\omega }_{\text{MJ}}^2(k)=\min [{\omega }_{\mathrm{J}}^2,{\omega }_{\mathrm{M}}^2]$, where MJ stands for magneto-Jeans. The critical wavelength of the system is therefore set as λ_M​J = min[λ_J, λ_M].

As can be inferred by the upper panel of Fig. 1, the collapse (the growth of the mass density) is significantly anisotropic since it depends on the ratio $\chi =k_z^2/k^2$. Anisotropy is also increased by differential rotation, as can be argued from the lower panel of the figure, which, moreover, pushes the critical wavelength to smaller scales (however, we note that ∂_{ln r} Ω ² must always be negative in order to have an MRI contribution Balbus & Hawley 1991). Moreover, as can be seen from the upper panel of Fig. 1, no dependence of instability critical scale on the value of χ emerges, except the discontinuity between χ = 0 and χ ≠ 0.

In Fig. 2, the critical length λ_M​J is plotted as a function of the ratio (ν_A/ν_s)², associated to the change of the magnetic field strength amplitude.

We note that MRI is valid for weak magnetic field only, but it implies no matter transport, and, if coupled with self-gravity, it extends the gravitational collapse to smaller scales, as clarified in the figure.

It is worth noting that the upper panel of Fig. 1 shows that, for a given k, the growth rate is larger for modes with χ = 1 (dotted line) than for modes with χ = 0 (solid line). Moreover, as clearly deduced from Fig. 2, at small values of k (large scales), the dominant instability is the Jeans one, whereas for large k (small scales), the dominant instability is associated to the MRI contribution. Below, we see how these considerations concern the analysis of the density-contrast evolution.

5. Time evolution of the density contrast

Let us now estimate the amount of the collapse anisotropy by studying the behavior of the density contrast defined as δ(r, t) = ρ₁/ρ₀ (where $\mathit{r}=(r{\widehat{\mathit{e}}}_r+z{\widehat{\mathit{e}}}_z)$ is the poloidal vector radius). We consider an initial overdense region of Gaussian form, i.e.,$$\begin{array}{c}\hfill \delta (\mathit{r},0)={\delta }_0{\mathrm{e}}^{-\frac{r^2+z^2}{2\sigma }},\end{array}$$(31)

whose Fourier transform is still Gaussian,$$\begin{array}{c}\hfill \tilde{\delta }(\mathit{k},0)={\tilde{\delta }}_0{\mathrm{e}}^{-\frac{k_r^2+k_z^2}{2\tilde{\sigma }}};\end{array}$$(32)

above, σ and $\tilde{\sigma }$ are the direct space and Fourier space variance, being reciprocal numbers (δ₀ and ${\tilde{\delta }}_0$ are assigned constants).

In the k-space, the single mode evolution is described by$$\begin{array}{c}\hfill \tilde{\delta }(\mathit{k},t)={\tilde{\delta }}_0{\mathrm{e}}^{-\frac{k^2}{2\tilde{\sigma }}+i\mathit{k}·\mathit{r}-i{\omega }_{MJ}t},\end{array}$$(33)

where ω_M​J = ω_M​J(k) is the MJ frequency derived in the previous section as a solution of the dispersion relation. Taking λ > λ_M​J, ω_M​J = iγ_M​J and perturbations grow exponentially. Transforming back Eq. (33) into the r-space, we can obtain the density contrast as a function of position until δ ≪ 1. In fact, when ρ₁ ≃ ρ₀, the linear analysis is clearly no longer valid and simulations of the nonlinear problem are required.

5.1. Evolution of the overdensity

The wave packet describing the evolution of the overdense region is now evaluated for different values of r. In this way, it is possible to investigate the dynamics of the collapse geometry by plotting the isolines of δ(r, t) on the r − z plane at different times (the system is invariant under translation along φ). The anisotropic behavior of the overdensity growth rate is evident from its dependence on the angle χ, and therefore we see the distortion of the isolines as time goes by, according to Fig. 3. Here, the numerical integration of the packet is performed taking ${\tilde{\delta }}_0=0.1\times e^{-k^2/2}$.

When t = 0, the isoline corresponds to a circle, but the shape is progressively altered as time increases, because the growth rate is larger along the direction of the background magnetic field. In fact, in Fig. 3, the isoline δ = 0.1 stretches along the z axis. Given a distance d, δ(r = d, z = 0, t)< δ(r = 0, z = d, t), i.e., the blob is more dense on the z axis and therefore it prevalently compresses in that direction. Furthermore, we note that the discontinuity in the growth rate derivative (see Fig. 1) has no visible effects on the density contrast, although this bend corresponds to a jump in the group velocity ∂_kω_M​J.

The shape of the isolines does not change if we choose another value for ${\tilde{\delta }}_0$, since it simply changes a constant that multiplies the Fourier transformation. We observe that we cannot overcome ${\tilde{\delta }}_0\sim 0.1$, otherwise we violate ρ₁ ≪ ρ₀ and the linear approximation is no longer valid.

5.2. Evaluation of the anisotropic growth

Let us now construct a coefficient in order to evaluate the anisotropy of the growing instability:$$\begin{array}{c}\hfill \alpha (t)=\frac{\delta (r=0,z=d,t)}{\delta (d,0,t)},\end{array}$$(34)

where d is a value of the distance from the center fixed with respect to the isoline δ = 0.1 on the z axis and at a given time, namely δ(0, d, 1.25) = 0.1. We can see that the quantity α can reach large values; for instance, at the chosen instant t = 1.25, we get α(1.25) ≃ 0.5. It must be stressed that α measures the anisotropy of the density contrast only, but it allows us to evaluate the whole density anisotropy. In fact, we can define an anisotropy coefficient for the whole density as$$\begin{array}{c}\hfill A(t)=1-\frac{\rho (d,0,t)}{\rho (0,d,t)},\end{array}$$(35)

and therefore this coefficient reaches 1 (and the anisotropy is of the order 100%) for ρ(d, 0, t)/ρ(0, d, t) → 0; it can be rewritten as function of α as$$\begin{array}{c}\hfill A(t)=1-\frac{1+\alpha \delta (0,d,t)}{1+\delta (0,d,t)}·\end{array}$$(36)

It is easy to evaluate that A(1.25) ≃ 5%, meaning that the onset collapse is just slightly anisotropic according to the considered linear regime (ρ₁ ≪ ρ₀).

It important to stress here that if we consider an initial perturbation where all modes are excited, the modes that grow faster are those with small k values, since the corresponding growth rate is larger (see Fig. 1). Here the MRI is essentially irrelevant and therefore the anisotropy feature, i.e., the difference between the dashed line and dotted line in Fig. 1, is only due to the fact that the Jeans mechanism exhibits different growth rates for modes with different χ values. We therefore deduce that magnetic tension is unimportant for sufficiently large scales and we are obtaining an equivalent behavior of the overdensity to the one predicted by Eq. (25) in the limit ω _A → 0 and k → 0. On the other hand, if we focus our attention on small scales, for which the Jeans instability is suppressed, then only the MRI is relevant, and, in this case, the growth rate is again anisotropic (clearly the MRI is strongly suppressed for χ → 0). Nonetheless, this effect is essentially uncoupled from the Jeans instability and the gravitation and magneto-rotational regions of the unstable spectrum simply coexist within the differentially rotating nebula profile.

6. Phenomenological considerations

Our investigation on the co-existence of MRI and Jeans instability and the main conclusion we reach about the noninteraction of the two mechanisms, which however offer a scenario for the generation of anisotropic structures, can be applied in different physical contexts, corresponding to different scales and systems across the Universe. The most natural implementation of the present study is in those astrophysical systems, such as filaments in the interstellar medium and nebula-like structures, which possess a sufficiently high level of ionization and self-gravity to be interpreted via the instability features, predicted by the dispersion relation (25).

Nonetheless, we aim to infer the validity of the present study on a cosmological level, in order to characterize the behavior of small Universe inhomogeneities in the evolution range between the hydrogen recombination (z ≃ 1100) and the formation of the most common structures at large scales (z ∼ 10), where z denotes the red-shift. However, in order to successfully address this characterization of Eq. (25), we must consider some subtle questions concerning the primordial Universe: 1) the Universe background is expanding, i.e., it is nonstationary, and it is also homogeneous in space; 2) only a very weak part of the Universe baryonic component is ionized, about a one part in one hundred thousand; 3) the Universe also contains dark energy (about 70% of its total energy density) and dark matter (about 25% of its energy density).

We now consider each of these questions separately point by point, arguing how the present study is de facto applicable to the early Universe, providing qualitative but reliable information about the dynamics of its inhomogeneities:

1. The nonstationarity of the Universe is relevant for the physics of the early cosmology, including the behavior of linear perturbations, only for spatial scales comparable with the Hubble size (L_H ≃ cH⁻¹(t), where H(t) is the expansion rate and c the speed of light). Such a scale roughly increases as L_H ∼ ct, t being the universal synchronous time, and it is well-known (Kolb & Turner 1990; Montani et al. 2011) that many cosmological scales, relevant for structure formation, become smaller than L_H simply because they increase slower in time, according to the cosmic scale factor behavior a(t) ∼ t^2/3. Therefore, when studying cosmological perturbations whose size is well below the Hubble length L_H, the effects of the expansion can be safely neglected and the background can be regarded as a steady one.

   An important feature introduced by the expansion, with respect to the steady case here considered, consists of a power-law growth in time of the perturbation, against an exponential instability. However, as well-known in the nonmagnetized case (Montani et al. 2011; Weinberg 1972) and also validated in the presence of a magnetic field in Lattanzi et al. (2012 and Pugliese et al. 2012; see also Vlahos et al. 2005), the concept of Jeans threshold scale can be defined in both cases (without or with expansion, respectively). Furthermore, it comes out that the value of the Jeans scale in the two cases differs for a numerical factor only. The presence of a non-stationary expansion therefore does not affect the physical content of our linear analysis, but it could affect the timescale of the considered processes.

   However, our request of a cylindrical symmetry, implied by the presence of differential rotation (we assume the existence of a privileged direction), suggests that our analysis must essentially concern the stability of subregions inside primordial structures, i.e., the background must be considered to be density contrast close to the nonlinear regime, almost bounded and disconnected by the Hubble flux. In this respect, our study can be regarded as relevant for the iteration of the so-called Jeans mechanism for the Universe fragmentation, especially around z ∼ 100, where the density contrast can be greater than unity and the dark matter role is not relevant to the jet (see below point 3).

2. Although the Universe is weakly ionized after the recombination (actually the Universe expansion prevents perfect hydrogen recombination), we can nonetheless reliably argue that the behavior of the ionized component is closely linked to the neutral matter behavior, meaning that our analysis can be applied to the baryonic neutral Universe as well, for a wide range of cosmological scales.

   Ions and neutral atoms actually interact via collisions, mediated by a drag coefficient γ_in ∼ 1.9 × 10⁻⁹ cm³ s ⁻¹. The relevance of this interaction process, commonly dubbed ambipolar diffusion (see, e.g., Sethi & Subramanian 2005; Sethi et al. 2008; Schleicher et al. 2009; Schleicher et al. 2009b), is properly characterized by the ambipolar Reynold number R_amb, defined as follows (Lattanzi et al. 2012) $$\begin{array}{c}\hfill R_{\mathrm{amb}}\equiv L/L_{\mathrm{amb}},L_{\mathrm{amb}}\equiv {\nu }_{\mathrm{A}}^2/\nu {\gamma }_{\mathrm{in}}{n}_{\mathrm{i}},\end{array}$$(37) where n_i is the ion number density and ν and L are a typical velocity and spatial scale of the considered system, respectively. In the present study, ν ∼ ν _s ≳ ν_A, since the MRI typically holds for high values of the plasma β parameter.

   It is well-known (Mestel & Spitzer 1956; Shu 1983; Banerjee & Jedamzik 2004; Li et al. 2006), that when R_amb ≪ 1, the ions and neutral atoms are very weakly interacting and we deal with two distinct components in the system; when R_amb ∼ 1, the two components are coupled and the ambipolar diffusion is relevant, being represented by a dissipative term in MHD (the single fluid representation holds, but a correction to ideal case must be taken into account); finally, for R_amb ≫ 1, that is, L ≫ L_amb, the coupling becomes very strong and the two species are so tightly evolving that the ideal MHD representation can be applied for the whole system, in agreement with the present study.

   It is possible to show (see Fig.1 in Lattanzi et al. 2012) that at the Universe background density for 10 < z < 1100, the Reynold number R_amb, the coupling becomes very strong and the two species are so tightly evolving that the ideal MHD representation can be applied for the whole system, in agreement with the present study.

   It is possible to show (see Fig.1 in Lattanzi et al. 2012) that at the Universe background density for 10 < z < 1100, the Reynold number R_amb remains much larger than unity for spatial scales containing a mass much greater than 10⁶ solar masses, i.e., for the most relevant cosmological scales. Since we intend to apply our analysis to the stability of overdense regions across the Universe, the situation is even slightly better, because L_amb is correspondingly smaller. Therefore, we can safely assume to be in a parameter region where the ambipolar diffusion dissipation term is actually negligible, but ions and neutrals are tightly bound forming a single fluid.

   In order to characterize the cosmological scales for which our analysis is predictive, in Fig. 4 we plot, as function of the red-shift (1 + z), the size associated to the minimal mass of the structure to deal with R_amb > 1 (we define the scale L_c which imply R_amb = 1), which corresponds to the value 3.5 of the ratio y_g ≡ 4πGρ₀/Ω², used in the plots of this paper and reasonable for a cylindrical configuration (cf. with Fig.1 in Lattanzi et al. 2012). The field intensity B₀(z) behaves as (1 + z)² and we set its present value as B₀(z = 0) = 10⁻⁹G. Since we are considering our background to be an overdensity of the expanding Universe, we consider ρ₀ = 3ρ_crit(1 + z)³ (ρ _crit being the Universe critical density): the factor 3 causes the considered region to be almost disconnect from the Hubble flux. We remark that L_c is defined by $L_{\mathrm{c}}^3=M_{\mathrm{c}}{y}_{\mathrm{g}}/4\pi {\rho }_0$, where M_c(z) denotes the critical mass contained within the scale L_amb above which the ambipolar diffusion can be safely neglected. Figure 4 clearly shows that the most relevant cosmological scales for structure formation are in the region where ideal MHD can be used to treat plasma coupled to neutral baryon matter.

   Fig. 4. Critical scale Lc ⇒ Ramb = 1 as function of the redshift for 4πGρ0 = 3.5Ω2. Marked region denotes scales implying Ramb > 1.

3. Despite the fact that dark energy is considered today to form about 70% of the present Universe, there is strong evidence (Peebles & Ratra 2003) that it is surely subdominant, or de facto negligible in any respect, for the Universe evolution at z ≳ 1. Therefore it does not concern the cosmological implementation for the dispersion relation (25). The situation is very different for dark matter, which forms approximately 80% of the total matter in the present Universe. Dark matter interacts with baryonic, neutral, and ionized matter, respectively, only through the gravitational interaction: its inhomogeneities are greater by a factor of about 20 than the baryonic perturbations, and therefore they constitute the gravitational skeleton of the structure formation (Kolb & Turner 1990; Montani et al. 2011). However, since the baryon to photon ratio is very small (and almost constant in time) in the background Universe (its value is ∼6.1 × 10⁻¹⁰), the baryons remain tightly coupled to photons, feeling the radiation pressure well after the hydrogen recombination; up to z ∼ 100 (Lattanzi et al. 2012; Kolb & Turner 1990). As a consequence, the neutral baryons are prevented from falling into the dark matter gravitational potential by the radiation pressure and therefore they are actually decoupled by the dark matter, but strongly coupled to the ionized and photon components of the Universe. In this situation, we are clearly neglecting here the dynamics of dark matter perturbations, but we expect to account for the presence of this component by including its contribution in the Universe mass density we introduced in the Jeans length scale (28). Furthermore, we observe that, even when the radiation pressure is suppressed by the Universe expansion (z < 100 ), the setting down of baryonic matter in the gravitational potential of dark matter requires a finite time, and therefore the decoupling of the two components remains valid roughly up to z ≳ 10. More specifically, it is worth noting that the decoupling between photons and baryons depends on the considered scale, since only structures larger than the photon mean free path can really be coupled to the radiation component: at the recombination age, this scale corresponds to a mass of ∼10¹¹ solar masses. From the point of view of the equation of state, the change regarding the baryon fluid pressure can be interpreted (see Weinberg 1972) as the passage of the polytropic index from 4/3 (baryon density behaves as the inverse of the volume, i.e., like a⁻³, a being the cosmic scale factor, while the radiation pressure goes as a⁻⁴), to the value 5/3, typical of a nonrelativistic fluid. However, to get a quantitative estimate of how the pressure decreases after the recombination and an idea of why it becomes negligible only after z ∼ 100, a kinetic evaluation of the sound speed velocity is necessary (see formula Lattanzi et al. 2012). Before recombination (for a scale greater than the photon mean free path), we can write the following expression for the sound speed. $$\begin{array}{c}\hfill {\nu }_{\mathrm{s}}^2=\frac{c^2}{3}\frac{k_{\mathrm{B}}{T}_{\mathrm{b}}}{m_{\mathrm{p}}{c}^2+1.5\times {10}^9}1.5\times {10}^9,\end{array}$$(38) where m_b and T_b denote the baryon mass and temperature, respectively (k_B being the Boltzmann constant). As long as the temperature of the baryon, coinciding with the photon one, remains sufficiently large, this value is close to that of an ultra-relativistic fluid, that is, ${\nu }_{\mathrm{s}}=c/\sqrt{3}$. After the recombination, the squared sound speed velocity reads $$\begin{array}{c}\hfill {\nu }_{\mathrm{s}}^2=\frac{5}{3}\frac{k_{\mathrm{B}}{T}_{\mathrm{b}}}{m_{\mathrm{b}}},\end{array}$$(39) which is a typical nonrelativistic value, some order of magnitude smaller than the estimate Eq. (38). The crucial point is that, up to z ∼ 100, the baryon and photon temperatures remain essentially equal, while for z < 100, baryons rapidly cool as a⁻² (instead as a⁻¹). Only after this time are baryons free particles, and they start to fall in the potential well of the dark matter. For z > 100, the residual pressure associated with Eq. (39) is responsible for acoustic oscillations of the baryon density and the Jeans scale selects which perturbations increase or oscillate. For z < 100, the sound velocity drastically decreases and the density by which the Jeans scale is calculated increases because the dark matter contribution must be included. As a result, the Jeans threshold scale significantly diminishes, while some perturbations enter the nonlinear dynamical regime.

Even considering the well-posed points above, the implementation of the present investigation of the anisotropic MRI-Jeans instability in a cosmological setting remains valid only on a rather qualitative level. Nonetheless, the obtained results encourage more careful analytical and numerical (maybe N-body) studies to clarify whether or not the formation of filament-like structures across the Universe can be explained with the intrinsic anisotropy of differentially rotating and small magnetized primordial sites.

In this respect, it is worth noting that the presence of the magnetic field introduces a privileged direction in space, which defines an intrinsic anisotropy of the perturbation dynamics, i.e., the angle between the background magnetic field and the perturbation wavenumber enters the dispersion relation. We have already outlined such a property in the case of a homogeneous (cosmological) background; see Lattanzi et al. (2012) and Pugliese et al. (2012). Clearly, the direction of propagation selects different contributions due to perturbed magnetic pressure and tension, respectively. In particular, the perturbed magnetic pressure depends on the angle mentioned above (it depends directly on the angle between B₀ and B₁, but the latter is orthogonal to the wavenumber because it has vanishing divergence) and this makes the response due to the magnetic field contrasting the gravitational force intrinsically anisotropic: perturbations propagating along the background magnetic field only provide tension, and, in that direction, the thermostatic pressure alone prevents gravitational collapse (we can speak of pure acoustic oscillations when stable modes are concerned). However, perturbed magnetic pressure has an intrinsic anisotropy and a positive or negative sign: it can support and contrast the ordinary pressure, altering the value of the Jeans scale considered here (for stable modes, we can speak of fast and slow magneto-acoustic oscillations, respectively).

Nonetheless, it is clear that, in the parameter region where MRI is suppressed in favor of the Jeans mode, the anisotropy due to the magnetic field essentially vanishes and the resulting anisotropic growth of the perturbations is due to the intrinsic anisotropy of the background profile, i.e., a privileged direction exists because a differential rotation of the system is taken along a given axis. In this respect, our background is not really a homogeneous one and we suggest that it must be cosmologically interpreted as a rotating primordial overdensity, on which we are studying the character of the Jeans instability. As already mentioned, the nonsteady behavior of the cosmological background is a minor feature here, because the considered overdensity can be close to the nonlinear regime, almost disconnected from the Universe expansion, which is removed by the internal bound energy, or negligible on a sufficiently small spatial sub-scale.

We conclude by observing that the aim of this section is only to set up the conceptual framework of a more rigorous and expectedly numerical analysis based on a real cosmological background, and discussing in some detail the linear dynamics of the coupled plasma and baryonic fluids. Since the magnetic field takes a very small value, as dictated by the cosmic microwave background radiation constraint B₀ ≲ 10⁻⁹G, and its energy decays like a radiation component, we could, in first approximation, reliably neglect its influence on the background metric of the expanding Universe (however, for an overdense region, the values taken by B₀ can be significantly larger; see Vlahos et al. 2005). Therefore, an interesting upgrade of our study in combination with a realistic cosmological setting could be ensured by the dynamics of a two-fluid system (plasma and neutral baryons), moving on a Tolmann-Bondi nonstationary and radially symmetric background (Montani et al. 2011). The presence of a magnetic field and ambipolar coupling would make this problem a subtle numerical study (maybe analytical under suitable simplifications) of general relativistic MHD. Nonetheless, we are confident that the qualitative scenario inferred here, together with the order of magnitude of our estimates, will survive in such a refinement of the problem. Clearly, when implementing the present scenario within an appropriate cosmological setting, it would be relevant to provide a precise characterization of the role played, after recombination, by the radiation pressure profile in the stability of the different cosmological scales.

7. Conclusions

We have analyzed the local stability of a self-gravitating and differentially rotating magnetized cloud in order to outline how the MHD effects influence the value of the Jeans length in this structure. In other words, we consider an axisymmetric and weakly magnetized system, in which the vertical shear is assumed to be small enough to neglect its stratification. This scheme is the basic paradigm with which the role played by MRI on the growth of gravitational instability is studied and it allows to properly estimate the contributions from magnetic pressure and tension in fixing the precise value of the Jeans critical length.

After a schematic characterization of the background system, we have written down the linear perturbation dynamics and determined the dispersion relation associated to the local linear evolution of the system. We have analyzed specific regimes and then we have outlined the dependence of the dispersion relation on the direction of the perturbation propagation with respect to that of the magnetic field background. We have subsequently studied such anisotropy properties of the dispersion relation concerning the evolution of a spherical overdense region. We have shown that the shape of the perturbation is deeply deformed during its evolution and that it becomes intrinsically anisotropic (up to order of unity), although the anisotropy of the density contrast cannot exceed a few percent, according to the considered linear regime. However, we have clarified how such anisotropy is essentially an intrinsic effect of the Jeans instability since the largest spatial scales, at which it appears to dominate, grow faster in time with respect to the smaller scales where MRI is relevant.

The idea of the Jeans mechanism for the structure fragmentation consists in the decreasing value of the Jeans length as the collapse proceeds and the structure average density increases. The obtained results suggest that such a mechanism is influenced by differential rotation and by the magnetic field too: as the fragmentation of the original background structure evolves, we can infer that the anisotropy degree increasingly affects the substructure shape, generating a filament-like class of subsystems.

The scenario traced above cannot be used to make generalizations regarding structure formation across the primordial Universe. Nonetheless the common presence of both differential rotation and weak magnetic fields in accreting cosmological structures leads us to infer that the generation of filaments from spherical unstable profiles could concern specific cosmological sites or even subclasses of structures (especially if connected in a nonlinear phase with hot dark matter gravitational skeletons) favoring anisotropic collapses.

References

All Figures