Energy Conservation in the Thin Layer Approximation: II. The Asymmetric Classic Case for Supernovae Remnant

Lorenzo Zaninetti

Physics Department, via P.Giuria 1, I-10125 Turin, Italy.

DOI: 10.4236/ijaa.2020.102009                

Abstract

Here we extend the conservation of energy in the framework of the thin layer approximation to the asymmetrical case. Four types of interstellar mediums are analysed, in which the density follows an inverse square profile, a power law profile, an exponential profile and a toroidal profile. An analytical solution for the radius as a function of time and the polar angle in spherical coordinates is derived in the case of the inverse square profile. The analytical and numerical results are applied to two supernova remnants: SN 1987A and SN 1006. The back reaction due to the radiative losses is evaluated in the case of the inverse square profile for the surrounding medium. Two models for the image formation are presented, which explain the triple ring visible in SN 1987A and the jet feature of SN 1006.

Keywords

Supernovae, General ISM, Supernova Remnants Supernovae, Individual (SN 1987A) Supernovae, Individual (SN 1006)

1. Introduction

The asymmetries observed in supernovae (SNs) and supernova remnants (SNRs) have been recently investigated by different approaches: on analysing a type Ia SN that exploded around CE 1900 [1] reported that it is asymmetric in the radio region but shows a bipolar morphology in the X region. The details of the deceleration of that SN in different directions were analysed by [2]. The symmetry of many SNRs in 6-cm and 20-cm Very Large Array images was measured by [3]. A three dimensional smoothed particle hydrodynamic simulation was used by [4] in order to explain the observed asymmetries in SNs. [5] implemented a three-dimensional MHD numerical code in which the initial mass distribution of the SN is asymmetric. The presence of an asymmetric circumstellar medium (CSM) responsible for the lack of narrow lines within the first two days of the explosion for the Type II-P supernova SN 2017gmr was suggested by [6]. The differences in the Si II line widths in type Ia supernovae with an asymmetric CSM were explained by [7]. An asymmetric explosion as seen from different viewing angles was suggested by [8] in order to explain some anomalies in the photospheric velocities. The model of energy conservation in the thin layer approximation for the expansion of SNRs was developed in [9] in the spherical case. Here we analyse the asymmetrical or non-spherical case. In order to accomplish this, Section 2 derives the basic differential equations which regulate the equation of motion for SNRs for four types of medium. Section 3 applies those analytical and numerical results to SN 1987A and SN 1006 and Section 4 analyses two models of the image formation for the two SNRs analysed.

2. Energy Conservation

A point in Cartesian coordinates is characterized by x, y and z, and the position of the origin is the centre of the explosive phenomena. The same point in spherical coordinates is characterized by the radial distance $r\in \left[\mathrm{0,}\infty \right]$, the polar angle $\theta \in \left[\mathrm{0,}\pi \right]$, and the azimuthal angle $\phi \in \left[\mathrm{0,2}\pi \right]$. In the following, the profiles of the density considered are independent of the azimuthal angle and are functions of the distance z from the $z=0$ plane. In spherical coordinates, the goal is to derive the instantaneous radius of expansion, r, as a function of the polar angle. The following approaches will model the CSM around the point of explosion of the SN which is characterized by a plane, $z=0$, which produces an up and down symmetry in the polar angle, $r\left(\theta \right)=r\left(\pi -\theta \right)$, and by a polar axis, $x=0,y=0$, which produces a right and left symmetry for the azimuthal angle $r\left(\phi \right)=r\left(\phi +\pi \right)$. The basic symmetries are now outlined

$r\left(\theta \right)=r\left(\pi -\theta \right)$ (1a)

$r\left(\phi \right)=r\left(\phi +\pi \right)$ (1b)

In other words, the numerical simulations will be developed in the first quadrant ( $\theta \in \left[\mathrm{0,}\frac{\pi }{2}\right]$, $\phi =0$ and then applied to the other three quadrants.

The conservation of kinetic energy in spherical coordinates along the solid angle $\Delta \Omega$ in the framework of the thin layer approximation is

$\frac{1}{2}{M}_0\left(r_0\mathrm{<mo>\; ;\; </mo>}\theta \right)v_0^2=\frac{1}{2}M\left(r\mathrm{<mo>\; ;\; </mo>}\theta \right)v^2\mathrm{,}$ (2)

where $M_0\left(r_0\mathrm{<mo>\; ;\; </mo>}\theta \right)$ and $M\left(r\mathrm{<mo>\; ;\; </mo>}\theta \right)$ are the swept masses at $r_0$ and r, while $v_0$ and v are the velocities of the thin layer at $r_0$ and r. The above conservation law when written as a differential equation is

$\frac{1}{2}M\left(r\right){\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2-\frac{1}{2}{M}_0{v}_0^2=0.$ (3)

The velocity as a function of the radius is

$v\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}\theta \right)=\frac{r_0^{3/2}{v}_0}{r^{3/2}}\mathrm{.}$ (4)

In the following, $r_0$ is the radius after which the density starts to decrease, r is the radius of expansion in spherical coordinates, ${\rho }_c$ is the density at $r=0$, $z=rcos\left(\theta \right)$ is the Cartesian coordinate along the Z axis and $cos\left(\theta \right)$ is the polar angle in spherical coordinates. The main astrophysical assumption adopted here is that the density of the CSM decreases as a function of the distance from the centre, due to the previous stellar winds. At the time of writing, there are no astronomical observations which outline such an effect, and therefore we test different profiles for the density and evaluate whether they are compatible with the observed sections of SNRs.

2.1. An Inverse Square Profile of the Density

We assume that the medium around the SN scales with the axial piecewise dependence

$\rho \left(r;r_0\right)=\{\begin{array}{l}{\rho }_c\text{if}r\le r_0\\ {\rho }_c{\left(\frac{r_0}{z}\right)}^2\text{if}r>r_0\end{array}$ (5)

The mass $M_0$ swept in the interval $\left[0,r_0\right]$ is

$M_0\left({\rho }_c\mathrm{,}{r}_0\right)=\frac{4}{3}{\rho }_c\pi r_0^3\mathrm{.}$

The total mass $M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{\rho }_c\right)$ swept in the interval [0, r] is

