Energy Conservation in the thin layer approximation: II. The asymmetric classic case for supernovae remnants

Here we extend the conservation of energy in the framework of the thin layer approximation to the asymmetrical case. Four types of interstellar medium 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.


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 were 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. down symmetry in the polar angle, r(θ) = r(π − θ), and by a polar axis, x = 0, y = 0, which produces a right and left symmetry for the azimuthal angle r(ϕ) = r(ϕ + π). The basic symmetries are now outlined In other words, the numerical simulations will be developed in the first quadrant (θ ∈ [0, π 2 ], ϕ = 0 and then applied to the other three quadrants.
The conservation of kinetic energy in spherical coordinates along the solid angle ∆Ω in the framework of the thin layer approximation is where M 0 (r 0 ; θ) and M (r; θ) 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 The velocity as a function of the radius is In the following, r 0 is the radius after which the density starts to decrease, r is the radius of expansion in spherical coordinates, ρ c is the density at r = 0, z = r cos(θ) is the Cartesian coordinate along the Z axis and cos(θ) 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.

An inverse square profile of the density
We assume that the medium around the SN scales with the axial piecewise dependence The mass M 0 swept in the interval [0,r 0 ] is The total mass M (r; r 0 , ρ c ) swept in the interval [0,r] is The positive solution of equation (2) gives the velocity as a function of the radius: v(r; r 0 , z 0 ) = r 0 3 (cos (θ)) 2 + 3 z 0 2 r − 3 z 0 2 r 0 r 0 cos (θ) v0 r 0 r 0 3 (cos (θ)) 2 + 3 z 0 2 r − 3 z 0 2 r 0 .
The differential equation which models the energy conservation is The analytical solution of the above differential equation is The above inverse square profile for density satisfies the basic symmetries as outlined in equations (1a) and (1b).

Back reaction for the inverse square profile
The radiative losses per unit length are assumed to scale as where is a constant and rho s is the density in the thin advancing layer which is 4 ρ. Inserting into the above equation the velocity to first order as given by equation (6), the radiative losses Q(r; r 0 , v 0 , z 0 , , θ) are The sum of the radiative losses between r 0 and r is given by the integral The conservation of energy in the presence of the back reaction due to the radiative losses is The analytical solution for the velocity to second order, v c (r; where The inclusion of the back reaction allows the evaluation of the SRS's maximum length, r back (r 0 , z 0 , θ, ), which can be derived by setting this velocity to zero: Figure 1 shows the SRS's maximum length as a function of the constant of conversion and the polar angle θ. The above figure shows that the SNR cannot reach an infinite length at an infinite time but has a finite lifetime and length.

A power law profile of the density
The medium is assumed to scale as where α is a positive real number. The total mass M (r; .
The differential equation which governs the motion is The above differential equation does not have an analytical solution; a Taylor expansion of order 3 covers a limited range in time The above power law profile for the density satisfies the basic symmetries as outlined in equations (1a) and (1b).

An exponential profile of the density
The medium is assumed to scale as In the absence of an analytical solution, we present the Taylor expansion of order 3 for the trajectory The above exponential profile for the density satisfies the basic symmetries as outlined in equations (1a) and (1b).

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 We now consider its intersection with the plane y = 0 and as a consequence the equation of the torus is now We now evaluate the intersection between the torus and a straight line of equation which crosses the centre, (z = 0, The above angle allows us to define the following three zones. 1. θ < θ crit : the line does not intersect the torus. 2. θ = θ crit : the line is tangent to the torus at one degenerate single point, with Cartesian coordinates (x crit , z crit ), see Figure 3. 3. θ > θ crit : the line intersects the torus at two points, with Cartesian coordinates (x 1 , z 1 ) and (x 2 , z 2 ), see Figure 4. The Cartesian coordinates at θ = θ crit are The Cartesian coordinates of the two intersections when θ > θ crit are We now assume that the density of matter is ρ c outside the torus and ρ 1 inside the torus. The mass swept in the third zone requires a careful analysis. When θ > θ crit , the mass swept along a line before the intersection with the torus, M I (r; ρ c ), is 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 (r; ρ c , ρ 1 , R T , r T ), is The mass swept when r is outside the torus in the third zone is The equation of motion is solved through the Euler method when the following recursive equations are solved where r n , v n , M n are the temporary radius, velocity, and total mass, respectively, ∆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.

Astrophysical applications
The complex structure of SN 1987A , see ithe 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, θ = π 2 . A useful resource for calibration is the geometric section of SN 1987A which is given as a sketch in Figure 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.45.0 keV band of Fig. 1 in [13]. The following observed radii can be extracted: R up = 11.69 pc in the polar direction and R eq = 8.7 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, obs , where r num is the theoretical radius of the SNR, r obs is the observed radius of the SNR, and the index j varies from 1 to the number of available observations.

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.

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 presentes 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.

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 .
The above results shows that an exponential profile for the density is comparable with the observed sections of the two SNRS here simulated.

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

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 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 < θ < π. Parameters as in Figure 10. The observer is at infinity of the x-axis and two lines of sight, s1 and s2, are marked. 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 , where A is the considered area, V the velocity and ρ the density, see formula (A28) in [17]. In our case A = R 2 ∆Ω, where ∆Ω is the considered solid angle along the chosen direction. This 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 observed luminosity along a given direction can be expressed as where 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.
-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 ∆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 stored in memory on M(i, j, k) with a variable integer number according to formula (37). -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 [18]. 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 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].

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(t; z0, r 0 , θ) of equation (8), and an upper boundary with radius given by the equation 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 (n) = r sup (n) sin(θ(n)) (42a) z sup (n) = r sup (n) cos(θ(n)) (42b) x inf (n) = r inf (n) sin(θ(n)) (42c) where n is the considered index. In order to build the first model we assume that z varies between z sup (1) and z sup (N ) 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. 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 < θ < π. Parameters as in Figure 10.

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 7 and γ-ray map 8, are theoretically modeled in Figure 25. The enigma of the jet, which is a cone, in a nearly spherical expansion is solved. 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 Φ=180 • , Θ=90 • and Ψ =-55 • . In this map Itr = Imax/1.2.

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).