A classical and a relativistic law of motion for SN1987A

In this paper we derive some first order differential equations which model the classical and the relativistic thin layer approximations in the presence of a circumstellar medium with a density which is decreasing in the distance $z$ from the equatorial plane. The circumstellar medium is assumed to follow a density profile with $z$ of hyperbolic type, power law type, exponential type or Gaussian type. The first order differential equations are solved analytically, or numerically, or by a series expansion, or by Pad\'e approximants. The initial conditions are chosen in order to model the temporal evolution of SN 1987A over 23 years. The free parameters of the theory are found by maximizing the observational reliability which is based on an observed section of SN 1987A.


Introduction
The theories of the expansion of supernovae (SN) in the circumstellar medium (CSM) are usually built in a spherical framework. Unfortunately, only a few SNs present a spherical expansion, such as SN 1993J , see Marcaide et al. (2009);Martí-Vidal et al. (2011). The more common observed morphologies are barrel or hourglass shapes, see Lopez (2014) for a classification. A possible classification for the asymmetries firstly identifies the center of the explosion and then defines the radius in the equatorial plane, R eq , then the radius in the downward direction, R down , and then the radius in the upward direction, R up , see Zaninetti (2000). The above classification allows introducing a symmetry: R down =R up means that the expansion from the equatorial plane along the two opposite polar directions is the same. A second symmetry can be introduced in the framework of spherical coordinates assuming independence from the azimuthal angle. The theories for asymmetric SNs or late supernova remnants (SNRs) have therefore been set up, we select some of them. Possible reasons for the distortion of SNRs have been extensively studied analytically and numerically by Chevalier and Gardner (1974); Bisnovatyi-Kogan et al. (1989); Igumenshchev et al. (1992); Arthur and Falle (1993); Maciejewski and Cox (1999). Two SNRs presenting a barrel morphology were observed and explained in Gaensler (1998). Numerical calculations of the interaction of an SN with an axisymmetric structure with a high density in the equatorial plane were carried out by Blondin et al. (1996).
New laws of motion, assuming that only a fraction of the mass which resides in the surrounding medium is accumulated in the advancing thin layer, were developed both in a classical framework in the presence of an exponential profile, see Zaninetti (2012) or an isothermal self-gravitating disk, see Zaninetti (2013), and both in a relativistic framework in the presence of an isothermal selfgravitating disk, see Zaninetti (2014). We now present some maximum observed velocities in SNs: the maximum velocity for Si II λ6355 vary in [15000,25000] km s −1 according to Figure 13 in Silverman et al. (2015) or in [13000,24000] km s −1 according to Figure 4 in Zhao et al. (2015). These high observed velocities demand a relativistic treatment of the theory. In this paper we introduce, in Section 2, four asymmetric density profiles, in Section 3 we derive the differential equations which model the thin layer approximation for an SN in the presence of four asymmetric types of medium and a relativistic treatment is carried out in Section 4 for two asymmetric types of medium.

Preliminaries
This section introduces the spherical coordinates and four density profiles with axial symmetry for the CSM: a hyperbolic profile, a power law profile, an exponential profile and a Gaussian profile.

Spherical coordinates
A point in Cartesian coordinates is characterized by x, y, and z. 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 1 presents a section of an asymmetric SN where there can be clearly seen the polar angle θ and the three observable radii R up , R down , and R eq .
The density ρ 0 can be obtained by introducing the number density expressed in particles cm −3 , n 0 , the mass of hydrogen, m H , and a multiplicative factor f , which is chosen to be 1.4, see McCray (1987), The astrophysical version of the total swept mass, expressed in solar mass units, M ⊙ , is therefore where z 0,pc , r 0,pc and r 0,pc are z 0 , r 0 and r expressed in pc units.

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 The astrophysical swept mass is

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; The astrophysical version expressed in solar masses is where b pc is the scale 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; where erf(x) is the error function, defined by The previous formula expressed is solar masses is M (r; r 0,pc , b pc , θ, n 0 ) = 1 (cos(θ)) 3 −1.024 10 −7 n 0 −1.4 10 6 r 0 ,pc 3 (cos(θ)) 3

The classical thin layer approximation
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 : The above differential equation is independent of the azimuthal angle ϕ. The 3D surface which represents the advancing shock of the SN consists of the rotation about the z-axis of the curve in the x − z plane defined by the analytical or numerical solution r(t); 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 Goldstein et al. (2002). 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 .

