Homologously Collapsing Stellar Cores

Review of Goldreich and Weber (1980)

This is principally a review of the dynamical model that Peter Goldreich & Stephen Weber (1980, ApJ, 238, 991) developed to describe the near-homologous collapse of stellar cores. As we began to study the Goldreich & Weber paper, it wasn't immediately obvious how the set of differential governing equations should be modified in order to accommodate a radially contracting (accelerating) coordinate system. I did not understand the transformed set of equations presented by Goldreich & Weber as equations (7) and (8), for example. At first, I turned to Poludnenko & Khokhlov (2007, Journal of Computational Physics, 220, 678) — hereafter, PK07 — for guidance. PK07 develop a very general set of governing equations that allows for coordinate rotation as well as expansion or contraction. Ultimately, the most helpful additional reference proved to be §19.11 (pp. 187 - 190) of Kippenhahn & Weigert [ KW94 ].

Governing Equations

Goldreich & Weber begin with the identical set of principal governing equations that serves as the foundation for all of the discussions throughout this H_Book. In particular, as is documented by their equation (1), their adopted equation of state is adiabatic/polytropic,

<math>~P = \kappa \rho^\gamma \, ,</math>

— where both <math>~\kappa</math> and <math>~\gamma</math> are constants — and therefore satisfies what we have referred to as the

Adiabatic Form of the
First Law of Thermodynamics
(Specific Entropy Conservation)

<math>~\frac{d\epsilon}{dt} + P \frac{d}{dt} \biggl(\frac{1}{\rho}\biggr) = 0</math> .

their equation (3) is what we have referred to as the

Euler Equation
in terms of the Vorticity,

<math>~\frac{\partial\vec{v}}{\partial t} + \vec\zeta \times \vec{v}= - \frac{1}{\rho} \nabla P - \nabla \biggl[\Phi + \frac{1}{2}v^2 \biggr] </math>

where, <math>~\vec\zeta \equiv \nabla\times \vec{v}</math> is the fluid vorticity; their equation (4) is the

Poisson Equation

<math>\nabla^2 \Phi = 4\pi G \rho</math>

and their equation (2) is what we have referred to as the

Eulerian Representation
or
Conservative Form
of the Continuity Equation,

