Emission lines due to ionizing radiation from a compact object in the remnant
of Supernova 1987A

C. Fransson,1∗ M. J. Barlow,2† P. J. Kavanagh,3,4† J. Larsson,5† O. C. Jones,6

B. Sargent, 7,8 M. Meixner,9 P. Bouchet,10 T. Temim,11 G. S. Wright,6

J. A. D. L. Blommaert,12 N. Habel,9 A. S. Hirschauer,7 J. Hjorth, 13 L.
Lenkić,14

T. Tikkanen,15 R. Wesson,16 A. Coulais,17,10 O. D. Fox,7 R. Gastaud,18

A. Glasse,6 J. Jaspers,3,4 O. Krause,19 R. M. Lau,20 O. Nayak,21

A. Rest,7,8 L. Colina,22 E. F. van Dishoeck,23,24 M. Güdel25,19,26 Th.
Henning,19

P.-O. Lagage,10 G. Östlin,1 T. P. Ray,4 B. Vandenbussche27

1Department of Astronomy, Stockholm University, The Oskar Klein Centre,
AlbaNova, SE-106 91 Stockholm, Sweden
2Department of Physics and Astronomy, University College London, London WC1E
6BT, UK
3Department of Experimental Physics, Maynooth University, Maynooth, Co
Kildare, Ireland
4Astronomy & Astrophyics Section, School of Cosmic Physics, Dublin Institute
for Advanced Studies, Dublin 2, Ireland
5Department of Physics, KTH Royal Institute of Technology, The Oskar Klein
Centre, AlbaNova, SE-106 91 Stockholm, Sweden
6 UK Astronomy Technology Centre, Royal Observatory, Blackford Hill,
Edinburgh, EH9 3HJ, UK
7 Space Telescope Science Institute, Baltimore, MD 21218, USA
8Department of Physics and Astronomy, The Johns Hopkins University, Baltimore,
MD 21218, USA
9 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA
91109, USA
10Université Paris-Saclay, Université Paris Cité, Commissariat à l’énergie
atomique et aux énergies alternatives (CEA), Centre National de la Recherche
Scientifique (CNRS), Astrophysique, Instrumentation, Modélisation (AIM),
L’Orme des Merisiers 91191 St Aubin, France
11Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA
12Astronomy and Astrophysics Research Group, Department of Physics and
Astrophysics,
Vrije Universiteit Brussel, B-1050 Brussels, Belgium
13Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, 2200
Copenhagen, Denmark
14Stratospheric Observatory for Infrared Astronomy Science Center,
Universities Space Research Association,
NASA Ames Research Center, Moffett Field, CA 94035, USA
15School of Physics & Astronomy, Space Research Centre, Space Park Leicester,
University of Leicester, Leicester LE4 5SP, UK
16School of Physics and Astronomy, Queen’s Buildings, Cardiff University,
Cardiff, CF24 3AA, UK
17Laboratoire d’Etudes du Rayonnement et de la Matière en Astrophysique et
Atmosphères (LERMA), Observatoire de Paris, Paris Sciences & Lettres (PSL)
Research University, French National Centre for Scientific Research (CNRS),
Sorbonne Université, Paris, France
18 Université Paris-Saclay, Commissariat à l’énergie atomique et aux énergies
alternatives (CEA),
Detectors, Electronics and Computing for Physics (DEDIP), 91191, Gif-sur-
Yvette, France
19Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg,
Germany
20National Science Foundation’s National Optical-Infrared Astronomy Research
Laboratory, Tucson, AZ 85719, USA
21NASA Goddard Space Flight Center, Greenbelt, MD 20770, USA
22Centro de Astrobiología, Consejo Superior de Investigaciones Científicas -
Instituto Nacional de Técnica Aeroespacial (CSIC-INTA), Torrejón de Ardoz,
E-28850, Madrid, Spain
23Max-Planck Institut für Extraterrestrische Physik, D-85748, Garching,
Germany
24Leiden Observatory, 2300 RA Leiden, The Netherlands
25Dept. of Astrophysics, University of Vienna, A-1180 Vienna, Austria
26Institute for Particle Physics and Astrophysics, Eidgenössische Technische
Hochschule (ETH) Zürich, 8093 Zürich, Switzerland
27Institute of Astronomy, Katholieke Universiteit Leuven, 3001 Leuven, Belgium
∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
†These authors contributed equally to this work.

> The nearby Supernova 1987A was accompanied by a burst of neutrino emission,
> which indicates that a compact object (a neutron star or black hole) was
> formed in the explosion. There has been no direct observation of this
> compact object. In this work, we observe the supernova remnant with JWST
> spectroscopy finding narrow infrared emission lines of argon and sulphur.
> The line emission is spatially unresolved and blueshifted in velocity
> relative to the supernova rest frame. We interpret the lines as gas
> illuminated by a source of ionizing photons located close to the center of
> the expanding ejecta. Photoionization models show that the line ratios are
> consistent with ionization by a cooling neutron star or pulsar wind nebula.
> The velocity shift could be evidence for a neutron star natal kick.

Supernova 1987A (SN 1987A) in the Large Magellanic Cloud is the most recent
supernova (SN) in the Local Group of galaxies, so it has been extensively
studied. It originated from the collapse of a blue supergiant star of 15 to 20
solar masses ($\rm~{}M_{\odot}$). Observations taken 35 years after the
explosion show that the remnant of SN 1987A contains several distinct
components (Fig. 1) (?). The inner ejecta containing the synthesized heavy
elements produced in the explosion, are now interacting with an equatorial
ring (ER) and the surrounding circumstellar medium (?). The ER is thought to
have originated $\sim 20,000$ years before the explosion, probably from a
binary merger experienced by the progenitor star (?, ?). The inner ejecta
produce broad emission lines with velocities of 2000 to 4000
$\rm~{}km~{}s^{-1}$, powered by radioactive decay of 44Ti nuclei and x-rays
generated by the interaction with the ER (?).

A compact object, such as a neutron star (NS) or black hole (BH), is expected
to be located in the central region, but has not been observed. A neutrino
burst was detected coinciding with SN 1987A (?, ?, ?); its duration indicates
that the compact object is a NS (?). Ionizing photons emitted by such a NS
have been predicted to excite emission lines from heavy elements in the ejecta
with velocities 100 to 1000 $\rm~{}km~{}s^{-1}$ (?). In this study we searched
for those emission lines.

## Infrared spectroscopy of SN1987A

We observed the SN 1987A remnant on 16 July 2022 (day 12,927 after the SN
explosion) using the JWST with the Mid-Infrared Instrument/Medium Resolution
Spectrometer
(MIRI/MRS) and the Near-Infrared Spectrograph/Integral Field Unit
(NIRSpec/IFU). Both instruments produce data cubes, which contain a spectrum
for each spatial pixel (?).

The data cubes exhibited an excess of emission in the central pixels of the
ejecta (Fig. 2) at observed peak wavelengths of $4.5293\pm 0.0003$ µm in the
NIRSpec data and at $6.98606\pm 0.00003$ µm in the MRS data, with the latter
line being much brighter. We identify these lines as due to [Ar vi] and [Ar
ii], which have rest wavelengths 4.52922 and $6.985274$ µm, respectively.
Figure 2, A to P shows velocity slices of the MRS data cube centered on the
wavelength of the [Ar ii] line and Fig. 2, Q to R shows slices of the NIRSpec
cube, containing the [Ar vi] line. The [Ar ii] line in SN 1987A has been
observed previously (?), but the spatial and spectral resolution were
insufficient to separate the emission from the ejecta and the ER. The JWST
observations show it arises from a central source, separated from the ER.

Fig. 3A shows the [Ar ii] line profile from the central source, which has a
narrow strong peak and a weaker broad extension to the red. We fitted this
asymmetric line profile with two Gaussian functions. Fig. S2 compares the line
from the center of the ejecta with the corresponding [Ar ii] line from the ER.
The narrow emission line component in the center is displaced by
$-259.6\rm~{}km~{}s^{-1}$ with respect to the rest frame of SN 1987A, which
has a heliocentric velocity +286.5 km s-1 (?), whereas in the spectrum of the
entire ER the [Ar ii] line shows a much smaller blueshift of $-12.2\pm
1.6\rm~{}km~{}s^{-1}$ (?). The velocity difference and spatial distribution
indicate that the emission line from the central source is not scattered light
from the ER.

The spectral resolution at the [Ar ii] wavelength corresponds to a full-width
at half-maximum (FWHM) of $87.7\pm 1.1\rm~{}km~{}s^{-1}$ (?). We measure that
the narrower ejecta line component has a FWHM $=121.9\pm
1.3\rm~{}km~{}s^{-1}$, so it is marginally spectrally resolved. The emission
from the entire ER has a FWHM$=286\rm~{}km~{}s^{-1}$, and the broader
component has a value of $=362.4\rm~{}km~{}s^{-1}$. Table 1 lists the velocity
offsets, FWHMs and luminosities that we measured for the emission lines from
the central ejecta.

Table 1: Parameters for emission line components from the central ejecta. Luminosities calculations adopt a distance to SN 1987A of 49.6 kiloparsecs (?). No broad component was required for the $[$Ar iii$]$, $[$Sr iv$]$, and $[$Sr iii$]$ lines (fig. S3). The $[$Ar vi$]$ lines were measured in the NIRSpec data; all other lines are from the MIRI data. Uncertainties are $1\sigma$ errors on the fitting parameters; upper limits are 3$\sigma$. Dashes indicate not applicable.  |  | Narrrow component | Broad component
---|---|---|---
Line | Wavelength | Velocity offset | FWHM | Luminosity | Velocity offset | FWHM | Luminosity
| µm | $\rm~{}km~{}s^{-1}$ | $\rm~{}km~{}s^{-1}$ | $10^{30}\rm~{}erg~{}s^{-1}$ | $\rm~{}km~{}s^{-1}$ | $\rm~{}km~{}s^{-1}$ | $10^{30}\rm~{}erg~{}s^{-1}$
$[$Ar vi$]$ | 4.529 |  -269$\pm$24 |  $<$377.1 |  3.04$\pm$0.42 |  15$\pm$ 33 |  $<$461.1 |  2.88$\pm$0.40
$[$Ar ii$]$ | 6.985 |  -259.6$\pm$0.4 |  121.9$\pm$1.3 |  142.9$\pm$2.7 |  -200.6$\pm$3.7 |  362.4$\pm$8.0 |  104.1$\pm$3.0
$[$Ar iii$]$ | 8.991 |  |  | $<3.5$ | $-$ | $-$ | $-$
$[$Sr iv$]$ | 10.510 |  -234$\pm$12.9 |  145$\pm$35 |  2.13$\pm$0.45 | $-$ | $-$ | $-$
$[$Sr iii$]$.1 | 18.713 |  -288$\pm$13 |  158$\pm$39 |  4.50$\pm$0.82 | $-$ | $-$ | $-$

The NIRSpec spectral resolution at the [Ar vi] $4.529$-µm line is $\sim
261\rm~{}km~{}s^{-1}$. Fitting a similar two-component model as for the [Ar
ii] line, we found that the peak of the [Ar vi] line has consistent velocity
to that of [Ar ii], but the broad components differ from the narrow
blueshifted component (Table 1). Both the [Ar ii] and [Ar vi] peaks are
consistent with being spatially unresolved [the spatial resolution FWHM is
$0.^{\prime\prime}350$ and $0.^{\prime\prime}212$ for [Ar ii] and [Ar vi],
respectively (?)]. The position of maximum emission in the [Ar ii] and [Ar vi]
lines spatially coincide (within the uncertainties) with each other and with
the previously determined center of the ER (Fig. 1), which is assumed to be
the explosion site (?).

In the MIRI data we also detected lines of [S iv] $10.51$ µm and [S iii]
$18.71$ µm (fig. S3). Both sulphur lines show weak components from the central
ejecta, displaced from the SN 1987A velocity by a similar amount to the argon
lines. We can only place an upper limit on the luminosity of an [Ar iii] line
at $8.991$ µm line (Table 1). In the NIRSpec data, we detected the [Ca iv]
line at $3.207$ µm (fig. S5) ; however, its velocity is shifted to the red and
its maximum emission is located slightly north of the [Ar vi] peak, so it
might be unrelated to the argon lines (?).

## Interpretation of the central emission

The only narrow emission components from the ejecta that we identify are from
Ar and S, which are both produced by the nuclear burning (fusion) of O and Si.
This indicates that the emission arises from the inner core of the ejecta
(fig. S7 ). Explosion models of SN 1987A (?) show that the innermost $\sim
0.07\rm~{}M_{\odot}$ of ejecta is dominated by He and 56Ni (the latter decays
into 56Fe). The next zone outward, with a mass of $\sim 0.3\rm~{}M_{\odot}$,
contains material produced by oxygen burning – Si, S, Ar, Ca, and, in its
outer region, also O (?).

The same elements have been observed in young oxygen-rich supernova remnants
(SNRs) from massive star explosions, including Cas A, SNR G54.1+0.3 and SNR
0540-69.3 (?, ?, ?, ?) (supplemenatry text). Both SNR G54.1+0.3 and SNR
0540-69.3 contain pulsar wind nebulae (PWNe), powered by the spin-down power
from a rapidly rotating NS with a strong magnetic field. By contrast, Cas A
contains a central, cooling NS (CNS), with a surface temperature $\sim
1.8\times 10^{6}$ K (?) without a PWN. The emission lines in Cas A are mainly
excited by a shock in the outer parts of the SN ejecta.

The narrow [Ar ii], [Ar vi], [S iii] and [S iv] lines that we observe from the
center of SN 1987A must be excited by a source of ionizing photons or a shock
wave. Potential sources include photons from a PWN generated by a NS, photons
directly from a CNS, accretion onto a compact object, or shock excitation by a
PWN shock. We also considered five other possibilities (supplementary text):
ionization by radioactive 44Ti; excitation produced by x-rays from the ejecta
and circumstellar medium interaction; a surviving companion star; reflection
of the narrow line emission from the ER collision by dust; or emission from an
ingoing, reverse shock. All five possibilities were excluded.

