Free Access
Issue
A&A
Volume 636, April 2020
Article Number A22
Number of page(s) 19
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201936718
Published online 08 April 2020

© ESO 2020

1. Introduction

The structure and morphology of supernova remnants (SNRs), which are the outcome of supernova (SN) explosions, reflect the properties of the parent SNe and the characteristics of their stellar progenitor systems. Thus, establishing firm connections between SNRs and their parent SNe and progenitor stars is essential for obtaining key information on many aspects of the explosion processes associated with SNe and of the final phases of stellar evolution. A unique opportunity to study these connections is offered by the remnant of SN 1987A in the Large Magellanic Cloud thanks to its youth (≈30 years) and proximity (≈50 kpc).

Observations of this remnant in different wavelength bands and at different epochs have revealed Doppler shifts in emission lines of heavy elements (e.g., [Fe II] and [Ni II]) up to velocities > 3000 km s−1 (e.g., Haas et al. 1990; Colgan et al. 1994; Utrobin et al. 1995; Larsson et al. 2016), suggesting an aspherical explosion and the production of high-velocity clumps of metal-rich ejecta in the SN. This scenario is also supported by the detection of emission lines from decay of 44Ti around day ≈10 000 (Boggs et al. 2015) which are redshifted with Doppler velocities of ≈700 ± 400 km s−1. All these lines of evidence point to large-scale asymmetry in the SN explosion, possibly a signature of the engine powering the burst.

Identifying the physical and geometrical properties of initial SN asymmetries is crucial to constraining explosion processes associated with SNe. Seminal studies (e.g., Patnaude et al. 2015; Wongwathanarat et al. 2015, 2017; Orlando et al. 2015, 2016; Ferrand et al. 2019) have shown that hydrodynamic/magnetohydrodynamic (HD/MHD) models can be very effective in studying the stellar progenitor-SN-SNR connection. Current models, however, treat the problem in a piecemeal way mainly due to the very different time and space scales involved, and the intrinsic three-dimensional (3D) nature of the phenomenon. Additional difficulty stems from the need to disentangle, in the remnant morphology, the effects of the SN explosion and of the structure of the stellar progenitor, from those of the interaction of the blast with the inhomogeneous ambient environment. As a result, many aspects of the stellar progenitor-SN-SNR connection remain uncertain.

Here we link the dynamical and radiative properties of the remnant of SN 1987A to the asymmetric SN explosion and to the structure of the progenitor star, through a comprehensive 3D hydrodynamic model which describes the evolution from the onset of the SN to the development of its remnant, accounting for the pre-SN structure of the progenitor star. To capture the large range in space scales (spanning 13 orders of magnitude) during the evolution, we coupled a 3D model of a core-collapse SN (Ono et al. 2020, in the following Paper I) with a 3D model of a SNR (Orlando et al. 2019a).

The paper is organized as follows. In Sect. 2, we describe the model, the numerical setup, and the synthesis of X-ray emission. In Sect. 3, we discuss the results and in Sect. 4, we draw our conclusions.

2. Problem description and numerical setup

We described the evolution of SN 1987A from the onset of the SN to the full-fledged remnant following a strategy analogous to that outlined in previous works (Orlando et al. 2015, 2016). First we performed 3D core-collapse SN simulations (see Sect. 2.2) which describe the onset of the SN, the propagation of the shock wave through the stellar interior, and the breakout of the shock wave at the stellar surface. Thus, the model accounts for the mixing and clumping of the stellar envelope and mantle, that is the ejecta, at the time of shock breakout. The SN simulations were initialized from pre-SN stellar models available in the literature (see Sect. 2.1) and they ended soon after the shock breakout, covering ≈20 h of evolution. Then the output of these simulations was used as an initial condition for the structure and chemical composition of ejecta in 3D MHD simulations which describe the subsequent evolution: the transition from the phase of SN to that of SNR and the interaction of the remnant with the inhomogeneous nebula around SN 1987A (see Sect. 2.3). The SNR simulations cover 50 years, thus describing the past evolution of SN 1987A until its current age and predicting its future evolution (until 2037).

2.1. Pre-supernova models

The progenitor star of SN 1987A, Sanduleak (Sk) −69° 202, was a blue supergiant (BSG) with an initial mass of ≈20 M (Hillebrandt et al. 1987). In the present study, we considered two progenitor star models of BSG (N16.3 and B18.3), which have been proposed in the literature to describe the average radial structure of Sk −69°  202 at collapse, and a model of red supergiant (RSG; model S19.8) selected to explore the evolution and structure of the remnant in the case of a progenitor star which differs significantly from the structure of a BSG. Figure 1 shows the radial profiles of ρr3 (where ρ is the mass density, and r is the radial distance from the center of the star) and the mass fractions of selected elements versus the interior mass, Mr, of the pre-SN models considered here. The radial profile of ρr3 provides useful information on the SN shock dynamics because its gradient determines where and when the SN blast wave decelerates (positive gradient) or accelerates (negative gradient).

thumbnail Fig. 1.

Radial profiles of ρr3 of the pre-SN models considered in this paper (upper left panel) and mass fractions of selected elements (see legend in each panel) versus the interior mass Mr in the pre-SN models N16.3 (upper right; Nomoto & Hashimoto 1988; Shigeyama & Nomoto 1990), B18.3 (lower left; Urushibata et al. 2018), and S19.8 (lower right; Woosley et al. 2002; Sukhbold et al. 2016).

Open with DEXTER

The first model, N16.3, denotes a 16.3 M BSG that was proposed early to describe Sk −69°  202 (Nomoto & Hashimoto 1988; Shigeyama & Nomoto 1990). The model assumes a main-sequence mass of 20 M which reduces to 16.3 M at the onset of the SN due to mass loss. For this reason, the model is also known as N20 in the literature (Wongwathanarat et al. 2015). At the time of the core-collapse, the model is characterized by a 6 M helium core and a 10.3 M compact hydrogen envelope with a radius of 3.4 × 1012 cm (≈49 R). In such a way, the BSG model satisfies the helium core mass, hydrogen envelope mass, luminosity and radius of Sk −69°  202 at collapse (Shigeyama & Nomoto 1990).

The second model, B18.3, was proposed recently (Urushibata et al. 2018) in the framework of the slow-merger scenario (Ivanova et al. 2002) and it denotes a 18.3 M BSG resulting from the merging of two massive stars with 14 M and 9 M, respectively. This model is analogous to the merger model previously proposed by Menon & Heger (2017). These models are particularly interesting because they reproduce most of the observational constraints of Sk−69°  202, such as the red-to-blue evolution, its lifetime (≈20 000 years), the total mass (≈18.3 M) and position in the Hertzsprung-Russell diagram (log Teff = 4.2, log L/L = 4.9) at collapse. Model B18.3 also predicts a surface chemical composition before envelope deflation of 0.139, 3.47 and 1.21 for He/H, N/C and N/O, respectively, in agreement with observations. Model B18.3 describes an external envelope which is quite similar to that of model N16.3 (see Fig. 1). Otherwise, the internal structure of B18.3 differs significantly from that of N16.3. In particular, model B18.3 is characterized by a C+O core with a smaller size and with a profile of ρr3 flatter than that of N16.3 (see Fig. 1). We expect, therefore, that the SNRs from these two models differ mainly for the structure and distribution of innermost ejecta rather than for outermost ejecta.

In addition to the above models, we examined the model S19.8 (Woosley et al. 2002; Sukhbold et al. 2016) which denotes a non-rotating, solar metallicity progenitor RSG with a main-sequence mass 19.8 M which reduces to 15.9 M at collapse. This model has a structure that is completely different from that of models N16.3 and B18.3 (see Fig. 1). However, it predicts helium core mass, hydrogen envelope mass, and luminosity consistent with those of Sk−69°  202 (although the stellar radius is much larger than that of Sk−69°  202). The C+O layer and the helium core mass are also similar to those of N16.3, although the ρr3 gradient of the helium layer is very different in the two models (Sukhbold et al. 2016) (see Fig. 1). Model S19.8 was considered in order to investigate how sensitive the remnant evolution is to the structure of the progenitor star.

2.2. Modeling the evolution of the supernova

We followed the evolution of the core-collapse SN by adopting the hydrodynamic model presented in Ono et al. (2013) but extended it to three dimensions. In Paper I, we describe in detail the 3D model and the corresponding numerical setup. In the following, we summarize the main features of the model.

The SN evolution was modeled by numerically solving the time-dependent hydrodynamic equations of mass, momentum, and energy conservation in a 3D Cartesian coordinate system (x, y, z). The model accounts for the effects of gravity (both self-gravity and gravitational effects of the central proto-neutron star), the fallback of material on the proto-neutron star, the explosive nucleosynthesis through a nuclear reaction network, the feedback of nuclear energy generation, and the energy deposition due to radioactive decays of isotopes synthesized in the SN explosion. The adopted equation of state depends on the different physical regimes simulated (see Paper I): at early phases, for 10−10 g cm−3 <  ρ <  1011 g cm−3 and 104 K< T <  1011 K, we adopted the Helmholtz equation of state (Timmes & Swesty 2000), which includes contributions from radiation, fully ionized nuclei, and degenerate or relativistic electrons and positrons; at later phases, we considered an equation of state that includes contributions from the radiation and ideal gas of elements (Ono et al. 2013) for ρ <  10−8 g cm−3, and a smooth transition between the latter equation of state and that of Helmholtz for 10−8 g cm−3 <  ρ <  10−7 g cm−3. The pressure from radiation was neglected in the optically thin regime by adopting an approximate approach described by Joggerst et al. (2010).

The SN explosion is initiated by injecting kinetic and thermal energies artificially around the composition interface of the iron core and the silicon-rich layer. We modeled aspherical explosions with clumpy structures to mimic the global anisotropies developing in magnetic jet-driven SNe (Khokhlov et al. 1999; Maeda & Nomoto 2003; Couch et al. 2009; Bear & Soker 2018) or, alternatively, in the context of neutrino heating mechanism, by convection in the neutrino heating layers and the standing accretion shock instability (SASI; e.g., Kifonidis et al. 2006; Marek & Janka 2009; Suwa et al. 2010; Nordhaus et al. 2010; Takiwaki et al. 2012; Wongwathanarat et al. 2015; Janka et al. 2016; Janka 2017; Walk et al. 2018). First we set the initial radial velocity as

(1)

where f(θ) is a function of θ in spherical coordinates selected in order to produce a bipolar-like explosion in the light of the observed elliptical shape of the inner ejecta distribution in SN 1987A (Larsson et al. 2013, 2016). More specifically, we explored four cases (see Paper I for details): a function with cosine, an exponential form, a power-law form, and an elliptical form. The degree of asymmetry of the explosion is regulated by the parameter β which is the ratio of the radial velocity on the polar axis, vpol, to that on the equatorial plane, veq, at a given radius (β = vpol/veq). Then we imposed asymmetry across the equatorial plane by changing the normalization of vr across this plane; this is regulated by a parameter α = vup/vdw, where vup and vdw are the initial radial velocities at a radius inside the shock at θ = 0° and θ = 180°, respectively. This kind of asymmetry may be introduced by the combined action of neutrino heating and SASI. Finally, we introduced large amplitude perturbations (30%) to the initial radial velocities by superposing functions of (θ, ϕ) in spherical coordinates, based on real spherical harmonics (see Paper I for details). We explored the space of parameters of the explosion (energy and parameters regulating the initial asymmetry); Table 1 reports the range of parameters explored and the parameters of the best-fit models N16.3, B18.3, and S19.8 (corresponding to models n16.3-high, b18.3-high, and s19.8-fid respectively, in Table 2 of Paper I).

Table 1.

Explored range of parameters of the SN models, and parameters of the best-fit models considered in this paper.

