format: text

text:
# The Quark Gluon Pion Plasma

## Abstract

While it is commonly believed that there is a direct transition from the hadronic to a quark gluon phase at high temperature, it would be prejudicial to rule out a sequence of dynamically generated intermediate scales. Using as guide, an effective lagrangian with unconfined gluons and constituent quarks, interacting with a chiral multiplet, we examine a scenario in which the system undergoes first-order transitions at $`T_{comp}`$, the compositeness scale of the pions, at $`T_\chi `$, the scale for spontaneous chiral symmetry breaking, and at $`T_c`$, the confinement temperature.

We find that at current energies, it is likely that the formation temperature of the plasma, $`T_0<T_{comp}`$, and that this is therefore a quark gluon pion plasma (QGPP) rather than the usual quark gluon plasma (QGP). We propose some dilepton-related signatures of this scenario.

We know that quarks and gluons are confined as hadrons and chiral symmetry is spontaneously broken. We do not however know if these two phenomena set in at an identical temperature. It is unlikely that the chiral symmetry restoration energy/temperature scale is lower than that of confinement, because if it were, hadrons would show parity doubling below the confinement but above the chiral breaking scale. This is not seen either experimentally or in finite temperature lattice simulations. Determining whether it is higher even in the simplest case of QCD with a two flavour $`SU_2(L)\times SU_2(R)`$ chiral symmetry is not straight-forward. Indeed, while there is a bonafide order parameter for chiral symmetry breaking — namely, the mass of the constituent quark — the Wilson loop in the presence of dynamical quarks, is not a valid order parameter for confinement — and there are no other known candidates. However, by looking at energy density or specific heat, we can get a fair idea of the change in the number of operational degrees of freedom or particle modes with temperature. The relevant lattice calculations indicate that the drop from the large number of degrees of freedom in the QGP/QGPP phase to a few degrees of freedom in the hadronic one takes place in one broad step in temperature, suggesting that the two transitions are close but not necessarily identical, i.e. that $`T_\chi \stackrel{\stackrel{}{>}}{\stackrel{}{}}T_c`$.

In the same spirit, we could ask whether pions necessarily come apart the moment chiral symmetry is restored, or equivalently, whether the compositeness scale for the pion, $`T_{comp}`$, coincides with $`T_\chi `$. The following analogy with superconductivity (SC) leads us to conclude that it may not:

In the usual BCS theory, the formation of Cooper pairs at $`T_{crit}`$ coincides with the appearance of a SC order parameter. This corresponds to $`T_{comp}=T_{crit}`$, as there are no Cooper pairs at scales above $`T_{crit}`$. Here, $`T_{comp}`$ is the compositeness scale for Cooper pairs, which are normally described as undergoing momentum space pairing, the interelectron distance being much less than the size of a Cooper pair. This is indeed what happens for a weak pairing interaction. For a much stronger interaction, Cooper pairs may actually form by real space pairing. The Cooper pair size is then smaller than the interelectron distance and it is possible that, even when SC is lost, i.e. $`T>T_{crit}`$, pairs continue to exist simply as bound states.

Evidence for the related precursor phenomena for the chiral transition had been suggested long ago . Recently, there have been several calculations of the QCD counterpart of the non-condensed, paired, pseudogap phase of high $`T_c`$ superconductivity , which suggest that the strong pairing situation is likely for the strong interactions.

We shall accordingly take $`T_{comp}>T_\chi `$, and argue that heavy-ion collisions are the ideal platform for testing this assumption. In doing so, we shall use the following $`SU_2(L)\times SU_2(R)`$ effective lagrangian, in which the quarks are coupled to a chiral multiplet, $`[\pi ,\sigma ]`$ , to keep track of particle modes at various scales:

$`L={\displaystyle \frac{1}{4}}G_{\mu v}^aG_{\mu v}^a{\displaystyle \overline{\psi }\left(D+g_y(\sigma +i\gamma _5\stackrel{}{\tau }\stackrel{}{\pi })\right)\psi }`$
$`{\displaystyle \frac{1}{2}}(_\mu \sigma )^2{\displaystyle \frac{1}{2}}(_\mu \stackrel{}{\pi })^2{\displaystyle \frac{1}{2}}\mu ^2(\sigma ^2+\stackrel{}{\pi }^2){\displaystyle \frac{\lambda }{4}}(\sigma ^2+\stackrel{}{\pi }^2)^2+\text{const}`$ (1)
This has chiral symmetry-breaking but no confinement. The masses of the scalar/pseudoscalars and fermions, obtained by minimizing the potentials above, are $`\mu ^2=\lambda <\sigma >^2`$ and $`m_\sigma ^2=2\lambda <\sigma >^2`$ respectively. This theory reproduces several broad features of the strong interaction at the mean field level , indicating that it is, perhpas, a reasonable guide to the physics at intermediate scales. More specifically,

1. It provides a nucleon, which is realized as a soliton with quarks bound in a skyrmion configuration . Such a nucleon:

1. Gives a natural explanation for the ‘Proton spin puzzle’: Quarks in the background fields are in a spin, ‘isospin’ singlet state in which the quark spin operator averages to zero. If we collectively quantize the soliton to get states of good spin and isospin, the quark spin operator picks up a small, non-zero expectation value .
    2. Seems to naturally produce the Gottfried sum rule .
    3. Yields, from first principles (albeit with some drastic QCD evolution), a set of structure functions for the nucleon which are close to the experimental ones .
2. As a finite temperature field theory, it yields screening masses that match with those obtained from the lattice simulation of finite temperature QCD with dynamical quarks .
3. It gives a consistent equation of state for strongly interacting matter at all densities .

This lagrangian and the above discussion on temperature scales lead us to consider the following scenario:

1. For $`T>T_{comp}`$, we have gluons and massless current quarks, in a pure QGP phase.
2. For $`T_{comp}>T>T_\chi `$, we have a QGPP phase for which the particle mode count is (i) 16 modes from 8 massless gluons, (ii) 24 from 2 flavors of massless current quarks and anti-quarks, and (iii) 3 degenerate pions, which are present, in this temperature range, in the form of non-Goldstone-boson bound states, which we assume to be massless. (We neglect the $`\sigma `$, which would otherwise contribute one more mode).

The operational degrees of freedom in this phase are then 40 — three more than in the QGP phase. This increase may be visible as a slight bump on an ‘Energy Density’ vs ‘Temperature’ plot, in a high resolution lattice simulation.
3. For $`T_\chi >T>T_c`$, chiral symmetry is spontaneously broken and the degrees of freedom are (i) gluons, (ii) constituent quarks, which acquire a constituent mass of $`300400`$ MeV from the spontaneous breaking of $`\chi `$-symmetry and are, as a result, Boltzman suppressed given that $`T_\chi T_c`$ is expected  to be $`O(200MeV)`$, (iii) pions as Goldstone bosons, and (iv) the $`\sigma `$, which also acquires a mass and may be neglected.

The effective degrees of freedom are thus 19, at this stage.
4. For $`T_c>T`$, quarks are confined, and the only massless modes left are the 3 Goldstone bosons/pions.

We can get an estimate of $`T_{comp}`$ by examining the behaviour of the above Lagrangian as a function of energy. We find that at a certain scale, both the wave function renormalization constant, $`Z`$, for the $`\pi `$ and $`\sigma `$, and the quartic scalar interaction, vanish, leaving us with a Yukawa and a scalar mass term. The vanishing of $`Z`$, or equivalently, of the kinetic term for the scalar and pseudoscalars, means that these are no longer dynamical degrees of freedom, i.e. they have ceased to exist as composite entities. Eliminating them, via their field equation, leaves behind a four-fermi term which gets weaker with increasing energy. While these results, gleaned using perturbative RNG, cannot be used all the way out to where Z goes to zero (since the Yukawa coupling has a Landau singlarity there), they can perhaps give a reasonable estimate for the compositeness scale. This works out to be around $`700800`$ MeV , well above the temperatures attained in the present-day heavy ion colliders. This implies that the initial state produced by these is very probably a QGPP.

To investigate this possibility further, we adapt to our scenario, the considerations of , which gives a very transparent and simple treatment of the underlying physics. In particular, following , we treat each of our transitions as first order, with each producing a mixed phase in which the number of degrees of freedom changes continuously, as one phase gives way to the other.

A plasma which thermalizes at a temperature, $`T_0`$, where $`T_{comp}<T_0<T_\chi `$, a proper time, $`t_0`$, after nuclear impact, then evolves as follows:

1. It Bjorken expands and cools , reaching $`T_\chi `$ at time, $`t_{x1}=(T_0/T_\chi )^3t_0`$.
2. At $`T_\chi `$, it undergoes a first-order transition from the quark-gluon-pion to the chirally broken, gluon-pion phase, from $`t_{x1}`$ to $`r_1t_{x1}`$, where $`r_1=40/19`$ is the ratio of the degrees of freedom in the initial and final phases.

N.B. If we assume that chiral restoration and deconfinement occur simultaneously at $`T_\chi =T_c`$, the quark phase lasts in the mixture from time $`t_1=t_{x1}`$ to $`t_2=rt_1`$ (where $`r=37/3`$), which is substantially longer.
3. Having passed to the chirally-broken phase, the plasma undergoes a second Bjorken expansion and thereby cools from $`T_\chi `$ to $`T_c`$. Since $`T_\chi T_c`$, this part of the evolution may, however, be neglected.
4. At $`T_c`$, we get a mixture of the gluon-pion and the purely hadronic confined phases, which lasts from $`r_1t_{x1}`$ to $`t_{x2}=r_2r_1t_{x1}`$, where $`r_2`$ is $`19/3`$.
5. Finally the pion phase expands from $`T_c`$ to $`T_f`$, the freezeout temperature , at which pions loose thermal contact with one another. This phase begins at time $`t_{x2}`$ and ends at $`t_f=(T_c/T_f)^3t_{x2}`$.

The basic changes in this scenario when compared to that of  are that, (i) the pionic phase starts much earlier, at the initial temperature, $`T_0`$, and not after confinement at $`T_c`$ and (ii) the quark phase is suppressed above $`T_c`$ and is thus shortened.

To see how these changes are reflected in the number of dileptons emitted by the plasma, we note that these come from $`\pi ^+\pi ^{}`$ and $`q\overline{q}`$ annihilations. The cross-section for $`q\overline{q}l^+l^{}`$, summed over the spin, flavour and color of all quarks is

$`\sigma _q[M]=F_q{\displaystyle \frac{4\pi }{3}}({\displaystyle \frac{\alpha }{M}})^2\sqrt{14({\displaystyle \frac{m_l}{M}})^2}\left(1+2({\displaystyle \frac{m_l}{M}})^2\right)`$ (2)
where the numerical factor, $`F_q=20/3`$, if we include just the $`u`$ and $`d`$ quarks in the flavor sum.

The cross-section for $`\pi ^+\pi ^{}l^+l^{}`$ likewise is,

$`\sigma _\pi [M]=F_\pi [M](4\pi /3)(\alpha /M)^2\sqrt{14(m_l/M)^2}(1+2(m_l/M)^2)\sqrt{14(m_\pi /M)^2}`$ (3)
where the pion form factor is given by,

