A new study on a type Iax stellar remnant and its probable association with SN 1181

High signal-to-noise spectroscopic data of both star and nebula were obtained with the OSIRIS instrument on the 10-m Gran Telescopio Canarias (GTC) telescope in late 2016 as part of our planetary nebula (PN) follow-up program. The faint nebula spectra from SN 1181 was presented in Paper I with only the [S ii] 6716 and 6731 Å doublet and [Ar iii] 7136 Å emission lines detectable. Here, we focus on the spectra we obtained for the stellar core.

OSIRIS covered the wavelength region from 3700 to 7000 Å. The slit width was 0.8″ (Paper I). The most prominent stellar feature is the O vi 3811/34 emission line of remarkable strength below 4000 Å (Figure 1), with additional broad spectral features beyond 4200 Å. The narrow absorption lines are known diffuse interstellar bands (DIBs) at 4428, 4502, 4727, 4762, 5418, 5456, 5779, 5796, 6204, 6281, and 6613 Å, identified using the NASA catalogue³³3https://leonid.arc.nasa.gov/DIBcatalog.html (Jenniskens & Desert, 1994).

We analyzed previously un-extracted spectra of the stellar core from the SparsePak fibre-based (fiber diameter 4.7″), pseudo-integral field unit on the 3.5-m WIYN telescope taken in late 2014. The spectral observations of the surrounding nebula Pa 30 are presented in Paper I.

Low resolution ($R\sim 540$) optical spectroscopy of Pa 30’s central star by our amateur collaborators with a 20-cm $F/D\sim 5$ telescope at Kermerrien observatory in Porspoder, France, brought its remarkable spectrum and its strong O vi 3811/34 emission line to our attention in October 2018. The low signal-to-noise spectra were not properly flux-calibrated and were scaled and smoothed to be compared against our OSIRIS/GTC spectrum in Figure 1.

Paper I adopted $A_{V}=2.4\,\pm\,0.1$ mag but field stars now imply this may be an underestimate. Figure 2 shows parallax and $A_{0}$ from Gaia DR3 field stars within a 5 arcmin radius (black points) from the central star of Pa 30. Red points are stars within 2 arcmin. The curve is the model extinction plot from Vergely et al. (2022) from the explore platform⁶⁶6https://explore-platform.eu/sdas at a spatial resolution of 25 pc (30 arcmin at the distance of Pa 30).

At the parallax of the central star (dashed line), the nominal extinction is $A_{0}=2.9$ mag, albeit with a large scatter between individual stars. The extinction is primarily due to clouds at around 400 pc and at 800 pc. The Vergely et al. (2022) curve ends at a distance of 3 kpc. Gaia DR3 finds some higher extinctions beyond 2.5 kpc but this is not included in the curve.

The data indicates an extinction $A_{0}=2.8\pm 0.4$ mag. The central value is the same as listed for our star in Gaia DR3, but that is fortuitous in view of the complexity of the spectrum. In this paper we will argue that the extinction for the central star is at the higher end of the range.

The central star of Pa 30 has a unique spectrum (Figure 1) with very broad, very high excitation emission lines and a striking discontinuity bluewards of 4000 Å. There is a complete absence of hydrogen, helium and nitrogen lines. The broad lines are mostly identified with very high ionization stages of oxygen (Figure 3). The sharp and significant emission line that cuts-on around 4000 Å is identified as O vi 3811/34 Å. This feature turns over near the atmospheric cut-off as confirmed in the recent spectrum of Garnavich et al. (2020). This doublet is the defining criteria for the WO spectral type in both Pop I stars (e.g., Crowther et al., 1998; Gvaramadze et al., 2019), as well as old low-mass central stars of planetary nebulae (designated as [WO]) before they enter the white dwarf phase (e.g., Smith & Aller, 1969; Sanduleak, 1971). The broadened emission at 5270/91 Å is identified as O vi. The central star of Pa 30 seems unique in arising from neither pathway if the double degenerate merger and accretion fed Type Iax SN are correct (Gvaramadze et al., 2019; Oskinova et al., 2020).

In addition to these strong, broad features, a number of weaker, individual lines are seen. Their profiles can be accurately fitted with Gaussians with FWHM of $0.015c$. The wings extend to velocities as high as 15,000 km s^-1 ($\approx 0.05c$). Several lines coincide with O viii and are tentatively identified as such. They are well centered on the systemic velocity.

The broad, prominent spectral features are indicative of a strong stellar wind. We modeled it using the non-LTE line-blanketed radiative transfer code cmfgen⁷⁷7http://kookaburra.phyast.pitt.edu/hillier/web/CMFGEN.htm adapted to WC and WO stars (Hillier & Miller, 1998; Hillier, 2012; Hillier & Dessart, 2012). Input parameters were our OSIRIS/GTC spectrum that includes the O vi doublet at 3811/44 Å with the interstellar extinction derived in Paper I, and a wind velocity of 15,000 km s^-1. A classic velocity law is used,

We explored a range of stellar luminosity and temperature, mass-loss rate, stellar radius, velocity law and chemical abundances, summarized in Table 1. Given the limited spectral coverage of our OSIRIS/GTC spectrum (and indeed of all available optical spectra- except that from Gaia, see above), we could not find a unique model that completely reproduces the observations (unlike Gvaramadze et al., 2019, although they also stated large uncertainties, e.g., Table 1). Nevertheless, many models could reproduce most spectral features. The results shown are best approximations and/or upper limits to the individual parameters. Most models tested could reproduce the visual and near-infrared photometry. Determining errors on the fits is tricky and we refrain from giving error bars as the errors are not easily transformed into a standard deviation and strong coupling exists between some parameters. Potentially four species can have mass fractions $>$0.1, so changing the mass fraction of one species necessitates a corresponding change of at least one other species. Therefore, below, we provide viable ranges for the stellar parameters and indicate the properties used to constrain them.