$M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{\rho }_c\mathrm{,}{z}_0\right)=\frac{1}{3}{\rho }_c{r}_0^3+\frac{{\rho }_c{z}_0^2\left(r-r_0\right)}{{\left(cos\left(\theta \right)\right)}^2}\mathrm{.}$

The positive solution of Equation (2) gives the velocity as a function of the radius:

$v\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{z}_0\right)=\frac{\sqrt{\left(r_0^3{\left(cos\left(\theta \right)\right)}^2+3z_0^{2}r-3z_0^2{r}_0\right)r_0}cos\left(\theta \right)v_0{r}_0}{r_0^3{\left(cos\left(\theta \right)\right)}^2+3z_0^{2}r-3z_0^2{r}_0}\mathrm{.}$ (6)

The differential equation which models the energy conservation is

$\frac{1}{2}\left(\frac{1}{3}{\rho }_c{r}_0^3+\frac{{\rho }_c{z}_0^2\left(r\left(t\right)-r_0\right)}{{\left(cos\left(\theta \right)\right)}^2}\right){\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2-\frac{1}{6}{\rho }_c{r}_0^3{v}_0^2=0.$ (7)

The analytical solution of the above differential equation is

$\begin{array}{l}r\left(t\mathrm{<mo>\; ;\; </mo>}{z}_0\mathrm{,}{r}_0\mathrm{,}\theta \right)\\ =\frac{-r_0\left(-{\left(cos\left(\theta \right)\right)}^{2/3}{\left(9z_0^2\left(t-t_0\right)v_0+2r_0^3{\left(cos\left(\theta \right)\right)}^2\right)}^{2/3}\sqrt[3]{2}+2r_0^2{\left(cos\left(\theta \right)\right)}^2-6z_0^2\right)}{6z_0^2}\mathrm{.}\end{array}$ (8)

The above inverse square profile for density satisfies the basic symmetries as outlined in Equations (1a) and (1b).

2.2. Back Reaction for the Inverse Square Profile

The radiative losses per unit length are assumed to scale as

$-4\epsilon {\rho }_s{v}^2{r}^2\mathrm{,}$ (9)

where $\epsilon$ is a constant and $rh{o}_s$ is the density in the thin advancing layer which is $4\rho$. Inserting into the above equation the velocity to first order as given by Equation (6), the radiative losses $Q\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}{z}_0\mathrm{,}\epsilon \mathrm{,}\theta \right)$ are

$Q\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}{z}_0\mathrm{,}\epsilon \mathrm{,}\theta \right)=-4\frac{\epsilon {\rho }_c{z}_0^2{r}_0^3{v}_0^2}{r_0^3{\left(cos\left(\theta \right)\right)}^2+3z_0^{2}r-3z_0^2{r}_0}\mathrm{.}$ (10)

The sum of the radiative losses between $r_0$ and r is given by the integral

$\begin{array}{l}L\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}{z}_0\mathrm{,}\epsilon \mathrm{,}\theta \right)={\displaystyle {\int }_{r_0}^{r}Q\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}{z}_0\mathrm{,}\epsilon \mathrm{,}\theta \right)\text{d}r}\\ =4/3\epsilon {\rho }_c{r}_0^3{v}_0^2\left(ln\left(r_0^3{\left(cos\left(\theta \right)\right)}^2+3z_0^2\left(-r_0+r\right)\right)-ln\left(r_0^3{\left(cos\left(\theta \right)\right)}^2\right)\right)\mathrm{.}\end{array}$ (11)

The conservation of energy in the presence of the back reaction due to the radiative losses is

$\begin{array}{l}1/2\frac{v^2\left(1/3r_0^3{\left(cos\left(\theta \right)\right)}^2+z_0^2\left(-r_0+r\right)\right){\rho }_c}{{\left(cos\left(\theta \right)\right)}^2}\\ +4/3\epsilon {\rho }_c{r}_0^3{v}_0^{2}ln\left(r_0^3{\left(cos\left(\theta \right)\right)}^2+3z_0^2\left(-r_0+r\right)\right)\\ -4/3\epsilon {\rho }_c{r}_0^3{v}_0^2\left(3ln\left(r_0\right)+ln\left({\left(cos\left(\theta \right)\right)}^2\right)\right)=1/6{\rho }_c{r}_0^3{v}_0^2\mathrm{.}\end{array}$ (12)

The analytical solution for the velocity to second order, $v_c\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}{z}_0\mathrm{,}\epsilon \mathrm{,}\theta \right)$, is

$\begin{array}{l}{v}_c\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}{z}_0\mathrm{,}\epsilon \mathrm{,}\theta \right)\\ =\frac{\sqrt{-A{r}_0\left(8ln\left(A\right)\epsilon -8ln\left(r_0^3{\left(cos\left(\theta \right)\right)}^2\right)\epsilon -1\right)}{v}_{0}cos\left(\theta \right)r_0}{A}\mathrm{,}\end{array}$ (13)

where

$A=r_0^3{\left(cos\left(\theta \right)\right)}^2+3z_0^{2}r-3z_0^2{r}_0\mathrm{.}$ (14)

The inclusion of the back reaction allows the evaluation of the SRS’s maximum length, $r_{back}\left(r_0\mathrm{,}{z}_0\mathrm{,}\theta \mathrm{,}\epsilon \right)$, which can be derived by setting this velocity to zero:

$r_{back}\left(r_0\mathrm{,}{z}_0\mathrm{,}\theta \mathrm{,}\epsilon \right)=\frac{{\text{e}}^{1/8\frac{8ln\left(r_0^3{\left(cos\left(\theta \right)\right)}^2\right)\epsilon +1}{\epsilon }}-r_0^3{\left(cos\left(\theta \right)\right)}^2+3z_0^2{r}_0}{3z_0^2}\mathrm{.}$ (15)

Figure 1 shows the SRS’s maximum length as a function of the constant of conversion $\epsilon$ and the polar angle $\theta$. The above figure shows that the SNR cannot reach an infinite length at an infinite time but has a finite lifetime and length.