$`F_\pi [M]=m_\rho ^4/((m_\rho ^2M^2)^2+m_\rho ^2\mathrm{\Gamma }_\rho ^2)+(1/4)m_\rho ^{}^4/((m_\rho ^{}^2M^2)^2+m_\rho ^{}^2\mathrm{\Gamma }_\rho ^{}^2)`$ (4)
The masses and decay widths of $`\rho `$ and $`\rho ^{}`$ resonances are $`m_\rho =775`$ MeV, $`m_\rho ^{}=1.6`$ GeV, $`\mathrm{\Gamma }_\rho =155`$ MeV, $`\mathrm{\Gamma }_\rho ^{}=260`$ MeV respectively.

The rate for producing pairs, per unit space-time, with invariant mass squared, $`M^2`$, and transverse energy, $`E_T=\sqrt{p_T^2+M^2}`$, where $`p_T`$ is the momentum transverse to the beam axis, is (for $`a=q,\pi `$),

$$\frac{dN_a}{d^4xdM^2dE_T}=\frac{\sigma _a[M]M^2}{4(2\pi )^4}(14\frac{m_a^2}{M^2})E_TK_0(E_T/T)=A_a[M]E_TK_0(E_T/T)$$
(5)
While using Bjorken’s model, it is expedient to write $`d^4x=d^2x_Tdytdt`$ where $`t`$ is the proper time and $`y`$, the rapidity, of the fluid element. For central collisions of equal mass nucleii $`d^2x_T=\pi R_A^2`$, where $`R_A`$ is the nuclear radius. We can further perform the $`t`$-integral to get:

$`{\displaystyle \frac{dN_q^{MT}}{dydM^2dE_T}}=\pi R_A^2A_q(M)\{3E_T^5T_0^6t_0^2(G[E_T/T_0]G[E_T/T_c])`$
$`+E_TK_0[E_T/T_c](1/2)(r_11)t_{x1}^2\},`$ (6)
for the quark channel in the multiple transitions (MT) scenario. The first term in braces comes from integrating along the cooling curve $`t=(T_0/T)^3t_0`$ from $`T_0`$ to $`T_c`$ and the second, from the coexistence line at $`T_\chi `$, which we approximate as being equal to $`T_c`$. The shortening of the quark phase is reflected in the presence here of $`r_1=40/19`$, as opposed to $`r=37/3`$ for the single transition (ST) at $`T_\chi =T_c`$ case.

The corresponding rate for $`\pi ^+\pi ^{}`$-annihilations is likewise,

$`{\displaystyle \frac{dN_\pi ^{MT}}{dydM^2dE_T}}=\pi R_A^2A_\pi (M)\{3E_T^5T_0^6t_0^2(G[E_T/T_0]G[E_T/T_c])`$
$`+3E_T^5T_c^6t_{x2}^2(G[E_T/T_c]G[E_T/T_f])+E_TK_0[E_T/T_c](t_{x2}^2t_{x1}^2)/2\}`$ (7)
where,

$`G[x]`$ $`=`$ $`x^3(8+x^2)K_3[x]`$ (8)
and $`K_i[x]`$ are the modified Bessel functions.

In this case, the first term in braces (absent altogether for the ST case) results from the cooling of pions in the QGPP, from $`T_0`$ and $`T_c`$, and the second, from their cooling in the hadronic phase, from $`T_c`$ to $`T_f`$. The latter coincides with the ST result if $`t_{x2}t_2`$. The last term comes from integrating over time, the volume fraction, $`f(t)`$, of the pion phase, on the coexistence line $`T_\chi T_c`$. Since pions are present both above and below each of the coexistence lines at $`T_\chi `$ and $`T_c`$, $`f(t)=1`$ for both these transitions. In the approximation, $`T_\chi =T_c`$, the integral is then simply,

$$_{t_{x1}}^{t_{x2}}t𝑑t=(t_{x2}^2t_{x1}^2)/2.$$
(9)
The analogous integral in the ST case is, by contrast, $`(r1)rt_1^2/2`$.

We note that these differential rates further integrate over $`E_T`$ to the closed-form expressions:

$`{\displaystyle \frac{dN_q^{MT}}{dydM^2}}=\pi R_A^2A_q(M)T_0^2t_0^2\{6(T_0/M)^4(H[M/T_0]H[M/T_c])`$
$`+(M/T_c)(T_0/M)^4K_1[M/T_c](r_11)\}`$ (10)
where,

$$H[x]=x^2(8+x^2)K_0[x]+4x(4+x^2)K_1[x]$$
(11)
and we have translated ratios of times into ratios of temperatures. Similarly,

$`{\displaystyle \frac{dN_\pi ^{MT}}{dydM^2}}=\pi R_A^2A_\pi (M)T_0^2t_0^2\{6(T_0/M)^4(H[M/T_0]H[M/T_c])+6(r_1r_2)^2(T_0/M)^4`$
$`(H[M/T_c]H[M/T_f])+(M/T_c)(T_0/T_c)^4K_1[M/T_c][(r_1r_2)^21]\}`$ (12)
Finally, the differential rates, $`dN_{tot}^{MT}/dydM^2dE_T`$ and $`dN_{tot}^{MT}/dydM^2`$ for the total number of pairs, are obtained by summing the right hand sides of Eqs.(The Quark Gluon Pion Plasma) and (The Quark Gluon Pion Plasma), and Eqs.(The Quark Gluon Pion Plasma) and (The Quark Gluon Pion Plasma) respectively.

We can now compare our scenario, in quantitative detail, with that of Ref.. For ease of comparison, we use the same parameter-values: $`t_0=1fm/c`$, and $`\pi R_A^2=127fm^2`$.

* In Fig.1, we plot $`dN_{tot}^{MT}/dydM^2dE_T`$, $`dN_q^{MT}/dydM^2dE_T`$ and $`dN_\pi ^{MT}/dydM^2dE_T`$ vs $`M`$, for the combinations of $`T_0`$ and $`T_c`$ used in Ref. but at values of $`E_T`$, for which (in the first two cases) the ST pion peaks, coming from annihilation via the $`\rho `$ and $`\rho ^{}`$ resonances, are largely extinguished. We note that in the MT scenario, the rate displayed reduces in magnitude as $`E_T`$ increases, but the pion peaks survive and stay well above the $`q\overline{q}`$ contribution. This results from pions being created when the plasma forms, as opposed to when it hadronizes. It is consistent with the observation in  that peak extinction at large $`E_T`$ and the initial temperature of the pions are distinctly correlated: the smaller this temperature, the smaller the $`E_T`$ at which the peaks get extinguished.
* We note that an enhancement in in the production of large $`E_T`$ dileptons, owing to pions being available at higher temperatures, will inevitably be accompanied by an enhancement in the production of large $`p_T`$ pions, a feature highlighted by Schucraft  long ago. We can add that pions present at the boundary of the initial plasma, will, being colourless, have no problem escaping from the interaction region. The availability of pions at higher temperatures is, incidentally, supported by the observation that the HBT size of the initial pion source is smaller than if pions formed only after hadronization .
* We would intuitively expect the difference between the two scenarios to depend (among other factors) on the time spent by the system in the QGPP phase. This would in turn be determined by the initial temperature, $`T_0`$. In Fig.2, we accordingly plot the ratio of $`dN_{tot}/dydM^2`$ for the MT and ST scenarios vs $`T_0`$ over the range $`T_\chi T_0T_{comp}`$, for several values of $`M`$. We find that the MT rate is $`2030\%`$ higher, when $`M`$ is not too far from the $`\rho `$, $`\rho ^{}`$ resonances, but drops sharply as $`M`$ increases and the drooping pion form factor begins to cut down the pion contribution.
* We have noted that quarks acquiring constituent masses larger than $`T_\chi `$ and decoupling at the chiral phase transition, reduces the lifetime of the quark phase above $`T_c`$. This, in turn, brings down the $`q\overline{q}l^+l^{}`$ rate for small $`M`$ and $`E_T`$. In Fig.1, the difference is visible but only when $`T_c=240MeV`$ and $`T_0=250MeV`$. In the total rate, the effect is swamped by the pion contribution for those values of $`M`$ for which the pion form factor is appreciable. However, from Fig.3, we see that for slightly larger values of $`M`$, the pion contribution disappears and the total differential rate $`dN_{tot}/dydM^2`$ reduces to just the quark contribution, for both the ST and MT cases. Determining the rate in this $`M`$-window would then give an indication of whether a chiral transition has occurred close to but distinctly above $`T_c`$.
* The behaviour of these curves as functions of $`T_0`$ can be further understood as follows: In the ST scenario, a higher $`T_0`$ does not change the pion contribution to dilepton production. However, it clearly increases the time spent by the quarks in the cooling phase and thus enhances the quark contribution to dileptons, pushing the crossover to the quark dominated regime to smaller M. For the MT scenario, both the quark and pion contrbutions to dileptons are enhanced and so the crossover to the quark-dominated regime is not expected to change much (though it may still move slightly towards smaller M, due to steep fall in the pion form factor). For both scenarios, the cooling term in the quark contribution increases in magnitude with $`T_0`$ while the Maxwell co-existence term (which accounts for the difference between the two scenarios) remains unchanged. Thus as $`T_0`$ increases, the ST and MT quark contributions approach each other.
* Finally, in Fig.4, we plot the ratio of $`dN_{tot}/dydM^2`$ for the MT and the No Transition (NoT) scenarios vs $`T_0`$. The latter corresponds to there being no phase transition, i.e. to the fireball’s cooling in the pure hadronic phase. It is seen that this ratio is enormous.
* We conclude that if the dilepton rates increase dramatically, a phase transition has probably occurred. If the pion peaks in $`dN_{tot}/dydM^2dE_T`$ vs $`M`$ survive to high values of $`E_T`$, the plasma has a good chance of being a QGPP. If the peaks do not survive, but rather reduce to a flat line, we are quite likely seeing a quark phase. The actual value of this constant rate, for the window in $`M`$ discussed above, will then decide whether the transition to the confined hadronic phase occurred directly or through a chirally broken phase.

In summary, we have argued that there are reasons to believe that a QPG may form in heavy-ion collisions, as a result of several, as opposed to a single, phase transition. We have explored this possibility on the basis of some simple assumptions. We have found that this scenario pushes to higher temperatures the advent of the QGP, but opens up, in return, the exciting possibilty of a QGPP, containing pions as non-Goldstone-boson bound states, well above $`T_\chi `$. We have further investigated experimental signatures based on dilepton production, which would help distinguish the single- from the multiple-transition scenario. Making this distinction would shed light not merely on how the QGP forms but on a range of assumptions which lie at the very core of our understanding of strong interaction dynamics.

Acknowledgement

VS would like to thank many colleagues and collaborators, George Ripka, Manoj Banerjee, Bojan Golli, Mike Birse, W. Broniowoski, J.P. Blaizot, N. D. Haridass, G. Baskaran, M.Rho and many others. VS thanks R. Rajaraman for suggesting that an experimental signature of the model used here should be presented in the context of heavy ion physics. We apologize for missing many relevant references — the list is too extensive.

Figure Captions

Figure 1 Plots of $`dN_{tot}^{MT}/dydM^2dE_T`$ vs $`M`$ for three $`(T_0,T_c)`$-combinations. For the values of $`E_T`$ in Fig.1.a and Fig.1.b, the pion peaks in the single transition rate are essentially non-existent. Also plotted are the quark contributions for both the MT (thick line) and ST (thin line) scenarios. These can be distinguished only for $`T_0=250MeV`$ and $`T_c=240MeV`$. The MT quark contribution is lower because of the shortening of the quark phase.

