{'fraction_non_alphanumeric': 0.038836560506035414, 'fraction_numerical': 0.02787061845703395, 'mean_word_length': 4.133134836239898, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The physics of asymmetric supernovae and supernovae remnants

We model the circumstellar medium with four density profiles: hyperbolic type, power law type, exponential type and Gaussian type. We solve analytically or numerically the four first-order differential equations which arise in the framework of the classical thin layer approximation. The non-cubic dependence of the swept mass with the advancing radius is also considered. We derive the equation of motion for the thin layer approximation in special relativity in two cases. The initial conditions are chosen in order to model the temporal evolution of SN 1987A over 23 years and of SN 1006 over 1000 years. We review the building blocks of the symmetrical and asymmetrical formations of the image.


Introduction
The term nova derives from the Latin novum, which means new. As an example, a 'Stella Nova', now SN 1604 , was observed both by Kepler in Italy, see [1], and Galileo, see [2,3] as well as by Korean astronomers, see [4,5].
More recent is the use of the term supernova (SN), see [6]. The expanding shell of matter from an SN, consisting of the supernova ejecta and swept-up matter, has been called an SN remnant (SNR), see [7]. Astronomers initially thought that the equation for an SN was symmetric, and the following two laws of motion, among others, were used. In the Sedov solution, the radius, r, scales as r ∝ t 0.4 , where t is the time, see [8,9]. In the momentum conservation for a thin layer approximation the radius scales as r ∝ t 0.25 , see [10,11]. Later the following two terms appeared: non-spherical SNR, see [12], and symmetry in SN explosions, see [13]. The observed asymmetries can be classified when the external surface is rotated in front of the observer. Often an equatorial or reflectional symmetry is visible in the astronomical images. This is the case with the bipolar planetary nebula He 2-104, see [14], the protoplanetary nebula M2-9, see [15], the young planetary nebula MyCn 18, see [16] and the nebula around Eta Carinae, see [17]. As an example of a rotation which brings the asymmetrical object face on, see Figure 1 in [18]. The magnetic fields around planetary nebulae have been observed by [19,20] and this observational fact has triggered simulations of the asymmetrical shape of planetary nebulae based on variations of the magnetic field, see [19,21,20]. In the following, the terms 'asymmetrical', 'asymmetric', 'non-symmetric', and 'aspherical', are supposed to have an identical meaning, i.e. they are synonyms. The SNRs can therefore be classified in light of the observed symmetry. A first example is SN 1993J , which presented a circular symmetry for 4000 days, see [22]. An example of a weak departure from circular symmetry is SN 1006 , in which a ratio of 1.2 between the maximum and minimum radius has been measured, see [23]. An example of axial symmetry is SN 1987A , in which three rings are symmetric with respect to a line which connects the centres, see [24]. We now briefly review some results on the expansion velocity of SNs. The spectropolarimetry (CA II IR triplet) of SN 2001el gives a maximum velocity of ≈ 26000 kms −1 , see [25]. The same triplet when searched in seven SNs of type Ia gives 10400 kms −1 ≤ v ≤ 17700 kms −1 , see Table I in [26]. A time series of eight spectra in SN 2009ig allows asserting that the velocity at the CA II line, for example, decreases in 12 days from 32000 kms −1 to 21500 kms −1 , see Fig. 9 in [27].
A recent analysis of 58 type Ia SNs in Si II gives 9660 kms −1 ≤ v ≤ 14820 kms −1 , see Table II in [28]. Other examples for the maximum velocity of expansion are ≈ 3700 kms −1 , see Fig. 20.21 in [29] and Fig. 6 in [30]. The previous analysis allows saying that the maximum velocity so far observed for an SN is v c ≈ 0.123, where c is the speed of light; this observational fact points to a relativistic equation of motion.
The above arguments leave some questions unanswered or only partially answered: • Is it possible to deduce an equation of motion in the framework of the thin layer approximation adopting different density profiles?
• Is it possible to deduce an equation of motion for an expanding shell assuming that only a fraction of the mass enclosed in the advancing sphere is absorbed in the thin layer?
• Is it possible to model the complex three-dimensional (3D) behaviour of the velocity field of an expanding nebula?
• Is it possible to make an evaluation of the reliability of the numerical results for the radius, compared to the observed values?
• Can we reproduce complicated features, such as the equatorial ring + two outer rings in SN 1987A , which are classified as 'Mysterious Rings', see [31]?
• Can we reproduce the appearance of the rim of SN 1006 , sometimes classified as 'A contrail from an alien spaceship' or 'A jet from a black hole', see [32].
• Is it possible to build cuts of the model intensity which can be compared with existing observations?

The density profiles
This section outlines the meaning of section for a 3D surface, the adopted symmetries, and then postulates the existence of four density profiles for the circumstellar medium (CSM) around the SNs: an hyperbolic profile, a power law profile, an exponential profile, and a Gaussian profile. The above profiles produce a reflectional symmetry about the Cartesian plane z = 0. We briefly recall that the density of the galactic H I scales with a Gaussian dependence as where z is the vertical distance in pc from the galactic plane, n H (0) = 1.11 particles cm −3 , b = 75.5 pc, and z < 1000 pc see [33]. In the following we will explore gradients in density with smaller values in respect to the scale h of the galactic H I and we will see which gradient produces the best results.

Geometrical section of a 3D surface
As an example, we take an ellipsoid in Cartesian coordinates which has the equation z 2 a 2 + where the lengths a, b, and c are the semi-axes. The intersection with the plane y = 0 is displayed in Figure 1 and Figure 2 displays the resulting section, which is an ellipse.   In the following, the meaning of geometrical section will be associated with the intersection between a complex 3D surface and the plane x = 0 or y = 0.

Spherical coordinates
A point in Cartesian coordinates is characterized by x, y, and z and the position of the origin is the center of the SN explosion. The same point in spherical coordinates is characterized by the radial distance r ∈ [0, ∞], the polar angle θ ∈ [0, π], and the azimuthal angle ϕ ∈ [0, 2π]. Figure 3 presents a geometrical section of an asymmetric SN in which is clearly seen the polar angle θ and the three observable radii R up , R down , and R eq . We now outline two methods to measure the asymmetries. The first method evaluates the ratio Rup Req which models the difference between the observed radius in the polar direction in respect to the equatorial direction. We can therefore speak of weak observed asymmetry when 1 < Rup Req < 2 and great observed asymmetry when Rup Req > 2. As an example, the efficiency along the equatorial direction is where R eq,theo is the theoretical radius given by the model along the equatorial direction.
The second method is connected with the availability of observed and theoretical geometrical sections of SN/SNR. An observational percentage reliability, obs , is introduced over the whole range of the polar angle θ, where r num is the theoretical radius given by a model, r obs is the observed radius, and the index j varies from 1 to the number of available observations. The above statistical method allows fixing the parameters of the theory in a scientific way rather than adopting an 'ad hoc' hypothesis.

A hyperbolic profile
The density is assumed to have the following dependence on z in Cartesian coordinates, where the parameter z 0 fixes the scale and ρ 0 is the density at z = z 0 . In spherical coordinates the dependence on the polar angle is Given a solid angle ∆Ω the mass M 0 swept in the interval [0, r 0 ] is The total mass swept, M (r; r 0 , z 0 , α, θ, ρ 0 ), in the interval [0, r] is M (r; r 0 , z 0 , α, θ, ρ 0 ) = 1 3 ρ 0 r 0 3 + 1 2 The density ρ 0 can be obtained by introducing the number density, n 0 , expressed in particles cm −3 , the mass of hydrogen, m H , and a multiplicative factor f , which is chosen to be 1.4, see [9], The astrophysical version of the total swept mass, expressed in solar mass units, M , is therefore M (r pc ; z 0,pc , n 0 , θ) ≈ 0.0172 n 0 z 0,pc r pc where z 0,pc , r 0,pc and r 0,pc are z 0 , r 0 and r expressed in pc.

A power law profile
The density is assumed to have the following dependence on z in Cartesian coordinates: where z 0 fixes the scale. In spherical coordinates, the dependence on the polar angle is The mass M 0 swept in the interval [0, r 0 ] in a given solid angle is The total mass swept, M (r; r 0 , α, θ, ρ 0 ), in the interval [0, r] is M (r; r 0 , α, z 0 , θ, ρ 0 ) = The astrophysical swept mass is M (r pc ; z 0,pc , α, n 0 , θ) ≈ 0.03444 n 0 z 0 ,pc where z 0,pc and r pc are z 0 and r expressed in pc.

An exponential profile
The density is assumed to have the following exponential dependence on z in Cartesian coordinates: where b represents the scale. In spherical coordinates, the density is The total mass swept, M (r; r 0 , b, θ, ρ 0 ), in the interval [0, r] is The astrophysical version expressed in solar masses is where r 0,pc , r pc and b pc are r 0 , r and b expressed in pc.

A Gaussian profile
The density is assumed to have the following Gaussian dependence on z in Cartesian coordinates: where b represents the standard deviation. In spherical coordinates, the density is The total mass swept, M (r; r 0 , b, θ, ρ 0 ), in the interval [0, r] is where erf(x) is the error function, defined by The previous formula expressed in solar masses is M (r pc ; r 0,pc , b pc , θ, n 0 ) = 1 (cos(θ)) 3 − −0.011481 r 0 ,pc 3 (cos (θ)) 3 +0.03444 e where r 0,pc , r pc and b pc are r 0 , r and b expressed in pc.

The classical thin layer approximation
This section reviews the standard equation of motion in the case of the thin layer approximation in the presence of a CSM and derives the equation of motion for the density profiles here analysed. One case of non-cubic dependence for the swept mass as a function of the radius is also considered.

Classical momentum conservation
The conservation of the classical momentum in spherical coordinates along the solid angle ∆Ω in the framework of the thin layer approximation states that where M 0 (r 0 ) and M (r) are the swept masses at r 0 and r, and v 0 and v are the velocities of the thin layer at r 0 and r. This conservation law can be expressed as a differential equation of the first order by inserting v = dr dt : A variant of the above equation introduces the non-cubic dependence (NCD), p, which 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 along the solid angle ∆Ω is 1 3 ρ 0 r 3 0 where ρ 0 is the central density of the ambient medium. The swept mass included in the thin layer which characterizes the expansion is Let us call M (r) the mass swept between 0 and r along the solid angle ∆Ω. The conservation of momentum in NCD gives The above two differential equations are independent of the azimuthal angle ϕ. The 3D surface which represents the advancing shock of the SN is generated by rotating the curve in the x − z plane defined by the analytical or numerical solution r(t) about the z-axis and this is the first symmetry. A second symmetry around the z = 0 plane allows building the two lobes of the advancing surface. The orientation of the 3D surface is characterized by the Euler angles (Θ, Φ, Ψ) and therefore by a total 3 × 3 rotation matrix, E, see [34].
The adopted astrophysical units are pc for length and yr for time; the initial velocity v 0 is expressed in pc yr −1 . The astronomical velocities are evaluated in km s −1 and therefore v 0 = 1.02 × 10 −6 v 1 where v 1 is the initial velocity expressed in km s −1 .

Motion with an hyperbolic profile
In the case of an hyperbolic density profile for the CSM, as given by equation (5), the differential equation which models momentum conservation is where the initial conditions are r = r 0 and v = v 0 when t = t 0 . The variables can be separated and the radius as a function of the time is The velocity as a function of the radius r is

Motion with a power law profile
In the case of a power-law density profile for the CSM as given by equation (11), the differential equation which models the momentum conservation is We now evaluate the integral which is The solution of the differential equation (34) can be found by solving numerically the following nonlinear equation

Motion with an exponential profile
In the case of an exponential density profile for the CSM as given by equation (16), the differential equation which models momentum conservation is There is no analytical solution. We present the following series solution of order 4 around t 0 A second approximate solution can be found by extracting the velocity from and Given a function f (r), the Padé approximant, after [35], is where the notation is the same as in [36]. The coefficients a i and b i are found through Wynn's cross rule, see [37,38], and our choice is o = 2 and q = 1. The choice of o and q is a compromise between precision, needing high values for o and q, and simplicity of the expressions to manage, needing low values for o and q. The inverse of the velocity expressed by the the Padè approximant is where The above result allows deducing a solution r 2,1 expressed through the Padè approximant where and

Motion with a Gaussian profile
In the case of a Gaussian density profile for the CSM as given by equation (20), the differential equation which models momentum conservation is There is no analytical solution.

Motion with an exponential profile, NCD
In the case of an exponential density profile for the CSM as given by equation (16), the differential equation which models momentum conservation in the NCD case, see equation (28), is with The above differential equation does not have an analytical solution and should be solved numerically.

The relativistic thin layer approximation
This section reviews the conservation of momentum in special relativity and subsequently introduces two laws of motion.

Relativistic momentum conservation
The conservation of relativistic momentum in spherical coordinates along the solid angle ∆Ω in the framework of the thin layer approximation gives Here, M 0 (r 0 ) and M (r) are the swept masses at r 0 and r, β = v c , β 0 = v 0 c , and v 0 and v are the velocities of the thin layer at r 0 and r, see [39,40] for more details. We have chosen as units pc for distances and yr for time, and therefore the speed of light is c = 0.306 pc yr −1 .

Relativistic motion with a hyperbolic profile
In the case of a hyperbolic density profile for the CSM as given by equation (5), the differential equation which models the relativistic momentum conservation is where the initial conditions are r = r 0 and v = v 0 when t = t 0 . The velocity expressed in terms of β can be derived from the above equation: where There is no analytical solution of (56), but we present the following series solution of order three around t 0 : More details can be found in [41].

Relativistic motion with an exponential profile
In the case of an exponential density profile for the CSM as given by equation (16), the differential equation which models the relativistic momentum conservation is where There is no analytical solution of (61), so we present the following series solution of order three around t 0 :

The simulated objects
This section reviews the astronomical observations of SN 1987A and SN 1006 .

Great asymmetry, SN 1987A
The complex structure of SN 1987A , see an Hubble Space Telescope (ST) image in Figure 4, can be classified as a torus only, a torus plus two lobes, and a torus plus 4 lobes, see [42,43]. The region connected with the radius of the advancing torus is here identified with our equatorial region, in spherical coordinates, θ = π 2 . The radius of the torus only as a function of time can be found in Table 2 of [42] or Fig.  3 of [44], see Figure 5 for a comparison of the two different techniques. The radius of the torus only as given by the method of counting pixels ( [44]) shows a more regular behaviour and we have calibrated our codes in the equatorial region in such a way that at time 23 years the radius is 0. 39 2 pc. Another useful resource for calibration is the geometrical section of SN 1987A reported as a sketch in Fig. 5 of [45]. This geometrical section was digitized and rotated in the x − z plane by −40 • , see Figure 6.
An hydrodynamical simulation of the torus only for SN 1987A is reported in Figure 4 of [46].

Weak asymmetry, SN 1006
This SN started to be visible in 1006 AD and currently has a radius of 12.2 pc, see [47,48]. The X-shape can be visualized in Figure 7 and 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 [23], and the X-map in the 0.45.0 keV band of Fig. 1 in [47]. 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 geometrical section of the above X-map was digitized and rotated in the x − z plane by −45 • , see Figure 9.
Many filament radial intensity profiles at 27 keV for SN 1006 can be found in Fig. 8 of [49]. According to [47], the velocity of the ejecta should be ≈ 3100 km s .

Astrophysical Results
This section reports the theoretical results compared with the astronomical observations for SN 1987A and SN 1006 .

The case SN 1987A
In the case of an hyperbolic density profile, see equation (6), we have an analytical solution for the motion, see equation (30). Figure 10      rotation around the z-axis of the previous geometrical section allows building a 3D surface, see Figure 11. In the case of a power law density profile, see equation (11), we evaluate a numerical solution for the motion, see the nonlinear equation (38). Figure 12 reports the numerical solution for the power law case as a cut of SN 1987A in the x − z plane.
In the case of an exponential density profile, see equation (16), the differential equation which regulates the motion is equation (39). Figure 13 compares the numerical solution, the approximate series solution, and the Padé approximant solution. The above figure clearly shows the limited range of validity of the power series solution. The good agreement between the Padé approximant solution and numerical solution, in Figure 13 the two solutions can not be distinguished, has a percentage error where r(t) is the numerical solution and r(t) 2,1 is the Padé approximant solution. Figure 14 shows the percentage error as a function of the polar angle θ. Figure 15 shows a cut of SN 1987A in the x − z plane evaluated with the numerical solution.
In the case of a Gaussian density profile, see equation (20), the differential equation (51). regulates the motion. Figure 16 shows a cut of SN 1987A in the x−z plane evaluated with a numerical solution for the Gaussian profile. In this        (1). In the case of an exponential density profile, see equation (16), the differential equation (61) regulates the relativistic motion. Figure 17 shows the numerical solution for SN 1987A in the x − z plane for the relativistic-exponential case.

The case SN 1006
In the case of an hyperbolic density profile, Figure 18 displays a cut of SN 1006 and Figure 19 a 3D surface.
In presence of an exponential density profile for the CSM as given by equation (16), in the NCD case we solve numerically the differential equation (54) with parameters as in Figure (1).

Symmetrical synchrotron radiation
This section introduces the transfer equation, the acceleration and losses for electrons, the mathematical diffusion in 3D, a simple model for the image formation, and a diffusive model.

The transfer equation
The transfer equation in the presence of emission only, see, for example, [50] or [51], is where I ν is the specific intensity, s is the line of sight, j ν is the emission coefficient, k ν is a mass absorption coefficient, ρ is the density of mass at position s, and the index ν denotes the frequency of emission. The intensity  Table 1, exponential profile and NCD case.  Table 1 and efficiency along the equatorial direction eq = 98.25%.
of radiation, i.e. the photon flux, is here identified with the counts at a given energy. The solution to Equation (64) is where τ ν is the optical depth at frequency ν: The framework of synchrotron emission, as described in sec. 4 of [52], is often used in order to explain an SNR, see for example [53,54,55,56,57]. The volume emissivity (power per unit frequency interval per unit volume per unit solid angle) of the ultrarelativistic electrons, according to [58], is where P (ν) is the total power radiated per unit frequency interval by one electron and N (E)dE is the number of electrons per unit volume, per unit solid angle along the line of sight, which are moving in the direction of the observer and whose energies lie in the range E to E + dE. In the case of a power law spectrum, where K is a constant. The value of the constant K can be found by assuming that the probability density function (PDF, in the following) for relativistic energy is of Pareto type as defined in [59]: with c > 0. In our case, c = γ − 1 and a = E min , where E min is the minimum energy. In the previous formula, we can extract where N 0 is the total number of relativistic electrons per unit volume, here assumed to be approximately equal to the matter number density. The previous formula can also be expressed as where m H is the mass of hydrogen. The emissivity of the ultrarelativistic electrons from a homogeneous and isotropic distribution of electrons whose N (E) is given by Equation (68) is, according to [58], ≈ 0.933 × 10 −23 α L (γ)KH (γ+1)/2 ⊥ 6.26 × 10 18 ν where ν is the frequency and α L (γ) is a slowly varying function of γ which is of the order of unity and is given by for γ ≥ 1 2 . We now continue to analyse a first case of an optically thin layer in which τ ν is very small (or k ν is very small) and where the density ρ is replaced by the concentration C(s) of relativistic electrons: where K e is a constant function of the energy power law index, magnetic field, and frequency of e.m. emission. The intensity is now The increase in brightness is proportional to the concentration integrated along the line of sight. A second case analyses a quadratic relationship between the emission coefficient and the concentration where K 2 is a constant function. This is true for • Free-free radiation from a thermal plasma, see formula (1.219) in [58] .
The intensity in the "thermal case" is

Acceleration and losses for electrons
An electron which loses energy due to synchrotron radiation has a lifetime of where E is the energy in ergs, H the magnetic field in gauss, and P r is the total radiated power, see Eq. 1.157 in [58]. The energy is connected to the critical frequency, see Eq. 1.154 in [58], by ν c = 6.266 × 10 18 HE 2 Hz .
The lifetime for synchrotron losses is A first astrophysical result is the evaluation of the distance, L syn , that a ultrarelativistic electron which radiates at a given frequency ν when is moving with the advancing shell at a given average velocity, v SN , where β SN = v SN c . As a practical example L syn = 3.84 10 −3 pc when β SN = 1/100, ν = 1.0 10 9 Hz and H = 1 Gauss.
Considering a given radius R SN of the advancing SN, we should speak of the 'in situ' acceleration of the electrons when In the framework of the 'in situ' acceleration of the electrons and following [60,61], the gain in energy in a continuous form for a particle which spirals around a line of force is proportional to its energy, E, where τ II is the typical time scale, where u is the velocity of the accelerating cloud belonging the advancing shell of the SN/SNR, c is the speed of light, and L II is the mean free path between clouds, see Eq. 4.439 in [58]. The mean free path between the accelerating clouds in the Fermi II mechanism can be found from the following inequality in time: which corresponds to the following inequality for the mean free path between scatterers The mean free path length for an SN/SNR which emits synchrotron emission in the radio region at 1.0 10 9 Hz, gives L < 5.4618 10 −4 u 3000 where u 3000 is the velocity of the accelerating cloud expressed in units of 3000 km/s and H 1 is the magnetic field expressed in units of 1 gauss. When this inequality is fulfilled, the direct conversion of the rate of kinetic energy into radiation can be adopted. Recall that the Fermi II mechanism produces an inverse power law spectrum in the energy of the type N (E) ∝ E −γ or an inverse power law in the observed frequencies N (ν) ∝ ν β with β = − γ−1 2 [58,62]. The strong shock accelerating mechanism referred to as Fermi I was introduced by [63] and [64]. The energy gain relative to a particle that is crossing the shock is where τ I is the typical time scale, where u is the velocity of the shock of the SN/SNR, c is the speed of light, and L I is the mean free path between scatterers. This process produces an energy spectrum of electrons of the type see [65]. The two mechanisms, Fermi I and Fermi II, produce the same results when More details can be found in [66].

3D diffusion from a spherical source
Once the number density, C, and the diffusion coefficient, D, are introduced, Fick's first equation changes its expression on the basis of the adopted environment, see for example equation (2.5) in [67]. In three dimensions, it is where t is the time and ∇ 2 is the Laplacian differential operator. In the presence of the steady state condition, The number density rises from 0 at r=a to a maximum value C m at r=b and then falls again to 0 at r=c, see Figure 22, which is adapted from Fig.  3.1 of [67].
The solution to equation (93) is where A and B are determined by the boundary conditions, and These solutions can be found in [67] or in [68].

Spherical Image
A first thermal model for the image is characterized by a constant temperature in the internal region of the advancing sphere. We recall that in the thermal model the temperature and the density are constant in the advancing shell, see equation (10.29) in [69]. We therefore assume that the number density C is constant in a sphere of radius a and then falls to 0. The line of sight, when the observer is situated at the infinity of the x-axis, is the locus parallel to the x-axis which crosses the position y in a Cartesian x − y plane and terminates at the external circle of radius a, see [70]. The length of this locus is The number density C m is constant in the sphere of radius a and therefore the intensity of the radiation is A comparison of the observed data of SN 1006 with the theoretical thermal intensity is displayed in Figure 23.
A second non-thermal model for the image is characterized by synchrotron emission by ultrarelativistic electrons in a thin layer around the advancing sphere. We therefore assume that the number density C is constant in the spherical thin layer, and in particular rises from 0 at r = a to a maximum value C m , remains constant up to r = b and then falls again to 0. The line of sight, when the observer is situated at the infinity of the x-axis, is the locus parallel to the x-axis which crosses the position y in a Cartesian x − y plane and terminates at the external circle of radius b, see [70]. The length of this locus is   The number density C m is constant between two spheres of radii a and b and therefore the intensity of radiation is The ratio between the theoretical intensity at the maximum (y = a) and at the minimum (y = 0) is given by A comparison of the observed data of SN 1993J and the theoretical nonthermal intensity is displayed in Figure 24.
The main result of this paragraph is that the intensity of the thermal model which has the maximum of the intensity at the centre of an SNR does not match the observed profiles. The observed profiles in the radio intensity have the maximum value at the rim, as predicted by the non-thermal model.

Diffusive model
The concentration rises from 0 at r=a to a maximum value C m at r=b and then falls again to 0 at r=c. The concentrations to be used are formulas (95) and (96) once r = x 2 + y 2 is imposed; these two concentrations are inserted in formula (65), which represents the transfer equation. The geometry of the phenomenon fixes three different zones (0 − a, a − b, b − c) for the variable y, see [71,70]; the first segment, I I (y), is The second segment, I II (y), is The third segment, I III (y), is The profile of I which is made by three segments (102), (103) and (105), can be calibrated against the real data of SN 1006 as in Fig. 8 (subfigure 14) of [49] and is reported in Figure 25.

Asymmetrical synchrotron radiation
This section reviews a numerical algorithm which allows building a theoretical map for the intensity of radiation, and then applies the algorithm to SN 1987A and SN 1006 .

The adopted algorithm
We assume two models for the luminosity: a first non thermal model and a second thermal model. The source of synchrotron luminosity in the first model is assumed here to be the flux of kinetic energy, L m , where A is the considered area, V the velocity and ρ the density, see formula (A28) in [72]. 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 is a constant of conversion from the mechanical luminosity to the observed luminosity in synchrotron emission. The above formula is derived in the framework of the 'in situ' acceleration of electrons, see formula (82). The numerical algorithm which allows us to build a complex image will now be outlined.
• 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 (106) in the non-thermal case or to formula (76) in the thermal case.
• 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 [34]. 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. The threshold intensity can be parametrized by I max , the maximum value of intensity characterizing the map.

The case of SN 1987A
The image of SN 1987A having polar axis along the z-direction is shown in Figure 26. The image for a realistically rotated SN 1987A is shown in Figure  27 and a comparison should be made with the observed triple ring system as reported in Fig. 1 of [24] where the outer lobes are explained by light-echo effects. More details can be found in [73].  Table 1 and non-thermal case. The image of SN 1006 is visible in different bands, such as radio, see [74,75], optical, see [76], and X-ray, see [77,78]. The 2D map in intensity of SN 1006 is visible in Figure 28. The intensity along the equatorial and polar directions of our image is reported in Figure 29; a comparison should be made with Fig.  4 in [77]. The projected flux as a function of the position angle is another interesting quantity to plot, see Figure 30, and a comparison should be made with Fig. 5 top right in [79].
After the previous graphs, it is more simple to present a characteristic feature such as the "jet appearance" visible in some maps, see our Figure  31. A comparison of the above figure should be made with the ASCA X-map in Figure 8, see nomenclature 'Rim', the X-map for the 6.33-6.53 keV band visible in Fig. 3b of [80], the X-map for the 2-7 keV band in Fig. 7 in [49], and the X-map for 2-4.5 keV in Fig. 2 in [81].  Table 1. The three Euler angles characterizing the orientation of the observer are Φ=90 • , Θ=-55 • and Ψ=-180 • . This combination of Euler angles corresponds to the observed image for the non-thermal case. In this map I tr = I max /1.1.

Conclusions
Type of medium: We have selected four density profiles, which decrease with the distance (zaxis) from the equatorial plane. The integral which evaluates the swept mass increases in complexity according to the following sequence of density profiles: hyperbolic, power law, exponential, and Gaussian.
Classical thin layer The application of the thin layer approximation with different profiles produces differential equations of the first order. The solution of the first-order differential equation can be analytical in the classical case, characterized by a hyperbolic density profile, see (30), and numerical in all other cases. We also evaluated the approximation of the solution as a power law series, see (40), or using the Padé approximant, see (48): the differences between the two approximations are outlined in Figure 13. In the case of an exponential density profile for the CSM as given by equation (16), the NCD approximation was implemented. The presence of a gradient for the density transforms a spherical symmetry into an axial symmetry and allows the appearance of the so called 'bipolar motion'.

Relativistic thin layer
The application of the thin layer approximation to the relativistic case produces first-order differential equations which can be solved only numerically or as a power series, see (59). and (62).

Symmetric image
The intensity of the image of a symmetrical SN or SNR in the case of an optically thin medium can be computed through an analytical evaluation of the lines of sight when the number density is constant between two spheres, see formula (100). In the case of a symmetrical diffusive process which is built in presence of three spheres, the intensity of emission is assumed to be proportional to the number density, see formulas ( 102), ( 103) and ( 105).
Asymmetric image A first model for the emissivity in the advancing layer assumes a proportionality to the flux of kinetic energy, see equation (105) where the density is assumed to be proportional to the swept material. A second model assumes an emissivity proportional to the square of the number density, see formula (76). This second allows simulating particular effects, such as the triple ring system of SN 1987A , see Figure 27. Another curious effect is the "jet appearance" visible in the weakly symmetric SN 1006 , see Fig. 31; this jet's effect is simulated in the framework of the first non thermal model. The jet/counter jet effect plays a significant role in current research, see the discussion in Section 5.2 of [82], where the jet appearance is tentatively explained by neutrino heating, see [83], or by the MHD jet, see [84]. Here in contrast we explain the appearance of the jet by the addition of three effects: • An asymmetric law of expansion due to a density gradient with respect to the equatorial plane. This produces an asymmetry in velocity.
• The direct conversion of the flux of kinetic energy into radiation.
• The image of the SNR as the composition of integrals along the line of sight.
A careful calibration of the various parameters involved can be done when cuts in intensity are available.