Figure 1. A 2D map of $r_{back}\left(r_0\mathrm{,}{z}_0\mathrm{,}\theta \mathrm{,}\epsilon \right)$ in pc. The parameters are $r_0=0.098\text{pc}$ and $z_0=1.5\times {10}^{-4}\text{pc}$.

2.3. A Power Law Profile of the Density

The medium is assumed to scale as

$\rho \left(r;r_0\right)=\{\begin{array}{l}{\rho }_c\text{if}r\le r_0\\ {\rho }_c{\left(\frac{r_0}{z}\right)}^{\alpha }\text{if}r>r_0\end{array}$ (16)

where $\alpha$ is a positive real number. The total mass $M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{\rho }_c,\alpha \right)$ swept in the interval $\left[0,r\right]$ is

$\begin{array}{l}M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{\rho }_c\mathrm{,}{z}_0\mathrm{,}\alpha \right)\\ =\frac{1}{3}{\rho }_c{r}_0^3-\frac{{\rho }_c\left(r^{-\alpha +3}{z}_0^{\alpha }{\left({\left(cos\left(\theta \right)\right)}^{-1}\right)}^{\alpha }-r_0^{-\alpha +3}{z}_0^{\alpha }{\left({\left(cos\left(\theta \right)\right)}^{-1}\right)}^{\alpha }\right)}{\alpha -3}\mathrm{.}\end{array}$

The velocity $v\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}\alpha \right)$ as a function of the radius is

$\begin{array}{l}v\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}\alpha \right)\\ =-\frac{\sqrt{-3r_0\left(\alpha -3\right)\left(-r_0^{-\alpha +3}{z}_0^{\alpha }{\left({\left(cos\left(\theta \right)\right)}^{-1}\right)}^{\alpha }+r^{-\alpha +3}{z}_0^{\alpha }{\left({\left(cos\left(\theta \right)\right)}^{-1}\right)}^{\alpha }-\frac{1}{3}{r}_0^3\left(\alpha -3\right)\right)}{v}_0{r}_0}{3r_0^{-\alpha +3}{z}_0^{\alpha }{\left({\left(cos\left(\theta \right)\right)}^{-1}\right)}^{\alpha }-3r^{-\alpha +3}{z}_0^{\alpha }{\left({\left(cos\left(\theta \right)\right)}^{-1}\right)}^{\alpha }+r_0^3\left(\alpha -3\right)}\mathrm{.}\end{array}$ (17)

The differential equation which governs the motion is

$\begin{array}{l}\frac{1}{2}\left(\frac{1}{3}{\rho }_c{r}_0^3-\frac{{\rho }_c\left({\left(r\left(t\right)\right)}^{-\alpha +3}{z}_0^{\alpha }{\left({\left(cos\left(\theta \right)\right)}^{-1}\right)}^{\alpha }-r_0^{-\alpha +3}{z}_0^{\alpha }{\left({\left(cos\left(\theta \right)\right)}^{-1}\right)}^{\alpha }\right)}{\alpha -3}\right){\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2\\ -\frac{1}{6}{\rho }_c{r}_0^3{v}_0^2=0.\end{array}$ (18)

The above differential equation does not have an analytical solution; a Taylor expansion of order 3 covers a limited range in time

$r\left(t\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}{z}_0\mathrm{,}{t}_0\mathrm{,}\alpha \right)=r_0+v_0\left(t-t_0\right)-\frac{3}{4}{r}_0^{-1-\alpha }{v}_0^2{z}_0^{\alpha }{\left({\left(cos\left(\theta \right)\right)}^{-1}\right)}^{\alpha }{\left(t-t_0\right)}^2\mathrm{.}$ (19)

The above power law profile for the density satisfies the basic symmetries as outlined in Equations (1a) and (1b).

2.4. An Exponential Profile of the Density

The medium is assumed to scale as

$\rho \left(r;r_0,z_0\right)=\{\begin{array}{l}{\rho }_c\text{if}r\le r_0\\ {\rho }_c\exp \left(-\frac{z}{z_0}\right)\text{if}r>r_0\end{array}$ (20)

The total mass $M\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{\rho }_c\mathrm{,}{z}_0\right)$ swept in the interval [0, r] is

$\begin{array}{c}M\left(r;r_0,{\rho }_c,z_0\right)=\frac{1}{3{\left(\cos \left(\theta \right)\right)}^3}{\rho }_c(r_0^3{\left(\cos \left(\theta \right)\right)}^3+3{\left(\cos \left(\theta \right)\right)}^2{\text{e}}^{-\frac{r_0\cos \left(\theta \right)}{z_0}}{r}_0^2{z}_0\\ -3{\left(\cos \left(\theta \right)\right)}^2{\text{e}}^{-\frac{r\cos \left(\theta \right)}{z_0}}{r}^2{z}_0+6\cos \left(\theta \right){\text{e}}^{-\frac{r_0\cos \left(\theta \right)}{z_0}}{r}_0{z}_0^2\\ -6\cos \left(\theta \right){\text{e}}^{-\frac{r\cos \left(\theta \right)}{z_0}}r{z}_0^2+6{\text{e}}^{-\frac{r_0\cos \left(\theta \right)}{z_0}}{z}_0^3-6{\text{e}}^{-\frac{r\cos \left(\theta \right)}{z_0}}{z}_0^3).\end{array}$ (21)

The velocity $v\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\right)$ as a function of the radius is

$\begin{array}{c}v\left(r\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\right)=\frac{1}{3{\left(\cos \left(\theta \right)\right)}^3}(r_0^3{\left(\cos \left(\theta \right)\right)}^3+3{\left(\cos \left(\theta \right)\right)}^2{\text{e}}^{-\frac{r_0\cos \left(\theta \right)}{z_0}}{r}_0^2{z}_0\\ -3{\left(\cos \left(\theta \right)\right)}^2{\text{e}}^{-\frac{r\cos \left(\theta \right)}{z_0}}{r}^2{z}_0+6\cos \left(\theta \right){\text{e}}^{-\frac{r_0\cos \left(\theta \right)}{z_0}}{r}_0{z}_0^2\\ -6\cos \left(\theta \right){\text{e}}^{-\frac{r\cos \left(\theta \right)}{z_0}}r{z}_0^2+6{\text{e}}^{-\frac{r_0\cos \left(\theta \right)}{z_0}}{z}_0^3-6{\text{e}}^{-\frac{r\cos \left(\theta \right)}{z_0}}{z}_0^3){\rho }_c.\end{array}$ (22)