The three primary modeling constraints come from matching the mean flux level, the strength of the O vi 3811/34 Å  doublet, and the O v 5572–5604 Å line strength. Only models with luminosities between 25,000 and 60,000 L_⊙ provide reasonable fits. Below the low luminosity limit either the O vi doublet is too weak and/or O v too strong, while at the upper limit either O vi is too strong and/or O v is too weak. At very high luminosities the O vi doublet strength is suppressed (but other O vi lines remain strong) since the upper levels preferentially decay directly to the ground state. We hence set an upper limit at $\leq$260,000 K⁸⁸8This is the classical effective temperature defined at a Rosseland optical depth of 2/3. Due to the dense wind this is somewhat lower (around 10%) than would be defined using the radius of the star at the sonic point..

The peak near 5300 Å is due to O vi offset to the red from observations, perhaps due to the overlapping O v line being too strong in the model (Figure 3). The model matches bumps seen at $\sim$4500 Å and $\sim$4659 Å, but fails to predict the feature near 4250 Å which is assigned to O viii. An alternative suggestion for this emission is Ne viii, which has approximately the same wavelength as the O viii line. In a model with a slightly lower mass loss, the Ne viii is clearly visible in predicted spectra though the overall fit to the spectrum is worse. As apparent from Figure 3, the broad emission lines blend together making it impossible to see the true continuum.

We provide some upper limits for mass fractions of key elements (He, C, O and Ne) assuming Solar abundance ratios for other elements, that is Na, Mg, Al, Si, S, Ar, Cl, Ca, Fe, and Ni⁹⁹9f For the following transitions: Na iv - Na x, Mg iv - Mg xi, Al iv - Al xii, Si iv - Si xi, S v - S xi, Ar vi - Ar x, Cl vi - Cl viii, Ca vi - Ca xii, Fe vi - Fe xvii, and Ni vi - Ni xvii. In these hot models helium has little influence on the spectra. From the apparent absence of the He ii 4686 Å line we place an upper limit of 0.4 on $X_{\rm He}$. The O vi features at 3811/34 Å and 5270/91 Å are clear indicators that the stellar atmosphere is oxygen-rich. Therefore we used an oxygen-rich model atmosphere with an oxygen mass fraction $<0.647$. This reproduces the O vi doublets, as well as O v lines. All the strongest features can be reproduced but not some weaker individual ones. The lines at 4500 and 4658 Å are sitting on a broader O vi feature whose strength depends on continuum choice, uncertain due to the density.

Carbon-free and carbon-enhanced models were tested. The C iv line can contribute to the 4568Å line but otherwise there is no easily distinguishable carbon line contribution. We set an upper limit for $X_{\rm C}\leq 0.261$. We explored a range of neon abundances (Table 1) and find an enhanced neon abundance could contribute to the spectral features at the 3600–4000 Å plateau (Ne vii), at 4341Å and at 6070Å (Ne viii). Inclusion of neon does not affect our ability to get a good fit to the oxygen lines. On the other hand, the 4341 Å line can also be identified with O viii (see below).

The fit shown in Figure 3 provides a ‘maximum’ model where as many features as possible are reproduced. Oxygen is well fitted, but carbon may be less certain and neon is likely significantly over predicted. The plotted model has $X_{\rm Ne}=0.436$ and may be considered as the upper limit. A more carbon-rich composition cannot reproduce the spectrum. The ranges of chemical abundance of carbon, oxygen and neon in Table 1 agree with those derived by Oskinova et al. (2020) from the central star’s X-ray spectrum.

The broad wind profile ($\approx$15,000 km s^-1) provides a slightly better fit to the broad spectral profiles for both the O vi doublets and O v feature at 5900 Å (Figure 3). This value is close to the measured value of Gvaramadze et al. (2019, 16,000 km s^-1). Since O v is present, this emerges from the outer regions, and so a higher terminal velocity is needed to fit both oxygen features simultaneously. As found by Gvaramadze et al. (2019) radiation pressure appears insufficient to drive the wind. Though the OSIRIS/GTC spectra do not extend below 3600 Å, the cmfgen models predict an O vi line at 3435 Å similar to Gvaramadze et al. (2019) and has been observed in the spectra of Garnavich et al. (2020).

After finalizing the model, we compared it to HST ultraviolet (UV) spectra of the stellar remnant (Sect. 2.4.4). To fit the UV data we needed to adjust the reddening to E(B-V)$=1.05$. Although we eventually updated our reddening, our initial models started off assuming the reddening from Paper I. With the new reddening the fluxes no longer matched, so we recomputed the model with scaled parameters that would yield the same spectral shape, and this model is shown in Figure 3.¹⁰¹⁰10To a good approximation an accurate distance is not needed for the spectral fitting. For stars with dense stellar winds, models with the same abundances, effective temperature, and value of $\dot{M}/R_{*}^{3/2}$ produce similar spectra – only differing in flux (Schmutz et al., 1989). The parameters for one distance can be scaled to a new distance using the following scaling relations: $L\propto d^{2}$, $R_{*}\propto d$ and $\dot{M}\propto d^{3/2}$. Given the O vi 3811/34 emission line’s strength, the estimated reddening should be more accurate than that obtained using optical fluxes over a limited wavelength range.

