1. Introduction
The observed universe is non-homogeneous at different scales. Perturbations of density lead to gravitational collapse of matter and formation of structures of different sizes: stars, galaxies and large-scale structures. Many astrophysical objects are essentially non-spherical, and their formation and evolution can be approximately described as a behaviour of self-gravitating ellipsoidal bodies in the frame of Newtonian gravity.
The theory of ellipsoidal figures was developed in classical works of Newton, Maclaurin, Jacobi, Poincaré and Chandrasekhar. The systematic analysis of ellipsoidal figures of equilibrium is presented in the book of Chandrasekhar (Reference Chandrasekhar1969). Different types of incompressible ellipsoids are investigated there on a basis of virial equations, as well as their dynamical and secular stability.
The shape of an ellipsoid is determined by the sizes of its three semi-axes. Non-rotating figures always take spherical form: all semi-axes are equal. In the presence of rotation a spheroid (axially symmetrical ellipsoid with two equal semi-axes) is formed (Maclaurin spheroid). If the angular momentum of rotation is large enough, development of different kinds of instabilities leads to formation of a three-axis ellipsoid, with rotation around the small axis (Jacobi ellipsoid). An incompressible spheroid can be secularly or dynamically unstable for transition to a three-axis ellipsoid. The points of the onset of the bar-mode dynamical and secular instabilities of the Maclaurin spheroid, for its transition into the Jacobi ellipsoid, are found in the book of Chandrasekhar (Reference Chandrasekhar1969).
An approach based on the ellipsoidal description allows one to investigate the time evolution of a body considering behaviour of its semi-axes. First of all, its initial collapse can be described, which is a central problem in formation of astrophysical objects. In the case of the collapse, compressible ellipsoidal objects are considered. Lynden-Bell (Reference Lynden-Bell1964, Reference Lynden-Bell1965) applied Chandrasekhar’s virial tensor method for the rotating, self-gravitating spheroid of pressure-free gas, and showed the growth of non-axisymmetric perturbations during the collapse. It was also discussed that the slowly shrinking Maclaurin spheroid will enter the Jacobi series if it shrinks slowly enough for the dissipative mechanisms to be operative. Zeldovich (Reference Zeldovich1964) considered the dynamics of rotating dust ellipsoids with linear dependence of velocities on coordinates, when the problem is solved exactly. There are numerical investigations of collapsing pressureless spheroids in papers of Lynden-Bell (Reference Lynden-Bell1964) and Lin et al. (Reference Lin, Mestel and Shu1965). In the paper of Rosensteel and Tran (Reference Rosensteel and Tran1991), the virial equations for rotating Riemann ellipsoids of incompressible fluid are demonstrated to form a Hamiltonian dynamical system. There is a detailed description of the ellipsoid model of rotating stars in papers of Lai et al. (Reference Lai, Rasio and Shapiro1993, Reference Lai, Rasio and Shapiro1994a ,Reference Lai, Rasio and Shapiro b ). Using a variational principle, they derived and investigated equations of the evolution of a compressible Riemann-S ellipsoid, incorporating viscous dissipation and gravitational radiation. Solutions of these approximate equations permitted them to obtain equilibrium models, and to investigate their dynamical and secular stability. The point of instability of compressible Maclaurin spheroids is found in Lai et al. (Reference Lai, Rasio and Shapiro1993) and Shapiro (Reference Shapiro2004). There is a wide-ranging review of this topic in the paper of Lynden-Bell (Reference Lynden-Bell1996) in memory of Chandrasekhar.
In the present paper we review the developed approximate model for consideration of dynamics of collapsing ellipsoidal figures. The present review is based on the papers of Bisnovatyi-Kogan (Reference Bisnovatyi-Kogan2004), Bisnovatyi-Kogan and Tsupko (Reference Bisnovatyi-Kogan and Tsupko2005) and Bisnovatyi-Kogan and Tsupko (Reference Bisnovatyi-Kogan and Tsupko2008a ). We derive and solve equations of the dynamical behaviour in a simple model of a compressible uniformly rotating three-axis ellipsoid. The ordinary differential equations for the ellipsoid semi-axes evolution in time are obtained with help of variation of the Lagrange function of the ellipsoid. In our model a motion along three axes takes place in the gravitational field of a uniform-density ellipsoid, with account of rotation, represented by centrifugal forces, and the isotropic pressure, represented in an approximate non-differential way.
Our model can be used for investigation of star dynamics, and also for the first stages of formation of large-scale structure consisting of dark matter, with appropriate approximations. Starting from § 2, we apply our model for investigation of dynamic stability of non-spherical stars. We show that deviation from the spherical symmetry in a non-rotating star with zero angular momentum leads to formation of a regularly or chaotically oscillating body, in which dynamical motion prevents the formation of the singularity. The non-spherical star without dissipative processes will never reach a singularity in the Newtonian dynamics. Therefore, collapse to a singularity is connected with a secular type of instability, even without rotation. Starting from § 7, we consider an approximate model for large scale structure formation. We estimate emission of very long gravitational waves (GWs) during a collapse of the dark matter, and discuss the possibility of gravitational lensing and polarization of the cosmic microwave background (CMB) by these waves.
Recently, Vandervoort (Reference Vandervoort2011) modelled the galaxy as a heterogeneous ellipsoid with an arbitrary stratification of the density distribution, and represented the mean motions of the stars in terms of a velocity field that sustains that density distribution consistently with the equation of continuity. The virial equations are reduced to a closed system of equations of motion that govern the semi-axes of the characteristic ellipsoid. The semi-axes of the ellipsoid characterizing the size and shape of the galaxy are functions of the time, whereas the stratification of the density does not change. In a second paper (Vandervoort Reference Vandervoort2014), the author described an application of the model to an investigation of chaotic behaviour in axisymmetric pulsations of non-rotating, spherically symmetric configurations and of rotating, axisymmetric configurations. See also papers of Sharif and Zaeem Ul Haq Bhatti (Reference Sharif and Zaeem Ul Haq Bhatti2013, Reference Sharif and Zaeem Ul Haq Bhatti2014) and Rodrigues (Reference Rodrigues2014).
2. Dynamic stability of non-spherical stars
Dynamic stability of spherical stars is determined by an average adiabatic power
$$\begin{eqnarray}{\it\gamma}=\left.\frac{\partial \log P}{\partial \log {\it\rho}}\right|_{S}.\end{eqnarray}$$
For a density distribution
$$\begin{eqnarray}{\it\rho}={\it\rho}_{0}\,{\it\varphi}\left(\frac{m}{M}\right),\end{eqnarray}$$
the star in the Newtonian gravity is stable against dynamical collapse when (Zeldovich and Novikov Reference Zeldovich and Novikov1967; Bisnovatyi-Kogan Reference Bisnovatyi-Kogan1989)
$$\begin{eqnarray}\int _{0}^{R}\left({\it\gamma}-\frac{4}{3}\right)P\frac{\text{d}m}{{\it\varphi}(m/M)}>0.\end{eqnarray}$$
Here
${\it\rho}_{0}$
is a central density,
$M$
is a stellar mass,
$m$
is the mass inside a Lagrangian radius
$r$
, so that
$m=4{\rm\pi}\int _{0}^{r}{\it\rho}r^{2}\,\text{d}r$
,
$M=m(R)$
,
$R$
is a stellar radius and
${\it\varphi}(m/M)$
is the Emden function corresponding to
${\it\gamma}=4/3$
.
Let us consider the equation of state
$P=K{\it\rho}^{{\it\gamma}}$
. A spherical star with such equation of state is stable against collapse at
${\it\gamma}>4/3$
and unstable at
${\it\gamma}<4/3$
, when the instability leads to a collapse to a singularity. At
${\it\gamma}=4/3$
, there is only one equilibrium state for a given mass, at a unique value of the constant
$K_{0}$
. This state is neutral to the stability, and may suffer contraction to a singularity at
$K<K_{0}$
or expansion to infinity at
$K>K_{0}$
.
Collapse of a spherical, initially unstable star may be stopped only by a stiffening of the equation of state, like neutron star formation at late stages of evolution or formation of a fully ionized stellar core with
${\it\gamma}=5/3$
at the collapse of clouds during star formation. Without such stiffening, a spherical star in the Newtonian theory would collapse into a point with an infinite density (singularity).
In the presence of a rotation, a star is becoming more dynamically stable against collapse. Due to the more rapid increase of centrifugal forces during contraction, in comparison with the Newtonian gravitational force, collapse of a rotating star will always be stopped at a finite density by centrifugal forces.
Here we show that deviation from the spherical symmetry in a non-rotating star with zero angular momentum leads to a similar stabilization, and a non-spherical star without dissipative processes will never reach a singularity. Therefore, collapse to a singularity is connected with a secular type of instability, even without rotation.
We calculate the dynamical behaviour of a non-spherical, non-rotating star after its loss of linear stability, and investigate nonlinear stages of contraction. We use an approximate system of dynamical equations, describing three degrees of freedom of a uniform self-gravitating compressible ellipsoidal body (Bisnovatyi-Kogan and Tsupko Reference Bisnovatyi-Kogan and Tsupko2008a
). We obtain that the development of an instability leads to the formation of a regularly or chaotically oscillating body, in which dynamical motion prevents the formation of the singularity. We find regions of chaotic and regular pulsations by constructing a Poincaré diagram for different values of the initial eccentricity and initial entropy. For simplicity, we restrict ourselves by calculating only spheroidal figures with
${\it\gamma}=4/3$
; see also results for
${\it\gamma}=6/5$
in Bisnovatyi-Kogan and Tsupko (Reference Bisnovatyi-Kogan and Tsupko2008a
).
3. Equations of motion for semi-axes of rotating three-axis ellipsoid
Here we derive equations of motion for the general case of a rotating three-axis ellipsoid and consider cases
${\it\gamma}=5/3$
and
${\it\gamma}=4/3$
.
Let us consider a compressible three-axis ellipsoid with semi-axes
$a\neq b\neq c$
$$\begin{eqnarray}\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1.\end{eqnarray}$$
We assume that the ellipsoid rotates uniformly with an angular velocity
${\it\Omega}$
around the axis
$z$
. Let us approximate the density of the matter
${\it\rho}$
in the ellipsoid as uniform. During collapse and posterior time evolution the sizes of the semi-axes
$a$
,
$b$
,
$c$
, angular velocity
${\it\Omega}$
and density
${\it\rho}$
are changed with time, but the ellipsoid keeps an ellipsoidal shape with the homogeneous distribution of density and uniform rotation.
The constant mass
$m$
and total angular momentum
$M$
of the uniform ellipsoid are connected with uniform density, angular velocity and semi-axes as (
$V$
is the volume of the ellipsoid)
$$\begin{eqnarray}m={\it\rho}\,V=\frac{4{\rm\pi}}{3}\,{\it\rho}\,abc,\quad M={\it\rho}\,{\it\Omega}\int _{V}(x^{2}+y^{2})\,\text{d}V=\frac{m}{5}\,{\it\Omega}(a^{2}+b^{2}).\end{eqnarray}$$
We assume a linear dependence of the velocities on the coordinates in the rotating frame:
$$\begin{eqnarray}{\it\upsilon}_{x}=\frac{{\dot{a}}x}{a},\quad {\it\upsilon}_{y}=\frac{{\dot{b}}y}{b},\quad {\it\upsilon}_{z}=\frac{{\dot{c}}z}{c}.\end{eqnarray}$$
If an homogeneous ellipsoidal body of inviscid fluid initially has velocities that are linear functions of position, then it will remain ellipsoidal and the internal velocities will remain linear functions of position for all time provided that the density remains spatially uniform (see, for example, Lynden-Bell Reference Lynden-Bell1996).
In absence of any dissipation, this ellipsoid is a conservative system. To derive the equations of motion (equations describing time evolutions of semi-axes), let us write a Lagrange function of the ellipsoid, which consists of kinetic energy
$U_{kin}$
and potential energy
$U_{pot}$
of the ellipsoid:
$$\begin{eqnarray}L=U_{kin}-U_{pot},\quad U_{pot}=U_{g}+E_{th}+U_{rot}.\end{eqnarray}$$
The kinetic energy is written as (using (3.3))
$$\begin{eqnarray}U_{kin}=\frac{1}{2}\,{\it\rho}\int _{V}({\it\upsilon}_{x}^{2}+{\it\upsilon}_{y}^{2}+{\it\upsilon}_{z}^{2})\,\text{d}V=\frac{m}{10}\,({\dot{a}}^{2}+{\dot{b}}^{2}+{\dot{c}}^{2}).\end{eqnarray}$$
The gravitational energy
$U_{g}$
of the uniform ellipsoid is defined as (Landau and Lifshitz Reference Landau and Lifshitz1993)
$$\begin{eqnarray}U_{g}=-\frac{3Gm^{2}}{10}\int _{0}^{\infty }\frac{\text{d}u}{\sqrt{(a^{2}+u)(b^{2}+u)(c^{2}+u)}}.\end{eqnarray}$$
In case of spheroids (
$a=b\neq c$
), this integral can be taken analytically; in case of ellipsoids it is expressed in elliptical integrals (see Bisnovatyi-Kogan Reference Bisnovatyi-Kogan2004; Bisnovatyi-Kogan and Tsupko Reference Bisnovatyi-Kogan and Tsupko2005).
The rotational energy is
$$\begin{eqnarray}U_{rot}=\frac{1}{2}\,{\it\rho}\int _{V}v_{rot}^{2}\,\text{d}V=\frac{1}{2}\,{\it\rho}\,{\it\Omega}^{2}\int _{V}(x^{2}+y^{2})\,\text{d}x\,\text{d}y\,\text{d}z=\frac{m}{10}\,{\it\Omega}^{2}(a^{2}+b^{2}).\end{eqnarray}$$
The angular velocity
${\it\Omega}$
is changed with time, and it is more convenient to express (3.7) via
$M$
, which is constant in absence of dissipation. Taking into account the second expression in (3.2), we obtain the relation
$$\begin{eqnarray}U_{rot}=\frac{5}{2}\,\frac{M^{2}}{m(a^{2}+b^{2})}.\end{eqnarray}$$
Different forms of thermal energy
$E_{th}$
are used depending on problem under consideration.
3.1. Rotating ellipsoid, case
${\it\gamma}=5/3$
First, we use the non-relativistic equation of state
$P=K{\it\rho}^{{\it\gamma}}$
with polytropic index
${\it\gamma}=5/3$
. It corresponds to the case of matter consisting of non-collisional non-relativistic dark matter particles.
The total thermal energy, with adiabatic index
$n=1.5~({\it\gamma}=5/3)$
, of non-relativistic dark matter particles in the ellipsoid is
$E_{th}\sim V^{-2/3}\sim (abc)^{-2/3}$
, and the relation between pressure
$P$
and thermal energy
$E_{th}$
is
$E_{th}=(3/2)PV$
. Let us write
$E_{th}$
using initial values
$E_{th,in}$
,
$a_{in}$
,
$b_{in}$
and
$c_{in}$
as
$$\begin{eqnarray}E_{th}=\frac{E_{th,in}(a_{in}b_{in}c_{in})^{2/3}}{(abc)^{2/3}}=\frac{{\it\varepsilon}}{(abc)^{2/3}}.\end{eqnarray}$$
We introduce here
${\it\varepsilon}$
, which we will call the entropy function. The entropy function
${\it\varepsilon}$
is constant in the conservative case, but it is variable in presence of dissipation.
By variation of the Lagrange function, we obtain the Lagrange equations of motion:
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=-\frac{3Gm}{2}\,a\int _{0}^{\infty }\frac{\text{d}u}{(a^{2}+u){\it\Delta}}+\,\frac{10}{3m}\,\frac{1}{a}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}+\frac{25M^{2}}{m^{2}}\,\frac{a}{(a^{2}+b^{2})^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{b}=-\frac{3Gm}{2}\,b\int _{0}^{\infty }\frac{\text{d}u}{(b^{2}+u){\it\Delta}}+\,\frac{10}{3m}\,\frac{1}{b}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}+\frac{25M^{2}}{m^{2}}\,\frac{b}{(a^{2}+b^{2})^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\ddot{c}=-\frac{3Gm}{2}\,c\int _{0}^{\infty }\frac{\text{d}u}{(c^{2}+u){\it\Delta}}+\,\frac{10}{3m}\,\frac{1}{c}\,\frac{{\it\varepsilon}}{(abc)^{2/3}},\quad {\it\Delta}^{2}=(a^{2}+u)(b^{2}+u)(c^{2}+u).\end{eqnarray}$$
At the right-hand side of these equations there are forces per unit ellipsoid mass. The first term is the gravitational force directed to the centre of the ellipsoid, the second term is the force of thermal pressure directed outward and the third term is the centrifugal force directed outward.
It is easy to check that equilibrium solutions of these equations are the Maclaurin spheroid and the Jacobi ellipsoid (Chandrasekhar Reference Chandrasekhar1969).
3.2. Non-rotating ellipsoid,
${\it\gamma}=4/3$
For investigation of dynamic stability of non-spherical bodies we consider a non-rotating compressible homogeneous ellipsoid with the equation of state
$P=K{\it\rho}^{{\it\gamma}}$
,
${\it\gamma}=4/3$
. The case
${\it\gamma}=5/3$
considered above is not interesting for investigation of dynamic stability, because an isentropic spherical star with
${\it\gamma}=5/3$
always stops contraction, and never suffers collapse to a singularity.
A spherical star with
${\it\gamma}=4/3$
collapses to a singularity at small enough
$K$
, and we show here how deviations from a spherical form prevent formation of any singularity. For
${\it\gamma}=4/3$
, the thermal energy of the ellipsoid is
$E_{th}\sim V^{-1/3}\sim (abc)^{-1/3}$
, and the value
$$\begin{eqnarray}{\it\varepsilon}=E_{th}(abc)^{1/3}=3\left(\frac{3m}{4{\rm\pi}}\right)^{1/3}K\end{eqnarray}$$
remains constant in time. A Lagrange function of the ellipsoid is written as
$$\begin{eqnarray}L=U_{kin}-U_{pot},\quad U_{pot}=U_{g}+E_{th},\quad E_{th}=\frac{{\it\varepsilon}}{(abc)^{1/3}}.\end{eqnarray}$$
Equations of motion describing the behaviour of the three semi-axes
$(a,b,c)$
are obtained from the Lagrange function (3.14a−c
) in the form
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=-\frac{3Gm}{2}\,a\int _{0}^{\infty }\frac{\text{d}u}{(a^{2}+u){\it\Delta}}+\,\frac{5}{3m}\,\frac{1}{a}\,\frac{{\it\varepsilon}}{(abc)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{b}=-\frac{3Gm}{2}\,b\int _{0}^{\infty }\frac{\text{d}u}{(b^{2}+u){\it\Delta}}+\,\frac{5}{3m}\,\frac{1}{b}\,\frac{{\it\varepsilon}}{(abc)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{c}=-\frac{3Gm}{2}\,c\int _{0}^{\infty }\frac{\text{d}u}{(c^{2}+u){\it\Delta}}+\,\frac{5}{3m}\,\frac{1}{c}\,\frac{{\it\varepsilon}}{(abc)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\Delta}^{2}=(a^{2}+u)(b^{2}+u)(c^{2}+u). & \displaystyle\end{eqnarray}$$
4. Dynamic stabilization of non-spherical bodies against unlimited collapse, numerical results
To obtain a numerical solution of (3.15)–(3.17), we write them in non-dimensional variables. Let us introduce the variables
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{t}=\frac{t}{t_{0}},\quad \tilde{a}=\frac{a}{a_{0}},\quad \tilde{b}=\frac{b}{a_{0}},\quad \tilde{c}=\frac{c}{a_{0}},\\ \displaystyle \tilde{m}=\frac{m}{m_{0}},\quad \tilde{{\it\rho}}=\frac{{\it\rho}}{{\it\rho}_{0}},\quad \tilde{U} =\frac{U}{U_{0}},\quad \tilde{E}_{th}=\frac{E_{th}}{U_{0}},\quad \tilde{{\it\varepsilon}}=\frac{{\it\varepsilon}}{{\it\varepsilon}_{0}}.\end{array}\!\right\}\end{eqnarray}$$
The scaling parameters
$t_{0},a_{0},m_{0},{\it\rho}_{0},U_{0},{\it\varepsilon}_{0}$
are connected by the following relations:
$$\begin{eqnarray}t_{0}^{2}=\frac{a_{0}^{3}}{Gm_{0}},\quad U_{0}=\frac{Gm_{0}^{2}}{a_{0}},\quad {\it\rho}_{0}=\frac{m_{0}}{a_{0}^{3}},\quad {\it\varepsilon}_{0}=U_{0}a_{0}.\end{eqnarray}$$
We have a system of non-dimensional equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=-\frac{3m}{2}\,a\int _{0}^{\infty }\frac{\text{d}u}{(a^{2}+u){\it\Delta}}+\,\frac{5}{3m}\,\frac{1}{a}\,\frac{{\it\varepsilon}}{(abc)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{b}=-\frac{3m}{2}\,b\int _{0}^{\infty }\frac{\text{d}u}{(b^{2}+u){\it\Delta}}+\,\frac{5}{3m}\,\frac{1}{b}\,\frac{{\it\varepsilon}}{(abc)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{c}=-\frac{3m}{2}\,c\int _{0}^{\infty }\frac{\text{d}u}{(c^{2}+u){\it\Delta}}+\,\frac{5}{3m}\,\frac{1}{c}\,\frac{{\it\varepsilon}}{(abc)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\Delta}^{2}=(a^{2}+u)(b^{2}+u)(c^{2}+u). & \displaystyle\end{eqnarray}$$
In (4.3)–(4.5), only non-dimensional variables are used, and the ‘tilde’ sign is omitted for simplicity in this section. The non-dimensional Hamiltonian (or non-dimensional total energy) is
$$\begin{eqnarray}\displaystyle H & = & \displaystyle U_{kin}+U_{g}+E_{th}=\frac{m}{10}\,({\dot{a}}^{2}+{\dot{b}}^{2}+{\dot{c}}^{2})\nonumber\\ \displaystyle & & \displaystyle -\,\frac{3m^{2}}{10}\int _{0}^{\infty }\frac{\text{d}u}{\sqrt{(a^{2}+u)(b^{2}+u)(c^{2}+u)}}+\frac{{\it\varepsilon}}{(abc)^{1/3}}.\end{eqnarray}$$
In case of the sphere (
$a=b=c$
,
${\dot{a}}={\dot{b}}={\dot{c}}$
), the non-dimensional Hamiltonian and the non-dimensional equations of motion reduce to
$$\begin{eqnarray}\displaystyle & \displaystyle H=\frac{3}{10}\,m{\dot{a}}^{2}-\frac{3}{5a}\left(m^{2}-\frac{5}{3}\,{\it\varepsilon}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=-\frac{1}{ma^{2}}\left(m^{2}-\frac{5}{3}\,{\it\varepsilon}\right). & \displaystyle\end{eqnarray}$$
${\it\varepsilon}$
$$\begin{eqnarray}{\it\varepsilon}_{eq}=\frac{3m^{2}}{5},\end{eqnarray}$$
at which the spherical star has zero total energy, and it may have an arbitrary radius. For smaller
${\it\varepsilon}<{\it\varepsilon}_{eq}$
, the sphere should contract to a singularity and, for
${\it\varepsilon}>{\it\varepsilon}_{eq}$
, there will be a total disruption of the star with an expansion to infinity. We solve here numerically the equations of motion for a spheroid with
$a=b\neq c$
, which, using (4.3)–(4.5), are written for the oblate spheroid with
$k=c/a<1$
as
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=\frac{3}{2}\frac{m}{a^{2}(1-k^{2})}\left[k-\frac{\arccos k}{\sqrt{1-k^{2}}}\right]+\,\frac{5}{3m}\,\frac{1}{a}\,\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{c}=-3\frac{m}{a^{2}(1-k^{2})}\left[1-\frac{k\arccos k}{\sqrt{1-k^{2}}}\right]+\frac{5}{3m}\,\frac{1}{c}\,\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}}; & \displaystyle\end{eqnarray}$$
$k=c/a>1$
as 
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=-\frac{3}{2}\frac{m}{a^{2}(k^{2}-1)}\left[k-\frac{\cosh ^{-1}k}{\sqrt{k^{2}-1}}\right]+\frac{5}{3m}\,\frac{1}{a}\,\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{c}=3\frac{m}{a^{2}(k^{2}-1)}\left[1-\frac{k\cosh ^{-1}k}{\sqrt{k^{2}-1}}\right]+\frac{5}{3m}\,\frac{1}{c}\,\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}{\it\varepsilon}_{\ast }=\frac{5}{3}\frac{{\it\varepsilon}}{m^{2}},\quad t_{\ast }=t\sqrt{m}.\end{eqnarray}$$
In these variables, (4.11)–(4.14) and (4.9) are written as (omitting the subscript ‘*’)
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=\frac{3}{2a^{2}(1-k^{2})}\left[k-\frac{\arccos k}{\sqrt{1-k^{2}}}\right]+\,\frac{1}{a}\,\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{c}=-\frac{3}{a^{2}(1-k^{2})}\left[1-\frac{k\arccos k}{\sqrt{1-k^{2}}}\right]+\frac{1}{c}\,\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}} & \displaystyle\end{eqnarray}$$
$k=c/a<1$
, 
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=-\frac{3}{2a^{2}(k^{2}-1)}\left[k-\frac{\cosh ^{-1}k}{\sqrt{k^{2}-1}}\right]+\frac{1}{a}\,\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{c}=\frac{3}{a^{2}(k^{2}-1)}\left[1-\frac{k\cosh ^{-1}k}{\sqrt{k^{2}-1}}\right]+\frac{1}{c}\,\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}} & \displaystyle\end{eqnarray}$$
$k=c/a>1$
and 
$$\begin{eqnarray}\ddot{a}=-\frac{1-{\it\varepsilon}}{a^{2}}\end{eqnarray}$$
for the sphere, where the equilibrium corresponds to
${\it\varepsilon}_{eq}=1$
. Near the spherical shape we should use expansions around
$k=1$
, which leads to equations of motion valid for both oblate and prolate cases:
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=-\frac{1-{\it\varepsilon}}{a^{2}}+\left(\frac{{\it\varepsilon}}{3}+\frac{3}{5}\right)\frac{1-k}{a^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{c}=-\frac{1-{\it\varepsilon}}{a^{2}}+\left(\frac{4{\it\varepsilon}}{3}-\frac{4}{5}\right)\frac{1-k}{a^{2}}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}H_{\ast }=\frac{H}{m^{2}},\end{eqnarray}$$
and omitting ‘*’ we have
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle H=\frac{{\dot{a}}^{2}}{5}+\frac{{\dot{c}}^{2}}{10}-\frac{3}{5a}\frac{\arccos k}{\sqrt{1-k^{2}}}+\frac{3}{5}\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}},\quad (\text{oblate})\\ \displaystyle H=\frac{{\dot{a}}^{2}}{5}+\frac{{\dot{c}}^{2}}{10}-\frac{3}{5a}\frac{\cosh ^{-1}k}{\sqrt{k^{2}-1}}+\frac{3}{5}\frac{{\it\varepsilon}}{(a^{2}c)^{1/3}},\quad (\text{prolate})\\ \displaystyle H=\frac{3{\dot{a}}^{2}}{10}-\frac{3}{5a}(1-{\it\varepsilon}),\quad (\text{sphere})\\ \displaystyle H=\frac{{\dot{a}}^{2}}{5}+\frac{{\dot{c}}^{2}}{10}-\frac{3}{5a}\left(1+\frac{{\it\delta}}{3}+\frac{2{\it\delta}^{2}}{15}\right)+\frac{3{\it\varepsilon}}{5a}\left(1+\frac{{\it\delta}}{3}+\frac{2{\it\delta}^{2}}{9}\right),\\ {\it\delta}=1-k,\quad (\text{around the sphere}),~|{\it\delta}|\ll 1.\end{array}\!\right\}\end{eqnarray}$$
The solution of the system of equations (4.11)–(4.14) was performed for initial conditions at
$t=0$
:
${\dot{c}}_{0}=0$
, different values of initial
$a_{0},{\dot{a}}_{0},k_{0}$
and different values of the constant parameter
${\it\varepsilon}$
. The results of our numerical calculations are:
-
(i) sphere (
$a=b=c$
):
${\it\varepsilon}=1$
– total energy equals zero (
$H=0$
), radius is arbitrary,
${\it\varepsilon}<1$
– spherical star collapses to a singularity,
${\it\varepsilon}>1$
– disruption of star with expansion to infinity; -
(ii) spheroid (
$a=b\neq c$
):at
${\it\varepsilon}\geqslant 1$
the total energy is
$H\geqslant 0$
– disruption of star with expansion to infinity,at
${\it\varepsilon}<1$
the total energy is
$H<0$
– the oscillatory regime is established.
In the oscillatory regime the dynamical motion prevents the formation of the singularity. The type of oscillatory regime depends on initial conditions, and may be represented either by regular periodic oscillations or by chaotic behaviour. Examples of two types of such oscillations are represented in figures 1 and 2 (regular and periodic) and in figure 3 (chaotic). For rigorous separation between these kinds of oscillations, we use a method developed by Poincaré (Lichtenberg and Lieberman Reference Lichtenberg and Lieberman1983). The same approach is used in recent papers of Vandervoort (Reference Vandervoort2011, Reference Vandervoort2014).
Example of regular motion of spheroid with
${\it\gamma}=4/3$
,
$H=-1/5$
and
${\it\varepsilon}=2/3$
. This motion corresponds to the full line in the Poincaré map in figure 4.
Example of regular motion of spheroid with
${\it\gamma}=4/3$
,
$H=-1/5$
and
${\it\varepsilon}=2/3$
. This motion corresponds to the point inside the regular region in the Poincaré map in figure 4.
Example of chaotic motion of spheroid with
${\it\gamma}=4/3$
,
$H=-1/5$
and
${\it\varepsilon}=2/3$
. This motion corresponds to grey points in the Poincaré map in figure 5.
5. Investigation of regular and chaotic oscillations via the Poincaré sections
To investigate regular and chaotic dynamics, we use the method of Poincaré section (Lichtenberg and Lieberman Reference Lichtenberg and Lieberman1983), and obtain the Poincaré map for different values of the total energy
$H$
. Let us consider a spheroid with semi-axes
$a=b\neq c$
. This system has two degrees of freedom. Therefore, in this case the phase space is four dimensional:
$a$
,
${\dot{a}}$
,
$c$
,
${\dot{c}}$
. If we choose, at given
${\it\varepsilon}$
, a value of the Hamiltonian
$H_{0}$
, we fix a three-dimensional energy surface
$H(a,{\dot{a}},c,{\dot{c}})=H_{0}$
. During the integration of (4.11)–(4.14), which preserve the constant
$H$
, we fix moments
$t_{i}$
when
${\dot{c}}=0$
. At these moments there are only two independent values (i.e.
$a$
and
${\dot{a}}$
), because the value of
$c$
is determined uniquely from the relation for the Hamiltonian at constant
$H$
.
For the same values of
$H$
and
${\it\varepsilon}$
, we solve the equations of motion (4.11)–(4.14) at initial
${\dot{c}}=0$
and different
$a,\,\,{\dot{a}}$
. For each integration we put the points in the plane
$(a,{\dot{a}})$
at the moments
$t_{i}$
. These points are the intersection points of the trajectories on the three-dimensional energy surface with a two-dimensional plane
${\dot{c}}=0$
, called a Poincaré section.
For each fixed combination of
${\it\varepsilon},H$
, we get the Poincaré map, represented in figures 4–7. The condition
${\dot{c}}=0$
splits in two cases, of a minimum and of a maximum of
$c$
. The Poincaré maps are drawn separately, either for the minimum or for the maximum of
$c$
, and both maps lead to identical results. The regular oscillations are represented by closed lines in the Poincaré map, and chaotic behaviour fills regions of a finite square with dots. These regions are separated from the regions of the regular oscillations by a separatrix line.
The Poincaré map for five regular and two chaotic trajectories in case of
${\it\gamma}=4/3$
,
$H=-1/5$
and
${\it\varepsilon}=2/3$
. The
$(a,{\dot{a}})$
values are taken in the minimum of
$c$
. Full black line is the bounding curve. The point inside the regular region corresponds to coherent oscillations with the same period for
$a$
and
$c$
values, represented in figure 2.
Zoom of previous figure.
The Poincaré map for two chaotic trajectories in case of
${\it\gamma}=4/3$
,
$H=-1/2$
and
${\it\varepsilon}=1/6$
. The
$(a,{\dot{a}})$
values are taken in the minimum of
$c$
. Full black line is the bounding curve.
The Poincaré map for six regular trajectories in the case of
${\it\gamma}=4/3$
,
$H=-3/50$
and
${\it\varepsilon}=9/10$
. The
$(a,{\dot{a}})$
values are taken in the minimum of
$c$
. Full black line is the bounding curve. The point inside the regular region corresponds to coherent oscillations with the same period for
$a$
and
$c$
values, similar to those represented in figure 2.
6. Dynamic stabilization of non-spherical bodies against unlimited collapse, discussion
The main result following from our calculations is the indication of a degenerate nature of formation of a singularity in unstable Newtonian self-gravitating gaseous bodies. Only pure spherical models can collapse to a singularity; any kind of non-sphericity leads to nonlinear stabilization of the collapse by a dynamic motion, and formation of a regularly or chaotically oscillating body. This conclusion is valid for all unstable equations of state, namely, for adiabatic ones with
${\it\gamma}<4/3$
. In addition to the case with
${\it\gamma}=4/3$
, we have calculated the dynamics of the model with
${\it\gamma}=6/5$
, and have obtained similar results.
Note that the region of chaotic behaviour in the Poincaré map is gradually increasing for
${\it\gamma}=4/3$
with decreasing of the entropy
${\it\varepsilon}$
and the total energy
$H$
. At
${\it\varepsilon}=9/10$
and
$H=-3/50$
, we have found only regular oscillations (figure 7), at
${\it\varepsilon}=2/3$
and
$H=-1/5$
both kinds of oscillations are present (figures 4 and 5) and only chaotic behaviour is found at
${\it\varepsilon}=1/6$
and
$H=-1/2$
(figure 6). We connect (Bisnovatyi-Kogan and Tsupko Reference Bisnovatyi-Kogan and Tsupko2005) this chaotic behaviour with development of an anisotropic instability, when radial velocities strongly exceed the transversal ones (Antonov Reference Antonov1973; Fridman and Polyachenko Reference Fridman and Polyachenko1985).
In reality, the presence of dissipation leads to damping of these oscillations, and to final collapse of the non-rotating model, when the total energy of the body is negative. In the case of a core-collapse supernova the main dissipation is due to emission of neutrinos. The time of the neutrino losses is much larger than the characteristic time of the collapse, so we may expect that the collapse leads to formation of a neutron star where non-spherical modes are excited and exist during several seconds after the collapse. In addition to the damping due to neutrino emission, shock waves will be generated, determining highly variable energy losses during the oscillations. Besides viscosity and radiation, which can damp the ellipsoid-like motions and allow collapse, there is the possibility that inertial interactions with higher-order modes that must be present in the real stratified bodies may cause an inertial cascade and drain the energy from the second-order modes in a non-secular way that does not depend on dissipative coefficients. It is very seductive to connect chaotic oscillations with highly variable emission observed in the prompt gamma ray emission of cosmic gamma ray bursts (Mazets et al. Reference Mazets, Golenetskii, Il’inskii, Panov, Aptekar, Gur’yan, Proskura, Sokolov, Sokolova and Kharitonova1981). The presence of rotation and magnetic field strongly complicates the picture of the core-collapse supernova explosion (Ardeljan et al. Reference Ardeljan, Bisnovatyi-Kogan and Moiseenko2005; Moiseenko et al. Reference Moiseenko, Bisnovatyi-Kogan and Ardeljan2006).
In this paper we consider in detail only spheroidal bodies. In reality, the spheroid will become a triaxial ellipsoid during the motion. In addition to spheroids we calculated many variants with triaxial figures (see Bisnovatyi-Kogan and Tsupko Reference Bisnovatyi-Kogan and Tsupko2005). Qualitatively we obtain the same results for ellipsoids: no singularity was reached for any
${\it\varepsilon}>0$
, and establishing an oscillatory (regular or chaotic) regime under negative total energy prevents the collapse. However, in the case of the ellipsoid with semi-axes
$a\neq b\neq c$
we have a system with three degrees of freedom and six-dimensional phase space. Therefore, we could not carry out a rigorous investigation of regular and chaotic types of motion by constructing the Poincaré map, as was done for a spheroid with two degrees of freedom, and restricted ourselves by a description of the spheroidal case.
In the frame of general relativity, dynamic stabilization against collapse by nonlinear non-spherical oscillations cannot be universal. When the size of the body approaches the gravitational radius, no stabilization is possible at any
${\it\gamma}$
. Nevertheless, the nonlinear stabilization may happen at larger radii, so after damping of the oscillations the star would collapse to a black hole. Due to development of non-spherical oscillations there is a possibility for emission of GWs during the collapse of non-rotating bodies with an intensity similar to rotating bodies, or even larger.
7. Zeldovich’s pancakes and large-scale structure of universe
One of the most important problems in modern cosmology is the formation of a large-scale structure of the universe. Observations of large-scale structures of the universe demonstrate strong inhomogeneity in the density distribution. Zeldovich (Reference Zeldovich1970a ,Reference Zeldovich b ) has investigated the fragmentation of a homogeneous medium under the action of gravitation. It has been shown that an arbitrary perturbation in a medium with zero pressure leads to formation of strongly non-spherical structures (discs), and only pure spherical perturbations lead to collapse into a point. These discs are now known as ‘Zeldovich’s pancakes’.
Depending on the medium properties, two types of pancakes can be formed.
-
(i) During contraction of the baryonic gas the collapse is stopped by the pressure; a shock wave is formed. Behind the shock front the cooling, fragmentation and formation of vortices occur, which is necessary for formation of spiral galaxies with large angular momentum.
-
(ii) During collapse of non-collisional dark matter the pancake is also formed, but no shocks are forming at the stage of maximal contraction. The medium is non-collisional, so, after the initial collapse, subsequent oscillations occur. The particles fly through each other many times, and mixing in the phase space occurs – so-called violent relaxation. The system reaches a quasi-equilibrium state via violent relaxation (Lynden-Bell Reference Lynden-Bell1967) and phase mixing during the most rapid dynamic stages. The established equilibrium state has a form depending on the angular momentum of rotation, and size depending on the final values of energy and entropy.
Subsequent nonlinear evolution and interaction between pancakes lead to formation of a very complicated large-scale structure of matter in the universe observed in the sky. It was shown by numerical simulations in two-dimensional space by Doroshkevich et al. (Reference Doroshkevich, Kotok, Polyudov, Shandarin, Sigov and Novikov1980). The first three-dimensional simulations were performed by Klypin and Shandarin (Reference Klypin and Shandarin1983), demonstrating the origin of filaments in the dark matter density field. For a review, see Shandarin and Zeldovich (Reference Shandarin and Zeldovich1989), where the developments of the original Zeldovich’s idea were summarized. The development of perturbations at the nonlinear stage and structure formation in the medium with baryonic gas and non-collisional cold dark matter were investigated by using numerical modelling by many groups (see for example Boily et al. Reference Boily, Clarke and Murray1999; Doroshkevich et al. Reference Doroshkevich, Müller, Retzlaff and Turchaninov1999; Sheth et al. Reference Sheth, Mo and Tormen2001; Boily et al. Reference Boily, Athanassoula and Kroupa2002). For numerical modelling, N-body simulations are used, which are very time consuming.
The first stage of large scale structure formation is a collapse of non-collisional dark matter with formation of a pancake and subsequent relaxation with formation of quasi-equilibrium structures. This stage can be investigated using an approximate model of large scale structure formation considered in works of Bisnovatyi-Kogan (Reference Bisnovatyi-Kogan2004) and Bisnovatyi-Kogan and Tsupko (Reference Bisnovatyi-Kogan and Tsupko2005). Forms of the collapsing matter have been approximated as a rotating homogeneous ellipsoid with two (Bisnovatyi-Kogan Reference Bisnovatyi-Kogan2004) and three (Bisnovatyi-Kogan and Tsupko Reference Bisnovatyi-Kogan and Tsupko2005) degrees of freedom. A correct description of pressure effects, attained by such an approach, and addition of a relaxation permit us to get the dynamics of motion without any numerical singularities. In our model a motion along three axes takes place in the gravitational field of a uniform-density ellipsoid, with account of the isotropic pressure, represented in an approximate non-differential way. The relaxation leads to a transformation of the kinetic energy of the ordered motion into the kinetic energy of the chaotic motion, and to an increase of the effective pressure and thermal energy. All losses are connected with runaway particles. The collapse of the rotating three-axis ellipsoid is approximated by a system of ordinary differential equations, where the relaxation and the losses of energy, mass and angular momentum are taken into account phenomenologically (Bisnovatyi-Kogan Reference Bisnovatyi-Kogan2004). The system is solved numerically for several parameters, characterizing the configuration.
8. Equations of motion with dissipation and numerical results
The equations of motion of a rotating dark matter ellipsoid without any dissipation are derived in § 3.1; see (3.10)–(3.12).
In reality, there is relaxation in the collisionless system, connected with phase mixing, which is called ‘violent relaxation’ (Lynden-Bell Reference Lynden-Bell1967). This relaxation leads to a transformation of the kinetic energy of the ordered motion into the kinetic energy of the chaotic motion and increases effective pressure and thermal energy.
The main transport process is an effective bulk viscosity. Therefore, there is a drag force, which is described phenomenologically by adding of the terms
$$\begin{eqnarray}-\!\frac{{\dot{a}}}{{\it\tau}_{rel}},\quad -\frac{{\dot{b}}}{{\it\tau}_{rel}},\quad -\frac{{\dot{c}}}{{\it\tau}_{rel}}\end{eqnarray}$$
in the right-hand parts of the equations of motion.
Dissipation leads to heat production. The rate of heat production leading to the growth of the entropy of the spheroidal matter was found in Bisnovatyi-Kogan (Reference Bisnovatyi-Kogan2004) from the equations of motion (3.10)–(3.12) with account of (3.4), allowing for the energy, mass and angular momentum losses.
Here we have scaled the relaxation time
${\it\tau}_{rel}$
by the Jeans characteristic time
$$\begin{eqnarray}{\it\tau}_{J}=\frac{2{\rm\pi}}{{\it\omega}_{J}}=\frac{2{\rm\pi}}{\sqrt{4{\rm\pi}G{\it\rho}}}=2{\rm\pi}\sqrt{\frac{abc}{3Gm}}\end{eqnarray}$$
with a constant value of
${\it\alpha}_{rel}$
$$\begin{eqnarray}{\it\tau}_{rel}=2{\rm\pi}\,{\it\alpha}_{rel}\sqrt{\frac{abc}{3Gm}}.\end{eqnarray}$$
The energy balance relation follows from the dynamic equations in the form
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}U_{tot}}{\text{d}t}=\frac{\partial U_{g}}{\partial m}\frac{\text{d}m}{\text{d}t}+\frac{\partial U_{rot}}{\partial m}\frac{\text{d}m}{\text{d}t}+\frac{\partial U_{rot}}{\partial M}\frac{\text{d}M}{\text{d}t}+\frac{\partial E_{th}}{\partial {\it\varepsilon}}\frac{\text{d}{\it\varepsilon}}{\text{d}t}-2\frac{U_{kin}}{{\it\tau}_{rel}}-\frac{U_{kin}}{m}\frac{\text{d}m}{\text{d}t}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & U_{tot}=U_{kin}+U_{g}+U_{rot}+E_{th}. & \displaystyle\end{eqnarray}$$
$U_{kin}$
; the characteristic times for mass, angular momentum and energy losses
${\it\tau}_{ml}$
,
${\it\tau}_{Ml}$
and
${\it\tau}_{el}$
are considerably greater than
${\it\tau}_{rel}$
. The equations describing different losses may be phenomenologically written as 
$$\begin{eqnarray}\frac{\text{d}U_{tot}}{\text{d}t}=-\frac{U_{kin}}{{\it\tau}_{ml}},\quad \frac{\text{d}m}{\text{d}t}=\frac{mU_{kin}}{U_{g}{\it\tau}_{ml}},\quad U_{g}<0,\quad \frac{M}{M_{in}}=\left(\frac{m}{m_{in}}\right)^{{\it\alpha}_{ml}/{\it\alpha}_{Ml}},\end{eqnarray}$$
where
$m_{in}$
and
$M_{in}$
are the initial values of the corresponding parameters. Scaling all the characteristic times by the Jeans value, we have
${\it\tau}_{el}={\it\alpha}_{el}{\it\tau}_{J},{\it\tau}_{ml}={\it\alpha}_{ml}{\it\tau}_{J},{\it\tau}_{Ml}={\it\alpha}_{Ml}{\it\tau}_{J},$
with constant values of
${\it\alpha}_{i}\,(i=el,ml,Ml)$
. The function
${\it\varepsilon}$
determines the entropy of the matter in the spheroid, so that, using (3.9), we have
$$\begin{eqnarray}\frac{\text{d}E_{th}}{\text{d}t}+P\frac{\text{d}V}{\text{d}t}=\frac{1}{(abc)^{2/3}}\frac{\text{d}{\it\varepsilon}}{\text{d}t},\quad E_{th}=\frac{{\it\varepsilon}}{(abc)^{2/3}},~V=\frac{4{\rm\pi}}{3}abc,~P=\frac{2}{3}\frac{E_{th}}{V}.\end{eqnarray}$$
Using (8.6) in (8.4), we obtain the equation for the entropy function
${\it\varepsilon}$
in presence of different losses in the form
$$\begin{eqnarray}\frac{\text{d}{\it\varepsilon}}{\text{d}t}=(abc)^{2/3}U_{kin}\left[\left(\frac{2}{{\it\tau}_{rel}}-\frac{1}{{\it\tau}_{el}}-\frac{2}{{\it\tau}_{ml}}\right)-\frac{U_{rot}}{U_{g}}\left(\frac{2}{{\it\tau}_{Ml}}-\frac{1}{{\it\tau}_{ml}}\right)+\frac{U_{kin}}{U_{g}{\it\tau}_{ml}}\right].\end{eqnarray}$$
The equations, describing approximately the dynamics of the formation of a stationary dark matter object, include the equations of motion (3.10)–(3.12) after adding terms with relaxation and the energy equation (8.8) with account of (3.5), (3.6), (3.8) and (8.6), describing the losses of the energy, mass and angular momentum.
To obtain a numerical solution of the equations, we write them in non-dimensional variables. Let us introduce the variables
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{t}=\frac{t}{t_{0}},\quad \tilde{a}=\frac{a}{a_{0}},\quad \tilde{b}=\frac{b}{a_{0}},\quad \tilde{c}=\frac{c}{a_{0}},\quad \tilde{m}=\frac{m}{m_{0}},\quad \tilde{M}=\frac{M}{M_{0}},\quad \tilde{{\it\rho}}=\frac{{\it\rho}}{{\it\rho}_{0}},\\ \displaystyle \tilde{U} =\frac{U}{U_{0}},\quad \tilde{E}_{th}=\frac{E_{th}}{U_{0}},\quad \tilde{{\it\varepsilon}}=\frac{{\it\varepsilon}}{{\it\varepsilon}_{0}},\quad \tilde{{\it\tau}_{i}}=\frac{{\it\tau}_{i}}{t_{0}}.\end{array}\!\right\}\end{eqnarray}$$
The scaling parameters
$t_{0},a_{0},m_{0},M_{0},{\it\rho}_{0},U_{0},{\it\Omega}_{0},{\it\varepsilon}_{0}$
are connected by the following relations:
$$\begin{eqnarray}t_{0}^{2}=\frac{a_{0}^{3}}{Gm_{0}},\quad M_{0}^{2}=Ga_{0}m_{0}^{3},\quad U_{0}=\frac{Gm_{0}^{2}}{a_{0}},\quad {\it\rho}_{0}=\frac{m_{0}}{a_{0}^{3}},\quad {\it\Omega}_{0}=\frac{M_{0}}{m_{0}a_{0}^{2}},\quad {\it\varepsilon}_{0}=U_{0}a_{0}^{2}.\end{eqnarray}$$
Here
$U_{0}$
is used for scaling of all kinds of energies. Hereinafter the ‘tilde’ sign is omitted for simplicity. In non-dimensional variables we have
$m=(4{\rm\pi}/3)\,{\it\rho}\,abc$
,
$M=(m/5)\,{\it\Omega}(a^{2}+b^{2})$
,
${\it\tau}_{i}=2{\rm\pi}\,{\it\alpha}_{i}\sqrt{(abc)/3m}$
.
Taking into account violent relaxation, total energy, mass and angular momentum losses, the dynamics of the system is described by the following non-dimensional system of equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=-\frac{{\dot{a}}}{m}\frac{\text{d}m}{\text{d}t}-\frac{3m}{2}\,a\int _{0}^{\infty }\frac{\text{d}u}{(a^{2}+u){\it\Delta}}+\frac{10}{3m}\,\frac{1}{a}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}+\frac{25M^{2}}{m^{2}}\,\frac{a}{(a^{2}+b^{2})^{2}}-\frac{{\dot{a}}}{{\it\tau}_{rel}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{b}=-\frac{{\dot{b}}}{m}\frac{\text{d}m}{\text{d}t}-\frac{3m}{2}\,b\int _{0}^{\infty }\frac{\text{d}u}{(b^{2}+u){\it\Delta}}+\frac{10}{3m}\,\frac{1}{b}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}+\frac{25M^{2}}{m^{2}}\,\frac{b}{(a^{2}+b^{2})^{2}}-\frac{{\dot{b}}}{{\it\tau}_{rel}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{c}=-\frac{{\dot{c}}}{m}\frac{\text{d}m}{\text{d}t}-\frac{3m}{2}\,c\int _{0}^{\infty }\frac{\text{d}u}{(c^{2}+u){\it\Delta}}+\frac{10}{3m}\,\frac{1}{c}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}-\frac{{\dot{c}}}{{\it\tau}_{rel}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{{\it\varepsilon}}=(abc)^{2/3}\,U_{kin}\left[\left(\frac{2}{{\it\tau}_{rel}}-\frac{1}{{\it\tau}_{el}}-\frac{2}{{\it\tau}_{ml}}\right)-\frac{U_{rot}}{U_{g}}\left(\frac{2}{{\it\tau}_{Ml}}-\frac{1}{{\it\tau}_{ml}}\right)+\frac{U_{kin}}{U_{g}\,{\it\tau}_{ml}}\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\begin{array}{@{}c@{}}{\dot{m}}=\displaystyle -\frac{1}{3{\it\tau}_{ml}}\,({\dot{a}}^{2}+{\dot{b}}^{2}+{\dot{c}}^{2})\bigg/\displaystyle \left(\int _{0}^{\infty }\frac{\text{d}u}{{\it\Delta}}\right),\quad \displaystyle \frac{M}{M_{in}}=\displaystyle \left(\frac{m}{m_{in}}\right)^{{\it\tau}_{ml}/{\it\tau}_{Ml}},\\ \quad {\it\Delta}^{2}=(a^{2}+u)(b^{2}+u)(c^{2}+u).\end{array}\end{eqnarray}$$
This system is solved numerically for several initial parameters. The case of a spheroid (
$a=b\neq c$
) is considered in the paper of Bisnovatyi-Kogan (Reference Bisnovatyi-Kogan2004); the case of an ellipsoid in the paper of Bisnovatyi-Kogan and Tsupko (Reference Bisnovatyi-Kogan and Tsupko2005). For investigation of the dynamical behaviour, we start the simulation from a spherical body of unit mass (
$m_{in}=1$
) with zero or small entropy
${\it\varepsilon}\ll 1$
. We also specify the parameters characterizing different dissipations. In all calculations with relaxation we use the following values:
${\it\alpha}_{rel}=3,{\it\alpha}_{el}={\it\alpha}_{ml}={\it\alpha}_{Ml}=15$
. We set non-zero initial velocities. Because of the rotation around the
$c$
-axis, the initial sphere transforms to a spheroid during the collapse if the initial parameters for the axes
$a$
and
$b$
are exactly equal (
$a_{in}=b_{in},{\dot{a}}_{in}={\dot{b}}_{in}$
). Note that during the motion spheroids may be not only oblate but prolate too (Bisnovatyi-Kogan Reference Bisnovatyi-Kogan2004). To study an appearance of three-axis figures the initial parameters for the axes
$a$
and
$b$
had been slightly different in all variants of the calculations.
In the first variant (figure 8), there is a large initial angular momentum
$M_{in}=0.5$
. The field of velocities is slightly perturbed by increasing
${\dot{b}}$
in comparison with
${\dot{a}}$
. Initially we observe the collapse and the formation of the pancake. During the motion the difference (a–b) increases due to development of a secular instability, and we obtain the transformation of the Maclaurin spheroid into the Jacobi ellipsoid in the dynamics. Because of the relaxation, the oscillating motion damps, and the configuration reaches the equilibrium state of the three-axis ellipsoid. The corresponding behaviours of total energy
$U_{tot}$
, entropy function
${\it\varepsilon}$
and mass
$m$
are represented in figures 9–11. We see that the main changes of these parameters take place at the stages of the first few oscillations: the entropy function
${\it\varepsilon}$
increases, the mass
$m$
and the total energy
$U_{tot}$
decrease. Then these parameters reach the equilibrium values.
Development of an instability at large angular momentum, and formation of a stationary triaxial figure.
The total energy evolution for the case considered in figure 8.
The entropy function evolution for the case considered in figure 8.
The mass evolution for the case considered in figure 8.
In the second variant of the calculations (figure 12), we set a small initial angular momentum
$M_{in}=0.1$
. After the collapse, there are temporary oblate and prolate spheroid appearances during relaxation. There is no secular instability, and the system reaches the equilibrium oblate spheroid of Maclaurin. Here we have an additional instability: there is increasing of the difference (a–b) at early stages of motion. However, this difference remains too small for visualizing in the figure. This is the instability characteristic to the system with purely radial trajectories (Antonov Reference Antonov1973; Fridman and Polyachenko Reference Fridman and Polyachenko1985). This instability takes place in all variants of calculations, but in most cases it is very small and always disappears during the relaxation and formation of the equilibrium figure. In the cases of large angular momentum this instability is not revealed because the difference (a–b) quickly increases as a result of the secular instability.
The relaxation to the spheroid at small angular momentum.
For the illustration of the radial instability, we consider an initial sphere without rotation (
$M_{in}=0$
) and without any relaxation process. Only three equations (8.11)–(8.13) with
${\it\tau}_{rel}=\infty$
have been used, at constant
$m=1$
,
$M=0$
and
${\it\varepsilon}$
. We have set small initial differences between the axes
$a$
–
$c$
, and have obtained an unstable behaviour evidently visible in figure 13.
Radial instability development in the non-rotating body, without dissipation.
For numerical investigation of the secular instability, we consider the equations of motion without dissipation, and take the equilibrium spheroid with small perturbation as the initial configuration. It is easy to obtain equilibrium parameters of the spheroid from the equations of motion with zero accelerations. The initial equilibrium configuration of the Maclaurin spheroid was slightly perturbed by increasing
$b$
. At small angular momentum we obtain the conservation of the initial configuration. There are small oscillations around the equilibrium spheroid, during which the difference (a–b) changes sign periodically. At large angular momentum we obtain the development of the secular instability and, after several low-amplitude oscillations, the Maclaurin spheroid is transformed into the Jacobi ellipsoid (see figure 14, where the variant with large angular momentum is represented). Because of absence of relaxation, the system does not reach the equilibrium configuration and oscillations continue. If we include relaxation in this variant, we should observe the final equilibrium state of the ellipsoid. Our calculations allow us to find approximately numerically the point of the onset of the secular instability of the compressible Maclaurin spheroid.
Development of the secular instability in a rapidly rotating body, without dissipation.
9. Equilibrium configurations and stability
The equilibrium of the uniformly rotating figure (spheroid or ellipsoid) is found from (8.11)–(8.13) with zero time derivatives:
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{3m}{2}\,a\int _{0}^{\infty }\frac{\text{d}u}{(a^{2}+u){\it\Delta}}+\frac{10}{3m}\,\frac{1}{a}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}+\frac{25M^{2}}{m^{2}}\,\frac{a}{(a^{2}+b^{2})^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{3m}{2}\,b\int _{0}^{\infty }\frac{\text{d}u}{(b^{2}+u){\it\Delta}}+\frac{10}{3m}\,\frac{1}{b}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}+\frac{25M^{2}}{m^{2}}\,\frac{b}{(a^{2}+b^{2})^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{3m}{2}\,c\int _{0}^{\infty }\frac{\text{d}u}{(c^{2}+u){\it\Delta}}+\frac{10}{3m}\,\frac{1}{c}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}x=\frac{u}{a^{2}},\quad k=\frac{c}{a},\quad k_{1}=\frac{b}{a},\quad j=\frac{M}{m},\quad {\it\epsilon}=\frac{{\it\varepsilon}}{m}\end{eqnarray}$$
and integrals
$$\begin{eqnarray}\displaystyle I_{1}(k,k_{1}) & = & \displaystyle \int _{0}^{\infty }\frac{\text{d}x}{(x+1)^{3/2}\sqrt{(k_{1}^{2}+x)(k^{2}+x)}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle I_{2}(k,k_{1}) & = & \displaystyle \int _{0}^{\infty }\frac{\text{d}x}{(x+k_{1}^{2})^{3/2}\sqrt{(1+x)(k^{2}+x)}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle I_{3}(k,k_{1}) & = & \displaystyle \int _{0}^{\infty }\frac{\text{d}x}{(x+k^{2})^{3/2}\sqrt{(1+x)(k_{1}^{2}+x)}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{3m}{2}\,a\,I_{1}(k,k_{1})+\frac{10}{3}\,\frac{{\it\varepsilon}}{(kk_{1})^{2/3}}+\frac{25j^{2}}{(1+k_{1}^{2})^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{3m}{2}\,ak_{1}\,I_{2}(k,k_{1})+\frac{10}{3}\,\frac{{\it\varepsilon}}{k_{1}(kk_{1})^{2/3}}+\frac{25j^{2}\,k_{1}}{(1+k_{1}^{2})^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{3m}{2}\,ak\,I_{3}(k,k_{1})+\frac{10}{3}\,\frac{{\it\varepsilon}}{k(kk_{1})^{2/3}}. & \displaystyle\end{eqnarray}$$
$m$
, we obtain 
$$\begin{eqnarray}\displaystyle 0=-\frac{1}{k^{2}}\frac{I_{1}(k,k_{1})}{I_{3}(k,k_{1})}+1+\frac{15j^{2}}{2{\it\epsilon}}\frac{(kk_{1})^{2/3}}{(1+k_{1}^{2})^{2}},\quad 0=-\frac{k_{1}^{2}}{k^{2}}\frac{I_{2}(k,k_{1})}{I_{3}(k,k_{1})}+1+k_{1}^{2}\frac{15j^{2}}{2{\it\epsilon}}\frac{(kk_{1})^{2/3}}{(1+k_{1}^{2})^{2}}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
At
$a=b$
,
$k_{1}=1$
, (9.9) are identical, and determine the equilibrium of the Maclaurin spheroid; see Bisnovatyi-Kogan (Reference Bisnovatyi-Kogan2004). For Jacobi ellipsoids with
$a\neq b\neq c$
, we obtain the following relation between
$k$
and
$k_{1}$
:
$$\begin{eqnarray}F(k,k_{1})=1-\frac{1}{k_{1}^{2}}+\frac{1}{k^{2}}\frac{I_{2}(k,k_{1})}{I_{3}(k,k_{1})}-\frac{1}{k^{2}}\frac{I_{1}(k,k_{1})}{I_{3}(k,k_{1})}=0.\end{eqnarray}$$
This equation has a trivial solution
$k_{1}=1$
at all
$k$
, corresponding to the Maclaurin spheroid. Let us find the bifurcation point of (9.10), at which non-trivial solutions appear. While
$k_{1}=1$
is always a root of (9.10), we may write
$F(k,k_{1})=(k_{1}-1)f(k,k_{1})$
. An additional root of the equation
$F(k,k_{1})=0$
appears when the root of the equation
$f(k,k_{1})=0$
appears at
$k_{1}=1$
. The root of the zero derivative equation
$F_{k_{1}}^{\prime }(k,k_{1})=f(k,k_{1})+(k_{1}-1)f_{k_{1}}^{\prime }(k,k_{1})=0$
at
$k_{1}=1$
coincides with the root of the equation
$f(k,1)=0$
(We are grateful to A. I. Neishtadt for useful discussion of this point.); therefore, the value of
$k$
at the bifurcation point is determined by the equation
$$\begin{eqnarray}\left.\frac{\partial F(k,k_{1})}{\partial k_{1}}\right|_{k_{1}=1}=0.\end{eqnarray}$$
Using (9.5) and (9.10), this equation is written in the form
$$\begin{eqnarray}2k^{2}\,I_{3}(k,k_{1})|_{k_{1}=1}+\left.\frac{\partial I_{2}(k,k_{1})}{\partial k_{1}}\right|_{k_{1}=1}-\left.\frac{\partial I_{1}(k,k_{1})}{\partial k_{1}}\right|_{k_{1}=1}=0.\end{eqnarray}$$
At
$k_{1}=1$
, we have the analytic expressions
$$\begin{eqnarray}\left.\frac{\partial I_{1}(k,k_{1})}{\partial k_{1}}\right|_{k_{1}=1}=-I_{0},\quad \left.\frac{\partial I_{2}(k,k_{1})}{\partial k_{1}}\right|_{k_{1}=1}=-3I_{0},\quad I_{0}=\int _{0}^{\infty }\frac{\text{d}x}{(1+x)^{3}\sqrt{k^{2}+x}},\end{eqnarray}$$
so that (9.12) is reduced to
$$\begin{eqnarray}k^{2}\,I_{3}=I_{0},\quad \text{where}~I_{3}=I_{3}(k,1)\quad \text{and}\quad I_{2}=I_{2}(k,1)=I_{1}=I_{1}(k,1).\end{eqnarray}$$
Taking analytically the integrals
$I_{0},I_{1}=I_{2}$
and
$I_{3}$
, we have
$$\begin{eqnarray}(1-k^{2})I_{0}=\frac{3}{4}I_{2}-\frac{k}{2},\quad I_{3}=\frac{2}{k(1-k^{2})}\left(1-\frac{k\arccos k}{\sqrt{1-k^{2}}}\right),\end{eqnarray}$$
$$\begin{eqnarray}I_{1}=I_{2}=-\frac{1}{k(1-k^{2})}\left(k^{2}-\frac{k\arccos k}{\sqrt{1-k^{2}}}\right).\end{eqnarray}$$
Using (9.15) in (9.14), we obtain the equation
$$\begin{eqnarray}\frac{\arccos k}{\sqrt{1-k^{2}}}=\frac{k(13-10k^{2})}{3+8k^{2}-8k^{4}},\end{eqnarray}$$
whose solution
$k=0.582724$
determines the bifurcation point at the sequence of the Maclaurin spheroids. For the uniform spheroid, the position of this point does not depend on the adiabatic index of the matter.
Above we have obtained the bifurcation point on the equilibrium curve of the Maclaurin spheroids using only the equilibrium relations for the Jacobi ellipsoids. The usual way for investigation of stability is to solve linearized equations of motion and to find the eigenfrequencies.
The non-dimensional equations of motion without relaxation are
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{a}=-\frac{3m}{2}\,a\int _{0}^{\infty }\frac{\text{d}u}{(a^{2}+u){\it\Delta}}+\,\frac{10}{3m}\,\frac{1}{a}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}+\frac{25M^{2}}{m^{2}}\,\frac{a}{(a^{2}+b^{2})^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{b}=-\frac{3m}{2}\,b\int _{0}^{\infty }\frac{\text{d}u}{(b^{2}+u){\it\Delta}}+\,\frac{10}{3m}\,\frac{1}{b}\,\frac{{\it\varepsilon}}{(abc)^{2/3}}+\frac{25M^{2}}{m^{2}}\,\frac{b}{(a^{2}+b^{2})^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\ddot{c}=-\frac{3m}{2}\,c\int _{0}^{\infty }\frac{\text{d}u}{(c^{2}+u){\it\Delta}}+\,\frac{10}{3m}\,\frac{1}{c}\,\frac{{\it\varepsilon}}{(abc)^{2/3}},\quad {\it\Delta}^{2}=(a^{2}+u)(b^{2}+u)(c^{2}+u).\end{eqnarray}$$
At given
$a,b,c$
, the equilibrium values of
${\it\varepsilon}$
and
$M$
(the equilibrium configuration of the ellipsoid) are obtained using the formulae
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\varepsilon}=\frac{9}{20}\,m^{2}\,c^{2}\,(abc)^{2/3}\int _{0}^{\infty }\frac{\text{d}u}{(c^{2}+u){\it\Delta}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{25M^{2}}{m^{2}}\,\frac{a^{2}}{(a^{2}+b^{2})^{2}}=\frac{3m}{2}\,a^{2}\int _{0}^{\infty }\frac{\text{d}u}{(a^{2}+u){\it\Delta}}-\frac{3m}{2}\,c^{2}\int _{0}^{\infty }\frac{\text{d}u}{(c^{2}+u){\it\Delta}}. & \displaystyle\end{eqnarray}$$
We linearize (9.18) and (9.19) around the equilibrium configuration of the Maclaurin spheroid (
$a=b,\;k=c/a$
). Subtracting (9.18) and (9.19) and using the equilibrium values (9.21) and (9.22), we obtain
$$\begin{eqnarray}{\it\delta}(\ddot{a}-\ddot{b})=3\frac{m}{a^{3}}\,\left[-k^{2}\int _{0}^{\infty }\frac{\text{d}u}{(1+u)(k^{2}+u)^{3/2}}+\int _{0}^{\infty }\frac{\text{d}u}{(1+u)^{3}(k^{2}+u)^{1/2}}\right]{\it\delta}(a-b).\end{eqnarray}$$
The losses for
${\it\varepsilon}$
,
$U_{tot}$
,
$m$
and
$M$
are quadratic to perturbations, so these values remain constant in linear approximation.
Taking
${\it\delta}(a-b)\sim \exp (-\text{i}{\it\omega}t)$
, we come to the characteristic equation
$$\begin{eqnarray}{\it\omega}^{2}=3\frac{m}{a^{3}}\,\left[k^{2}\int _{0}^{\infty }\frac{\text{d}u}{(1+u)(k^{2}+u)^{3/2}}-\int _{0}^{\infty }\frac{\text{d}u}{(1+u)^{3}(k^{2}+u)^{1/2}}\,\right].\end{eqnarray}$$
These integrals are expressed in analytical functions:
$$\begin{eqnarray}\displaystyle {\it\omega}^{2}=\frac{3m}{a^{3}}\,\left[k^{2}\left(\frac{2}{k(1-k^{2})}-\frac{2\arccos k}{(1-k^{2})^{3/2}}\!\right)-\left(\!\frac{3\arccos k}{4(1-k^{2})^{5/2}}-\frac{3k}{4(1-k^{2})^{2}}-\frac{k}{2(1\!-\!k^{2})}\!\right)\right]. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
For the border of stability
${\it\omega}^{2}=0$
, we again obtain (9.17).
The spheroid loses its stability for the transformation into a three-axis ellipsoid, at the bifurcation point
$k=0.58272$
,
$e=\sqrt{1-c^{2}/a^{2}}=0.81267$
. Our approximate equilibrium equations, even without dissipation, describe the uniformly rotating ellipsoids (spheroids), which are not connected by adiabatic relations, and contain a ‘hidden’ non-conservation of the local angular momentum, which preserves the uniformity. Therefore, in presence of this ‘hidden’ non-conservation, the loss of stability takes place exactly in the bifurcation point. The loss of stability is not connected with the relaxation and follows from our equations even in the case when all global integrals remain constant. In the exact approach the instability in this point happens only in the direct presence of dissipative terms (Chandrasekhar Reference Chandrasekhar1969). The pure adiabatic spheroidal system preserves stability until
$e=0.952887$
, where it becomes unstable via a vibrational mode.
According to the hypothesis of Ostriker and Peebles (Reference Ostriker and Peebles1973), the stability of an isolated axially symmetric system is determined by the ratio
$U_{rot}/|U_{g}|$
. They determined from numerical experiments the critical value for various configurations as
$0.14\pm 0.03$
. It was found in Lai et al. (Reference Lai, Rasio and Shapiro1993) and Shapiro (Reference Shapiro2004) that compressible spheroids become secularly unstable to triaxial deformations at the bifurcation point, where
$U_{rot}/|U_{g}|=0.1375$
, independent of the adiabatic index
$n$
,
${\it\gamma}=1+1/n$
. Our formula gives exactly the same result, which is also confirmed by our numerical simulations.
10. Emission of very long gravitational waves and gravitational lensing by gravitational waves
Gravitational radiation is produced during the collapse of a non-spherical body. Gravitational radiation during the formation of a pancake was estimated by Thuan and Ostriker (Reference Thuan and Ostriker1974), and was improved by Novikov (Reference Novikov1975), who took into account the most important stage of the radiation during a bounce. The formula for the estimation of the total energy emitted during the collapse and bounce is
$$\begin{eqnarray}U_{GW}\approx 0.01\left(\frac{r_{g}}{a_{eq}}\right)^{7/2}\left(\frac{a_{eq}}{c_{min}}\right)\,Mc^{2}.\end{eqnarray}$$
Here
$a_{eq}$
is the equilibrium equatorial radius of the pancake,
$c_{min}$
is its minimal half-thickness during the bounce and
$r_{g}=2GM/c^{2}$
is the Schwarzschild gravitational radius of the body. From our calculations we have
$a_{eq}/c_{min}\leqslant 100$
. The value of
$r_{g}/a_{eq}$
we estimate using the observed average velocity of a galaxy in the rich cluster
$v_{g}\sim 3000~\text{km}~\text{s}^{-1}$
, and taking
$$\begin{eqnarray}\frac{r_{g}}{a_{eq}}\sim \left(\frac{v_{g}}{c}\right)^{2}\approx 10^{-4}.\end{eqnarray}$$
Then the energy carried away by the GW may be estimated as
$U_{GW}\approx 10^{-14}Mc^{2}$
. If all dark matter had passed through the stage of a pancake formation, then very long GWs with a wavelength of the order of the size of the galactic cluster have an average energy density in the universe of
$$\begin{eqnarray}{\it\varepsilon}_{GW}\approx 10^{-14}{\it\rho}_{dm}c^{2}\approx 3\times 10^{-23}~\text{erg}~\text{cm}^{-3}.\end{eqnarray}$$
Here we have used for estimation the average dark matter density
${\it\rho}_{dm}=3\times 10^{-30}~\text{g}~\text{cm}^{-3}$
. The average strength
$E_{GW}$
of the very long GW may be estimated, taking the relation (Landau and Lifshitz Reference Landau and Lifshitz1993)
$$\begin{eqnarray}{\it\varepsilon}_{GW}=\frac{c^{2}}{16{\rm\pi}G}{\dot{h}}^{2},\end{eqnarray}$$
where
$h$
is the metric tensor perturbation (non-dimensional) connected with the GW; we consider only the scalar, having in mind the averaged value of this perturbation. Taking into account
${\dot{h}}={\it\omega}h=2{\rm\pi}ch/{\it\lambda}$
,
${\it\lambda}\sim 10~\text{Mpc}$
is the wavelength of the GW of the order of the size of the cluster. From the comparison of (10.3) and (10.4), we obtain the averaged metric perturbation due to a very long GW in the form (Bisnovatyi-Kogan Reference Bisnovatyi-Kogan2004)
$$\begin{eqnarray}\bar{h}=\frac{2{\it\lambda}}{c^{2}}\left(\frac{G{\it\varepsilon}_{GW}}{{\rm\pi}}\right)^{1/2}\approx 6\times 10^{-11}\end{eqnarray}$$
for the values of
${\it\lambda}$
and
${\it\varepsilon}_{GW}$
mentioned above. Incidentally, we may expect a 10 times larger amplitude of the GW than the averaged over the volume value.
According to general relativity, any gravitational field can change the trajectory of photons or, in other words, deflect light rays. Hence, the gravitational field may act as a gravitational lens.
Gravitational lensing by gravitational waves in different cases was considered by many authors (see Braginsky et al. Reference Braginsky, Kardashev, Polnarev and Novikov1990; Faraoni Reference Faraoni1992; Damour and Esposito-Farèse Reference Damour and Esposito-Farèse1998 and references therein). It was found that the deflection angle vanishes for any localized GW packet because of transversality of GWs (Damour and Esposito-Farèse Reference Damour and Esposito-Farèse1998). Thus, if the photon passes through the finite GW pulse its deflection due to this wave is equal to zero. In Bisnovatyi-Kogan and Tsupko (Reference Bisnovatyi-Kogan and Tsupko2008b ), we notice that displacement between trajectories of the photon before and after passing the wave may occur. On the basis of this result, we obtain an approximate formula for the estimation of observational effects.
We confirm analytically the vanishing of the deflection angle for plane-wave pulses. However, we have found that the GW changes the photon propagation in another way, simply shifting its whole trajectory after passing through the GW (see figure 15). This displacement is found analytically for the photon passing through the plane GW. Displacement takes place mainly in the case of isolated wave pulses, which have a form similar to the top part of the sinusoid or when it has non-symmetrical top and bottom parts of the wave profile, and may, in principle, vanish for a periodic wave of a long duration.
The observational effect of the displacement in the trajectory of a photon.
The directions of photons passing through the GW packet do not change; therefore, any focusing of rays does not occur in this case. Thus, the displacement in trajectories does not lead to any magnification effect. But the displacement leads to a change of the angular position of an object for a distant observer. The change of the angular position due to passing of the light ray through the GW pulse
${\rm\Delta}{\it\alpha}_{d}$
is (see figure 15)
$$\begin{eqnarray}{\rm\Delta}{\it\alpha}_{d}=\frac{{\rm\Delta}y}{D_{s}}\simeq \frac{h{\it\delta}}{D_{s}},\end{eqnarray}$$
where
$h$
is the amplitude of the GW pulse,
${\it\delta}$
is its thickness and
$D_{s}$
is the distance between the source and the observer.
Let us estimate the change of the angular position for the GW pulses produced during formation of large-scale structure of the universe in dark matter using the following parameters:
$h=10^{-11}$
,
${\it\delta}=\text{Mpc}$
and
$D_{s}=100~\text{Mpc}$
. Then we obtain
$$\begin{eqnarray}{\rm\Delta}{\it\alpha}_{d}\simeq 2\times 10^{-8}~\text{arcsec}.\end{eqnarray}$$
Another possible manifestation of GWs is the influence of them on the polarization of the cosmic microwave background (Polnarev Reference Polnarev1985). It is known that primordial GWs lead to formation of B-mode polarization (Seljak and Zaldarriaga Reference Seljak and Zaldarriaga1997); see also the textbook of Gorbunov and Rubakov (Reference Gorbunov and Rubakov2011). This effect is now being hunted for by different observational projects; see e.g. BICEP2 (http://bicepkeck.org/). We suppose that B-mode polarization formation can be explained not only by primordial GWs but also by GWs appearing in the epoch of the formation of the large-scale structure of the universe. The main distinguishing feature of GWs formed during large scale structure formation is anisotropy in the direction to rich clusters of galaxies.
11. Conclusions
We present here an approximate way of investigation of collapse and dynamic stability based on a model of compressible homogeneous three-axis ellipsoids. The investigation of the dynamics is reduced to numerical solution of systems of ordinary differential equations of semi-axis time evolution. In the presence of dissipation processes there are also equations for mass, angular momentum, energy and entropy evolution. Derivation of the equations of motion has been performed by variation of the Lagrange function of the ellipsoid. Different equations of state can be easily considered in this model. The right-hand sides of the equations have an analytical form for the spheroid case and can be written in terms of elliptical integrals for the three-axis ellipsoid case. Equilibrium configurations can also be investigated by taking zero accelerations in all equations. We apply this method for investigation of dynamics of non-spherical stars and dark matter large scale structure formation.
Dynamic stability of non-spherical stars was investigated. We have solved the equations that describe, in a simplified way, the Newtonian dynamics of a self-gravitating non-rotating spheroidal body after loss of stability. We have found that contraction to a singularity occurs only in a pure spherical collapse, and deviations from spherical symmetry stop the contraction through the stabilizing action of nonlinear non-spherical oscillations. A real collapse occurs after damping of the oscillations because of energy losses, shock wave formation or viscosity.
We have investigated the dynamics of a three-axis dark matter ellipsoid. Equations of motion for the axes of a uniform compressible ellipsoid have been obtained by variation of the Lagrange function, in which the violent relaxation and losses of matter, energy and angular momentum have been included phenomenologically.
The system was solved numerically, until the formation of stationary rotating figures in presence of the relaxation. For lower angular momentum
$M$
we have the formation of a compressed spheroid, while at larger
$M$
we follow the development of three-axial instability and formation of a three-axial ellipsoid. The instability in this approximation happens at the bifurcation point of the sequence of Maclaurin spheroids, where the Jacobi ellipsoidal system starts.
The bifurcation point coinciding with the point of loss of stability is found analytically in the form of a simple formula, by static and dynamic approaches. Numerical and analytical considerations give identical results. Development of instability, connected with radial orbits, is obtained for slowly rotating collapsing bodies.
The weak but very long GWs emitted mainly in the stages of the collapse and pancake formation form a long-wave GW background. The existence of such a long GW in the space between the source and the observer may be registered as an action of the gravitational lens. The effect of lensing can be due to displacement of the whole trajectory if a light ray is passing through the GW pulse. The displacement leads to a change of the angular position of an object for a distant observer.
The challenging idea is explanation of B-mode polarization formation not only by primordial GWs but also by gravitational waves appearing in the epoch of the formation of the large-scale structure of the universe.
Acknowledgements
The work of GSBK and OYuT was partly supported by the Russian Science Support Foundation, Russian Basic Research Foundation Grant No. 14-02-00728, and the Russian Federation President Grant for Support of Leading Scientific Schools, Grant No. NSh-261.2014.2.
The work of GSBK was also partially supported by the Russian Foundation for Basic Research Grant No. OFI-M 14-29-06045.