The differential equation which governs the motion is

$\begin{array}{l}\frac{1}{6{\left(cos\left(\theta \right)\right)}^3}3(-2\left(\frac{1}{2}{\left(cos\left(\theta \right)\right)}^2{\left(r\left(t\right)\right)}^2+r\left(t\right)cos\left(\theta \right)z_0+z_0^2\right)\\ \times {\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2{z}_0{\text{e}}^{-\frac{r\left(t\right)cos\left(\theta \right)}{z_0}}\\ +z_0{\left(\frac{\text{d}}{\text{d}t}r\left(t\right)\right)}^2\left({\left(cos\left(\theta \right)\right)}^2{r}_0^2+2r_0{z}_{0}cos\left(\theta \right)+2z_0^2\right){\text{e}}^{-\frac{r_{0}cos\left(\theta \right)}{z_0}}\\ -\frac{1}{3}{r}_0^3{\left(cos\left(\theta \right)\right)}^3\left(v_0-\frac{\text{d}}{\text{d}t}r\left(t\right)\right)\left(v_0+\frac{\text{d}}{\text{d}t}r\left(t\right)\right)){\rho }_c=0.\end{array}$ (23)

In the absence of an analytical solution, we present the Taylor expansion of order 3 for the trajectory.

$r\left(t\mathrm{<mo>\; ;\; </mo>}{r}_0\mathrm{,}{v}_0\mathrm{,}{z}_0\mathrm{,}{t}_0\right)=r_0+v_0\left(t-t_0\right)-\frac{3}{4}\frac{v_0^2{\left(t-t_0\right)}^2}{r_0}{\text{e}}^{-\frac{r_{0}cos\left(\theta \right)}{z_0}}\mathrm{.}$ (24)

The above exponential profile for the density satisfies the basic symmetries as outlined in Equations (1a) and (1b).

2.5. A Toroidal Profile of the Density

A torus is swept out by revolving a small circle of radius $r_T$ about an axis lying in the same plane as the circle but outside it, say at a distance of $R_T$, see Figure 2. The above toroidal profile of the density, once the axis of revolution is identified with the polar axis, satisfies the basic symmetries as outlined in Equations (1a) and (1b). In Cartesian coordinates, the torus satisfies the equation

${\left(R_T-\sqrt{x^2+y^2}\right)}^2+z^2-r_T^2=0.$ (25)

We now consider its intersection with the plane $y=0$ and as a consequence the equation of the torus is now

${\left(R_T-x\right)}^2+z^2-r_T^2=0.$ (26)

We now evaluate the intersection between the torus and a straight line of equation which crosses the centre, ( $z=0,x=0$ ),

$z=xcot\left(\theta \right)\mathrm{,}$ (27)

where $\theta$ is the polar angle in the spherical coordinates for the circle which represents the torus in the first quadrant. The line touches the torus at a single point at the critical value of the polar angle, ${\theta }_{crit}$ :

Figure 2. A 3D view of a torus.

${\theta }_{crit}=arctan\left(\frac{\sqrt{R_T^2-r_T^2}}{r_T}\right)\mathrm{.}$ (28)

The above angle allows us to define the following three zones.

1) $\theta <{\theta }_{crit}$ : the line does not intersect the torus.

2) $\theta ={\theta }_{crit}$ : the line is tangent to the torus at one degenerate single point, with Cartesian coordinates ( $x_{crit}\mathrm{,}{z}_{crit}$ ), see Figure 3.

3) $\theta >{\theta }_{crit}$ : the line intersects the torus at two points, with Cartesian coordinates ( $x_1\mathrm{,}{z}_1$ ) and ( $x_2\mathrm{,}{z}_2$ ), see Figure 4.

The Cartesian coordinates at $\theta ={\theta }_{crit}$ are

$x_{crit}=\frac{R_T^2-r_T^2}{R_T}$ (29a)

$z_{crit}=\frac{\sqrt{R_T^2-r_T^2}{r}_T}{R_T}\mathrm{.}$ (29b)

The Cartesian coordinates of the two intersections when $\theta >{\theta }_{crit}$ are

$x_1=\left(R_{T}sin\left(\theta \right)-\sqrt{-R_T^2{\left(cos\left(\theta \right)\right)}^2+r_T^2}\right)sin\left(\theta \right)$ (30a)

$z_1=\left(R_T\sin \left(\theta \right)-\sqrt{-R_T^2{\left(\cos \left(\theta \right)\right)}^2+r_T^2}\right)\cos \left(\theta \right)$ (30b)

Figure 3. The straight line which intersects the centre (red) and is tangent to the torus (blue) in the first quadrant.

Figure 4. The intersections between a straight line which intersects the centre (red), the torus (blue) to the tangent point (asterisk).

$x_2=\left(R_T\sin \left(\theta \right)+\sqrt{-R_T^2{\left(\cos \left(\theta \right)\right)}^2+r_T^2}\right)\sin \left(\theta \right)$ (30c)

$z_2=\left(R_T\sin \left(\theta \right)+\sqrt{-R_T^2{\left(\cos \left(\theta \right)\right)}^2+r_T^2}\right)\cos \left(\theta \right)$ (30d)

We now assume that the density of matter is ${\rho }_c$ outside the torus and ${\rho }_1$ inside the torus. The mass swept in the third zone requires a careful analysis. When $\theta >{\theta }_{crit}$, the mass swept along a line before the intersection with the torus, $M_I\left(r\mathrm{<mo>\; ;\; </mo>}{\rho }_c\right)$, is

$M_I\left(r\mathrm{<mo>\; ;\; </mo>}{\rho }_c\right)=\frac{1}{3}{\rho }_c{r}^3\mathrm{,}$ (31)

where r is the momentary radius of expansion in spherical coordinates.