Figure 2 Plots of the ratio of the $`dN_{tot}/dydM^2`$ rate vs $`T_0`$ for the MT and ST scenarios. For the lower $`M`$-values considered in these graphs, the pion form factor is appreciable. It seen that the MT rate is distinctly higher than the ST rate.

Figure 3 Plots of the ratio of the $`dN_{tot}/dydM^2`$ rates for the MT and ST scenarios vs $`T_0`$. The pion contribution is essentially absent for the higher $`M`$-values chosen. Note that the MT rate is now lower than the ST rate. The dashed lines represent the ratio of just the quark contributions. Quarks contribute fewer dileptons in the MT scenario, where the quark phase is shortened by the occurrence of the chiral transition.

Figure 4 The ratios of the MT rate to the rate if no transition occurred at all.


text:
# RELATIVISTIC GRAVITATIONAL COLLAPSE IN NON–COMOVING COORDINATES: THE POST-QUASISTATIC APPROXIMATION

## I Introduction

The problem of general relativistic gravitational collapse has attracted the attention of researchers since the seminal paper by Oppenheimer and Snyder . The motivation for such interest is easily understood: the gravitational collapse of massive stars represents one of the few observable phenomena, where general relativity is expected to play a relevant role. Ever since that work, much was written by researchers trying to provide models of evolving gravitating spheres. However this endeavour proved to be difficult and uncertain. Different kinds of advantages and obstacles appear, depending on the approach adopted for the modeling.

Thus, numerical methods (see  and references therein) are enabling researchers to investigate systems which are extremely difficult to handle analytically. In the case of General Relativity, numerical models have proved valuable for investigations of strong field scenarios and have been crucial to reveal unexpected phenomena . Even specific difficulties associated with numerical solutions of partial differential equations in presence of shocks are being overpassed . By these days what seems to be the main limitation for numerical relativity is the computational demands for 3D evolution, prohibitive in some cases . Nevertheless, purely numerical solutions usually hinder to catch general, qualitative, aspects of the process.

On the other hand, analytical solutions although more suitable for a general discussion (see  and references therein), are found, either for too simplistic equations of state and/or under additional heuristic assumptions whose justification is usually uncertain.

Therefore it seems useful to consider nonstatic models which are relatively simple to analyze but still contain some of the essential features of a realistic situation.

Accordingly, it is our purpose in this work to present an approach for modeling the evolution of self–gravitating spheres, which may be regarded as a “compromise” between the two approaches mentioned above (analytical and numerical).

Indeed, the proposed method, starting from any interior (analytical) static spherically symmetric (“seed”) solution to Einstein equations, leads to a system of ordinary differential equations for quantities evaluated at the boundary surface of the fluid distribution, whose solution (numerical), allows for modeling the dynamics of self–gravitating spheres, whose static limit is the original “seed” solution.

The approach is based on the introduction of a set of conveniently defined “effective” variables, which are effective pressure and energy density, and an heuristic ansatzs on the latter , whose rationale and justification become intelligible within the context of the post–quasistatic appproximation defined below. In the quasistatic approximation (see next section), the effective variables coincide with the corresponding physical variables (pressure and density) and therefore the method may be regarded as an iterative method with each consecutive step corresponding to a stronger departure from equilibrium. In this work we shall restrain ourselves to the post–quasistatic level (see next section for details).

At this point it is important to stress a crucial difference between this method and the one proposed many years ago with a similar structure, but based on radiative Bondi coordinates (see  and references therein): in the latter the introduced effective variables do not coincide with the corresponding physical variables in the quasistatic approximation (they do coincide in the static limit), and accordingly the ansatz on those variables remains as an heuristic assumption, only justified by the eventual suitability of the obtained models.

The fluid distribution under consideration will be assumed to be dissipative. Indeed, dissipation due to the emission of massless particles (photons and/or neutrinos) is a characteristic process in the evolution of massive stars. In fact, it seems that the only plausible mechanism to carry away the bulk of the binding energy of the collapsing star, leading to a neutron star or black hole is neutrino emission . Consequently, in this paper, the matter distribution forming the selfgravitating object will be described as a dissipative fluid.

In the diffusion approximation, it is assumed that the energy flux of radiation (as that of thermal conduction) is proportional to the gradient of temperature. This assumption is in general very sensible, since the mean free path of particles responsibles for the propagation of energy in stellar interiors is in general very small as compared with the typical length of the object. Thus, for a main sequence star as the sun, the mean free path of photons at the centre, is of the order of $`2cm`$. Also, the mean free path of trapped neutrinos in compact cores of densities about $`10^{12}g.cm.^3`$ becomes smaller than the size of the stellar core .

Furthermore, the observational data collected from supernovae 1987A indicates that the regime of radiation transport prevailing during the emission process, is closer to the diffusion approximation than to the streaming out limit .

However in many other circumstances, the mean free path of particles transporting energy may be large enough as to justify the free streaming approximation. Therefore our formalism will include simultaneously both limiting cases of radiative transport (diffusion and streaming out), allowing for describing a wide range situations.

Besides the usual physical variables (energy density, pressure, velocity, heat flow, etc.) we shall also incorporate into discussion other quantities which are expected to play an important role in the evolution of evolving self–gravitating systems, such as the Weyl tensor, the shear of the fluid and the Tolman mass. Therefore these quantities will be calculated and used in the process of modeling. It is also worth mentioning that although the most common method of solving Einstein’s equations is to use commoving coordinates (e.g. , ), we shall use noncomoving coordinates, which implies that the velocity of any fluid element (defined with respect to a conveniently chosen set of observers) has to be considered as a relevant physical variable ().

The plan of the paper is as follows. In Section 2 we define the conventions and give the field equations and expressions for the kinematical and physical variables we shall use, in noncomoving coordinates. The proposed approach is presented and explained in Section 3. In Section 4 we illustrate the method by means of three examples. Finally a discussion of results is presented in Section 5.

## II Relevant Equations and Conventions

### A The field equations

We consider spherically symmetric distributions of collapsing fluid, which for sake of completeness we assume to be locally anisotropic, undergoing dissipation in the form of heat flow and/or free streaming radiation, bounded by a spherical surface $`\mathrm{\Sigma }`$.

The line element is given in Schwarzschild–like coordinates by

$$ds^2=e^\nu dt^2e^\lambda dr^2r^2\left(d\theta ^2+sin^2\theta d\varphi ^2\right),$$
(1)
where $`\nu (t,r)`$ and $`\lambda (t,r)`$ are functions of their arguments. We number the coordinates: $`x^0=t;x^1=r;x^2=\theta ;x^3=\varphi `$.

The metric (1) has to satisfy Einstein field equations

$$G_\mu ^\nu =8\pi T_\mu ^\nu ,$$
(2)
which in our case read :

$$8\pi T_0^0=\frac{1}{r^2}+e^\lambda \left(\frac{1}{r^2}\frac{\lambda ^{}}{r}\right),$$
(3)
$$8\pi T_1^1=\frac{1}{r^2}+e^\lambda \left(\frac{1}{r^2}+\frac{\nu ^{}}{r}\right),$$
(4)
$`8\pi T_2^2=8\pi T_3^3=`$ $``$ $`{\displaystyle \frac{e^\nu }{4}}\left(2\ddot{\lambda }+\dot{\lambda }(\dot{\lambda }\dot{\nu })\right)`$ (5)
$`+`$ $`{\displaystyle \frac{e^\lambda }{4}}\left(2\nu ^{\prime \prime }+\nu ^2\lambda ^{}\nu ^{}+2{\displaystyle \frac{\nu ^{}\lambda ^{}}{r}}\right),`$ (6)
$$8\pi T_{01}=\frac{\dot{\lambda }}{r},$$
(7)
where dots and primes stand for partial differentiation with respect to $`t`$ and $`r`$, respectively.

In order to give physical significance to the $`T_\nu ^\mu `$ components we apply the Bondi approach .

Thus, following Bondi, let us introduce purely locally Minkowski coordinates ($`\tau ,x,y,z`$)

$$d\tau =e^{\nu /2}dt;dx=e^{\lambda /2}dr;dy=rd\theta ;dz=rsin\theta d\varphi .$$

Then, denoting the Minkowski components of the energy tensor by a bar, we have

$$\overline{T}_0^0=T_0^0;\overline{T}_1^1=T_1^1;\overline{T}_2^2=T_2^2;\overline{T}_3^3=T_3^3;\overline{T}_{01}=e^{(\nu +\lambda )/2}T_{01}.$$

Next, we suppose that when viewed by an observer moving relative to these coordinates with proper velocity $`\omega `$ in the radial direction, the physical content of space consists of an anisotropic fluid of energy density $`\rho `$, radial pressure $`P_r`$, tangential pressure $`P_{}`$, radial heat flux $`\widehat{q}`$ and unpolarized radiation of energy density $`\widehat{ϵ}`$ traveling in the radial direction. Thus, when viewed by this moving observer the covariant tensor in Minkowski coordinates is

$`\left(\begin{array}{cccc}\rho +\widehat{ϵ}& \widehat{q}\widehat{ϵ}& 0& 0\\ \widehat{q}\widehat{ϵ}& P_r+\widehat{ϵ}& 0& 0\\ 0& 0& P_{}& 0\\ 0& 0& 0& P_{}\end{array}\right).`$

Then a Lorentz transformation readily shows that

$$T_0^0=\overline{T}_0^0=\frac{\rho +P_r\omega ^2}{1\omega ^2}+\frac{2Q\omega e^{\lambda /2}}{(1\omega ^2)^{1/2}}+ϵ,$$
(8)
$$T_1^1=\overline{T}_1^1=\frac{P_r+\rho \omega ^2}{1\omega ^2}\frac{2Q\omega e^{\lambda /2}}{(1\omega ^2)^{1/2}}ϵ,$$
(9)
$$T_2^2=T_3^3=\overline{T}_2^2=\overline{T}_3^3=P_{},$$
(10)
$$T_{01}=e^{(\nu +\lambda )/2}\overline{T}_{01}=\frac{(\rho +P_r)\omega e^{(\nu +\lambda )/2}}{1\omega ^2}\frac{Qe^{\nu /2}e^\lambda }{(1\omega ^2)^{1/2}}(1+\omega ^2)e^{(\nu +\lambda )/2}ϵ,$$
(11)
with

$$Q\frac{\widehat{q}e^{\lambda /2}}{(1\omega ^2)^{1/2}}$$
(12)
and

$$ϵ\widehat{ϵ}\frac{(1+\omega )}{(1\omega )}.$$
(13)
Note that the coordinate velocity in the ($`t,r,\theta ,\varphi `$) system, $`dr/dt`$, is related to $`\omega `$ by

$$\omega =\frac{dr}{dt}e^{(\lambda \nu )/2}.$$
(14)
Feeding back (811) into (37), we get the field equations in the form

