© The Authors 2023

1 Introduction

Over the past three decades and more, supernova (SN) 1987A has proven to be one of the most important transient events ever observed. Its location only ~50 kpc from us (Pietrzyński et al. 2019), in the Large Magellanic Cloud, combined with the high spatial resolution of the Hubble Space Telescope (HST) and adaptive optics (AO) facilities, has allowed us to study this SN and its transition into a SN remnant in a uniquely detailed manner. The insights gained from this follow-up concern topics such as nucleosynthesis (e.g., Jerkstrand et al. 2011), the late-time evolution of the progenitor system, and the different power sources of the SN, from radioactive ⁵⁶Ni and ⁴⁴Ti decay to interaction (for reviews, see, e.g., McCray & Fransson 2016; Milisavljevic & Fesen 2017). Previously, SN 1987A had been observed in the near- and mid-infrared (NIR and MIR, respectively) using ground-based telescopes, the Stratospheric Observatory for Infrared Astronomy (SOFIA), and the Spitzer Space Telescope (e.g., Bouchet et al. 2006; Dwek et al. 2010; Arendt et al. 2016; Fransson et al. 2016; Larsson et al. 2016; Matsuura et al. 2022), with which it has been possible to study dust, molecules, and ejecta morphology. With the recent launch of the James Webb Space Telescope (JWST), it is now possible to increase the detail of the follow-up even further in the infrared (the first JWST study of SN 1987A was presented by Larsson et al. 2023). With that in mind, an analysis of the time evolution of SN 1987A in the NIR is now timely for setting the baseline for JWST studies.

At most wavelengths, the evolving remnant of SN 1987A is dominated by the interaction between the SN ejecta and its circumstellar medium. The circumstellar medium has a well-known triple-ring structure, with an inner equatorial ring (ER) and two outer rings in an hourglass-like shape inclined by ~40º (e.g., Tziamtzis et al. 2011). The triple-ring structure is well replicated by a model in which the matter was ejected 20 000 yr ago (Crotts & Heathcote 2000), following a merger of two stars in a binary system, with initial masses of 14–15 M_⊙ and 7–9 M_⊙ (Menon & Heger 2017; Urushibata et al. 2018), creating the known immediate progenitor, a blue supergiant of ~20 M_⊙. The optical and NIR brightness of the ER has been dominated since about the turn of the millennium by hotspots created by the forward shock of the SN crashing into dense clumps in the ER (Borkowski et al. 1997; Sonneborn et al. 1998; Lawrence et al. 2000). The hotspots emit especially strongly in optical and NIR lines (Kjær et al. 2007; Gröningsson et al. 2008). Since about 2009, these hotspots have been fading (Fransson et al. 2015; Arendt et al. 2016, 2020) as the fastest ejecta have cleared the ER and now interact with matter outside the ring (Larsson et al. 2019; Kangas et al. 2022), creating new hotspots there. The fading and destruction of the ER will, however, take some time (e.g., Orlando et al. 2019).

Dust has a large influence on the infrared emission of the system. Cool dust formed in the ejecta (on the order of 0.5 M_⊙ of silicate and carbon dust combined, at ~20 K; Indebetouw et al. 2014; Matsuura et al. 2015; Cigan et al. 2019)¹ has been argued to be responsible for most of the far-infrared (FIR) luminosity of the ejecta, which also provides most of its bolometric luminosity (McCray & Fransson 2016). In the ER, cool dust is less important, but warmer dust dominates the MIR flux (while the inner ejecta contributes very little in the MIR; Arendt et al. 2020) – the spectral energy distribution is consistent with silicate dust at 180 K (Bouchet et al. 2006; Dwek et al. 2010), with a weaker, hotter component at shorter wavelengths. Amorphous carbon dust at 525 K was deemed consistent with the hotter component by Arendt et al. (2016). Later observations revealed emerging emission between 30 and 70 µm at 10 000 days, indicative of a dust mass of >3 × 10⁻⁴ M_⊙ in the ER (Matsuura et al. 2019), compared to ~10⁻⁵ M_⊙ at 8000 days. This suggests ongoing dust formation in the ER or possibly a contribution from the warmest ejecta dust. Compared to the gas mass of the ER (≳6 × 10⁻² M_⊙; Mattila et al. 2010), this would imply a dust-to-gas mass ratio evolving from ≲10⁻³ (affected by the destruction of preexisting dust by SN radiation) to ≲10⁻². The dust responsible for the hotter component, on the other hand, is being destroyed faster than the gas cools (Arendt et al. 2016); either small grains are evaporated by ultraviolet photons, or large grains are destroyed in grain-grain collisions. The properties of the dust created by SN 1987A can shed light on the importance of SNe as dust producers in galaxies (Indebetouw et al. 2014).

Any contribution of even hotter ER dust in the NIR has not been well studied, and other sources also contribute to the NIR emission, such as thermal free-free emission from shocked gas. Much as in the optical (e.g., Fransson et al. 2013), numerous narrow emission lines can be seen in the ER spectrum (e.g., Kjær et al. 2007), dominated by the shocked and unshocked dense clumps that form the hotspots. Broader emission lines originate from collisionally excited ejecta material crossing the reverse shock. Synchrotron emission from material shocked by the ejecta-ER interaction is dominant in the radio (e.g., Indebetouw et al. 2014; Zanardo et al. 2014; Cendes et al. 2018), and while the FIR-to-MIR flux of SN 1987A is dominated by dust, a synchrotron spectrum continuing to the NIR might contribute to the continuum there, and is allowed by reported FIR upper limits (Cigan et al. 2019). These components can be disentangled by studying the evolution of the morphology and spectrum of the NIR emission.

In this paper, we study the morphology and flux evolution of the ER at NIR wavelengths. We use image subtraction to compare the morphology not only between epochs, but also between wavelengths. We especially focus on the NIR continuum emission, which we isolate from and compare to the line emission as well. We describe the data we use in Sect. 2 and our analysis methods and results in Sect. 3. We discuss our findings in Sect. 4 and present our conclusions in Sect. 5. All uncertainties in the paper correspond to 1σ unless otherwise noted.

2 Observations and data reduction

2.1 Imaging observations

We used NIR imaging data of SN 1987A taken using NAOS-CONICA (NACO; Lenzen et al. 2003; Rousset et al. 2003) on the Very Large Telescope (VLT) at the European Southern Observatory (ESO) from 2006 to 2017. These observations are listed in Table 1. The NACO images were taken using AO in the J, H, and Ks bands. The NACO data from 2006 are from Program 078.D-0322(A), (PI Danziger), obtained through the public ESO Science Archive²; the data after 2006 are from Programs 086.D-0713(D), 088.D-0638(D), 090.D-0645(E), 094.D-0505(D), and 096.D-0882(D), (PI Fransson). The individual exposures were taken using a dithering (jitter) pattern of random offsets within a 7-arcsecond box around SN 1987A (8″ for the 2006 images). Due to problems with the NACO detector causing correlated noise in two quadrants, which resulted in strong stripes in the images, a dithering pattern was selected for the 2017 observations to align the SN with the two better quadrants. Star 2 (Walker & Suntzeff 1990; Walborn et al. 1993) was used as a natural guide star for the AO observations. The pixel scale of the NACO images is 0.013″ pixel⁻¹ and the field of view 14″ × 14″. The full width half maximum (FWHM) angular resolution is listed in Table 1 for each observation, but we note that the point spread function (PSF) is not well modeled by a single Gaussian due to extended “wings” common to AO observations.