The remaining potential explanations all include the presence of a young NS or
a BH in the center. We consider a BH unlikely because the progenitor star of
SN 1987A is thought to have had a total mass of $\lesssim 20\rm~{}M_{\odot}$
and an Fe core mass of $\lesssim 2\rm~{}M_{\odot}$ (?). The explosion ejected
at least $0.07\rm~{}M_{\odot}$ of 56Ni (?). That leaves less core mass
remaining than is required to form a BH, which is $\gtrsim 2.2\rm~{}M_{\odot}$
(?). Accretion of ejecta bound to the NS could be a possibility. The estimated
luminosity from this, $\sim 2\times 10^{32}\rm~{}erg~{}s^{-1}$ at 35 years
(?), is too low compared with the observed luminosity (Table 1), but the
relevant physical parameters are uncertain.

A hot CNS would have a temperature of $\gtrsim 10^{6}$ K and luminosity
$\gtrsim 10^{34}\rm~{}erg~{}s^{-1}$ (?), which is sufficient to ionize part of
the inner ejecta (?). An analogous case is the NS in Cas A, although it is
older and fainter than what is expected for SN 1987A.

The energy loss from a pulsar is in the form of electrons, positrons and
magnetic fields, which produce an expanding high-pressure PWN. High-frequency
synchrotron radiation from the relativistic particles can ionize the ejecta
(?). The luminosities and spectra of the PWNs discussed above are not
consistent with the bolometric and monochromatic luminosity in SN 1987A,
whereas a slowly rotating and/or weakly magnetized NS is compatible with the
observations (?). We expect an expanding PWN bubble to sweep up a dense shell
of ejecta, separated from the PWN by a shock, which could also excite the
lines we observed. The velocity of the dense shell would depend on the energy
input from the pulsar along with the density and expansion velocity of the
ejecta (?). These parameters are all uncertain, but estimates for SN 1987A
produce a shell velocity of 150 to 300 $\rm~{}km~{}s^{-1}$ and a PWN shock
velocity of 30 to 50 $\rm~{}km~{}s^{-1}$ (?).

## Photoionization models

To investigate these scenarios in more detail, we have updated a modelling
code (?) and applied it to model the photoionization from a PWN and from a CNS
(?). The resulting emission line spectrum is primarily determined by the
ionization parameter, $\xi=L_{\rm ion}/n_{\rm ion}r^{2}$, where $L_{\rm ion}$
is the ionizing luminosity, $n_{\rm ion}$ is the number density of ions and
$r$ is the distance to the ionizing source; it is affected to a lesser extent
by the spectral form, $L(\nu)$, where $\nu$ is the frequency (?). For the
ionizing spectrum of a PWN we assume a power law,
$L(\nu)\propto\nu^{-\alpha}$, with $\alpha=1.1$. This matches the
ultraviolet–to–x-ray spectrum of the Crab Nebula (?), and is within the range
of $0.5<\alpha<1.1$ previously found for x-rays in young PWNs (?, ?). We also
investigated using steeper power laws and found that it has little effect
(supplementary text). Previous work has constrained $L_{\rm ion}\lesssim
10^{36}\rm~{}erg~{}s^{-1}$ at energies above 13.6 eV (?).

For the young CNS scenario, previous models have found an upper limit to the
thermal luminosity of $\lesssim 4\times 10^{35}\rm~{}erg~{}s^{-1}$ and a
surface temperature of $\sim(1.5$ to $3)\times 10^{6}$ K at 35 years after the
explosion (?, ?, ?). Central compact objects, like the NS in Cas A, are
thought to have carbon-dominated atmospheres (?), although alternative models
producing a hard spectrum have been discussed (?). We adopted a carbon-
dominated CNS spectrum (?) with $L_{\rm ion}=3\times
10^{35}\rm~{}erg~{}s^{-1}$.

The models depend on the density in the ejecta core, which is uncertain. After
explosion, the mass inside the oxygen core of a 19 $\rm~{}M_{\odot}$
progenitor is $\sim 3.0\rm~{}M_{\odot}$, excluding the material that forms the
NS (?). The line profiles indicate a core velocity of $\sim
2200\rm~{}km~{}s^{-1}$ (?). For a uniform core density $n_{\rm ion}\approx
2.6\times 10^{3}(A/22)^{-1}\ \rm~{}cm^{-3}$ after 35 years , where $A$ is the
mean atomic weight in the O-Si-S-Ar core. However, both observations (?, ?, ?,
?) and simulations (?), indicate that clumping and filaments produce regions
with higher densities in this zone, whereas the Fe-rich regions can have lower
density as a result of early heating by decay of 56Ni (?). If a PWN is present
it produces a dense shell (as discussed above), which is likely to be unstable
and clumpy (?, ?). We therefore assume a density of $n_{\rm ion}=2.6\times
10^{4}$ cm-3, which is 10 times the average estimated above, in both models.

Figure 4 shows the lines predicted by the models of the PWN and CNS scenarios
for an ejecta composition (fig. S7) typical of explosive oxygen burning in a
SN with a $19-\rm~{}M_{\odot}$ progenitor star (?, ?). The blueshifted line
velocities and dust effects indicate that only a fraction of the total
ionizing luminosity is converted to the observed emission lines (supplementary
text). We therefore normalize the [Ar ii] luminosity by a covering factor,
chosen to match the observed luminosity, and scale the other lines to their
relative model luminosities (?). For consistency, we include only the
luminosities of the blueshifted narrow component.

The relative luminosities of the [Ar ii] $6.985$-µm and the [Ar vi] $4.529$-µm
lines are reproduced in both models (Fig. 4), with this assumed density and
$\xi\approx 0.3$ (?). Both models predict other high-ionization lines,
including [Ca iv] $3.207$-µm, [Ca v] $4.159$-µm, [Ar v] $7.902$-µm, [Ar iii]
$8.991$-µm, [S iv] $10.51$-µm, [Ar v] $13.10$-µm, [S iii] $18.71$-µm, [Ar iii]
$21.83$-µm, [O iv] $25.89$-µm and [Fe ii] $25.99$-µm (Fig. 4). Most of these
lines have previously been observed in young, oxygen-rich SNRs (supplementary
text). However, in SN 1987A the lines at wavelengths higher than $\sim 8$ µm
are much weaker than predicted (relative to the [Ar ii] line) or absent
entirely, although we do detect faint emission for [S iv] $10.51$-µm, [S iii]
$18.71$-µm (fig. S3).

## Dust effects and the observed emission

Infrared observations have shown that a large amount of dust was formed in the
ejecta of SN1987A (?, ?), consisting of a mixture of silicates and graphite
grains (?, ?). However, the composition of the dust close to the center is
likely to differ from that of the outer regions (?). Dust nucleation models
have predicted that silicates dominate over carbon-rich dust in the central
regions, (?, ?, ?, ?), of which Mg2SiO4 (forsterite) is expected to have the
highest abundance. Observations of dust in young SNRs have shown that
silicates are the dominant type of dust (?).

Silicate dust grains strongly absorb wavelengths longer than $\sim 8$ µm (fig.
S8)). This is a potential explanation for why we do not detect the long-
wavelength lines at the luminosities predicted by the photoionization models.

We are unable to estimate the optical depth $\tau_{\rm abs}$ at those
wavelengths, because of the unknown composition, grain size, and clumping.
However, $\tau_{\rm abs}$ is likely to be large (?). Figure 4 also shows the
dust absorption expected for an assumed optical depth $\tau_{\rm abs}=6.5$ at
$10$ µm and amorphous forsterite composition (?) and the line luminosities
after including this dust absorption. Althopugh the predicted luminosities of
lines shorter than $8$ µm are only slightly affected, those at longer
wavelengths are substantially suppressed. For [S iv] $10.51$-µm and [S iii]
$18.71$-µm lines, this correction brings the model predictions close to the
observed line luminosities. The observed upper limit on the [Ar iii]
$8.991$-µm line is consistent with the model.

The CNS model (Fig. 4A) predicts strong [Si i] lines at 1.607 and 1.645 µm.
These are not detected in our observations, perhaps because they are blended
with broad lines from the ejecta. Dust scattering, dust depletion and
ionization of silicon by the ultraviolet flux from the ejecta and ring could
all reduce the predicted luminosities of these lines (?).

We have also modeled the expected line luminosities from a PWN shock (?). Most
of the observed lines are well reproduced (supplementary text) (fig. S12),
without requiring dust effects. However, the [S iii] $18.71$-µm line is
underproduced by a factor $\sim 25$.

## The nature of the compact object

Given the uncertain model assumptions about the ionizing spectrum, density,
elemental abundances and dust properties, we conclude that the observed lines
can be explained by either the CNS or PWN photoionization models. Both require
the presence of a NS. Because thermal emission from the CNS should always be
present at some level, a combination of both models is also possible. For the
model parameters adopted above, the PWN shock model is less likely but not
excluded (?).

There have been previous indications of emission from a central object in SN
1987A. Sub-millimeter dust maps of the ejecta show a peak in the temperature
close to the central region, which might be due to extra energy input from a
central source (?). The peak temperature position does not coincide with the
location of maximum [Ar vi] emission in our observations (?); however, this
does not rule out additional heating from a compact object (?). X-ray
observations of SN 1987A have detected nonthermal emission at energies above
10 keV, which has been interpreted as being due to a central source (?, ?);
this would support the PWN scenario. However, the same observations have
alternatively been interpreted as thermal emission from the shocks in the ER
(?).

The position and velocity of the argon emission lines constrain any natal NS
kick velocity - a velocity imparted on the NS during its formation by
anisotropies in the explosion mechanism. The offset in position of the narrow
emission components corresponds to a velocity $416\pm 206\rm~{}km~{}s^{-1}$
given the time since the explosion (supplementary text). However, converting
this to a kick velocity depends on the separation between the NS and the line-
emitting region, which differs between the CNS and PWN models.

Our observations of high ionization emission lines in the centre of the SN
1987A remnant indicate the presence of a source of ionizing photons, probably
a NS. We cannot determine which scenario is more likely – young CNS or PWN –
but both require the presence of a NS.

Figure 1:  Optical image of SN 1987A taken in 2022, 35 years after explosion.
Hubble Space Telescope (HST) in the F625W filter (?), which is dominated by
H$\alpha$ emission. The freely expanding inner ejecta and ER are labelled.
White contours mark the [Ar vi] 4.529-µm emission observed with NIRSpec (40
and 70% of the maximal surface brightness in Fig. 2S). The white star denotes
the center of the ER (?).

Figure 2: JWST observations of the [Ar ii] and [Ar vi] lines in SN 1987A. (A
to O) Velocity slices of the MIRI/MRS data cube around the [Ar ii] line at
$6.985$ µm. Labels in each panel indicate the corresponding velocity range,
measured with respect to the systemic velocity of SN 1987A (?). Emission is
apparent from a central source, which has its maximum intensity offset in
velocity from the ER by $\sim-250\rm~{}km~{}s^{-1}$. The white circle in (F)
indicates the region used to extract the spectrum in Fig. 3A. (P) Residuals
between (F) and (A). (Q) NIRSpec observations of the [Ar vi] $4.529$-µm line,
combining all velocity slices between $-450$ and $+380\rm~{}km~{}s^{-1}$. (R)
The NIRSpec data stacked outside that velocity range. (S) Residuals between
panels (Q) and (R), separating the emission from the central source. The
bright point in the lower left areas of (Q) and (R) is an unrelated star.
Color bars are in units of megajansky (MJy) per steradian.  Figure 3: Velocity
profiles of the [Ar ii] 6.985-µm and [Ar vi] 4.529-µm lines from the central
source. (A) The [Ar ii] line in the MRS data (black histogram), extracted from
the region indicated in Fig. 2F. The orange curve is a model fitted to the
data, consisting of two Gaussian components (blue and grey dashed curves). (B)
The same [Ar ii] data (black histogram), compared to the [Ar vi] line from the
NIRSpec data (red histogram). The [Ar vi] data have a lower spectral
resolution and have been scaled by a factor of 76 for display.

Figure 4: Photoionization models compared to observed lines. (A) Model of the
oxygen-burning layers expected in a SN ionized by a young CNS. (B) Model of
the same material ionized by a PWN. The optical depth due to silicate dust
(magenta curve, on the right axis), assumes Mg2SiO4 composition and $\tau_{\rm
abs}=6.5$ at 10 µm in both cases. Circles show the model-predicted line
luminosities (left axis)before (red) and after (blue) including dust
absorption, which are connected by dashed vertical lines. Line identifications
are labelled near each red circle. Black stars are the observed luminosities
(Table 1) for the [S iii], [S iv], [Ar ii] and [Ar vi] lines, with error bars
showing 1$\sigma$ uncertainties, which are mostly smaller than the symbol
size. The upper limit for the [Ar iii] line (black triangle) is 3$\sigma$.