The simulations were performed using the adaptive mesh hydrodynamic code FLASH (Fryxell et al. 2000). The code consists of inter-operable modules that can be combined to generate different astrophysical applications and it was designed to make efficient use of massive parallel computers using the message passing interface (MPI) library. For the 3D SN simulations, we used the directionally split Eulerian version of the piecewise parabolic method (PPM; Colella & Woodward 1984), which provides a formal second-order accuracy in both space and time. Also a hybrid Riemann solver which switches to an HLLE solver inside shocks was used to avoid an odd-even instability (decoupling) that can develop from shocks aligned with the grid (Quirk et al. 1997).

The gravitational effects are included by adopting a spherically symmetric approximation1: first we calculated angle-averaged radial profiles of density and, then, we estimated local gravitational potentials from enclosed masses at each radius. A point source gravity from the mass of the proto-neutron star was also included; the total mass that accretes onto the proto-neutron star during the evolution was added to the point mass (see Paper I for details).

The nuclear reaction network includes 19 species: neutron (n), proton (p), 1H, 3He, 4He, 12C, 14N, 16O, 20Ne, 24Mg, 28Si, 32S, 36Ar, 40Ca, 44Ti, 48Cr, 52Fe, 54Fe, and 56Ni (see Weaver et al. 1978 for the network chain). The calculations were performed by using the MA28 sparse matrix package (Duff et al. 1986) and a time integration scheme based on the Bader-Deuflhard method (Bader & Deuflhard 1983), both available in FLASH. The distribution and evolution of elements were followed by adopting a multiple fluids approach. Each fluid is associated to one of the heavy elements of the nuclear reaction network; the fluid evolution is calculated by an advection equation which is solved in addition to the hydrodynamic equations. The energy depositions due to radioactive decays of 56Ni to 56Co to 56Fe are calculated by adopting the method described by Joggerst et al. (2009) (see also Ono et al. 2013 and Paper I for more details).

In order to follow the large physical scales spanned during the remnant expansion, we followed a mesh strategy similar to that proposed by Ono et al. (2013). The initial computational domain is a Cartesian box extending between ≈ − ​5000 km and ≈5000 km in the x, y, and z directions, with the SN explosion assumed to sit at the origin of the coordinate system (x0, y0, z0) = (0, 0, 0). In such a way, the domain includes the inner regions of the oxygen-rich layer of the progenitor star at collapse. In this mesh, we exploited the adaptive mesh refinement (AMR) capabilities of FLASH (through the PARAMESH library; MacNeice et al. 2000) to reach a maximum spatial resolution of 4.9 km in low-resolution simulations and 2.4 km in high-resolution simulations. In Paper I, we describe in detail the AMR strategy adopted for our SN simulations. Then the computational domain was gradually extended as the forward shock propagates: when the forward shock was close to one of the boundaries of the Cartesian box, the physical size of the computational domain was extended by a factor of 1.2 in all directions. The physical quantities of the old domain were re-mapped in the new domain; the physical quantities in the extended region were set to the values of the pre-SN star and of the wind (after the shock breakout). This approach is possible because the propagation of the forward shock is supersonic and does not introduce errors larger than 0.1% after 40 re-mappings (Ono et al. 2013). We found that about 75 re-mappings were necessary to follow the propagation of the SN blast wave through the star until the breakout of the shock wave at the stellar surface, about 20 h after the core-collapse. The final domain extends between ≈ − ​2 × 1014 cm and ≈2 × 1014 cm in the x, y, and z directions, leading to a maximum spatial resolution of ≈2 × 1011 cm.

2.3. Modeling the evolution of the supernova remnant

The output of the SN simulations about 20 h after the core-collapse was used as initial condition for 3D MHD simulations describing the interaction of the blast wave with the ambient medium. The numerical setup for the SNR model has been described in detail in Orlando et al. (2019a). Briefly, the evolution and transition from the SN phase to the SNR phase were modeled by numerically solving the time-dependent MHD equations of mass, momentum, energy, and magnetic flux conservation in a 3D Cartesian coordinate system (x, y, z).

The calculations were performed using PLUTO (Mignone et al. 2007, 2012), an advanced multi-dimensional Godunov-type code for astrophysical plasmas and high Mach number flows. The code was extended by additional computational modules to calculate the deviations from temperature-equilibration between electrons and ions, and the deviations from equilibrium of ionization of the most abundant ions (see Orlando et al. 2019a for the details of the implementation). The former were implemented by including the almost instantaneous heating of electrons at shock fronts up to kT ≈ 0.3 keV by lower hybrid waves (Ghavamian et al. 2007) and the effects of Coulomb collisions for the calculation of ion and electron temperatures in the post-shock plasma (Orlando et al. 2015). The deviations from equilibrium of ionization were calculated through the computation of the maximum ionization age in each cell of the spatial domain (Orlando et al. 2015).

We followed the chemical evolution of the ejecta by adopting a multiple fluids approach (Orlando et al. 2016). Each fluid is associated to one of the species modeled in the SN simulations and initialized with the corresponding abundance in the output of the SN simulations. As for the composition of the circumstellar medium (CSM) around SN 1987A, we adopted the abundances inferred from the analysis of X-ray observations (Zhekov et al. 2009). This allowed us to follow the chemical evolution of the ejecta and to map the spatial distribution of heavy elements at different epochs during the evolution of SN 1987A. The different fluids mix together during the remnant evolution: the density of a specific element in a fluid cell at time t is calculated as ρel = ρCel, where Cel is the mass fraction of each element and the index “el” refers to a different element.

We modeled the CSM as in previous works (Orlando et al. 2015, 2019a). The immediate surrounding of the SN event was modeled as a spherically symmetric wind characterized by a gas density proportional to r−2 (where r is the radial distance from the center of explosion), a mass-loss rate of w = 10−7 M yr−1, and a wind velocity uw = 500 km s−1 (Morris & Podsiadlowski 2007). The termination shock of the wind is located at rw = 0.05 pc. The circumstellar nebula consists of an extended ionized H II region, a dense inhomogeneous equatorial ring, and two less dense rings lying in planes almost parallel to the equatorial one, but displaced by about 0.4 pc above and below the central ring. The H II region is uniform with density nHII; its inner edge in the equatorial plane is at distance rHII from the center of explosion. The central ring is composed of a uniform smooth component and Ncl high-density spherical clumps mostly located in its inner portion (Chevalier & Dwarkadas 1995; Sugerman et al. 2005). The uniform component has density nrg, radius rrg, and an elliptical cross section with the major axis wrg lying on the equatorial plane and height hrg. The high-density clumps have a diameter wcl and their plasma density and radial distance from SN 1987A are randomly distributed around the values ⟨ncl⟩ and ⟨rcl⟩ respectively. We explored the space of parameters of the CSM around the fiducial values derived in a previous work (Orlando et al. 2015). Table 2 reports the range of parameters explored and the best-fit parameters found for the three progenitor stars considered (N16.3, B18.3, and S19.8).

Table 2.

Explored and best-fit parameters of the CSM for the models of SN 1987A.

In order to trace the evolution of the different plasma components (ejecta, H II region, and ring material), we introduced passive tracers, each associated with a different component. Each tracer is initialized to one in cells belonging to the plasma component and to zero elsewhere. The continuity equations of the tracers are solved in addition to our set of MHD equations. Some tracers are also used to store information on the shocked plasma (time, shock velocity, and shock position, namely Lagrangian coordinates, when a cell of the mesh is shocked either by the forward or by the reverse shock) which are required to synthesize the thermal X-ray emission (see Sect. 2.4). These last tracers are initialized to zero.

The ambient magnetic field (most likely that originating from the progenitor star) is modeled as in Orlando et al. (2019a) as a “Parker spiral”, in analogy with the spiral-shaped magnetic field on the interplanetary medium of the solar system (Parker 1958). This is the simplest field configuration for a rotating magnetized progenitor star, resulting from the rotation of the star and from the expanding stellar wind. For our purpose, we considered a pre-SN magnetic field characterized by an average strength at the stellar surface B0 ≈ 3 kG (well within the range inferred for magnetic massive stars, e.g., Donati & Landstreet 2009) and a strength at the inner edge of the nebula (at a distance from the center of explosion r ≈ 0.08 pc) B1 ≈ 1 μG. This pre-SN magnetic field configuration leads to field strength within the nebula in agreement with observations (Zanardo et al. 2018). As found in Orlando et al. (2019a), the field plays a role in preserving inhomogeneities of the CSM (as the rings and the clumps) from complete fragmentation by limiting the growth of hydrodynamic instabilities that would develop at their boundaries (Orlando et al. 2008), whereas it does not influence the overall expansion and evolution of the blast wave.

We adopted a mesh strategy similar to that used to model the evolution of the SN explosion to follow the large physical scales spanned during the remnant expansion (see Sect. 2.2 and Orlando et al. 2019a). The initial computational domain is a Cartesian box extending between −3 × 1014 cm and 3 × 1014 cm in the x, y, and z directions, thus including the final domain of the SN simulations (see Sect. 2.4). For the evolution of the SNR, the box was covered by a non-uniform grid with the smallest mesh size (and the highest spatial resolution) in the remnant interior. The grid is made of 12803 zones and consists of the following: (i) a uniform grid patch with 5123 points, extending between −7.5 × 1013 cm and 7.5 × 1013 cm in the x, y, and z directions (including metal-rich ejecta), and leading to a maximum resolution of ≈2.9 × 1011 cm, and (ii) a non-uniform grid patch in the rest of the domain with resolution ranging between ≈2.9 × 1011 cm close to the high-resolution patch and ≈5.8 × 1011 cm close to the border of the domain. This non-uniform mesh allows us to describe in more detail the distribution of metal-rich ejecta in the remnant interior. Then, as for SN simulations, we gradually extended the computational domain by a factor 1.2 as the forward shock propagates outward, and remapped the physical values to the new domains. All the physical quantities in the extended region are set to the values of the pre-SN CSM. We found that about 49 re-mappings were necessary to follow the interaction of the blast wave with the CSM during 50 years of evolution. The final domain extends between −1.3 × 1018 cm and 1.3 × 1018 cm in the x, y, and z directions, leading to a maximum spatial resolution of ≈1.3 × 1015 cm in the remnant interior and of ≈2.6 × 1015 cm in the outer layers of the remnant. All physical quantities were set to the values of the pre-SN CSM at all boundaries.

2.4. Synthesis of X-ray emission

From the model results, we synthesized the thermal X-ray emission originating from the impact of the blast wave with the nebula, by following the approach outlined in previous studies (Orlando et al. 2015; Miceli et al. 2019). The system was rotated about the three axes to fit the orientation of the ring in the plane of the sky inferred from the analysis of optical data (Sugerman et al. 2005): ix = 41°, iy = −8°, and iz = −9°. Then, for each jth cell of the spatial domain, we derived: the maximum ionization age τj = nejΔtj (where nej is the particle number density in the j-th domain cell and Δtj is the time since the plasma in the cell was shocked), the emission measure em (where is the hydrogen number density in the cell and Vj is the cell volume), and the electron temperature Tej (derived from the ion temperature, plasma density, and Δtj by assuming Coulomb collisions and starting from an electron temperature at the shock front kT = 0.3 keV, which is assumed to be the same at any time as a result of instantaneous heating by lower hybrid waves; Ghavamian et al. 2007; Orlando et al. 2015). The X-ray emission was synthesized in the [0.1, 10] keV band from the values of τj, emj, and Tej, assuming that the system is at a distance D = 51.4 kpc (Panagia et al. 1999), and by using the non-equilibrium of ionization (NEI) emission model VPSHOCK available in the XSPEC package along with the NEI atomic data from ATOMDB (Smith et al. 2001). For the CSM, we adopted the metal abundances derived from the analysis of deep Chandra observations of SN 1987A (Zhekov et al. 2009); for the ejecta, we used the metal abundances calculated during the ejecta evolution after the core-collapse. Finally, the X-ray spectrum from each cell was filtered through the photoelectric absorption by the interstellar medium, assuming a column density NH = 2.35 × 1021 cm−2 (Park et al. 2006), and folded through the instrumental response of either XMM-Newton-EPIC or Chandra-ACIS. X-ray lightcurves were derived by integrating the spectra from the cells in spectral bands of interest and in the whole spatial domain; X-ray images were derived by integrating the spectra in bands of interest and along the line-of-sight (LoS).