We also used public pipeline-reduced NIR images from the Near Infrared Camera and Multi-Object Spectrometer (NICMOS) and Wide-Field Camera 3 (WFC3) instruments aboard the HST between 1997 and 2011 (GO 7434, 7821, 9428, 10549, 10867, 11653, and 12241; PI Kirshner). The HST/NICMOS images were taken using the filters F110W, F160W, and F205W and obtained from the Mikulski Archive for Space Telescopes (MAST) Portal³. The WFC3 images were taken using F110W and F160W. The WFC3 images were drizzled (Fruchter & Hook 2002) with a pixel scale of 0.05″ pixel⁻¹ in the final drizzled image, while in NICMOS images the pixel scale is 0.025″ pixel⁻¹ in F110W and F160W and 0.05″ pixel⁻¹ in F205W. In order to compare the NIR morphology to that in the optical, we used images at epochs close to the NACO observations, taken with the Advanced Camera for Surveys (ACS) and WFC3 on the HST. These images were taken using the F625W filter between 2006 and 2017 (GO 10867, 12241, 13181, and 13810, PI Kirshner; and GO 14753, PI Fransson) and have a final pixel scale of 0.025″ pixel⁻¹ in the drizzled images. The HST images used in this study are listed in Table 1 as well.

2.2 NACO image reduction

The NACO images were reduced with a combination of ESO pipelines⁴ and Image Reduction and Analysis Facility (IRAF)⁵ tasks as described below. In some of the epochs, the images exhibited a faint horizontal striping pattern, which was first removed from each raw image with a custom IRAF script. For this purpose, the script creates a one-dimensional image by taking the median of the pixel values along each image row. The 400 brightest pixels in each row were excluded from the median to reduce the influence of bright stars. The one-dimensional image was then subtracted from the columns of the original image. The stripe removal was skipped for images from 2015 as the final image quality was not improved.

After the stripe removal, the exposures were sky-subtracted in sets of three consecutive exposures to remove time-varying artifacts in the exposures. For N exposures numbered 1,…,N, the sets were made out of exposures [1,2,3];[2,3,4]; [N − 2, N − 1, N], yielding N − 2 sets of three exposures. The sky subtraction was carried out one set at a time using the ESO pipeline recipe naco_img_jitter, which subtracts a sky frame (the median of the three object frames) from each object frame, aligns the exposures, and stacks them. The recipe yielded one sky-subtracted image for every set of three exposures, N − 2 in total. The sky-subtracted images were aligned based on the centroid coordinates of a manually selected star in each image. The aligned images were stacked as the median of the images. The sky-subtracted images were visually inspected and images were rejected from the stacking at this step if either individual hotspots or the ring as a whole were not visually discernible. Finally, we ran the stripe removal script on the stacked image to remove any bands or stripes left over by the previous steps. Figure 1 shows a section of the final H-band image for the 2010 epoch as an example, including the ER, the SN ejecta and nearby bright stars – Stars 2 and 3 as in Walker & Suntzeff (1990) and Walborn et al. (1993), and a third star, referred to here as Star A. We also show the sequences of ER evolution in all the NACO and NICMOS images used in this study in Fig. 2. The WFC3 images have a worse image quality and the evolution is not as apparent from them, but for completeness we show them in Fig. 3.

Fig. 1 Section of the reduced NACO image for the 2010 epoch in the H band, showing the hotspot-dominated ER, the inner SN ejecta in the middle, and nearby bright stars.

2.3 Integral field spectroscopy

We also studied the ER using integral field unit data cubes from the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI; Bonnet et al. 2003; Eisenhauer et al. 2003) on the VLT. The SINFONI data we have used here, taken in the J band in 2005 and in the HK bands between 2005 and 2017, have been previously published in several papers (Kjær et al. 2010; Fransson et al. 2016; Larsson et al. 2013, 2016, 2019). We summarize the basic relevant information in Table 2. The data cubes were reduced (sky-subtracted, wavelength- and flux-calibrated) using ESO pipelines and custom software. For details on the observations and the reduction process, see the aforementioned papers. The systematic uncertainty of SINFONI flux calibration is roughly 10%.

Fig. 2 Evolution of the morphology of the ER of SN 1987A in NICMOS and NACO images from 1997 to 2017. In each image, north is up and east is to the left. The field of view is 2.5″. The color scale in each image is logarithmic.

Fig. 3 WFC3 images from 2009 and 2011. In each image, north is up and east is to the left. The field of view is 2.5″. The color scale in each image is logarithmic.

3 Analysis and results

3.1 Photometry of the ER and continuum contribution

3.1.1 Measurement and corrections

The flux from the ER was measured from images as a difference of fluxes in two elliptical apertures. The larger, outer aperture contained both the ER and the SN ejecta in the middle, and the smaller, inner aperture contained the ejecta in order to remove its contribution. The position angles and sizes for both apertures and the eccentricity for the inner aperture were selected to make them match the shapes of the ER and the ejecta as closely as possible⁶. The apertures are shown in Fig. 4. For the NACO images, the outer aperture had an eccentricity of 0.63, a semimajor axis of ~2.2″, and a position angle of 7 degrees. The inner aperture had an eccentricity of 0.7, a semimajor axis of ~0.44″, and a position angle of 77 degrees. The eccentricity of the outer aperture corresponds to an axial ratio of 0.78, which was chosen to match the ratio found by Fransson et al. (2015). The photometry was carried out in Starlink GAIA. For the NICMOS and WFC3 NIR images, the larger, outer aperture was also an ellipse with an eccentricity of 0.63 and position angle of 7 degrees, but a semimajor axis of~1.1″ as the PSF does not suffer from AO effects (bottom panel of Fig. 4). The inner aperture had an eccentricity of 0.6, a position angle of 77 degrees, and a semimajor axis of ~0.4″.

The fluxes measured from the NACO images were calibrated using known magnitudes of Star 2 (see Fig. 1) from Walborn et al. (1993), assumed to be constant, whereas the HST images were flux-calibrated by the pipeline. Errors are not reported in Walborn et al. (1993), and neither are individual NIR measurements of Star 2, but those of the variable Star 3 are. The standard error of the mean⁷ of the Star 3 magnitudes – after subtracting a linear fit to its decline – is between 0.015 and 0.033 mag, depending on band; as Star 2 is of similar brightness, the precision of its magnitude should also be similar. Considering that the decline of Star 3 may not be exactly linear, we adopted 0.05 mag as the uncertainty for Star 2 (i.e., J = 15.09 ± 0.05, H = 15.08 ± 0.05, and K = 15.06 ± 0.05, Vega). Star 2 might, however, conceivably vary in brightness on a longer timescale than that probed by Walborn et al. (1993).