## References and Notes

  * 1. R. McCray, C. Fransson, Annu. Rev. Astron. Astrophys 54, 19 (2016).
  * 2. J. Larsson, et al., Astrophys. J. Lett. 949, L27 (2023).
  * 3. T. Morris, P. Podsiadlowski, Science 315, 1103 (2007).
  * 4. T. Morris, P. Podsiadlowski, Mon. Not. R. Astron. Soc. 399, 515 (2009).
  * 5. E. N. Alekseev, L. N. Alekseeva, V. I. Volchenko, I. V. Krivosheina, Soviet Journal of Experimental and Theoretical Physics Letters 45, 589 (1987).
  * 6. R. M. Bionta, et al., Phys. Rev. Lett.  58, 1494 (1987).
  * 7. K. Hirata, et al., Phys. Rev. Lett.  58, 1490 (1987).
  * 8. A. Burrows, Astrophys. J. 334, 891 (1988).
  * 9. R. A. Chevalier, C. Fransson, Astrophys. J. 395, 540 (1992).
  * 10. Materials and methods are available as supplementary materials .
  * 11. R. G. Arendt, et al., Astron. J. 151, 62 (2016).
  * 12. P. Gröningsson, et al., Astron. Astrophys. 492, 481 (2008).
  * 13. O. C. Jones, et al., Mon. Not. R. Astron. Soc. 523, 2519 (2023).
  * 14. G. Pietrzyński, et al., Nature 567, 200 (2019).
  * 15. D. Alp, et al., Astrophys. J. 864, 174 (2018).
  * 16. S. E. Woosley, A. Heger, Phys. Rep. 442, 269 (2007).
  * 17. R. A. Chevalier, R. P. Kirshner, Astrophys. J. 233, 154 (1979).
  * 18. J. D. T. Smith, et al., Astrophys. J. 693, 713 (2009).
  * 19. T. Temim, P. Slane, S. P. Reynolds, J. C. Raymond, K. J. Borkowski, Astrophys. J. 710, 309 (2010).
  * 20. B. J. Williams, et al., Astrophys. J. 687, 1054 (2008).
  * 21. W. C. G. Ho, C. O. Heinke, Nature 462, 71 (2009).
  * 22. T. Sukhbold, S. E. Woosley, A. Heger, Astrophys. J. 860, 93 (2018).
  * 23. P. Bouchet, et al., Astron. Astrophys. 245, 490 (1991).
  * 24. C. Musolino, C. Ecker, L. Rezzolla, arXiv e-prints p. arXiv:2307.03225 (2023).
  * 25. R. A. Chevalier, Phys. Rep. 256, 95 (1995).
  * 26. M. V. Beznogov, A. Y. Potekhin, D. G. Yakovlev, Phys. Rep. 919, 1 (2021).
  * 27. M. Lyutikov, et al., Mon. Not. R. Astron. Soc. 489, 2403 (2019).
  * 28. X.-H. Li, F.-J. Lu, Z. Li, Astrophys. J. 682, 1166 (2008).
  * 29. E. V. Gotthelf, Astrophys. J. 591, 361 (2003).
  * 30. M. V. Beznogov, J. Novak, D. Page, A. R. Raduta, Astrophys. J. 942, 72 (2023).
  * 31. D. Page, et al., Astrophys. J. 898, 125 (2020).
  * 32. J. A. J. Alford, J. P. Halpern, Astrophys. J. 944, 36 (2023).
  * 33. J. Spyromilio, P. A. Pinto, European Southern Observatory Conference and Workshop Proceedings, Volume 37, I. J. Danziger, K. Kjär, eds. (European Southern Observatory, Garching, 1991), p. 423.
  * 34. H. Li, R. McCray, Astrophys. J. 387, 309 (1992).
  * 35. J. Larsson, et al., Astrophys. J. 768, 89 (2013).
  * 36. M. Matsuura, et al., Mon. Not. R. Astron. Soc. 469, 3347 (2017).
  * 37. M. Gabler, A. Wongwathanarat, H.-T. Janka, Mon. Not. R. Astron. Soc. 502, 3264 (2021).
  * 38. J. M. Blondin, R. A. Chevalier, Astrophys. J. 845, 139 (2017).
  * 39. D. H. Wooden, et al., Astrophys. J. Suppl. Ser. 88, 477 (1993).
  * 40. M. Matsuura, et al., Science 333, 1258 (2011).
  * 41. M. Matsuura, et al., Astrophys. J. 800, 50 (2015).
  * 42. R. Wesson, M. J. Barlow, M. Matsuura, B. Ercolano, Mon. Not. R. Astron. Soc. 446, 2089 (2015).
  * 43. A. Sarangi, I. Cherchneff, Astron. Astrophys. 575, A95 (2015).
  * 44. A. Sarangi, M. Matsuura, E. R. Micelotta, Space Sci. Rev. 214, 63 (2018).
  * 45. A. Sluder, M. Milosavljević, M. H. Montgomery, Mon. Not. R. Astron. Soc. 480, 5580 (2018).
  * 46. A. Sarangi, Astron. Astrophys. 668, A57 (2022).
  * 47. P. Cigan, et al., Astrophys. J. 886, 51 (2019).
  * 48. H.-P. Gail, A. Tamanai, A. Pucci, R. Dohmen, Astron. Astrophys. 644, A139 (2020).
  * 49. E. Greco, et al., Astrophys. J. Lett. 908, L45 (2021).
  * 50. E. Greco, et al., Astrophys. J. 931, 132 (2022).
  * 51. D. Alp, J. Larsson, C. Fransson, Astrophys. J. 916, 76 (2021).
  * 52. O. Jones, JWST MRS SN 1987A, Mikulski Archive for Space Telescopes https://doi.org/10.17909/pz4v-ej33 (2023) .
  * 53. J. Larsson, JWST NIRSpec SN 1987A, Mikulski Archive for Space Telescopes https://doi.org/10.17909/175h-7x33 (2023) .
  * 54. P. Kavangh, JWST MRS 10 Lac, Mikulski Archive for Space Telescopes https://doi.org/10.17909/7an5-yp54 (2023) .
  * 55. J. Larsson, HST/WFC3 F625W & F657N, 2022, SN 1987A, Mikulski Archive for Space Telescopes https://doi.org/10.17909/53qa-pd67 (2023) .
  * 56. C. Fransson, Photoionization code: Zenodo, claesob/SNcodes
https://doi.org/10.5281/zenodo.10476657 (2024) .

  * 57. C. Fransson, Shock code: Zenodo, claesob/SN87A_shock
https://doi.org/10.5281/zenodo.10476710 (2024) .

  * 58. O. C. Jones, et al., Astrophys. J. 958, 95 (2023).
  * 59. M. Wells, et al., Publ. Astron. Soc. Pac. 127, 646 (2015).
  * 60. P. Jakobsen, et al., Astron. Astrophys. 661, A80 (2022).
  * 61. S. Gonzaga, W. Hack, A. Fruchter, J. Mack, The DrizzlePac Handbook (2012).
  * 62. Gaia Collaboration, et al., Astron. Astrophys. 674, A1 (2023).
  * 63. K. M. Pontoppidan, et al., Astrophys. J. Lett. 936, L14 (2022).
  * 64. J. Larsson, et al., Astrophys. J. 873, 15 (2019).
  * 65. I. Argyriou, et al., Astron. Astrophys. 675, A111 (2023).
  * 66. I. D. Howarth, J. Murray, D. Mills, D. S. Berry, DIPSO: Spectrum analysis code, Astrophysics Source Code Library, record ascl:1405.016 (2014).
  * 67. P. Virtanen, et al., Nature Methods 17, 261 (2020).
  * 68. C. B. Tarter, W. H. Tucker, E. E. Salpeter, Astrophys. J. 156, 943 (1969).
  * 69. N. R. Badnell, Astrophys. J. Suppl. Ser. 167, 334 (2006).
  * 70. Z. Altun, A. Yumak, N. R. Badnell, S. D. Loch, M. S. Pindzola, Astron. Astrophys. 447, 1165 (2006).
  * 71. Z. Altun, et al., Astron. Astrophys. 474, 1051 (2007).
  * 72. S. A. Abdel-Naby, D. Nikolić, T. W. Gorczyca, K. T. Korista, N. R. Badnell, Astron. Astrophys. 537, A40 (2012).
  * 73. E. A. Bleda, Z. Altun, N. R. Badnell, Astron. Astrophys. 668, A72 (2022).
  * 74. D. A. Verner, G. J. Ferland, K. T. Korista, D. G. Yakovlev, Astrophys. J. 465, 487 (1996).
  * 75. W. Cunto, C. Mendoza, F. Ochsenbein, C. J. Zeippen, Astron. Astrophys. 275, L5 (1993).
  * 76. J. S. Kaastra, R. Mewe, Astron. Astrophys. Suppl. Ser. 97, 443 (1993).
  * 77. T. R. Kallman, R. McCray, Astrophys. J. Suppl. Ser. 50, 263 (1982).
  * 78. G. Del Zanna, K. P. Dere, P. R. Young, E. Landi, Astrophys. J. 909, 38 (2021).
  * 79. S. Nahar, Atoms 8, 68 (2020).
  * 80. NIST Atomic Spectra Database https://physics.nist.gov/asd .
  * 81. P. A. M. van Hoof, Galaxies 6, 63 (2018).
  * 82. J. Pelan, K. A. Berrington, Astron. Astrophys. Suppl. Ser. 110, 209 (1995).
  * 83. H. E. Saraph, P. J. Storey, Astron. Astrophys. Suppl. Ser. 115, 151 (1996).
  * 84. S. N. Nahar, B. Shafique, European Physical Journal D 77, 45 (2023).
  * 85. K. Kjær, B. Leibundgut, C. Fransson, A. Jerkstrand, J. Spyromilio, Astron. Astrophys. 517, A51 (2010).
  * 86. C. Fransson, et al., Astrophys. J. 768, 88 (2013).
  * 87. H. Li, R. McCray, R. A. Sunyaev, Astrophys. J. 419, 824 (1993).
  * 88. C. Kozma, C. Fransson, Astrophys. J. 497, 431 (1998).
  * 89. J. Spyromilio, W. P. S. Meikle, D. A. Allen, Mon. Not. R. Astron. Soc. 242, 669 (1990).
  * 90. M. R. Haas, et al., Astrophys. J. 360, 257 (1990).
  * 91. J. P. Hughes, C. E. Rakowski, D. N. Burrows, P. O. Slane, Astrophys. J. Lett. 528, L109 (2000).
  * 92. D. Milisavljevic, R. A. Fesen, Astrophys. J. 772, 134 (2013).
  * 93. R. A. Chevalier, C. Fransson, Astrophys. J. 420, 268 (1994).
  * 94. T. K. Nymark, C. Fransson, C. Kozma, Astron. Astrophys. 449, 171 (2006).
  * 95. R. Strickland, J. M. Blondin, Astrophys. J. 449, 727 (1995).
  * 96. R. S. Sutherland, G. V. Bicknell, M. A. Dopita, Astrophys. J. 591, 238 (2003).
  * 97. C. Fransson, C. Kozma, New Astron. Rev. 46, 487 (2002).
  * 98. A. Jerkstrand, C. Fransson, C. Kozma, Astron. Astrophys. 530, A45 (2011).
  * 99. S. E. Boggs, et al., Science 348, 670 (2015).
  * 100. C. Kozma, C. Fransson, Astrophys. J. 390, 602 (1992).
  * 101. N. N. Chugai, R. A. Chevalier, R. P. Kirshner, P. M. Challis, Astrophys. J. 483, 925 (1997).
  * 102. K. K. Madsen, et al., Astrophys. J. 801, 66 (2015).
  * 103. A. Bamba, et al., Publ. Astron. Soc. Jpn. 74, 1186 (2022).
  * 104. R. A. Chevalier, Astrophys. J. Lett. 539, L45 (2000).
  * 105. J. J. Hester, Annu. Rev. Astron. Astrophys 46, 127 (2008).
  * 106. T. Temim, et al., Astrophys. J. 836, 129 (2017).
  * 107. R. G. Arendt, E. Dwek, G. Kober, J. Rho, U. Hwang, Astrophys. J. 786, 55 (2014).
  * 108. J. Rho, et al., Mon. Not. R. Astron. Soc. 479, 5101 (2018).
  * 109. F. Molster, C. Kemper, Space Sci. Rev. 119, 3 (2005).
  * 110. T. Henning, Annu. Rev. Astron. Astrophys 48, 21 (2010).
  * 111. M. Min, J. W. Hovenier, A. de Koter, Astron. Astrophys. 404, 35 (2003).
  * 112. B. T. Draine, B. S. Hensley, Astrophys. J. 910, 47 (2021).
  * 113. A. Tamanai, A. Pucci, R. Dohmen, H.-P. Gail, Astrophys. J. 845, 6 (2017).
  * 114. C. Jäger, J. Dorschner, H. Mutschke, T. Posch, T. Henning, Astron. Astrophys. 408, 193 (2003).
  * 115. J. Dorschner, B. Begemann, T. Henning, C. Jaeger, H. Mutschke, Astron. Astrophys. 300, 503 (1995).
  * 116. C. Jaeger, H. Mutschke, B. Begemann, J. Dorschner, T. Henning, Astron. Astrophys. 292, 641 (1994).
  * 117. L. B. Lucy, I. J. Danziger, C. Gouiffes, P. Bouchet, IAU Colloq. 120: Structure and Dynamics of the Interstellar Medium, G. Tenorio-Tagle, M. Moles, J. Melnick, eds. (Springer-Verlag, Berlin, Lecture Notes in Physics, 1989), vol. 350, p. 164.
  * 118. T. Kangas, et al., Mon. Not. R. Astron. Soc. 511, 2977 (2022).
  * 119. B. Müller, Living Reviews in Computational Astrophysics 6, 3 (2020).
  * 120. F. J. M. Rietmeijer, A. Pun, Y. Kimura, J. A. Nuth, Icarus 195, 493 (2008).
  * 121. H.-T. Janka, Astrophys. J. 837, 84 (2017).
  * 122. M. Ono, et al., Astrophys. J. 888, 111 (2020).
  * 123. S. Orlando, et al., Astron. Astrophys. 636, A22 (2020).
  * 124. A. De Luca, Journal of Physics Conference Series (2017), vol. 932 of Journal of Physics Conference Series, p. 012006.

ACKNOWLEDGEMENTS: We thank R. Chevalier, T. Janka, S. Woosley and the
anonymous referees for constructive comments on the paper and useful
discussions. MIRI development was supported by NASA, the European Space Agency
(ESA), the Belgian Science Policy Office (BELSPO), the Centre Nationale
d’Etudes Spatiales (CNES), the Danish National Space Centre, Deutsches Zentrum
fur Luft-und Raumfahrt (DLR), Enterprise Ireland, Ministerio De Economiá y
Competividad, Netherlands Research School for Astronomy (NOVA), Netherlands
Organisation for Scientific Research (NWO), Science and Technology Facilities
Council, Swiss Space Office, Swedish National Space Agency, and UK Space
Agency This study is based on observations made with the NASA/ESA/CSA JWST.
The data were obtained from the Mikulski Archive for Space Telescopes at the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127
for JWST.

Funding: C.F., J.L. and G.Ö. acknowledge support from the Swedish National
Space Agency. J.L. acknowledges support from the Knut & Alice Wallenberg
Foundation. M.J.B. acknowledges support from European Research Council
Advanced Grant ERC-2015-AdG-694520-SNDUST. P.J.K. and J.J. acknowledge support
from Science Foundation Ireland/Irish Research Council Pathway programme under
Grant Number 21/PATH-S/9360. O.C.J. acknowledges support from an STFC Webb
fellowship. M.M. and N.H. acknowledge support through NASA/JWST grant
80NSSC22K0025. M.M. and L.L. acknowledge support from the NSF through grant
2054178. J.H. was supported by a VILLUM FONDEN Investigator grant (project
number 16599). O.N. acknowledges support from STScI Director’s Discretionary
Fund. T.Tikkanen acknowledges financial support from the UK Science and
Technology Facilities Council, and the UK Space Agency. Part of the research
of M.M. and N.H. was carried out at the Jet Propulsion Laboratory, California
Institute of Technology, under a contract with the National Aeronautics and
Space Administration (80NM0018D0004). L.C. acknowledges support by grant
PIB2021-127718NB-100 from the Spanish Ministry of Science and Innovation/State
Agency of Research MCIN/AEI/10.13039/501100011033. J.A.D.L.B. and B.V. thank
the Belgian Federal Science Policy Office (BELSPO) for the provision of
financial support in the framework of the PRODEX Programme of the European
Space Agency (ESA).