The mass swept when r is inside the torus in the third zone, $M_{II}\left(r\mathrm{<mo>\; ;\; </mo>}{\rho }_c\mathrm{,}{\rho }_1\mathrm{,}{R}_T\mathrm{,}{r}_T\right)$, is

$\begin{array}{l}{M}_{II}\left(r\mathrm{<mo>\; ;\; </mo>}{\rho }_c\mathrm{,}{\rho }_1\mathrm{,}{R}_T\mathrm{,}{r}_T\right)\\ =\frac{1}{3}{\rho }_c{\left({\left(R_{T}sin\left(\theta \right)-\sqrt{-R_T^2{\left(cos\left(\theta \right)\right)}^2+r_T^2}\right)}^2\right)}^{3/2}\\ +\frac{1}{3}{\rho }_1\left(r^3-{\left({\left(R_{T}sin\left(\theta \right)-\sqrt{-R_T^2{\left(cos\left(\theta \right)\right)}^2+r_T^2}\right)}^2\right)}^{3/2}\right)\mathrm{.}\end{array}$ (32)

The mass swept when r is outside the torus in the third zone is

$\begin{array}{l}{M}_{III}\left(r\mathrm{<mo>\; ;\; </mo>}{\rho }_c\mathrm{,}{\rho }_1\mathrm{,}{R}_T\mathrm{,}{r}_T\right)\\ =\frac{1}{3}{\rho }_c{\left({\left(R_{T}sin\left(\theta \right)-\sqrt{-R_T^2{\left(cos\left(\theta \right)\right)}^2+r_T^2}\right)}^2\right)}^{3/2}\\ +\frac{1}{3}{\rho }_1{\left({\left(R_{T}sin\left(\theta \right)+\sqrt{-R_T^2{\left(cos\left(\theta \right)\right)}^2+r_T^2}\right)}^2\right)}^{3/2}\\ -{\left({\left(R_{T}sin\left(\theta \right)-\sqrt{-R_T^2{\left(cos\left(\theta \right)\right)}^2+r_T^2}\right)}^2\right)}^{3/2}\mathrm{<mo\; stretchy="false">\; )\; </mo>}\\ +\frac{1}{3}{\rho }_c\left(r^3-{\left({\left(R_{T}sin\left(\theta \right)+\sqrt{-R_T^2{\left(cos\left(\theta \right)\right)}^2+r_T^2}\right)}^2\right)}^{3/2}\right)\mathrm{.}\end{array}$ (33)

The equation of motion is solved through the Euler method when the following recursive equations are solved

$r_{n+1}=r_n+v_n\Delta t$ (34a)

$v_{n+1}=v_n{\left(\frac{M_n\left(r_n\right)}{M_{n+1}\left(r_{n+1}\right)}\right)}^{1/2}\mathrm{,}$ (34b)

where $r_n$, $v_n$, $M_n$ are the temporary radius, velocity, and total mass, respectively, $\Delta t$ is the time step and n is the index. Due to the fact that the velocity is continuously updated, see Equation (34b), the method turns out to be stable.

3. Astrophysical Applications

The complex structure of SN 1987A, see the Hubble Space Telescope (ST) image in Figure 5, can be classified as a torus only, a torus plus two lobes, and a torus plus 4 lobes, see [10] [11]; we therefore speak of a strong asymmetry. The region connected with the radius of the advancing torus is here identified with our

equatorial region, in spherical coordinates, $\theta =\frac{\pi }{2}$. A useful resource for calibration

for calibration is the geometric section of SN 1987A which is given as a sketch in Fig. 5 of [12]. This geometric section was digitized and rotated in the x-z plane by −40˚, see Figure 6; it will be the astronomical section of reference in order to test the simulations.

SN 1006 started to be visible in 1006 AD and currently has a radius of 12.2 pc, see [13] [14]. The X-shape is shown in Figure 7 and the γ-ray shape (100 GeV) is shown in Figure 8.

More precisely, on referring to the above X-image, it can be observed that the radius is greatest in the north-east direction, see also the radio map of SN 1006 at 1370 MHz by [15], and the X-map in the 0.4 - 5.0 keV band of Fig. 1 in [13]. The following observed radii can be extracted: $R_{up}=11.69\text{pc}$ in the polar direction and $R_{eq}=8.7\text{pc}$ in the equatorial direction. A geometric section of the above X-map was digitized and rotated in the x-z plane by −45˚, see Figure 9; it will be the test for the simulations based on the previous data we can speak of a weak asymmetry for SN 1006. The reliability of the models is evaluated in terms of the percentage reliability, ${\epsilon }_{\text{obs}}$,

Figure 5. An ST image of SN 1987A in the year 1997. Credit is given to the Hubble Space Telescope.

Figure 6. Geometric section of SN 1987A in the x-z plane adapted by the author from Fig. of [12].

Figure 7. A Chandra X-ray (Red, Green, Blue) image of SN 1006 in the year 2013. Credit is given to the Chandra X-ray Observatory.

Figure 8. An high energy stereoscopic system (HESS) γ-ray image of SN 1006. Credit is given to the Max-Planck-Institut für Kernphysik.

Figure 9. Geometric section of SN 1006 in the x-z plane adapted by the author from our Figure 7.

${\epsilon }_{\text{obs}}=100\left(1-\frac{{\displaystyle {\sum }_j{\left|r_{\text{obs}}-r_{\text{num}}\right|}_j}}{{\displaystyle {\sum }_j{r}_{\text{obs}}{}_{\mathrm{,}j}}}\right)\mathrm{,}$ (35)

where $r_{\text{num}}$ is the theoretical radius of the SNR, $r_{\text{obs}}$ is the observed radius of the SNR, and the index j varies from 1 to the number of available observations.

3.1. Results for the Inverse Square Profile

In the case of an inverse square profile of the density we have an analytical solution as given by Equation (8). Figure 10 displays a cut of SN 1987A in the x-z plane.

A rotation around the z-axis of the previous geometric section allows building a 3D surface, see Figure 11.

Figure 10 displays a cut of SN 1987A in the x-z plane.

Figure 12 displays a cut of SN 1006 in the x-z plane and Figure 13 a 3D display.