No further model modifications were made and there is excellent agreement between the STIS/HST fluxes and the OSIRIS/GTC spectrum. Figure 3 shows the cmfgen model correctly predicts the O v and Ne viii blend at $\sim$2800 Å and the strong O v emission at $\sim$1370 Å. The model also predicts strong O vi resonance emission (1032, 1038 Å). Absorption shortward of 1500 Å is from the C vi resonance transition (1548, 1551 Å) being too strong, possibly indicating a lower carbon abundance than in the model. The strong 2200 Å band and its Galactic variation make it difficult to judge agreement between model and observation for features in the band.

We computed additional models to better constrain the parameters and potentially improve the fit, but this did not lead to a change in the adopted parameters. One possible exception is the neon abundance – a model with the abundance reduced by a factor of 2 might give a better fit to the UV spectrum. However, a better understanding of the interstellar 2200 Å band is required to draw firm conclusions.

In a recent study of two WO stars in the Large Magellanic Cloud by Aadland et al. (2022) it was found that X(C)/X(O) was greater than 5 despite the presence of strong O vi features at 3811/34Å. However this star is descended from a massive O star and has a very different evolutionary history than the central star of Pa30. Further on, the spectra are quite distinct - in the WO stars, for example, the C iv spectrum is still very prevalent and neon lines are relatively weak and can be matched by a neon abundance similar to that expected from the processing of CNO into ¹⁴N in the CNO cycle, and the subsequent processing into ²²Ne. The high temperatures used in our models will weaken the C iv lines, but this alone cannot explain the weakness of the carbon spectrum and the strength of the oxygen spectrum.

The X-ray spectrum indicates an additional high temperature gas component (Oskinova et al., 2020) that is strong compared to what is typical for OB stars. The unabsorbed X-ray to bolometric luminosity ratio, $\log(L_{\rm X}/L_{\rm bol})$, is in the range of $-4.8$ to $-6.1$, one to two orders of magnitude larger than OB stars, i.e., $\log(L_{\rm X}/L_{\rm bol})=-6.912\,\pm\,0.153$ (Sana et al., 2006). Oskinova et al. (2020) attribute the bulk of the remnant’s thermal X-ray emission to the wind’s outer regions.

The wind model describes well the shape and many features of the optical spectrum. However, some weaker lines are not fitted and the 4340 Å feature requires Ne viii at a high neon abundance. This line coincides with a high-excitation O viii line. We therefore explore a spectral contribution from the high-temperature gas component which is not included in the wind model.

As a hydrogenic ion, O viii has a restricted number of lines. The 4340 Å line coincides with the O viii 8–9 transition and the 6069 Å line with the O viii 9–10 transition. These transitions have a wavelength ratio of 5/7 and the two line profiles overplot accurately when wavelengths are scaled by this ratio. The observed intensity ratio of the two lines (after subtraction of an estimated continuum) is reproduced by these O viii lines for an extinction in agreement with the foreground extinction. This strengthens identification as O viii, although other lines may contribute to these features. The temperature required for O viii is $>$1 MK, far above the photospheric or wind temperature.

We ran exploratory photo-ionisation models using the cloudy code (ver. 17.03 Ferland et al., 2017). In the ‘coronal mode’ gas temperature is pre-defined and all lines are assumed collisionally excited. There is no input radiation field. A single shell model was used, with a constant density and temperature and without considering optical depth effects. Models were run for a range of temperatures and abundances. The hydrogen density was set to $\log n=7.5$. All lines were assumed to have a Gaussian profile with sigma of $v/c=\Delta\lambda/\lambda=0.0092$ as found from line fitting. The cloudy line output list was convolved with this Gaussian profile to obtain model spectra. The model spectra were multiplied by the extinction curve of Cardelli et al. (1988) and scaled to the peak intensity of the 4340 Å line. The ratio of the integrated line strength of the 4337 Å and 6064 Å lines was used to derive the extinction $A_{V}$, where all contributions to these features were taken into account.

The broad wind features at 3800 Å and 5200 Å are not reproduced in this model: the temperatures are too high for significant O vi emission. However, the weaker lines absent from the wind model could be reproduced with this hot gas.

The model fit is shown in Figure 4 along with the observed spectra for a gas temperature of $T=4$ MK. The best fit gave abundances (number ratios) of He/C/N/O/Ne = 0.73/0.07/0.02/0.16/0.01, where helium must be considered as an upper limit. These abundances are within the wind model range. The brightest predicted lines between 3000 Å and 7000 Å are listed in Table 2.

An estimated constant-value continuum was subtracted from the observed spectrum across the full wavelength range with no attempt to correct for any continuum slope. The line strengths depend on the continuum choice. This introduces an uncertainty, as the dominant broad wind features obliterate most of the continuum. The subtraction of a constant level from the WIYN spectrum brought both the blue and red ends to zero. For the GTC spectrum, a constant subtraction left an offset between the ends. We used the WIYN spectrum for the fit. The procedure indicated an extinction of $A_{V}=2.9\pm 0.2$ mag from the WIYN data, in excellent agreement with that found in Sect. 2.5 and with the UV extinction indicated by the wind model. We adopted $A_{V}=3.0$ mag for the final model. The model is shown with the continuum-subtracted WIYN spectrum in Figure 4.

The cloudy model fits the 4337, 4655, and 6064 Å  with the O viii 8-9, 10-12, and 9-10 transitions, respectively, though the observed 4655 Å line is stronger than predicted. Another O viii line at 5790 Å coincides with a peak inside the broad wind band. This identification is not as secure because of the unknown optical depth in this wind band. The 4800 Å feature is reproduced with Ne ix. The He ii line at 4686 Å contributes to the observed line at this wavelength together with the O viii line. A good fit requires a helium abundance ratio of He/O = 4 (by number). However, there is also a C vi and an Ne x line at this wavelength which could contribute. Hence, this helium abundance is considered as an upper limit. The 4940 Å feature is not reproduced. A weaker model line predicted at 5290 Å is due to C vi but it is likely obscured by the broad wind feature.