Authors contributions: C.F. led the project, performed the modeling, wrote
most of the text and contributed to the data analysis. M.J.B. and P.B. wrote
the proposal science case, designed the observations and coordinated the SN
1987A Guaranteed Time Observations (GTO) Program #1232. M.J.B. also
contributed to the data analysis and writing. P.J.K. performed the MIRI MRS
data reduction and contributed to the analysis and writing. J.L. contributed
to the data analysis and writing. O.C.J. designed the observations and was the
science lead for the SN 1987A observations program. B.S. performed the NIRSpec
data reduction. M.M. is the co-PI of the SN 1987A observations program,
provided three of the nine hours of GTO time used and contributed to the data
analysis. T.Te. contributed to the interpretation of the data. A.S.H.
optimized the observational program. N.H. and A.R. assisted with data
processing. A.C., A.G., J.A.D.L.B., J.H., J.J., L.L., O.K., O.N., O.D.F.,
R.G., R.M.L., R.W. and T.Ti. reviewed the manuscript and provided comments for
incorporation. G.S.W., L.C., E.F.v.D., M.G., Th.H., P.-OL., G.Ö., T.P.R., B.V.
are MIRI PI and co-PIs, who reviewed the manuscript and provided comments for
incorporation.

Competing interests: The authors declare that they have no competing
interests.

Data and materials availability: The observations are available at the
Mikulski Archive for Space Telescopes https://mast.stsci.edu/ under proposal
IDs 1232 and 1524 for JWST, and 16789 for HST. The specific MRS and NIRSpec
SN1987A observations we used are archived (?, ?), as are the 10 Lac
observation used for calibration (?) and the HST observations (?).

The source code used for the PWN and CNS models is available at
https://github.com/claesob/SN-codes and archived at Zenodo (?), while the
source code for the shock models is available at
https://github.com/claesob/SN87A_shock and archived at Zenodo (?).

Supplementary Materials

  * Materials and Methods

  * Supplementary Text

  * Figs. S1 to S12

  * Tables S1 to S3

  * References (58-125)

Supplementary Materials for

Emission lines due to ionizing radiation from a compact object in the remnant
of Supernova 1987A

C. Fransson∗, M. J. Barlow, P. J. Kavanagh, J. Larsson, O. C. Jones, B.
Sargent, M. Meixner, P. Bouchet, T. Temim, G. S. Wright, J. A. D. L.
Blommaert, N. Habel , A. S. Hirschauer, J. Hjorth, L. Lenkić, T. Tikkanen, R.
Wesson, A. Coulais, O. D. Fox, R. Gastaud, A. Glasse, J. Jaspers, O. Krause,
R. M. Lau, O. Nayak, A. Rest, L. Colina, E. F. van Dishoeck, M. Güdel Th.
Henning, P.-O. Lagage, G. Östlin, T. P. Ray, B. Vandenbussche

∗ Corresponding author<EMAIL_ADDRESS>
This PDF file includes:

Materials and Methods
Supplementary Text
Figs. S1 – S12
Tables S1 – S3
References 58 – 125

##  1 Materials and methods

###  1.1 Analysis of the narrow emission lines from the ejecta

Detailed descriptions of the MIRI/MRS and NIRSpec/IFU observations (Proposal
ID 1232, Principal Investigator: G. Wright) and data reduction can be found in
sect. 2 of (?) and sect. 2 of (?), respectively. The MRS observations were
obtained on 2022 July 16 and the IFUs provided complete spectral coverage from
4.9 to 27.9 µm at medium spectral resolving powers
$R=\lambda/\Delta\lambda\sim$ 4000–1500 (?). Here $\Delta\lambda$ is the FWHM
of of a spectral line at wavelength $\lambda$. The MRS fields of view range
from $3.^{\prime\prime}2\times 3.^{\prime\prime}7$ in channel 1, the shortest
wavelength channel, to $6.^{\prime\prime}6\times 7.^{\prime\prime}7$ in the
long wavelength channel 4 (?). The total MRS integration time was 9396 s,
corresponding to approximately 3.75 hr when including all overheads. The
NIRSpec/IFU observations were obtained on the same date, using the
G140M/F100LP, G235M/F170LP, and G395M/F290LP grating/filter combinations. The
gratings covered the wavelength ranges 0.97–1.88 µm (G140M), 1.66–3.15 µm
(G235M), and 2.87–5.20 µm (G235M), with a spectral resolving power that
increased from $R\sim 700$ at the shortest wavelengths to $R\sim 1300$ at the
longest wavelengths in each grating (?). The total exposure times were 1751 s
for each of G140M and G235M, and 1225 s for G395M. The NIRSpec/IFU field of
view was $3.^{\prime\prime}1\times 3.^{\prime\prime}2$ at all wavelengths. The
fields of view of both instruments covered both the ER and the ejecta.

We used a recent HST Wide Field Camera 3 (WFC3) observation in the F657N
filter (obtained on 2022 September 6, program ID 16789) to perform the
alignment between the MRS and NIRSpec data. Another HST WFC3 image from this
program taken in the F625W filter (obtained on 2022 September 5) was used for
illustrating the main emission components of the system in Fig. 1. Both
observations consisted of four dithered exposures, with total exposure times
of 2600 s (F657N) and 1080 s (F625W). The exposures were combined using
DrizzlePac (?), which also applies distortion corrections and removes cosmic
rays. The resulting images were subsequently aligned with Gaia data release 3
(?) to improve the astrometry.

####  1.1.1 Spatial distribution of the Ar emission

Astrometric uncertainties can be introduced into JWST data by guide star
catalogue errors and roll uncertainty (?). Therefore, the MRS and NIRSpec data
had to be spatially aligned to extract spectra from identical regions and
compare the spatial distribution of the emission. For this we used the HST
WFC3 obseervations above.

For NIRSpec, two stars within the field of view (FOV), as well as several
hotspots in the ER, were used to perform the alignment with HST (?). However,
no suitable stars were available within the FOV of MRS. The ER of SN 1987A was
therefore used to perform the alignment for MRS. We produced a continuum-
subtracted image of the H i 7.460 µm line, which we expect to originate from
the same hotspots in the ER as the H$\alpha$ emission in the HST F657N image
(?). The limited spatial resolution of the MRS means that most of the hotspots
are blended into a diffuse structure. We fitted three bright hotspots with
Gaussians, from which we determined the offset of their centroids from the
corresponding hotspots in the HST image. In addition, we fitted the whole ER
with an elliptical model with diffuse emission (described in (?)) and compared
the offsets of the centers of the ellipses, which gave consistent results for
the offset. The uncertainty in the alignment between MRS and HST estimated
from the hotspot fits was $\sim 0.^{\prime\prime}03$, which dominates the
uncertainty in the comparison between MRS and NIRSpec.

To assess the spatial distribution of the [Ar ii] narrow component and whether
it is spatially resolved, we compared the radial profile to an empirical MRS
PSF. While Star 3 (RA=05:35:28.26, Dec=-69:16:11.89) is seen in the MRS
spectral band 1C datacube, its location of at the very edge makes it
unsuitable for such a comparison. Instead we made use of a JWST calibration
observation of a spatially unresolved source, the O-star 10 Lac
(RA=22:39:15.67, Dec=+39:03:00.97). We reprocessed these data (PID 1524, PI:
D. Law) following the same method described in (?) to produce an MRS band 1C
datacube. To enable a direct comparison between the two, we constructed
datacubes for both 10 Lac and SN1987A in the MRS IFU coordinate system so the
broadening in the along-slice direction is in the same dimension (?).

From the 10 Lac data, we isolated and summed the velocity slices around the
[Ar ii] line and determined the source centroid by fitting a two-dimensional
Gaussian model to the PSF peak. We defined annular bins of 1 pixel width
centered on the centroid and determined the surface brightness in each bin
before normalising by the innermost bin. The resulting 10 Lac radial profile
is shown in Fig. S1B.

To extract the radial profile of the [Ar ii] narrow component we isolated the
narrow component from the [Ar ii] broad component, the [Ar ii] ER emission,
and the dust continuum in the ring. This was performed by first producing a
datacube in a narrow velocity range around the [Ar ii] line. Datacubes were
also produced from narrow velocity ranges to the low and high end of the [Ar
ii] line and used to measure and subtract the continuum in the [Ar ii]
datacube, removing the dust continuum. We fitted the spectra extracted from
each pixel of the continuum-subtracted datacube with a multi-component
Gaussian model. The central wavelengths and widths of the two components
representing the broad and narrow [Ar ii] line were fixed to the best fitting
values in Table. 1, which were determined from the spectrum extracted from the
entire ejecta region. Additional Gaussian components were included to account
for the [Ar ii] emission from the ER. The resulting models were subtracted to
leave only the [Ar ii] narrow component in the datacube. As with the 10 Lac
cubes, the frames were then summed, the centroid determined, and the radial
profile extracted using annular bins of 1 pixel width. The resulting radial
profile and comparison with 10 Lac is shown in Fig. S1B. We fitted the
profiles with Gaussian models to determine the FWHMs of the [Ar ii] narrow
component and 10 Lac to be $0.^{\prime\prime}353~{}(\pm 0.^{\prime\prime}006)$
and $0.^{\prime\prime}350~{}(\pm 0.^{\prime\prime}005)$, respectively. The
best fits are shown in Fig. S1-right. We therefore conclude that the [Ar ii]
narrow component is spatially unresolved.

Figure S1:  Radial profiles of the [Ar vi] and [Ar ii] lines. (A) Profile of
the [Ar vi] line (blue circles) and ‘Star 3’ (orange stars). The enhanced
emission in the [Ar vi] profile compared to Star 3 outside
$0.^{\prime\prime}2$ may be due to residual continuum emission from the
ejecta. (B) Profile of the narrow component of the [Ar ii] line (blue
circles). The flux profile of 10 Lac (orange stars) is included as a
comparison of an unresolved points source in MRS band 1C. The best-fit
Gaussian models are shown by the overlapping dashed blue and solid orange
lines for the [Ar ii] source and 10 Lac, respectively.

We performed an equivalent analysis of the radial profile of the [Ar vi]
emission compared with the radial profile of Star 3 (Fig. S1A), which
represents an unresolved point source. The radial profile of the [Ar vi]
source was obtained from the continuum subtracted image (Fig. 2S), while that
of the star was obtained from the original image at the same wavelength (Fig.
2Q). The lower signal and lower spectral resolution in the [Ar vi] line
compared to [Ar II] makes the separation of the line emission into two
components uncertain, so we show the profile of the full line in Fig. S1A. We
find that the peak of the emission is unresolved, while there is some evidence
of more extended emission outside 0.”2. However this could be due to imperfect
continuum subtraction.

The centroid of the [Ar vi] source in the background-subtracted image is
located $38\pm 22$ milliarcseconds (mas) east and $31\pm 22$ mas south of the
geometric center of the ER determined from HST optical images (?). For
comparison, the centroid of the blue component of the [Ar ii] emission is
located $65\pm 26$ mas east and $51\pm 30$ mas south of the geometric center.
The [Ar ii] and [Ar vi] positions are consistent with each other and within 3
sigma of the geometric center. The uncertainty in the [Ar vi] position is
dominated by statistical uncertainties, while that of [Ar ii] is dominated by
the systematic astrometric uncertainties (see above).

####  1.1.2 Spectral line fitting and search for additional lines

To model the emission lines in the MRS spectra we used the Gaussian emission
line fitting (elf) routines within the dipso spectral analysis package (?).
For each line component, elf returns an estimate for the best-fitting central
wavelength or velocity, simultaneously with 1$\sigma$ uncertainties on each
line’s FWHM and total flux The resulting models are shown in Figs. 3, S2 and
S3.

Figure S2: Velocity profiles and extraction regions for the [Ar ii] 6.985 µm line. (A) Velocity profile of the central source (black histogram) together with velocity profiles extracted from different parts of the ER (blue, orange, brown and grey histograms were extracted from the north, east, south and west ER, respectively). (B) Extraction regions for the velocity profiles in A, superposed on an integrated image of the [Ar ii] emission. There is a velocity shift between the central component and the different ER extractions. The offset of the stronger ejecta component of the [Ar ii] line is $-259.6\pm 0.4\rm~{}km~{}s^{-1}$ and the FWHM is $121.9\pm 1.3\rm~{}km~{}s^{-1}$, narrower than the FWHM of the ring extractions (see Table S1). There is a red extension of the line from the central source.  Figure S3: Velocity profiles of the [Ar vi] 4.529 µm, [S iii] 18.71 µm and [S iv] 10.51 µm lines from the central region. (A) The [Ar vi] line (black histogram). The orange curve is a model fitted to the data, consisting of two Gaussians (blue and red dashed lines). (B) The [S iii] line (black histogram). The orange curve is a model fitted to the data, consisting of two Gaussians (blue and red dashed lines). (C) The [S iv] line (black histogram) fitted by a Gaussian model (orange curve). In the case of [Ar vi], the two Gaussian components are associated with the central source, with a minimal contribution from the ER (cf. Fig. S4). For [S iii], the red component is dominated by scattered light from the ER. For all three lines, the blue component has a similar velocity.  Table S1: Parameters for the [Ar ii] line in different parts of the ER. The four extraction regions are shown in Fig. S2B. The lines were fitted with a Gaussian model.  Region | Offset | FWHM
---|---|---
| $\rm~{}km~{}s^{-1}$ | $\rm~{}km~{}s^{-1}$
North | $-153.3\pm 2.0$ | $253.0\pm 4.7$
West | $-26.6\pm 0.8$ | $225.8\pm 1.8$
South | $127.7\pm 3.2$ | $263.1\pm 8.1$
East | $22.2\pm 3.4$ | $231.3\pm 9.6$

