Blandford-McKee 1976 Reading Note

This post is to take some reading notes of Blandford-McKee 1976 self-similar solution, and derive the main equations in detail.

Relativistic Shock Jump Condition

Define the dimensionless 4-velocity \(u = ds/cd\tau\), therefore spacial component of u is \(\gamma \beta\), and time component is \(1\).

We further define the relative dimensionless velocity \(\beta_{ij}\), e.g. \(\beta_{1s}\) the relative velocity between shock and upstream unshocked plasma. From SR, we easily get the following relations

\[\begin{equation} \begin{aligned} &\beta_{1s} = \frac{\beta_{2s}+\beta_{21}}{1+\beta_{2s}\beta_{21}}\\ &\beta_{2s} = \frac{\beta_{1s}-\beta_{21}}{1-\beta_{1s}\beta_{21}}\\ &\beta_{21} = \frac{\beta_{1s}-\beta_{2s}}{1-\beta_{1s}\beta_{2s}}\\ &\gamma_{1s} = \gamma_{2s}\gamma_{21}(1+\beta_{2s}\beta_{21})\\ &\gamma_{2s} = \gamma_{1s}\gamma_{21}(1-\beta_{1s}\beta_{21})\\ &\gamma_{21} = \gamma_{1s}\gamma_{2s}(1-\beta_{1s}\beta_{2s})\\ &u_{1s} = \gamma_{2s}\gamma_{21}(\beta_{21}+\beta_{2s})\\ &u_{2s} = \gamma_{1s}\gamma_{21}(\beta_{1s}-\beta_{21})\\ &u_{21} = \gamma_{1s}\gamma_{2s}(\beta_{1s}-\beta_{2s}) \end{aligned} \end{equation}\]

This is useful to obtain jump conditions in upstream rest frame. The shock needs to satisfy mass, momentum and energy conservation. In shock frame, the conservation condition could be derived from \(\partial_x T^{1\nu} = 0\). Obviously we must first have the particle number conservation

\[\begin{equation} \langle nu\rangle = \langle n\gamma \beta \rangle = 0 \end{equation}\]

Then from \(\partial_x T^{00} = 0\) follows the energy conservation condition:

\[\begin{equation} \langle \gamma^2 h \beta\rangle = \langle \gamma h u\rangle = 0 .\end{equation}\]

From \(\partial_x T^{01} = 0\), we have the momentum conservation:

\[\begin{equation} \langle \gamma^2 \beta^2 h +p\rangle = \langle hu^2+p\rangle = 0 \end{equation}\]

here \(h\) is the enthalpy in unit volume in fluid rest frame. We need to note that the \(n, p, e , h\) are all quantities in the fluid rest frame. The relativistic transformation is contained in the dimensionless 4-speed.

Rankine-Hugoniot Condition for Cold upstream

Cold upstream means that \(e_1 = \rho_1 c^2, p_1= 0\). Thus from the equation of state \(p = (\hat{\gamma} - 1)(e - \rho c^2)\), there is \(\frac{p}{n} = (\hat{\gamma}-1)(\frac{e}{n}-mc^2)\). If the internal energy at upstream \(e_1-\rho c^2=0\), since the \(e/n\) is Lorentz invariant, we have

\[\begin{equation} \frac{e_2}{n_2} = \gamma_{21} m c^2 \end{equation}\]

Therefore,

\(\frac{p_2}{n_2} = (\hat{\gamma}-1)(\gamma-1) mc^2\).

From the energy conservation \(\gamma_{1s} h_1/n_1=\gamma_{2s} h_2/n_2\), we have \(\begin{equation} \frac{\gamma_{1s}}{\gamma_{2s}} = \hat{\gamma}(\gamma_{21} -1)+1 \end{equation}\)

Keep this in mind and substitute it into

\[\gamma_{2s} = \gamma_{1s}\gamma_{21}(1-\beta_{1s}\beta_{21}),\]