$$\frac{\rho +P_r\omega ^2}{1\omega ^2}+\frac{2Q\omega e^{\lambda /2}}{(1\omega ^2)^{1/2}}+ϵ=\frac{1}{8\pi }\left\{\frac{1}{r^2}+e^\lambda \left(\frac{1}{r^2}\frac{\lambda ^{}}{r}\right)\right\},$$
(15)
$$\frac{P_r+\rho \omega ^2}{1\omega ^2}+\frac{2Q\omega e^{\lambda /2}}{(1\omega ^2)^{1/2}}+ϵ=\frac{1}{8\pi }\left\{\frac{1}{r^2}e^\lambda \left(\frac{1}{r^2}+\frac{\nu ^{}}{r}\right)\right\},$$
(16)
$`P_{}={\displaystyle \frac{1}{8\pi }}\{{\displaystyle \frac{e^\nu }{4}}(2\ddot{\lambda }+\dot{\lambda }(\dot{\lambda }\dot{\nu }))`$ (17)
$`{\displaystyle \frac{e^\lambda }{4}}(2\nu ^{\prime \prime }+\nu ^2\lambda ^{}\nu ^{}+2{\displaystyle \frac{\nu ^{}\lambda ^{}}{r}})\},`$ (18)
$$\frac{(\rho +P_r)\omega e^{(\nu +\lambda )/2}}{1\omega ^2}+\frac{Qe^{\nu /2}e^\lambda }{(1\omega ^2)^{1/2}}(1+\omega ^2)+e^{(\nu +\lambda )/2}ϵ=\frac{\dot{\lambda }}{8\pi r}.$$
(19)
Observe that if $`\nu `$ and $`\lambda `$ are fully specified, then (1519) becomes a system of algebraic equations for the physical variables $`\rho `$, $`P_r`$, $`P_{}`$, $`\omega `$, $`Q`$ and $`ϵ`$. Obviuosly, in the most general case when all these variables are non–vanishing, the system is underdetermined, and two equations of state should be given. In general, whenever $`Q0`$ a transport equation has to be assumed. In the case originally considered by Bondi  (locally isotropic fluid and free streaming regime, $`Q=0`$) the system is closed. For the adiabatic ($`ϵ=Q=0`$), and locally isotropic fluid ($`P_r=P_{}`$) the system is overdetermined, and a constraint on the physical variables appears.

At the outside of the fluid distribution, the spacetime is that of Vaidya, given by

$$ds^2=\left(1\frac{2M(u)}{}\right)du^2+2dud^2\left(d\theta ^2+sin^2\theta d\varphi ^2\right),$$
(20)
where $`u`$ is a coordinate related to the retarded time, such that $`u=constant`$ is (asymptotically) a null cone open to the future and $``$ is a null coordinate ($`g_{}=0`$). It should be remarked, however, that strictly speaking, the radiation can be considered in radial free streaming only at radial infinity.

The two coordinate systems ($`t,r,\theta ,\varphi `$) and ($`u,,\theta ,\varphi `$) are related at the boundary surface and outside it by

$$u=tr2Mln\left(\frac{r}{2M}1\right),$$
(21)
$$=r.$$
(22)
In order to match smoothly the two metrics above on the boundary surface $`r=r_\mathrm{\Sigma }(t)`$, we first require the continuity of the first fundamental form across that surface. Then

$$\left[e^{\nu _\mathrm{\Sigma }}e^{\lambda _\mathrm{\Sigma }}\dot{r}_\mathrm{\Sigma }^2\right]dt^2=\left[1\frac{2M}{R_\mathrm{\Sigma }}+2\frac{dR_\mathrm{\Sigma }}{du}\right]du^2,$$
(23)
where $`R=R_\mathrm{\Sigma }(u)`$ is the equation of the boundary surface in ($`u,R,\theta ,\varphi `$) coordinates.

From (23), using (14), (21) and (20) it follows

$$e^{\nu _\mathrm{\Sigma }}=1\frac{2M}{R_\mathrm{\Sigma }},$$
(24)
$$e^{\lambda _\mathrm{\Sigma }}=1\frac{2M}{R_\mathrm{\Sigma }}.$$
(25)
Where, from now on, subscript $`\mathrm{\Sigma }`$ indicates that the quantity is evaluated at the boundary surface $`\mathrm{\Sigma }`$.

Next, the unit vector $`n_\mu `$, normal to the boundary surface, has components

$$n_\mu ^{(+)}=(\beta \frac{dR_\mathrm{\Sigma }}{du},\beta ,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}),$$
(26)
where $`+`$ indicates that the components are evaluated from the outside of $`\mathrm{\Sigma }`$, and $`\beta `$ is given by

$$\beta =\left(1\frac{2M(u)}{R_\mathrm{\Sigma }}+2\frac{dR_\mathrm{\Sigma }}{du}\right)^{1/2}.$$
(27)
The unit vector normal to $`\mathrm{\Sigma }`$, evaluated from the inside, is given by

$$n_\mu ^{()}=(\dot{r}_\mathrm{\Sigma }\gamma ,\gamma ,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}),$$
(28)
with

$$\gamma =\left(e^{\lambda _\mathrm{\Sigma }}\dot{r}_\mathrm{\Sigma }^2e^{\nu _\mathrm{\Sigma }}\right)^{1/2}.$$
(29)
Let us now define a time–like vector $`v^\mu `$ such that

$$v^{\mu (+)}=\beta \delta _u^\mu +\beta \frac{dR_\mathrm{\Sigma }}{du}\delta _R^\mu $$
(30)
and

$$v^{\mu ()}=\frac{e^{\nu _\mathrm{\Sigma }/2}}{\left(1\omega _\mathrm{\Sigma }^2\right)^{1/2}}\delta _t^\mu +\frac{\omega _\mathrm{\Sigma }e^{\lambda _\mathrm{\Sigma }/2}}{\left(1\omega _\mathrm{\Sigma }^2\right)^{1/2}}\delta _r^\mu .$$
(31)
Then, junction conditions across $`\mathrm{\Sigma }`$, require (besides (23))

$$\left(T_{\mu \nu }n^\mu n^\nu \right)_\mathrm{\Sigma }^{(+)}=\left(T_{\mu \nu }n^\mu n^\nu \right)_\mathrm{\Sigma }^{()},$$
(32)
$$\left(T_{\mu \nu }n^\mu v^\nu \right)_\mathrm{\Sigma }^{(+)}=\left(T_{\mu \nu }n^\mu v^\nu \right)_\mathrm{\Sigma }^{()},$$
(33)
where the expressions for the energy momentum tensor at both sides of the boundary surface are

$$T_{\mu \nu }^{()}=\left(\rho +P_{}\right)u_\mu u_\nu P_{}g_{\mu \nu }+\left(P_rP_{}\right)s_\mu s_\nu +q_\mu u_\nu +q_\nu u_\mu +ϵl_\nu l_\mu $$
(34)
and

$$T_{\mu \nu }^{(+)}=\frac{1}{4\pi R^2}\frac{dM}{du}\delta _\mu ^0\delta _\nu ^0,$$
(35)
with

$$u^\mu =(\frac{e^{\nu /2}}{\left(1\omega ^2\right)^{1/2}},\frac{\omega e^{\lambda /2}}{\left(1\omega ^2\right)^{1/2}},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}),$$
(36)
$$s^\mu =(\frac{\omega e^{\nu /2}}{\left(1\omega ^2\right)^{1/2}},\frac{e^{\lambda /2}}{\left(1\omega ^2\right)^{1/2}},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}),$$
(37)
$$l^\mu =(e^{\nu /2},e^{\lambda /2},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}),$$
(38)
where $`u^\mu `$ denotes the four velocity of the fluid, $`s^\mu `$ is a radially directed space–like vector orthogonal to $`u^\mu `$, $`l^\mu `$ is a null outgoing vector, and

$$q^\mu =Q(\omega e^{(\lambda \nu )/2},\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}).$$
(39)
Then it follows from (32) and (33)

$$[P_r+\widehat{ϵ}]_\mathrm{\Sigma }=\left[\frac{1}{4\pi R^2}\frac{dM}{du}\beta ^2\right]_\mathrm{\Sigma },$$
(40)
$$[Qe^{\lambda /2}(1\omega ^2)^{1/2}+\widehat{ϵ}]_\mathrm{\Sigma }=\left[\frac{1}{4\pi R^2}\frac{dM}{du}\beta ^2\right]_\mathrm{\Sigma }.$$
(41)
Eqs. (23), (40) and (41) are the necessary and sufficient conditions for a smooth matching of the two metrics (1) and (20) on $`\mathrm{\Sigma }`$. Combining (40) and (41) we get

$$\left[P_r\right]_\mathrm{\Sigma }=\left[Qe^{\lambda /2}\left(1\omega ^2\right)^{1/2}\right]_\mathrm{\Sigma },$$
(42)
expressing the discontinuity of the radial pressure in the presence of heat flow, which is a well known result .

Next, it will be useful to calculate the radial component of the conservation law

$$T_{\nu ;\mu }^\mu =0.$$
(43)
After tedious but simple calculations we get

$$\left(8\pi T_1^1\right)^{}=\frac{16\pi }{r}\left(T_1^1T_2^2\right)+4\pi \nu ^{}\left(T_1^1T_0^0\right)+\frac{e^\nu }{r}\left(\ddot{\lambda }+\frac{\dot{\lambda }^2}{2}\frac{\dot{\lambda }\dot{\nu }}{2}\right),$$
(44)
which in the static case becomes

$$P_r^{}=\frac{\nu ^{}}{2}\left(\rho +P_r\right)+\frac{2\left(P_{}P_r\right)}{r},$$
(45)
representing the generalization of the Tolman–Oppenheimer–Volkof equation for anisotropic fluids .

### B The kinematical variables

The components of the shear tensor are defined by

$$\sigma _{\mu \nu }=u_{\mu ;\nu }+u_{\nu ;\mu }u_\mu a_\nu u_\nu a_\mu \frac{2}{3}\mathrm{\Theta }P_{\mu \nu },$$
(46)
where

$$P_{\mu \nu }=g_{\mu \nu }u{}_{\mu }{}^{}u_{\nu }^{};\mathrm{\Theta }=u_{;\mu }^\mu ;a_\mu =u^\nu u_{\mu ;\nu },$$
(47)
denote the projector onto the three space orthogonal to $`u^\mu `$, the expansion and the four acceleration, respectively.

A simple calculation gives

$$\mathrm{\Theta }=\frac{e^{\nu /2}}{2\left(1\omega ^2\right)^{1/2}}\left(\dot{\lambda }+\frac{2\omega \dot{\omega }}{1\omega ^2}\right)+\frac{e^{\lambda /2}}{2\left(1\omega ^2\right)^{1/2}}\left(\omega \nu ^{}+2\omega ^{}+\frac{2\omega ^2\omega ^{}}{1\omega ^2}+\frac{4\omega }{r}\right),$$
(48)
$$\sigma _{11}=\frac{2}{3\left(1\omega ^2\right)^{3/2}}\left[e^\lambda e^{\nu /2}\left(\dot{\lambda }+\frac{2\omega \dot{\omega }}{1\omega ^2}\right)+e^{\lambda /2}\left(\omega \nu ^{}+\frac{2\omega ^{}}{1\omega ^2}\frac{2\omega }{r}\right)\right],$$
(49)
$$\sigma _{22}=\frac{e^\lambda r^2\left(1\omega ^2\right)}{2}\sigma _{11},$$
(50)
$$\sigma _{33}=\frac{e^\lambda r^2\left(1\omega ^2\right)}{2}\mathrm{sin}^2\theta \sigma _{11},$$
(51)
$$\sigma _{00}=\omega ^2e^\lambda e^\nu \sigma _{11},$$
(52)
$$\sigma _{01}=\omega e^{\left(\nu \lambda \right)/2}\sigma _{11},$$
(53)
$$a_0=\frac{1}{1\omega ^2}\left[\left(\frac{\omega \dot{\omega }}{1\omega ^2}+\frac{\omega ^2\dot{\lambda }}{2}\right)+e^{\nu /2}e^{\lambda /2}\left(\frac{\omega \nu ^{}}{2}+\frac{\omega ^2\omega ^{}}{1\omega ^2}\right)\right],$$
(54)
$$a_1=\frac{1}{1\omega ^2}\left[\left(\frac{\omega \omega ^{}}{1\omega ^2}+\frac{\nu ^{}}{2}\right)+e^{\nu /2}e^{\lambda /2}\left(\frac{\omega \dot{\lambda }}{2}+\frac{\dot{\omega }}{1\omega ^2}\right)\right],$$
(55)
and for the shear scalar $`\sigma `$