The NIRSpec spectra were fitted using the curve_fit function of scipy.optimize
(?), which gave consistent results with elf. The models were constrained by
requiring the FWHM of the Gaussians to be no smaller than the spectral
resolution at the relevant wavelengths. For the [Ar vi] line, which is well
modelled with two Gaussians (Fig. S3), the uncertainty ranges on the FWHM of
both components reached the resolution limit, so only upper limits on the FWHM
are provided in Table 1.

While the [Ar ii] $\lambda~{}6.985$ µm, [Ar vi] $\lambda~{}4.529$ µm lines are
in regions of low dust continuum emission, the comparison with the mid-
infrared (MIR) spectra of young supernova remnants (see below) led us to
search for additional lines at longer wavelengths. The strong background from
the dust continuum, as well as the decreasing spatial resolution with
wavelength makes it increasingly more difficult to separate lines emitted by
the central ejecta from the strong ER and interstellar medium (ISM) emission.
However, the velocity offset seen in especially the [Ar ii] $\lambda~{}6.985$
µm helps in separating narrow ejecta components from the ER and ISM
contributions. In this way, we detected two more weak lines with similar
velocity offsets as the [Ar ii] line. Figure S3 shows the [S iii] $18.71$ µm
and [S iv] $10.51$ µm line profiles from the central ejecta. The velocity
offsets and luminosities are summarized in Table 1. These velocities are
consistent with the offset in the [Ar ii] line, and are therefore likely to
originate from the same source. We are only able to place a 3$\sigma$ upper
limit on the flux of [Ar iii] $8.991$ µm, corresponding to a line luminosity
upper limit of $<3.5\times 10^{30}\rm~{}erg~{}s^{-1}$.

Figure S4: Spectra from the central source and the ER at 4.30-4.60 µm. The
black and red histograms show spectra extracted from the central region and
ER, respectively. A H i line is present at 4.377 µm, which is stronger than
the [Ar vi] line in the ER, but which does not appear in the central
extraction region. This indicates that scattered light from the ER makes a
minimal contribution to the [Ar vi] profile extracted from the central region.

Figure S5: Line profile and spectral map of the [Ca iv] 3.207 µm line. (A)
Velocity profile extracted from the central region (black histogram) fitted
with a Gaussian model (orange curve). (B) Velocity slice integrated over the
[Ca iv] line in panel A. The white contours show the [Ar vi] 4.529 µm emission
(Fig. 2S), while the cross shows the geometric center of the ER (?).  Figure
S6: Continuum emission in NIRSpec and the central [Ar vi] source. The image of
the continuum (color scale) was produced by integrating over the 3.42–3.58 µm
wavelength interval. The black contours show the central [Ar vi] 4.529 µm
emission (Fig. 2S). There is a weak enhancement in the continuum just north of
the [Ar vi] source.

In the NIRSpec data we identify two additional lines which could be related to
the central source, [Ca iv] 3.207 µm and [Ca v] 4.159 µm. The [Ca v] line,
however, is likely blended with an H2 line at 4.181 µm (?), so its properties
are uncertain. We model the [Ca iv] line with a single Gaussian with FWHM
$854\pm 40\rm~{}km~{}s^{-1}$ and a peak velocity shifted to the red by $676\pm
17\rm~{}km~{}s^{-1}$ (Fig. S5A), substantially different from the Ar and S
lines. The [Ca iv] emission is located $3\pm 10$ mas east and $114\pm 15$ mas
north of the center (S5B), which implies a significant offset to the north
compared to the [Ar vi] emission. Expressed in terms of an equivalent velocity
in the plane of the sky for the freely expanding ejecta, (with velocity
proportional to the distance from the center), the [Ca iv] emission is offset
by $757\pm 96\rm~{}km~{}s^{-1}$ to the north of the center. We therefore
interpret the [Ca iv] line as not directly related to the [Ar ii] emission. A
connection with the compact object is possible, because the emission region
and line width are much smaller than observed for the general ejecta, which
reaches velocities $\sim 5000\rm~{}km~{}s^{-1}$ close to the ER, while the
ionisation state is higher (?). The region just north of the center also shows
enhanced continuum emission at $3.5$ µm (Fig. S6), whereas there is no point
source in the continuum at the location of the Ar emission.

###  1.2 Spectral modeling

####  1.2.1 Photoionization model

To model the photoionization by the radiation from a PWN or CNS we update an
earlier photoionisation code (?). We now include all important ionization
states of H, He, C, N, O , Ne, Na, Mg, Al, Si, S, Ar, Ca, and Fe. In (?) Ar
was not included in the models, and the Si-S-Ar-Ca zone was not discussed.

The ionization structure is set (to first order) by the the ionization
parameter (?),

$\xi=\frac{L_{\rm ion}}{n_{\rm ion}r^{2}}$ (S1)

which is proportional to the ratio of the number density of ionizing photons
to the number density of atoms. $L_{\rm ion}$ is the ionizing luminosity,
$n_{\rm ion}$ the ion number density and $r$ the distance to the source. For
discriminating between different ionizing sources the spectral shape is also
important for the fractions of the different ionization stages and the
temperature, as illustrated by the different models.

For the ionization equilibrium we include photoionization, collisional
ionization, radiative and dielectronic recombination, including updated
recombination data (?, ?, ?, ?, ?), photoionization cross sections were
updated from fits in (?) and the TOP base (?), and Auger ionization from the
inner shells (?). Charge transfer between elements heavier than H and He are
lacking, with a few exceptions, and introduce an uncertainty.

The ionizing flux is calculated in the ‘outward only’ approximation, including
diffuse emission from continuum emission and line emission (?), while the
radiative transfer in the lines is treated with the Sobolev approximation. The
line emission is calculated using detailed multilevel atoms in non-local
thermodynamical equilibriumnon-local (NLTE), with atomic data from the Chianti
Atomic Database (?), NORAD-ATOMIC-DATA base (?), NIST Atomic Spectra Database
(?) and the Atomic Line List (?). We adopt updated collision strengths for the
important fine-structure lines of [Ar ii] (?) and [Ar vi] (?), and for [Ca IV]
(?). The calculation is stopped when the temperature reaches 100 K, where no
lines below 30 $\mu$m are excited. The updated code is written in Fortran and
available on Zenodo (?), together with the output files.

####  1.2.2 Abundances and structure of the core

Models of the progenitor and nucleosynthesis of SN 1987A indicate a zero age
main sequence (ZAMS) mass of $15-20\rm~{}M_{\odot}$ (?). Fig. S7 shows the
predicted composition in the central region of an explosion of a
$19\rm~{}M_{\odot}$ ZAMS progenitor (?), adopted for the modeling. Other
models in the 15 to 20 $\rm~{}M_{\odot}$ range have similar structure and
mainly differ in the varying masses of the different zones. Close to the NS
the predicted composition is dominated by He and 56Ni, decaying into 56Fe.

Figure S7: Fractional abundances by number as a function of mass from the
center for the core of a $19\rm~{}M_{\odot}$ model (?). Ar is most abundant in
the zone between $1.8-2.1\rm~{}M_{\odot}$ from the center. The region interior
to $\sim 1.65\rm~{}M_{\odot}$ (dashed vertical line) becomes the NS.

For the O – Si – S – Ar – Ca zone between $1.9-2.1\rm~{}M_{\odot}$ in mass
coordinates, we have taken an average fractional abundances over the zone
(Table S2). Within the zone the abundances of the most abundant elements O,
Si, S, Ar and Ca, vary by approximately a factor of two (Fig. S7). For
comparison we model the line emission from from the O – Ne – C – Mg zone in
Sect. 2.6.3 , with averaged abundances between $2.10-4.0\rm~{}M_{\odot}$
(Table S2).

Table S2: Averaged fractional abundances used in the calculations. Element | Abundance by number
---|---
| O–Si–S–Ar–Ca | O–Ne–C–Mg
H | $2.06\times 10^{-7}$ | $6.78\times 10^{-9}$
He | $2.30\times 10^{-5}$ | $1.10\times 10^{-4}$
C | $1.58\times 10^{-4}$ | $4.04\times 10^{-2}$
N | $1.92\times 10^{-5}$ | $3.42\times 10^{-5}$
O | $5.80\times 10^{-1}$ | $7.45\times 10^{-1}$
Ne | $1.52\times 10^{-4}$ | $1.67\times 10^{-1}$
Na | $1.00\times 10^{-9}$ | $1.84\times 10^{-6}$
Mg | $6.46\times 10^{-4}$ | $3.12\times 10^{-2}$
Al | $1.85\times 10^{-4}$ | $2.83\times 10^{-3}$
Si | $2.26\times 10^{-1}$ | $6.79\times 10^{-3}$
S | $1.56\times 10^{-1}$ | $1.17\times 10^{-3}$
Ar | $2.74\times 10^{-2}$ | $1.47\times 10^{-4}$
Ca | $8.29\times 10^{-3}$ | $3.51\times 10^{-5}$
Fe | $1.40\times 10^{-3}$ | $4.34\times 10^{-4}$

The density in the SN core is uncertain.The inner ejecta in SN 1987A is both
clumpy and anisotropic (?, ?, ?). To get an estimate of the density we
nevertheless approximate this with a constant density sphere. For constant
emissivity the line profile is then parabolic, with flux density $F(V)\propto
1-(V/V_{\rm core})^{2}$, where $V$ is the velocity from the line center and
$V_{\rm core}$ the maximum velocity of the core. The full width velocity at
half maximum, $V_{\rm FWHM}$, is related to this by $V_{\rm core}=V_{\rm
FWHM}/\sqrt{2}$.

Observations of the [O i] $\lambda\lambda$6300 , 6364 lines during the first
500 days after explosion were modeled with a maximum velocity $\sim
1700\rm~{}km~{}s^{-1}$ (?). Similarly, observations of these and other lines
in 1995 were modeled with a constant emissivity profile to $\sim
1700\rm~{}km~{}s^{-1}$, steeply decreasing at higher velocities. A fit with a
parabolic line profile to the [O i] $\lambda$6300 line from 1999 (?) gives a
good fit up to $\sim 2200\rm~{}km~{}s^{-1}$, but shows an extension above the
parabolic profile at higher velocity, as expected from the 3D observations
above. From this model we find $V_{\rm FWHM}\approx 3500\rm~{}km~{}s^{-1}$,
implying $V_{\rm core}\approx 2500\rm~{}km~{}s^{-1}$. A similar fit to the [Fe
i] $\lambda~{}\ 1.443$ µm, [Fe ii] $\lambda~{}1.534$ µm, and [Fe ii]
$\lambda~{}\ 25.98$ µm lines (?) give $V_{\rm core}\approx
2700\rm~{}km~{}s^{-1}$. As an average we adopt $V_{\rm core}=2200\pm
500\rm~{}km~{}s^{-1}$.

The ejecta mass inside the oxygen core of a 19 $\rm~{}M_{\odot}$ progenitor is
$\sim 3.0\rm~{}M_{\odot}$, excluding the NS (?). For a uniform oxygen core
density with the above mass and $V_{\rm core}=2200\rm~{}km~{}s^{-1}$, this
corresponds to a density of $\sim 9.66\times 10^{-20}\ \rm g~{}cm^{-3}$ at 35
years, or a number density of ions of $n_{\rm ion}=2.63\times
10^{3}(A/22)^{-1}\ \rm~{}cm^{-3}$, where $A$ is the mean atomic weight, $A\sim
30$ in the inner Si zone, $A\sim 22$ in the O – Si –S –Ar –Ca zone and $A\sim
16.8$ in the O – Ne – C – Mg zone.

However, both observations and hydrodynamical multidimensional simulations (?)
have shown that clumping and filaments produce considerably higher density in
the emitting gas. The line ratio of the [O i] $\lambda\lambda~{}6300,6364$
doublet in SN 1987A during the first year indicated a volume filling factor of
$\sim 0.1$ in the oxygen core (?, ?), implying a factor ten higher density
compared to the constant density estimate above. For the iron-dominated region
there is evidence of expansion from the radioactive heating during the first
months after explosion (?). Modeling of the [Fe ii] optical and infrared (IR)
lines indicated filling factors in the range $0.2-0.5$, despite constituting
only $2.4\%$ of the core mass (?, ?). The low density of the Fe dominated
regions in combination with the strong emission lines from the radioactively
powered ejecta could be responsible for the weak narrow [Fe ii] lines. As
mentioned above, the 3D IFU observations of the [Si i]/[Fe ii]
$\lambda~{}1.64$ µm line (?, ?), the [Fe i] $\lambda~{}\ 1.44$ µm line (?),
and lines of CO and SiO (?), indicate a highly asymmetric and clumpy
distribution. This agrees qualitatively with hydrodynamical simulations of the
explosion (?), where the combination of an early reverse shock and 56Ni
heating produces a clumpy structure with high-density regions mixed with low-
density cavities. If a PWN is present, this would produce further clumping.
Both analytical calculations (?) and hydrodynamical simulations (?) predict
that the shell swept up by the pulsar wind is unstable and can lead to
filamentation. The strongest effects are in the break-out phase when the
bubble shell accelerates down the steep density gradient. Based on the
observed filling factor above, $\sim 0.1$ , for the oxygen core above, we
adopt a density ten times the uniform core density in Fig. 4. Models with
higher filling factors are discussed in Sect. 2.6.1.

There is also evidence for mixing in velocity: H is mixed in from the envelope
to the core (?, ?), while Fe is mixed in the opposite direction (?, ?). An
example of a clumpy, anisotropic and strongly velocity mixed ejecta is Cas A,
where knots and filaments formed in the explosion from different
nucleosynthesis zones are completely mixed in velocity (?, ?).

The inhomogeneous, anisotropic structure close to the region where the central
emission originates, and also the presence of dust, means that only a fraction
of the ionizing luminosity is converted into optical and IR emission lines.
The rest either escapes the center and is absorbed at larger radii in the
ejecta, or is converted to thermal radiation by dust. To model this we
introduce a parameter representing the covering factor, $f$, of the line
emitting regions as seen from the ionizing source, $f=\Omega/4\pi$, where
$\Omega$ is the solid angle. The predicted line luminosities are a factor $f$
lower compared to the spherical case. $f$ is a free parameter that we fit to
the data so that the [Ar ii] $6.985$ µm line agrees with the observations.
This determines the luminosities of the other lines in the model from the
relative line fluxes.

####  1.2.3 Shock models