eliminate \(\gamma_{2s}\) and after a bit complicated formulation we get

\[\begin{equation} \gamma_{1s}^2 = \frac{(\gamma_{21} + 1)[\hat{\gamma}(\gamma_{21}-1)+1]^2}{\hat{\gamma}(2-\hat{\gamma} )(\gamma_{21}-1)+2}. \end{equation}\]

From Lorentz factor, we can directly write the dimensionless 4-velocity \(u^2 = \gamma^2 \beta^2 = \gamma^2-1\)

\[\begin{equation} \begin{aligned} u_{1s}^2 = &\frac{(\gamma_{21} - 1)(\gamma_{21}\hat{\gamma}+1)^2}{(\gamma_{21}-1)+2}\\ u_{2s}^2 =&\frac{(\gamma_{21} - 1)(\hat{\gamma}-1)^2}{(\gamma_{21}-1)+2} \end{aligned} \end{equation}\]

Therefore, from \(u_{2s} n_{2} = u_{1s}n_1\), the compression ratio is

\[\begin{equation} \frac{n_2}{n_1} = \frac{\hat{\gamma}\gamma_{21}+1}{\hat{\gamma}-1} \end{equation}\]

These are the Rankine-Hugoniot conditions of relativistic shock, see eq. (3-5) in BM76.

Energy-Momentum Tensor

The energy-momentum tensor could be expressed in the following form

\[\begin{equation} T^{\mu\nu} = (e+p) u^\mu u^\nu + p \eta ^{\mu\nu}, \end{equation}\]

where \(e\) is the energy density (internal energy + rest mass), \(p\) is the pressure, \(u = \gamma( 1,\beta_x,\beta_y,\beta_z)\) the dimensionless 4-velocity, \(\eta = (-1,1,1,1)\) the Minkowski metric.

Specifically, if a relativistic fluid is moving with \(\beta\) in \(x\) direction, the complete matrix form of EM Tensor is

\[\begin{equation} T^{\mu\nu} = \left(\begin{matrix} \gamma^2h-p & \gamma^2\beta h & 0&0\\ \gamma^2\beta h & \gamma^2 \beta^2 h +p & 0 & 0 \\ 0&0&p&0\\ 0&0&0&p \end{matrix}\right) \end{equation}\]

Eq. (11) and (12) are in the spherical coordinate with angular symmetry. Thus \(\eta = (-1,1,r^2,r^2\sin^2\theta)\) and \(u = \gamma(1,\beta,0,0)\). We can now calculate the derivative of the EM Tensor,

\[\begin{aligned} &\partial_0T^{00} = \frac{\partial (e+p)\gamma^2}{\partial t} - \frac{\partial p}{\partial t}\\ &\partial_0T^{0r} = \frac{\partial (e+p)\gamma^2 \beta}{\partial t}\\ &\partial_r T^{r0} = \frac{1}{r^2 } \frac{\partial}{\partial r } \left[r^2 (e+p)\gamma^2\beta\right]\\ &\partial_r T^{rr} = \frac{1}{r^2 } \frac{\partial}{\partial r } \left[r^2 (e+p)\gamma^2\beta^2\right] + \frac{1}{r^2 } \frac{\partial}{\partial r } \left(r^2p\right) \end{aligned}\]

From the energy and momentum conservation, i.e. \(\nabla _{\mu} T^{\mu\nu} = 0\). Sum up the corresponding derivatives, for \(\nu = 0\), we get

\[\begin{equation} \frac{\partial}{\partial t } \frac{e+p\beta^2}{1-\beta^2}+\frac{1}{r^2 } \frac{\partial}{\partial r }\frac{r^2(e+p)\beta}{1-\beta^2}=0, \end{equation}\]

same as eq. (11).

For \(\nu=1\),

