Revision as of 19:49, 7 November 2014

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 continuity equation, the Euler equation, and the Poisson equation become, respectively,

<math>~\frac{1}{\rho} \frac{d\rho}{dt} </math>	<math>~=</math>	<math>~-~ a^{-2} \nabla_\mathfrak{x}^2 \psi \, ;</math>
<math>~\frac{d\psi}{dt} </math>	<math>~=</math>	<math>~\frac{1}{2} a^{-2} ( \nabla_\mathfrak{x} \psi )^2 - H - \Phi \, ;</math>
<math>~a^{-2}\nabla_\mathfrak{x}^2 \Phi </math>	<math>~=</math>	<math>~4\pi G \rho \, .</math>

Reconciling with Goldreich & Weber

The set of three principal governing equations, as just derived, are intended to match equations (7) - (9) of Goldreich & Weber (1980). The following is a framed image of equations (7) - (9) as they appear in the Goldreich & Weber publication:

Principal Governing Equations from Goldreich & Weber (1980)

For discussion purposes, next we will retype this set of equations, altering only the variable names and notation to correspond with ours. Assuming that we have interpreted their typeset expressions correctly, the governing equations, as derived by Goldreich & Weber, are,

<math>~\frac{1}{\rho} \frac{\partial\rho}{\partial t} ~+ a^{-1}(a^{-1}\nabla_\mathfrak{x}\psi - \dot{a}\mathfrak{x})\cdot \frac{\nabla_\mathfrak{x}\rho}{\rho}+~ a^{-2} \nabla_\mathfrak{x}^2 \psi </math>	<math>~=</math>	<math>0 \, ;</math>
<math>~\frac{\partial\psi}{\partial t} - \frac{\dot{a} \mathfrak{x}}{a} \cdot \nabla_\mathfrak{x} \psi~+ \frac{1}{2} a^{-2}( \nabla_\mathfrak{x} \psi )^2 + H + \Phi </math>	<math>~=</math>	<math>0 \, ;</math>
<math>~a^{-2}\nabla_\mathfrak{x}^2 \Phi - 4\pi G \rho</math>	<math>~=</math>	<math>~0 \, .</math>

Notice that our expression for the Poisson equation matches the expression presented by Goldreich & Weber, but it isn't immediately obvious whether or not the other two pairs of equations match. Let's rearrange the terms in Goldreich & Weber's continuity equation and in their Euler equation to emphasize overlap with ours:

<math>~\frac{1}{\rho} \biggl[ \frac{\partial\rho}{\partial t} + (a^{-1}\nabla_\mathfrak{x}\psi - \dot{a}\mathfrak{x})\cdot a^{-1}\nabla_\mathfrak{x}\rho \biggr] </math>	<math>~=</math>	<math>~-~ a^{-2} \nabla_\mathfrak{x}^2 \psi \, ;</math>
<math>~\frac{\partial\psi}{\partial t} +(a^{-1}\nabla_\mathfrak{x}\psi - \dot{a}\mathfrak{x})\cdot a^{-1}\nabla_\mathfrak{x}\psi</math>	<math>~=</math>	<math>\frac{1}{2} a^{-2}( \nabla_\mathfrak{x} \psi )^2 - H - \Phi </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>.

@@ Line 262: / Line 262: @@
    </td>
    <td align="left">
-<math>~\frac{1}{2a^2} ( \nabla_\mathfrak{x} \psi )^2 - H - \Phi  \, ;</math>
+<math>~\frac{1}{2} a^{-2} ( \nabla_\mathfrak{x} \psi )^2 - H - \Phi  \, ;</math>
    </td>
 </tr>

Difference between revisions of "User:Tohline/Apps/GoldreichWeber80"