The JHKs magnitudes of the star in the Two-Micron All Sky Survey⁸ (Skrutskie et al. 2006) are consistent within 1σ with the Walborn et al. (1993) magnitude. We also measured the AB magnitudes of Star 2 in the HST filters; these are listed in Table 3. With our assumed uncertainty, the AB magnitudes of Star 2 based on Walborn et al. (1993) are J_ab = 16.00 ± 0.05, H_AB = 16.47 ± 0.05, and Ks_AB = 16.91 ± 0.05; the measured magnitudes in F110W, F160W, and F205W, respectively, stay close to these values from 1997 to 2011, but do show some variation on the order of 0.1 mag over multiple-year timescales (with one outlying epoch in 2006 in F160W), indicating the possibility of similar scatter in NACO JHKs fluxes.

The NACO images have PSFs characteristic to AO, where the narrow, corrected PSF core is surrounded by a fainter, uncorrected halo with broad “wings.” The outer aperture was thus spread far out from the ring to compensate for the PSF. To be able to use a large aperture to measure the ER flux, nearby stars (mainly Star 3 and Star A in Fig. 1) needed to be removed from the images. A PSF fit (using SNOoPY⁹) was made on the three bright stars closest to the ring. The PSF model was then subtracted from the stars in the image, leaving some residuals at the locations of the removed stars. The residual from Star 3 was removed by replacing the values of the residual pixels with the mean value of the nearby background, measured from an annulus around the residual, using the imedit task in IRAF.

In addition to the total fluxes from the ER, we estimated the continuum fluxes. We determined the fraction of continuum emission contributing to the total flux at each epoch of the SINFONI data, which we in turn interpolated to correct the NACO fluxes. We extracted the spectrum at each SINFONI epoch and band using QFitsView¹⁰, with an aperture covering the entire ER slightly farther out than in the NICMOS photometry (excluding the ejecta region in a similar way) and a sky extraction outside the ER. The semimajor axis was approximately 1.3″ and the semiminor axis 1.1″, but as the extraction region is set by eye in QFitsView, there is slight variation in the size of the region, on the order of 10% (affecting only the faint emission in the wings of the SINFONI PSF). In each spectrum, we manually removed each clearly detected emission line, down to the surrounding noise level. We then interpolated the remaining continuum-dominated spectrum using a Gaussian process (GP) regression algorithm (Rasmussen & Williams 2006). We used the Python-based george package (Ambikasaran et al. 2015), which implements various different kernel functions. Matérn kernels with v parameter of 3/2 or 5/2, common in Gaussian process fitting, were used. We show an example of a SINFONI H-band spectrum before and after this process in Fig. 5. The NACO JHKs filters were integrated over the GP-interpolated continuum spectrum and the original, and the ratio between the resulting fluxes is the fraction of continuum flux. The continuum fractions (CFs) at each epoch and band are listed in Table 4.

As is clear from Table 4, the CF does not change much between 2005 and 2017, and the linear interpolation of continuum fluxes at the NACO epochs is straightforward. No extrapolation is needed as the NACO epochs all fall within the 6830–11269 days range of the SINFONI observations. For the J-band observations, we assume the same fraction as estimated for day 6830 at all epochs; this assumption is unlikely to affect the results much considering the nearly constant CF at the HKs bands. No continuum correction is performed for NICMOS F110W and F205W data, as these filters are considerably wider than the NACO ones, extending outside the SINFONI coverage and thus presumably resulting in a different CF that cannot be determined from SINFONI data in the same way. We illustrate this problem in Fig. 6. F160W is closer to H, but no SINFONI data exist for interpolating the CF before 6000 days. While the change of the CF over time is slow or nonexistent over the NACO observations, assuming this for the NICMOS epochs as well is unfounded, as the ER flux evolves faster over these epochs.

We also used the SINFONI spectra to determine the total and continuum flux density at those epochs. We integrated the normalized NACO JHKs filter functions over the GP-ŕìtted SINFONI continuum spectra described above. An assumed 10% uncertainty from the flux calibration and another 10% from extraction and line removal was included in quadrature. The NACO and SINFONI fluxes do not fully agree, as can be seen in Fig. 7. A multiplication factor of ~1.1, ~1.2, or ~0.9 must be included in order for the light curves to match over the full span of the observations, depending on the band. The SINFONI data are affected by similar concerns as NACO, namely large PSF wings, which may have an effect in the extraction of the spectrum, but this is not likely to explain the difference completely, and systematic flux calibration effects may be more important. There are also differences between fluxes measured from NICMOS and NACO images around 7200 days – the NACO flux is lower in J and higher in H. This is affected by the different filters, however, so the numbers especially in J and F110W are not directly comparable, and the discrepancy between NACO and NICMOS is similar to that between NACO and SINFONI only in H.

Fig. 4 Apertures used for measuring the fluxes of the ER and the SN ejecta in the NACO (top) and NICMOS images (bottom). The residuals of two nearby stars from the subtraction of a PSF model are visible in the NACO image. The residual of Star 3 has been manually masked out. North is up and east is to the left.

Fig. 5 Example of a SINFONI H-band spectrum from 2017 before (solid gray line) and after line removal and GP fitting (dotted black line). The extent of the NACO H-band filter is overplotted (dashed red line). Savitzky-Golay smoothing has been applied to the original spectrum.

Fig. 6 Wavelength ranges of NACO and NICMOS filters overlaid on an arbitrarily scaled, Savitzky-Golay-smoothed SINFONI JHK spectrum from 2005 (Kjær et al. 2010).

3.1.2 Photometry results

We applied the multiplication factors of 1.1 (J), 1.2 (H), and 0.85 (Ks) to fluxes measured from SINFONI, as in Fig. 7, to match the light curves as described above. Assuming this correction, we list the ER flux densities in Fig. 5. We also show the light curve and the evolution of the continuum flux density in the NIR compared to the MIR (Arendt et al. 2016, 2020) in Fig. 8. The JHKs light curves peaked between 8000 and 9000 days (i.e., roughly 2010–2012) and declined by about a third in the next 2500 days in H and Ks. The optical light curve in R peaked around day 8000 (2010) and in B slightly later; they also declined by roughly a third in the 3000 days after that. The rise toward the peak is also similar in shape. The evolution of the NIR continuum brightness after 7000 days resembles that in the MIR as well, especially in H and Ks.