\[\begin{equation} \frac{\partial }{\partial t}\frac{(e+p)\beta}{1-\beta^2}+ \frac{1}{r^2 } \frac{\partial}{\partial r } \frac{r^2 (e+p)\beta^2}{1-\beta^2} + \frac{\partial p}{\partial r} + \frac{2p}{r}. \end{equation}\]

There seems to be an additional \({2p}/{r}\) compared to eq. (12) of BM76.

With \(e = 3p\) the ultra-relativistic approximation for equation of state, and using \(\gamma = 1/\sqrt{1-\beta^2}\), we expand eq (11-12,14-15).

\[\begin{equation} \begin{aligned} (4\gamma^2-1)\frac{\partial p}{\partial t} + 8\gamma p \frac{\partial \gamma}{\partial t } + 4\gamma^2\beta\frac{\partial p }{\partial r} + 4p(2\gamma\beta + \frac{1}{\gamma \beta})\frac{\partial \gamma }{\partial r }+\frac{8p\gamma^2 \beta}{r}= 0&\hspace{2em}\text{(eq. 11)}\\ 4 \gamma^2 \beta \frac{\partial p }{\partial t } + 4p(2\gamma \beta + \frac{1}{\gamma \beta})\frac{\partial \gamma }{\partial t } + (4 \gamma^2 -3 )\frac{\partial p}{ \partial r } + 8p\gamma \frac{\partial \gamma }{\partial r} + \frac{8p(\gamma^2-1)}{r } = 0&\hspace{2em}\text{(eq. 12)}\\ \\ (\gamma^2-1)\frac{\partial p}{\partial t} + 4 \gamma p \frac{\partial \gamma}{\partial t } + 4 \beta p\gamma \frac{\partial \gamma}{\partial r} + \beta \gamma^2 \frac{\partial p }{\partial r }= 0&\hspace{2em}\text{(eq. 14)}\\ 3\gamma^2 \frac{\partial p}{\partial t} + 4 \gamma p \frac{\partial \gamma }{\partial t } + 3\beta \gamma^2\frac{\partial p}{\partial r } + (4\gamma \beta p +\frac{4p}{\gamma\beta})\frac{\partial \gamma}{\partial r} + \frac{8\gamma^2 p\beta}{r} = 0&\hspace{2em}\text{(eq. 15)} \end{aligned} \end{equation}\]

No additional approximations are made in the expansion above. Eq. 14 could be derived by combining eq. 11-12 and eliminate the \(p/r\) term. Eq. 15 is obtained by adding eq. 14 with eq. 11.

Therefore,

\[\begin{equation} \begin{aligned} \frac{d}{dt}(p\gamma^4) = &\gamma^2 \frac{\partial p }{\partial t}&\hspace{2em}\text{(eq. 14)}\\ \frac{d}{dt}\ln(p^3\gamma^4) = &-\frac{4}{r^2} \frac{\partial}{\partial r}(r^2 \beta)&\hspace{2em}\text{(eq. 15)} \end{aligned} \end{equation}\]