3. Results

3.1. Early distribution and mixing of ejecta in the remnant

We ran several 3D high-resolution simulations of the SN, searching for the parameters (explosion energy, Eexp, and parameters for the initial asymmetry, α and β; see Sect. 2.2) best reproducing the observed shifts and broadening of [Fe II] lines (Haas et al. 1990; Colgan et al. 1994; Utrobin et al. 1995). Table 1 reports the range of parameters explored and the best-fit parameters for each pre-SN models adopted.

In Paper I, we describe in detail the evolution of the SN from the energy release soon after the core-collapse to the breakout of the shock wave at the stellar surface, covering about 20 h of evolution. In all the cases considered, the initial large scale asymmetry leads to an expanding bipolar structure. During the early phases, the bipolar structure is subject to development of Kelvin–Helmholtz (KH) and Rayleigh–Taylor (RT) instabilities. The growth of these instabilities depends on the degree of asymmetry of the explosion: the higher the asymmetry, the faster their growth. The morphology of the bipolar structure depends on the progenitor star models: it is the narrowest in model B18.3 due to the flatter ρr3 profile and smaller size of the C+O core (see Fig. 1), which allow for a rapid expansion of the shock wave. The initial bipolar structure is roughly kept up to the shock breakout even if its morphology can be significantly modified during the shock propagation due to structure of the stellar interior and the growth of hydrodynamic instabilities. In Paper I, we found that an elliptical form for the function f(θ) in Eq. (1) is required to reproduce the bipolar-like explosion of SN 1987A; Fig. 2 shows the corresponding angle dependence of the initial radial velocity. Thus, in the present paper, we focused the analysis only on this case.

thumbnail Fig. 2.

Angle dependence of the initial radial velocity in the [x, z] plane of the Cartesian coordinate, adopting an elliptical form for the function f(θ) in Eq. (1). The different curves are for different values of β = vpol/veq and α = vup/vdw = 1.5.

Open with DEXTER

After the shock breakout, the blast wave expands through the wind of the progenitor star before the impact onto the nebula. In this phase, the metal-rich ejecta expand almost homologously, carrying the fingerprints of the asymmetric explosion and of the structure of the progenitor star. The X-ray emission arising from the shocked wind and from the shocked ejecta is very faint due to the low values of wind density. In fact, no significant X-ray emission was detected during the first two years of evolution of SN 1987A. Indeed, in this phase, the expanding ejecta have been detected in the lines of [Fe II] 26 μm and 18 μm (e.g., Haas et al. 1990; Colgan et al. 1994; Utrobin et al. 1995). Thus, from the models, we synthesized these lines, taking into account the Doppler shift due to bulk motion of material along the LoS.

Following Haas et al. (1990), we considered several approximations and simplifications to synthesize the [Fe II] lines. In particular, we assumed that the temperature of the emitting iron is isothermal and that the level populations of iron depend only on the temperature (Haas et al. 1990). In fact, our models neglect the heating of ejecta due to the radioactive decay of 56Co and their cooling due to optical emission lines, and show small variations of the temperature of unshocked innermost ejecta at the epoch of observations. Nevertheless, even in the presence of these heating and cooling mechanisms, our assumption is reasonable because the heating and cooling rates scale similarly with the iron density (Colgan & Hollenbach 1988). We assumed also that the lines of [Fe II] 26 μm and 18 μm are produced by optically thin plasma. This assumption is supported by the analysis of Haas et al. (1990) which shows that most of the measured emission of [Fe II] lines comes from optically thin regions. Another assumption is that the emission of [Fe II] lines is proportional to the Fe density (Haas et al. 1990). In the light of the above assumptions, we derived the line profiles by simply calculating the distribution of iron mass vs. the LoS velocity and, then, convolved the profiles with a Gaussian of size 400 km s−1, to approximate the resolution of the observed spectra (Haas et al. 1990). Finally, we explored different explosion energies and asymmetries, along with different orientations of the asymmetry with respect to the LoS, to identify the physical and geometrical properties of the SN explosion which led to the shifts and broadening of [Fe II] lines in early observations (e.g., Haas et al. 1990; Colgan et al. 1994; Utrobin et al. 1995).

We found that a high energy of the explosion is required to reproduce the observations: at the shock breakout, our best-fit models have an explosion energy of ≈2 × 1051 erg (see Table 1). The explosion energy was constrained from the comparison between the maximum iron velocities reached in our models and those inferred from the analysis of [Fe II] lines during the first two years of evolution, and from the comparison between the X-ray lightcurves synthesized from the models and those observed during the first 30 years of evolution (see Sect. 3.2). Lower (higher) energies produce lower (higher) iron velocities and lower (higher) X-ray fluxes at odds with the observations. We note that the explosion energies found are at the upper end of the range of values, Eexp = 0.8 − 2.0 × 1051 erg, proposed in the literature for SN 1987A (see Table 1 in Handy et al. 2014). From bolometric lightcurve fitting what can be constrained is not the explosion energy itself but the ratio between the explosion energy, Eexp, and the mass of hydrogen envelope, Menv (e.g., Shigeyama & Nomoto 1990). Thus, the estimates of explosion energies from bolometric lightcurve fitting depend on the model adopted for the progenitor star. Studies that considered progenitor BSGs that evolved as single stars suggest explosion energies in the range 1.2 − 1.5 × 1051 erg (e.g., Utrobin et al. 2019). Studies that modeled explosions of BSGs from binary mergers (analogous to model B18.3) found slightly higher values, Eexp ≈ 1.7 × 1051 erg (Menon et al. 2019). The difference between merger models and single star models is basically due to the fact that the helium core mass in the former is rather smaller than that in the latter and that the envelope composition is significantly different in the two cases.

The two parameters α and β in Table 1 determine the level of asymmetry of the explosion and, therefore, the shifts and broadening of [Fe II] lines. The α parameter regulates the asymmetry across the equatorial plane so that its main effect is to shift the centroid of the lines. The β parameter is the ratio of the radial velocity on the polar axis to that on the equatorial plane. In this case, the main effect is to determine the broadening of Fe II lines. We started the exploration of the parameter space with low values of β and values of α in the range [1 − 2], and progressively increased β up to when the observed line broadening was reproduced by one of our models. We found that the observations require β = 16, indicating a high degree of asymmetry in the explosion with the radial velocity on the polar axis, vpol, 16 times larger than the radial velocity on the equatorial plane, veq, and with the radial velocity at θ = 0°, vup, 1.5 times larger than the radial velocity at θ = 180°, vdw. This kind of asymmetry may result from a strong sloshing of the stalled shock in polar directions, due to the SASI (e.g., Blondin et al. 2003), or, alternatively, in magnetic jet-driven SNe (e.g., Khokhlov et al. 1999), or in single-lobe SNe (e.g., Hungerford et al. 2005). As a consequence, most of the energy is released along the polar axis.

We note that, among the forms considered for the function f(θ) in Eq. (1), the elliptical is the only one reproducing the observed line profiles for the range of β explored (see Paper I for details). Values of β >  16 would produce lines broader than observed. We cannot exclude that other forms of f(θ) might reproduce the observed line profiles for β >  16. This, however, would indicate an even more extreme asymmetry in the explosion which seems to be unlikely. Finally, we note that a higher (lower) explosion energy would require a lower (higher) value of β. Explosion energies > 2 × 1051 erg would be out of the range of values inferred from the analysis of observed bolometric lightcurves and, as mentioned before, values of β >  16 are unlikely. Thus, the values selected (namely β = 16 and Eexp ≈ 2 × 1051 erg) for the favorite models are a compromise between very energetic and highly asymmetric explosions.

As for the orientation of the asymmetry with respect to the LoS, we found that the polar axis should be almost lying in the plane of the central ring of the nebula, roughly in the direction of the LoS but with an offset of 40°, and with the most energetic lobe propagating away from the observer in the south (see Fig. 3). The projection of this axis in the plane of the sky is offset by 15° from the northsouth axis.

thumbnail Fig. 3.

Schematic view of the orientation of the initial asymmetry (lightblue and red arrows) with respect to the dense inner ring (lightblue circle) and to the LoS (gray arrow along negative y-axis). The negative y-axis points toward the observer. The yardstick indicating the length scale is on the left.

Open with DEXTER

As a result of the asymmetry, after the shock breakout, most of the mass of 56Ni resides in two main clumps which propagate in opposite directions with an angle of 40° from the LoS: a massive one moves away from the observer in the south (red arrow in Fig. 3); the other (less massive) moves toward the observer in the north (lightblue arrow in the figure). These clumps may be those responsible for the so-called Bochum event observed at days 20–100, when the Hα profile exhibited two features nearly symmetric on both the blue and red wings of the Hα line (e.g., Hanuschik et al. 1988; Phillips & Heathcote 1989). The analysis of observations has suggested that these features are probably due to two 56Ni clumps ejected during the explosion (e.g., Utrobin et al. 1995; Nagataki et al. 1997; Nagataki 2000) and propagating in opposite directions with an angle ranging between 31° and 45° from the LoS (e.g., Utrobin et al. 1995; Wang et al. 2002). Our modeled 56Ni clumps fit well in this scenario.

The distribution of iron-rich ejecta at day 730 (the time of the first output of SNR simulations), after the decay of 56Ni in 56Fe, reflects the initial large-scale asymmetry (see Fig. 4; see also online Movie 1). The figure also compares the line profiles of [Fe II] 26 μm and 18 μm synthesized from the models with those observed (Haas et al. 1990). The only model able to reproduce the observed redshift and broadening of [Fe II] lines and, in particular, the wings of the profiles extending up to velocities ≈3000 km s−1, is B18.3, although it underestimates the mass of iron with low velocity (the dip at the center of the line) and does not reproduce the feature at v ≈ 4000 km s−1, due to a fast isolated clump of iron not modeled by our simulations. This suggests that the initial large-scale asymmetry was more complex than modeled here. The other models do not account for iron-rich ejecta with velocities > 2200 km s−1, for all the explosion energies and asymmetries explored (see Table 1).

thumbnail Fig. 4.

Isosurfaces of the iron distribution for models N16.3 (upper left panel), B18.3 (lower left), and S19.8 (lower right) at day 730. The semi-transparent isosurfaces correspond to a value of the iron density which is at 10% of the peak density for each model (maximum iron density is 4.9 × 10−16 g cm−3 for N16.3, 6.0 × 10−16 g cm−3 for B18.3, and 1.1 × 10−15 g cm−3 for S19.8). The colors give the velocity along the LoS in units of 1000 km s−1 on the isosurface; the color coding is defined at the bottom of each panel. The yardstick indicating the length scale is on the left. The upper right panel shows the line profiles of [Fe II] 26 μm and 18 μm of SN 1987A observed during the first 2 years of evolution (symbols; Haas et al. 1990) and synthesized from models N16.3, B18.3, and S19.8 at day 730 (solid lines). See online Movie 1 for an animation of these data; a navigable 3D graphic of the iron distribution for the three models is available at https://skfb.ly/6P7sH.

Open with DEXTER