We calculated models of line emission from shocks, with a setup corresponding
to a PWN shock (Sect. 2.2). The code is a major update of (?, ?), similar to
(?). We solve the steady state hydrodynamic equations, including cooling and
radiative emission in the same way as for the photoionization models, but
including the time dependent terms in the ionization and thermal balance. The
atomic data is the same as used for the photoionization models. We do not
include any ionization precursor in the model, as in (?), but have tested the
sensitivity to the assumed state (here assumed to be in the third ionization
stage, e.g. [O iii] etc.), finding that different assumptions have limited
effects on the resulting line luminosities. An inherent limitation to these
one-dimensional models is that they are subject to instabilities (?, ?). If
the real geometry is not spherical or planar, a range of shock velocities
would be produced by oblique shocks. Therefore we regard these models as
indicative only. The code is written in Fortran and available on Zenodo (?),
together with the output files.

###  1.3 Relation to the dust peak

Dust could thermalize emission from a compact object if the optical depth is
high at optical-IR wavelengths (?, ?, ?). A submillimetre map of the dust
temperature in SN 1987A (?) shows an increase in the temperature at a position
slightly north of the center. The position of the [Ar vi] centroid is located
at $38\pm 22$ mas east and $31\pm 22$ mas south of the center (section 1.1.1).
The location of the dust peak was only quoted as the position of the brightest
pixel, without an uncertainty (?) (though the offset of this pixel of 72 mas
east and 44 mas north of the center was stated to be 3-5 times the total
alignment uncertainty). This prevents us from assessing the significance of
the offset between the dust temperature peak and Ar emission site. If the dust
peak is due to heating by a CNS, the CNS must be very close to the location of
the peak (?). The offset between that position and the Ar emission we observe
is difficult to explain in this scenario.

The centroid of the [Ca iv] line at $3\pm 10$ mas east and $114\pm 15$ mas
north from the center (section 1.1) is closer to the dust peak than the [Ar
vi] centroid, but any possible association between [Ca iv] and the dust
emission cannot be quantified due to the lack of uncertainty on the dust peak.

##  2 Supplementary text

###  2.1 Ejecta lines and comparison to other young supernova remnants

We only detect these the narrow emission line components from S, Ar and
possibly Ca. The chemical composition in the center of the ejecta is likely
dominated by elements resulting from oxygen and silicon burning. At the low
temperature of the high metallicity inner ejecta ($\lesssim 2000$ K, see
below), ions which have low excitation temperature near infrared (NIR) and MIR
transitions are predicted to dominate the cooling and therefore the observed
spectrum.

Both theoretical models and observations of other SNRs support this
interpretation (see below). The models predict that Si is the most abundant
element in this region; the strongest of the lines from this element are the
[Si i] 1.607, 1.645 µm and [Si ii] 34.81 µm. However, Si i is easily ionized,
while the [Si ii] line is outside our observed wavelength range. Although the
[S iii] 18.71 µm and [S iv] 10.51 µm lines dominate the cooling of the Si – S
– Ar – Ca zone, their wavelengths are in a region of high dust absorption
(main text) and are therefore expected to be weak.

Several other young supernova remnants (SNRs) have been observed to show
strong high-ionization lines with abundances characteristic of the inner
ejecta, powered by synchrotron radiation from a PWN or thermal emission from
shocks. The most extensively studied case is SNR Cas A, where fast-moving
knots are dominated by [S ii], [Ar iii] and [Ca ii], mixed with O lines with
varying strengths, indicating ejecta that has undergone different degrees of O
burning (?). MIR Spitzer observations of SNR Cas A have shown strong lines of
[Ni ii] $6.636$, [Ar ii] $6.985$, [Ar iii] $8.991$, [S iv] $10.51$, [Ne ii]
$12.81$, [Ne v] $14.32$, [Ne iii] $15.56$, [Fe ii] $17.94$, [S iii] $18.71$,
[O iv] $25.89$, [Fe ii] $25.99$, [S iii] $33.48$, and [Si ii] $34.82$ µm (?).
Because SNR Cas A does not contain a PWN, these lines must be excited by
either shocks or the ionizing radiation produced by the reverse shock.

Both SNR G54.1+0.3 and SNR 0540-69.3 contain a PWN, with strong synchrotron
emission. In the case of SNR G54.1+0.3 MIR spectra show ejecta lines of H2
12.3, 17.0 µm and atomic lines from [Ar ii] $6.985$, [S iv] $10.51$, [Ne ii]
$12.81$, [Cl ii] $14.37$, [S iii] $18.71$, [Fe ii] $25.99$ , [S iii] $33.48$,
and [Si ii] $34.82$ µm (?). The MIR $10-37$ µm spectrum of SNR 0540-69.3 shows
similar high-ionization lines of [O iv], [Ne ii-iii], [Si ii], [S iii-iv] and
[Fe ii] (?).

In all three cases the same set of atomic lines are seen in these remnants,
with few variations. All three have a strong [Ar ii] line and strong lines
above $\sim 8$ µm from highly ionized heavy elements. While also SN 1987A has
a strong [Ar ii] line, it differs from these SNRs by having weak or abscent
lines above $\sim 8$ µm. We discuss a possible explanation for this in the
main text and in Sect. 2.4.

###  2.2 Origins of the narrow ejecta lines

The observed displacement, $\sim-250\rm~{}km~{}s^{-1}$, of the [Ar ii] and [Ar
vi] lines from the systemic velocity of SN 1987A, is much smaller than the
general ejecta velocity, which reaches $\sim 5000\rm~{}km~{}s^{-1}$ close to
the ER and $\sim 10,000\rm~{}km~{}s^{-1}$ at the reverse shock out of the ER
plane (?). Together with the spatial distribution, this indicates that the Ar
emission must have an origin close to the center. If the emission was coming
from a spherically symmetric shell, this would have resulted in a boxy, flat-
topped line profile, with a minimum and maximum velocity equal to the velocity
of the shell, inconsistent with the observations. Instead we observe a
displaced line with a FWHM width $\sim 122\rm~{}km~{}s^{-1}$ of the [Ar ii]
line.

There are several possible scenarios for the excitation of the narrow ejecta
lines, including emission from a PWN produced by a NS in the center, thermal
emission directly from a CNS, accretion onto a compact object, shock
excitation in the inner ejecta, excitation by radioactive 44Ti, a surviving
companion star, dust reflection of the narrow line emission from the ER
collision or emission from a reverse shock. We investigate each of these
possibilities below.

Dust reflection from the ER is possible because [Ar ii] is seen both in the ER
and central ejecta. Many silicates scatter light mainly below $\sim 8$ µm, but
are strongly absorbing above this wavelength (Sect. 2.4). However, if this was
the origin of the [Ar ii] lines we would expect similar reflections for the
other strong lines from the ER below $\sim 8$ µm, which we do not observe.
Fig. S4 shows that a nearby H i line at $\sim 4.376$ µm is stronger than the
[Ar vi] line in the ER spectrum but does not show a central emission.
Equivalent situations occur for many other strong lines in the NIRSpec range.
In addition, the spatial FWHM of the central component [Ar vi] line is
$\lesssim 50\%$ of that of the ER and is blueshifted by over 200 km s-1 from
the ER emission.

More than $\sim 1500$ days after explosion, radioactive 44Ti dominates energy
input to the inner ejecta (?). Modeling of the optical spectra (?) and
observations of hard X-rays (?, ?) have shown that $\sim 1.5\times
10^{-4}\rm~{}M_{\odot}$ of 44Ti was produced in the explosion, which undergoes
radioactive decay emitting gamma-rays and positrons. Most of the gamma-rays
escape, while positrons are mostly trapped in the ejecta, where their
annihilation provides an energy input of $\sim 1.2\times
10^{36}\rm~{}erg~{}s^{-1}$ (?). If a large fraction of this energy was
deposited in the Ar region of the ejecta, it could power the observed emission
lines. However, the narrow line width indicates the emission must arise close
to the center (see above). This is unlike the broad lines from the ejecta,
such as the strong [Fe ii] $\lambda~{}1.644$ µm and $\lambda~{}25.99$ µm
lines, which have a velocity of at least 2500 $\rm~{}km~{}s^{-1}$ (?, ?).
Because of the low energy of the secondary electrons produced by the positron
thermalization (?), this process is unlikely to produce ions higher than
singly ionized ions from the non-thermal positrons (?, ?), so would produce
very little [Ar vi] emission. While the broad components seen in the [Fe i]
and [Fe ii] lines arise from 44Ti and X-ray input from the ER collision (?),
we therefore conclude the 44Ti mechanism is unlikely to be the source of the
narrow Ar components. This was also the conclusion from a detailed study of
the heating of the submillimter dust ’blob’ (?).

Another possibility is that the X-ray input from the ejecta and circumstellar
medium interaction could be responsible for the high ionization lines from the
center. Models predict that a small fraction of the hard X-rays penetrate all
the way to the central regions (?). However, most of the luminosity in the
X-rays from the shocks is emitted at energies below $\sim 1$ keV, so absorbed
in the outer parts of the ejecta. The effect of the X-ray input can be seen in
the [Fe ii] 5.340 µm line, which is strong especially the south-west region [
(?), figures 15 and 16]. However, there is little emission in this line from
the central region, indicating that X-rays from the ER do not penetrate to the
centre of the ejecta. It is therefore unlikely that external X-ray flux could
excite the narrow high-ionisation emission lines we observe in the center,
although it could be responsible for some [Ar ii] emission from the region
close to the ER.

A reverse shock is known to be present in SN 1987A, but only close to the ER,
with a velocity of $\gtrsim 4000\rm~{}km~{}s^{-1}$ (?). The inner ejecta must
have been traversed by a reverse shock, but only within a few days of the
explosion. Simulations have predicted the formation of an additional reflected
shock, but it is expected to dissipate within a year (?). We therefore exclude
the possibility of reverse shock excitation.

A surviving companion star is constrained to $\lesssim 3.7\rm~{}M_{\odot}$ if
on the main sequence (?), which corresponds to an A to F star. These produce
few ionising photons and have soft spectra, insufficient to explain the [S
iii], [S iv] and [Ar vi] emission.

###  2.3 Scenarios with a compact object

Assuming no pulsar activity, the minimum bolometric luminosity of a central NS
is that of a young CNS, which has a temperature of $\gtrsim 10^{6}$ K and
luminosity $\gtrsim 10^{34}\rm~{}erg~{}s^{-1}$ (?), sufficient to ionize a
substantial region of the inner core. We discussed this possibility in the
main text and give more details below.

For the scenario with a PWN, discussed in the main text, we assume the PWN
produces X-ray emission with a power law spectrum with
$L(\nu)\propto\nu^{-\alpha}$ with power law index $\alpha=1.1$, similar to the
Crab PWN, normalized to a total luminosity above 13.6 eV less than the
inferred upper limit of $10^{36}\rm~{}erg~{}s^{-1}$ (?). From broad band X-ray
observations of several PWN spectra there is a significant break in the power
law at energies $4-14$ keV (?, ?), consistent with a cooling break (?).
However this has little effect on the line emission (Sect. 2.6.2).

If a PWN is present, it produces synchrotron emission and a slow, dynamic
shock that propagates into the ejecta (?). This shock is likely to be cooling,
emitting line emission in the IR to soft X-ray range, and can therefore
contribute to the emission from the inner ejecta. Because the pressure is
approximately constant behind the shock, the cooling produce a high density,
thin shell between the ejecta and PWN. The expansion velocity of this shell,
$V_{\rm shell}$, depends on the energy input by electromagnetic radiation,
relativistic particles and magnetic fields from the pulsar to the PWN. We
expect that to be close to the spin-down luminosity of the pulsar, $L_{\rm
sd}$.

To estimate $L_{\rm sd}$ we note that the ratio of the synchrotron luminosity
to the spin-down luminosity is $\sim 28\%$ for the Crab PWN (?). An upper
limit to the bolometric luminosity from a compact object in SN 1987A is $\sim
10^{36}\rm~{}erg~{}s^{-1}$ (?). Using this, we estimate a total spin down
luminosity of $L_{\rm sd}\lesssim 3\times 10^{36}\rm~{}erg~{}s^{-1}$. For an
oxygen core mass $M_{\rm core}\approx 3\rm~{}M_{\odot}$ and a core velocity
$V_{\rm core}\approx 2200\rm~{}km~{}s^{-1}$ the kinetic energy of the core is
$\sim 7.2\times 10^{49}$ ergs. For a constant density core we then estimate
$V_{\rm shell}$ [ (?), their equation 2.10]

$V_{\rm shell}=275\left(\frac{L_{\rm sd}}{3\times
10^{36}\rm~{}erg~{}s^{-1}}\right)^{0.2}\left(\frac{V_{\rm core}}{2200\
\rm~{}km~{}s^{-1}}\right)^{0.3}\left(\frac{M_{\rm
core}}{3\rm~{}M_{\odot}}\right)^{-0.2}\left(\frac{t}{35\ \rm
years}\right)^{0.2}\rm~{}km~{}s^{-1}$ (S2)

where $t$ is the time after the SN explosion. Alternatively, for a core with a
density decreasing with radius, $r$ as $\rho\propto r^{-1}$:

$V_{\rm shell}=148\left(\frac{L_{\rm sd}}{3\times
10^{36}\rm~{}erg~{}s^{-1}}\right)^{0.25}\left(\frac{V_{\rm core}}{2500\
\rm~{}km~{}s^{-1}}\right)^{0.5}\left(\frac{M_{\rm
core}}{3\rm~{}M_{\odot}}\right)^{-0.25}\left(\frac{t}{35\rm\
years}\right)^{0.25}\rm~{}km~{}s^{-1}.$ (S3)

The PWN shock velocity into the expanding ejecta is in the two cases $V_{\rm
shell}/6$ and $V_{\rm shell}/5$, respectively. The resulting shock velocities
are then $46\rm~{}km~{}s^{-1}$ and $30\rm~{}km~{}s^{-1}$, respectively. The
low shock velocity is because of the relative velocity of the shell and
expanding ejecta. The uncertainties in the parameters are large, in particular
these estimates assume a spherically-symmetric uniform core, which is a
unlikely to be the case (?). Nevertheless, we conclude that a PWN shell can
produce line emission at low velocity.