Our cloudy hot-gas model is able to reproduce most lines not present in the wind model, using abundances equivalent to that of the wind model, for a gas temperature of $T=3.5$–$4$ MK. Importantly, all lines which are predicted to be detectable are seen in the spectrum, with a sole exception the Ne x line at 6900 Å, which is affected by a telluric feature.

The carbon abundance is poorly constrained because of a lack of bright lines though several fainter carbon lines do coincide with features in the spectra. Figure 5 shows a model with the carbon abundance (by number) increased to C/O = 4.5, and with a lower helium abundance. This produces a good fit to many of the lines but adds a strong component to the 5290 Å line which is already well fitted in the wind model without this component. There is no other evidence for such a high carbon abundance and the wind model needs oxygen-rich gas. We therefore do not adopt these abundances.

Figure 6 shows the same fit, but now the observed spectra are shown with the wind model subtracted. Only wavelength regions free of the broad wind features are shown. The wind model accounts for almost all of the continuum present at these wavelengths. A small excess was subtracted shortward of 5000 Å. This shows that the combination of the wind and hot gas models can explain much of the observed spectra.

The abundances in the cloudy model were adopted from the wind model. The nitrogen abundance is not constrained in this model because of a lack of bright lines. There are potential C vi lines at 6195 Å and at 4683 Å, which can fit lines at these locations, however these lines need a C/O ratio well above unity to fit the corresponding features, while the wind model does not support C/O $>1$. In addition, the 5290 Å C vi line becomes too strong at higher carbon abundance. The temperature is mainly constrained by O viii lines which affect the 6064 Å line at temperatures below 3 MK, and the Ne x line at 6894 Å which becomes too strong for $T>4$ MK. Because of the model’s simplicity these constraints should be viewed with caution, along with the assumption of uniform temperature and density, and the fact that other wind lines may be present.

The line widths and temperature indicate that the hot gas is closely related to the stellar wind. The line widths are similar or slightly smaller to that of the wind. The assumption of a gaussian profile is not identical to the wind profile of Eq (1). The symmetric line profiles favor a physical location embedded in the outer wind. The gas can be heated by shocks as regions with different speeds collide. This is a common situation in OB stellar winds. A smaller line width for the hot gas could argue for a location further inside the wind but this remains to be investigated.

Figure 7 shows the Gaia BP/RP spectrum. It covers a much larger wavelength range to the red. The sensitivity is less than that of our other spectra, and there are calibration issues as mentioned before. The top panel shows the wind model overplotted with the Gaia spectra. It shows that the wind feature at 7500–8000Å is well reproduced. We use this ‘above atmosphere’ spectrum to discuss the neon abundance.

Table 2 lists a number of expected neon lines. The chosen neon abundance was based on lines in our optical spectra. The predicted Ne x line at 6894 Å is too strong but this coincides with a strong telluric feature that was removed from the ground-based spectra.

The bottom panel shows the Gaia spectrum with the wind model (red, dashed; shifted down for clarity). The lower red line shows the hot-gas model. The blue-dashed line shows the same model but with three times higher neon abundance. The neon lines now disagree with the observed spectrum. From this we deduce an upper limit of Ne/O $<0.15$ by number.

The wind spectrum is seen to overpredict the line at 8200 Å. In the wind model this is a Ne vii complex, while the hot-gas model has an O viii line at this location. The model that is shown has a high neon abundance. Again this points at a lower neon abundance for the remnant.

The presence of strong magnetic fields has been proposed by Gvaramadze et al. (2019). This is based on predictions of WD merger models which yield $B\sim 200$ MG (Ji et al., 2013), and on the expectation that the stellar wind is magnetically driven from a very rapidly rotating star. For the SN 1181 stellar remnant, Kashiyama et al. (2019) carried out magneto-hydrodynamic simulations that argue for a WD remnant with a strong magnetic field (20 – 50 MG) but still an order of magnitude lower than Gvaramadze et al. (2019).

Magnetic fields have been proposed for hot ($T>60$ kK), hydrogen-poor (DO) white dwarfs with O viii lines and other ultra-high excitation lines (UHE) in absorption. Reindl et al. (2019) attributes these coronal lines to gas trapped in an equatorial, magnetically-confined, shock-heated magnetosphere. The confinement requires

where $B$ is the magnetic field strength, $R_{*}$ the stellar radius, $v_{w}$ the wind speed, and $\dot{M}$ the mass-loss rate. For this star, the mass-loss rate from the cmfgen models require $B>0.1$ MG for $\eta>1$, a modest magnetic field for such stars.

The hydrogenic O viii lines can provide strong constraints on magnetic field strength due to Zeeman splitting of the line into $\pi$ and $\sigma$ components. For moderate fields, the central $\pi$ component remains unchanged in wavelength, whilst the two $\sigma$ components shift to either side. The relative strength of the $\pi$ component depends on the magnetic field’s orientation angle: it is absent when the magnetic field points at the observer and it is equal in intensity to the sum of the $\sigma$ lines if the magnetic field is perpendicular to the line of sight. Each of the three Zeeman components can be considered as a Gaussian of width derived from the temperature. The observed line is the sum of these three Gaussians. For small shifts much less than the width of each line, this sum is itself a Gaussian widened by the shift. For larger shifts the line profile deviates.