The differences in the line profiles synthesized from the models can be further investigated by deriving the mass distributions of selected elements versus the radial velocity, vrad, at day 730. The left panels of Fig. 5 shows these distributions for models N16.3, B18.3, and S19.8. In the figure, ΔMi is the mass of the i-th element in the velocity range [v; v + Δv], and Mi is the total mass of the i-th element. The velocity is binned with Δv = 200 km s−1. At day 730, we assumed that all 56Ni has already decayed in 56Fe; thus, the figure reports only the distribution of the latter. The figure shows that the metal distribution can be strongly affected by the internal structure (density and temperature) of the progenitor star. In particular, the model adopting a progenitor RSG (S19.8) shows strong differences with respect to the other two cases which adopt a progenitor BSG (N16.3 and B18.3). This is a direct consequence of the distinct structure of the progenitor stars (see Fig. 1). In model S19.8, the various species (and in particular 1H and 4He) cannot reach velocities larger than ≈6000 km s−1 at variance with the other two models. This is mainly due to the structure of the extended hydrogen envelope of the RSG; in this case the blast wave is continuously decelerated during its propagation through the hydrogen envelope (see Paper I). As a consequence, the velocity of the forward shock after the breakout at the stellar surface is significantly lower in model S19.8 than in models N16.3 and B18.3, and ≈60% of the total explosion energy in S19.8 is kinetic (whilst, in the other two models, more than 90% of the explosion energy is kinetic).

thumbnail Fig. 5.

Mass distributions of 1H, 4He, 56Fe, 16O, 28Si, and 44Ti versus radial velocity vrad (left panels) and velocity component along the LoS, vLoS (right panels) at day 730 for models N16.3 (top), B18.3 (middle), and S19.8 (bottom).

Open with DEXTER

In the case of the two progenitor BSGs (N16.3 and B18.3), although their distributions are similar for all nuclear species, some evident differences appear at the high-velocity tail of the 56Fe, 16O, 28Si, and 44Ti distributions. In particular, in model B18.3, these species can reach higher velocities and they have more mass in the high-velocity tail of the distributions than in model N16.3. For instance, the fastest 56Fe (product of decay from 56Ni) mixes up to velocities of ≈5000 km s−1 in model B18.3, while it is about 10% slower in model N16.3. This is basically due to the fact that the velocity of the expanding nickel and, then, of its product of decay, iron, depends on the structure of the overlying C+O and He layers (see the gradients of ρr3 radial profiles in Fig. 1). High velocities can be obtained if the nickel is able to reach the hydrogen layer before the development of a strong reverse shock during the propagation of the blast wave through the He layer. In model B18.3, C+O and He layers are less extended and less dense than in model N16.3 (see Fig. 1). In this way, the nickel can penetrate efficiently through the helium and hydrogen layers. As a result, the metal-rich ejecta in model B18.3 can expand faster than in model N16.3, resulting in a larger matter mixing. Some differences between models N16.3 and B18.3 are evident also in the high velocity tails of 1H and 4He where, in N16.3, these species have higher velocities and more mass than in B18.3. One of the consequences is that the shock after the breakout at the stellar surface is slightly faster in model N16.3.

The lines of [Fe II] 26 μm and 18 μm (shown in Fig. 4) are expected to reflect the differences in the mass distributions versus the velocity component along the LoS, vLoS, through different Doppler shifts and broadenings. The right panels in Fig. 5 shows the mass distributions versus vLoS, assuming the orientation of the initial asymmetry required to fit the observations (see Fig. 3). The differences noted in the distributions versus vrad are also evident in the distributions versus vLoS. As a consequence, the line profiles reflect the distinct structure of the corresponding progenitor stars and, besides to constrain the explosion energy and asymmetry, enabled us to identify the progenitor star most appropriate for SN 1987A. The comparison of simulations with observations of SN 1987A (e.g., Haas et al. 1990) showed that the model best reproducing the observations is B18.3 (see Fig. 4).

3.2. Evolution of the X-ray emitting remnant

We evolved the models for 50 years, describing the interaction of the blast wave with the nebula. We ran several 3D high-resolution simulations of the SNR to explore the model parameter space, synthesizing (remote sensing) observables to compare the model results with available observations (see Sect. 2.4). In particular, we explored several density structures of the nebula within and around the range of parameters discussed in the literature (e.g., Sugerman et al. 2005), searching for a set of values that best reproduces the X-ray lightcurves and morphology of SN 1987A as observed in the last 30 years (e.g., Orlando et al. 2015). The range of parameters explored and the best-fit parameters found for the three models (N16.3, B18.3, and S19.8) are shown in Table 2.

We found that the remnant evolution is similar in models which consider a progenitor BSG (N16.3 and B18.3), whilst the evolution in the case of a progenitor RSG (model S19.8) largely differs from the others. The remnant hits the inhomogeneous nebula around year 3 in models N16.3 and B18.3 and around year 11 in model S19.8. In the latter case, the forward shock is slower than in the other models due to the early deceleration suffered by the shock during its propagation through the extended hydrogen envelope of the RSG (see Sect. 3.1). The impact of the blast wave onto the nebula produces a sudden increase of X-ray flux either due to significant X-ray emission arising from the shocked gas of the H II region (in models N16.3 and B18.3) or because of large X-ray emission produced by the shocked outermost ejecta (in model S19.8). Figure 6 shows the X-ray lightcurves synthesized from the three models. Online Movies 2, 3, and 4 show the evolution of the 3D volume rendering of particle density of the shocked plasma (left panel), and the corresponding synthetic X-ray emission maps in the [0.5, 2] keV and [3, 10] keV bands (upper right panels), and three-color composite representations (lower right panels), for models N16.3, B18.3, and S19.8, respectively; the synthetic X-ray emission maps are also shown in Fig. 7.

thumbnail Fig. 6.

X-ray lightcurves (solid lines) in the [0.5, 2] keV (upper panels) and [3, 10] keV (lower panels) bands synthesized from models N16.3 (left panels), B18.3 (center panels), and S19.8 (right panels) compared to the lightcurves of SN 1987A observed with Rosat (black triangles; Haberl et al. 2006), ASCA (green diamonds; Orlando et al. 2015), Chandra (magenta circles; Helder et al. 2013; Frank et al. 2016) and XMM-Newton (cyan circles; Haberl et al. 2006; Maggi et al. 2012). Dotted, dashed, and dot-dashed lines indicate the contribution to emission from the shocked ejecta, the shocked plasma from the H II region, and the shocked plasma from the ring, respectively.

Open with DEXTER

thumbnail Fig. 7.

Synthetic X-ray emission maps from models N16.3 (top), B18.3 (middle), and S19.8 (bottom) at the labeled times. Each panel consists of 4 quadrants: emission maps in the [0.5, 2] keV (upper left) and [3, 10] keV (upper right) bands, and corresponding three-color composite representations (lower left and right). Each image is normalized to its maximum for visibility; the X-ray maps were convolved with a Gaussian of size 0.15 arcsec to approximate the spatial resolution of Chandra; the three-color composite images were smoothed with a Gaussian of 0.025 arcsec; the colors in the composite show the emission from shocked ejecta (green), shocked plasma from the ring (red), and from the H II region (blue). See online Movies 2, 3, and 4 for animations of these data.

Open with DEXTER

During the remnant-nebula interaction, the evolution follows three phases as also found in previous studies (Orlando et al. 2015, 2019a): in the first phase, the blast wave travels through the H II region; in the second, it interacts with the dense equatorial ring; finally, in the third phase, it leaves the ring and starts to travel through a less dense environment. In models N16.3 and B18.3, the first phase occurs between years 3 and 15 when the X-ray emission is dominated by shocked material from the H II region and from the outermost ejecta (see also upper and middle panels in Fig. 7 and online Movies 2, and 3). In model S19.8, this phase lasts between year 11 (when the blast hits the H II region) and year 22 (when the blast reaches the ring) and, at variance with the other cases, the X-ray emission is always largely dominated by shocked ejecta (see right panels of Fig. 6, lower panels of Fig. 7, and online Movie 4).

The remnant enters the second phase when the blast wave hits the dense equatorial ring. This occurs around year 15 in models N16.3 and B18.3 and around year 22 in model S19.8. In this phase the forward shock travels through the ring, and the X-ray emission from the shocked material of the ring becomes the dominant component in models N16.3 and B18.3 (see Figs. 6 and 7, and online Movies 2, and 3). Conversely, the X-ray emission continues to be largely dominated by shocked ejecta in model S19.8 (see Figs. 6 and 7, and online Movie 4) and this makes this model totally different from the others. In this phase, the ambient magnetic field preserves the ring from erosion by hydrodynamic instabilities and the ring survives the passage of the blast (Orlando et al. 2019a). From year 26, the ring is fully shocked in models N16.3 and B18.3 and its X-ray emission gradually fades away. In model S19.8, the ring is fully shocked around year 42.

During the whole evolution, the reverse shock travels through the inner envelope of the SNR and the X-ray emission from shocked ejecta gradually increases as the shock encounters higher and higher densities of ejecta. In models N16.3 and B18.3, the third phase of evolution begins around year 34, when the blast wave, after leaving the ring, travels through a less dense environment and the X-ray emission becomes dominated by shocked ejecta (Figs. 6 and 7, and online Movies 2, and 3). In model S19.8, the blast leaves the ring around year 42, but the emission continues to be dominated by shocked ejecta (see online Movie 4).

By comparing the synthetic X-ray lightcurves with those observed, the models best representing the X-ray observations are those considering a progenitor BSG (N16.3 and B18.3). In the case of model S19.8 (with a progenitor RSG), we did not find a set of parameters of the CSM compatible with that inferred from observations (e.g., Sugerman et al. 2005) which is able to reproduce the observed X-ray lightcurves; this is mainly due to the fact that, in this model, the X-ray emission is largely dominated by shocked ejecta and this flux contribution is well above the observed values during the whole interaction of the remnant with the nebula. The synthetic lightcurves in this case are almost insensitive to the parameters of the CSM. Clearly the X-ray emitting remnant of SN 1987A keeps memory of the structure of the progenitor star. As for models N16.3 and B18.3, both are able to reproduce in detail the observed lightcurves with similar parameters of the CSM (compatible with those inferred from observations; e.g., Sugerman et al. 2005). In fact, the structure of the outermost ejecta (namely those interacting with the nebula in the period analyzed) is similar in these two models (see Fig. 1).

To further investigate which model between N16.3 and B18.3 better fits the X-ray observations, we compared the size of the expanding X-ray emitting torus visible in synthetic images (see Fig. 7) with the size observed in actual X-ray observations. To this end, we analyzed the synthetic X-ray images in the [0.5 − 10] keV band, adopting a methodology analogous to that used for actual X-ray observations (Frank et al. 2016), and assuming a distance D = 51.4 kpc from SN 1987A (Panagia et al. 1999), the same used for the synthesis of X-ray emission. Figure 8 shows the best-fit radius of the torus vs time inferred from the analysis of actual observations (symbols; Frank et al. 2016) and those derived from the synthetic images (solid lines). The observed torus radii at different epochs are nicely reproduced by model B18.3 which is able to describe the rapid expansion of the torus between years 10 and 15 (when the blast was propagating through the H II region) and the slower expansion of the torus between years 15 and 30 (when the blast was propagating through the denser material of the equatorial ring). Then the model predicts a faster expansion of the torus after year 32, namely after the blast has left the ring and started to travel through a less dense environment.

thumbnail Fig. 8.

Best-fit radius of the X-ray emitting torus vs time inferred from the analysis of actual X-ray observations (symbols; Frank et al. 2016) and that derived from the analysis of synthetic X-ray images for models B18.3 (red solid line), N16.3 (blue), and S19.8 (green).

Open with DEXTER

The other two models fail to reproduce the observed expansion of the X-ray emitting torus: model N16.3 roughly reproduces the expansion when the blast wave was traveling through the H II region, but it overestimates the torus radius after year 15, namely when the blast wave was traveling through the ring; model S19.8 does not reproduce the observed expansion of the torus and largely underestimates its radius during the whole evolution. In fact, model N16.3 requires a size of the equatorial ring slightly larger than in B18.3 to fit the X-ray lightcurves (see Table 2 and Fig. 6); this is necessary because, in N16.3, the forward shock after the breakout at the stellar surface is faster than in B18.3 (due to the structure of the progenitor star; see Sect. 3.1). As a consequence, the radius of the X-ray emitting torus in N16.3 is overestimated when the blast is traveling through the equatorial ring. As it concerns model S19.8, it largely underestimates the torus size during the whole evolution due to a significant lower velocity of the forward shock caused by the early deceleration suffered by the shock during its propagation through the extended hydrogen envelope of the progenitor RSG (see Sect. 3.1).