MIR observations of SNR G54.1+0.3 show fine-structure emission lines from [Ar
ii] and other lines similar to other SNRs (?, ?). This emission was
interpreted as coming from a slow shock, $\sim 25\rm~{}km~{}s^{-1}$, from the
expanding PWN into the ejecta, producing a density enhancement. A similar
model could in principle also explain our observations. An expanding
anisotropic PWN, as observed in SNR G54.1+0.3, referred to as a ’jet’ (?),
could explain the fact that we only observe narrow blueshifted emission lines.
This can, however, also be due to internal ejecta dust absorption (see below).

###  2.4 Dust absorption and scattering

Observations of dust in young, oxygen-rich SNRs, including Cas A (?) and SNR
G54.1+0.3 (?, ?) have been interpreted as showing different dust properties.
Modelling of the dust SED used mainly Mg0.7SiO2.7 and a secondary component of
other silicates, carbon or alumina dust (?). Alternatively, it has been
proposed that the $20$ µm peak in SNR G54.1+0.3 is due to SiO2 grains, while
the $10$ µm peak is produced by SiC grains and polycyclic aromatic
hydrocarbons (PAHs) (?). The latter composition also matches the dust emission
in Cas A (?).

Other observations of Cas A have indicated different dust properties in the
different abundance regions in the ejecta (?). The Ar ii and Ar iii regions in
the ejecta have strong silicate peaks, which were attributed to Mg0.7SiO2.7
(?).

Although there is degeneracy in the models and conflicting interpretations of
the dust spectra, a common result is that silicates are invoked to explain the
$10-20$ µm region. Because dust formation occurs at high temperatures in the
cooling ejecta, a large fraction of the silicates can be in crystalline dust,
as is seen in e.g., the outflows of asymptotic giant branch (AGB) stars and in
comets (?, ?). Whether the silicates are crystalline or amorphous depends on
the cooling time scale compared to the expansion time scale. Fast cooling
produces amorphous dust while very slow cooling favours crystalline.

Independent of the exact dust composition, size and structure of the grains, a
common feature of silicates is a large drop of the absorption shortwards of
$8-9$ µm. Depending on the composition and crystallinity of the dust this can
be more or less sharp. We invoke this property to explain the different MIR
spectra of SN 1987A compared to other young SNRs, discussed in Sect. 2.1, as
well as the models discussed in this section.

Fig. S8 shows the optical depth as a function of wavelength for several
silicates, both in amorphous and crystalline forms, produced in laboratory
experiments by different methods and therefore having different structures.
The absorption is calculated in the Rayleigh limit, assuming grains are small
compared to the wavelength and with the continuous ellipsoidal distribution
(CDE2) (?, ?). For comparison, we have normalized all optical depth curves to
$\tau_{\rm abs}=10$ at $10$ µm. Optical constants for amorphous Mg2SiO4 have
been obtained using different techniques: laser ablation (?, ?) was used to
produce the sample, or the sol-gel method (?). Both methods lead to similar
absorption curves above $\sim 15$ µm, but below the 10 µm peak the sol-gel
method gives a flatter slope. Therefore even for the amorphous case the
optical properties depend on the structure for the same composition. A similar
behaviour is seen for enstatite, MgSiO3, where the MgSiO3 data are from glass
(?) and sol-gel data (?). Observations of several AGB stars show a steeper
slope below the 10 µm peak than amorphous silicates produced by the sol-gel
method (?). Crystalline enstatite (?) has a similar steep slope below 10 µm,
but with considerable structure in the absorption above the peak.

Figure S8: Optical depth of absorption as function of wavelength for different
silicate compositions and structures. The optical depths have been normalized
to $\tau_{\rm abs}=10$ at the $10$ µm peak. There is a drop in $\tau_{\rm
abs}$ below $\sim 8$ µm in all cases. The wavelength of the [Ar ii] $6.985$ µm
line is shown by the vertical, dashed line. Optical constants of amorphous
Mg2SiO4 obtained by laser ablation are from (?, ?) (solid line) and by the
sol-gel method from (?) (black dashed line). For enstatite, MgSiO3, data for
glass are from (?), for crystalline from (?) (blue dashed line) and by the
sol-gel method from (?) (magenta dash-dotted line).

At wavelengths shorter than $\sim 8$ µm dust scattering usually dominates
absorption by silicates (?), especially at the shorter wavelengths. The
scattering and absorption properties depend strongly on the size of the
grains, as well as on composition, especially if Fe is included, e.g. MgFeSiO4
(?). The grain size is also important, with larger grains increasing the
scattering. Theoretical calculations indicate a distribution that peaks at
$\sim 0.1$ µm, extending to $\sim 1$ µm (?).

Pure scattering does not thermalize the photons, but produces broader lines
from the expanding ejecta and more spatially extended emission that is
difficult to separate from the general ejecta emission. If there is optically
thick dust in the core, mixed with the line emitting clumps, emission from the
far side relative to the observer (redder velocities) could be absorbed, as
was first proposed for SN 1987A as an explanation of the blue shift of the
emission lines during the first years (?). This, together with clumping, could
explain the blueshift of the Ar lines we observe. This requires the dust to
have some absorption below 8 µm, which is the case for olivine, MgFeSiO4 (?).

###  2.5 Detailed results for the photoionization models in the main text

Figure 4 shows the results of our ’standard’ photoionization models for the
PWN and CNS cases. For these models with $n_{\rm ion}=2.6\times
10^{4}\rm~{}cm^{-3}$, we find that $\xi\approx 0.5$ gives the best fit to the
line ratios for both models. Specifically, we used for the PWN model $L_{\rm
ion}=5\times 10^{34}\rm~{}erg~{}s^{-1}$ and $r=2.8\times 10^{15}$ cm and
$L_{\rm ion}=3\times 10^{35}\rm~{}erg~{}s^{-1}$ and $r=5.6\times 10^{15}$ cm
for the CNS model. However, only the flux, $L_{\rm ion}/4\pi r^{2}$, affects
the results.

####  2.5.1 Temperature and ionization structure

Figure S9 shows the temperature and ionization structure of selected elements
for the two cases. The ionization structure of the CNS and PWN models are
quite different as a result of the different ionizing spectra. The blackbody
like CNS spectrum leads to diffuse zones with a large range of ionization
stages at the same radius, while the steeper PWN spectrum produces more
discrete zones. The high metallicity leads to a general lower temperature
compared to a gas with solar abundances, $\lesssim 10^{4}$ K in the high
ionization zone, decreasing to a few hundred K in the partially ionized zone.
A consequence of this is that most of the cooling is dominated by NIR and MIR
fine structure lines of O i-iii, Si i-ii, S i-iv, Ar ii-v, with low excitation
temperatures.

Figure S9: Photoionisation models of the CNS (A) and a PWN (B) cases. In the
top panel the temperature, $T_{\rm e}$, is shown with the red curve and
electron fraction, $X_{\rm e}$, with the blue. The different ionization
fractions, $X$, are color and line coded as shown by the legend. There are
distinct ionization zones for the PWN, which overlap in the CNS case.

The thickness of the ionized region is determined by the number of ionizing
photons, electron density and recombination rates of the dominant ions, and
the ionization parameter, so $\Delta r\propto\xi$. Figure S9 shows the
thicknesses in both models are $\lesssim 10^{16}$ cm, which is too small to
resolve in our observations.

Compared to photoionization regions with solar abundances, the higher
abundance of heavy elements, which are efficient coolants, results in a lower
temperature, $\lesssim 1000$ K, in most of the line emitting region, even if
the ionization is higher. The result is that few of the usually strong UV
resonance and semi-forbidden lines seen in objects with solar abundances can
be excited, and the cooling and emission is instead dominated by MIR and far-
IR lines from heavy elements. Only close to the inner boundary, and for a high
ionization parameter, is there a zone with temperature $\gtrsim 10^{4}$ K
where UV lines are excited.

In addition to the lines discussed in the main text, the photoionisation
models predict other strong lines at shorter and longer wavelengths. Longwards
of the NIRSpec and MRS range the models predict strong fine structure lines,
including [Si ii] 34.81, [S iii] 33.49, [O iii] 51.81, 88.33 µm. However, no
observations are available in this range. In the optical, the models predict
strong [O iii] 4959, 5007 Å, [Ar iv] 4712, 4741 Å, [Ar v] 6435, 7006 Å, [Ca
ii] 7291, 7323 Å, [Ar iii] 7136, 7751 Å , with line ratios that differ between
the models. However, dust scattering could affect these predictions.

####  2.5.2 Lines in the NIR

In the NIRSpec range [Si i] $1.099,1.607,1.645$ µm lines, [S i] $1.033$ µm and
[S iii] $0.9533$ µm lines are predicted. These lines are stronger in the CNS
models, due to the extended Si i zone with a temperature $\gtrsim 10^{3}$ K
(Fig. S9), which is high enough for collisional excitation of these lines. The
PWN models have more well-defined ionization zones, and a much lower
temperature in the Si i zone, resulting in weak [Si i] lines. We do not find
any narrow lines at these wavelengths in the NIRSpec spectrum, which may speak
in favour of the PWN model. The limited spectral resolution of NIRSpec and
overlapping broad ejecta lines, however, makes it difficult to distinguish
these. Dust scattering can also broaden the lines both in wavelength and
spatially. The low ionization potential of Si i, 8.15 eV, means that any
radiation below 1520 Å can ionize Si i. This effect is seen in the PWN models,
which have a strong continuum at these wavelengths. The strong UV radiation
from the ER shocks, including Lyman $\alpha$ (?), will have the same effect,
unless dust extinction to the center absorbs the radiation from the ER. There
is also a strong UV field at these energies from the ejecta itself in the form
of Lyman $\alpha$ and two-photon emission from hydrogen and helium powered by
44Ti decay (?).

The main difference between the models and observations in the NIRSpec range
are the strong [Ca iv] 3.207 µm and [Ca vi] 4.159 µm lines. These lines are in
the models an order of magnitude stronger than the [Ar vi] line in the same
wavelength range. The [Ca iv] 3.207 µm is observed from the central region
(Fig. S5), but is coming from a different location and is redshifted. This
discrepancy could be caused by uncertainties in atomic data, in particular
photoionization cross sections, collision strengths and recombination rates.
The Ca iv and Ca v ionic fractions are at least an order of magnitude larger
than the Ar vi abundance (Fig. S9). These effects could lead to overestimates
of the Ca line strengths. The relative abundances between calcium and argon
depend on the ZAMS mass (?) and vary within the O–Si–S–Ar–Ca zones (Fig. S7),
but only within a factor of about two. The relative Ar to Ca abundances in
this zone for models between 15 to 19 $\rm~{}M_{\odot}$, vary by a similar
amount (?). Adopting a 15 $\rm~{}M_{\odot}$ model gives similar results as for
the 19 $\rm~{}M_{\odot}$ model. However, there could be systematic
uncertainties in the treatment of both the pre-SN burning and the explosive
oxygen burning in these one-dimensional models, which could produce incorrect
different abundance ratios (?).

An alternative explanation for the high Ca line luminosities in the model
could be the neglect of depletion of Ca to dust. Ca is not included in the
dust models in (?) nor (?), so no quantitative predictions are available.
However, calciosilia, like wollastonite, CaSiO3 and diopside, MgCaSi2O6 (?),
could be produced by the dust formation in the inner ejecta, , lowering the Ca
abundance in the gas phase. The same applies to the NIR [Si i] lines, which
are predicted to be strong in the CNS model but not observed. The large mass
of silicates indicated by the dust observations, and formation of molecular
SiO (?, ?), could lower the gas-phase abundance of atomic Si.

####  2.5.3 The covering factor

The ionization structure, and resulting line ratios, are determined by the
ionization parameter, $\xi=L_{\rm ion}/n_{\rm ion}r^{2}$, for the same
spectral shape, $L(\nu)$. Because we do not know how large the solid angle of
the emitting gas is from the ionizing source we normalize the [Ar ii] $6.985$
µm luminosity in the model to the observed value by a covering factor
$f=L_{\rm obs}({\rm[Ar~{}\textsc{ ii}]})/L_{\rm
model}({\rm[Ar~{}{\textsc{ii}}]})$. In the models $\sim 20\%$ of $L_{\rm ion}$
emerges in the [Ar ii] line, so $f\approx 5\ L_{\rm obs}$([Ar ii])$/L_{\rm
ion}$. The ionization parameter is determined by the relative line
luminosities, primarily the [Ar i vi]/[Ar ii] ratio. Therefore we can write
$f\approx 5\ L_{\rm obs}$([Ar ii])$/\xi n_{\rm ion}r^{2}$, with $\xi\approx
0.2-0.4$. This means that $f$ is determined by the distance between the
ionizing source and emitting gas which is not known, except for an upper
limit.

The luminosity, distance to the ionizing source and density of the PWN model
in Fig. 4 were chosen to be $L_{\rm ion}=5\times 10^{34}\rm~{}erg~{}s^{-1}$,
$r=2.8\times 10^{15}$ cm and $n_{\rm ion}=2.6\times 10^{4}\rm~{}cm^{-3}$ for
the PWN model, resulting in a covering factor $f=0.022$. To produce the same
ionization parameter for the CNS model with the same density and a luminosity
$L_{\rm ion}=3\times 10^{35}\rm~{}erg~{}s^{-1}$, indicated by the CNS models
(?), the distance has to be $5.6\times 10^{15}$ cm, which results in
$f=0.005$. For a lower NS luminosity the scaling is the same as above,
$f\propto L_{\rm ion}^{-1}$, but the luminosity is more constrained to
$\gtrsim 3\times 10^{34}\rm~{}erg~{}s^{-1}$ (?) than for the PWN case.

These low covering factors indicate that only a fraction of the total ionizing
luminosity is converted into the observed lines, and/or a lower luminosity and
distance. A small covering factor is consistent with other aspects of our
observations. The observed lines observed are shifted with respect to the rest
velocity and also narrow in velocity, which shows that only a small solid
angle is responsible for the observed luminosity. The total luminosity of the
line emission is less than $\sim 10^{33}$ erg s-1, which is at least two
orders of magnitude less than expected for a CNS (?, ?). The difference may be
due to line emission, similar to the observed, but absorbed by optically thick
dust in the NIR. Alternatively, the continuum flux from the ionizing source
could be directly thermalized by dust or absorbed by bound-free absorption in
the ejecta. This would result in broad, low ionization lines, which would not
be distinguished from the radioactively and X-ray powered lines, like the [Si
i] and [Fe ii] lines in the NIRSpec and MRS ranges (?, ?). In the estimate of
the observed total line emission above we have only included the narrow
component of the lines in Fig. 3, while the broader component also contributes
to the total [Ar ii] luminosity.