$$\sigma =\sqrt{3}\left(\frac{\mathrm{\Theta }}{3}\frac{e^{\lambda /2}}{r}\frac{\omega }{\sqrt{1\omega ^2}}\right).$$
(56)
### C The Tolman mass

The Tolman mass for a spherically symmetric distribution of matter is given by (eq.(24) in )

$`m_T=`$ $`4\pi {\displaystyle _0^{r_\mathrm{\Sigma }}}r^2e^{(\nu +\lambda )/2}\left(T_0^0T_1^12T_2^2\right)𝑑r`$ (57)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^{r_\mathrm{\Sigma }}}r^2e^{(\nu +\lambda )/2}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{\mathrm{\pounds }}{\left[\left(g^{\alpha \beta }\sqrt{g}\right)/t\right]}}\right)g^{\alpha \beta }𝑑r,`$ (58)
where $`\mathrm{\pounds }`$ denotes the usual gravitational lagrangian density (eq.(10) in ). Although Tolman’s formula was introduced as a measure of the total energy of the system, with no commitment to its localization, we shall define the mass within a sphere of radius $`r`$, completely inside $`\mathrm{\Sigma }`$, as

$`m_T=`$ $`4\pi {\displaystyle _0^r}r^2e^{(\nu +\lambda )/2}\left(T_0^0T_1^12T_2^2\right)𝑑r`$ (59)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^r}r^2e^{(\nu +\lambda )/2}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{\mathrm{\pounds }}{\left[\left(g^{\alpha \beta }\sqrt{g}\right)/t\right]}}\right)g^{\alpha \beta }𝑑r.`$ (60)
This extension of the global concept of energy to a local level  is suggested by the conspicuous role played by $`m_T`$ as the “effective gravitational mass”, which will be exhibited below. Even though Tolman’s definition is not without its problems , we shall see that $`m_T`$, as defined by (60), is a good measure of the active gravitational mass, at least for the systems under consideration.

Let us now evaluate expression (60). The first integral in that expression

$$I4\pi _0^rr^2e^{(\nu +\lambda )/2}\left(T_0^0T_1^12T_2^2\right)𝑑r,$$
(61)
may be transformed to give (see  for details)

$`I=`$ $`e^{(\nu +\lambda )/2}\left[m(r,t){\displaystyle \frac{4\pi }{3}}r^3T_1^1\right]`$ (62)
$``$ $`{\displaystyle _0^r}e^{(\lambda \nu )/2}{\displaystyle \frac{r^2}{2}}\left(\ddot{\lambda }+{\displaystyle \frac{\dot{\lambda }^2}{2}}{\displaystyle \frac{\dot{\lambda }\dot{\nu }}{2}}\right)𝑑r,`$ (63)
where the mass function $`m`$, as usually is defined by

$$e^{\lambda (r,t)}=12m(r,t)/r.$$
(64)
Next, from (eq.(13) in )

$$\frac{}{t}\left(\frac{\mathrm{\pounds }}{\left[\left(g^{\alpha \beta }\sqrt{g}\right)/t\right]}\right)=\mathrm{\Gamma }_{\alpha \beta }^0+\frac{1}{2}\delta _\alpha ^0\mathrm{\Gamma }_{\beta \sigma }^\sigma +\frac{1}{2}\delta _\beta ^0\mathrm{\Gamma }_{\alpha \sigma }^\sigma ,$$
(65)
and so the second integral $`(II)`$ in (60) may be expressed as

$$II=\frac{1}{2}_0^rr^2e^{(\lambda \nu )/2}\left(\ddot{\lambda }+\frac{\dot{\lambda }^2}{2}\frac{\dot{\lambda }\dot{\nu }}{2}\right)𝑑r.$$
(66)
Thus

$$m_TI+II=e^{(\nu +\lambda )/2}\left[m(r,t)4\pi r^3T_1^1\right].$$
(67)
This is, formally, the same expression for $`m_T`$ in terms of $`m`$ and $`T_1^1`$ that appears in the static (or quasi-static) case (eq.(25) in ).

Replacing $`T_1^1`$ by (4) and $`m`$ by (64), one also finds

$$m_T=e^{(\nu \lambda )/2}\nu ^{}\frac{r^2}{2}.$$
(68)
This last equation brings out the physical meaning of $`m_T`$ as the active gravitational mass. Indeed, it can be easily shown  that the gravitational acceleration $`a`$ of a test particle, instantaneously at rest in a static gravitational field, as measured with standard rods and coordinate clock is given by

$$a=\frac{e^{(\nu \lambda )/2}\nu ^{}}{2}=\frac{m_T}{r^2}.$$
(69)
A similar conclusion can be obtained by inspection of Eq. (45) (valid only in the static or quasi-static case) . In fact, the first term on the right side of this equation (the “gravitational force” term) is a product of the “passive” gravitational mass density $`(\rho +P_r)`$ and a term proportional to $`m_T/r^2`$.

### D The Weyl tensor

Since the publication of Penrose‘s work , there has been an increasing interest in studying the possible role of Weyl tensor (or some function of it) in the evolution of self-gravitating systems . This interest is reinforced by the fact that for spherically symmetric distribution of fluid, the Weyl tensor may be defined exclusively in terms of the density contrast and the local anisotropy of the pressure (see below), which in turn are known to affect the fate of gravitational collapse .

Now, using Maple V, it is found that all non–vanishing components of the Weyl tensor are proportional to

$`W{\displaystyle \frac{r}{2}}C_{232}^3`$ $`=`$ $`W_{(s)}+{\displaystyle \frac{r^3e^\nu }{12}}\left(\ddot{\lambda }+{\displaystyle \frac{\dot{\lambda }^2}{2}}{\displaystyle \frac{\dot{\lambda }\dot{\nu }}{2}}\right)`$ (70)
where

$$W_{(s)}=\frac{r^3e^\lambda }{6}\left(\frac{e^\lambda }{r^2}\frac{1}{r^2}+\frac{\nu ^{}\lambda ^{}}{4}\frac{\nu ^2}{4}\frac{\nu ^{\prime \prime }}{2}\frac{\lambda ^{}}{2r}+\frac{\nu ^{}}{2r}\right),$$
(71)
corresponds to the contribution in the static case.

Also, from the field equations and the definition of the Weyl tensor it can be easily shown that (see  for details)

$$W=\frac{4\pi }{3}_0^rr^3\left(T_0^0\right)^{}𝑑r+\frac{4\pi }{3}r^3\left(T_2^2T_1^1\right).$$
(72)
## III The method

We have now available all the ingredients required to present our method, however before doing so some general considerations will be necessary.

### A Equilibrium and departures from equilibrium

The simplest situation, when dealing with self–gravitating spheres, is that of equilibrium (static case). In our notation that means that $`\omega =ϵ=Q=0`$, all time derivatives vanishes, and we obtain the generalized Tolman–Oppenheimer–Volkof equation (45).

Next, we have the quasistatic regime. By this we mean that the sphere changes slowly, on a time scale that is very long compared to the typical time in which the sphere reacts to a slight perturbation of hydrostatic equilibrium, this typical time scale is called hydrostatic time scale  (sometimes this time scale is also referred to as dynamical time scale, e.g. see the third reference in ). Thus, in this regime the system is always very close to hydrostatic equilibrium and its evolution may be regarded as a sequence of static models linked by (19). This assumption is very sensible because the hydrostatic time scale is very small for many phases of the life of the star. It is of the order of $`27`$ minutes for the Sun, $`4.5`$ seconds for a white dwarf and $`10^4`$ seconds for a neutron star of one solar mass and $`10`$ Km radius. It is well known that any of the stellar configurations mentioned above, generally, change on a time scale that is very long compared to their respective hydrostatic time scales. Let us now translate this assumption in conditions to $`\omega `$ and metric functions.

First of all, slow contraction means that the radial velocity $`\omega `$ as measured by the Minkowski observer is always much smaller than the velocity of light ($`\omega 1`$). Therefore we have to neglect terms of the order $`O(\omega ^2)`$.

Then (44) yields

$$\ddot{\lambda }+\frac{\dot{\lambda }^2}{2}\frac{\dot{\nu }\dot{\lambda }}{2}=8\pi re^\nu \left[(P_r+ϵ)^{}+\left(\rho +P_r+2ϵ\right)\frac{\nu ^{}}{2}2\frac{P_{}P_rϵ}{r}\right].$$
(73)
(observe the contribution of $`ϵ`$ to, both, $`P_r`$ and $`\rho `$ and the fact that $`ϵ`$, $`\omega `$ and $`Q`$ are of the same order of smallness, in this approximation).

Since by assumption, in this regime the system is always (not only at a given time $`t`$) in equilibrium (or very close to), (45) and (73) imply then, for an arbitrary slowly evolving configuration

$$\ddot{\lambda }\dot{\nu }\dot{\lambda }\dot{\lambda }^20,$$
(74)
and of course, time derivatives of any order of the left hand side of the hydrostatic equilibrium equation must also vanish, for otherwise the system will deviate from equilibrium. This condition implies, in particular, that we must demand in this regime

$$\ddot{\nu }0.$$

Finally, from the time derivative of (7), and using (11), it follows that

$$\dot{\omega }O(\ddot{\lambda },\dot{\lambda }\omega ,\dot{\nu }\omega ).$$
(75)
which implies that we have also to neglect terms linear in the acceleration. On purely physical considerations, it is obvious that the vanishing of $`\dot{\omega }`$ is required to keep the system always in equilibrium.

Thus, in the quasistatic regime we have to assume

$$O(\omega ^2)=\dot{\lambda }^2=\dot{\nu }^2=\dot{\lambda }\dot{\nu }=\ddot{\lambda }=\ddot{\nu }=0,$$
(76)
implying that the system remains in (or very close to) equilibrium. However, during their evolution, self–gravitating objects may pass through phases of intense dynamical activity, with time scales of the order of magnitude of (or even smaller than) the hydrostatic time scale, and for which the quasi–static approximation is clearly not reliable (e.g.,the collapse of very massive stars  and the quick collapse phase preceding neutron star formation, see for example  and references therein). In these cases it is mandatory to take into account terms which describe departure from equilibrium.

### B The effective variables and the post–quasistatic approximation

Let us now define the following effective variables:

$$\stackrel{~}{\rho }=T_0^0=\frac{\rho +P_r\omega ^2}{1\omega ^2}+\frac{2Q\omega e^{\lambda /2}}{(1\omega ^2)^{1/2}}+ϵ,$$
(77)
$$\stackrel{~}{P}=T_1^1=\frac{P_r+\rho \omega ^2}{1\omega ^2}+\frac{2Q\omega e^{\lambda /2}}{(1\omega ^2)^{1/2}}+ϵ.$$
(78)
In the quasistatic regime the effective variables satisfy the same equation (45) as the corresponding physical variables (taking into account the contribution of $`ϵ`$ to the “total” energy density and radial pressure, whenever the free streaming approximation is being used). Therefore in the quasistatic situation (and obviously in the static too), effective and physical variables share the same radial dependence. Next, feeding back (77) and (78) into (15) and (16), these two equations may be formally integrated , to obtain:

$$m=4\pi _0^rr^2\stackrel{~}{\rho }𝑑r,$$
(79)
$$\nu =\nu _\mathrm{\Sigma }+_{r_\mathrm{\Sigma }}^r\frac{2(4\pi r^3\stackrel{~}{P}+m)}{r(r2m)}𝑑r.$$
(80)
From where it is obvious that for a given radial dependence of the effective variables, the radial dependence of metric functions become completely determined.

With this last comment in mind, we shall define the post–quasistatic regime as that corresponding to a system out of equilibrium (or quasiequilibrium) but whose effective variables share the same radial dependence as the corrresponding physical variables in the state of equilibrium (or quasiequilibrium). Alternatively it may be said that the system in the post–quasistatic regime is characterized by metric functions whose radial dependence is the same as the metric functions corresponding to the static (quasistatic) regime. The rationale behind this definition is not difficult to grasp: we look for a regime which although out of equilibrium, represents the closest possible situation to a quasistatic evolution (see more on this point in the last Section).

### C The algorithm

Let us now outline the approach that we propose:

1. Take an interior solution to Einstein equations, representing a fluid distribution of matter in equilibrium, with a given

$$\rho _{st}=\rho (r)P_{rst}=P_r(r)$$
2. Assume that the $`r`$ dependence of $`\stackrel{~}{P}`$ and $`\stackrel{~}{\rho }`$ is the same as that of $`P_{rst}`$ and $`\rho _{st}`$, respectively.
3. Using equations (80) and (79), with the $`r`$ dependence of $`\stackrel{~}{P}`$ and $`\stackrel{~}{\rho }`$, one gets $`m`$ and $`\nu `$ up to some functions of $`t`$, which will be specified below.
4. For these functions of $`t`$ one has three ordinary differential equations (hereafter referred to as surface equations), namely:

1. Equation (14) evaluated on $`r=r_\mathrm{\Sigma }`$.
    2. Equation (44) evaluated on $`r=r_\mathrm{\Sigma }`$.
    3. The equation relating the total mass loss rate with the energy flux through the boundary surface.
5. Depending on the kind of matter under consideration, the system of surface equations described above may be closed with the additional information provided by the transport equation and/or the equation of state for the anisotropic pressure and/or additional information about some of the physical variables evaluated on the boundary surface (e.g. the luminosity).
6. Once the system of surface equations is closed, it may be integrated for any particular initial data.
7. Feeding back the result of integration in the expressions for $`m`$ and $`\nu `$, these two functions are completely determined.
8. With the input from the point 7 above, and using field equations, together with the equations of state and/or transport equation, all physical variables may be found for any piece of matter distribution.

### D The Surface equations

As it should be clear from the above the crucial point in the algorithm is the system of surface equations. So, let us specify them now.

Introducing the dimenssionles variables

$$A=r_\mathrm{\Sigma }/m_\mathrm{\Sigma }(0),$$

$$F=12M/A,$$

$$M=m_\mathrm{\Sigma }/m_\mathrm{\Sigma }(0),$$

$$\mathrm{\Omega }=\omega _\mathrm{\Sigma },$$

$$\alpha =t/m_\mathrm{\Sigma }(0),$$

where $`m_\mathrm{\Sigma }(0)`$ denotes the total initial mass, we obtain the first surface equation by evaluating (14) at $`r=r_\mathrm{\Sigma }`$, one gets

$$\frac{dA}{d\alpha }=F\mathrm{\Omega }.$$
(81)
Next, using junction conditions, one obtains from (64), (15) and (19) evaluating at $`r=r_\mathrm{\Sigma }`$, that

$$\frac{dM}{d\alpha }=F(1+\mathrm{\Omega })\widehat{E},$$
(82)
with

$$\widehat{E}=4\pi r_\mathrm{\Sigma }^2(\widehat{ϵ}_\mathrm{\Sigma }+\widehat{q}_\mathrm{\Sigma }),$$
(83)
where the first and second term on the right of (82), represent the gravitational redshift and the Doppler shift corrections, respectively.

Then, defining the luminosity perceived by an observer at infinity as

$$L=\frac{dM}{d\alpha }.$$

we obtain the second surface equation in the form

$$\frac{dF}{d\alpha }=\frac{F}{A}(1F)\mathrm{\Omega }+2L/A.$$
(84)
The third surface equation may be obtained by evaluating at the boundary surface, the conservation law $`T_{1;\mu }^\mu =0`$, which reads

$$\stackrel{~}{P}^{^{}}+\frac{(\stackrel{~}{\rho }+\stackrel{~}{P})(4\pi r^3\stackrel{~}{P}+m)}{r(r2m)}=\frac{e^\nu }{4\pi r(r2m)}\left(\ddot{m}+\frac{3\dot{m}^2}{r2m}\frac{\dot{m}\dot{\nu }}{2}\right)+\frac{2}{r}(P_{}\stackrel{~}{P}).$$
(85)
Now, in the following section we consider two relatively simple models with a separable effective density, i.e., $`\stackrel{~}{\rho }=f(t)h(r)`$; thus equation (85) evaluated at the boundary surface leads to

$$\frac{d\mathrm{\Omega }}{d\alpha }=\mathrm{\Omega }^2\left[\frac{8F}{A}+2Fk(r_\mathrm{\Sigma })+4\pi \stackrel{~}{\rho }_\mathrm{\Sigma }A(3\mathrm{\Omega }^2)\right]\frac{F}{\stackrel{~}{\rho }_\mathrm{\Sigma }}\left[R\frac{2}{A}\left(P_\mathrm{\Sigma }\stackrel{~}{\rho }_\mathrm{\Sigma }\mathrm{\Omega }^2\frac{\overline{E}(1+\mathrm{\Omega })}{4\pi r_\mathrm{\Sigma }^2}\right)\right],$$
(86)
where

$$R=\left[\stackrel{~}{P}^{^{}}+\frac{\stackrel{~}{P}+\stackrel{~}{\rho }}{12m/r}(4\pi r\stackrel{~}{P}+\frac{m}{r^2})\right]_\mathrm{\Sigma },$$
(87)
$$\overline{E}=\widehat{E}(1+\mathrm{\Omega }),$$
(88)
and

$$k(r_\mathrm{\Sigma })=\frac{d}{dr_\mathrm{\Sigma }}ln\left(\frac{1}{r_\mathrm{\Sigma }}_0^{r_\mathrm{\Sigma }}𝑑rr^2h(r)/h(r_\mathrm{\Sigma })\right).$$
(89)
Before analyzing specific models, some interesting conclusions can be obtained at this level of generality. One of these conclusions concerns the condition of bouncing at the surface which, of course, is related to the occurrence of a minimum radius $`A`$. According to (81) this requires $`\mathrm{\Omega }=0`$, and we have

$$\frac{d^2A}{d\alpha ^2}=F\frac{d\mathrm{\Omega }}{d\alpha },$$
(90)
or using (86)

$$\frac{d\mathrm{\Omega }}{d\alpha }(\mathrm{\Omega }=0)=\frac{F}{\stackrel{~}{\rho }_\mathrm{\Sigma }}\left[R\frac{2}{A}\left(P_\mathrm{\Sigma }\frac{\widehat{E}}{4\pi r_{\mathrm{\Sigma }}^{}{}_{}{}^{2}}\right)\right].$$
(91)
Observe that a positive energy flux ($`\widehat{E}`$) tends to decrease the radius of the sphere, i.e., it favors the compactification of the object, which is easily understandable. The same happens when $`R>0`$ or $`P_\mathrm{\Sigma }<0`$. The opposite effect occurs when these quantities have the opposite signs. Now, for a positive energy flux the sphere can only bounce at its surface when

$$\frac{d\mathrm{\Omega }}{d\alpha }(\mathrm{\Omega }=0)0.$$

According to (91) this is equivalent to

$$R(\mathrm{\Omega }=0)+\frac{2P_\mathrm{\Sigma }}{A}0.$$
(92)
A physical meaning can be associated to this equation as follows. For non–radiating, static configuration, $`R`$ as defined by (87) consists of two parts. The first term which together with $`[2(P_{}P_r)/r]_\mathrm{\Sigma }`$ represents the hydrodynamical force (see (45)) and the second which is of course the gravitational force. The resulting force in the sense of increasing $`r`$ is precisely $`R+[2(P_{}P_r)/r]_\mathrm{\Sigma }`$, if this is positive a net outward acceleration occurs and vice–versa. Equation (92) is the natural generalization of this result for general non–static configurations.

As mentioned before, besides the surface equations, in some cases (depending on the type of matter under consideration) further information has to be provided in the form of equation of state for the tangential stresses and/or transport equation. In the next Section we shall illustrate our method with three examples, one of which refers to an anisotropic fluid, and for which we shall further assume the equation of state (see ,)

$$P_{}P_r=\frac{C(\stackrel{~}{P}+\stackrel{~}{\rho })(4\pi r^3\stackrel{~}{P}+m)}{(r2m)},$$
(93)
where $`C`$ is a constant.

## IV Examples

The only purpose of the present Section is to illustrate the proposed method. For simplification we shall consider only the adiabatic case ($`ϵ=Q=0`$). For all these models we shall calculate the physical and geometrical variables for any piece of matter, as function of the timelike coordinate. In spite of the simplicity of the models, some interesting conclusions about the physical meaning of different variables may be reached.

One of the models has as the “seed” solution the well known Schwarzshild interior solution, whose properties have been extensively discussed in the literature. The second example is based on an anisotropic fluid without radial pressure. Models of this kind have also been discussed extensively since the original Einstein paper (see ). Finally, the third example represents the dynamic version of the Tolman VI static solution (), whose equation of state, as is well known, approaches that for highly compressed Fermi gas.

### A Schwarzschild–type model

This model is inspired in the well known interior Schwarzschild solution. Accordingly we take

$$\stackrel{~}{\rho }=f(t),$$
(94)
where $`f`$ is an arbitrary function of $`t`$. And the expression for $`\stackrel{~}{P}`$ is

$$\frac{\stackrel{~}{P}+\frac{1}{3}\stackrel{~}{\rho }}{\stackrel{~}{P}+\stackrel{~}{\rho }}=\left(1\frac{8\pi }{3}\stackrel{~}{\rho }r^2\right)^{1/2}k(t),$$
(95)
where $`k`$ is a function of $`t`$ to be determined from junction conditions (42), which in terms of effective variables, becomes

