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  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 . The symmetry of many SNRs in 6-cm and 20-cm Very Large Array images was measured by . A three dimensional smoothed particle hydrodynamic simulation was used by  in order to explain the observed asymmetries in SNs.  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 . The differences in the Si II line widths in type Ia supernovae with an asymmetric CSM were explained by . An asymmetric explosion as seen from different viewing angles was suggested by  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  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 , the polar angle , and the azimuthal angle . In the following, the profiles of the density considered are independent of the azimuthal angle and are functions of the distance z from the 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, , which produces an up and down symmetry in the polar angle, , and by a polar axis, , which produces a right and left symmetry for the azimuthal angle . The basic symmetries are now outlined
In other words, the numerical simulations will be developed in the first quadrant ( , 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 and are the swept masses at and r, while and v are the velocities of the thin layer at 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, is the radius after which the density starts to decrease, r is the radius of expansion in spherical coordinates, is the density at , is the Cartesian coordinate along the Z axis and 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
The mass swept in the interval is
The total mass swept in the interval [0, r] is
The positive solution of Equation (2) gives the velocity as a function of the radius:
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).
2.2. Back Reaction for the Inverse Square Profile
The radiative losses per unit length are assumed to scale as
where is a constant and is the density in the thin advancing layer which is . Inserting into the above equation the velocity to first order as given by Equation (6), the radiative losses are
The sum of the radiative losses between 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, , is
The inclusion of the back reaction allows the evaluation of the SRS’s maximum length, , 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.
Figure 1. A 2D map of in pc. The parameters are and .
2.3. A Power Law Profile of the Density
The medium is assumed to scale as
where is a positive real number. The total mass swept in the interval is
The velocity as a function of the radius is
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).
2.4. An Exponential Profile of the Density
The medium is assumed to scale as
The total mass swept in the interval [0, r] is
The velocity as a function of the radius is
The differential equation which governs the motion is
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).
2.5. A Toroidal Profile of the Density
A torus is swept out by revolving a small circle of radius about an axis lying in the same plane as the circle but outside it, say at a distance of , 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 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, ( ),
where 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, :
Figure 2. A 3D view of a torus.
The above angle allows us to define the following three zones.
1) : the line does not intersect the torus.
2) : the line is tangent to the torus at one degenerate single point, with Cartesian coordinates ( ), see Figure 3.
3) : the line intersects the torus at two points, with Cartesian coordinates ( ) and ( ), see Figure 4.
The Cartesian coordinates at are
The Cartesian coordinates of the two intersections when are
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).
We now assume that the density of matter is outside the torus and inside the torus. The mass swept in the third zone requires a careful analysis. When , the mass swept along a line before the intersection with the torus, , 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, , 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 , , are the temporary radius, velocity, and total mass, respectively, 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  ; 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, . A useful resource for calibration
for calibration is the geometric section of SN 1987A which is given as a sketch in Fig. 5 of . 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.
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 , and the X-map in the 0.4 - 5.0 keV band of Fig. 1 in . The following observed radii can be extracted: in the polar direction and 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, ,
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 .
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.
where is the theoretical radius of the SNR, 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.
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 , , , and give .
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 , and .
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 , , , and give .
Figure 13. 3D surface of SN 1006 with parameters as in Figure 12, inverse square profile. The three Euler angles are , and .
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 , , , , and give .
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 , , , , and give .
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 , , , and give .
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 , , , and give .
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 .
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 , , , , , , , and give .
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 , , , , , , , and give .
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, ,
where A is the considered area, V the velocity and the density, see Formula (A28) in . In our case , 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
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 . 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 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 which contains 4003 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 from the first surface. As an example, we fixed , 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 with a variable integer number according to Formula (43).
· Each point of has spatial coordinates which can be repre-sented by the following matrix, A,
The orientation of the object is characterized by the Euler angles and therefore by a total rotation matrix, E, see . The matrix point is represented by the following 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, , given by the observational techniques, is now analysed. In both models the threshold intensity can be parametrized by , the maximum value of intensity characterizing the map.
More details on the theory of the image can be found in .
4.2. Results for the First Model
The radiation of the SN is supposed to originate from between a lower boundary of radius , given as an example by 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
where n is the considered index. In order to build the first model we assume that z varies between and 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 . 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 .
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 , and . In this map .
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 , and . In this map .
Table 1. Synoptic percentage reliability for the best model of SNRs with different profiles for the density.
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.
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.
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).
At the time of writing, an animated version of Figure 11 is visible at http://personalpages.to.infn.it/~zaninett/image/sn1987a_animation.gif.
 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.
 Stafford, J. and Lopez, L.A. (2017) Measuring the Symmetry of Supernova Remnants in the Radio. American Astronomical Society Meeting Abstracts #229, 148-155.
 Collier, A., Bachrach, H., Fryer, C. and Ellinger, C. (2017) Asymmetry in Supernovae. American Astronomical Society Meeting Abstracts #229, 434-438.
 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.
 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.
 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.
 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.
 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.
 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.
 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.