**International Journal of Astronomy and Astrophysics** Vol.3 No.2(2013), Article ID:32351,6 pages DOI:10.4236/ijaa.2013.32010

Three Dimensional Evolution of SN 1987A in a Self-Gravitating Disk

Dipartimento di Fisica, via P. Giuria 1, Turin, Italy

Email: zaninetti@ph.unito.it

Copyright © 2013 L. Zaninetti. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received February 16, 2013; revised March 19, 2013; accepted March 28, 2013

**Keywords:** Supernovae; General Supernovae; Individual (SN 1987A) ISM; Supernova Remnants

ABSTRACT

The introduction of an exponential or power law gradient in the interstellar medium (ISM) allows to produce an asymmetric evolution of the supernova remnant (SNR) when the framework of the thin layer approximation is adopted. Unfortunately both the exponential and power law gradients for the ISM do not have a well defined physical meaning. The physics conversely is well represented by an isothermal self-gravitating disk of particles whose velocity is everywhere Maxwellian. We derived a law of motion in the framework of the thin layer approximation with a control parameter of the swept mass. The photon’s losses, which are often neglected in the thin layer approximation, are modeled trough velocity dependence. The developed framework is applied to SNR 1987A and the three observed rings are simulated.

1. Introduction

The expansion of supernova remnant (SNR) is usually explained by the theory of the shock waves, as an example [1,2], which analyzes self-similar solutions with varying inverse power law exponents for the density profile of the advancing matter, , and ambient medium,. The SNRs conversely often present a weak departure from the circular symmetry, as an example, in SN 1006 1006 a ratio of 1.2 between maximum and minimum radius which has been measured, see [3] or great asymmetry as the triple ring of SN 1987A 1987a, see [4,5].

The asymmetries of SNRs are hardly explained by the shock theory which is usually developed in spherical or planar symmetry. The presence of gradients in the density of the medium allows to explain such asymmetries when the conservation of the momentum is adopted. As an example an exponential decay allows to model SN1006 1006 and SN 1987A 1987a, see [6]. At the same time the various gradients that can be adopted such as Gaussian, power law and exponential do not have a precise physical basis. The auto-gravitating Maxwellian medium provides a gradient given by the secant law. In this paper we will answer the following questions.

• Is it possible to deduce an analytical formula for the temporal evolution of an SN in the presence of a vertical profile density as given by an ISD?

• Is it possible to deduce an equation of motion for SN with an adjustable parameter that can be found from a numerical analysis of the diameter—time relationship?

• Is it possible to simulate the appearance of the triple ring system of SN 1987A 1987a?

2. A Modified Thin Layer Approximation

The thin layer approximation with non cubic dependence (NCD), , in classical physics assumes that only a fraction of the total mass enclosed in the volume of the expansion accumulates in a thin shell just after the shock front. The global mass between 0 and is

where is the density of the ambient medium. The swept mass included in the thin layer which characterizes the expansion is

(1)

The mass swept between 0 and is

(2)

The conservation of radial momentum requires that, after the initial radius,

(3)

where and are the radius and velocity of the advancing shock. The velocity is

(4)

The velocity as a function of the radius when the expansion starts at the origin of the coordinates is:

(5)

The velocity can be expressed as

(6)

and the integration between, and, gives

(7)

We solve the previous equation for and we obtain

(8)

The velocity is:

(9)

where

Equation (8) can also be solved with a solution of type, being a constant, and the asymptotic result is:

(10)

The asymptotic solution for the velocity is

(11)

3. Asymmetrical Law of Motion with NCD

Given the Cartesian coordinate system, the plane will be called the equatorial plane, , where is the latitude angle which has range, and is the distance from the origin. In our framework, the polar angle of the spherical coordinate system is. Here we adopt the density profile of a thin self-gravitating disk of gas which is characterized by a Maxwellian distribution in velocity and distribution which varies only in the z-direction (ISD). The number density distribution is