The above results can be easily reproduced because we have an analytical expression for the trajectory.

3.2. Results for a Power Profile

In the case of a power law profile of the density we obtained a numerical solution for the differential Equation (18). Figure 14 presents the results for SN 1987A in the x-z plane and Figure 15 those for SN 1006.

The above results shows that a variable power law profile can model the observed sections of the two SNRS here simulated.

3.3. Results for an Exponential Profile

In the case of an exponential profile of the density we obtained a numerical solution for the differential Equation (23). Figure 16 presents the results for SN 1987A in the x-z plane and Figure 17 those for SN 1006.

Figure 10. Geometric section of SN 1987A in the x-z plane with an inverse square profile (green points) and the observed profile (red stars). The parameters $r_0=0.001355\text{pc}$, $z_0=0.0001\text{pc}$, $t=27.7\text{yr}$, $t_0=0.053\text{yr}$ and $v_0=25000\text{km}\cdot {\text{s}}^{-1}$ give ${\epsilon }_{\text{obs}}=\mathrm{91.92\; <mi>\; \%\; </mi>}$.

Figure 11. 3D surface of SN 1987A with parameters as in Figure 10 in the framework of an inverse square profile. The three Euler angles are $\Theta =40$, $\Phi =60$ and $\Psi =60$.

Figure 12. Geometric section of SN 1006 in the x-z plane with an inverse square profile (green points) and the observed profile (red stars). The parameters $r_0=0.098\text{pc}$, $z_0=1.5\times {10}^{-4}\text{pc}$, $t=1000\text{yr}$, $t_0=9.057\text{yr}$ and $v_0=10600\text{km}\cdot {\text{s}}^{-1}$ give ${\epsilon }_{\text{obs}}=\mathrm{94.34\; <mi>\; \%\; </mi>}$.

Figure 13. 3D surface of SN 1006 with parameters as in Figure 12, inverse square profile. The three Euler angles are $\Theta =70$, $\Phi =70$ and $\Psi =70$.

Figure 14. Geometric section of SN 1987A in the x-z plane with a power law profile (green points) and the observed profile (red stars). The parameters $\alpha =1.5$, $r_0=0.0038\text{pc}$, $z_0=0.0001\text{pc}$, $t=27.7\text{yr}$, $t_0=0.148\text{yr}$ and $v_0=25000\text{km}\cdot {\text{s}}^{-1}$ give ${\epsilon }_{\text{obs}}=\mathrm{91.92\; <mi>\; \%\; </mi>}$.

Figure 15. Geometric section of SN 1006 in the x-z plane with a power law profile (green points) and the observed profile (red stars). The parameters $\alpha =1.5$, $r_0=0.098\text{pc}$, $z_0=0.00015\text{pc}$, $t=1000\text{yr}$, $t_0=9.057\text{yr}$ and $v_0=10600\text{km}\cdot {\text{s}}^{-1}$ give ${\epsilon }_{\text{obs}}=\mathrm{93.562\; <mi>\; \%\; </mi>}$.

Figure 16. Geometric section of SN 1987A in the x-z plane with an exponential profile (green points) and the observed profile (red stars). The parameters $r_0=0.04\text{pc}$, $z_0=0.033\text{pc}$, $t=27.7\text{yr}$, $t_0=0.0528\text{yr}$ and $v_0=25000\text{km}\cdot {\text{s}}^{-1}$ give ${\epsilon }_{\text{obs}}=\mathrm{92.45\; <mi>\; \%\; </mi>}$.

Figure 17. Geometric section of SN 1006 in the x-z plane with an exponential profile (green points) and the observed profile (red stars). The parameters $r_0=2.45\text{pc}$, $z_0=4.9\text{pc}$, $t=1000\text{yr}$, $t_0=0.052\text{yr}$ and $v_0=25000\text{km}\cdot {\text{s}}^{-1}$ give ${\epsilon }_{\text{obs}}=\mathrm{94.26\; <mi>\; \%\; </mi>}$.

The above results shows that an exponential profile for the density is comparable with the observed sections of the two SNRS here simulated.

3.4. Results for a Toroidal Profile

In the case of a toroidal profile of the density, we obtained a numerical solution for the two recursive Equations (34a) and (34b). Figure 18 presents the results for SN 1987A in the x-z plane and Figure 19 for SN 1006.

A toroidal region around the point of the explosion with density greater than the surrounding medium produces theoretical sections which are comparable with the observed sections of the two SNRS here simulated.

4. Theory of the Image

The astronomical images are given by cut or 2D maps for the intensity of emission. The shape of the intensity of emission is a combination of different processes which are as follows.

1) An assumption on the transfer equation: we adopt an optically thin medium.

2) Type of emission: we choose non thermal emission.

3) Geometry of the radiative zone: the zone of emission resides in the thin advancing shell, which in our case is not spherical.

More details on the above assumptions can be found in Section 7 of [16].

Figure 18. Geometric section of SN 1987A in the x-z plane for a toroidal profile (green points) and the observed profile (red stars). The parameters $r_0=0.15\text{pc}$, $R_T=1.02r_0\text{pc}$, $r_T=0.9R_T\text{pc}$, ${\rho }_c=1$, ${\rho }_1=10{\rho }_c$, $t=27.7\text{yr}$, $t_0=5.65\text{yr}$, $v_0=26000\text{km}\cdot {\text{s}}^{-1}$ and $h=0.01\text{yr}$ give ${\epsilon }_{\text{obs}}=\mathrm{86.28\; <mi>\; \%\; </mi>}$.

Figure 19. Geometric section of SN 1006 in the x-z plane for a toroidal profile (green points) and the observed profile (red stars). The parameters $r_0=3.6\text{pc}$, $R_T=1.02r_0\text{pc}$, $r_T=0.86R_T\text{pc}$, ${\rho }_c=1$, ${\rho }_1=3{\rho }_c$, $t=1000\text{yr}$, $t_0=135.64\text{yr}$, $v_0=26000\text{km}\cdot {\text{s}}^{-1}$ and $h=1\text{yr}$ give ${\epsilon }_{\text{obs}}=\mathrm{94.41\; <mi>\; \%\; </mi>}$.