The difference between NACO and HST filter functions is apparent from the wildly different light-curve evolution in F110W and J around 2010 (~8500 days) – the FllOW flux reaches values similar to the R band, a factor of 2 higher than in J¹¹. A difference between these filters can be caused by strong emission lines not included in the J filter (i.e., between 8000 and 11000 Å; see Fig. 6); the strongest of these is expected to be He I at 1.083 µm (as the He I 2.058 µm line is the strongest emission line we see in the K-band spectra), only a small fraction of which is included in the NACO J filter. We do note that He I lines measured from our SINFONI spectra (see below) do not increase in strength as drastically as the F110W flux does, but the 1.083 µm line was previously noted to be much stronger than expected based on other He lines, owing to different excitation mechanisms (Mattila et al. 2010). Recent JWST observations (Larsson et al. 2023), with an extremely strong 1.083 µm line from the ER and the RS, corroborate this. Both the soft X-ray (Frank et al. 2016) and Hα emission (Fransson et al. 2015) rose by a similar factor as F110W from ~6000 to ~8000 days. We do not have enough data to be able to say whether the F110W flux kept rising past 9000 days, as did the hard X-rays (Alp et al. 2021). The H and F160W filters are much closer in wavelength coverage, and the light curves in these filters are unsurprisingly much more similar as well.

Fig. 7 Light curves of estimated continuum flux density based on SINFONI and NACO data.

3.2 Emission line fluxes

We measured the fluxes of several narrow emission lines in the SINFONI spectra – lines originating in both shocked and unshocked matter in the ER (which we cannot distinguish between, but the former become dominant over time). As we lack J-band spectra after 2005, we concentrate on lines in the other bands. We list the measured line fluxes in Table 6. In each case, the uncertainty is dominated by the assumed systematic errors we ascribe to the flux calibration and the extraction of the spectrum. We show the evolution of the line fluxes compared to the SINFONI continuum flux density in Fig. 9.

The line fluxes mostly evolved similarly to the continuum flux. There are slight variations in evolution between lines, with the largest difference between the [Fe II] and Br 10 lines in the H band. The flux in all measured lines peaked in the 8708-day spectrum, consistently with the continuum light curve. This is, again, slightly later than in the optical: fluxes of Hα, [O III] and [Fe XIV] lines from the shocked ring reported by Fransson et al. (2015) peaked around 8000 days or, in the case of [O III], around 7000 days. The lines from unshocked matter were disentangled from their high-resolution spectra, but their contribution to the total flux is small and they peaked earlier; thus they do not account for the delay between the NIR and optical peak.

3.3 Expansion velocity of the ER

We used the NACO images to determine the size of the ER, and its expansion velocity, in the Ks band where the hotspots can be resolved most consistently. In 2017, on the other hand, the NACO Ks image is visibly of considerably worse quality than in H with fewer measurable hotspots, and we use H instead¹². We also used F160W on NICMOS, as these images have a better angular resolution than F205W or the WFC3 NIR images (see Table 1). The time baseline from our NACO images corresponds to that in Larsson et al. (2019); they note that at earlier epochs the size evolution is less linear, and thus epochs before 7000 days are not used in our velocity fit. At each epoch, we measured the centroids of all visible hotspots and fitted an ellipse to the centroid coordinates. The semimajor axis of the ellipse from this fit was then converted to a physical size using a distance of 49.6 kpc to the Large Magellanic Cloud (Pietrzyński et al. 2019).

We obtain a best-fit expansion velocity of 690±190 km s⁻¹ since 2006 (~7000 days); we show this fit in Fig. 10. The size evolution appears to have changed around this time in the optical. The measured sizes and the expansion velocity since 2006 are both in excellent agreement with the values measured from optical hotspots by Larsson et al. (2019), who obtained a velocity of 680±50 km s⁻¹. This velocity corresponds to the velocity of shocks propagating through dense clumps in the ER. Before 2006 the optical hotspots may be located slightly further out than in the NIR, but any difference between the two is not significant. In the MIR, on the other hand, Matsuura et al. (2022) measured an expansion velocity of 3920 km s⁻¹, with an initially smaller extent of the MIR ring, but after ~10 000 days larger than in the optical. This velocity tracks a more diffuse emission than the optical/NIR hotspots, however.

Fig. 8 NIR light curve of the ER. Top: total ER flux in the NICMOS and NACO bands, compared to the R band (Larsson et al. 2019). Bottom: continuum flux evolution in the NIR, scaled to match at ~7000 days, compared to the MIR light curve, also dominated by the continuum (Arendt et al. 2016, 2020).

Fig. 9 Emission line flux evolution over time, compared to the continuum flux density in the H and K bands as measured from the SINFONI spectra. All fluxes are scaled to match in the first SINFONI spectrum.

Fig. 10 ER sizes in the NIR and optical as a function of time (points; optical from Larsson et al. 2019), and the best expansion velocity fit to the NIR size (line). Only epochs of >7000 days are included in the fit. After 7000 days, the ER sizes and expansion velocities determined from optical and NIR hotspots are practically identical, and even before that the sizes are consistent within ~1σ.

Fig. 11 Subtractions between selected epochs in the NICMOS F110W, F160W, and F205W bands. Each subtraction is scaled to the ER peak counts in the later image used, and the color scale is constant. North is up and east is to the left. Positive pixels indicate brightening between epochs. The field of view is 2.5″ × 2.1″.

3.4 Difference images

3.4.1 Evolution of the NIR morphology

NIR difference images were constructed to visualize the evolution of the SN ejecta and the ER between different epochs. For all of the images, an earlier image was subtracted from a later one in the same band. All of the NACO images were aligned to the 2006 epoch, and the NICMOS images to the 1997 epoch, with the IRAF tasks geomap and geotran prior to the subtraction. The alignment was done using several manually selected bright stars in both images. A general geometry that takes into account shift, rotation, skew, and distortion between the images was used. The image subtraction was carried out in a modified version of ISIS¹³ using the Optimal Image Subtraction method and manual stamp selection. In this method, the image with the better seeing is convolved with a kernel derived from comparison with the image with the poorer seeing, yielding two images with similar PSFs and background and intensity levels, allowing for a meaningful subtraction of the images (Alard & Lupton 1998; Alard 2000). Image counts are scaled to match the later image in each case. For each subtraction, we used several stars in the field to determine the kernel and the scaling factor, the exact number depending on image quality (i.e., how many stars were visible). We show the resulting subtraction images in Figs. 11 and 12; not all epochs are included here, depending on the time between epochs and image quality.

The NICMOS subtractions mainly show the considerable fading of the central ejecta between 1997 and 2002 and the brightening of the ER. The minimum of the optical light curve of the ejecta in fact occurred between these epochs, and it then brightened until at least 2019 (see, e.g., Larsson et al. 2019); however, the central peak of the ejecta still fades over time (and the edges brighten) as the ejecta expands and becomes less centrally peaked (Fig. 2). The central fading can be seen later as well, including the NACO subtractions, but less clearly, as the relative change in ejecta morphology slows down over time.