$$\stackrel{~}{P}_\mathrm{\Sigma }=\stackrel{~}{\rho }_\mathrm{\Sigma }\mathrm{\Omega }^2.$$
(96)
Thus, using (95) and (96) we have for the effective variables

$$\stackrel{~}{\rho }=\frac{3(1F)}{8\pi r_\mathrm{\Sigma }^2},$$
(97)
$$\stackrel{~}{P}=\frac{\stackrel{~}{\rho }}{3}\left\{\frac{\chi F^{1/2}3\psi \xi }{\psi \xi \chi F^{1/2}}\right\},$$
(98)
with

$$\xi =[1(1F)(r/r_\mathrm{\Sigma })^2]^{1/2}$$

and

$$\chi =3(\mathrm{\Omega }^2+1)(1F),$$

$$\psi =(3\mathrm{\Omega }^2+1)(1F).$$

And for the metric functions $`m`$ and $`\nu `$ we get, using (79) and (80)

$$m=m_\mathrm{\Sigma }(r/r_\mathrm{\Sigma })^3,$$
(99)
$$e^\nu =\left\{\frac{\chi F^{1/2}\psi \xi }{2(1F)}\right\}^2.$$
(100)
The third surface equation for this model becomes

$$\frac{d\mathrm{\Omega }}{d\alpha }=\frac{\mathrm{\Omega }^2}{2A}(73\mathrm{\Omega }^2+3F(\mathrm{\Omega }^21)).$$
(101)
This equation together with (81) and (84) form the set of surface equations for this model. We have integrated it numerically and from this integration, all physical variables are found for any piece of the fluid distribution, following the algorithm described above.

Figures (1)–(5) exhibit the behaviour of $`\rho `$, $`P`$, $`\omega `$, $`\sigma `$ and $`W`$ for an initially contracting configuration, as function of $`\alpha `$ and different pieces of matter.

Figure (6) shows the profile of $`\omega `$ as function of $`r/r_\mathrm{\Sigma }`$ for $`\alpha =10`$.

### B Lemaitre–Florides–type model

This model has as the “seed” solution a configuration with homogeneous energy density and vanishing radial pressure. Configurations of this kind were suggested by the first time by Lemaitre ().

The corresponding effective variables now are:

$$\stackrel{~}{\rho }=f(t)$$
(102)
and

$$\stackrel{~}{P}=0.$$
(103)
Observe that in this case, because of (102) and (103), it follows from (77) and (78) that the radial pressure is discontinuos at the boundary surface, with

$$P_{r\mathrm{\Sigma }}=3(1F)\mathrm{\Omega }^2/8\pi r_\mathrm{\Sigma }^2,$$
(104)
for otherwise either $`\rho _\mathrm{\Sigma }`$ or $`\mathrm{\Omega }`$ should vanish at $`\mathrm{\Sigma }`$. Therefore the only way to “dynamize” this model is by relaxing boundary conditions, allowing for the presence of a kind of surface tension.

Once the effective variables are defined, we only need the value of the tangential pressure at the boundary to close the system of surface equations. This is obtained by evaluating (93) at $`\mathrm{\Sigma }`$.

Next, following the algorithm, all physical variables may be found for any piece of material as functions of the timelike coordinate. Although we are not going to exhibit them here, because the graphics are not particularly illuminating, we wanted to present an example which, besides the fact that implies an anisotropic fluid, requires the introduction of a surface tension, for allowing the application of the algorithm.

### C Tolman VI–type model

Our last example is based on the Tolman VI solution. Accordingly the effective variables for this model will be

$$\stackrel{~}{\rho }=\frac{3g(t)}{r^2}$$
(105)
and

$$\stackrel{~}{P}=\frac{g(9bK(r/r_\mathrm{\Sigma }))}{(9b(r/r_\mathrm{\Sigma }))r^2},$$
(106)
where $`g`$ and $`b`$ are functions of $`\alpha `$, to be obtained from (96). Then,

$$\stackrel{~}{\rho }=\frac{3(1F)}{24\pi r^2}.$$
(107)
Using (79) and (80) we get

$$m=m_\mathrm{\Sigma }r/r_\mathrm{\Sigma },$$
(108)
$$\nu =lnF+\frac{8\pi g}{F}\left\{4ln(r/r_\mathrm{\Sigma })+8ln\left(\frac{b(r/r_\mathrm{\Sigma })K}{b9}\right)\right\}.$$
(109)
Finally, solving the surface equations for this model, $`m`$ and $`\nu `$ are completely determined and all physical variables can thereby be calculated. Besides the intrinsic physical interest of the equation of state of this “seed” model mentioned before, it is interesting because of the fact that the static limit of the model (unlike the previous ones) is “unstable”, in the sense that it requires a specific value of the gravitational potential at the boundary, namely $`m_\mathrm{\Sigma }(0)/r_\mathrm{\Sigma }=3/14`$. For values above (below) this, the sphere starts to collapse (expand).

Figures (7) and (8) display the evolution of velocity ($`\omega `$) for different regions of the sphere, and for initial values of $`F`$ corresponding to values of $`m_\mathrm{\Sigma }(0)/r_\mathrm{\Sigma }`$, above and below the equilibrium value, respectively. Figures (9) and (10), represent the evolution of the Weyl tensor ($`W`$), for some internal region and the boundary surface respectively, and initials values of $`F`$ corresponding to values of $`m_\mathrm{\Sigma }(0)/r_\mathrm{\Sigma }`$, above and below equilibrium. Finally, figures (11) and (12) exhibit the behaviour of the shear ($`\sigma `$) for different regions and initial values of $`F`$ corresponding to values of $`m_\mathrm{\Sigma }(0)/r_\mathrm{\Sigma }`$, above and below equilibrium.

We shall comment on these graphics in the next Section.

## V Conclusions

A method has been presented, which allows for the description of radiating selfgravitating relativistic spheres. In its most general form, the approach incorporates the two limiting cases of radiation transport (free streaming and diffusion) as well as the possibility of dealing with anisotropic fluids.

The cornerstone of the algorithm is an ansatz, based on an specific definition of the post–quasistatic approximation, namely: considering different degrees of departure from equilibrium, the post–quasistatic regime (i.e. the next step after the quasistatic situation) is defined as that, characterized by metric functions whose radial dependence is the same as that of the quasistatic regime. This in turn implies, that the effective variables defined above, share the same radial dependence as the correspondig physical variables of the quasistatic regime. The rationale behind this definition seems intelligible when it is remembered that in the latter case (the quasistatic) the effective variables share the same radial dependence as that of the physical variables in the static regime. Thus, starting with a static configuration, the first “level” off equilibrium, beyond the quasitatic situation, is represented by the post–quasistatic regime.

Once the static (“seed”) solution has been selected, then the definition of the effective variables together with surface equations, allows for determination of metric functions, which in turn lead to the full description of physical variables as functions of the timelike coordinate, for any region of the sphere. In this process, depending on the kind of matter and/or the prevailing transport approximation, additional equations of state and/or transport equations and/or some of the surface variables (e.g. the luminosity) have to be specified.

All physical variables having been found (particularly the energy density and the radial pressure) then we may, in principle, go to the next step assuming that the effective variables now share the same radial dependence as that of the physical variables just obtained. In this sense the algorithm may be regarded as an iterative approach. For obvious reasons we have restrained ourselves to the first step of the process. It remains to be seen if available physical evidence justifies to go through the complexities associated with the “post–post–quasitatic” approximation.

In order to illustrate the method, and without the pretension of modeling specific astrophysical scenarios, we have presented three examples, in the simplest (adiabatic) case.

In the first model, the profiles of the shear and the Weyl tensor, clearly illustrate the “dynamics” of the model, tending to zero in the static limit. The fact that these two quantities vanish in the quasitatic regime (for this specific model), brings out further their relevance in the treatment of situations off equilibrium. On the other hand however, the velocity profiles show almost no difference between the two regimes. Deviations from homology contraction due to relativistic gravitational effect are also indicated.

The purpose of the second example was to illustrate the implementation of the algorithm, for anisotropic fluids. The very particular form of the “seed” equation of state of this model, imposses discontinuity (surface tension) of radial pressure at the boundary. Of course such discontinuity vanishes at the static (or quasistatic) regime.

Finally, a model based on the Tolman VI solution was presented. This static solution, as was already mentioned, requires a specific value of $`m_\mathrm{\Sigma }(0)/r_\mathrm{\Sigma }`$, accordingly any deviation from this value leads to deviations from the static regime (observe that the quasitatic regime is incompatible with this solution). The velocity profiles indicate that all regions either expand or contract, and therefore, cracking (different signs of the velocity for different regions of the sphere) will not occur . This is consistent with the established fact that cracking only occurs for anisotropic fluids or isotropic fluids with outgoing radiation in the free streaming approximation.

Also, the profiles of the Weyl tensor and the shear, clearly diverging from the initial values as time proceeds and the evolution becomes more and more “dynamic”, stress once again their roles in describing departures off equilibrium.

###### Acknowledgements.

WB was benefited from research support from the Consejo de Investigación, Universidad de Oriente, under Grant CI-5-1002-1054/01 and from FONACIT under grant S1–98003270 to the Universidad de Oriente. NOS acknowledges financial support of CNPq. LH acknowledges financial assistance under grant BFM2000-1322 (M.C.T. Spain) and from Cátedra-FONACIT, under grant 2001001789.


text:
# 2D spectroscopy of double-barred galaxies

## 1. Double-barred galaxies: history of the problem.

Quarter of century ago, de Vaucouleurs (1975) found that an inner bar-like structure was nested in a large-scale bar at the optical image of NGC 1291. The first systematic observational study of this phenomenon was made by Buta & Crocker (1993) who published a list of 13 galaxies, in which the secondary bar was arbitrarily oriented in relation to the primary one. In the following years the isophotal analysis technique had allowed to detect secondary bars at the optical (Wozniak et al., 1995; Erwin & Sparke, 1998) and at the NIR images of barred galaxies (Friedli et al., 1996; Laine et al., 2001). Up to now, we have found, in total, 71 candidates to double-barred galaxies in the literature (Moiseev, 2001). Also numerous theoretical papers exist, but secondary bar dynamics is still far from being understood. Shlosman, Frank & Begeman (1989) have shown that a new circumnuclear bar may be formed as a result of gas instabilities within a primary bar. Maciejewski & Sparke (2000) found some families of closed orbital loops which can support both bars. Dynamics of stellar-gaseous secondary bars have been studied in numerical simulations (Pfenninger & Norman, 1990; Friedli & Martinet, 1993). New hydrodynamical simulations of the gas behaviour in double bars are presented by Heller, Shlosman & Englmaier (2001), Maciejewski et al. (2002) and Shlosman & Heller (2001). To test these theoretical predictions, new observational data is required.

In this paper we try to find an answer to the question – Are the secondary bars dynamically decoupled systems?

## 2. Necessity of 2D spectroscopy

Mostly of the information about double-barred objects has been obtained from galaxy imaging (optical or NIR). However the photometric structure revealed for them can be explained in a less exotic manner, without involvement of the secondary bars. For example, the following explanations of the turn of the inner isophote major axis are possible:

* a dynamically independent secondary bar ($`\mathrm{\Omega }_P\mathrm{\Omega }_S`$, where $`\mathrm{\Omega }_P`$ and $`\mathrm{\Omega }_S`$ are angular rotation speeds of the primary and secondary bars respectively);
* dynamically coupled bars $`(\mathrm{\Omega }_P=\mathrm{\Omega }_S)`$, or $`x_2`$-orbits between two inner Lindblad resonances (ILRs) in the primary bar;
* an elliptical ring in the plane of the disk, at the ILR resonance of the primary bar;
* a polar disk (or an inner polar ring) within the primary bar (similar structures have been found in a number of ordinary early-type disk galaxies, see Afanasiev & Sil’chenko, 2000; Sil’chenko, this proceeding);
* a the projection of the central part of an oblate bulge inside the primary bar onto the sky plane;
* a complex distribution of dust and star formation regions in the circumnuclear region.

It seems to us that additional observations of the stellar and gas kinematics in these strange objects might prove if the secondary bars are new dynamically isolated structures in barred galaxies. Nevertheless, such observations are very rare (Wozniak, 1999; Emsellem et al., 2001). Emsellem et al. (2001) have presented long-slit stellar kinematics of 4 double-barred objects. However two-dimensional distributions of velocities and velocity dispersions can give us a more detailed information because the objects observed are non-axisymmetrical.

## 3.  Observations

In 2000 we observed a sample of the northern double-barred galaxies at the 6m telescope of the SAO RAS by using various 2D spectroscopy methods. The circumnuclear regions were investigated with the integral-field spectrograph MPFS (Afanasiev et al., this Proceeding). The velocity fields of the ionized gas in the \[OIII\]$`\lambda 5007`$ and $`H_\beta `$ emission lines, the stellar line-of-sight velocity and velocity dispersion (hereafter $`\sigma `$) fields were constructed over $`16{}_{}{}^{\prime \prime }\times 15^{\prime \prime }`$ field of view. For some galaxies we obtained large-scale ionized-gas velocity fields in the \[NII\]$`\lambda 6583`$ and/or $`H_\alpha `$ emission lines by using the scanning Fabry-Perot Interferometer (IFP). The analysis of the morphology of the galaxies was done by means of JHK’ images obtained with CAMILA (Cruz-Gonzalez et al., 1994), at the 2.1m telescope (OAN, Mexico). Optical images from the 6m telescope and high-resolution images from the HST Archive were also used for the same purpose.

## 4. Results

Fig. 1 shows our data for NGC 2950. It is an SB0 galaxy, at the images of which the inner bar-like structure is clearly seen that is turned by approximately $`60^{}`$ with respect to the primary bar. This “secondary bar” was firstly described by Wozniak et al. (1995) and Friedli et al. (1996). Surprisingly, the position angle (hereafter PA) of the kinematical axis defined from the circumnuclear stellar velocity field, $`PA_{KIN}=125\pm 3^{}`$, is not significantly different from the orientation of the outer isophotes of the disk, $`PA_{DISK}=120\pm 5^{}`$. This fact is unexpected111Our present value of the $`PA_{DISK}`$ obtained from the deep R-band image disagrees with the $`PA_{DISK}=110^{}`$ from Friedli et al.(1996). The latter value was also used by Moiseev (2001) that has resulted in a erroneous conclusion on the significant difference between $`PA_{KIN}`$ and $`PA_{DISK}`$ in NGC 2950. because the secondary bar-like structure is sharply distinguished at the galactic images and must affect the stellar orbital motions if it is triaxial.

The central ellipsoidal distribution ($`r5^{\prime \prime }`$) is seen at the velocity dispersion map (see Fig. 1), with the major axis of this distribution being coincident with the outer (not with the inner!) bar orientation. This situation is typical for single-barred galaxies, where distributions of line-of-sight $`\sigma _{}`$ in the bar would be symmetric to the bar axis but not to the major axis of the disk (Miller & Smith, 1979; Vauterin & Dejonghe, 1997). In other objects observed we have found more complex velocity dispersion distributions than regular elliptical peaks (see Moiseev, 2001). These are peaks shifting from the dynamical centers or central velocity dispersion minima. The latter may be connected with a presence of a cold dynamical component (Emsellem, 2001) or with the evolution of the stellar bars (see simulations by Khoperskov, Moiseev & Chulanova, 2002).

We have considered some statistical properties of the our sample. If a secondary bar is truly decoupled then the galaxy has three main directions: $`PA_{DISK}`$ of a disk major axis, and also $`PA_{PRIMARY}`$ and $`PA_{SECONDARY}`$ that are traced by both bars in the sky plane. Fig. 2 shows the $`PA`$ of the observed kinematical features versus the $`PA`$ of these main directions. Here the kinematical features are $`PA_\sigma _{}`$ – the orientation of the elliptical $`\sigma _{}`$-peak (if it is present) and $`PA_{KIN}`$ – orientation of the central stellar kinematics axis. One can clearly see that the stellar velocity dispersion distributions correlate only with the primary bars. The secondary bars have not any influence on the stellar kinematics.

The bottom plot in Fig. 2 reveals two important facts. Firstly, the deviations of $`PA_{KIN}`$ and $`PAs`$ of the both bars from the galactic major axis do not correlate. Secondly, one can see that the circumnuclear kinematical axes, deduced from the velocity fields, correlate with the global disk orientations. Therefore both gas and stars orbital motions are mostly axisymmetrical in the nuclei of the observed galaxies. The residual deviations $`PA(KIN)`$ from $`PA_{DISK}`$ could be connected with non-circular gas motions. A detailed consideration of the ionized gas kinematics in individual objects reveals various types of such motions. For example, the large-scale (primary) bar in NGC 470 provokes non-circular gas flows in the IFP velocity field. At the secondary bar distance ($`r<8^{\prime \prime }`$) the isovelocity turn may be explained by a radial gas inflow in the bar or by a polar disk appearance within the primary bar. The comparison with the isophotal analysis results allows us to choose the disk instead of the bar.

In the table 1 we summarize all the circumnuclear properties that were found in our sample. Although various kinematical features are present, only in the case of NGC 3368 the gas kinematics may be explained by a dynamically decoupled secondary bar.

## 5. Do mini-spirals nest in the primary bars?

In order to study the residual brightness morphology of the galaxies, we have subtracted elliptical isophotal models from the optical and NIR images. Interestingly, about half of all galaxies have mini-spirals or pseudo-spirals around their nuclei. Three examples of nuclear spirals are shown in Fig. 3. The multi-armed flocculent spiral in NGC 4736 is studied also by Elmegreen, Elmegreen & Eberwein (2002). The two-armed spiral in NGC 5566 resembles model calculations of the gas behaviour in a stellar bar by Englmaier & Shlosman (2000). Pérez et al. (2000) found a flocculent mini-spiral structure on the HST sharp-divided images inside $`r<6^{\prime \prime }`$ in NGC 6951. We suggest that this spiral causes a unusual residual (non-circular) velocity picture that was observed at our Fabry-Perot velocity field of NGC 6951.

In the recent years, the interest for the nuclear spirals has considerably grown, because these features have been found in the high-resolution HST images of some galaxies(Van Dokkum & Franx (1995), Elmegreen et al.(2002) and references there). The origin of these structures is yet unknown, different models are proposed in the literature. In our sample a clear connection is seen between the large-scale bars and mini-spiral structure within them.

## 6.  Conclusions

A sample of double-barred galaxies has been observed at the 6m telescope by the 2D-spectroscopy methods. In order to know about the morphology of the sample galaxies, we also carried out NIR imaging observation with CAMILA, at the 2.1m telescope (OAN, México) The velocity fields of the ionized gas and stars and the line-of-sight stellar velocity dispersion maps are constructed and analyzed. The following features are found:

* The line-of-sight velocity dispersion distributions of stars are connected with the primary bars only. Also the circumnuclear stellar kinematical axes do not show significant misalignment with the $`PA`$’s of the main discs. Therefore the secondary bars are not dynamically decoupled stellar systems.
* Remarkable non-circular motions of the ionized gas are found in $`60\%`$ of all objects. These features may be associated with radial gas inflow in the primary bar, or with non-circular motions between ILRs of the primary bar, or with polar (inclined) circumnuclear discs, or with gaseous motions in the circumnuclear mini-spirals.
* The analysis of the optical and NIR images shows that mini-spiral structures exist at the secondary bar radial scales in 6 out of 13 objects.

Thus we suggest that a dynamically independent secondary bar is a much more rare phenomenon than it follows from the isophotal analysis of images itself.

### Acknowledgments.

We would like to thank V.L. Afanasiev and O.K. Sil’chenko for their interest and numerous discussions. This work was partially supported by the RFBR grant (Project No. 01-02-17597) and by the Federal program ”Astronomy” (Project 1.2.3.1). This work is based on the observations from the 6m telescope of the Special Astrophysical Observatory of the Russian Academy of Sciences operated under the financial support of the Science Department of Russia (registration number 01-43). The authors would like to thank all the staff of the Observatorio Astronómico Nacional, in San Pedro Martir, México for their assistance during the NIR observations

## References

Afanasiev, V.L. & Sil’chenko, O.K. 2000, AJ, 119, 126

Buta, R. & Crocker, D.A. 1993, AJ, 105, 1344

Cruz-Gonzalez, I., Carrasco, L., Ruiz, E., Salas, L., Skrutskie, M., et al. 1994, Rev. Mex. Astron. Astrofis., 29, 197

Elmegreen, D.M., Elmegreen B.G., & Eberwein, K.S. 2002, astro-ph/0112527

Englmaier, P., & Shlosman, I., 2000, ApJ, 528, 677

Emsellem, E. 2001, astro-ph/0107299

Emsellem, E., Greusard, D., Combes, F., et al. 2001, A&A, 368, 52

Erwin, P. & Sparke, L.S. 1998, (astro-ph/9811345)

de Vaucouleurs, G. 1975, ApJS, 29, 193

Friedli, D. & Martinet, L. 1993, A&A, 277, 27

Friedli, D., Wozniak, H., Rieke, M., Martinet, L., Bratschi, P. 1996, A&A, 118, 461

Heller, C., Shlosman, I. & Englmaier, P. 2001, ApJ, 553, 661

Khoperskov, A.V., Moiseev, A.V., & Chulanova, E.A. 2002, Bull. SAO, 52, 135

Laine, S., Shlosman, I., Knapen, J.H., Peletier R.F. 2001, to appear in ApJ, (astro-ph/0108029)

Maciejewski, W. & Sparke, L.S. 2000, MNRAS, 313, 745

Maciejewski, W., Teuben, P.J., Sparke, L.S., Stone, J.M. 2002, MNRAS, 329, 502

Miller, R.H. & Smith, B.F. 1979, ApJ, 227, 785

Moiseev, A.V. 2001, Bull. SAO, 51, 140, (astro-ph/0111220)

Pérez, E., Márquez, I., Durret, F. et al. 2001, A&A, 353, 893

Pfenninger, D. & Norman, C.A. 1990, Apj, 363, 391

Shlosman, I., Frank, J. & Begeman, M.C. 1989, Nature, 338, 45

Shlosman, I., & Heller, C. 2001, astro-ph/0109536

van Dokkum, P.G., & Franx, M. 1995, AJ, 110, 2027

Vauterin, P., & Dejonghe, H. 1997, MNRAS, 286, 812

Wozniak, H., Friedli, D., Martinet, L., et al. 1995, A&AS, 111, 115