We note that a smaller ring size in model N16.3 (necessary to fit the observed radii of the torus) would have produced a steepening of the soft X-ray lightcurve well before that observed when the blast hits the ring. In other words, model N16.3 can reproduce either the X-ray lightcurves or the size of the X-ray emitting torus, whilst model B18.3 is able to naturally reproduce both at the same time.

3.3. Distribution of 44Ti in the evolved remnant

During the interaction of the blast wave with the inhomogeneous nebula, the metal-rich ejecta in the remnant interior expand and form a bipolar structure which lies almost in the plane of the ring and which reflects the initial large-scale asymmetry (see Fig. 3). The early 56Ni clumps and the distribution of 56Fe (Fig. 4), all lie along this structure. The projection of the bipolar structure in the plane of the sky is offset by 15° from the northsouth axis and coincides with an elongated region, evident in maps of [Si I]+[Fe II] lines of near-infrared and optical observations of SN 1987A (Larsson et al. 2016), that lies roughly in the direction of the northsouth axis but with an offset by 10° −20°. The ejecta probed by these emission lines are dominated by radioactive input from 44Ti decay.

Figure 9 and the online Movie 5 show the spatial distribution of 44Ti for the three models around day ≈10 000, when NuSTAR observations have resolved the emission lines from decay of 44Ti at 67.87 and 78.32 keV (Boggs et al. 2015). As found for 56Ni and 56Fe, again most of the mass of 44Ti resides in two clumps moving in opposite directions along the bipolar structure (Fig. 9, and online Movie 5): the most massive clump propagates away from the observer in the south at an angle of 40° from the LoS. Thus, the ejecta powered by radiative decay of 44Ti are expected to lie preferentially along the bipolar structure, consistently with the observed maps of [Si I]+[Fe II] lines (Larsson et al. 2016).

thumbnail Fig. 9.

Isosurfaces of the titanium distribution for models N16.3 (upper left), B18.3 (lower left), and S19.8 (lower right) at day 9875. The semi-transparent isosurface corresponds to a value of the titanium density which is at 10% of the peak (maximum titanium density is 8.6 × 10−22 g cm−3 for N16.3, 6.9 × 10−22 g cm−3 for B18.3, and 2.3 × 10−21 g cm−3 for S19.8). The colors give the velocity along the LoS in units of 1000 km s−1 on the isosurface; the color coding is defined at the bottom of each panel. The ring indicates the position of the dense inner ring. The yardstick on the left indicates the length scale. Upper right panel: best-fit Gaussian for the 44Ti lines observed between 2012 and 2014 (black line; Boggs et al. 2015) and lines profile synthesized from models N16.3 (blue line), B18.3 (red line), and S19.8 (green line) at day 9875. The shaded region represents the 90% confidence area for the position of the peak and for the width of the Gaussian. See online Movie 5 for an animation of these data; a navigable 3D graphic of the titanium distribution for the three models is available at https://skfb.ly/6P7sI.

Open with DEXTER

The evident differences in the distributions of the models originate again from the distinct structure of the corresponding progenitor stars. In model N16.3, most of the mass of 44Ti is concentrated in a single clump elongated along the bipolar structure. The portion of the clump expanding away from the observer includes ≈61% of the total mass of 44Ti (i.e., ≈5.4 × 10−4M). In the other two models (B18.3 and S19.8), most of the mass resides in two clumps moving in opposite directions along the bipolar structure. In both cases, the most massive clump propagates away from the observer in the south with ≈64% of the titanium mass (total mass ≈5.7 × 10−4M) in B18.3, and ≈60% of the titanium mass (total mass ≈8.9 × 10−4M) in S19.8. Considering that the half life of 44Ti is 60 years, the total 44Ti mass decreases by about 27% after ≈10 000 days, leading to ≈3.9 × 10−4M for N16.3, ≈4.2 × 10−4M for B18.3, and ≈6.5 × 10−4M for S19.8.

We note that all the models overestimate the total mass of titanium inferred from observations which, in general, fall in the range [0.5 − 2]×10−4M (Seitenzahl et al. 2014; Boggs et al. 2015). This is basically due to the limited number of nuclei included in the adopted nuclear reaction network. In fact, the mass of titanium can be overestimated by a factor ≈3 if compared with the mass calculated by a larger nuclear reaction network (464 nuclei; Mao et al. 2015, and Paper I). Indeed, taking into account this factor, the models predicting a mass of 44Ti consistent with observations are N16.3 and B18.3 (total mass ≈1.3 × 10−4M and ≈1.4 × 10−4M, respectively).

NuSTAR observations around day ≈10 000 have revealed 44Ti lines redshifted with Doppler velocities of ≈700 ± 400 km s−1 (black line in the upper right panel of Fig. 9; Boggs et al. 2015). From the models, we synthesized the profile of 44Ti lines assuming a spectral resolution of 900 eV in the range of wavelengths of interest (around [60 − 90] keV), and taking into account the Doppler shift and broadening due to motion of titanium along the LoS. We fit the synthetic profiles with a Gaussian function, following the analysis of NuSTAR observations (see upper right panel of Fig. 9). Table 3 reports the velocity of the peak of the Gaussian, vrs, corresponding to the line redshift, and the standard deviation, σ, corresponding to the Doppler broadening after accounting for the NuSTAR spectral resolution at these energies. In all the models, the synthetic 44Ti lines are redshifted. However, model N16.3 produces a line profile (blue line in the figure) with redshift much lower than observed (≈260 km s−1) and which is not compatible with the observations. The other two models produce lines with redshift velocities of ≈370 km s−1 in B18.3 (red line in the figure) and ≈840 km s−1 in S19.8 (green line) which are both consistent with the observations within the 90% confidence area (shaded gray region in the figure) for the position and width of the Gaussian best-fitting the observations. In the three models, the redshift originates from the initial asymmetry which assumes a parameter α = vup/vdw = 1.5 (see Table 1); this means that the lobe propagating away from the observer is more energetic than that propagating toward the observer.

Table 3.

Parameters of the Gaussian functions best fitting the observed and synthetic lines of 44Ti.

3.4. Molecular structures in the evolved remnant

The structure of innermost ejecta can be a direct probe of anisotropies developed in the immediate aftermath of the SN. Recently, this structure has been resolved by ALMA infrared observations of cold molecular gas (Abellán et al. 2017) that have shown that the distributions of carbon and silicon monoxide (CO and SiO) are characterized by a torus-like structure perpendicular to the equatorial ring.

Our models do not describe the formation of molecules. Nevertheless, we can evaluate the size scales and spatial distributions of CO and SiO following the approach described in Abellán et al. (2017). More specifically, we used the modeled atomic density distributions of 12C, 16O, and 28Si to derive the square root of the products (12C × 16O) and (28Si × 16O). These can be considered as proxies of CO and SiO. In a forthcoming paper (Ono et al., in prep.), we will calculate accurately the molecule formation in SN 1987A using the MHD models described here. Figure 10 and the online Movie 6 show the distributions of (in the following C×O) and (in the following Si×O) at day 9875. Models N16.3 and B18.3 show a dense toroidal structure for both C×O and Si×O at 50% of the peak density located around the center of the explosion. The torus forms because the fast outflows from the bipolar explosion pierce the outer shells of C, O, and Si; then this material is swept out by the outflows and progressively accumulates at their sides, producing the torus-like feature. The torus forms soon after the core-collapse; then it expands following the evolution of the remnant. In model S19.8, this structure is not clearly evident and incomplete shells of C×O and Si×O develop around the center of the explosion.

thumbnail Fig. 10.

Isosurfaces of the distributions of C×O (red) and Si×O (green) at day 9875 from selected view angles for models N16.3 (left panels), B18.3 (center panels), and S19.8 (right panels). The semi-transparent isosurfaces correspond to a value of density which is at 50% of the respective peak density (maximum density is 1.3 × 10−19 g cm−3, 3.2 × 10−19 g cm−3, and 1.9 × 10−19 g cm−3 for Si×O and 9.9 × 10−20 g cm−3, 4.3 × 10−19 g cm−3, and 1.3 × 10−19 g cm−3 for C×O, for models N16.3, B18.3, and S19.8 respectively). The gray arrow lying on the y axis indicate the position of the observer. The ring indicates the position of the dense inner ring. See online Movie 6 for an animation of these data; a navigable 3D graphic of the C×O and Si×O distributions for the three models is available at https://skfb.ly/6P7sI.

Open with DEXTER

The axis of the torus (when present) lies in the plane of the equatorial ring in the direction of the initial large-scale asymmetry: the clumps of early 56Ni, 56Fe, and 44Ti, all travel along this axis. Therefore, the models show that the initial explosion asymmetry required to fit the line profiles of [Fe II] 26 μm and 18 μm few years after the SN predicts a toroidal structure of C×O and Si×O, in models N16.3 and B18.3, which is oriented in the same way as the dense structure of CO and SiO inferred from observations (Abellán et al. 2017). We note that the dense structures of CO and SiO observed with ALMA are quite irregular and their torus-like shape (as suggested by observations) relies just on few plumes of dense CO and SiO. Nevertheless, the observations show that these plumes lie on a preferential plane whose orientation coincides with that of the torus-like feature of models N16.3 and B18.3.

Figure 11 shows the distributions of C×O and Si×O for isosurfaces corresponding to densities at 30% and 70% of the peak. At 30% of the peak density, the central deficit of density is fully enclosed into a shell of C×O and Si×O in all models. At 70% of the peak, the torus-like structure evident at 50% of the peak in models N16.3 and B18.3 is still clearly visible, whilst an irregular shell characterizes the distribution of C×O and Si×O in model S19.8.

thumbnail Fig. 11.

As in Fig. 10 but for semi-transparent isosurfaces corresponding to a value of density which is at 30% (upper panels) and 70% (lower panels) of the respective peak.

Open with DEXTER

A more quantitative comparison between the modeled and the observed structures can be done by comparing the typical size scales inferred from observations and obtained with the three models around day ≈10 000. Figure 12 shows maps of C×O and Si×O derived by integrating their density distributions along the axis of the toroidal structure at day 9875 for the three models. The modeled distributions are more regular than those observed; in particular both C×O and Si×O exhibit clear and regular ring-like distributions at odds with the observations that show more clumpy and irregular ring-like morphologies for CO and SiO distributions. This is most likely due to the idealized initial asymmetric explosion prescribed in our simulations. Nevertheless, the large-scale morphology of the modeled structures and their typical length scales are very similar to those observed.

thumbnail Fig. 12.

Maps of C×O (upper panels) and Si×O (lower panels) integrated along the axis of the torus-like structure at day 9875 for models N16.3, B18.3, and S19.8. The yardstick in the upper panels indicates the length scale in cm and in arcsec, assuming a distance of 51.4 kpc.

Open with DEXTER

From Fig. 12 we note that, in all models, C×O presents a maximum extension larger than Si×O as also found in the observed distributions of CO and SiO (Abellán et al. 2017). This is more clear by inspecting Fig. 13 which shows the angle-averaged radial distributions of C×O and Si×O derived from the maps in Fig. 12. In all the models, C×O and Si×O are spatially distinct as in the observations. Models N16.3 and B18.3 predict a size of C×O and Si×O (see Table 4) of the same order of magnitude of the observed maximum extension of CO and SiO (≈ ± 1.53 × 1017 cm and ≈ ± 1.17 × 1017 cm, respectively; Abellán et al. 2017). Model S19.8 predicts a size of C×O and Si×O significantly larger than observed.

thumbnail Fig. 13.