Between 1997 and 2002, the brightening of the ER occurs mostly on the north and east sides, and between 2002 and 2006, more strongly on the west side. In the 2010–2006 NACO subtraction, meanwhile, a stronger change is seen in the ER morphology. The western half of the ER brightens considerably between 2006 and 2010, while the eastern half mostly declines. As is apparent especially in the H band, however, there is also brightening indicating expansion at the easternmost edge. The emission extends roughly equally far from the hotspot peaks on both edges. In the following epochs, the J-band image quality is much worse, but in the H and Ks bands a steady decline can be seen all over the ER, except for the ongoing expansion of the ER resulting in a slight brightening toward the advancing edges. A change in the east-west balance as noticeable as that between 2010 and 2006 is not seen afterward, and the ER stays skewed to the west side, which can also be seen in Fig. 2. The subtractions between each epoch and 2006 show a similar timeline: by 2017, the ER is only brighter than in 2006 at the western and, to a lesser extent, eastern edge, indicating expansion. Fransson et al. (2015) noted a similar evolution in the optical.

Fig. 12 Subtractions between selected epochs in the NACO JHKs bands. Each subtraction is scaled to the ER peak counts in the later image used, and the color scale is constant across each row. North is up and east is to the left. Positive pixels indicate brightening between epochs. The field of view is 2. 5″ × 2.1″.

3.4.2 Line and continuum maps

We constructed maps of the [Fe II] 1.644 µm line, the He I 2.058 µm line and Brγ. These were created from the data cubes using QFitsView, subtracting the continuum immediately surrounding the line. In order to study the spatial distribution of the line flux compared to the continuum, we also extracted maps of the continuum redward of each line. The line maps were constructed by including a range of wavelengths large enough to include the different velocities in different parts of the ER (after correcting for the recession velocity of 287 km⁻¹; Gröningsson et al. 2008); in the case of the [Fe II] line we included 1.644–1.647 µm (±350 km⁻¹), in He I 2.0579–2.0640 µm (−350 to 550 km⁻¹), and in Brγ 2.1657–2.1718 µm (−350 to 550 km⁻¹). The continuum map was constructed by, respectively, only including the wavelength ranges 1.652–1.655, 2.0677–2.0751, and 2.1743–2.1804 µm.

In order to study whether the distribution of the continuum and line emission is the same, in each case we normalized the maps by matching the average surface brightness of all visible hotspots across the ER between them, then subtracted the continuum from the line map. While the contribution of the continuum to the line map itself has already been subtracted, this comparison clarifies the difference between the spatial distributions – if the line and continuum originate from the same locations, the subtraction map is practically flat, but any significant positive or negative differences indicate a different spatial distribution. We show the resulting maps in Fig. 13. The continuum flux gradually becomes stronger relative to the lines slightly outside the ER hotspots, as can be seen via the ring of negative flux in the normalized and subtracted images by 2014, especially in the K-band subtractions. Positive hotspots can be seen inside the negative ring in 2014 and 2017, indicating a relatively high line emission contribution there compared to the region slightly outside the ER. Before this, the extent of the lines is closer to that of the continuum flux, resulting in this structure (a negative ring and positive hotspots) being less clear.

Fig. 13 Line minus continuum difference maps of three strong lines at each epoch, normalized using the average hotspot brightness before subtraction. The contours illustrate the extent of the line emission, and positive fluxes indicate that the normalized line emission is stronger than the continuum. The influence of the ejecta and Star 3 is strongest at the weakest ER line of the three, [Fe II] 1.64 µm. Over time, the continuum flux outside the hotspots strengthens compared to the line flux, resulting in a clear negative ring slightly outside the line emission in the later epochs. North is up and east is to the left. The color scale is constant across each row. The field of view is 2.5″ × 2.1″.

Fig. 14 F625W-Ks maps at each epoch, illustrating the stronger diffuse emission in Ks since around 2013. The images were convolved to the same seeing using stars, then normalized to the average hotspot brightness in the F625W image before subtraction in each case. Positive pixels indicate stronger normalized emission in the F625W band than in the Ks band. The location of a star coinciding with the ER is marked with a plus sign. North is up and east is to the left. The field of view is 2.5″ × 2.1″.

Fig. 15 F625W contours overlaid on the Ks maps at each epoch after convolving to the same PSF (see text). The location of a star coinciding with the ER is marked with a plus sign. North is up and east is to the left. The field of view is 2.5″ × 2.1″.

3.4.3 Differences with optical

In order to study the differences between optical and NIR emission, we used the HST F625W maps (mostly Hα line emission) together with the Ks-band maps from NACO. We first performed sub-pixel alignment, as previously described, between the NIR images and the F625W images closest in time to them. We omitted the 2011 NACO epoch here, as the F625W images closest to it in time match the 2010 and 2012 epochs much better and were used in comparisons to them instead (see Table 1). We were unable to use ISIS for subtraction, as the colors of the stars in the image can be expected to vary considerably compared to the ER, and they are thus unsuitable for scaling. We did, however, use ISIS to convolve the image with the better seeing in a similar fashion using field stars. We then used the peak counts of the hotspots in the seeing-matched, aligned images to scale the images to each other, then finally subtracted the Ks image from the F625W image at a similar epoch. The resulting maps are shown in Fig. 14. We also plot the Ks-band images overlaid with contours corresponding to the F625W images in Fig. 15.

A difference is visible in the Ks-band outer emission compared to F625W, especially starting in 2013. A ring of Ks-band emission is located outside the ER in the subtraction, similar to those in the line-continuum subtractions. The difference is less clear in 2011 and especially 2006; the outer ER emission thus gains strength in the Ks band over time, suggesting a faster expansion. The expansion velocities measured from ellipses fitted to hotspots are consistent (Fig. 10), that is to say, the difference is caused by diffuse emission outside the hotspots. The difference is not readily visible from the comparison of the Ks images with the F625W contour maps, however – only the eastern difference can barely be seen without the subtraction. The eastern ER is less dominated by hotspots than the west, making differences caused by other components more visible.

3.5 East versus west fluxes

We repeated our photometry for the east and west halves of the ER separately in order to quantify the morphological evolution seen in the difference images (Fig. 12). We masked out the other half in each NACO image, split at the geometric center of the ER as measured using the flux peaks of the ER at the western and eastern extremes, and repeated the measurement as in Sect. 3.1. The resulting flux densities are listed in Table 7. The east- and west-side light curves and east-west flux density ratios are shown in Fig. 16; east-west ratios are additionally listed in Table 8. It is apparent that the east-west ratio changes considerably between 2006 and 2010, similarly to what can be seen in the difference images, but stays fairly constant after 2010. Both halves of the ER decline similarly after the light curve peak.

4 Discussion

4.1 NIR evolution versus other wavelengths

Most aspects of the NIR evolution of the ER of SN 1987A follow the optical evolution. The ER light curves in the optical and NIR (and MIR when available) are quite similar, peaking close to the same epoch (~8500 days, i.e., ~2010; Fransson et al. 2015; Larsson et al. 2019) in J and slightly later in HKs and declining at a similar rate after that (Fig. 8), indicating the same underlying power source: the ejecta-ER interaction. The evolution of the ER morphology resembles that in the optical as well (Fransson et al. 2015), with the east–west asymmetry appearing around the time of the light curve peak – the east side has been declining in the optical since 7000–7500 days (~2006–2007), and apart from a slight brightening at the eastern edge of the ER due to expansion, our subtractions show the same trend (Fig. 12). Furthermore, the expansion of the ER as a whole also proceeds at the same rate as in the optical when measured from the centroids of hotspots (Fig. 10); in other words, within the hotspots themselves, the locations of optical and NIR emission are practically identical.