The shift $\Delta\lambda$ between the two $\sigma$ lines can be parametrized as $\Delta\lambda=\alpha B$, where $\alpha$ is a constant for each ion and transition. The value for $\alpha$ is tabulated by Blom & Jupén (2002, their Table 1). For transitions considered here, that is O viii 8-9, O viii 9-10, and C vi 7-8, and for field strengths in Gauss, the $\alpha$ parameter is $1.8\times 10^{-5}$, $3.4\times 10^{-5}$, and $2.6\times 10^{-5}$, respectively. The 6064 Å O viii line has a FWHM of 93 Å. Assuming a separation between the $\sigma$ lines equal to the FWHM would have strongly distorted the shape, we constrain $B<2.7$ MG. The 4340 Å line with FWHM of 69 Å gives a conservative upper limit of 3.8 MG.

We ran models adding the Zeeman components assuming all three Zeeman lines have equal intensity. This ratio corresponds to a magnetic field orientation to the line of sight of 55 degrees. We then fit the profile of the 6064 Å line, where the intrinsic width is a fitting parameter. For $B=1.2$ MG an acceptable fit to the line shape is found. For $B=1.5$ MG the line profile deviation was clearly visible at 6064 Å, and developed a flat peak. Running the same test on the 4337 Å line shows it is indeed a bit less sensitive. We also ran the test with the central $\pi$ line twice as strong as each accompanying $\sigma$ line corresponding to a magnetic field oriented along the line of sight. Here the line remains centrally peaked which makes it easier to fit the profile. A notable deviation from the observed profile of the 6064 Å line was found for a field of 2.5 MG. At this field strength, half the line width is due to Zeeman splitting.

The O viii line shapes give an upper magnetic field limit of $B<2.5$ MG, easily within the requirement of the magnetically-confined toroidal structure proposed by Reindl et al. (2019). Based on the available data we can not rule out the alternative of such a structure.

This new field strength limit is far below that proposed by Gvaramadze et al. (2019, 200 MG) and of typical low-mass, magnetic WDs (10 MG). There is some caution here. If the hot gas is located in the wind, then the surface magnetic field may be underestimated because the field strength of a magnetic dipole will go down as $1/r^{2}$ with distance $r$. Also, the star has a very high luminosity and temperature, and so a larger radius than a typical WD. If it evolves as a WD it will contract significantly, and the field would strengthen. Hence, our limits are not inconsistent with known field strengths of magnetic WDs, which are of order 10 MG but they do not support values as high as 200 MG previously proposed.

The stellar remnant (Gaia DR3 526287189172936320) has been observed by various terrestrial and space-based all-sky photometric surveys over many decades. Some offer short cadences to detect transient events (e.g., TESS 2-min and 30-min cadence), while others cover decades of sparse photometric sampling with long exposures (e.g, DASCH). Many have variable depth and sensitivity and large effective detector pixel sizes and large photometric apertures. This can result in blending and nearby source contamination.

A faint, foreground star (Gaia DR3 526287189166341504) is 2.3″ eastwards of the stellar remnant (Figure 8). Though this star is clearly distinguished in short-exposure IPHAS (Drew et al., 2005) images, it is $\sim 4.5$ mag fainter in Gaia filters and well below limits of earlier photographic imagery so should not affect recorded optical photometry. A very bright star (Gaia DR3 526287257892412928) is $\sim$21″ directly west. Because recent all-sky surveys use large photometric apertures (e.g., the ASAS-SN PSF is $\sim$15″, or single-pixel aperture $\sim$17″ in OMC), the bright star ($G\sim 11$ mag) can affect adjacent pixels and influence photometric results. However, it is well separated from the stellar remnant in scanned data from all older photographic plates recorded in DASCH and so it is unlikely any variability for the stellar remnant is affected by the brighter star in 20^th century observations.

Here, we show $B$ band results from various surveys as this filter has the best temporal coverage. The data include photographic photometry from DASCH, the Ukrainian FON astrographic catalogue (Andruk et al., 2016) and the Palomar Observatory Sky Surveys (POSS-I and POSS-II); CCD data from the AAVSO Photometric All-Sky Survey (APASS DR10, Henden et al., 2018), the INT Galactic Plane Survey (IGAPS/UVEX; $g,r$; Monguió et al., 2020), the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS1, $g,r,i$; Tonry et al., 2012), the Asteroid Terrestrial-impact Last Alert System survey (ATLAS, $o,c$; Tonry et al., 2018), the Transiting Exoplanet Survey Satellite (TESS; Ricker et al., 2015), and the Near-Earth Object Wide-field Infrared Survey Explorer Reactivation Mission (NEOWISE, $3.6,4.5\rm\mu m$; Mainzer et al., 2014). The single epoch IGAPS/UVEX $g,r$ photometry from Pan-STARRS1, and the Pan-STARRS1 $g,r,i$ photometry, were converted to $B$ using the transformations of Tonry et al. (2012). The optical light curves are shown in Figure 9 and the infrared light curve is shown in Figure 10.

The stellar remnant is faint and the optical light curve sparsely populated over much of its recorded history (Figure 9) but evidence of dimming is seen. The photographic $B_{j}$ magnitude estimates, calibrated via the GSC 2.3.2 (Lasker et al., 2008) provided by DASCH for measurements between 1924.8 and 1950.7, show a mean $B_{j}$ magnitude of $15.6\pm 0.2$ mag, and maximum magnitude variation of 0.85 from 37 data points with a dimming slope of 0.02 magnitudes/year and an $R^{2}=0.38$ (which indicates a weak though positive correlation). Exposures where the recorded plate limit was within 0.3 mags of the star’s value are omitted as unreliable. Seven nearby stars, both fainter and brighter (including the bright star to the west), have much lower $R^{2}$ values (average of 0.047 with rms of 0.064) and are flat (average slope of 0.0007 mag/yr with rms also of 0.007). See the top panel in Figure 9. Overall, the temporal coverage is rather sparse apart from the ATLAS data¹¹¹¹11Which has a pixel size of $\sim 2$″ in two bands: the “cyan” (c) band (4200-6500Å) and the “orange” (o) band (5600-8200Å).. For the ATLAS dataset (before binning), no significant peaks were found in Lomb-Scargle periodograms so the source appears stable on short timescales in optical bands.