(12)

where is the density at, is a scaling parameter and is the hyperbolic secant ([7,8]). On assuming that the expansion starts at, we can write and therefore

(13)

where is the radius of the advancing shell.

The 3D expansion that starts at the origin of the coordinates will be characterized by the following properties

• Dependence of the momentary radius of the shell on the polar angle that has a range.

• Independence of the momentary radius of the shell from, the azimuthal angle in the x-y plane, that has a range.

The mass swept, , along the solid angle, between 0 and is

(14)

where

(15)

where is the initial radius and the mass of the hydrogen. The integral is

(16)

(17)

where

(18)

is the dilogarithm, see [9-11].

The conservation of the momentum gives

(19)

where is the velocity at and is the initial velocity at where is the NCD parameter. This means that only a fraction of the total mass enclosed in the volume of the expansion accumulates in a thin shell just after the shock front. The emission of photons produces losses in the momentum carried by the advancing shell here assumed to be proportional to. The generalized conservation law is now

(20)

where is a constant. At the moment of writing is not possible to deduce the parameter which controls the quantity of momentum carried away by the photons due to the presence of the parameter which controls the quantity of swept mass during the expansion. According to the previous generalized Equation (20) an analytical expression for can be found but it is complex and therefore we omit it. In this differential equation of the first order in the variable can be separated and the integration term by term gives

(21)

where is the time and the time at. We therefore have an equation of the type

(22)

where has an analytical but complex form. In our case in order to find the root of the nonlinear Equation (22) the FORTRAN SUBROUTINE ZRIDDR from [12] has been used.

The physical units have not yet been specified, for length and for time are perhaps an acceptable astrophysical choice. On adopting the previous units, the initial velocity is expressed in pc/yr and should be converted into; this means that V_{0} = 1.020738 × 10^{−6}V_{1} where is the initial velocity expressed in.

4. Application to SN 1987A

4.1. The Simulation

The SN 1987A exploded in the Large Magellanic Cloud in 1987. The observed image is complex and we will follow the nomenclature of [13] which distinguish between torus only, torus +2 lobes and torus +4 lobes. In particular we concentrate on the torus which is characterized by a distance from the center of the tube and the radius of the tube. Table 2 of [13] reports the relationship between distance of the torus in arcsec and time since the explosion in days. Figure 3 in [14] reports diameter in pc versus year when the counting pixels is adopted. Table 1 reports the data used in the simulation, Figures 1 and 2 report the diameter and the velocity in the equatorial plane respectively as well the observational data as found by the counting pixels method, see [14]. Figure 3 reports the asymmetric expansion in a section crossing the center, Figures 4 and 5 report the radius and the velocity respectively as a function of the position angle θ.

The swept mass in our model is not constant but varies with the value of latitude, see Figure 6.

4.2. The Image

We assume that the emissivity is proportional to the number density.

(23)

where is the emission coefficient and is a constant function. This can be the case of synchrotron radiation in presence of an isotropic distribution of electrons with a power law distribution in energy, ,

Table 1. Numerical value of the parameters of the simulation for SN 1987A.

Figure 1. Diameter as a function of time for a self-gravitating medium (full line) when θ = 0.001 (equatorial plane) and astronomical data of torus only as extracted from Figure 3 in [14] (dashed line). Physical parameters as in Table 1.

Figure 2. Velocity as a function of time for a self-gravitating medium (full line) when θ = 0.001 (equatorial plane) and astronomical data of torus only as extracted from Figure 3 in [14] (empty stars). Physical parameters as in Table 1.

Figure 3. Section of SN 1987A in the x-z plane. The horizontal and vertical axis are in pc. Physical parameters as in Table 1.

Figure 4. Radius in pc of SN 1987A as a function of the position angle in degrees for a self-gravitating medium. Physical parameters as in Table 1.