From the particle conservation \(\frac{\partial n'}{\partial t} = -\frac{1}{r^2}\frac{\partial }{\partial r }(r^2 n' \beta),\) we could yield \(\frac{dn'}{dt} = -\frac{n'}{r^2}\frac{\partial r^2\beta}{\partial r} = \frac{n'}{4}\frac{d}{dt}\ln(p^3\gamma^4)\)

Thus

\[\frac{4}{n'}\frac{dn'}{dt} = \frac{n'}{4}\frac{d}{dt}\ln(p^3\gamma^4)\]

which turns into \(\begin{equation} \frac{d}{dt}\left(\frac{p^3\gamma^4}{n'^4}\right) = \frac{d}{dt}\left(\frac{p}{n^{4/3}} \right) = 0 \end{equation}\)

Self-Similar Solution of Adiabatic Blast Wave

The energy between R0 and R1 is \(E = \int_{R_0}^{R_1} 16 \pi p\gamma^2 r^2 dr\). In the coefficient,\(4\pi r^2\) is due to spherical coordinate, \(4\gamma^2p\) is the shocked energy density. Assuming an adiabatic expanding blast wave, the total energy is conserved \(E \propto R^3\Gamma^2\).

Further assume \(\Gamma \propto t^{-m}\) with \(m > -1\), first expand \(\beta\) with \(\Gamma\) to \(O(\Gamma^{-2})\)

\[R = \int \beta dt = t\left[1-\frac{1}{2(m+1)\Gamma^2}\right].\]

BM76 chose the self similar variable as

\[\begin{equation} \chi = [1+2(m+1) \Gamma^2]\left(1-\frac{r}{t}\right) \end{equation}\]

The post-shock fluid could write as

\[\begin{equation} \begin{aligned} p =& \frac{2}{3}h_1\Gamma^2 f(\chi)\\ \gamma_{21}^2 =& \frac{1}{2}\Gamma g(\chi)\\ n' =& 2n_1\Gamma^2 h(\chi) \end{aligned}, \end{equation}\]

where \(\Gamma = \gamma_{1s}\) is the shock Lorentz factor and \(\Gamma^2 = 2 \gamma_{21}^2\) in the ultra-relativistic limit. \(n'\) is seen from the upstream rest frame and \(n'= \gamma n_2\).

Then we can substitute the old variables \(r,t\) with new variables \(\Gamma^2, \chi\). Be careful here that \(\Gamma^2 \propto t^{-m}\) and is independent upon \(r\). Then we have the elements of Jacobi Matrix,

\[\begin{aligned} t\frac{\partial \Gamma^2}{\partial t} =& -m\Gamma^2\\ t\frac{\partial \Gamma^2 }{\partial r} = &0\\ t\frac{\partial\chi}{\partial t} =& (m+1)(2\Gamma^2-\chi)+1\\ t\frac{\partial\chi}{\partial r} =& -[1+2(m+1)\Gamma^2] \end{aligned}\]

Specially, we have substituted \(\Gamma^2 =A t^{-m}\) to compute \(t\frac{\partial\chi}{\partial t}\). The old derivatives are now

\[\begin{equation} \begin{aligned} &\frac{\partial}{\partial\ln t} = -m\frac{\partial}{\partial \ln \Gamma^2}+ [(m+1)(2\Gamma^2-\chi)+1]\frac{\partial }{\partial\chi}\\ &t\frac{\partial}{\partial r} = -[1+2(m+1)\Gamma^2]\frac{\partial }{\partial\chi}\\ &\frac{d}{d\ln t} = -m\frac{\partial}{\partial \ln \Gamma^2}+ (m+1)\left(\frac{2}{g}-\chi\right)\frac{\partial }{\partial\chi} \end{aligned} \end{equation}\]

The third equation is from convective derivative \(d/dt = \partial_t + \beta_{21}\partial_r\). There is the function \(g\) because of the \(\beta_{21}\).

Now we have the new derivatives. Just substitute them into (14-16), and only keep the terms with highest order of \(\Gamma^2\). Solve the linear equations and we get the self similar solution of adiabatic blast wave

\[\begin{equation} \begin{aligned} \frac{1}{g}\frac{d \ln f}{d\chi} =& \frac{8(m-1)-(m-4)g\chi}{(m+1)(4-8g\chi+g^2\chi^2)}\\ \frac{1}{g}\frac{d \ln g}{d\chi} =& \frac{(7m-4)-(m+2)g\chi}{(m+1)(4-8g\chi+g^2\chi^2)}\\ \frac{1}{g}\frac{d \ln h}{d\chi} =& \frac{2(9m-8)-2(5m-6)g\chi + (m-2) g^2\chi^2}{(m+1)(4-8g\chi+g^2\chi^2)(2-g\chi)} \end{aligned} \end{equation}\]

There is no \(h\) in eq (14,15), so we could first align these two equations and solve derivative of \(f,g\).