Average radial distribution of C×O (red curve) and Si×O (green) from the axis of the torus-like structure at day 9875 for models N16.3, B18.3, and S19.8. Symbols indicate the extension of the C×O (red stars) and Si×O (green stars) structures at 50% of the peak of the profiles; vertical dashed lines show the maximum extension of CO (red) and SiO (green) inferred from ALMA observations (Abellán et al. 2017).

Open with DEXTER

Table 4.

Maximum extension of CO and SiO inferred from the observations and derived from the models.

Despite the idealization of the initial asymmetry of the SN explosion, we found that in models N16.3 and B18.3 the asymmetry and its orientation, which are necessary to reproduce the redshift and broadening of observed [Fe II] lines during the first two years after the SN (Haas et al. 1990) and the redshift of 44Ti observed at day ≈10 000, naturally produce a toroidal dense structure in the distributions of C×O and Si×O (proxies of CO and SiO molecules) at day ≈10 000 with position, orientation, and size similar to those of observed CO and SiO distributions (Abellán et al. 2017). We conclude that all ejecta anisotropies discussed above originate from the same large-scale asymmetry in the SN explosion.

3.5. Neutron star kick

The highly asymmetric explosion is expected to cause momentum transfer to a newly formed compact object possibly produced by SN 1987A (e.g., Wongwathanarat et al. 2013): the proto-neutron star and the center of gravity of the ejecta are expected to move in opposite directions due to momentum conservation. In the case of the initial asymmetry required by our simulations to reproduce the observations of the SN and the SNR, the kick is mainly regulated by the parameter α of the asymmetry which also determines the shift of the centroid of [Fe II] lines about 1–2 years after the SN event. The parameter is well constrained by the comparison between the modeled and the observed shifts (see Sect. 3.1). Thus, the kick was estimated by assuming momentum conservation as in Wongwathanarat et al. (2013): it is calculated by imposing the initial total momentum equal zero (see also Paper I). We found that all the three cases investigated predict that the compact object would have received a kick velocity of ≈300 km s−1 toward the observer in the north (see the lightblue arrow in Fig. 3; see also Paper I).

We note that the initial asymmetry may have been more complex than that modeled here, possibly producing plumes and fingers in different directions. In fact, as discussed in Sect. 3.1, our models do not describe, for instance, the isolated fast clump of iron responsible for the spectral feature redshifted at v ≈ 4000 km s−1 in [Fe II] lines (see upper right panel of Fig. 4). That feature may be a signature of a structure of the initial asymmetry more complex than modeled here. One may ask, for instance, if the explosion was not strictly bipolar but had dominant outflows separated by less than 180°. In our models, we considered simple idealized forms for the asymmetry, with a bipolar structure. The fact that one of our models is able to broadly reproduce the observed line profiles suggests that two dominant outflows were produced by the explosion, a prominent one redshifted and the other blueshifted, separated by an angle which was not much lower than 180°. Other smaller scale structures may have been produced at the collapse, but our simulations show that the bulk of the explosion kinematics was dominated by the bipolar outflows. Nevertheless, we note that the predicted kick velocity might be a lower limit to the actual velocity. In fact, the high velocity spectral feature in the red wing of [Fe II] lines suggests that the bulk of ejecta (not only iron) might present more mass traveling away from the observer. In this case, our models might underestimate the momentum of ejecta traveling away from the observer and, consequently, the momentum of the compact object traveling toward the observer.

Our estimate for the neutron star kick is in good agreement with previous studies (Nagataki 2000; Wongwathanarat et al. 2013; Janka et al. 2017). Early studies pointed out that a highly asymmetric explosion required for SN 1987A would have produced a kick velocity of the same order of magnitude as found here (Nagataki 2000), although the authors suggested that the proto-neutron star should move toward the observer in the south because they have assumed that the initial asymmetry was aligned with the axis of the equatorial ring. More recently, the comparison of the iron distribution derived in a 3D neutrino-driven explosion model for a 15 M star appropriate for SN 1987A (Janka et al. 2017) with 3D maps of silicon and iron derived from spectral and imaging observations of SN 1987A with HST-STIS and VLT-SINFONI (Larsson et al. 2016) suggested a kick velocity of the same order of magnitude toward the observer in the north, but with an angle with respect to the LoS slightly smaller than in our models (compare Fig. 4 in Janka et al. 2017 with Fig. 3 in this paper).

Finally, our model predictions are also in good agreement with recent findings from the analysis of ALMA observations of dust in the SN 1987A ejecta (Cigan et al. 2019). The analysis has shown a dust blob in ALMA images which is slightly offset along the north-east direction from the estimated position of Sk −69°  202. Cigan et al. (2019) have interpreted the blob as dust heated by the compact object and estimated its kick velocity to be ≈700 km s−1. The estimate of the kick was done by considering the age of SN 1987A and the offset between the brightest part of the blob (interpreted as an early pulsar wind nebula) and the position of Sk −69°  202. However, in the Crab Nebula, the position of the compact object and that of the brightest part of the pulsar wind nebula are not coincident (Weisskopf et al. 2000). Thus, some uncertainty in the determination of the offset and of the kick velocity should be considered. Nevertheless, the predicted position of the neutron star and the lower limit on its kick velocity derived with our favorite model B18.3 are both fully consistent with the interpretation proposed by Cigan et al. (2019).

4. Summary and conclusions

We investigated the evolution of the remnant of SN 1987A with the aim of unveiling the link between the radiative and dynamical properties of the remnant and the properties of the parent SN explosion and structure of the progenitor star. To this end, we developed a comprehensive 3D hydrodynamic model which describes the long-term evolution of SN 1987A from the onset of the SN to the development of its remnant, accounting for the pre-SN structure of the progenitor star. The model was realized by coupling a 3D model of a core-collapse SN (Ono et al. 2013 and Paper I) with a 3D model of a SNR (Orlando et al. 2015, 2019a).

The SN model has the relevant physics required to describe the evolution in the immediate aftermath of the core-collapse (see Sect. 2.2): the explosive nucleosynthesis through a nuclear reaction network, the feedback of nuclear energy generation, the energy deposition due to radioactive decays of isotopes synthesized in the SN, the effects of gravity (both self-gravity and gravitational effects of the central proto-neutron star), and the fallback of material on the proto-neutron star. Aspherical explosions (as suggested by observations; Larsson et al. 2016) were simulated to mimic the global anisotropies developing in jet-driven SNe, or by convection in the neutrino heating layers and the SASI (Paper I).

The 3D simulations of the SN were initialized from pre-SN stellar models available in literature (see Sect. 2.1): a 16.3 M BSG (model N16.3; Shigeyama & Nomoto 1990), a 18.3 M BSG (B18.3; Urushibata et al. 2018), and a 15.9 M RSG (S19.8; Sukhbold et al. 2016). Model B18.3 results from the merging of two massive stars and reproduces most of the observational constraints of the progenitor star of SN 1987A: the red-to-blue evolution, its lifetime, the total mass and the position in the Hertzsprung–Russell diagram at collapse (Urushibata et al. 2018).

We then used the output from the SN simulations (about 20 h after the core-collapse) to start 3D MHD simulations, describing the subsequent expansion of the ejecta through the inhomogeneous pre-SN environment. The SNR model includes the relevant physics to describe the evolution of the remnant and its interaction with the inhomogeneous ambient environment (see Sect. 2.3): a plausible configuration of the ambient magnetic field, the deviations from temperature-equilibration between electrons and ions, and the deviations from equilibrium of ionization of the most abundant ions. The inhomogeneous pre-SN ambient environment was described as a nebula consisting of a dense ring surrounded by a H II region with 3D geometry compatible with observations (e.g., Sugerman et al. 2005; Smith et al. 2013). The simulations cover 50 years, thus describing the evolution of SN 1987A since the SN and providing some hints on its evolution until 2037. Current and future observations will help in unveiling the structure and geometry of the medium beyond the dense equatorial ring (e.g., Larsson et al. 2019) and will provide important constraints for the models.

We explored different explosion energies and asymmetries, and different orientations of the asymmetry with respect to the LoS, to identify the physical and geometrical properties of the SN explosion which led to the shifts and broadening of [Fe II] lines in early observations (e.g., Haas et al. 1990; Colgan et al. 1994; Utrobin et al. 1995). Furthermore, we explored several density structures of the nebula around the parameters found in previous works (Orlando et al. 2015, 2019a) and within the range of values inferred from observations (e.g., Dewey et al. 2012; Sugerman et al. 2005), searching for parameters which best reproduce X-ray lightcurves and morphology of the remnant observed in the last 30 years (e.g., Haberl et al. 2006; Maggi et al. 2012; Helder et al. 2013; Orlando et al. 2015; Frank et al. 2016), and the observed late distribution of metal rich ejecta (e.g., Boggs et al. 2015; Abellán et al. 2017).

The best match with observations was found for models with a highly asymmetric SN explosion with a total energy at the shock breakout of ≈2 × 1051 erg (see Table 1). The energy is slightly higher than that estimated in a previous work (Orlando et al. 2015), but still within the range of values expected for SN 1987A. For all the progenitor stars investigated, the observations require that a significant fraction of energy was channeled along an axis almost lying in the plane of the central ring, roughly in the direction of the LoS but with an offset of 40° (see Fig. 3); the projection of this axis in the plane of the sky is offset by 15° from the northsouth axis. Such large asymmetries may result from a combination of low-mode convection in the neutrino heating layers and strong axial sloshing of the shock due to the SASI (e.g., Blondin et al. 2003); alternatively, extreme asymmetries may develop in magnetic, jet-driven SNe (e.g., Khokhlov et al. 1999) or in single-lobe SNe (e.g., Hungerford et al. 2005).

We found that the structure of the progenitor star and the initial asymmetry of the SN both play a central role in determining the mixing of 56Ni to speeds exceeding 3000 km s−1, and the observed high redshifts of [Fe II] lines around day 400 and 44Ti lines around day 10 000. Among the models investigated, we found that B18.3 is the only one satisfying the observational constraints for the progenitor star, the SN, and the SNR (see Table 5). In fact, this model is the only one able to reproduce the mixing of ejecta which leads to the observed velocity distributions of 56Ni and 56Fe. This is possible because the C+O and He layers in this model are significantly less extended and less dense than in the others. As a result, nickel (and heavy elements formed in the innermost layers of the remnant) can penetrate efficiently through the helium and hydrogen layers producing the observed mixing and velocities.

Table 5.

Comparison of model results and observations.

The imprint of the progenitor star is also evident in the X-ray observations of the SNR. The X-ray lightcurves of SN 1987A can be reproduced by the two BSG (N16.3 and B18.3) but not by the RSG progenitor (S19.8). This is due to the structure of the extended hydrogen envelope of the RSG which causes a significant deceleration of the blast wave before the shock breakout, so that, when the blast travels through the CSM, it has a velocity much smaller than in the other two models examined. Furthermore, the differences in the structure of the progenitor stars make B18.3 the model best reproducing also the size and morphology of the X-ray emitting remnant during its evolution, whilst the other two models fail in reproducing at the same time fluxes and size of the X-ray source. We conclude that the dynamical and radiative properties of the remnant require a progenitor BSG that resulted from the merging of two massive stars before the collapse.

The observed distribution of cold molecular gas (e.g., CO and SiO) is also consistent with the same progenitor star and explosion asymmetry needed to reproduce the high redshifts and broadenings of [Fe II] and 44Ti lines, and the X-ray lightcurves and morphology of the remnant. This suggests that all these features are signatures of the same process associated with the SN.

As a result of the highly asymmetric explosion, we estimated that a compact object possibly produced by SN 1987A would have received a kick velocity of ≈300 km s−1 toward the observer in the north (lightblue arrow in Fig. 3; see Sect. 3.5). Currently the newly formed neutron star (if any) has not been detected yet. This is probably due to heavy absorption of its emission caused by a high local absorbing column density along the LoS (Orlando et al. 2015; Alp et al. 2018). Future detection of the neutron star in SN 1987A will provide a strong constraint on the initial large-scale asymmetry and on the mechanism of neutron star kick.