We extracted archival photometry from the Wide-field Infrared Survey Explorer (WISE) 3-band Cryo data release (Wright et al., 2010), and the NEOWISE latest data release. Figure 10 shows the averaged photometry per epoch in bands $W1$ (3.4µm) and $W2$ (4.6µm). Each epoch covers about six months. Frame inspection showed there may be bright star contamination within the standard photometric aperture of NEOWISE (radius 8.25″). There is no nebula contamination as this is not visible at those wavelengths. The $W1$ and $W2$ light curves show a modest decline of $\sim$0.01 mag/year over 11 years, half the rate seen at earlier epochs in the optical $B$-band.

The wind model indicates some weak but no strong emission features in the WISE $W1$ and $W2$ filters: the emission is dominated by continuum. Since the Pa 30 nebula is not visible at these wavelengths, we can exclude nebular lines such as [Fe ii] and [Ar vi] or the 3.3µm PAH band. Therefore the potential decline would be likely due to the continuum.

Kashiyama et al. (2019) propose that the central star (WD J005311 in their paper) is a highly-magnetic WD spinning at an angular frequency of 0.2 to 0.5 s^-1. This would suggest very short light-curve variability at a period of 12.6 to 31.4 seconds. The current data does not trace this time scale.

In Section 3.1 we hinted at variability in the stellar wind. The SparsePak/WIYN data from 2014, our 2016 OSIRIS/GTC spectra, the 2017 spectrum by Gvaramadze et al. (2019), and the averaged spectra in 2020 by Garnavich et al. (2020) are examined for long-term spectral variability (Figure 11). The SparsePak/WIYN fibre is 4.7″ and should include all stellar flux assuming the fibre is well placed. Comparison of the synthetic photometry from this spectrum (using the Python package pyphot¹²¹²12https://github.com/mfouesneau/pyphot) to the photometric measurements provided by the individual surveys justifies this assumption. Hence, all slit spectra are scaled to match the SparsePak/WIYN spectrum. To improve the signal-to-noise ratio the SparsePak/WIYN and Garnavich et al. (2020) spectra have been smoothed with a mean filter of length of 7 and 9 pixels respectively.

The most striking differences are intensity changes in the O vi 3811/34 emission line, as seen in panels (a) and (b) of Figure 11. The discontinuity in the Gvaramadze et al. (2019) spectrum at 3850 Å and the different spectral shape further to the blue are to be taken with caution, since flux calibration is less reliable in the far blue. Garnavich et al. (2020) found that the O vi 3811/34 line shows short-term $<2$hour variability, which was attributed to clumpiness in the stellar wind. Calibration differences among the different instruments do not allow a direct comparison of wind variability within the O vi 3811/34 line, but the OSIRIS/GTC spectra overall shape appears consistent with the averaged spectra of Garnavich et al. (2020). A decrease in the peak emission of the O vi 5300 Å line between 2014 and 2017 (panel (e) in Figure 11) is also seen, perhaps associated with a change in the stellar wind. Small changes in continuum slope beyond 5500 Å can be due to instrumental and calibration effects between epochs, as the photometry appears to have remained constant throughout (cf. Sect. 4). Future photometric monitoring observations should be obtained at cadences of 10 seconds or less in the visual and at angular resolutions sufficient to allow separation of the target from the adjacent bright star (e.g., $\leq 5$″).

Although no hydrogen lines were detected in the nebular spectrum from large aperture telescopes (Paper I), the outer Pa 30 nebula is tentatively detected in deep, $\approx$1 Rayleigh sensitivity¹³¹³13Compared to the $\simeq$3 Rayleighs for IPHAS where only the central star is visible. H$\alpha$ imaging from the low angular resolution (pixel size 1.6′) Virginia Tech “H$\alpha$” Spectral-line Survey (VTSS; Dennison et al., 1998) that revealed a single enhanced pixel at Pa 30’s position. The 5-$\sigma$ excursion within a 15′$\times$15′ region indicates a low-surface brightness hydrogen (H$\alpha$) shell of similar size to the [O iii] shell found in Paper I. The VTSS survey limit corresponds to $5.661\times 10^{-18}$ erg s^-1 cm^-2 arcsec^-2 at H$\alpha$. The continuum-corrected surface brightness within a 1 pc radius of Pa 30 is $\approx 5.6$ Rayleigh $\approx 8.8\times 10^{-13}$ erg s^-1 cm^-2. Adopting the extinction law of Howarth (1983) for H$\alpha$, $c_{H\alpha}=0.99\,E(B-V)$, then the dereddened H$\alpha$ brightness becomes $7.4\times 10^{-12}$ erg s^-1 cm^-2 for $A_{V}=2.9$, which translates to $L_{\rm H\alpha}\leq 1.2\,\rm L_{\odot}$. Deeper narrow-band H$\alpha$ imagery is needed to confirm this tentative detection.

Pa 30 is visible at mid- and far-infrared wavelengths though not detected in the IRAS 12µm maps. In the far-infrared, the shell is round as indicated by AKARI 65 to 140µm maps, extending from 88″ up to 128″($\pm 5$″) from the central star, or else a radius of about 0.98 – 1.42 pc at the adopted distance. Foreground interstellar emission confuses the images beyond 100µm. We attribute the far-infrared wavelengths to dust emission and the mid-infrared emission to line emission.

Given the infrared shell sizes and the uncertainties in archival point-source photometry, we have re-estimated surface brightnesses for all IRAS and AKARI images. Images were downloaded from IRSA/IPAC with the same angular resolution (15″/pixel): HIRES maps for IRAS, and Far-Infrared Surveyor maps for AKARI. An 100″ radius aperture was selected for all maps. We use the conversion factor of Ueta et al. (2019) to derive flux densities. The final photometry beyond 100µm was corrected for interstellar contamination (Table 3). IRAS 12µm and AKARI 160µm flux densities are upper limits (non detection). The new measurements are given in Figure 12, along with archival photometry from Pan-STARRS1, Gaia, IPHAS, 2MASS, and WISE. The IRAS and AKARI fluxes integrated over the nebula indicate the shell flux peaks at $\lambda\sim$90µm, while the foreground emission is colder. Using the aperture photometry results for the nebula and the VOSA tool, we derive a shell blackbody temperature of 60 K (cf. Figure 12).

We ran a cloudy model for the shell assuming it is ionized gas with normal ISM abundances. We assume a blackbody star with $L_{*}=3\times 10^{4}\,\rm L_{\odot}$, $T_{*}=2.3\times 10^{5}\,\rm K$, surrounded by a shell with inner and outer radius of 0.8 and 1.2 pc, and density $n=17\,\rm cm^{-3}$. Standard ISM abundances and dust content were used as defined in cloudy, with a range of grain size and a dust-to-gas ratio of $3.8\times 10^{-3}$. This gives a total gas mass of $\sim 2.1\,\rm M_{\odot}$ and a dust mass of $8\times 10^{-3}\,\rm M_{\odot}$. The model predicts broadband fluxes listed in Table 3, in good agreement with the infrared measurements. At shorter wavelengths, the model flux is dominated by a few strong emission lines: [Ne vi] at 7.6µm, [Ne v] at 14.7µm and 25µm, and [O iv] at 25.9µm. The longer wavelength emission cannot be explained by emission lines and is from dust. The predicted H$\alpha$ flux is $2.2\times 10^{-11}\,\rm erg\,s^{-1}\,cm^{-2}$, which, after extinction correction, gives 13 Rayleighs, 2.3$\times$ the VTSS H$\alpha$ surface brightness. The model explains the infrared photometry but overpredicts the apparent H$\alpha$ surface brightness.

Alternatively, it can be assumed that the H$\alpha$ emission comes from the accelerated gas seen in the [S ii] emission reported in Paper I. The average density of this component, derived from the [S ii] lines, is $N_{\rm e}=120\,\rm cm^{-3}$. We assume a temperature $T_{\rm e}=10^{4}$ K and use the relation (Pottasch, 1984):

where $F({\rm H}\beta)$ is in units of $10^{-11}$ erg cm^-2 s^-1, distance $d$ in kpc, electron temperature $T_{e}$ in K, and electron density $N_{e}$ in cm^-3. The intrinsic H$\beta$ flux, derived from the H$\alpha$ flux assuming recombination Case B for a theoretical H$\alpha$ to H$\beta$ ratio of 2.85, is $F(\rm H\beta)\approx 7.2\times 10^{-12}$ erg cm^-2 s^-1. Hence, an ionized mass of $\approx 0.3\,\rm M_{\odot}$ is derived.

The mass is therefore dominated by the outer shell, rather the accelerated gas seen in [S ii]. The overprediction of H$\alpha$ for the outer shell may indicate that this gas is hydrogen-poor, in which case the mass is less than the derived 2.1 M_⊙.

The presence of dust in the circumstellar shell is strongly supported. The mass is dependent on composition and grain size. The dust origin remains open: it can be supernova dust, interstellar dust heated by the star, or dust from previous mass loss. The circumstellar shell could represent swept-up ISM or matter ejected in a mass-loss episode before the supernova explosion, such as a relic planetary nebula.

Oskinova et al. (2020) estimate the amount of hydrogen-free, X-ray-emitting nebula gas to be $0.1~{}\rm M_{\odot}$. Our ionized mass estimate from the tentative H$\alpha$ detection is an additional 0.3 to 2.1 $\rm M_{\odot}$ (Sec. 5.2). However, the relation of this gas to the supernova is unclear (see next subsection) and its kinematics are not known.

Consider the kinetic energy of the H-poor ejecta, $E_{\rm kin}=(1/2)\,M_{\rm ej}v^{2}_{\rm exp}$, where $v_{\rm exp}$ is the ejecta expansion velocity. In Paper I, we find the shocked emission has a radial velocity of 1100 km s^-1. If the extended halo-like feature in the AKARI 90µm data (radius $128\pm 5$ arcsec) is the farthest extent of the nebula, then the ejecta velocity would be $\sim$1700 km s^-1  not much different to the shocked emission. Both values are consistent with the general decline of ejecta velocities more than 100 days post-eruption in Type Iax SNe (e.g., Kawabata et al., 2018; Srivastav et al., 2021).

Using the above ejecta mass and velocities, the kinetic energy of the ejecta $E_{\rm kin}$ ranges from $1\times 10^{48}$ to $3\times 10^{48}$ erg. These are conservative estimates since this is a lower limit for ejected mass, while the expansion velocity at the time of eruption could have been higher.

We compare these tentative estimates against known extragalactic Type Iax SNe in Figure 13. These include theoretical models by Lach et al. (2022, open diamonds), the proposed range of values for Type Iax by Foley et al. (2013, black dashed bar), and a few examples of extragalactic sources: 2008ha and 2010ae (gray; Stritzinger et al., 2014), 2007gd (blue; McClelland et al., 2010), 2014ck (orange; Tomasella et al., 2016), 2019gsc (magenta; Srivastav et al., 2020), and 2020kg (green; Srivastav et al., 2021). Both the low $E_{\rm kin}$ and ejecta mass estimate are consistent with observed ranges for Type Iax sources. However, they do not agree with the model predictions for pure deflagrations of Chandrasekhar-mass CO WDs of Lach et al. (2022).

The stellar luminosity and temperature place the star on a location in the HR diagram similar to post-AGB stars. In these stars, the stellar luminosity comes from residual hydrogen burning and follows the core mass–luminosity relation, originally proposed by Paczyński (1970). Using the star’s luminosity determined as $L_{*}=(40\pm 10)\times 10^{3}\,\rm L_{\odot}$, this relation gives a core mass (which is essentially the mass of the remnant star) of $M_{*}=1.20\pm 0.17\,\rm M_{\odot}$. The relation from Vassiliadis & Wood (1994) gives 0.1 M_⊙ less.

The application of these models to post-merger evolution can be queried. In the absence of detectable hydrogen, helium burning may be expected. The models of Vassiliadis & Wood (1994) find a similar luminosity for stable helium burning, however this helium very quickly runs out in post-AGB stars and the luminosity subsequently declines.

Merger models published by Schwab et al. (2016) and Schwab (2021) reach a similar location in the HR diagram as carbon burning stars, but take 10 to 20 kyr to reach this. They reach the post-AGB phases at 30% lower luminosity than predicted by the Vassiliadis & Wood (1994) for their mass. Lower mass merger models by Schwab & Bauer (2021) take 170 kyr to reach the post-AGB phase and similar models from Wu et al. (2022) reach the right luminosity but do so as cool red giant stars. All these merger models spend a long time as red giants. The models assume a red giant mass loss recipe. The evolution could conceivably be much accelerated if an instantaneous mass loss is assumed, either due to the explosion or akin to common envelope systems, thereby skipping the red giant phase and the early post-AGB phase. Accretion may also affect the evolution in the post-AGB region of the HR diagram.

Assuming the change in luminosity by 30%, the mass of the star would become $1.45\pm 0.2\,\rm M_{\odot}$. The mass range does not require the star to have super-Chandrasekhar mass as proposed by Gvaramadze et al. (2019), although this is not ruled out either.

The range of possible masses fits with WD merger products. The CO WD mass distribution peaks at approximately 0.6 M_⊙, while helium WDs are approximately 0.45 M_⊙. Two merging CO WDs or a CO+He WD can therefore reach the minimum remnant mass. O-Ne WDs are more massive and a merger involving one (Gvaramadze et al., 2019) would lead to a more massive remnant. The moderate neon abundance may not support involvement of such a WD.

The mass ejected in the merger should be included in the progenitor masses. The minimum ejecta mass is that derived from the X-ray emission (Oskinova et al., 2020) of 0.1 M_⊙. The outer shell may contain further mass ejected during the merger. The hydrogen mass in the shell cannot come from the merger as WDs have very thin hydrogen layers and this should not be be counted in the ejecta mass. The hydrogen is seen only in the VTSS detection and needs confirmation. If the outer shell is hydrogen-free, then it may add another 0.1 M_⊙. The merger ejecta are therefore estimated at $M_{\rm ej}=0.15\pm 0.05\,\rm M_{\odot}$.

This leaves the question of where the outer shell originates. If there is hydrogen in the outer shell, it can be swept-up ISM or from mass lost by the WD progenitor. If so the faint, outer shell is a relic planetary nebula, re-ionized by the luminous, hot post-merger star.

Pa 30 is only the second detection of dust emission in Type Iax SNe and the only known case of cold dust. The other example is the recent Type Iax SN 2014dt in M 61. Spitzer observations by Fox et al. (2016) found excess mid-infrared emission (3.6 and 4.5$\rm\mu m$) one year post-eruption in 2014dt, that could be a pre-existing dusty circumstellar matter or newly formed dust ($\sim 10^{-5}\,\rm M_{\odot}$ of carbonaceous material). Coincidentally, 2014dt has been suggested as a possibly bound remnant (Kawabata et al., 2018). We note that cold dust, as in Pa 30, would not be detectable in extragalactic SNe.

Theoretical models of Type Iax explosions predict a kicked stellar remnant (e.g., Jordan et al., 2012; Kashyap et al., 2018). However, the radial velocity of Pa 30’s stellar remnant, as well as the mean radial velocity of its surrounding nebula, is nearly zero, while the tangential velocity is low ($\sim 30$ km s^-1). The Gaia DR3 proper motion indicates the central star could have moved by at most $\sim 2$″ since 1181 AD. It clearly remains at the center of Pa 30. So these data are already a challenge to existing models. The nearby faint star (Figure 8) cannot be a kicked companion as this star is in the foreground ($1.1^{+0.4}_{-0.3}$ kpc; Bailer-Jones et al., 2021) with a much higher proper motion (21.285 mas yr^-1)

Lach et al. (2022) estimate kick velocities between 6.9 km s^-1 and 369.8 km s^-1 for bound remnants in Type Iax explosions. These could be in the ballpark for Pa 30.

The (only) CO+ONe WD merger model is by Kashyap et al. (2018) with a predicted kick velocity of $\sim 90$ km s^-1. However, an ONe WD may not be required here because of the remnant mass and the moderate neon abundance.