The optical east-west ratio in 2006 was ~0.9, but by 2010 it had decreased to ~0.6, and by 2017 below 0.5 (Larsson et al. 2019). In the NIR, the east-west ratio decreases from 0.7–0.8 in 2006 to 0.5–0.6 since 2010, depending on wavelength. In the soft X-rays (Frank et al. 2016), the ratio behaves somewhat similarly to its optical and NIR equivalent: ~1.1 in 2006, ~0.8 in 2010 and ~0.6 in 2015. This east-west evolution is also evident in the MIR (Matsuura et al. 2022), with an increasing fraction of the emission originating in the west over time after roughly 6500 days (2005), albeit going from a ratio of ~3 to ~0.8 instead. This is in stark contrast with the radio (e.g., Zanardo et al. 2018), millimeter (Indebetouw et al. 2014) and hard X-ray emission (Frank et al. 2016), where the east side is the brighter one, especially in terms of polarized radio flux (however, according to Cendes et al. 2018, the radio asymmetry decreased over time and may have reversed after ~12000 days).

However, we see some differences to F625W and MIR emission as well. The diffuse emission outside the hotspots in the Ks band is seen to gradually become stronger (relative to the hotspots themselves) than in F625W, as shown by the F625W–Ks maps (Fig. 14). This only becomes apparent around 2013. As the expansion velocity we measure is based on an elliptical fit to hotspots, this difference has little effect on the velocity measurement. A large fraction of the emission in both bands originates in the hotspots, which results in little visible difference without the subtraction (Fig. 15). The continuum contribution in the Ks band is around 75% in our NACO images, whereas in F625W, the flux is dominated by the narrow and broad components of Hα (Fransson et al. 2013). This contributes to the observed difference if the extended emission in Ks is dominated by continuum, as Fig. 13 indicates.

The expansion velocity of the diffuse ring in the MIR is several times higher than the hotspot expansion velocity in the NIR and optical (~3900 km s⁻¹ as opposed to ~700 km s⁻¹; Matsuura et al. 2022), which tracks shocks propagating through dense ER clumps. This is nearly identical to the velocity in the radio (Zanardo et al. 2013), while the X-rays run the gamut from ~3100 km s⁻¹ in >2 keV X-rays (Frank et al. 2016), through ~1900 km s⁻¹ between 0.5 and 2 keV, to consistent with zero (or even with slow shrinking) in the 0.3–0.8 keV band. The MIR ring was initially smaller than in the NIR before about 10000 days (~2014), when it caught up to the extent of the hotspots, and larger afterward. However, the more diffuse NIR emission component outside the hotspots, responsible for the seemingly larger extent of the ER compared to the optical, may track this faster expansion. Even so, despite the similar light curve, most of the NIR emission thus seems to originate in a different location than the 180 K dust that dominates the MIR (or the X-ray-emitting shocked gas, except for the softest X-rays).

In the MIR, the brightness difference is associated with a shift in temperature: the dust on the east side was hotter in 2005 and on the west side in 2017 (Matsuura et al. 2022), possibly because of a density difference, which could also explain the stronger NIR/optical emission from the interaction on the west side. Larsson et al. (2023) also found that the extent of the reverse shock is largest on the southeast side, the faintest section of the ER with the fewest visible hotspots in our images, indicating a lower density there.

Fig. 16 Comparison between the east and west halves of the ER. Top: evolution of the east-west flux density ratio over time in each band. The clearest change occurs between 2006 and 2010. Bottom: light curves of the east (open symbols and dashed lines) and west (solid symbols and lines) halves. The strongest difference in evolution again occurs between 2006 and 2010, and mainly in the west. Afterward, the evolution is similar in both halves.

4.2 The properties of the NIR continuum

The continuum consistently forms a majority of the total NIR emission from the ER between 6800 and 11 300 days (2005–2017), and possibly earlier as well. We have shown that the spatial extent of the continuum emission is not identical to the line emission. Hotspots in the ER are responsible for much of both the lines and the continuum emission, as is apparent from our images, the lack of evolution in the CF and hence the similarity of the continuum and line emission light curves (see Fig. 9). However, the fraction of the continuum emission outside the hotspots is larger than that of the line emission Fig. 13. This, together with the larger CF in the NIR than in the optical as mentioned above, indicates that the difference between F625W and NIR, with NIR emission similarly being stronger beyond the hotspots, is specifically caused by the continuum emission.

The MIR continuum of the ER is dominated by 180 K silicate dust (Dwek et al. 2010) that is collisionally heated by gas, which itself is heated by the ejecta-ER interaction and subsequently emits X-rays. Emission from such collisionally heated dust should thus follow the X-ray evolution, but the MIR/X-ray flux ratio (and, similarly, the NIR/X-ray ratio) lessens over time, indicating destruction of dust by sputtering (i.e., dust-gas collisions within shocked gas that return the dust grain material into the gas phase). An additional hotter component is also seen with T ≳ 350 K. Arendt et al. (2016) fit this component with 525-K amorphous carbon dust, which is destroyed over time faster than the gas cools – either by ultraviolet photons or by grain–grain collisions, depending on the grain size. These components are not responsible for most of the JHKs emission, however – nor are they spatially coincident with it.

We performed blackbody fits to the continuum flux densities in the MIR (Arendt et al. 2020) and NIR to investigate an additional component that may be present in the NIR – assuming only one dominant continuum source. This is likely not true, as in addition to blackbody emission, thermal bremsstrahlung and nonthermal emission can contribute. We fit two blackbody components, the cooler of which is fixed at 525 K (Arendt et al. 2016). We show these fits and their temperature uncertainties in Fig. 17. We did not include the J-band fluxes, as it is clear from the figure that one blackbody is not sufficient to fit all NIR bands. We used both NACO and SINFONI epochs, but included a correction factor in the latter as described in Sect. 3.1. Finally, we introduced a 10% uncertainty in the MIR fluxes, which otherwise would be weighted much more heavily than the NIR fluxes in the fit. The temperature of the hotter blackbody component in the resulting fits ranges between 1800 and 2800 K. The individual photometric points show some scatter (see, e.g., Fig. 8), and some epochs (e.g., around 9000 d) are affected by a small dip in Ks-band fluxes (possibly due to slight color variation in Star 2), resulting in the highest blackbody fit temperatures. The rest cluster at 1800–2400 K, close to the evaporation temperature of graphite grains (Guhathakurta & Draine 1989).