Our long-term simulations provided a way to link the remnant of SN 1987A to its parent SN explosion and to the internal structure of its progenitor star. Extending the present work to other young SNRs will permit advances in resolving pending questions about the physical processes associated with SNe and the final stages of stellar evolution.


1

This approximation allowed us to reduce the computational cost needed if, instead, the Poisson equation for self-gravity is solved.

Acknowledgments

We thank the anonymous referee for useful suggestions that have allowed us to improve the paper. We acknowledge that the results of this research have been achieved using the PRACE Research Infrastructure resource Marconi based in Italy at CINECA (PRACE Award N.2016153460). Additional numerical computations were carried out complementarily on XC40 (YITP, Kyoto University), Cray XC50 (Center for Computational Astrophysics, National Astronomical Observatory of Japan), HOKUSAI (RIKEN). This research also used computational resources of the K computer provided by the RIKEN Center for Computational Science through the HPCI System Research project (Project ID:hp180281). The FLASH code, used in this work, is developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago (USA). The PLUTO code is developed at the Turin Astronomical Observatory (Italy) in collaboration with the Department of General Physics of Turin University (Italy) and the SCAI Department of CINECA (Italy). SO would like to thank H.-T. Janka for helpfull discussions on this work and A. Mignone for his support with the PLUTO code. SO, MM, GP, FB acknowledge financial contribution from the agreement ASI-INAF n.2017-14-H.O, and partial financial support by the PRIN INAF 2016 grant “Probing particle acceleration and γ-ray propagation with CTA and its precursors”. This work is supported by JSPS Grants-in-Aid for Scientific Research “KAKENHI” Grant Numbers JP26800141 and JP19H00693. SN and MO wish to acknowledge the support from the Program of Interdisciplinary Theoretical & Mathematical Sciences (iTHEMS) at RIKEN (Japan). SN also acknowledges the support from Pioneering Program of RIKEN for Evolution of Matter in the Universe (r-EMU). The navigable 3D graphics have been developed in the framework of the project 3DMAP-VR (3-Dimensional Modeling of Astrophysical Phenomena in Virtual Reality; Orlando et al. 2019b) at INAF-Osservatorio astronomico di Palermo (http://cerere.astropa.unipa.it/progetti_ricerca/HPC/3dmap_vr.htm) and uploaded on the Sketchfab platform.

References

  1. Abellán, F. J., Indebetouw, R., Marcaide, J. M., et al. 2017, ApJ, 842, L24 [NASA ADS] [CrossRef] [Google Scholar]
  2. Alp, D., Larsson, J., Fransson, C., et al. 2018, ApJ, 864, 174 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bader, G., & Deuflhard, P. 1983, NuMat, 41, 373 [Google Scholar]
  4. Bear, E., & Soker, N. 2018, MNRAS, 478, 682 [NASA ADS] [Google Scholar]
  5. Blondin, J. M., Mezzacappa, A., & DeMarino, C. 2003, ApJ, 584, 971 [NASA ADS] [CrossRef] [Google Scholar]
  6. Boggs, S. E., Harrison, F. A., Miyasaka, H., et al. 2015, Science, 348, 670 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  7. Chevalier, R. A., & Dwarkadas, V. V. 1995, ApJ, 452, L45 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cigan, P., Matsuura, M., Gomez, H. L., et al. 2019, ApJ, 886, 51 [NASA ADS] [CrossRef] [Google Scholar]
  9. Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  10. Colgan, S. W. J., & Hollenbach, D. J. 1988, ApJ, 329, L25 [NASA ADS] [CrossRef] [Google Scholar]
  11. Colgan, S. W. J., Haas, M. R., Erickson, E. F., Lord, S. D., & Hollenbach, D. J. 1994, ApJ, 427, 874 [NASA ADS] [CrossRef] [Google Scholar]
  12. Couch, S. M., Wheeler, J. C., & Milosavljević, M. 2009, ApJ, 696, 953 [NASA ADS] [CrossRef] [Google Scholar]
  13. Dewey, D., Dwarkadas, V. V., Haberl, F., Sturm, R., & Canizares, C. R. 2012, ApJ, 752, 103 [NASA ADS] [CrossRef] [Google Scholar]
  14. Donati, J.-F., & Landstreet, J. D. 2009, ARA&A, 47, 333 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  15. Duff, I. S., Erisman, A. M., & Reid, J. K. 1986, Direct Methods for Sparse Matrices (Oxford: Clarendon) [Google Scholar]
  16. Ferrand, G., Warren, D. C., Ono, M., et al. 2019, ApJ, 877, 136 [NASA ADS] [CrossRef] [Google Scholar]
  17. Frank, K. A., Zhekov, S. A., Park, S., et al. 2016, ApJ, 829, 40 [NASA ADS] [CrossRef] [Google Scholar]
  18. Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
  19. Ghavamian, P., Laming, J. M., & Rakowski, C. E. 2007, ApJ, 654, L69 [NASA ADS] [CrossRef] [Google Scholar]
  20. Haas, M. R., Colgan, S. W. J., Erickson, E. F., et al. 1990, ApJ, 360, 257 [NASA ADS] [CrossRef] [Google Scholar]
  21. Haberl, F., Geppert, U., Aschenbach, B., & Hasinger, G. 2006, A&A, 460, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Handy, T., Plewa, T., & Odrzywołek, A. 2014, ApJ, 783, 125 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hanuschik, R. W., Thimm, G., & Dachs, J. 1988, MNRAS, 234, 41P [NASA ADS] [CrossRef] [Google Scholar]
  24. Helder, E. A., Broos, P. S., Dewey, D., et al. 2013, ApJ, 764, 11 [NASA ADS] [CrossRef] [Google Scholar]
  25. Hillebrandt, W., Hoeflich, P., Weiss, A., & Truran, J. W. 1987, Nature, 327, 597 [NASA ADS] [CrossRef] [Google Scholar]
  26. Hungerford, A. L., Fryer, C. L., & Rockefeller, G. 2005, ApJ, 635, 487 [NASA ADS] [CrossRef] [Google Scholar]
  27. Ivanova, N., Podsiadlowski, P., & Spruit, H. 2002, MNRAS, 334, 819 [NASA ADS] [CrossRef] [Google Scholar]
  28. Janka, H. T. 2017, Neutrino-Driven Explosions (Cham: Springer International Publishing AG), 1095 [Google Scholar]
  29. Janka, H.-T., Melson, T., & Summa, A. 2016, Ann. Rev. Nucl. Part. Sci., 66, 341 [NASA ADS] [CrossRef] [Google Scholar]
  30. Janka, H. T., Gabler, M., & Wongwathanarat, A. 2017, in Supernova 1987A:30 years later - Cosmic Rays and Nuclei from Supernovae and their Aftermaths, eds. A. Marcowith, M. Renaud, G. Dubner, A. Ray, & A. Bykov, IAU Symp., 331, 148 [NASA ADS] [Google Scholar]
  31. Joggerst, C. C., Woosley, S. E., & Heger, A. 2009, ApJ, 693, 1780 [NASA ADS] [CrossRef] [Google Scholar]
  32. Joggerst, C. C., Almgren, A., Bell, J., et al. 2010, ApJ, 709, 11 [NASA ADS] [CrossRef] [Google Scholar]
  33. Khokhlov, A. M., Höflich, P. A., Oran, E. S., et al. 1999, ApJ, 524, L107 [NASA ADS] [CrossRef] [Google Scholar]
  34. Kifonidis, K., Plewa, T., Scheck, L., Janka, H.-T., & Müller, E. 2006, A&A, 453, 661 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Larsson, J., Fransson, C., Kjaer, K., et al. 2013, ApJ, 768, 89 [NASA ADS] [CrossRef] [Google Scholar]
  36. Larsson, J., Fransson, C., Spyromilio, J., et al. 2016, ApJ, 833, 147 [NASA ADS] [CrossRef] [Google Scholar]
  37. Larsson, J., Fransson, C., Alp, D., et al. 2019, ApJ, 886, 147 [NASA ADS] [CrossRef] [Google Scholar]
  38. MacNeice, P., Olson, K. M., Mobarry, C., de Fainchtein, R., & Packer, C. 2000, Comput. Phys. Commun., 126, 330 [NASA ADS] [CrossRef] [Google Scholar]
  39. Maeda, K., & Nomoto, K. 2003, ApJ, 598, 1163 [NASA ADS] [CrossRef] [Google Scholar]
  40. Maggi, P., Haberl, F., Sturm, R., & Dewey, D. 2012, A&A, 548, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Mao, J., Ono, M., Nagataki, S., et al. 2015, ApJ, 808, 164 [NASA ADS] [CrossRef] [Google Scholar]
  42. Marek, A., & Janka, H.-T. 2009, ApJ, 694, 664 [NASA ADS] [CrossRef] [Google Scholar]
  43. Menon, A., & Heger, A. 2017, MNRAS, 469, 4649 [NASA ADS] [Google Scholar]
  44. Menon, A., Utrobin, V., & Heger, A. 2019, MNRAS, 482, 438 [NASA ADS] [CrossRef] [Google Scholar]
  45. Miceli, M., Orlando, S., Burrows, D. N., et al. 2019, Nat. Astron., 3, 236 [NASA ADS] [CrossRef] [Google Scholar]
  46. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
  47. Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7 [NASA ADS] [CrossRef] [Google Scholar]
  48. Morris, T., & Podsiadlowski, P. 2007, Science, 315, 1103 [NASA ADS] [CrossRef] [Google Scholar]
  49. Nagataki, S. 2000, ApJS, 127, 141 [NASA ADS] [CrossRef] [Google Scholar]
  50. Nagataki, S., Hashimoto, M., Sato, K., & Yamada, S. 1997, ApJ, 486, 1026 [NASA ADS] [CrossRef] [Google Scholar]
  51. Nomoto, K., & Hashimoto, M. 1988, Phys. Rep., 163, 13 [NASA ADS] [CrossRef] [Google Scholar]
  52. Nordhaus, J., Burrows, A., Almgren, A., & Bell, J. 2010, ApJ, 720, 694 [NASA ADS] [CrossRef] [Google Scholar]
  53. Ono, M., Nagataki, S., Ito, H., et al. 2013, ApJ, 773, 161 [NASA ADS] [CrossRef] [Google Scholar]
  54. Ono, M., Nagataki, S., Ferrand, G., et al. 2020, ApJ, 888, 111 [NASA ADS] [CrossRef] [Google Scholar]
  55. Orlando, S., Bocchino, F., Reale, F., Peres, G., & Pagano, P. 2008, ApJ, 678, 274 [NASA ADS] [CrossRef] [Google Scholar]
  56. Orlando, S., Miceli, M., Pumo, M. L., & Bocchino, F. 2015, ApJ, 810, 168 [NASA ADS] [CrossRef] [Google Scholar]
  57. Orlando, S., Miceli, M., Pumo, M. L., & Bocchino, F. 2016, ApJ, 822, 22 [NASA ADS] [CrossRef] [Google Scholar]
  58. Orlando, S., Miceli, M., Petruk, O., et al. 2019a, A&A, 622, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Orlando, S., Pillitteri, I., Bocchino, F., Daricello, L., & Leonardi, L. 2019b, Res. Am. Astron. Soc., 3, 176 [NASA ADS] [Google Scholar]
  60. Panagia, N. 1999, in New Views of the Magellanic Clouds, eds. Y. H. Chu, N. Suntzeff, J. Hesser, & D. Bohlender, IAU Symp., 190, 549 [NASA ADS] [Google Scholar]
  61. Park, S., Zhekov, S. A., Burrows, D. N., et al. 2006, ApJ, 646, 1001 [NASA ADS] [CrossRef] [Google Scholar]
  62. Parker, E. N. 1958, ApJ, 128, 664 [Google Scholar]
  63. Patnaude, D. J., Lee, S.-H., Slane, P. O., et al. 2015, ApJ, 803, 101 [NASA ADS] [CrossRef] [Google Scholar]
  64. Phillips, M. M., & Heathcote, S. R. 1989, PASP, 101, 137 [NASA ADS] [CrossRef] [Google Scholar]
  65. Quirk, J. J. 1997, in Upwind and High-Resolution Schemes, eds. M. Yousuff Hussaini, B. van Leer, & J. Van Rosendale (Berlin: Springer), 550 [Google Scholar]
  66. Seitenzahl, I. R., Timmes, F. X., & Magkotsios, G. 2014, ApJ, 792, 10 [NASA ADS] [CrossRef] [Google Scholar]
  67. Shigeyama, T., & Nomoto, K. 1990, ApJ, 360, 242 [NASA ADS] [CrossRef] [Google Scholar]
  68. Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91 [NASA ADS] [CrossRef] [Google Scholar]
  69. Smith, N., Arnett, W. D., Bally, J., Ginsburg, A., & Filippenko, A. V. 2013, MNRAS, 429, 1324 [NASA ADS] [CrossRef] [Google Scholar]
  70. Sugerman, B. E. K., Crotts, A. P. S., Kunkel, W. E., Heathcote, S. R., & Lawrence, S. S. 2005, ApJS, 159, 60 [NASA ADS] [CrossRef] [Google Scholar]
  71. Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H.-T. 2016, ApJ, 821, 38 [NASA ADS] [CrossRef] [Google Scholar]
  72. Suwa, Y., Kotake, K., Takiwaki, T., et al. 2010, PASJ, 62, L49 [NASA ADS] [CrossRef] [Google Scholar]
  73. Takiwaki, T., Kotake, K., & Suwa, Y. 2012, ApJ, 749, 98 [NASA ADS] [CrossRef] [Google Scholar]
  74. Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 [NASA ADS] [CrossRef] [Google Scholar]
  75. Urushibata, T., Takahashi, K., Umeda, H., & Yoshida, T. 2018, MNRAS, 473, L101 [NASA ADS] [CrossRef] [Google Scholar]
  76. Utrobin, V. P., Chugai, N. N., & Andronova, A. A. 1995, A&A, 295, 129 [NASA ADS] [Google Scholar]
  77. Utrobin, V. P., Wongwathanarat, A., Janka, H. T., et al. 2019, A&A, 624, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Walk, L., Tamborra, I., Janka, H.-T., & Summa, A. 2018, Phys. Rev. D, 98, 123001 [NASA ADS] [CrossRef] [Google Scholar]
  79. Wang, L., Wheeler, J. C., Höflich, P., et al. 2002, ApJ, 579, 671 [NASA ADS] [CrossRef] [Google Scholar]
  80. Weaver, T. A., Zimmerman, G. B., & Woosley, S. E. 1978, ApJ, 225, 1021 [NASA ADS] [CrossRef] [Google Scholar]
  81. Weisskopf, M. C., Hester, J. J., Tennant, A. F., et al. 2000, ApJ, 536, L81 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  82. Wongwathanarat, A., Janka, H. T., & Müller, E. 2013, A&A, 552, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Wongwathanarat, A., Müller, E., & Janka, H.-T. 2015, A&A, 577, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Wongwathanarat, A., Janka, H.-T., Müller, E., Pllumbi, E., & Wanajo, S. 2017, ApJ, 842, 13 [NASA ADS] [CrossRef] [Google Scholar]
  85. Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Modern Phys., 74, 1015 [NASA ADS] [CrossRef] [Google Scholar]
  86. Zanardo, G., Staveley-Smith, L., Gaensler, B. M., et al. 2018, ApJ, 861, L9 [NASA ADS] [CrossRef] [Google Scholar]
  87. Zhekov, S. A., McCray, R., Dewey, D., et al. 2009, ApJ, 692, 1190 [NASA ADS] [CrossRef] [Google Scholar]

Movies

Movie of Fig. 4 (Access here)

Movie of Fig. 7 (Access here)

Movie of Fig. 7 (Access here)

Movie of Fig. 7 (Access here)

Movie of Fig. 9 (Access here)

Movie of Fig. 10 (Access here)

All Tables

Table 1.

Explored range of parameters of the SN models, and parameters of the best-fit models considered in this paper.

Table 2.

Explored and best-fit parameters of the CSM for the models of SN 1987A.

Table 3.

Parameters of the Gaussian functions best fitting the observed and synthetic lines of 44Ti.

Table 4.

Maximum extension of CO and SiO inferred from the observations and derived from the models.

Table 5.

Comparison of model results and observations.

All Figures

thumbnail Fig. 1.

Radial profiles of ρr3 of the pre-SN models considered in this paper (upper left panel) and mass fractions of selected elements (see legend in each panel) versus the interior mass Mr in the pre-SN models N16.3 (upper right; Nomoto & Hashimoto 1988; Shigeyama & Nomoto 1990), B18.3 (lower left; Urushibata et al. 2018), and S19.8 (lower right; Woosley et al. 2002; Sukhbold et al. 2016).

Open with DEXTER
In the text
thumbnail Fig. 2.

Angle dependence of the initial radial velocity in the [x, z] plane of the Cartesian coordinate, adopting an elliptical form for the function f(θ) in Eq. (1). The different curves are for different values of β = vpol/veq and α = vup/vdw = 1.5.

Open with DEXTER
In the text
thumbnail Fig. 3.

Schematic view of the orientation of the initial asymmetry (lightblue and red arrows) with respect to the dense inner ring (lightblue circle) and to the LoS (gray arrow along negative y-axis). The negative y-axis points toward the observer. The yardstick indicating the length scale is on the left.

Open with DEXTER
In the text
thumbnail Fig. 4.

Isosurfaces of the iron distribution for models N16.3 (upper left panel), B18.3 (lower left), and S19.8 (lower right) at day 730. The semi-transparent isosurfaces correspond to a value of the iron density which is at 10% of the peak density for each model (maximum iron density is 4.9 × 10−16 g cm−3 for N16.3, 6.0 × 10−16 g cm−3 for B18.3, and 1.1 × 10−15 g cm−3 for S19.8). The colors give the velocity along the LoS in units of 1000 km s−1 on the isosurface; the color coding is defined at the bottom of each panel. The yardstick indicating the length scale is on the left. The upper right panel shows the line profiles of [Fe II] 26 μm and 18 μm of SN 1987A observed during the first 2 years of evolution (symbols; Haas et al. 1990) and synthesized from models N16.3, B18.3, and S19.8 at day 730 (solid lines). See online Movie 1 for an animation of these data; a navigable 3D graphic of the iron distribution for the three models is available at https://skfb.ly/6P7sH.

Open with DEXTER
In the text
thumbnail Fig. 5.

Mass distributions of 1H, 4He, 56Fe, 16O, 28Si, and 44Ti versus radial velocity vrad (left panels) and velocity component along the LoS, vLoS (right panels) at day 730 for models N16.3 (top), B18.3 (middle), and S19.8 (bottom).

Open with DEXTER
In the text
thumbnail Fig. 6.

X-ray lightcurves (solid lines) in the [0.5, 2] keV (upper panels) and [3, 10] keV (lower panels) bands synthesized from models N16.3 (left panels), B18.3 (center panels), and S19.8 (right panels) compared to the lightcurves of SN 1987A observed with Rosat (black triangles; Haberl et al. 2006), ASCA (green diamonds; Orlando et al. 2015), Chandra (magenta circles; Helder et al. 2013; Frank et al. 2016) and XMM-Newton (cyan circles; Haberl et al. 2006; Maggi et al. 2012). Dotted, dashed, and dot-dashed lines indicate the contribution to emission from the shocked ejecta, the shocked plasma from the H II region, and the shocked plasma from the ring, respectively.

Open with DEXTER
In the text
thumbnail Fig. 7.

Synthetic X-ray emission maps from models N16.3 (top), B18.3 (middle), and S19.8 (bottom) at the labeled times. Each panel consists of 4 quadrants: emission maps in the [0.5, 2] keV (upper left) and [3, 10] keV (upper right) bands, and corresponding three-color composite representations (lower left and right). Each image is normalized to its maximum for visibility; the X-ray maps were convolved with a Gaussian of size 0.15 arcsec to approximate the spatial resolution of Chandra; the three-color composite images were smoothed with a Gaussian of 0.025 arcsec; the colors in the composite show the emission from shocked ejecta (green), shocked plasma from the ring (red), and from the H II region (blue). See online Movies 2, 3, and 4 for animations of these data.

Open with DEXTER
In the text
thumbnail Fig. 8.

Best-fit radius of the X-ray emitting torus vs time inferred from the analysis of actual X-ray observations (symbols; Frank et al. 2016) and that derived from the analysis of synthetic X-ray images for models B18.3 (red solid line), N16.3 (blue), and S19.8 (green).

Open with DEXTER
In the text
thumbnail Fig. 9.

Isosurfaces of the titanium distribution for models N16.3 (upper left), B18.3 (lower left), and S19.8 (lower right) at day 9875. The semi-transparent isosurface corresponds to a value of the titanium density which is at 10% of the peak (maximum titanium density is 8.6 × 10−22 g cm−3 for N16.3, 6.9 × 10−22 g cm−3 for B18.3, and 2.3 × 10−21 g cm−3 for S19.8). The colors give the velocity along the LoS in units of 1000 km s−1 on the isosurface; the color coding is defined at the bottom of each panel. The ring indicates the position of the dense inner ring. The yardstick on the left indicates the length scale. Upper right panel: best-fit Gaussian for the 44Ti lines observed between 2012 and 2014 (black line; Boggs et al. 2015) and lines profile synthesized from models N16.3 (blue line), B18.3 (red line), and S19.8 (green line) at day 9875. The shaded region represents the 90% confidence area for the position of the peak and for the width of the Gaussian. See online Movie 5 for an animation of these data; a navigable 3D graphic of the titanium distribution for the three models is available at https://skfb.ly/6P7sI.

Open with DEXTER
In the text
thumbnail Fig. 10.

Isosurfaces of the distributions of C×O (red) and Si×O (green) at day 9875 from selected view angles for models N16.3 (left panels), B18.3 (center panels), and S19.8 (right panels). The semi-transparent isosurfaces correspond to a value of density which is at 50% of the respective peak density (maximum density is 1.3 × 10−19 g cm−3, 3.2 × 10−19 g cm−3, and 1.9 × 10−19 g cm−3 for Si×O and 9.9 × 10−20 g cm−3, 4.3 × 10−19 g cm−3, and 1.3 × 10−19 g cm−3 for C×O, for models N16.3, B18.3, and S19.8 respectively). The gray arrow lying on the y axis indicate the position of the observer. The ring indicates the position of the dense inner ring. See online Movie 6 for an animation of these data; a navigable 3D graphic of the C×O and Si×O distributions for the three models is available at https://skfb.ly/6P7sI.

Open with DEXTER
In the text
thumbnail Fig. 11.

As in Fig. 10 but for semi-transparent isosurfaces corresponding to a value of density which is at 30% (upper panels) and 70% (lower panels) of the respective peak.

Open with DEXTER
In the text
thumbnail Fig. 12.

Maps of C×O (upper panels) and Si×O (lower panels) integrated along the axis of the torus-like structure at day 9875 for models N16.3, B18.3, and S19.8. The yardstick in the upper panels indicates the length scale in cm and in arcsec, assuming a distance of 51.4 kpc.

Open with DEXTER
In the text
thumbnail Fig. 13.

Average radial distribution of C×O (red curve) and Si×O (green) from the axis of the torus-like structure at day 9875 for models N16.3, B18.3, and S19.8. Symbols indicate the extension of the C×O (red stars) and Si×O (green stars) structures at 50% of the peak of the profiles; vertical dashed lines show the maximum extension of CO (red) and SiO (green) inferred from ALMA observations (Abellán et al. 2017).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.