###  2.6 Alternative PWN models

####  2.6.1 Sensitivity to the core density

As discussed above, the density is uncertain because of clumping and
filaments. Our models in Fig. 4 assumed a factor of ten higher density than
the average in the core. To check the sensitivity to this we calculated a PWN
model with the average density, $n_{\rm ion}=2.6\times 10^{3}\rm~{}cm^{-3}$.
To match the observed [Ar vi] /[Ar ii] ratio, the ionization parameter, $\xi$,
has to be within $\sim 50\%$ of that of the the models with the higher
density. This requires the ionizing luminosity to be reduced by the same
factor as the density, or alternatively the distance increased by the square
root of the density, or $1\times 10^{34}\rm~{}erg~{}s^{-1}$. The results of
this model are shown in Fig. S10.

Figure S10: PWN photoionization model without clumping. Same as Fig. 4B but
assuming a ten times lower density, corresponding to a constant density core.

Comparing Figures 4 and S10, we find that the [Ar vi] /[Ar ii] ratio is the
same, as was expected. While qualitatively most of the MIR line luminosities
relative to the [Ar ii] line are similar, quantitatively the [S iii] $18.71$
µm line is considerably stronger in the model, caused by the more extended S
iii zone in the low density model. The fairly strong [S iii] lines at
$0.90-0.95$ µm in the high density model are absent in the low density model,
which is a result of a lower temperature in the S iii zone. Similarily, the
weaker [S iv] $10.52$ µm line is a result of the lower temperature in the S iv
zone.

The main effect of the density in these models is that collisional de-
excitation of the NIR and MIR fine-structure lines becomes important as the
density approaches the critical density of these lines, $n_{\rm crit}=A_{\rm
u,l}/C_{\rm u,l}$, where $A_{\rm u,l}$ is the radiative transition rate and
$C_{\rm u,l}$ is the collisional de-excitation rate between upper level $u$
and lower level $l$. Typically $n_{\rm crit}\approx 10^{4}\rm~{}cm^{-3}$ for
these lines. Collisional de-excitaion results in a lower cooling rate and
therefore a higher temperature as the density increases towards $n_{\rm
crit}$, explaining the weaker optical and NIR lines, with high excitation
temperatures, in the high density models.

####  2.6.2 PWN models with different ionizing spectra

X-ray observations have been interpreted (?) as indicating a non-thermal steep
power law with $L(\nu)\propto\nu^{-\alpha}$ with $\alpha=1.8$ in hard X-rays
from a PWN [however the same observations have also been interpreted as
indicating thermal emission (?)]. A bolometric limit is 138 solar luminosities
or $5.4\times 10^{35}\rm~{}erg~{}s^{-1}$ for a compact object in SN 1987A (?).
In (?) X-ray observations between 10-20 keV are fitted by a steep spectrum
with $\alpha=1.8$ and 10-20 keV luminosity $3.0\times
10^{34}\rm~{}erg~{}s^{-1}$. Assuming that the $\alpha=1.8$ spectrum applies to
13.6 eV, the luminosity between 13.6 eV to 20 keV is $1.4\times
10^{37}\rm~{}erg~{}s^{-1}$, i.e., a factor of more than twenty higher than the
limit above. Extending into the UV and optical would increase this
discrepancy. To comply with the power law and luminosity from (?) and the
limit above, there must therefore be a break to a flatter spectrum at soft
X-ray energies.

We calculated two PWN models with this power law and luminosity between
$10-20$ keV (?). The density and ionization parameter are the same as the PWN
model in Fig. 4. To produce a luminosity of $3.0\times
10^{34}\rm~{}erg~{}s^{-1}$ between 10-20 keV and a continuous spectrum
requires $L_{\rm ion}\gtrsim 5\times 10^{35}\rm~{}erg~{}s^{-1}$ for
$\alpha=1.1$ below the break. With $L_{\rm ion}=1\times
10^{36}\rm~{}erg~{}s^{-1}$ the break to the $\alpha=1.8$ spectrum occurs at
$E_{\rm b}=4.1$ keV, while $L_{\rm ion}=2\times 10^{36}\rm~{}erg~{}s^{-1}$
gives $E_{\rm b}=1.1$ keV.

Table S3: Luminosities relative to the $[$Ar ii$]$ 6.985 µm line for different power law breaks. No correction for dust absorption has been applied.  Ion |  1Wavelength1234 |  No break111 |  $E_{\rm b}=4.1$ keV |  $E_{\rm b}=1.1$ keV
---|---|---|---|---
| µm |  |  |
$[$S iii$]$ |  10.9530 | 0.9584 | 0.8841 | 0.9398
$[$C i$]$ |  10.9824 | 0.2444 | 0.2239 | 0.2400
$[$Ca iv$]$ |  13.2068 | 0.0966 | 0.0930 | 0.0995
$[$Ca v$]$ |  14.1585 | 0.0136 | 0.0133 | 0.0139
$[$Ar vi$]$ |  14.5292 | 0.0132 | 0.0129 | 0.0137
$[$Ar ii$]$ |  16.9850 | 1.0000 | 1.0000 | 1.0000
$[$Ar v$]$ |  17.9141 | 0.0213 | 0.0207 | 0.0221
$[$Ar iii$]$ |  18.9914 | 0.9455 | 0.9209 | 0.9698
$[$S iv$]$ | 10.5077 | 0.7319 | 0.7110 | 0.7574
$[$Ar v$]$ | 13.0680 | 0.0081 | 0.0079 | 0.0084
$[$S iii$]$ | 18.6992 | 0.7328 | 0.7160 | 0.7533
$[$Ar iii$]$ | 21.8315 | 0.0346 | 0.0340 | 0.0354
$[$Fe iii$]$ | 22.9192 | 0.0057 | 0.0056 | 0.0061
$[$Ar iv$]$ | 25.8800 | 0.0299 | 0.0292 | 0.0314
$[$Fe ii$]$ | 25.9814 | 0.0344 | 0.0324 | 0.0320
$[$S ii$]$ | 34.8100 | 0.9410 | 0.8834 | 0.8965

In Table S3 we show the line luminosities relative to the $[$Ar ii$]$ 6.9850
µm line for these two models, together with the model with no break (Fig. 4).
The change in the luminosities for these lines are less than $10\%$, and shows
that a break above $\sim 1$ keV has only marginal significance for the IR
emission line spectrum. This is because the ionization structure is mainly
determined by the spectrum below $\sim 1$ keV and the harder X-rays have
little effect for $\xi\sim 0.5$.

####  2.6.3 PWN models from the O – Ne – C– Mg zone

Figure S11: PWN photoionization model for the O-Ne-C-Mg zone. Same symbols and
colors as in Fig. 4

We calculated models with abundances from the next nuclear burning zone,
mainly O, Ne, C and Mg (Fig. S7). The cooling is dominated by C i, O i \- O
iii and Ne ii \- Ne v. The predicted spectrum (Fig. S11) is therefore very
different compared to Fig. 4, with weaker lines of Ar, but instead strong
lines of C, O, Ne, S due to the differing abundances (Table S2). This is
consistent with previous models (?). Although the relative strengths of the
[Ar ii] and [Ar ii] lines can be reproduced, this model can from the more than
two order of magnitude stronger [Ne ii], [Ne iii] and [Ne v] lines be excluded
as the origin of the observed narrow lines.

###  2.7 Results for shock models

Figure S12: Structure (A) and predicted line luminosities (B) from a 75 km s-1
PWN shock. Line styles and colours indicate ionisation states (see legend). No
dust absorption has been applied. The symbols are the same as in Fig. 4.

We calculated a model with parameters close to estimates for a PWN shock
(Sect. 2.2), $V_{\rm s}=50\rm~{}km~{}s^{-1}$ and pre-shock density $\rho_{\rm
pre-shock}=2\times 10^{3}\rm~{}cm^{-3}$. This model, however, predicts an [Ar
vi]/[Ar ii] ratio more than an order of magnitude lower than observed, and
also underpredicts the other high ionization lines. To get better agreement a
higher velocity is required. As discussed in Sect 1.2.2, the ejecta is
expected to be clumpy. If the pre-shock density is lower than the average used
above, the shock velocity will be higher, $V_{\rm s}\propto\rho^{-1/2}_{\rm
pre-shock}$ for the same bubble pressure. Fig. S12 shows the result from a
model with density, $\sim 1\times 10^{3}\rm~{}cm^{-3}$ and $V_{\rm
s}=75\rm~{}km~{}s^{-1}$, which more closely matches the observations,
especially the high ionization lines. This model does not include any ejecta
absorption, because it would completely extinguish the already low
luminosities of lines above $\sim 8\ \mu$m relative to the [Ar ii] line.

Fig. S12A shows the resulting temperature, density and ionization for this
model as a function of the column density of the ions, $N_{\rm ion}$. The
model predicts an almost constant temperature, $\sim 10^{6}$ K, in the post-
shock zone until the abrupt drop at $N_{\rm ion}\approx 10^{15}$ cm-2, which
is triggered by a thermal instability. Because the ionization time scale is
similar to the flow time scale for many ions, the ionization increases from
the pre-shock ionization behind the shock until the thermal instability. The
temperature here drops by a factor of $\gtrsim 10^{2}$ until it is stabilized
at $1000-1500$ K, where cooling is balanced by photoionization heating from
the hot region. At the same time the degree of ionization decreases and this
region is dominated by neutral elements, although there is also a substantial
fraction of singly and doubly ionized ions. Because the whole post-shock
region is close to pressure equilibrium, the density varies roughly the
opposite way to the temperature (Fig. S12A), with an ion density $\gtrsim
10^{6}\rm~{}cm^{-3}$. A magnetic pressure reduces this compression. The
thickness of the shocked region is $\sim 10^{11}$ cm, which would produce a
narrow velocity line if the shock has a limited extent in solid angle.

This model reproduces the observed relative luminosities of the [Ar vi], [Ar
iii], and [S iv] lines, within the uncertainties. However, the predicted [S
iii] 18.70 µm line is too faint by more than an order of magnitude, caused by
the rapid recombination to lower ionization stages after the thermal
instability. This may be a problem for this model. In addition, the model
predicts fairly strong [Si i] 1.607, 1.645 µm lines as a result of the high
temperature in the Si i region (Fig. S12).

This shock model predicts strong lines in the UV, including O v] 1218 Å, Si iv
1394, 1403 Å, O iv] 1401, 1408 Å, and S iv] 1405-1424 Å, and in the optical [O
iii] 4959, 5007 Å, [Ar v] 6435, 7006 Å, [Ca ii] 7291, 7323 Å, [Ar iii] 7136,
7751 Å. Dust scattering and absorption could decrease the fluxes.
Unfortunately, the latest spectral HST observations from 2017–2018 (?), do not
have sufficient S/N or spectral resolution for a meaningful comparison.

In summary, the PWN shock models have both pros and cons. They can explain the
velocity offset and narrow FWHM of the Ar lines, in particular if emission is
from an anisotropic PWN bubble [as in SNR G54.1+0.3 (?)]. The PWN scenario
does not require dust absorption to explain the lack of lines above $\sim 8\
\mu$m, although dust is known to be present in SN 1987A (?, ?). However the
observations indicate a lower density and higher shock velocities compared to
the predictions of PWN model, although a clumpy medium could explain that
difference. The most serious discrepance is the faint MIR [S iii] 18.70 µm
line in the model. In addition, the predicted [Si i] 1.607, 1.645 µm lines in
the NIR range are not seen, although this is similar to the CNS case (section
2.5), and the same caveats, connected to the observations, ionization from the
ER and ejecta and dust depletion, applies here. Given the uncertainties in the
atomic data, model parameters, geometry etc. we cannot rule out the PWN shock
scenario completely.

###  2.8 Constraints on the NS kick velocity

The position and velocity of the high ionization emission can in some
scenarios provide a limit on the NS kick, which has been predicted by core
collapse simulations (?).

The centroid of the [Ar vi] position is $38\pm 22$ mas east and $31\pm 22$ mas
south of the geometric center of the ER (Sect. 1.1.1). If we assume the ER
center was the explosion site (?), over 35 years that translates to $249\pm
146\rm~{}km~{}s^{-1}$ east and $208\pm 146\rm~{}km~{}s^{-1}$ south of the ER
center, assuming this to be the explosion site, or a combined velocity of
$324\pm 206\rm~{}km~{}s^{-1}$ on the plane of the sky. Adding the radial
velocity of $259.6\pm 0.4\rm~{}km~{}s^{-1}$ (Table1) in quadrature we get a
total velocity of $416\pm 206\rm~{}km~{}s^{-1}$ relative to the explosion
position for the [Ar ii, vi] emitting gas.

We compare this position and radial velocity to simulations of SN 1987A (?,
?), which predicted a kick velocity of $\sim-300\rm~{}km~{}s^{-1}$ in the line
of sight, and to the north of the center. While the radial velocity of the [Ar
ii] line is similar to the predictions, the position of the [Ar vi] emission
differs substantially.

The kick velocity of the NS depends on the relative positions of the NS and
the emitting gas in the different scenarios. In the PWN case the emission
would be exterior to or in the shocked shell, separating the bubble of
relativistic particles and the ejecta (Sect. 2.2). The observed radial
velocity of $\sim-260\rm~{}km~{}s^{-1}$, as well as sky position, are
therefore indicative of the line emitting site, rather than that of the NS.

For a low magnetic field and/or a long pulsar period the non-thermal energy
input from the NS would be low. In particular, the class of NSs known as
central compact objects (CCOs) in young SNRs, which represent a large fraction
of young pulsars, have magnetic fields, $10^{10}-10^{11}$ G, and spin periods
$0.3\pm 0.15$ s (?). In that case, the dynamic effects on the surrounding gas
of the NS would be small, so the emitting gas could be closer to the NS, so
the transverse and radial velocity above might approximate the kick velocity.
The FWHM of the [Ar ii] line, $\sim 122\rm~{}km~{}s^{-1}$, might then indicate
the size of the region ionized by the CNS, $\sim 1.3\times 10^{16}$ cm, which
is within an order of magnitude of the size of the ionized region in the CNS
model. This scenario could explain why only a blue shift is observed, rather
than a more boxy line profile, as would be expected for a spherically
symmetric shell without absorption.