4.1. How to Build an Image

The first model assumes that the intensity of the image is proportional to the length along the line of sight within the emitting region. The derivation of this length can be complicated, due to the complex morphology of the analysed object and to the infinite points of view of the observer. Figure 20 presents a sketch of the emitting layer, some lines of sight, and the position of the observer. This sketch allows building a geometric model for the theoretical image formed between the two layers. The second model for the intensity of the observed radiation assumes a synchrotron source for luminosity proportional to the flux of kinetic energy, $L_m$,

$L_m=\frac{1}{2}\rho A{V}^3\mathrm{,}$ (36)

where A is the considered area, V the velocity and $\rho$ the density, see Formula (A28) in [17]. In our case $A=R^2\Delta \Omega$, where $\Delta \Omega$ is the considered solid angle along the chosen direction. This means

$L_m=\frac{1}{2}\rho \Delta \Omega R^2{V}^3\mathrm{,}$ (37)

where R is the instantaneous radius of the SNR and $\rho$ is the density in the advancing layer in which the synchrotron emission takes place. The observed luminosity along a given direction can be expressed as

$L={\epsilon }_m{L}_m\mathrm{,}$ (38)

Figure 20. First quadrant for the layer between the upper (red curve) and the lower (green curve) section of SN 1987A in the x-z plane with an inverse square profile when $0<\theta <\pi$. Parameters as in Figure 10. The observer is at infinity of the x-axis and two lines of sight, s1 and s2, are marked.

where ${\epsilon }_m$ is a constant of conversion from the mechanical luminosity to the observed luminosity in the considered band. The numerical algorithm which allows us to build a complex image is now outlined.

· An empty (value = 0) memory grid $M\left(i\mathrm{,}j\mathrm{,}k\right)$ which contains 400³ pixels is considered.

· We first generate an internal 3D surface by rotating the ideal image of 180˚ around the polar direction and a second external surface at a fixed distance $\Delta R$ from the first surface. As an example, we fixed $\Delta R=R/12$, where R is the momentary radius of expansion. The points on the memory grid which lie between the internal and external surfaces are stored in memory on $M\left(i\mathrm{,}j\mathrm{,}k\right)$ with a variable integer number according to Formula (43).

· Each point of $M\left(i\mathrm{,}j\mathrm{,}k\right)$ has spatial coordinates $x\mathrm{,}y\mathrm{,}z$ which can be repre-sented by the following $1\times 3$ matrix, A,

$A=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]$ (39)

The orientation of the object is characterized by the Euler angles $\left(\Phi \mathrm{,}\Theta \mathrm{,}\Psi \right)$ and therefore by a total $3\times 3$ rotation matrix, E, see [18]. The matrix point is represented by the following $1\times 3$ matrix, B,

$B=E\cdot A$ (40)

· The intensity map is obtained by summing the points of the rotated images along a particular direction.

· The effect of the insertion of a threshold intensity, $I_{tr}$, given by the observational techniques, is now analysed. In both models the threshold intensity can be parametrized by $I_{max}$, the maximum value of intensity characterizing the map.

More details on the theory of the image can be found in [16].

4.2. Results for the First Model

The radiation of the SN is supposed to originate from between a lower boundary of radius $r_{inf}$, given as an example by $r\left(t\mathrm{<mo>\; ;\; </mo>}{z}_0\mathrm{,}{r}_0\mathrm{,}\theta \right)$ of Equation (8), and an upper boundary with radius given by the equation

$r_{sup}=r_{inf}+\frac{r_{inf}}{12}\mathrm{.}$ (41)

The above layer is visualized in Figure 21 for SN 1987A when an inverse square profile is adopted.

The image of SN 1987A is formed between the two layers. The data of the two layers are stored in memory as Cartesian coordinates

$x_{sup}\left(n\right)=r_{sup}\left(n\right)sin\left(\theta \left(n\right)\right)$ (42a)

$z_{sup}\left(n\right)=r_{sup}\left(n\right)cos\left(\theta \left(n\right)\right)$ (42b)

$x_{inf}\left(n\right)=r_{inf}\left(n\right)sin\left(\theta \left(n\right)\right)$ (42c)

$z_{inf}\left(n\right)=r_{inf}\left(n\right)cos\left(\theta \left(n\right)\right)$ (42d)

where n is the considered index. In order to build the first model we assume that z varies between $z_{sup}\left(1\right)$ and $z_{sup}\left(N\right)$ where N represents the number of points. The image is assumed to be proportional to the distance between the x coordinate on the upper layer and the x coordinate on the lower layer as a function of the position z on the polar axis, see Figure 22. In principle, astronomers should produce these observational cuts in order to compare the observations with the theory here developed.

The same algorithm allows building a map of the theoretical intensity for the rotated image of SN 1987A, see Figure 23. The observed image of SN 1987A should be digitized in order to make a comparison with the theory.

4.3. Results for the Second Model

This second model allows simulating particular effects, such as the triple ring system of SN 1987A, see Figure 24, where three zones or holes with theoretical intensity under the threshold value are visible. The enigma of the three holes is solved. Characteristic features, such as the “jet appearance” visible in some maps for SN 1006, see the X-ray map 3 and γ-ray map 3, are theoretically modeled in Figure 25. The enigma of the jet, which is a cone, in a nearly spherical expansion is solved.

Figure 21. The layer between the upper (green points) and the lower (red points) section of SN 1987A in the x-z plane with an inverse square profile when $0<\theta <\pi$. Parameters as in Figure 10.

Figure 22. Cut of the mathematical intensity I, as a function of the polar z-axis for SN 1987A with an inverse square profile. Parameters as in Figure 10 and $N=240$.

Figure 23. Map of the theoretical intensity I, for SN 1987A with an inverse square profile made by 400 × 400 pixels. Parameters as in Figure 10.

Figure 24. Model map of SN 1987A rotated in accordance with the observations, for an inverse square medium with parameters as in Figure 10. The three Euler angles characterizing the orientation are $\Phi ={105}^{\circ }$, $\Theta ={45}^{\circ }$ and $\Psi =-{165}^{\circ }$. In this map $I_{tr}=I_{max}/2.2$.