Small-grain dust that is being destroyed over time by photoe-vaporation could be evident in the NIR data through a blackbody component at a relatively constant (and high) temperature. Grain destruction within the ER could also increase the relative contribution of the outer emission. Despite the expanding diffuse component, most of the NIR continuum is coming from the hotspots (i.e., dense ER clumps) and expanding slower than the MIR continuum emission. Both the warm and hypothetical hot dust would likely have been formed in the same mass loss episode(s) – however, differences in their location could arise after the SN as graphite grains could survive in hotter locations, and the collisional heating of the dust becomes more efficient in the high-density clumps (e.g., Matsuura et al. 2022). This could result in a higher relative prevalence of hot graphite dust in the clumps compared to cooler silicate dust.

Assuming that the NIR is dominated by dust, we can estimate the dust mass similarly to Fox et al. (2010) using their Eq. (2), from which we obtain M_d = F_νD²/(B_vκ), where M_d is the dust mass, F_ν the continuum flux density, B_ν the fitted hot blackbody component and κ a mass absorption coefficient, a function of the grain size and the frequency. With a grain size of 0.1 µm (below ~0.1 µm the size has a minimal effect on κ) we estimate that a hot dust mass on the order of only ~1–4 × 10⁻¹² M_⊙ could account for the entire HKs continuum, depending on the epoch – this range corresponds to the scatter in blackbody temperatures. For comparison, Arendt et al. (2016) estimated the stronger 525 K component to originate from roughly 0.5% of the total dust mass in the ER, which, assuming a total of 10⁻⁵ M_⊙ at 8000 days (Matsuura et al. 2019), amounts to ~5 × 10⁻⁸ M_⊙ – we obtain a value of ~2 × 10⁻⁸ M_⊙ from the 525 K graphite component in our blackbody fit around the same epoch. Matsuura et al. (2022) obtained a similar estimate of ~10⁻⁸ M_⊙. A low mass of graphite dust is not surprising based on the low carbon abundance in the CNO-processed ER (Lundqvist & Fransson 1996). Since the 525 K component is both stronger and cooler than the NIR continuum (therefore having a much larger emitting area) and has a lower κ (Fox et al. 2010), the ~2000 K dust mass is orders of magnitude lower still.

However, the J-band continuum flux is clearly affected by another component as well; this component can also extend into the H band. The component is already visible in 2006, when the F110W/J ratio is lower than at later times, making it unlikely to be simply the He I 1.083 µm line at the edge of the J filter. The temperature fits above assume that there is only one continuum component in H and Ks and no contribution from, for example, blended weak lines, thermal bremsstrahlung from gas either at the ER or the reverse shock, or a broad-line component from the latter. The blackbody fit is affected in reality by these additional components. The H band is more densely filled with lines than the Ks band, which makes it more difficult to disentangle broad reverse-shock lines from the continuum; thus the effect of reverse-shock lines would be expected to be higher in H, which can increase the apparent blackbody temperature. It is likely that the continuum includes a reverse-shock + weak-line component as well. The reverse shock flux in optical lines is a non-negligible fraction of the total line flux (Fransson et al. 2013, 2015). Thermal bremsstrahlung from gas at ≳ 10 000 K, typical for the dense ER clumps (e.g., Orlando et al. 2019), on the other hand, could have a flat spectrum across the HKs bands, which would also distort the temperature from a blackbody fit by increasing the relative contribution in the weaker band (and thus usually the temperature, as this band tends to be H). This component may contribute significantly to the total NIR flux, albeit less so in J. As a result, the dust mass we obtain is an upper limit only.

Synchrotron emission from the forward or reverse shock may contribute as well. It plays an important role in the radio and, to a lesser extent, FIR (e.g., Cendes et al. 2018; Cigan et al. 2019). Its extension into NIR (which would be strongest in Ks) is a candidate for a relatively fast-expanding diffuse continuum component. Extrapolating the synchrotron spectrum from ALMA frequencies to the NIR with a ν^−0.7 power law (Cigan et al. 2019), one obtains a flux density of ~0.1 mJy (strongest in the Ks band: ~0.13 mJy) at 9280 d, around the peak of the NIR light curve. This is roughly a third of the continuum flux in the NIR at that epoch. The radio fluxes have increased until at least 11000 days (2017), implying that the synchrotron component at the post-peak epochs can constitute even more than a third of the NIR continuum flux. The effect of a synchrotron contribution would be largest in Ks, decreasing the apparent blackbody temperature (i.e., possibly offsetting the effects of bremsstrahlung and the weak/broad lines on it). This, however, assumes no break frequency in the spectrum between 1 and ~800 µm. The ER is not visible in the FIR images of Indebetouw et al. (2014), but the upper limits shown in Cigan et al. (2019) do not constrain whether such a break is present.

Lundqvist et al. (2020) used an observed break frequency and the age of the SN remnant SNR 0540-69.3 to estimate its magnetic field; we can similarly estimate that a break frequency ν_c requires a magnetic field B = (τ/6 × 10¹¹ s)^−2/3ν_c^−1/3 (Pacholczyk 1970), where τ is the lifetime of the synchrotron-emitting electrons. With τ ≈ 15 yr (at 9280 d, since about the beginning of the ER-ejecta interaction) and ν_c > 3 × 10¹⁴ Hz, corresponding to 1 µm, we obtain B ≲ 1.8 mG. Berezhko et al. (2011) argued that ≲10 mG was required to produce the observed radio and X-ray spectrum at 20 yr. This would translate into ν_c ~ 6 × 10¹¹ Hz (~500 µm, roughly the lowest value allowed by Cigan et al. 2019) and an insignificant synchrotron component in the total NIR emission. However, Petruk et al. (2023), based on the Orlando et al. (2019) models, place the magnetic field strength much lower than 1 mG in the hotspots and beyond, allowing the synchrotron spectrum to continue unbroken through the NIR. Even in the Berezhko et al. (2011) scenario, regions of low magnetic field may still contribute, and we note that the continuum “ring” outside the hotspots clearly does not dominate the total continuum emission either.

Fig. 17 Two-component blackbody fits to the MIR and NIR flux densities, with the cooler dust fixed at a temperature of 525 K (Arendt et al. 2016). The J-band points clearly do not fit the same blackbody function as the HKs points in most cases, and are not included in the fits. The SINFONI epochs (with the correction factor included) are plotted with triangles and the NACO epochs with circles. For clarity, all epochs except 6832 days are multiplied by arbitrary constants.

5 Conclusions

We examined the morphology of the ER of the remnant of SN 1987A in the NIR, its evolution over time, and differences between the NIR and other wavelengths. We also attempted to isolate the continuum and line emission from the ER, and examined differences in morphology between them. Our conclusions are as follows:

• The evolution of the NIR morphology broadly follows that seen in the optical, with a similar light curve, expansion speed as measured from hotspots (~700 km s⁻¹), and a decreasing east-west brightness ratio. However, the NIR emission at later epochs is stronger outside the ring defined by the hotspots than in F625W (but still weaker than the hotspots), hinting at a faster-expanding diffuse component;