Figure 5. Velocity in km/s of SN 1987A as a function of the position angle in degrees for a self-gravitating medium. Physical parameters as in Table 1.

Figure 6. Swept mass of SN 1987A in arbitrary units as a function of the position angle in degrees, for a self-gravitating medium. Physical parameters as in Table 1.

(24)

where is a constant. The source of synchrotron luminosity is assumed here to be the flux of kinetic energy, ,

(25)

where is the considered area, see formula (A28) in [15]. In our case, which means

(26)

where is the instantaneous radius of the SNR and is the density in the advancing layer in which the synchrotron emission takes place. The total observed luminosity can be expressed as

(27)

where is a constant of conversion from the mechanical luminosity to the total observed luminosity in synchrotron emission.

The numerical algorithm which allows us to build a complex image is now outlined.

• An empty (value = 0) memory grid which contains pixels is considered

• We first generate an internal 3D surface by rotating the ideal image around the polar direction and a second external surface at a fixed distance from the first surface. As an example, we fixed, where is the momentary radius of expansion. The points on the memory grid which lie between the internal and external surfaces are memorized on with a variable integer number according to formula (26) and density proportional to the swept mass, see Figure 6.

• Each point of has spatial coordinates which can be represented by the following matrix, ,

(28)

The orientation of the object is characterized by the Euler angles and therefore by a total rotation matrix, , see [16]. The matrix point is represented by the following matrix, ,

(29)

• 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, , given by the observational techniques, is now analyzed. The threshold intensity can be parametrized to, the maximum value of intensity characterizing the map.

The image of SN 1987A having polar axis along the z-direction, is shown in Figure 7. The image for a realistically rotated SN 1987A is shown in Figure 8 and a comparison should be done with the observed triple ring system as reported in Figure 1 of [5] where the outer lobes are explained by light-echo effects.

Two cuts of intensity are reported in Figure 9, and is interesting to note that the polar cut is characterized by four spikes. Here conversely the two outer lobes are explained by the image theory applied to a thin advancing shell.

5. Conclusions

5.1. Law of Motion

The presence of a thin self-gravitating disk of gas which is characterized by a Maxwellian distribution in velocity

Figure 7. Map of the theoretical intensity of SN 1987A in the presence of a self-gravitating medium. Physical parameters as in Table 1. The three Euler angles characterizing the orientation are Φ = 180˚, Θ = 90˚ and Ψ = 0˚. This combination of Euler angles corresponds to the rotated image with the polar axis along the z-axis.

Figure 8. Model map of SN 1987A 1987a rotated in accordance with the observations, for a self-gravitating medium. Physical parameters as in Table 1. the orientation of the observer are Φ = 105˚, Θ = 55˚ and Ψ = −165˚. This combination of Euler angles corresponds to the observed image. In this map.

Figure 9. Two cuts of the model intensity across the center of the realistically rotated SN 1987A: equatorial cut (full line) and polar cut (dotted line).

and distribution of density which varies only in the zdirection produces an axial symmetry in the temporal evolution of a SN. The resulting Equation (22) has a form which can be solved numerically. The results depend from the chosen latitude and more precisely the radius as function of the time is maximum at the poles and minimum at the equator. A comparison between the observed radius and the velocity at the equator, see [14], gives an agreement of 92% for the theoretical versus observed radius and of 62% for the theoretical versus observed velocity. This is the first result that confirms the rationality of the method here presented. The momentum carried away by the photons can be modeled by introducing a reasonable form for the photon’s losses, see Equation (20).

5.2. Images

The combination of different processes such as the bipolar 3D shape , the thin layer approximation (which means advancing layer having thickness of the momentary radius), the conversion of the rate of kinetic energy into luminosity and the observer’s point of view allows to simulate the observed triple ring morphology of SN 1987A 1987a. The framework here adopted is similar to that adopted to explain the asymmetric shape of the nebula around -Carinae (Homunculus), see [17,18]. We therefore suggest that the astronomers should map the non rotated radius and velocity of SN 1987A 1987a following the outlines of Table 1 in [19]. These maps, if realized, will produce a clearer framework for the observed triple ring. The suggested similarity of the triple ring system of SN 1987A 1987a with -Carinae is the second result which ensures the rationality of the method here presented.

