Three dimensional evolution of SN 1987A in a self-gravitating disk

The introduction of a 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 ,that are often neglected in the thin layer approximation, are modeled trough a velocity dependence. The developed framework is applied to SNR 1987A and the three observed rings are simulated.


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

A modified thin layer approximation
The thin layer approximation with non cubic dependence (NCD) , p , 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 R 0 is 4 3 πρR 3 0 where ρ is the density of the ambient medium. The swept mass included in the thin layer which characterizes the expansion is The mass swept between 0 and R is The conservation of radial momentum requires that, after the initial radius R 0 , where R and V are the radius and velocity of the advancing shock. The velocity is The velocity as a function of the radius when the expansion starts at the origin of the coordinates is: The velocity can be expressed as and the integration between R , t and R 0 , t 0 gives We solve the previous equation for R and we obtain The velocity is: Equation (8) can also be solved with a solution of type R = K(t − t 0 ) α , k being a constant, and the asymptotic result is: The asymptotic solution for the velocity is

Asymmetrical law of motion with NCD
Given the Cartesian coordinate system (x, y, z), the plane z = 0 will be called the equatorial plane, z = R sin(θ), where θ is the latitude angle which has range [−90 • ↔ +90 • ], and R is the distance from the origin. In our framework, the polar angle of the spherical coordinate system is 90 − θ. 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 where n 0 is the density at z = 0, h is a scaling parameter and sech is the hyperbolic secant ( [7,8]). On assuming that the expansion starts at z ≈ 0 , we can write z = R sin(θ) and therefore where R 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 [−90 • ↔ +90 • ].
• Independence of the momentary radius of the shell from φ , the azimuthal angle in the x-y plane, that has a range [0 The mass swept, M , along the solid angle ∆ Ω, between 0 and R is where where R 0 is the initial radius and m H the mass of the hydrogen. The integral is where is the dilogarithm, see [9,10,11]. The conservation of the momentum gives whereṘ is the velocity at R andṘ 0 is the initial velocity at R = R 0 where p 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 R 2Ṙ . The generalized conservation law is now where C l is a constant. At the moment of writing is not possible to deduce the parameter C l which controls the quantity of momentum carried away by the photons due to the presence of the parameter p which controls the quantity of swept mass during the expansion. According to the previous generalized equation (20) an analytical expression for V (R) can be found but it is complex and therefore we omit it. In this differential equation of the first order in R the variable can be separated and the integration term by term gives where t is the time and t 0 the time at R 0 . We therefore have an equation of the type where F(R, R 0 , h) N L 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, pc for length and yr for time are perhaps an acceptable astrophysical choice. On adopting the previous units, the initial velocity V 0 is expressed in pc yr and should be converted into km s ; this means that V 0 = 1.020738 10 −6 V 1 where V 1 is the initial velocity expressed in km s .

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]. Fig. 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.

The image
We assume that the emissivity is proportional to the number density.
where j ν is the emission coefficient and K 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, N (E),  Table 1.
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. where K s is a constant. The source of synchrotron luminosity is assumed here to be the flux of kinetic energy, L m , where A is the considered area, see formula (A28) in [15]. In our case A = 4πR 2 , which means where R 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 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.  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.
• An empty (value=0) memory grid M(i, j, k) which contains N DIM 3 pixels is considered • We first generate an internal 3D surface by rotating the ideal image 180 • around the polar direction and a second external surface at a fixed distance ∆R from the first surface. As an example, we fixed ∆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 memorized on M(i, j, k) with a variable integer number according to formula (26) and density ρ proportional to the swept mass, see Fig.  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.
• Each point of M(i, j, k) has spatial coordinates x, y, z which can be represented by the following 1 × 3 matrix, A, The orientation of the object is characterized by the Euler angles (Φ, Θ, Ψ) and therefore by a total 3 × 3 rotation matrix, E, see [16]. The matrix point is represented by the following 1 × 3 matrix, B, • 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 analyzed. The threshold intensity can be parametrized to I max , the maximum value of intensity characterizing the map.
The image of SN 1987A having polar axis along the z-direction, is shown in Fig. 7. The image for a realistically rotated SN 1987A is shown in Fig. 8 and a   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 Itr = Imax/4 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.

Conclusions
Law of motion The presence of a thin self-gravitating disk of gas which is characterized by a Maxwellian distribution in velocity and distribution of density which varies only in the z-direction produces an axial symmetry in the temporal evolution of a SN. The resulting Eq. (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 eqn. ( 20).

Images
The combination of different processes such as the bipolar 3D shape , the thin layer approximation (which means advancing layer having thickness ≈ 1/10 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 . 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 following the outlines of Table 1 in [19]. These maps , if realized , will produce a more clear framework for the observed triple ring. The suggested similarity of the triple ring system of SN 1987A with η-Carinae is the second result which ensures the rationality of the method here presented.