{'fraction_non_alphanumeric': 0.03509556275846652, 'fraction_numerical': 0.018143884356022502, 'mean_word_length': 4.533292104720676, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Application of the Galerkin method to the problem of stellar stability, gravitational collapse and black hole formation

Approximate approach based on the Galerkin method is suggested for the investigation of equilibrium stellar models, a relativistic collapse problem and black hole formation. Some results of its simplified version - energetic method- are reviewed, and equations for general Gelerkin method are presented.


Introduction
Relativistic collapse leading to a black hole formation is a final stage of the evolution of the sufficiently massive stars with a mass M > (30 − 50)M ⊙ . To describe these collapse dynamical equations in general relativity (GR) are solved numerically (Baumgarte et al., 1995(Baumgarte et al., , 1996. In order to simplify the problem for obtaining a rapid approximate answer for such a questions, like collapse behaviour of different stellar masses, comparative neutrino light curves and their spectra, a simplified method is suggested. It is based on the classical Galerkin method, widely used in different mechanical and physical problems, connected with numerical computations (Fletcher, 1984). For the problem of stellar collapse this method may be considered as a generalization of the well known energetic method, which was applied for calculations of stellar stability (Zeldovich andNovikov, 1965, 1967) and initial stages of stellar collapse (Bisnovatyi-Kogan, 1968) in post-newtonian approximation (PN). Here we formulate the Galerkin method for the problems of relativistic collapse in GR.
2 Energetic method in post-newtonian approximation (PN) Consider first the formulation of the energetic method in PN. The main idea of this method is an averaging of the static or dynamic equations over the volume, using a given prescribed density distribution inside the star. We consider here only spherically symmetric stars, where the density is prescribed by polytropic distribution with a central density, obtained from algebraic (static) or ordinary differential (dynamics) equation. Consider polytropic equation of state P = Kρ γ . Combining the equilibrium and continuity equations, we get the equation for ρ(r) d dr Where r is the Euler radial coordinate, ρ is the density. Transforming to dimensionless quantities ( Chandrasekhar, 1939 ) θ and ξ, defined as ρ = ρ c θ n , r = αξ, α = (n + 1)K 4πG ρ 1 n −1 c where n is an polytrope index, ρ c is the central density, gives the equilibrium Lane-Emden equation with the boundary conditions θ = 1, dθ dξ = 0 at ξ = 0. The stellar boundary corresponds to ξ = ξ 1 so that θ(ξ 1 ) = 0. The stellar mass M expressed in terms of variables (2) becomes M = 4π and hence from (3), Obviously, at n = 3, γ = 4/3 the stellar mass is independent of ρ c and is exactly determined by the constant K in the equation of state. The stellar mass increases at γ > 4/3 with increasing ρ c , and at γ < 4/3 decreases. If the polytropic power coincides with the adiabatic power, γ = γ ad , then the star is stable at γ > 4/3 and unstable at γ < 4/3, and γ = 4/3 corresponds to the boundary case and represents the indifferent equilibrium. Real stars are not polytropes, but the condition γ = 4/3 is approximately valid at the boundary of stability if γ is treated as an adiabatic power properly averaged over the star. A strict derivation of stability conditions, including GR, is made by use of the variational method (Harrison et al., 1965). Equating the first variation to zero yields the equilibrium equation, while the stability condition requires the second variation to be positive. For an isentropic polytrope with γ = 4/3, ρ c is arbitrary, whereas the density distribution θ(ξ) is invariant against homologous contraction or expansion. Let us treat these properties as valid also for the case where is γ = 4/3 only in average. We then derive the equilibrium and stability conditions, using the simplified variational method based on the assumption of homology and conservation of stellar structure at density variations (Zeldovich and Novikov, 1967), usually called energetic method. Write down the total energy of an instantaneously static star analogous to the potential energy of consevative mechanical system: The first term here represents the internal energy ε i , the second the newtonian gravitational energy ε G , and the third, ε GR , a small GR correction rg r = 2Gm c 2 r ≪ 1 is the small parameter evaluated by Zeldovich and Novikov (1967) for the matter distribution over a n = 3 polytrope. The term containing the Newtonian gravitational energy of equilibrium star may be explicitely evaluated for arbitrary polytropic equation of state (Landau and Lifshits, 1976).
For an adiabat with γ = γ ad , we have E = n P ρ and for a star in equilibrium takes place Where ε N is the total energy of a newtonian star. The total energy of a stable star is negative, therfore the stability requires that n < 3, γ > 4/3. The radius of a polytrope is, using (2), (4) and (5) Here are used the values for a polytrope of n = 3: ξ 1 = 6.89685, M 3 = 2.01824. The ratio of ρ c to the average density ρ M = 4π 3 ρR 3 is ρc ρ = 4π 54.18. From (9), (7) and (6) we have Only one parameter, ρ 1/3 c or R, varies with homologous variations: is an invariant function (11) and hence, the energy variations reduce to ordinary derivatives. Using (10)- (11) and taking the entropy to be constant at variations, we get from (6) the equilibrium condition The second derivative of the energy at S = const turns into zero on the boundary of stability: We use here the thermodynamic relations Equations (12) and (13) have been obtained by Bisnovatyi-Kogan (1966) and used for determining the boundary of stability for various stellar models.

Supermassive stars with a hot dark matter
For M > 10 4 M ⊙ the main reason of instability are GR effects. The entropy of such supermassive stars in critical state is so large that the pressure is determined mainly by the radiation with a small admixture of plasma, important for stability, but giving a very short time until the onset of instability. A common way to overcome this instability is to consider rotating superstars, what may postpone the moment of collapse to 3 × 10 4 years for solid body rotation with angular momentum and mass losses (Bisnovatyi-Kogan, Zeldovich and Novikov, 1967), and much longer for a differentially rotating star evolving with almost constant angular momentum (Fowler, 1966;Bisnovatyi-Kogan and Ruzmaikin, 1973). Formation of supermassive stars on early stages of the Universe expansion, their loss of stability with subsequent collapse or explosion (Bisnovatyi-Kogan, 1968;Fricke, 1973;Fuller et al, 1986) could be important for early formation of heavy elements, observed in spectra of the most distant objects with red shift ∼ 5, creation of perturbations for large scale structure formation, influence on small scale fluctuations of microwave background radiation (Peebles, 1987;Cen et al, 1993). A necessity of a presence of a dark matter in modern cosmological models makes it important to include it into stability analysis of supermassive stars. This was done by McLaughlin and Fuller (1996), who dealed with nonrotating superstars. The same problem for rotating superstars, using energetic method, was solved by Bisnovatyi-Kogan (1998). The rotational effects occure to be more important for realistic choice of parameters.

Stability analysis
In supermassive stars with equation of state P = P r + P g = aT 4 3 + ρRT , ( where R is a gas constant, a is a constant of a radiative energy density ) there is P r ≫ P g due to high entropy of such stars. Besides, such stars are fully convective and entropy is uniform over them, so the spatial structure is well described by a polytropic distribution, corresponding to γ = 4/3. The influence of a hot dark matter, which density does not change during perturbations, should be taken by account of a newtonian gravitational energy of the star in the dark matter potential, because GR effects of a dark matter are of a higher order of magnitude (McLaughlin, Fuller, 1996). For radiation dominated plasma there is a following expression for the adiabatic index, determining the stability to a collapse where β = Pg P = 4R S . In the radiationaly dominated supermassive star there is a unique connection between its mass M and entropy per unit mass S (Zeldovich and Novikov, 1967) where a is a constant of the radiation energy density, and numerical coefficient is related to the polytropic density distribution with γ = 4/3. At the point of a loss of stability the critical value of an average adiabatic index < γ > in selfgravitating nonrotating star with account of post-newtonian corrections is determined by a relation (Zeldovich and Novikov, 1967) Here averaging is done according to (13). From comparison between (15) and (17) we get a well known relation for a critical central density of a supermassive star stabilized by plasma Here and below we consider for simplicity a pure hydrogen plasma. Newtonian energy of a superstar ε nd in the gravitational field of uniformly distributed dark matter with a density ρ d is written as The gravitational potential of a uniform dark matter ϕ d is written as where R d is much larger then stellar radius R, and M d is a total mass of the dark matter halo. Stability does not depend on normalization of the gravitational potential so we shall omit the constant value in (19). It follows from (19) and . Critical value of the adiabatic index is obtained from (13), and for nonrotating superstar in presence of a dark matter the critical value of an average adiabatic index < γ > crnrot is determined by The relation for a critical density in presence of a dark matter is obtained by comparison of (15) and (21). We get Here relations (9) is used for R, and R 0 ρr 4 dr = 6.95 × 10 −4 ρ c R 5 , based on Emden polytropic distribution with γ = 4/3, was used (see Bisnovatyi-Kogan, 1989). Relation (22) is reduced to Solution of (23) is presented in Fig.1.

Stability of rotating stars
Consider a rigid rotation, when its energy is a small correction to the energy of radiation and the energetic method is a good approach. When losses of an angular momentum during evolution are negligible we distinguish between rapidly rotating (RR) and slowly rotating (SR) superstars. In RR case a superstar reaches the state of rotational equatorial breaking before loosing its dynamical instability, and in SR case instability comes first. If a superstar has an angular momentum J, then its rotational energy ε rot ≈ 1.25J 2 ρ 2/3 c M −5/3 , and a ratio ε rot /ε GR remains constant during evolution (Bisnovatyi-Kogan, Zeldovich, Novikov, 1967). In presence of rotation and dark matter the critical value of the adiabatic index < γ > crrot is written as and the relations for determination of a critical central density, instead of (23), is written as As follows from (25), a superstar does dot loose its stability when ε rot > |ε GR |. This qualitative result, obtained in the post-newtonian approximation, remains to be valid in a strong gravitational field and reflects a presence of a limiting specific angular momentum a lim = GM/c, so that a black holes with a Kerr metric may exist only at a < a lim (Misner, Thorne, Wheeler, 1973). A RR superstar in a course of the evolution reaches instead a limit of a rotational instability, and equatorial mass shedding begins, leading to a loss of an angular momentum. Such star will loose the stability when the anuglar momentum will become less then the limiting value. The stage of a mass loss was examined by (Bisnovatyi-Kogan, Zeldovich, Novikov, 1967), where it was shown that this stage may last about 10 times longer, then a maximum evolution time to approach the rotational instability point. RR star reaches the stage of a rotational instability at different central densities, depending on J, but the ratio of rotatonal and Newtonan gravitational energy on the mass-shedding curve remains constant (Bisnovatyi-Kogan, Zeldovich, Novikov, 1967) 1 The energy of a rotating supertar in equilibrium in presence of a hot dark matter may be written as In the main term for a superstar in equilibrium a relation is valid Taking into account (26), (28), we get an expression for an equilibrium energy along the mass-shedding curve (with variable J) The curve ε eq (ρ c ) has a minimum at the central density, determined by a relation The relation (26) is used for determination of an angular momentum of the superstar J = J 0 with ρ c from (30) and J = J 1 < J 0 with ρ c from (31). Solution of this equation is also given in fig.1. As may be seen from fig.1, the stabilizing effect of rotation on the mass-shedding curve at J = J 1 is more important, then stabilization by a hot dark matter, which influence decreases at increasing of a central density. With account of a longer state of the evolution with mass loss until reaching the point of the loss of stability at larger J, this conclusion becomes even stronger.

Collapse and explosions of supermassive stars
To study dynamical processes by the energetic method we use, instead of the energy variation, the energy conservation law in the form Figure 1: The correction terms δ GR and δ dm + δ GR , the quantities β/6 (line c), β/6+(ε rot /ε N ) sh /3 (line b), and 2β/6+ (ε rot /ε N ) sh /3 (line a), as functions of the central density of a supermassive star with M = 10 6 M ⊙ , and dark matter density of 10 −5 g/cm 3 . The instability points for nonrotating star occure at intersection of the correction term curves with the line c. Mass shedding in the stable star with angular momentum J 0 (see text) occures at intersection of correction term curves with the line b, and critical point on the mass-shedding curve is determined by a corresponding intersection with the line a.
Here ε k = 1 2 M 0 v 2 dm is the kinetic energy of the superstar. Using thermodynamic relation we get From the mass conservation law, in presence of homologycal motion with the fixed density distribution in space (11), we obtain a space velocity distribution. Considering r as a Lagrangian radius, corresponding to the mass m ≤ M , r(M ) = R, we may write Using definitions of the velocity in Lagrangian coordinates, we get its homolog- With account of the equality 1 0 (r 2 /R 2 )ϕ(ν)dν = 0.108, we get a dynamical equation in the form 0.597 On dynamical stages of the superstar evolution, after its loss of stability the radiative transfer losses (−Q r ) are less important, then neutrino losses (−Q ν ) and energy production in nuclear reactions Q n . The equation determining entropy changes is averaged over the star with a weight function (T 4 ∼ Eρ) for radiation dominated superstar, which entropy is taken homogenous Nuclear reactions of pp and CN O hydrogen burning and 3α helium burning have been included in the calculations of Bisnovatyi-Kogan (1968). The equaions (36),(37), together with equations for averaged composition of hydrogen X, helium Y , and equation of state with account of e + e − pairs creation and relativistic relations for electrons, had been solved numerically. Primodial chemical composition with only hydrogen and helium was taken as initial condition. In the process of contraction after a loss of stability, 3α reaction produces 12 C, which initiates a CN O hydrogen burning, pp reaction remaining always unimportant. It was obtained that in stars with M < 1.5 × 10 5 M ⊙ collapse is reversed, and they explode, enriching the intergalactic and interstellar gas with heavy elements. Such explosions could happen on stages, preceding the epoch of a galaxy formation. That is one of the way to explain high metallicity in the distant quasars and intergalactic gas in galaxy clusters. Similar calculations made for rotating superstars, and for normal (solar) composition shift the boundary between collapsing and exploding superstars to higher masses (Fricke,1973).

Formulation of Galerkin method
The Galerkin method allows to find an approximate solution of differential equations. In this method the solution of the partial or ordinary differential equation is represented in the form (Fletcher, 1984) where ϕ i are known analytical functions, and coefficients (or functions) α i are to be determined. Functions ϕ i are assumed to be not necessary orthogonal. When substituted to the original differential equation the coefficients α i will satisfy a set of algebraic or ordinary differential equations, obtained by minimizing the corresponding functional. Consider how the Galerkin method can be applied to the problem of the relativistic star collapse. The functional which should be minimized, is a total stellar energy. It contains a kinetic energy (32), which in PN is written explicitly. In GR the total energy cannot be splitted into different components, and the collapsing body should be described by means of a stress energy tensor T ij . Equating to zero the variation of the total energy gives dynamical equations of the collapse. In the process of contraction the collapse accelerates and finally approaches a free fall which is described by a Tolman solution (Landau and Lifshits, 1962). The approximate method should be formulated in such a way, that its solution approaches the Tolman one at final stages.
We may use the Galerkin method to approximate the system of equations in partial derivatives by the system of ordinary differential equations. The density, as a function of a lagrangian coordinate a and time t, ρ(a, t) is represented as a sum of α i (t)ϕ i (a). Since ϕ(a) is invariant under the lagrangian transformations, the partial time derivative of this product will transform to dαi dt ϕ i (a). Thus, the Galerkin method allows to separate variables t, and a. The same separation has to be done for a velocity function v(a, t).

Stability conditions in PN
In PN approximation the density in Galerkin method is written as a following sum As for a function ϕ 0 (a), it is convenient to take corresponding Emden profile for one of polytropic indices. For other functions we may chose ϕ k = cos 1+2K 2 πa, which have increasing number of nodes.
Then satisfaction of the boundary conditions: will be provided. The minimization of the energy functional (6) for finding an equilibrium model is reduced to zero partial derivatives leading in the static case of constant α i to a set of N algebraic equations for finding equilibrium α eq i . Stability of a model is found from an evaluation of the second variation δ 2 ε, which in the energetic method is reduced to the algebraic equation (13). In the Galerkin method with several scaling functions ϕ i (a), the second variation δ 2 ε is represented by a quadratic form In the energetic method the condition of stability reduces to one equation ∂ 2 ε/∂α 2 > 0, given in (13). In the Galerkin method the stability is related to positive definiteness of the quadratic form (42), what is provided (Smirnov, 1958) by the positiveness of the determinant and all its main minors. Remind, that main minors are determinants, obtaind from the main determinant (43) after eliminating lines and columns, intersecting on the main diagonal at i = k. For two functions in (39) the positiveness of the main determinant (43), and two partial derivatives ∂ 2 ε/∂α 2 1 > 0 and ∂ 2 ε/∂α 2 2 > 0 are enough for stellar stability. Loss of stability happens close before the point where the main determinant, or one of its main minors becomes zero.
In approximate presentation of the trial function in the Galerkin method, the minimal value of the second energy variation is larger, then its value for a real trial function. So zero values of the determinant (43), or one of its main minors, guarantees the onset of instability. Their positiveness is not an exact guarantee of the stability, but comparison of the energetic method with an exact stability analysis shows a good presicion of this approximate approach in most realistic cases (Bisnovatyi-Kogan, 1989). Energetic method corresponds to a homologeous trial function for displacement δr ∼ r. In the Galerkin method the trial function may be determined with a better precision. In fact, the coefficients δα i for the trial function of a density are determined as an eigenvector of a set of uniform linear equations The eigenvector δα e i is used for obtaining an approximate eigenfunction (44), and eigenvalues λ p are related to the square eigenfrequencies of the stellar model. The positive definiteness of the quadratic form (42) coincides with the positiveness of all eigenvalues λ p . So, Galerkin method permits to investigate a stability by finding approximate eigenvalues and eigenfunctions of a linear set of algebraic equations, instead of finding the same values from a second order differential equation exactly.

Relations for total energy and barion number in GR
In GR a barion number density n is used in (39) instead of the mass density ρ, which in this case is not presenting itself a conserved value. To obtain GR expression for the total energy of a spherically symmetric radially moving body take a stress-energy tensor for the perfect gas where ǫ is an internal energy per baryon, u i is a four velocity of the matter. Inside the spherically symmetric star the Schwarzschild type metric is used where dΩ 2 = dθ 2 + sin 2 θdϕ 2 , and λ is expressed as a functional of T 0 0 (Landau and Lifshits, 1962) To avoid confusion, we determine r and t to be independent Euler variables and v to be a velocity field of the medium. A number of baryons a inside a Lagrangian radius r(a) is considered as a Lagrangian independent coordinate Then, the world line of each shell is r(t, a). There exist a set of transformations of the Euler type metric (47) to its Lagrangian representation: where τ scales time of the co-moving observer. We determine a physical velocity of the shell as a displacement dl, measured by the Euler rest observer, with respect to Euler time, measured by the same observer: Where φ to be determined from (47). In a case of a test particle it gives a well known result for the free falling particle in a Schwarzschild metric Covariant components v α and v 2 are given by the following relations (Landau and Lifshits, 1962) In a case of a diagonal metric, the dimensionless 3 -metric tensor is connected with the 4 -metric tensor as γ αβ = −g αβ . The mass-energy functional for the spherically symmetrical relativistic star is written as in (48) Four velocity of the fluid can be expressed in terms of the tree velocity (Landau and Lifshits, 1962) Then, T 0 0 component of the stress-energy tensor (46) is written as From (48),(49) making use of (56) we get

Velocity distribution over the star in GR
Assuming the profile law for n, let us obtain the profile law for the velocity. To separate variables a and t we need a relation between the pertubations of the boundary (∂r/∂t) a=N , and velocity of the inner shell (∂r/∂t) a . Differentiating (58) gives what contains terms withė. After differentiation, equation (57) contains a second order derivativer. Going on this iterative process, we obtain an infinite system of differential equations: . . ., completed by a similar set of equations obtained from differentiation of (59). In PN the velocity distribution (35) is found from the continuity equation, and than is used in derivation of the equation of motion (36). In full GR description the infinite set of equations (60)-(62) should be used for obtaining the velocity distribution, instead of (35). Such complication is connected with an influence of the 'kinetic' energy on the geometry of the star interior. After substituting of expansion on n, similar to (39), we obtain ordinary differential equations instead of the partial ones. In practice the system (60) -(62) has to be cutted off on some step. Final stages of the gravitational collapse are close to a free fall of a dust spherical cloud. Thus usingė from the Tolman solution should be good for late stages of the collapse. For the zero-pressure dust collapse without thermal processes the energy inside a given lagrangian shell does not change, so in (59) the relationė(a) = (∂e/∂t) a = 0 could be approximately used, making unnesessary to solve the equations (60) -(62) for obtaining the profile v(a).

Equations of motion
The equation determining the change of the entropy of the star is written similar to (37). In that case the entropy distribution over the star should not change its form in the process of the collapse, or be uniform. The integration should be done over the lagrangian coordinate da = 4πnr 2 dr 1− 2Ge c 4 r . To derive equation of motion we should, like in (36), to take derivatives from all thermodynamic functions at constant entropy distribution, as follows from (32),(33). When only one term is is left in (39), as in the energetic method, the equation (60) is used also for derivation of the equation of motion. In the case of several terms in (39), the equation (60) should be written for each coefficient α k , taking constant other α i , i = k. This the number of equations will be equal to the number of unknown functions α i . This procedure is equivalent to a minimization of the energy functional over all coefficients α i (t), or over all values of α i for the equilibrium solution.
It is important to stress, that when using (60) -(62) for derivation of the velocity profile, all coefficients α i (t) must be differentiated together.

Conclusion
We suggest an approximate method for solving the problem of stellar stability, and dynamical stages of relativistic stellar collapse, using Galerkin method. It should be mentioned, that this paper, together with the general review of the contemporary state of art of the problem contains the description of the essentially new results. The stability analysis of the rotating stars with the presence of the surrounding hot dark matter background, as well as the application of the energetic method together with the Galerkin method to the problem of relativistic gravitational collapse are such new results.
In the case of stability investigation we get a set of algebraic equations and search of the eigenvalues of matrics, instead of solution of eigenvalue problem for differential equation of a second order. In the case of collapse set of ordinary differential equations should be solved, instead of partial ones. It is evident, that solution of ordinary equations is much simpler, what could permit to investigate a wide set of equation of states, stellar masses, and to get different neutrino light curves during a black hole formation.
Several examples of application of the simplest version of Galerkin method, a well known energetic method, are presented. A simple solution obtained by this method for complicated problems show explicitely the advantage of such approximate approach.