<math>~\frac{\partial\rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0</math>

although, for the derivation, below, we prefer to start with what we have referred to as the

Standard Lagrangian Representation
of the Continuity Equation,

<math>\frac{d\rho}{dt} + \rho \nabla \cdot \vec{v} = 0</math>

Tweaking the set of principal governing equations, as we have written them, to even more precisely match equations (1) - (4) in Goldreich & Weber (1980), we should replace the state variable <math>~P</math> (pressure) with <math>~H</math> (enthalpy), keeping in mind that, <math>~\gamma = 1 + 1/n</math>, and, as presented in our introductory discussion of barotropic supplemental relations,

<math>~H = \biggl( \frac{\gamma}{\gamma-1} \biggr) \kappa \rho^{\gamma-1} \, ,</math>

and,

<math>~\nabla H = \frac{\nabla P}{\rho} \, .</math>

Imposed Constraints

Goldreich & Weber (1980) specifically choose to examine the spherically symmetric collapse of a <math>~\gamma = 4/3</math> fluid. With this choice of adiabatic index, the equation of state becomes,

<math>~H = 4 \kappa \rho^{1/3} \, .</math>

And because a strictly radial flow-field exhibits no vorticity (i.e., <math>\vec\zeta = 0</math>), the Euler equation can be rewritten as,

<math>~\frac{\partial v_r}{\partial t} </math>

<math>~-~ \nabla_r \biggl[ H + \Phi + \frac{1}{2}v_r^2 \biggr] \, .</math>

Goldreich & Weber also realize that, because the flow is vorticity free, the velocity can be obtained from a stream function, <math>~\psi</math>, via the relation,

<math>~\vec{v} = \nabla\psi ~~~~~\Rightarrow~~~~~v_r = \nabla_r\psi </math> and <math>~\nabla\cdot \vec{v} = \nabla_r^2 \psi \, .</math>

Hence, the continuity equation becomes,

<math>~-~ \nabla_r^2 \psi \, ,</math>

and the Euler equation becomes,

<math>~\frac{\partial }{\partial t} \biggl[ \nabla_r \psi \biggr]</math>

<math>~-~ \nabla_r \biggl[ H + \Phi + \frac{1}{2}(\nabla_r \psi)^2 \biggr] \, .</math>

Since we are, up to this point in the discussion, still referencing the inertial-frame radial coordinate, the <math>~\nabla_r</math> operator can be moved outside of the partial time-derivative on the lefthand side of this equation to give,

<math>~\nabla_r \biggl[ \frac{\partial \psi}{\partial t} + H + \Phi + \frac{1}{2}(\nabla_r \psi)^2 \biggr]</math>

This means that the terms inside the square brackets must sum to a constant that is independent of spatial position. Following the lead of Goldreich & Weber, this "integration constant" will be incorporated into the potential, in which case we have,

<math>~\frac{\partial \psi}{\partial t} </math>

<math>~-~ \biggl[ H + \Phi + \frac{1}{2} ( \nabla_r \psi )^2 \biggr] \, ,</math>

which matches equation (5) of Goldreich & Weber (1980).

Now, because it is more readily integrable, we ultimately would like to work with a differential equation that contains the total, rather than partial, time derivative of <math>~\psi</math>. So we will take this opportunity to shift from an Eulerian representation of the Euler equation to a Lagrangian representation, invoking the same (familiar to fluid dynamicists) operator transformation as we have used in our general discussion of the Euler equation, namely,

<math>~\frac{\partial\psi}{\partial t} ~~ \rightarrow ~~ \frac{d\psi}{dt} - \vec{v}\cdot \nabla\psi \, .</math>

In the context of Goldreich & Weber's model, we are dealing with a one-dimension (spherically symmetric), radial flow, so,

<math>\vec{v}\cdot \nabla\psi = v_r \nabla_r \psi \, .</math>

But, given that we have adopted a stream-function representation of the flow in which <math>~v_r = \nabla_r\psi</math>, we appreciate that this term can either be written as <math>~v_r^2</math> or <math>~(\nabla_r\psi)^2</math>. We choose the latter representation, so the Euler equation becomes,

<math>~\frac{d\psi}{dt} - (\nabla_r\psi)^2</math>

<math>~-~ \biggl[ H + \Phi + \frac{1}{2} ( \nabla_r \psi )^2 \biggr] \, ,</math>

or, combining like terms on the left and right,

<math>~\frac{1}{2} ( \nabla_r \psi )^2 - H - \Phi \, .</math>

Dimensionless and Time-Dependent Normalization

Length

In their investigation, Goldreich & Weber (1980) chose the same length scale for normalization that is used in deriving the Lane-Emden equation, which governs the hydrostatic structure of a polytrope of index <math>~n</math>, that is,

<math> a_\mathrm{n} \equiv \biggl[\frac{1}{4\pi G}~ \biggl( \frac{H_c}{\rho_c} \biggr)\biggr]^{1/2} \, , </math>

where the subscript, "c", denotes central values. In this case <math>~(n = 3)</math>, substitution of the equation of state expression for <math>~H_c</math> leads to,

<math> a = \rho_c^{-1/3} \biggl(\frac{\kappa}{\pi G}\biggr)^{1/2} \, . </math>

Most significantly, Goldreich & Weber (see their equation 6) allow the normalizing scale length to vary with time in order for the governing equations to accommodate a self-similar dynamical solution. In doing this, they effectively adopted an accelerating coordinate system with a time-dependent dimensionless radial coordinate,

<math>~\vec\mathfrak{x} \equiv \frac{1}{a(t)} \vec{r} \, .</math>

This, in turn, will mean that either the central density varies with time, or the specific entropy of all fluid elements (captured by the value of <math>~\kappa</math>) varies with time, or both. In practice, Goldreich & Weber assume that <math>~\kappa</math> is held fixed, so the time-variation in the scale length, <math>~a</math>, reflects a time-varying central density; specifically,

<math> \rho_c = \biggl(\frac{\kappa}{\pi G}\biggr)^{3/2} [a(t)]^{-3} \, . </math>

Given the newly adopted dimensionless radial coordinate, the following replacements for the spatial operators should be made, as appropriate, throughout the set of governing equations:

<math>~\nabla_r ~\rightarrow~ a^{-1} \nabla_\mathfrak{x}</math> and <math>~\nabla_r^2 ~\rightarrow~ a^{-2} \nabla_\mathfrak{x}^2 \, .</math>

Specifically, the Poisson equation becomes,

<math>\nabla_\mathfrak{x}^2 \Phi = 4\pi G a^2 \rho \, ;</math>

the Euler equation becomes,

<math>~\frac{1}{2a^2} ( \nabla_\mathfrak{x} \psi )^2 - H - \Phi \, ;</math>

and the continuity equation becomes,

<math>~-~ a^{-2} \nabla_\mathfrak{x}^2 \psi \, .</math>

Mass-Density and Speed

Next, Goldreich & Weber (1980) (see their equation 10) choose to normalize the density by the central density, specifically defining a dimensionless function,

<math>f \equiv \biggl( \frac{\rho}{\rho_c} \biggr)^{1/3} \, .</math>

Keeping in mind that <math>~n = 3</math>, this is also in line with the formulation and evaluation of the Lane-Emden equation, where the primary dependent structural variable is the dimensionless polytropic enthalpy,

<math>\Theta_H \equiv \biggl( \frac{\rho}{\rho_c} \biggr)^{1/n} \, .</math>

Also, Goldreich & Weber (1980) (see their equation 11) normalize the gravitational potential to the square of the central sound speed,

<math>c_s^2 = \frac{\gamma P_c}{\rho_c} = \frac{4}{3} \kappa \rho_c^{1/3} = \frac{4}{3}\biggl(\frac{\kappa^3}{\pi G}\biggr)^{1/2} [a(t)]^{-1} \, .</math>

Specifically, their dimensionless gravitational potential is,

<math>~\sigma</math>

<math>~\equiv</math>

<math>~\frac{\Phi}{c_s^2} = \biggl[ \frac{3}{4} \biggl( \frac{\pi G}{\kappa^3} \biggr)^{1/2} a(t) \biggr] \Phi \, ,</math>

and the similarly normalized enthalpy may be written as,

<math>~\frac{H}{c_s^2} </math>	<math>~=</math>	<math>~\biggl[ \frac{3}{4} \biggl( \frac{\pi G}{\kappa^3} \biggr)^{1/2} a(t) \biggr] 4\kappa \rho^{1/3} </math>
	<math>~=</math>	<math>~3 \biggl( \frac{\rho}{\rho_c} \biggr)^{1/3} </math>
	<math>~=</math>	<math>~3f \, .</math>

With these additional scalings, the continuity equation becomes,

<math>~-~ a^{-2} \nabla_\mathfrak{x}^2 \psi \, ,</math>

the Euler equation becomes,

<math>~ \biggl[ \frac{3}{4} \biggl( \frac{\pi G}{\kappa^3} \biggr)^{1/2} a(t) \biggr] \biggl[ \frac{d\psi}{dt} - \frac{1}{2a^2} ( \nabla_\mathfrak{x} \psi )^2 \biggr] </math>

<math>~ - 3 f - \sigma \, ;</math>

and the Poisson equation becomes,

<math>\nabla_\mathfrak{x}^2 \sigma = 3f^3 \, .</math>

Homologous Solution

Goldreich & Weber (1980) discovered that the governing equations admit to an homologous, self-similar solution if they adopted a stream function of the form,

<math>~\frac{1}{2}a \dot{a} \mathfrak{x}^2 \, ,</math>

which, when acted upon by the various relevant operators, gives,

<math>~\nabla_\mathfrak{x}\psi</math>	<math>~=</math>	<math>~a \dot{a} \mathfrak{x} \, ,</math>
<math>~\nabla^2_\mathfrak{x}\psi</math>	<math>~=</math>	<math>~ \biggl( \frac{1}{2}a \dot{a} \biggr) \frac{1}{\mathfrak{x}^2} \frac{d}{d\mathfrak{x}} \biggl[\mathfrak{x}^2 \frac{d}{d\mathfrak{x}} \mathfrak{x}^2 \biggr] = 3 a \dot{a} \, , </math>
<math>~\frac{d\psi}{dt}</math>	<math>~=</math>	<math>~\mathfrak{x}^2 \biggl[ \frac{1}{2}\dot{a}^2 + \frac{1}{2}a\ddot{a} \biggr] \, .</math>

Hence, the radial velocity profile is,

<math>~v_r = a^{-1}\nabla_\mathfrak{x} \psi</math>

<math>~\dot{a}\mathfrak{x} \, ; </math>

the continuity equation gives,

<math>~3 ~\frac{d\ln f}{dt} </math>	<math>~=</math>	<math>-~ \frac{3\dot{a}}{a} </math>
<math>\Rightarrow~~~~\frac{d\ln f}{dt} + \frac{d\ln a}{dt} </math>	<math>~=</math>	<math>~0 \, ,</math>

which means that the product, <math>~af</math>, is independent of time; and the Euler equation becomes,

<math>~ - 3 f - \sigma </math>	<math>~=</math>	<math>~ \biggl[ \frac{3}{4} \biggl( \frac{\pi G}{\kappa^3} \biggr)^{1/2} a(t) \biggr] \biggl\{ \mathfrak{x}^2 \biggl[ \frac{1}{2}\dot{a}^2 + \frac{1}{2}a\ddot{a} \biggr] - \frac{1}{2} ( \dot{a} \mathfrak{x} )^2 \biggr\} </math>
	<math>~=</math>	<math>~ \frac{3}{8} \biggl( \frac{\pi G}{\kappa^3} \biggr)^{1/2} (a \mathfrak{x})^2 \ddot{a} </math>
<math>~\Rightarrow~~~~ \frac{(f + \sigma/3)}{\mathfrak{x}^2} </math>	<math>~=</math>	<math>~ - \frac{1}{8} \biggl( \frac{\pi G}{\kappa^3} \biggr)^{1/2} a^2 \ddot{a} \, , </math>

which matches equation (12) of Goldreich & Weber (1980).

Because everything on the lefthand side of Goldreich & Weber's scaled Euler equation depends only on the dimensionless spatial coordinate, <math>~\mathfrak{x}</math>, while everything on the righthand side depends only on time — via the parameter, <math>~a(t)</math> — both expressions must equal the same (dimensionless) constant. Goldreich & Weber (1980) (see their equation 12) call this constant, <math>~\lambda/6</math>. From the terms on the lefthand side, they conclude (see their equation 13) that the dimensionless gravitational potential is,

<math>~\sigma</math>

<math>~\frac{1}{2} \lambda ~\mathfrak{x}^2 - 3f \, .</math>

From the terms on the righthand side they conclude, furthermore, that the nonlinear differential equation governing the time-dependent variation of the scale length, <math>~a</math>, is,

<math>~-~\frac{4\lambda}{3} \biggl( \frac{\kappa^3}{\pi G} \biggr)^{1/2} \, .</math>

As Goldreich & Weber (1980) point out, this nonlinear differential equation can be integrated twice to produce an algebraic relationship between <math>~a</math> and time, <math>~t</math>. The required mathematical steps are identical to the steps used to analytically solve the classic, spherically symmetric free-fall collapse problem. First, rewrite the equation as,

where,

<math> ~B \equiv \frac{8\lambda}{3} \biggl( \frac{\kappa^3}{\pi G} \biggr)^{1/2} \, , </math>

has the same dimensions as the product, <math>~GM</math> (see the free-fall collapse problem), that is, the dimensions of "length-cubed per unit time-squared." Then, multiply both sides by <math>~2\dot{a} = 2(da/dt)</math> to obtain,

<math>~ 2\dot{a} \frac{d\dot{a}}{dt} </math>	<math>~=</math>	<math>~-B \biggl( a^{-2} \frac{da}{dt} \biggr) </math>
<math>~\Rightarrow~~~~ \frac{d\dot{a}^2}{dt} </math>	<math>~=</math>	<math>~B \frac{d}{dt} \biggl( \frac{1}{a} \biggr) \, ,</math>

which integrates once to give,

or,

<math> ~dt = \biggl( \frac{B}{a} + C \biggr)^{-1/2} da \, . </math>

For the case, <math>~C = 0</math>, this differential equation can be integrated straightforwardly to give (see Goldreich & Weber's equation 15),

For the cases when <math>~C \ne 0</math>, Wolfram Mathematica's online integrator can be called upon to integrate this equation and provide the following closed-form solution,

<math> \frac{a}{C} \biggl( \frac{B}{a} + C \biggr)^{1/2} - \frac{B}{2C^{3/2}} \ln \biggl[2aC^{1/2} \biggl( \frac{B}{a} + C \biggr)^{1/2} + B + 2aC \biggr] \, . </math>

Related Discussions

Material that appears after this point in our presentation is under development and therefore
may contain incorrect mathematical equations and/or physical misinterpretations.
| Go Home |

As Goldreich & Weber (1980) point out, because all terms in this equation are inside the gradient operator, the sum of the terms inside the square brackets must equal a constant — that is, the sum must be independent of spatial position throughout the spherically symmetric configuration. If, following Goldreich & Weber's lead, we simply fold this integration constant into the potential, the Euler equation becomes (see their equation 8),

<math>~\frac{\partial \psi}{\partial t} - \biggl( \frac{\dot{a}}{a} \biggr)\psi + H + \Phi + \frac{1}{2}\biggl(\frac{1}{a} \nabla_x\psi \biggr)^2 </math>

<math>~\frac{\partial \rho}{\partial t} + \rho \nabla_r \cdot \vec{v} + \vec{v}\cdot \nabla_r \rho</math>	<math>~=</math>	<math>~0</math>
<math>~\Rightarrow ~~~~ \frac{1}{\rho} \frac{\partial \rho}{\partial t} + \nabla_r \cdot \vec{v} + \vec{v}\cdot \frac{\nabla_r \rho}{\rho}</math>	<math>~=</math>	<math>~0</math>
<math>~\Rightarrow ~~~~ \frac{1}{\rho} \frac{\partial \rho}{\partial t} + a^{-1} \nabla_x \cdot \biggl[ a^{-1} \nabla_x \psi \biggr] + a^{-1} \nabla_x \psi \cdot \frac{a^{-1}\nabla_x \rho}{\rho}</math>	<math>~=</math>	<math>~0</math>
<math>~\frac{1}{\rho} \frac{\partial \rho}{\partial t} + a^{-1}(a^{-1} \nabla_x\psi - \dot{a} \vec{x}) \cdot \frac{\nabla_x\rho}{\rho} + a^{-2} \nabla_x^2\psi </math>	<math>~=</math>	<math>~0</math>

Goldreich & Weber's (1980) Governing Equations After Initial Length Scaling (yet to be demonstrated)

<math>~\frac{1}{\rho} \frac{\partial \rho}{\partial t} + a^{-1}(a^{-1} \nabla_x\psi - \dot{a} \vec{x}) \cdot \frac{\nabla_x\rho}{\rho} + a^{-2} \nabla_x^2\psi </math>	<math>~=</math>	<math>~0</math>
<math>~\frac{\partial \psi}{\partial t} - \frac{\dot{a}}{a} \vec{x}\cdot \nabla_x\psi + \frac{1}{2} a^{-2} \| \nabla_x\psi\|^2 + H + \Phi</math>	<math>~=</math>	<math>~0</math>
<math>~ a^{-2} \nabla_x^2\Phi - 4\pi G \rho </math>	<math>~=</math>	<math>~0</math>
where, <math>~\vec{x} \equiv \frac{\vec{r}}{a} \, ,</math> and it is understood that derivatives in the <math>~\nabla_x</math> and <math>~\nabla_x^2</math> operators are taken with respect to the dimensionless radial coordinate, <math>~x</math>.

User:Tohline/Apps/GoldreichWeber80

Contents