REFERENCES

- R. A. Chevalier, “Self-Similar Solutions for the Interaction of Stellar Ejecta with an External Medium ApJ 1982258, 790,” 1982.
- R. A. Chevalier, “The Radio and X-Ray Emission from Type II Supernovae ApJ 259, 302,” 1982.
- S. P. Reynolds and D. M. Gilmore, “Radio Observations of the Remnant of the Supernova of A.D. 1006. I—Total Intensity Observations AJ 92, 1138,” 1986.
- T. Morris and P. Podsiadlowski, “The Triple-Ring Nebula around SN 1987A: Fingerprint of a Binary Merger Science 315, 1103 (Preprint arXiv:astro-ph/0703317) Three Dimensional Evolution of SN 1987A in a Self-Gravitating Disk,” 2007.
- A. Tziamtzis, P. Lundqvist, P. Gröningsson and S. Nasoudi-Shoar, “The Outer Rings of SN 1987A A&A 527 A35 (Preprint 1008.3387),” 2011.
- L. Zaninetti, “On the Spherical-Axial Transition in Supernova Remnants Astrophysics and Space Science 337, 581 (Preprint 1109.4012),” 2012.
- G. Bertin, “Dynamics of Galaxies,” Cambridge University Press, Cambridge, 2000.
- P. Padmanabhan, “Theoretical Astrophysics. Vol. III: Galaxies and Cosmology,” Cambridge University Press, Cambridge, 2002.
- C. J. Hill, “Ueber Die Integration Logarithmisch-Rationaler Differentiale,” Journal für die Reine und Angewandte Mathematik, Vol. 1828, No. 3, 2009, pp. 101-159. doi:10.1515/crll.1828.3.101
- L. Lewin, “Polylogarithms and Associated Functions,” New York, North Holland, 1981.
- F. W. J. Olver, D. W. Lozier, R. F. Boisvert and C. W. Clark, “NIST Handbook of Mathematical Functions,” Cambridge University Press, Cambridge, 2010.
- W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, “Numerical Recipes in FORTRAN. The Art of Scientific Computing,” Cambridge University Press, Cambridge, 1992.
- J. L. Racusin, S. Park, S. Zhekov, D. N. Burrows, G. P. Garmire and R. McCray, “X-Ray Evolution of SNR 1987A: The Radial Expansion ApJ 703, 1752 (Preprint 0908. 2097),” 2009.
- B. T. Chiad, L. M. Karim and L. T. Ali, “Study the Radial Expansion of SN 1987A Using Counting Pixels Method,” International Journal of Astronomy and Astrophysics, Vol. 2, 2012, pp. 199-203.
- D. S. De Young, “The Physics of Extragalactic Radio Sources,” University of Chicago Press, Chicago, 2002.
- H. Goldstein, C. Poole and J. Safko, “Classical Mechanics,” Addison-Wesley, San Francisco, 2002.
- L. Zaninetti, “On the 3D Structure of the Nebula around Carinae Advances in Space Research 46, 1341 (Preprint 1010.2364),” 2010.
- L. Zaninetti, “Interaction of Planetary Nebulae, Eta-Carinae and Supernova Remnants with the Interstellar Medium,” In: B. Cancellieri and V Mamedov, Eds., Interstellar Medium: New Research, Nova Publishers, New York, 2012, pp. 91-146.
- N. Smith, “The Structure of the Homunculus. I. Shape and Latitude Dependence from H2 and Fe II Velocity Maps of eta Carinae ApJ 644, 1151 (Preprint arXiv:astro-ph/0602464),” 2006.