• Although the line and continuum emission follow similar light curves, the continuum emission is similarly stronger than the line emission outside the hotspots in the later epochs. This indicates that the continuum is responsible for the outer emission not seen in the optical;

• Continuum emission constitutes a majority of the total flux in the NIR in all bands and at all epochs, but not all bands can be simultaneously fitted using a blackbody function. If the HKs-band continuum is dominated by dust, its temperature is in the vicinity of 2000 K from 2005 to 2017. The location of this hypothetical hot dust would also be different from the warm dust dominating in the MIR, which expands much faster and is more diffuse in morphology. This dust does not dominate in the J band, however, and thermal bremsstrahlung is also expected to contribute; thus, the mass of any hot dust is not more than a few times 10⁻¹² M_⊙ – orders of magnitude below the warm dust mass;

• Synchrotron emission can contribute significantly to the NIR continuum in regions where the magnetic field is below ~2 mG, especially in the Ks band. It cannot dominate the full NIR continuum emission in the ER, which mostly originates in the hotspots, but can be the source of the outer continuum emission.

With the recent launch of the JWST with its superior spatial resolution and sensitivity, the community is now able to study the ER and other regions in the remnant of SN 1987A in greater detail. These observations will no doubt cast more light on the different emission components in the NIR and allow for better constraints to be placed on their origins.

Acknowledgements

We thank the anonymous referee for their suggestions. We thank John Danziger, Jason Spyromilio, Jesper Sollerman, Anders Jerkstrand, Roger Chevalier, and Rubina Kotak for contributing to the NACO observations and for helpful discussion. This work was supported by the Swedish National Space Agency and the Swedish Research Council. S.M. was funded by the Academy of Finland project 350458. This work is partially based on observations made with the NASA/ESA Hubble Space Telescope (programmes GO 7434, 7821, 9248, 10549, 10867, 11653, 12241, 13181, 13810, and 14753), obtained through the data archive at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555.

References

All Tables

All Figures

	Fig. 1 Section of the reduced NACO image for the 2010 epoch in the H band, showing the hotspot-dominated ER, the inner SN ejecta in the middle, and nearby bright stars.
In the text

	Fig. 2 Evolution of the morphology of the ER of SN 1987A in NICMOS and NACO images from 1997 to 2017. In each image, north is up and east is to the left. The field of view is 2.5″. The color scale in each image is logarithmic.
In the text

	Fig. 3 WFC3 images from 2009 and 2011. In each image, north is up and east is to the left. The field of view is 2.5″. The color scale in each image is logarithmic.
In the text

	Fig. 4 Apertures used for measuring the fluxes of the ER and the SN ejecta in the NACO (top) and NICMOS images (bottom). The residuals of two nearby stars from the subtraction of a PSF model are visible in the NACO image. The residual of Star 3 has been manually masked out. North is up and east is to the left.
In the text

	Fig. 5 Example of a SINFONI H-band spectrum from 2017 before (solid gray line) and after line removal and GP fitting (dotted black line). The extent of the NACO H-band filter is overplotted (dashed red line). Savitzky-Golay smoothing has been applied to the original spectrum.
In the text

	Fig. 6 Wavelength ranges of NACO and NICMOS filters overlaid on an arbitrarily scaled, Savitzky-Golay-smoothed SINFONI JHK spectrum from 2005 (Kjær et al. 2010).
In the text

	Fig. 7 Light curves of estimated continuum flux density based on SINFONI and NACO data.
In the text

	Fig. 8 NIR light curve of the ER. Top: total ER flux in the NICMOS and NACO bands, compared to the R band (Larsson et al. 2019). Bottom: continuum flux evolution in the NIR, scaled to match at ~7000 days, compared to the MIR light curve, also dominated by the continuum (Arendt et al. 2016, 2020).
In the text

	Fig. 9 Emission line flux evolution over time, compared to the continuum flux density in the H and K bands as measured from the SINFONI spectra. All fluxes are scaled to match in the first SINFONI spectrum.
In the text

	Fig. 10 ER sizes in the NIR and optical as a function of time (points; optical from Larsson et al. 2019), and the best expansion velocity fit to the NIR size (line). Only epochs of >7000 days are included in the fit. After 7000 days, the ER sizes and expansion velocities determined from optical and NIR hotspots are practically identical, and even before that the sizes are consistent within ~1σ.
In the text

	Fig. 11 Subtractions between selected epochs in the NICMOS F110W, F160W, and F205W bands. Each subtraction is scaled to the ER peak counts in the later image used, and the color scale is constant. North is up and east is to the left. Positive pixels indicate brightening between epochs. The field of view is 2.5″ × 2.1″.
In the text

	Fig. 12 Subtractions between selected epochs in the NACO JHKs bands. Each subtraction is scaled to the ER peak counts in the later image used, and the color scale is constant across each row. North is up and east is to the left. Positive pixels indicate brightening between epochs. The field of view is 2. 5″ × 2.1″.
In the text

Fig. 13 Line minus continuum difference maps of three strong lines at each epoch, normalized using the average hotspot brightness before subtraction. The contours illustrate the extent of the line emission, and positive fluxes indicate that the normalized line emission is stronger than the continuum. The influence of the ejecta and Star 3 is strongest at the weakest ER line of the three, [Fe II] 1.64 µm. Over time, the continuum flux outside the hotspots strengthens compared to the line flux, resulting in a clear negative ring slightly outside the line emission in the later epochs. North is up and east is to the left. The color scale is constant across each row. The field of view is 2.5″ × 2.1″.

In the text

Fig. 14 F625W-Ks maps at each epoch, illustrating the stronger diffuse emission in Ks since around 2013. The images were convolved to the same seeing using stars, then normalized to the average hotspot brightness in the F625W image before subtraction in each case. Positive pixels indicate stronger normalized emission in the F625W band than in the Ks band. The location of a star coinciding with the ER is marked with a plus sign. North is up and east is to the left. The field of view is 2.5″ × 2.1″.

In the text

	Fig. 15 F625W contours overlaid on the Ks maps at each epoch after convolving to the same PSF (see text). The location of a star coinciding with the ER is marked with a plus sign. North is up and east is to the left. The field of view is 2.5″ × 2.1″.
In the text

Fig. 16 Comparison between the east and west halves of the ER. Top: evolution of the east-west flux density ratio over time in each band. The clearest change occurs between 2006 and 2010. Bottom: light curves of the east (open symbols and dashed lines) and west (solid symbols and lines) halves. The strongest difference in evolution again occurs between 2006 and 2010, and mainly in the west. Afterward, the evolution is similar in both halves.

In the text

Fig. 17 Two-component blackbody fits to the MIR and NIR flux densities, with the cooler dust fixed at a temperature of 525 K (Arendt et al. 2016). The J-band points clearly do not fit the same blackbody function as the HKs points in most cases, and are not included in the fits. The SINFONI epochs (with the correction factor included) are plotted with triangles and the NACO epochs with circles. For clarity, all epochs except 6832 days are multiplied by arbitrary constants.

In the text