Figure 25. Model map of SN 1006 rotated in accordance with the γ-ray observations with HESS, for a toroidal medium. Physical parameters as in Figure 19. The three Euler angles characterizing the orientation of the observer are $\Phi ={180}^{\circ }$, $\Theta ={90}^{\circ }$ and $\Psi =-{55}^{\circ }$. In this map $I_{tr}=I_{max}/1.2$.

Table 1. Synoptic percentage reliability for the best model of SNRs with different profiles for the density.

5. Conclusions

Models

We have adopted four models in the framework of the conservation of energy for the thin layer approximation and applied them to two SNRs; the results for the percentage reliability, see Formula (35) are shown in Table 1. The best fit for the asymmetric section of SN 1987A is represented by the exponential model for decreasing density in the polar direction and by the toroidal profile of the density for SN 1006.

Image theory

Two models for the intensity of the image in an SNR have been presented, which both require an analytical or numerical function for the advancing radius in terms of the polar angle. The appearance of the triple ring in SN 1987A, see Figure 24, and the jet-feature in SN 1006, see Figure 25, have been simulated.

Back reaction

The derivation of an approximate form for the losses, see Equation (10), allows deriving a finite rather than infinite length for an SNR, see Equation (15).

Acknowledgements

At the time of writing, an animated version of Figure 11 is visible at http://personalpages.to.infn.it/~zaninett/image/sn1987a_animation.gif.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Borkowski, K.J., Gwynne, P., Reynolds, S.P., et al. (2017) Asymmetric Expansion of the Youngest Galactic Supernova Remnant G1.9+0.3. The Astrophysical Journal Letters, 837, L7. https://doi.org/10.3847/2041-8213/aa618c
[2]	Reynolds, S.P. (2016) Asymmetric Expansion of the Youngest Galactic Supernova Remnant G1.9+0.3. Supernova Remnants: An Odyssey in Space after Stellar Death, Chania, 6-11 June 2016, 147-157.
[3]	Stafford, J. and Lopez, L.A. (2017) Measuring the Symmetry of Supernova Remnants in the Radio. American Astronomical Society Meeting Abstracts #229, 148-155. https://ui.adsabs.harvard.edu/abs/2017AAS...22914807S/abstract
[4]	Collier, A., Bachrach, H., Fryer, C. and Ellinger, C. (2017) Asymmetry in Supernovae. American Astronomical Society Meeting Abstracts #229, 434-438. https://ui.adsabs.harvard.edu/abs/2017AAS...22943408C/abstract
[5]	Moranchel Basurto, A., Velázquez, P.F., Schneiter, E.M. and Esquivel, A. (2019) Asymmetries in the Emission from Young Supernova Remnants: The Case of Tycho. Supernova Remnants: An Odyssey in Space after Stellar Death II, Chania, 3-8 June 2019, 190-200.
[6]	Andrews, J.E., Sand, D.J., Valenti, S., et al. (2019) SN 2017gmr: An Energetic Type II-P Supernova with Asymmetries. The Astrophysical Journal, 885, 43.
[7]	Livneh, R. and Katz, B. (2019) An Asymmetric Explosion Mechanism May Explain the Diversity of Si II Line Widths in Type Ia Supernovae. Monthly Notices of the Royal Astronomical Society, 494, 5811-5824. https://doi.org/10.1093/mnras/staa974
[8]	Strasburger, E., Siebert, M. and Foley, R. (2020) Constraining the Asymmetry of Type Ia Supernovae Using Late-Time Spectra. Bulletin of the American Astronomical Society, 52, 276-277.
[9]	Zaninetti, L. (2020) Energy Conservation in the Thin Layer Approximation: I. The Spherical Classic Case for Supernovae Remnants. International Journal of Astronomy and Astrophysics, 10, 71-88. https://doi.org/10.4236/ijaa.2020.102006
[10]	Racusin, J.L., Park, S., Zhekov, S., Burrows, D.N., Garmire, G.P. and McCray, R. (2009) X-Ray Evolution of SNR 1987A: The Radial Expansion. The Astrophysical Journal, 703, 1752. https://doi.org/10.1088/0004-637X/703/2/1752
[11]	McCray, R. (2017) The Physics of Supernova 1987A. In: Alsabti, A.W. and Murdin, P., Eds., Handbook of Supernovae, Springer International Publishing, Cham, 2181-2210. https://doi.org/10.1007/978-3-319-21846-5_96
[12]	France, K., McCray, R., Fransson, C., Larsson, J., Frank, K.A., Burrows, D.N., Challis, P., Kirshner, R.P., Chevalier, R.A., Garnavich, P., Heng, K., Lawrence, S.S., Lundqvist, P., Smith, N. and Sonneborn, G. (2015) Mapping High-Velocity Hα and Lyα Emission from Supernova 1987A. The Astrophysical Journal Letters, 801, L16. https://doi.org/10.1088/2041-8205/801/1/L16
[13]	Uchida, H., Yamaguchi, H. and Koyama, K. (2013) Asymmetric Ejecta Distribution in SN 1006. The Astrophysical Journal, 771, 56. https://doi.org/10.1088/0004-637X/771/1/56
[14]	Katsuda, S. (2017) Supernova of 1006 (G327.6+14.6). In: Alsabti, A.W. and Murdin, P., Eds., Handbook of Supernovae, Springer International Publishing, Cham, 63-81.
[15]	Reynolds, S.P. and Gilmore, D.M. (1986) Radio Observations of the Remnant of the Supernova of A.D. 1006. I. Total Intensity Observations. Astronomical Journal, 92, 1138-1144. https://doi.org/10.1086/114244
[16]	Zaninetti, L. (2018) The Physics of Asymmetric Supernovae and Supernovae Remnants. Advances in Applied Physics, 6, 25-73. https://doi.org/10.12988/aap.2018.897
[17]	De Young, D.S. (2002) The Physics of Extragalactic Radio Sources. University of Chicago Press, Chicago.
[18]	Goldstein, H., Poole, C. and Safko, J. (2002) Classical Mechanics. Addison-Wesley, San Francisco.