The case of SN 1987A
The complex structure of SN 1987A can be classified as a torus only, a torus plus two lobes, and a torus plus 4 lobes, see Racusin et al. (2009). 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 Racusin et al. (2009) or Figure 3 of Chiad et al. (2012), see Figure 2 for a comparison of the two different techniques. The radius of the torus only as given by the counting pixels method (Chiad et al. (2012)) shows a more regular behavior and we calibrate 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 a section of SN 1987A reported as a sketch in Figure 5 of France et al. (2015). This section was digitized and rotated in the x − z plane by −40 • , see Figure 3. The above approximate section allows introducing an observational percentage reliability, ǫ obs , over the whole range of the polar angle θ, Figure 3: Section of SN 1987A in the x − z plane adapted by the author from Figure 5 in France et al. (2015).
where r num is the theoretical radius, r obs is the observed radius, and the index j varies from 1 to the number of available observations, in our case 81. The above statistical method allows fixing the parameters of the theory in a scientific way rather than adopting an "ad hoc" hypothesis.

Motion with an hyperbolic profile
In the case of a hyperbolic density profile for the CSM as given by Eq.
(1), 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 and The velocity as a function of the radius r is v(t) = 2 r 0 3 v 0 cos (θ) 2 r 0 3 cos (θ) − 3 r 0 2 z 0 + 3 r 2 z 0 .

Motion with a power law profile
In the case of a power-law density profile for the CSM as given by Eq. (7), the differential equation which models the momentum conservation is We now evaluate the following integral which is The solution of the differential equation (30) can be found solving numerically the following nonlinear equation More precisely we used the FORTRAN SUBROUTINE ZBRENT from Press et al. (1992) and Figure 6 reports the numerical solution as a cut of SN 1987A in the x − z plane.

Motion with an exponential profile
In the case of an exponential density profile for the CSM as given by Eq. (12), the differential equation which models momentum conservation is An analytical solution does not exist and we present the following series solution of order 4 around t 0 A second approximate solution can be found by deriving the velocity from (35): and Given a function f (r), the Padé approximant, after Padé (1892), is f (r) = a 0 + a 1 r + · · · + a p r p b 0 + b 1 r + · · · + b q r q , where the notation is the same of Olver et al. (2010). The coefficients a i and b i are found through Wynn's cross rule, see Baker (1975); Baker and Graves-Morris (1996) and our choice is p = 2 and q = 1. The choice of p and q is a compromise between precision, high values for p and q, and simplicity of the expressions to manage, low values for p and q. The inverse of the velocity expressed by the the Padè approximant is where N 21 = (r − r 0 )(9 e − r 0 cos(θ) b br − 9 e − r 0 cos(θ) b br 0 + 2 cos(θ)rr 0 − 2 cos(θ)r 0 2 − 4 br + 10 r 0 b) (42) and D21 = 2 v 0 cos (θ) rr 0 − cos (θ) r 0 2 − 2 br + 5 r 0 b The above result allows deducing a solution r 2,1 expressed through the Padè approximant where where r(t) is the numerical solution and r(t) 2,1 is the Padé approximant solution. Figure 8 shows the percentage error as a function of the polar angle θ. Figure  9 shows a cut of SN 1987A in the x − z plane evaluated with the numerical solution.

Relativistic motion with an exponential profile
In the case of an exponential density profile for the CSM as given by Eq. (12), the differential equation which models the relativistic momentum conservation is where An analytical solution of (56) does not exist, so we present the following series solution of order three around t 0 : (57) Figure 12 shows the numerical solution for SN 1987A in the x − z plane.

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 (25), and numerical in all other cases. We also evaluated the approximation of the solution as a power law series, see (36), or using the Pade approximant, see (44): the differences between the two approximations are outlined in Figure 7.
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 (57) and (54).
The astrophysical case The application of the theory here developed is connected with a clear definition of the advancing SN's surface in 3D. We have concentrated the analysis on SN 1987A with the cuts of the advancing surface in the x − z plane when the z axis is in front of the observer. Another choice of the point of view of Figure 12: Section of SN 1987A in the x − z plane with an exponential profile: relativistic case (green points) and observed profile (red stars). The parameters r 0 =0.1 pc, b = 0.02 pc, t = 21.86 yr, t 0 = 0.1 yr, v 0 = 70000 km s −1 , and β 0 = 0.233 give ǫ obs = 90.23%. the observer would complicate the situation, and a comparison between theory and observation then requires the introduction of the three Euler angles which characterize the observer, see the rotated advancing surface of SN 1987A shown in Figure 5. Is interesting note that the rotation of the observed image with the polar axis aligned with the z-direction has been done for the Homunculus nebula, see Figure 4 in Smith (2006), but only an approximate section of the H-α imaging connected with SN 1987A has been already reported, see Figure  5 in France et al. (2015). A more precise definition of the section of SN 1987A will help the theoretical determination of the parameters maximizing the observational reliability, ǫ obs , see 23. Is important note that the observational reliability gives already an acceptable result and ǫ obs lies within the interval [88%-92%] for the six models, four classical and two relativistic, here analyzed.