# Multi-scale PIIFD for Registration of Multi-source Remote Sensing Images

Chenzhong Gao, Wei Li  C. Gao and W. Li are with the School of Information and
Electronics, Beijing Institute Technology, Beijing 100081 China (e-mail:
liwei089@ieee.org).

###### Abstract

This paper aims at providing multi-source remote sensing images registered in
geometric space for image fusion. Focusing on the characteristics and
differences of multi-source remote sensing images, a feature-based
registration algorithm is implemented. The key technologies include image
scale-space for implementing multi-scale properties, Harris corner detection
for keypoints extraction, and partial intensity invariant feature descriptor
(PIIFD) for keypoints description. Eventually, a multi-scale Harris-PIIFD
image registration algorithm framework is proposed. The experimental results
of four sets of representative real data show that the algorithm has
excellent, stable performance in multi-source remote sensing image
registration, and can achieve accurate spatial alignment, which has strong
practical application value and certain generalization ability.

###### Index Terms:

Image Registration, Multi-Source, Remote Sensing, Scale-Space, Harris Corner,
PIIFD

##  I Introduction

Image registration has always been one of the most challenging works in image
processing, especially for multi-source images. In order to achieve in-depth
analysis, it is necessary to take advantage of the images obtained by
different sensors [1, 2, 3, 4], and the accurate registration of these images
is the fundamental premise. Multi-source image registration is an essential
step to implement multi-source data fusion, multi-source collaborative
classification, and joint analysis.

In general, satellite remote sensing data provide spatial calibration
information (spatial reference), such as latitude and longitude, map grid,
etc. But in practice, it is found that there can be problems of spatial
information loss. Besides, the spatial reference of different satellites or
sensors often has the problem of inconformity, with a difference of a few
pixels or even large deviations. It leads to inaccurate ground object
correspondence, which can cause a decline in the accuracy of image processing
results in data fusion, collaborative classification, change detection, image
mosaic, and other applications.

In practical work, remote sensing image registration is usually carried out
manually, by using ground control points, or through observation of the human
eye, which is time-consuming, laborious, and inefficient. The existing
automatic registration algorithm cannot meet the requirements of multi-source
remote sensing image registration in accuracy, stability, efficiency, and
adaptability. Many mature algorithms are dependent on the single-mode
properties of images [5, 6, 7, 8, 9]. However, multi-source remote sensing
image registration is a more complex problem in which image pairs are
generally obtained from different sensors, times, and viewpoints. Thus,
differences between the imaging principle, resolution, condition, and many
other factors can cause a multi-source remote sensing image pair to contain
changes in scale, blur, rotation, geometric and radiance distortion, content
details, and noise, etc, which invalidate the algorithms.

Multi-source remote sensing images are rich in types and have different
characteristics, bringing great difficulty to realize universal registration
algorithm. For example, hyperspectral image (HSI) and multispectral image
(MSI), which provide wealthy spectral information, and light detection and
ranging (LiDAR) data, which provide elevation information, are used together
for joint classification [1, 2, 3]. Correspondence analysis with optical and
synthetic aperture radar (SAR) images is investigated in depth [4]. In
addition, the fusion of infrared and optical images is also one of the most
important and commonly used technical means [10]. Panchromatic(PAN) images
with high spatial resolution and texture information are often used for fusion
with HSI and MSI. Therefore, it is in great need of effective, universal, and
robust algorithms for multi-source remote sensing image registration.

In this paper, focusing on the characteristics and differences of multi-source
remote sensing images, a robust feature-based image registration algorithm is
implemented, in which we use Harris corner detector for keypoints extraction,
scale-space to implement multi-scale properties, and PIIFD for key points
description. Eventually, a multi-scale PIIFD (MS-PIIFD) image registration
algorithm framework is proposed, which can solve the main problems and
difficulties: 1) different scales and blur, 2) radiance distortion, 3)
geometric distortion, 4) different content details. The experimental results
of 15 scenes of representative data show that the algorithm has excellent,
stable performance in multi-source remote sensing image registration, which
can achieve accurate spatial alignment, and has good universality and strong
practical application value.

##  II Related Works

Figure 1: The proposed multi-source remote sensing image registration
framework of MS-PIIFD.

Existing automatic image registration algorithms can be generally divided into
two categories, the traditional algorithms and the algorithms based on deep
learning. Traditional registration algorithms are systematically summarized
and classified into area-based methods and feature-based methods according to
technical means [11].

###  II-A Area-Based Methods

The area-based registration methods, also called intensity-based methods,
register images by directly using correlation of pixel values, mutual
information, or transform domain of the image pairs. Usually, the global image
or a predefined local window is used to search for matching [12]. By
establishing the similarity relationship of the original pixel values [13,
14], statistical properties [15, 16, 17], or transform domain, such as Fourier
transformation [18, 19], and wavelet transformation [20], the optimal spatial
transformation model between the image pairs is found through parameter
optimization methods. The advantage of the area-based methods lies in high
registration accuracy, but the drawbacks are particularly obvious. The
drawbacks of the area-based methods mainly lie in three aspects.

The first drawback is that this type of method only works for simple spatial
changes, such as translation and rotation, which are called rigid
transformations. If the image pairs contain more complex spatial changes, such
as scale difference, affine or perspective transformation, viewpoint changes,
and various nonlinear geographic distortions, it is almost impossible to
calculate the transformation parameters by comparing the similarity of the
information in a window, which will make the methods fail.

The second drawback, which is the most serious, is that measuring similarity
based on pixel values makes such algorithms too dependent on the modality of
the images. In other words, area-based methods cannot handle intensity
differences. For multi-source images, especially multi-sensor images, because
of the different purpose and principle of image collection, their pixel values
represent different physical information, so that there are multi-modal
properties between them. Take optical and thermal infrared image registration
as an example. The former not only has higher spatial resolution and clearer
edge of the objects but also represents a lot of texture information. The
latter has a lower resolution, blurry edges, and no complex content. The pixel
values between the image pairs may have non-linear differences or be totally
different. For multi-temporal images, the image content is inconsistent,
resulting in the similarity measurement being meaningless. Furthermore, noise
also has a great influence on these methods.

Another drawback is that, when using a local window to search for matching,
the content in the window should be weakly textured. If the area contained in
the window is relatively smooth, there will be a lack of obvious features for
matching. This problem also leads to a lot of redundant computations, making
the algorithms less efficient.

The area-based methods are not good choices for multi-source remote sensing
image registration.

###  II-B Feature-Based Methods

The feature-based registration methods detect images’ salient structure to
extract features as sparse representation, which can greatly reduce the
computational complexity. Moreover, the feature-based methods are more robust
to the scene difference and noise interference of the images, being relatively
efficient.

The commonly used features in image registration include point, line, edge
profile, and area features. Among them, compared with other features, the
coordinate position of point features can be directly used to calculate the
parameters of the transformation model, while other features often need to be
transformed into control points. For example, the line feature is represented
by the endpoints. Therefore, the point feature is an ideal matching primitive
[21]. The feature-based, especially point feature-based algorithm has always
been the focus of image registration research and is relatively mature at
present. The point feature-based registration algorithm involves the following
steps: 1) feature points detection, 2) feature description, 3) feature
matching, 4) mismatch removal, 5) parameter estimation and transformation.

Scale-invariant feature transform (SIFT) [22, 23, 5] is one of the most
classic, effective, and commonly used feature point matching methods. It first
constructs a Gaussian scale space, supported by scale-space theory [24], and
extracts local extreme points in a DoG (Difference of Gaussian) pyramid as the
feature points. Then, feature points are described by statistical gradient
histograms. SIFT has a good scale and rotation invariance. Many improved
algorithms based on SIFT were proposed in recent years. [25, 26, 6, 27, 28,
29, 30] Harris et al [31] proposed a feature point detection method, which
detects the corner structure that changes sharply in multiple directions in an
image by calculating corner response according to the gradient. It is an
effective and stable feature points detection algorithm, which is widely
applied in feature-based image registration. The oriented FAST and rotated
BRIEF (ORB) [7] uses improved accelerated segment test (FAST) [32] and
directional binary robust independent elementary features (BRIEF) [33]
respectively for feature points detection and description. The algorithm has
low time complexity, which is suitable for real-time registration. Binary
robust invariant scalable keypoints (BRISK) [8] and fast retina keypoint
(FREAK) [9] are the improved algorithms of BRIEF, which are scale invariant,
rotation invariant, and robust to noise. KAZE feature[34] uses the nonlinear
diffusion filtering method to establish the scale-space without blurring the
image, different from SIFT.

However, most algorithms are only applicable to single-modal image
registration. With the increasing requirements of multi-source image
registration, several algorithms for multi-modal image registration were
proposed. Based on SIFT algorithm, Chen et al [9] proposed a partial intensity
invariant feature descriptor (PIIFD) for multi-source retinal image
registration, which has intensity and rotation invariance but is unable to
deal with scale differences. PSO-SIFT [35] introduces a new gradient
definition and an enhanced feature matching method for the registration of
remote sensing images with intensity differences. Radiation-variation
insensitive feature transform (RIFT) [36] proposes a maximum index map (MIM)
for feature description, which has high intensity invariance, but low rotation
invariance, and is not scale invariant. Most of the research on multimodal
image registration is focused on medical images. However, except for some
differences in image characteristics and details, the multi-source remote
sensing image registration has little difference in essence with medical image
registration, so the algorithms also have universality.

###  II-C Deep Learning-Based Methods

In recent years, research has attempted to achieve image registration by using
deep learning methods. In general, relevant methods can be divided into two
categories [37]: using deep learning networks to estimate the similarity
measure of two images to drive iterative optimization, and directly using deep
regression networks to predict transformation parameters. Deep learning-based
methods first achieved good results in medical image registration. More and
more research on remote sensing image registration began to shift its focus to
the neural network and deep learning. Although the algorithms based on deep
learning are relatively advanced, they often need plenty of annotated samples
for model training, in which the model is highly targeted. So the deep
learning-based methods cannot be widely applied to all kinds of multi-source
data, and are not very mature at present.

##  III Proposed Registration Method

Figure 2: The results of three feature point detection methods in an optical-
infrared image pair.

Figure 3: The Harris corner detection results with and without LNMS in an
optical-infrared image pair.

The process of the MS-PIIFD registration algorithm is shown in Fig. 1. The
input is the image pair to be registered. The preprocessing is carried out
first, including simple denoising and data normalization. The images used to
obtain features are single-band grayscale images. For multi-band images, such
as HSI and MSI, since the proposed algorithm is robust to radiation
distortion, it is optional to select one of the bands as the input. Harris
corner point detection is carried out on the pre-processed image. Then the
scale-space of each image is created, in which a Gaussian pyramid is built.
The PIIFD descriptor of each detected Harris corner point is established in
the scale space. The feature points are matched according to the descriptors,
and then the mismatch removal is carried out. Finally, the spatial
correspondence of the matching pairs is used to calculate (estimate) the
parameters of the transformation model between the image pair.

###  III-A Feature Point Extraction

To ensure the accuracy of registration, the feature points should be highly
stable and repetitive between the image pair. For multi-source remote sensing
images, the difference between scale, blur, intensity, or content details can
lead to instability of local features. The feature point detection in SIFT and
FAST algorithms is too sensitive to the intensity. Harris corner detection
[31] focuses on locating obvious ’corner’ structure, which is distinct and
stable even in multi-model images. Fig. 2 shows the feature point detection
results of the three methods in an optical-infrared image pair. By comparison,
the feature points obtained by the Harris corner detector are the most highly
repetitive. Harris corner detector also has the advantages of strong
robustness to local noise and high computation efficiency. Therefore, the
Harris corner detector is adopted to extract feature points from multi-source
remote sensing images.

Harris corner point is detected by calculating the corner response of each
pixel:

$cornerness=\frac{{{\rm{det}}(\textbf{\emph{M}})}}{{{\rm{tr}}(\textbf{\emph{M}})}}$
(1)

$\textbf{\emph{M}}=\begin{bmatrix}\sum\limits_{(u,v)\in\bf{W}}\bf{I}_{x}(u,v)^{2}&\sum\limits_{(u,v)\in\bf{W}}\bf{I}_{x}(u,v)\bf{I}_{y}(u,v)\\\
\sum\limits_{(u,v)\in\bf{W}}\bf{I}_{x}(u,v)\bf{I}_{y}(u,v)&\sum\limits_{(u,v)\in\bf{W}}\bf{I}_{y}(u,v)^{2}\end{bmatrix}$
(2)

where $\rm{det}(\textbf{\emph{M}})$ is the determinant of M,
$\rm{tr}(\textbf{\emph{M}})$ is the trace of M, ${{{\bf{I}}_{x}}(u,v)}$ and
${{{\bf{I}}_{y}}(u,v)}$ are the gradient of the image in the $x$ and $y$
directions, ${\bf{W}}$ is a gaussian window.

One key point worth noting is that multi-source remote sensing images may have
serious intensity distortion or scale difference, which results in the
nonuniform distribution and low repeatability of feature points in the image
pair. This problem is visually illustrated in Fig. 3, which is the Harris
corner detection result of an optical-infrared image pair with a large scale
difference. The result of directly filtering with a threshold of the optical
image is shown in Fig. 3(a), in which most feature points are concentrated in
the left part of the image. However, the approximate position of the infrared
image in the visible image is indicated by the red box, and only a small
number of feature points are distributed in it, which are very nonuniform. The
proposed algorithm adopts local non-maximum suppression (LNMS) to solve this
problem. From Fig. 3(c), it can be seen that, with LNMS, the feature points
are far more uniformly distributed, and the repeatability is higher.

###  III-B Multi-scale PIIFD

PIIFD [38] is an improved algorithm based on SIFT [5] to solve the problem of
the full inverse (light and dark flipping) and nonlinear intensity difference
in multi-source retinal images. Similar problems are evident in multi-source
remote sensing images, which can be seen through the image pairs in Fig. 2 and
Fig. 3. So in this paper, the PIIFD algorithm is used to describe the feature
points.

In the original PIIFD [38], at first, each feature point is given a main
orientation (0 $\sim\pi$), by calculating the local second-order gradient,
which ensures rotation invariance. Then the neighborhood of each feature
point, the size of which is fixed at $40*40$ pixels, is used for description.
Taking the main orientation as the datum (0∘), the gradient within the
neighborhood is counted. Finally, a descriptor vector of length $128$ is
generated for feature point description. However, for multi-source remote
sensing images with different scales, the information in the local
neighborhood can have a large difference, which will lead to different
descriptors of the same point.

Figure 4: Structure of the scale space (gauss pyramid) used in MS-PIIFD.

In order to realize multi-scale PIIFD, local information of feature points
needs to be extracted in the scale space of the images. The method of building
image’s Gaussian pyramids based on the scale-space theory [24] is adopted. The
schematic diagram of establishing Gaussian pyramid of image in the proposed
algorithm is shown in Fig. 4. The original image is sampled down step by step
to obtain a series of images with different scales, that is, the first layer
in each octave. Then a series of Gaussian blurs are performed on each octave:

${\bf{L}}={\bf{G}}*{\bf{I}}$ (3)
${\bf{G}}=\frac{1}{{\sqrt{2\pi{\sigma^{2}}}}}{e^{\frac{{-({x^{2}}+{y^{2}})}}{{2{\sigma^{2}}}}}}$
(4)

where ${\bf{I}}$ is the original image, ${\bf{G}}$ is a Gaussian kernel with
standard deviation of $\sigma$, and ${\bf{L}}$ is the Gaussian blur image.

After the scale-space of the images is established, for each Harris corner
point, the multi-scale PIIFD descriptor is calculated by obtaining the
neighborhood information at the corresponding location in the scale space. The
algorithm flow is shown in Algorithm 1, where ${\bf{G}}_{\bf{I}}(x,y,O,L)$ is
the Gaussian scale space of image $\bf{I}$ established through down-sampling
and Gaussian blur, $O$ is the number of octave with different scale, $L$ is
the number of layers with different blurring degree. The algorithm outputs the
feature point descriptor set $D_{P}$, which contains 128-dimensional
descriptors for each feature point at each scale. The feature described in the
scale-space can reduce the influence of scale difference of multi-source
images and obtain the best matching result of feature points.

Algorithm 1 Proposed multi-scale PIIFD for registration

Input: feature point set $P_{\bf{I}}$, image scale space
${\bf{G}}_{\bf{I}}(x,y,O,L)$, PIIFD window scale $S$
Output: feature descriptor set $D^{P}$

1:In each layer of ${\bf{G}}_{\bf{I}}(x,y,O,L)$:

2: Calculate the second-order gradient of the image

3: Calculate the first-order gradient of the image

4: For each feature point $P$ in $P_{\bf{I}}$:

5: Calculate the corresponding position

6: Calculate the main orientation using second-order gradient

7: Taking the main orientation as the reference direction ($0^{\circ}$),
establish a square window with size (side length) of S, and divide it into 16
small squares

8: Statistics the first-order gradient direction within each small square to
obtain feature descriptor $D(P,O,L)$

9:Output feature descriptor set $D_{P_{\bf{I}}}$

###  III-C Matching and Transformation

After all the descriptors of each feature point are created, the best-bin-
first (BBF) [39] matching method is adopted to obtain the initial matches of
the feature points, when the cosine similarity is used as the measurement. For
multi-source images, there is likely to be a large number of mismatches in
initial matches, which causes common outlier removal algorithms to fail, such
as RANSAC (random sample consensus). In the proposed algorithm, the
consistency of the main orientation and spatial distribution of the feature
points [38] are used for mismatch removal. This is because although there may
be spatial distortion in the multi-source remote sensing image pairs, there
will not be too much deformation on the whole. Then, assuming that there is
only overall rotation or scaling between the image pair. So the angle between
the main orientation of each pair of correctly matched feature points should
be consistent, and their spatial relationship with other correctly matched
points in the two images should also be consistent. Therefore, this method is
suitable for multi-source remote sensing images.

Finally, three transformation models, similarity, affine, and projective
transformation, are selected for parameter estimation and image transformation
to deal with the different spatial distortion that may exist between the image
pairs, where the least square method is adopted to estimate parameters.

##  IV Experimental Results and Discussion

Figure 5: Registration results of the proposed MS-Harris-PIIFD in three forms, including matching results (left), fusion form (middle), and checkerboard form (right). TABLE I: Number of correct matching pairs comparison by five registration methods. Method | Scene
---|---
a | b | c | d | e | f | g | h | i | j | k | l | m | n | o
SIFT [5] | 22 | 14 | 112 | 5 | 0 | 0 | 12 | 108 | 6 | 9 | 10 | 0 | 16 | 0 | 0
PIIFD [38] | 0 | 0 | 0 | 0 | 0 | 0 | 48 | 0 | 0 | 0 | 12 | 180 | 0 | 18 | 17
PSO-SIFT [35] | 0 | 21 | 58 | 0 | 6 | 18 | 22 | 107 | 16 | 21 | 8 | 127 | 0 | 8 | 0
RIFT+ [36] | 0 | 7 | 7 | 3 | 2 | 2 | 190 | 0 | 158 | 504 | 76 | 0 | 0 | 410 | 214
Proposed | 181 | 16 | 383 | 90 | 63 | 114 | 175 | 76 | 55 | 271 | 155 | 293 | 112 | 82 | 113

To validate the effectiveness of the proposed MS-PIIFD, a distinct set of
multi-source image data are used to test the robustness and universality of
the proposed algorithm, and four state-of-the-art multi-source registration
algorithms including SIFT [5], PIIFD [38], PSO-SIFT [35] and RIFT+ [36] are
selected for comparison. Except for the feature point detection and
description method, all processing is the same, including image normalization,
noise reduction, feature point matching, and mismatch removal, etc. The
experiments are all implemented using MATLAB2021a on a platform with Intel
Core I7-8700 CPU.

###  IV-A Experimental Data

The data used in experiments are multi-source image pairs of 15 scenes,
labeled a$\sim$o, in which there are both multi-sensor and multi-temporal
images, including HSI, MSI, LiDAR-derived DSM (digital surface model), SAR,
optical images, infrared images, depth map, and artificially produced
rasterized map. These selected experimental data basically contain all the
problems in multi-source image registration, including large resolution
differences, serious radiance distortions, spatial distortion, rotation,
content details differences, and image noise. To further test the universality
of the proposed algorithm, not only satellite data but also UAV data and
ground scene data are selected.

###  IV-B Parameter Setting

The proposed MS-PIIFD algorithm contains four main parameters that can greatly
affect the registration, which are the size of the LNMS window and the
threshold of the number of feature points in the Harris corner detector, and
the total number of octaves and layers of the Gaussian pyramid.

In Harris corner points detection, since the positional relationship and scale
difference of the image pairs are unknown in advance, we hope that the feature
points are distributed as uniformly as possible in the image pair. So, the
ratio of the size of the LNMS window depends on the ratio of the pixel number
of the two images:

$ratio=\sqrt{\frac{{M\times N}}{{m\times{\rm{n}}}}}$ (5)

where $M,N$ and $m,n$ are the length and width of the two images. In addition,
the minimum unit size in the established PIIFD descriptor is $10\times 10$.
Thus, the minimum size of the LNMS window is set to $10\times 10$, to prevent
poor separability of the descriptors due to the feature points being
excessively close. High repetition of feature points is required to provide
sufficient matching candidates, so it is important to provide sufficient
points to be matched. But too many feature points not only increase the
computational burden but also bring a lot of outliers. Through experiments,
the threshold of the Harris corner points number is set between 500 and 1000.
The specific number is determined according to the size of the image. In this
experiment, the maximum number of points was set at 1000.

To ensure the efficiency of registration, it is necessary to set the suitable
total number of octaves and layers, when too large scale space will increases
the computational burden a lot. Through experiments, the total numbers of
octaves and layers are set to 3 and 4, which can basically handle most multi-
source image pairs.

###  IV-C Visual Results

The visual registration results of the proposed algorithm are shown in Fig. 3.
The results of each scene are respectively displayed in three forms: matching
results of feature points, fusion form, and checkerboard form of the image
pairs. Although the characteristics of scene contents of multi-source images
are not consistent, the alignment of the images can be evaluated by roughly
comparing the coincidence of objects’ edges.

By artificial comparison, there is no obvious deviation or ’virtual shadow’ in
the images, except for scene f. This is because the imaging quality of the
infrared image is low, and the image has a large spatial distortion,
especially at the edge of the image. Simply using similarity, affine, or
projection transformations will not align the image in all parts. In general,
some functional interpolation and fitting methods will be applied in this case
to achieve a better alignment effect, but it is not the focus of this paper,
and further research will pay attention to these issues. It can be found that
the spatial difference of images is within 1$\sim$2 pixels. Apart from the
edge positions of scene f, there is no obvious registration failure in any of
the scenes, and the registration algorithm is effective.

###  IV-D Numerical Results

The registration accuracy is usually measured by root mean squared error
(RMSE) and other similar methods, through comparing the similarity of the
image pair after registration. However, since the intensity of multi-source,
especially multi-modal images, is different, the methods of measuring
similarity do not have physical significance. Therefore, the registration
effect is evaluated by directly comparing the number of correct matches of the
feature points, because the number of correct control points is positively
correlated with the accuracy of the transformation model. When the number of
matched feature points is less than 3 pairs, the parameters of the spatial
transformation model cannot be solved, and the registration is considered to
be failed.

The numbers of correct matching pairs of the five registration methods are
shown in Table 1. In scene b, the proposed algorithm gets relatively few
matches, and 5 less than PSO-SIFT. Through analysis, it is because that the
texture information of the SAR image is relatively complex and chaotic, which
reduces the performance of the proposed algorithm, but still obtains enough
correct matches. In the scene g, h, i, j, n, o, the proposed algorithm does
not get the most matches, and in contrast, the Rift+ is doing better. It is
because there is basically no scale difference between the selected images of
these scenes. In particular, the scene i, j, n, o are image pairs that have
been artificially modified to the same scale. However, the proposed MS-PIIFD
still obtains a considerable number of correct matches. Among all the
competitive methods, the proposed registration algorithm is the only one that
does not fail in any scene and has significant advantages in most scenes.

##  V Conclusions

In this paper, an effective image registration algorithm based on MS-PIIFD was
proposed for multi-source remote sensing image registration. Through
experiments and analysis, the algorithm can well achieve the registration of
multi-source remote sensing images and performs well in practical work. The
image registration error can reach within 1$\sim$2 pixels, and the objects in
the scene can be aligned accurately. Through tests of various scene data
registration, it is concluded that the algorithm has good robustness,
stability, and universality, and has strong practical application value.

## References

  * [1] M. Zhang, W. Li, Q. Du, L. Gao, and B. Zhang, “Feature extraction for classification of hyperspectral and lidar data using patch-to-patch cnn,” _IEEE Transactions on Cybernetics_ , vol. 50, no. 1, pp. 100–111, 2018.
  * [2] X. Zhao, R. Tao, W. Li, H.-C. Li, Q. Du, W. Liao, and W. Philips, “Joint classification of hyperspectral and lidar data using hierarchical random walk and deep cnn architecture,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 58, no. 10, pp. 7355–7370, 2020.
  * [3] C. Liu, R. Tao, W. Li, M. Zhang, W. Sun, and Q. Du, “Joint classification of hyperspectral and multispectral images for mapping coastal wetlands,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , 2020.
  * [4] M. Zhang, W. Li, R. Tao, and S. Wang, “Transfer learning for optical and sar data correspondence identification with limited training labels,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , 2020\.
  * [5] D. G. Lowe, “Distinctive image features from scale-invariant keypoints,” _International Journal of Computer Vision_ , vol. 60, no. 2, pp. 91–110, 2004\.
  * [6] H. Bay, T. Tuytelaars, and L. Van Gool, “Surf: Speeded up robust features,” in _European conference on computer vision_. Springer, 2006, pp. 404–417.
  * [7] E. Rublee, V. Rabaud, K. Konolige, and G. Bradski, “Orb: An efficient alternative to sift or surf,” in _2011 International conference on computer vision_. Ieee, 2011, pp. 2564–2571.
  * [8] S. Leutenegger, M. Chli, and R. Y. Siegwart, “Brisk: Binary robust invariant scalable keypoints,” in _2011 International conference on computer vision_. Ieee, 2011, pp. 2548–2555.
  * [9] A. Alahi, R. Ortiz, and P. Vandergheynst, “Freak: Fast retina keypoint,” in _2012 IEEE Conference on Computer Vision and Pattern Recognition_. Ieee, 2012, pp. 510–517.
  * [10] J. Ma, Y. Ma, and C. Li, “Infrared and visible image fusion methods and applications: A survey,” _Information Fusion_ , vol. 45, pp. 153–178, 2019\.
  * [11] B. Zitova and J. Flusser, “Image registration methods: a survey,” _Image and vision computing_ , vol. 21, no. 11, pp. 977–1000, 2003.
  * [12] A. Gruen, “Development and status of image matching in photogrammetry,” _The Photogrammetric Record_ , vol. 27, no. 137, pp. 36–57, 2012.
  * [13] W. K. Pratt, “Correlation techniques of image registration,” _IEEE Transactions on Aerospace and Electronic Systems_ , no. 3, pp. 353–358, 1974.
  * [14] A. Mahmood and S. Khan, “Correlation-coefficient-based fast template matching through partial elimination,” _IEEE Transactions on image processing_ , vol. 21, no. 4, pp. 2099–2108, 2011.
  * [15] P. Viola and W. M. Wells III, “Alignment by maximization of mutual information,” _International journal of computer vision_ , vol. 24, no. 2, pp. 137–154, 1997.
  * [16] J. P. Pluim, J. A. Maintz, and M. A. Viergever, “Mutual-information-based registration of medical images: a survey,” _IEEE Transactions on Medical Imaging_ , vol. 22, no. 8, pp. 986–1004, 2003.
  * [17] F. P. Oliveira and J. M. R. Tavares, “Medical image registration: a review,” _Computer methods in biomechanics and biomedical engineering_ , vol. 17, no. 2, pp. 73–93, 2014.
  * [18] E. De Castro and C. Morandi, “Registration of translated and rotated images using finite fourier transforms,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , no. 5, pp. 700–703, 1987.
  * [19] X. Tong, Z. Ye, Y. Xu, S. Liu, L. Li, H. Xie, and T. Li, “A novel subpixel phase correlation method using singular value decomposition and unified random sample consensus,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 53, no. 8, pp. 4143–4156, 2015.
  * [20] J. Le Moigne, W. J. Campbell, and R. F. Cromp, “An automated parallel image registration technique based on the correlation of wavelet features,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 40, no. 8, pp. 1849–1864, 2002.
  * [21] A. A. Goshtasby, _2-D and 3-D image registration: for medical, remote sensing, and industrial applications_. John Wiley & Sons, 2005.
  * [22] D. G. Lowe, “Object recognition from local scale-invariant features,” in _Proceedings of the seventh IEEE international conference on computer vision_ , vol. 2. Ieee, 1999, pp. 1150–1157.
  * [23] D. G. Lowe, “Local feature view clustering for 3d object recognition,” in _Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. CVPR 2001_ , vol. 1. IEEE, 2001, pp. I–I.
  * [24] T. Lindeberg, “Scale-space theory: A basic tool for analyzing structures at different scales,” _Journal of applied statistics_ , vol. 21, no. 1-2, pp. 225–270, 1994.
  * [25] Y. Ke and R. Sukthankar, “Pca-sift: A more distinctive representation for local image descriptors,” in _Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004._ , vol. 2. IEEE, 2004, pp. II–II.
  * [26] K. Mikolajczyk and C. Schmid, “A performance evaluation of local descriptors,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 27, no. 10, pp. 1615–1630, 2005.
  * [27] J.-M. Morel and G. Yu, “Asift: A new framework for fully affine invariant image comparison,” _SIAM journal on imaging sciences_ , vol. 2, no. 2, pp. 438–469, 2009.
  * [28] A. Sedaghat, M. Mokhtarzade, and H. Ebadi, “Uniform robust scale-invariant feature matching for optical remote sensing images,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 49, no. 11, pp. 4516–4527, 2011.
  * [29] F. Dellinger, J. Delon, Y. Gousseau, J. Michel, and F. Tupin, “Sar-sift: a sift-like algorithm for sar images,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 53, no. 1, pp. 453–466, 2014.
  * [30] A. Sedaghat and H. Ebadi, “Remote sensing image matching based on adaptive binning sift descriptor,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 53, no. 10, pp. 5283–5293, 2015.
  * [31] C. G. Harris, M. Stephens, _et al._ , “A combined corner and edge detector.” in _Alvey vision conference_ , vol. 15, no. 50. Citeseer, 1988, pp. 10–5244.
  * [32] E. Rosten, R. Porter, and T. Drummond, “Faster and better: A machine learning approach to corner detection,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 32, no. 1, pp. 105–119, 2008.
  * [33] M. Calonder, V. Lepetit, M. Ozuysal, T. Trzcinski, C. Strecha, and P. Fua, “Brief: Computing a local binary descriptor very fast,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 34, no. 7, pp. 1281–1298, 2011.
  * [34] P. F. Alcantarilla, A. Bartoli, and A. J. Davison, “Kaze features,” in _European conference on computer vision_. Springer, 2012, pp. 214–227.
  * [35] W. Ma, Z. Wen, Y. Wu, L. Jiao, M. Gong, Y. Zheng, and L. Liu, “Remote sensing image registration with modified sift and enhanced feature matching,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 14, no. 1, pp. 3–7, 2016\.
  * [36] J. Li, Q. Hu, and M. Ai, “Rift: Multi-modal image matching based on radiation-variation insensitive feature transform,” _IEEE Transactions on Image Processing_ , vol. 29, pp. 3296–3310, 2019.
  * [37] G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, and C. I. Sánchez, “A survey on deep learning in medical image analysis,” _Medical image analysis_ , vol. 42, pp. 60–88, 2017.
  * [38] J. Chen, J. Tian, N. Lee, J. Zheng, R. T. Smith, and A. F. Laine, “A partial intensity invariant feature descriptor for multimodal retinal image registration,” _IEEE Transactions on Biomedical Engineering_ , vol. 57, no. 7, pp. 1707–1718, 2010.
  * [39] J. S. Beis and D. G. Lowe, “Shape indexing using approximate nearest-neighbour search in high-dimensional spaces,” in _Proceedings of IEEE computer society conference on computer vision and pattern recognition_. IEEE, 1997, pp. 1000–1006.


# MUSE observations of the optical nebula surrounding the central compact
object in the Vela Junior Supernova Remnant

Janette Suherli,1 Samar Safi-Harb,1 Ivo R. Seitenzahl,2 Parviz Ghavamian,3
Wynn C. G. Ho,4 Chuan-Jui Li,5 Ashley J. Ruiter,2 Ralph S. Sutherland,6 and
Frédéric P. A. Vogt7
1Department of Physics and Astronomy, University of Manitoba, Winnipeg,
Manitoba, R3T 2N2, Canada
2School of Science, University of New South Wales, Australian Defence Force
Academy, Canberra, ACT 2600, Australia
3Department of Physics, Astronomy and Geosciences, Towson University, Towson,
MD 21252, USA
4Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue,
Haverford, PA 19041, USA
5Institute of Astronomy and Astrophysics, Academia Sinica, No. 1, Sec. 4,
Roosevelt Rd. Taipei 10617, Taiwan
6Research School of Astronomy and Astrophysics, Australian National
University, Cotter Road, Weston Creek, ACT 2611, Australia
7Federal Office of Meteorology and Climatology – MeteoSwiss, Chemin de
l’Aérologie 1, 1530 Payerne, Switzerland E-mail<EMAIL_ADDRESS>
(Accepted 2023 November 27. Received 2023 November 27; in original form 2023
October 24)

###### Abstract

Central Compact Objects (CCOs), neutron stars found near the centre of some
Supernova Remnants (SNRs), have been almost exclusively studied in X-rays and
are thought to lack the wind nebulae typically seen around young, rotation-
powered pulsars. We present the first, spatially-resolved, morphological and
spectroscopic study of the optical nebula observed at the location of CXOU
J085201.4–461753, the CCO in the heart of the Vela Junior SNR. It is currently
the only Galactic CCO with a spatially coincident nebula detected at optical
wavelengths, whose exact nature remains uncertain. New MUSE integral field
spectroscopy data confirm that the nebula, shaped like a smooth blob extending
8$\arcsec$ in diameter, is dominated by [N ii]λλ6548,6583 emission. The data
reveals a distinct and previously unobserved morphology of the Hα emission,
exhibiting an arc-like shape reminiscent of a bow shock nebula. We observe a
significantly strong [N ii] emission relative to Hα, with the [N
ii]λλ6548,6583 up to 34 times the intensity of the Hα emission within the
optical nebula environment. Notably, the [N ii] and Hα structures are not
spatially coincident, with the [N ii] nebula concentrated to the south of the
CCO and delimited by the Hα arc-like structure. We detect additional emission
in [N i], He i, [S ii], [Ar iii], [Fe ii], and [S iii]. We discuss our
findings in the light of a photoionization or Wolf-Rayet nebula, pointing to a
very massive progenitor and further suggesting that very massive stars do not
necessarily make black holes.

###### keywords:

techniques: imaging spectroscopy – ISM: individual objects: RX J0852.0–4622 –
stars: neutron

††pubyear: 2023††pagerange: MUSE observations of the optical nebula
surrounding the central compact object in the Vela Junior Supernova
Remnant–MUSE observations of the optical nebula surrounding the central
compact object in the Vela Junior Supernova Remnant

##  1 Introduction

The nature and environment of Central Compact Objects (CCOs) constitute a
crucial and fundamental observational test of our understanding of the
evolution of massive stars and the theoretical modeling of the supernova (SN)
event process in core-collapse explosions. Typified by the CCO first
discovered with the Chandra X-ray Observatory in the Cas A supernova remnant
(SNR), CCOs have been identified near the centres of young core-collapse SNRs
with no pulsar wind nebulae (PWNe) associations, and have been studied almost
exclusively in the X-ray band. Their X-ray spectra are dominated by thermal
emission, with blackbody temperatures of 0.2–0.5 keV and typical luminosities
of $L_{X}\sim 10^{33}-10^{34}$ erg s-1 (Pavlov et al., 2004; see De Luca,
2017; Safi-Harb, 2017, for reviews).

Among the 15 known CCOs or CCO candidates111To access the complete list of
Galactic SNRs that host CCOs, users can apply a global filter by typing ”CCO”
on the SNRcat page: http://snrcat.physics.umanitoba.ca/SNRtable.php (Ferrand &
Safi-Harb, 2012)., three are detected as pulsars with inferred surface dipole
magnetic fields ($\sim$ 1010-11 G) and characteristic ages that are several
orders of magnitude larger than that of their associated SNRs. As such, CCOs
have been designated as anti-magnetars (Gotthelf et al., 2013) given their
relatively low surface dipole magnetic fields compared to the 1014-15 G
inferred fields of magnetars. Their unique properties have been attributed to
the invoked fall-back from the SN burying their magnetic field and the
magnetic field evolution in neutron stars (Ho, 2011; Rogers & Safi-Harb,
2016). One of the key open questions about CCOs is their SN progenitors. A
recent systematic X-ray study of the CCO-hosting ejecta-dominated SNRs showed
evidence for low-mass progenitors or a relatively low explosion energy (Braun
et al., 2019; Braun et al., 2023). In contrast, CXOU J085201.4–461753, the CCO
that is investigated in this article, has been suggested to be the product of
a massive progenitor (Slane et al., 2001; Allen et al., 2015).

CXOU J085201.4–461753 sits near the heart of Vela Junior (Vela Jr. or
G266.2–1.2 or RX J0852.0–4622; Slane et al., 2001; Mereghetti, 2001; Pavlov et
al., 2001), a Galactic SNR that was discovered 25 years ago through hard X-ray
observations (E $>$ 1.3 keV; Aschenbach, 1998). It is predominantly
characterized by non-thermal X-ray radiation and has been detected in high-
energy gamma-rays (E $\gtrsim$ 100 GeV; e.g. Aharonian et al., 2005, 2007; H.
E. S. S. Collaboration et al., 2018). Vela Jr. remains, thus far, one of the
iconic sources in our Galaxy for studying particle acceleration in SNRs. Its
age and distance were initially estimated to be $\sim$680 year and $\sim$200
pc, respectively, based on the direct detection of $\gamma$-ray emission
attributed to the decay of 44Ti (Iyudin et al., 1998), which remains
controversial (see e.g. Renaud et al., 2006). However, an age of 0.5-1 kyr and
a distance of 700 $\pm$ 200 pc, derived from the expansion rate of the remnant
(Allen et al., 2015), are now favoured in the literature.

A particularly enigmatic property of CCOs is their exclusive detection in
X-rays, with no confirmed counterparts detected in any other wavelength. This
has made the serendipitous discovery of an optical ring of emission in the
Small Magellanic Cloud SNR 1E0102.2–7219 surrounding a CCO candidate (Vogt et
al., 2018) especially remarkable. The true nature of this source however
remains uncertain (Hebbar et al., 2020; Long et al., 2020; Banovetz et al.,
2021; Li et al., 2021). To date, the Vela Jr. CCO stands as the only Galactic
CCO known to have an optical nebular counterpart. The optical nebula was
initially identified by Pellizzoni et al. (2002) and was then thought to be Hα
emission. It has a roughly circular shape with a diameter of $\sim$6″(Figure 1
– top left), and was proposed to be either a bow-shock nebula produced by the
CCO that moves supersonically through the interstellar medium (ISM), or a
photo-ionization nebula powered by the radiation from the CCO. Subsequent
observations using a broad R-band filter revealed that the shape of the nebula
was not exactly spherical, but rather resembles a kidney-bean or an arc-like
structure (Mignani et al., 2007, Figure 1 – top middle).

Figure 1: The comparison of 15″$\times$ 15″Hα images of the Vela Jr. CCO
optical nebula from previous observations and our MUSE observations in Hα and
[N ii]. The red circle marks the CCO location on each panel, and west of the
nebula is a field star, Star Z. Each frame colour corresponds to the filter
bandpass line colour in Figure 7. Top left: Digitized image of photographic
plate Hα image (Parker & Phillipps, 1998) from the 3.9 m UK Schmidt Telescope
(UKST) at the Australian Astronomical Observation (AAO) showing a nebula with
a diameter of roughly 6″(Pellizzoni et al., 2002). Top middle: R-band image
from VLT/Focal Reducer Spectrograph (FORS1) showing a kidney-bean shape nebula
(Mignani et al., 2007). Top right: HST/WFPC2 F656N observations in which the
optical nebula was not detected (Mignani et al., 2009). Bottom left: A
serendipitous OmegaCAM (Arnaboldi et al., 1998; Kuijken et al., 2002; Kuijken,
2011) observation using the narrow band H α filter NB_659 (Mignani et al.,
2019). Bottom middle: Hα integrated image generated from the non-continuum-
subtracted MUSE datacube, showing a clear arc-like shape structure. Bottom
right: [N ii] MUSE integrated image, showing a smooth blob-shaped morphology.

Narrow-band imaging observations were conducted using the Wide Field Planetary
Camera 2 (WFPC2) on the Hubble Space Telescope (HST) but these failed to
detect the optical nebula, down to a 3σ flux limit (Mignani et al., 2009,
Figure 1 – top right). This non-detection was attributed to either the F656N’s
narrower and bluer bandwidth compared to those of the filters used in previous
ground-based observations, or the possibility that the nebula is not dominated
by Hα emission after all. Follow-up spectroscopy observations revealed that
the spectrum of the nebula is instead dominated by strong and narrow [N
ii]λλ6548,6583, line emissions with only a minor contribution from Hα. This
led Mignani et al. (2019) to rule out the bow-shock nebula scenario and
proposed a plausible association between the nebula and a bright planetary
nebula candidate, Wray 16-30 (see Reynoso et al., 2006), located west of the
CCO, based on the similarities in their spectra.

In this article, we present the results from our study of the Vela Jr. CCO
field with Multi Unit Spectroscopic Explorer (MUSE; Bacon et al., 2010), the
integral field spectrograph mounted on the 8.2-m Unit Telescope 4 Very Large
Telescope (UT4 VLT) at the Cerro Paranal Observatory in Chile. The MUSE
datacubes allowed us to construct spatially-resolved maps of the optical
nebula’s morphology with sub-arcsecond resolution and derive diagnostic
measurements from the spectra. Section 2 describes our MUSE observations and
touches on the data reduction and analysis processes. Our results and
characterization of the optical nebula is presented in Section 3. We discuss
our findings in Section 4, and conclude the article in Sect 5.

##  2 Observations, Data Reduction, and Post-Processing

###  2.1 MUSE observations

We observed the central 2′ field of Vela Jr. SNR using MUSE in Wide Field Mode
with Ground Layer Adaptive Optics (WFM-AO) configuration (Arsenault et al.,
2008; Ströbele et al., 2012) in Service Mode (Program ID 0104.D-0092(B); P.I.:
F.P.A. Vogt). MUSE-WFM (nominal) mode covers a spectral range from 4800 Å to
9300 Å and samples the wavelength in 1.25 Å bins at spectral resolutions
ranging from R $\sim$ 1770 on the blue side to R $\sim$ 3590 on the red side.
It has a field of view of 1 $\times$ 1 arcmin-1 with a spatial sampling of 0.2
arcsec pixel-1. We arranged our MUSE fields to optimally cover the CCO (RA =
08${}^{\text{h}}$52${}^{\text{m}}$01$\aas@@fstack{s}$37, Dec =
-46°17′53$\aas@@fstack{\prime\prime}$50; Mignani et al., 2007) optical nebula
and most of the structures extending from Wray 16-30 to the northwest, while
also avoiding the bright star in the northeast.

To prevent any background diffuse nebular contamination in the sky frames, we
selected an empty sky field, away from the Galactic plane, at
07${}^{\text{h}}$57${}^{\text{m}}$15$\aas@@fstack{s}$5;
-62°41′00$\aas@@fstack{\prime\prime}$00\. We also applied some position angle
changes to each object exposure to minimize the residual background artefacts
from the MUSE’s 24 identical Integral Field Unit modules. Our observations
were carried out in 8 object-sky-pair observing blocks (OBs), for a total of
16 OBs, over 8 nights between 2020 (UT February 25, March 21 and 22222The
program was affected by the COVID-19 closure of the VLT with no science
operation on UT4 from March 23, 2020 to November 12, 2020, December 14, 15,
and 17) and 2021 (UT January 10 and 17). Each OB pair consisted of 3 $\times$
740 s exposure on the object and a 165 s sky exposure.

Figure 2: Top: The archival white-light image of the central 2′ field of Vela
Jr. SNR, obtained with OmegaCAM using the narrow-band Hα filter NB_659. The
Vela Jr. CCO, optical nebula, star Z of Pavlov et al. (2001), and Wray 16-30
are labeled. The white square indicates the 10″$\times$ 10″ cutout centred on
the CCO, corresponding to the area of the bottom right image. Bottom left: The
continuum-subtracted image from the MUSE WFM datacube, combining Hα (cyan), [N
ii] (magenta), and [S ii] (yellow). The image quality of the MUSE datacube is
0.46 arcsec at 5000 Å. Note that the components of the [N ii] and [S ii]
doublets are resolved spectrally in the MUSE data, and we use only the
emission of the brighter components for this figure. Bottom right: The
continuum-subtracted 10″$\times$ 10″ region, centred on Vela Jr. CCO position
(indicated by the red circle), highlighting the spatially resolved optical
nebula.  Figure 3: Narrow-band emission line intensity maps of all the bright
lines detected in the Vela Jr. CCO optical nebula (from the continuum-
subtracted datacube, and the intensity is in units of 10-20 erg Å-1 cm-2 s-1).
For [N ii] and [S ii] doublets, that are both resolved at the R$\sim$3000
spectral resolution of MUSE, only the maps of the brighter lines are shown.
Each map represents the 10″$\times$ 10″region centred on the CCO, with north
is up, and is overlaid with [N ii]λ6583 intensity contours (white dotted
lines). The red circles in all maps mark the positions of the Vela Jr. CCO,
and the black dashdot lines indicate the estimated angular distances from the
CCO location to the peak of each emission structure (see Table 1). The grey
filled-circles in all images are the subtracted region of Star Z, a field star
in the area (Pavlov et al., 2001).

###  2.2 Standard data reduction, calibrations, and post-processing

Our MUSE data were reduced individually using the automated MUSE pipeline
(v2.8.7; Weilbacher et al., 2016; Weilbacher et al., 2020) on EsoReflex (ESO
Recipe Flexible Execution Workbench v2.11.5; Freudling et al., 2013). We
followed the standard reduction procedure for the basic calibration, as well
as the sky background subtraction. The output of the MUSE data reduction
pipeline is a datacube (with RA, DEC, and wavelength as the three axes) in
FITS format. The final datacube has a size of 607 $\times$ 623 spatial pixels
(spaxels) = 121.4 $\times$ 124.6 arcsec, which corresponds to a physical size
of 0.59 $\times$ 0.60 pc, assuming an upper limit on the distance to the Vela
Jr. SNR of 1 kpc (Allen et al., 2015). The image quality of the final datacube
is 0.46 arcsec at 5000 Å, sampled with 0.2$\times$0.2 square arcsec spaxels.
For the purpose of our analysis in this work, we are focusing on the optical
nebula’s extent within a 10″$\times$ 10″ area centred on the location of the
Vela Jr. CCO.

We corrected the World Coordinate System (WCS) solution of the MUSE datacube
using the Gaia Data Release 3 (DR3) catalogue (Gaia Collaboration et al.,
2021). In order to isolate the emission lines from the underlying stellar or
nebular continuum in our data and to accurately measure the line properties,
we subtracted the continuum from our spectra. First, we fitted the continuum
within each spaxel in the datacube using the Locally Weighted Scatterplot
Smoothing algorithm (LOWESS; Cleveland, 1979; the application of LOWESS is
described in Vogt et al., 2017), and then performed the spaxel-by-spaxel
continuum subtractions. We utilized brutifus
2019.08.3333https://brutifus.github.io/index.html (Vogt, 2019), a Python
module designed to do post-processing on datacubes from integral fields
spectrographs, to execute the WCS correction and continuum subtraction.

##  3 Results and Characterization

The top panel of Figure 2 displays the white-light image of the central 2′
region of Vela Jr. SNR. This archival image was obtained with OmegaCAM at the
VLT Survey Telescope (VST Arnaboldi et al., 1998; Kuijken et al., 2002;
Kuijken, 2011) in 2015 as part of the VST Photometric H-Alpha Survey of the
Southern Galactic Plane and Bulge (VPHAS+ Drew et al., 2014), using the
narrow-band Hα filter NB_659. The Vela Jr. CCO position and the optical nebula
are labeled, and the white square indicates the 10″$\times$ 10″ cutout centred
on the CCO, corresponding to the bottom right image of Figure 2. Star Z is a
field star mentioned in Pavlov et al. (2001), and Wray 16-30 is a diffuse
structure initially designated as a Be star but later reclassified as a
planetary nebula (Wray, 1966; Reynoso et al., 2006) and is associated with the
star Ve 7-27 (this object was several times improperly referenced as Ve 2-27;
Gaia Collaboration et al., 2021).

The full field of our continuum-subtracted MUSE data as well as the
10″$\times$ 10″ cutout centred on Vela Jr. CCO are shown in the bottom row of
Figure 2, with colours representing the emission from Hα (cyan), [N ii]λ6583
(magenta), and [S ii]λ6716 (yellow). The red circle on the bottom right image
indicates the location of the CCO, with a 0.7″ radius corresponding to the
recomputed Chandra position uncertainty (Mignani et al., 2007). As can be seen
from the 10″ cutout, the optical nebula is unambiguously dominated by [N ii]
emission and is spatially coincident with the Vela Jr. CCO position. There is
a subtle asymmetry toward the south, where the brightness appears slightly
enhanced. Additionally, Hα emission is also evident in a distinct arc-like
shape, located approximately 3″ south of the CCO. The [N ii] and Hα structures
do not spatially align, with the [N ii] nebula is delimited by the Hα arc-like
structure. The [S ii] emission displays a notably different morphology
compared to both the [N ii] and Hα emissions, appearing as a circular blob
structure with a slight off-centre concentration relative to the [N ii]
emission.

###  3.1 Nebular morphology in the emission lines

Our MUSE datacube provides the very first spatially resolved view of the
optical nebula coincident with the location of Vela Jr. CCO. It is evident
that the nebula is predominantly composed of [N ii] emission. The [N ii]
nebula has a diameter of approximately 8″, which corresponds to about 0.04 pc,
assuming a distance of 1 kpc to the Vela Jr. SNR (Allen et al., 2015). In
addition to the bright emissions from Hα λ6563 and [S ii]λλ6716,6731, we also
detected emissions from [N i]λλ5198,5200, [N ii]λ5755, He i λ6678, λ7065444The
He i λ5876 emission line is not observable in our MUSE datacube due to the
blocking implemented around NaD lines in AO modes observations. This blocking
is necessary to prevent contamination from the sodium light used in the laser
guide star system., [Ar iii]λ7136, [Fe ii]λ8617, and [S iii]λ9069.

In Figure 3, we present the narrow-band emission line intensity maps (in erg
Å-1 cm-2 s-1) of the lines detected in the Vela Jr. CCO optical nebula,
extracted from the continuum-subtracted MUSE datacube. Each map represents a
10″$\times$ 10″ area, centred on the CCO and was generated by summing the
channel images over which the optical nebula structure is seen. The overlaid
contour represents the intensity levels of [N ii]λ6583 emission. The red
circle in each map marks the CCO position and the black dashdot line indicates
the estimated angular distance from the CCO to the peak of emission structure
(see column 5 in Table 1). These maps showcase the two distinct morphologies
of the optical nebula:

  1. 1.

An arc-like shape resembling a bow-shock nebula, observed in Hα.

  2. 2.

Smooth blob structures, observed in [N ii], [S ii], [Ar iii], [Fe ii], and [S
iii] emissions, and kidney-bean shaped blob seen in [N i] and He i emissions.

Figure 4: Pseudo-RGB images of different emission lines combinations, with
north is up, illustrating the spatial distribution of the emissions detected
in the Vela Jr. CCO optical nebula: Top left: [N ii]λ6583, [N ii]λ5755, and [N
i]λλ5198,5200 composite; Top right: [S iii]λ9069, [S ii]λ6716, and Hα
composite; Bottom left: He iλ6678, [N ii]λ6583, and Hα composite; Bottom
right: [S iii]λ9069, [Fe ii]λ8617, and [Ar iii]λ7136 composite. The red circle
in each image marks the position of Vela Jr. CCO.

We created pseudo-RGB images from the continuum-subtracted MUSE datacube to
illustrate the spatial distribution of different emission lines detected in
the optical nebula, shown in Figure 4.

  1. 1.

The first composite, [N ii]λλ6548,6583, [N ii]λ5755, and [N i]λλ5198,5200,
highlights the different ionization stratification of nitrogen lines and
reveals the elongated, kidney-bean-shaped blob of the [N i] emission tracing
the outer edge of [N ii]. The [N ii] emissions appear just slightly off-centre
from the CCO (Figure 4–top left).

  2. 2.

The second composite combining [S iii]λ9069, [S ii]λ6716, and Hα, shows that
both singly- and doubly-ionized sulfur emissions have diffuse circular blob
shapes that are near the apex of the Hα emission (Figure 4–top right).

  3. 3.

The third composite image, which combines He iλ6678 and Hα with respect to [N
ii]λ6583\. The He i emission appears to have a kidney-bean shaped blob and
traces the south outer edge of [N ii] emission (Figure 4–bottom left).

  4. 4.

Lastly, the fourth composite image, combining the emission lines of the
products of advanced stages of stellar evolution, [S iii]λ9069, [Fe ii]λ8617,
and [Ar iii]λ7136, shows similar circular shapes concentrated near the CCO,
offset slightly to the southwest (Figure 4–bottom right).

We noted that the MUSE datacube reveals a clear absence of structures from the
Hβ and [O iii]λλ4959,5007 emissions within the Vela Jr. CCO optical nebula.

###  3.2 Spectral analysis of the optical nebula

Figure 5: Hα and [N ii]λ6583 line profiles from 6555 Å to 6590 Å, obtained by
integrating over 5 $\times$ 5 spaxels boxes (1″$\times$ 1″) and overplotted on
Hα (cyan dotted-line) and [N ii] (magenta dotted-line) intensity contour maps.
Within each box, the blue and black lines respectively represent the spectra
from MUSE data and the Hα emission background-corrected spectra. The red
rectangle indicates the region from which the line profile spectra were used
for the correction. On the right, a magnified snippet shows the full extent of
the spectral lines strengths.

Figure 5 shows the intensity contour maps of the optical nebula morphology in
Hα (cyan) and [N ii] (magenta), superimposed with the spectral line profiles
from 6540 Å to 6590 Å, covering the Hα and [N ii]λ6583 lines (the blue and
black lines). These line profiles were generated by integrating the spectra
over 5 $\times$ 5 spaxels boxes (1″$\times$ 1″). The line profile of Hα
appears fairly constant in strength across the spaxel boxes, shown as the blue
lines in each box. This suggests the presence of some consistent background Hα
contribution across the field. To correct the Hα line profile for our
analysis, we selected an area outside the optical nebula where there is no
significant emission from the Hα nebula, as indicated by the red rectangle in
Figure 5. We then used the average line profile from this region for the
correction. The black lines represent the corrected Hα line profile, and the
snippet on the right provides a magnified view of the line profiles, showing
the full extent of the [N ii]λ6583 line strength.

We fitted the integrated spectrum of each emission line detected in the
optical nebula individually with a single-component Gaussian function, with
the exception for the [N ii]λλ6548,6583 and [S ii]λλ6716,6731 where each was
fitted jointly. We imposed the same velocity and line width for both emission
line doublets and constrained the flux ratio [N ii]λ6583/λ6548 to be equal to
the theoretical value of $\sim$3 (Osterbrock & De Robertis, 1985). The line
fitting was conducted on the integrated spectra using our own custom Python
routine. It starts by determining the spatial region for spectra integration,
which is chosen based on the intensity level of each emission line’s nebular
structure. A non-linear least square optimization method was then employed to
fit a single-component Gaussian to the integrated spectrum. The uncertainties
of the fitted parameters were measured as one standard deviation errors. The
resulting measurements, which include the total observed line flux, radial
velocity, and the corrected full width at half maximum (FWHM) velocity
dispersion, accounting for the instrument line broadening of MUSE, were
derived for all prominent lines and are presented in Table 1. Additionally, we
estimated the angular distance from the CCO location to the peak emission
location of each line emission through visual inspection (the last column of
Table 1 and indicated with the black dashdot lines in Figure 3).

The derived radial velocities of most of the emission lines indicate that the
nebular structures are moving toward us, with the exception for the Hα and [S
iii]. From the integrated spectra, we measured an unusually high [N
ii]6548$+$6583/Hα value of 33.8 for the optical nebula structure.

Using the measurements from the integrated spectra presented in Table 1, we
determined the electron temperature and electron density in the optical nebula
utilizing the [N ii] (λ6548+λ6583)/λ5755 and the [S ii]λ6716/λ6731 flux
ratios. We derived these values using the PyNeb code (Luridiana et al., 2015)
which solves the equilibrium equations for an n-level atom and works
iteratively with the ratio of the [N ii] lines, as well as of the [S ii] lines
to obtain more accurate values. The [N ii] lines provide an estimate for the
mean electron temperature of 7800 K for the optical nebula, and the [S ii]
lines give an electron density value of 35 cm-3.

Table 1: Observed total (integrated) line fluxes, radial velocities, derived full width at half maximum (FWHM) velocities that are corrected for the MUSE’s instrument line broadening, and the angular distances estimated visually from the Vela Jr. CCO location to the peak emission of each line emission. The uncertainties of the fitted parameters were measured as one standard deviation error. NOTE – We did not perform the fitting for [Fe ii]λ8617, and He i λ7065 due to their low integrated flux level. Line |  f${}_{\textnormal{\lambda}}$,obs | v | FWHM (corrected) | Angular Distance
---|---|---|---|---
| (10-16 erg cm-2 s-1) | (km s-1) | (km s-1) | (″)
[N i]λ5197 | 2.76 | -2.6 | 44.5 | 2.24
[N i]λ5200 | 3.39 | -2.6 | 44.3 | 2.24
[N ii]λ5755 | 1.50 | -12.1 | 61.9 | 1.02
[N ii]λ6548 | 66.36 | -9.0 | 54.6 | 1.26
Hα λ6563 | 8.17 | 3.3 | 81.7 | 3.26
[N ii]λ6583 | 209.84 | -9.0 | 54.3 | 1.26
He i λ6678 | 1.14 | -2.6 | 87.1 | 3.82
[S ii]λ6716 | 13.76 | -23.3 | 163.6 | 1.6
[S ii]λ6731 | 9.70 | -23.3 | 163.3 | 1.6
[Ar iii]λ7136 | 1.51 | -0.3 | 136.7 | 0.89
[S iii]λ9069 | 3.79 | 17.8 | 71.2 | 1.20

###  3.3 A CCO-ionized optical nebula

Figures 2–4 of the inner 4–6 arcseconds optical nebula show unusual emission,
including very strong [N ii] emission relative to Hα, as well as diffuse [S
ii] and [S iii]. The undetected (or weak) Hβ is also suggestive of a non-
standard Hα/Hβ ratio, which makes the common Balmer decrement reddening
estimation problematic. However, it is often the case in photoionized nebulae
modeling, such as fast shock precursors and other partially ionized plasmas
(Sutherland & Dopita, 2017; Dopita & Sutherland, 2017) for the Hα/Hβ ratio to
be intrinsically greater than the canonical 2.86 Case B (optically thick
nebula where the recombinations of H atom to the ground state generate photons
that are absorbed locally) ratio and Hα can be enhanced.

We performed a preliminary photoionization modeling with MAPPINGS V
v5.2.0555https://mappings.anu.edu.au (Sutherland et al., 2018) using a Wolf-
Rayet (WR) star progenitor – the late phase of massive star evolution
identified by high mass-loss rate and fast stellar winds. As a hypothesis, we
took a WN star (a subtype of WR stars, characterized by strong helium and
nitrogen lines) surface composition of a star with an initial mass of 100
M${\sun}$ from Roy et al. (2020) and evolved it past the late-type WN (WNL)
stage to the full He core burning WN of 50–60 M☉. By this late stage of
evolution, the atmosphere and surface composition are insensitive to the
initial abundances of the star, dominated by helium and nitrogen, and are
heavily CNO-processed. Under the assumption that these compositions from the
last 500,000 years of the stellar lifetime could constitute a late stage WN
wind environment for any hot CCO remnant, we constructed a simple spherical
isobaric X-ray nebula model using MAPPINGS V to test the general feasibility.

Using the Slane et al. (2001) estimated 0.47 keV blackbody temperature for the
Vela Jr. CCO’s thermal X-ray emission, a $5.45\times 10^{6}$ K blackbody as a
uniform 10 km radius source was placed in the WN stage with gas composition
listed in Table 2. The pressure of the nebula was iterated until it yielded
the densities that satisfy the observed [S ii]λ6731/λ6716 density ratio. The
result of our test model is summarized in Figure 6. The resulting nebula has a
$\sim 2.5\times 10^{17}$ cm radius in spherical symmetry and displays a high
[N ii]/Hα ratio, an extended [S ii], and no [S iii]λ6312 line (as the nebula
was largely too cool) in its spectra, along with a strong [N i] line, as shown
in the top panels of Figure 6. The inner $5\times 10^{16}$ cm of this nebula
is very hot and highly ionized, producing no observable optical nebular lines.
The nebula rapidly cools in an extended region spanning 0.5–2.5$\times
10^{17}$ cm, with temperatures from 3 $\times 10^{3}$ to 1.5$\times 10^{4}$ K,
and emits observable optical lines (Figure 6 bottom panels).

Table 2: Gas composition of 60 M☉ WN star progenitor for MAPPPINGS V test model | Number | Mass |  | Number | Mass
---|---|---|---|---|---
Element | $\log(n_{H})=0.0$ | Fraction | Element | $\log(n_{H})=0.0$ | Fraction
H | 0.00 | 2.26E-01 | Al | -4.90 | 7.59E-05
He | -0.07 | 7.64E-01 | Si | -4.12 | 4.81E-04
C | -4.58 | 7.03E-05 | Si | -4.50 | 2.30E-04
N | -2.76 | 5.53E-03 | Cl | -5.25 | 4.49E-05
O | -5.06 | 3.14E-05 | Ar | -4.35 | 4.02E-04
Ne | -3.74 | 8.32E-04 | Ca | -5.25 | 5.03E-05
Na | -4.46 | 1.77E-04 | Fe | -4.29 | 6.39E-04
Mg | -4.14 | 3.95E-04 | Ni | -4.65 | 2.92E-04
Figure 6: MAPPINGS V test model of a nebula with 60 M☉ WN progenitor. Top: The
simulated spectra of the nebula. Bottom: The intrinsic variation of electron
temperatures for a spherical model of the nebula as a function of the nebular
radius and the spherical MAPPINGS V model of the ionization structures within
the nebula.

This preliminary result is circumstantial, and more detailed modeling that
includes partial CCO surface heat and non-thermal components will be presented
in future work. Nevertheless, it supports the existence of a CCO, and suggests
that at least some of the nebulosity is attributed to an X-ray blackbody
object, in addition to the widespread shock processes. The general ionization
structure of the nebula from the test model broadly agrees with the
observations, specifically in reproducing the strong widespread of nitrogen
and the extended sulfur emission in particular, which distinguished this from
a typical H II region. With further exploration and modeling, this hypothesis
may resolve if the inner 4″is a neutron-star-photoionized nebula or bubble,
and the bright Hα arc to the south is an interaction or bow shock structure.
If it is true that the CCO is a $\sim 5.5$ MK, with about 10 km radius source,
and able to photoionize a WN composition circum–progenitor environment this
lends support to the Allen et al. (2015) suggestion that Vela Jr has a massive
progenitor.

##  4 Discussion

###  4.1 Comparison of MUSE observations with previous Hα observations

We compiled previously published Hα observations of the optical nebula from
Pellizzoni et al. (2002), Mignani et al. (2007), Mignani et al. (2009), and
Mignani et al. (2019), alongside the MUSE images in Hα and [N ii] from the
non-continuum-subtracted datacube in Figure 1. Each frame colour corresponds
to the line colour of the respective filter bandpass range in Figure 7. Since
its discovery in 2002, the optical nebula spatially coincident at the position
of the CCO in Vela Jr. SNR was initially believed to be Hα emission, with a
kidney-bean-like shape. MUSE data unequivocally reveal that the observed shape
of the optical nebula is primarily comprised of [N ii] emission rather than
Hα, and the Hα emission exhibits a distinct arc-like morphology. This could be
attributed to the fact that we observe a very high [N ii]/Hα intensity ratio
for the optical nebula and also that the Hα filters used in the most of the
previous observations also cover the [N ii]λλ6548,6583 lines within their
bandpass response ranges (see Figure 7 for details).

Figure 7: The filter bandpass response range of all the filters used in the
previous optical observations of the Vela Jr. CCO optical nebula is
represented by coloured lines, that match the frame colour of the images in
Figure 1. The MUSE wavelength ranges from which each MUSE image is integrated
are shown as yellow shaded areas. Left: The filter bandpasses spanning from
5400 Å to 8400 Å, to show the extended range of the R Bessel filter used for
the VLT/FORS1 observations. Right: The filter bandpasses range from 6520 Å to
6650 Å, to emphasize the narrow-band filters. Star symbols indicate the
central wavelength of each filter

###  4.2 The nature of the optical nebula

The optical nebula observed at the location of the CCO in Vela Jr. SNR has
sparked debate regarding its nature since its discovery over 20 years ago.
Extensive radio and X-ray studies have not detected any PWN associated with
the CCO (e.g. Pavlov et al., 2001; Reynoso et al., 2006; Maxted et al., 2018).
PWNe are expected typically to be powered by the rotational energy loss of an
underlying pulsar or neutron star and, when young and still inside their SNR,
they are confined by the pressure of the surrounding SNR material. They are
also typically discovered in radio and X-rays. The lack of a PWN around the
Vela Jr. CCO is consistent with the lack of PWNe in the CCO class.

Our MUSE observations have revealed that the morphology of the optical nebula,
particularly in Hα, is reminiscent of a bow shock nebula, which is typically
observed when pulsars move supersonically in their surrounding medium creating
a nebula that is confined by the ram pressure in the ISM (see e.g. Gaensler &
Slane, 2006). The bow-shock-like structure is at approximately 3″away from the
neutron star, which at a distance of 1 kpc, corresponds to about 0.015 pc.
This distance is comparable to termination shocks seen in other PWNe. Such
nebulae are characterized by being non-radiative and dominated by Hα Balmer
line, with no detectable lines from any other element in the optical band.
Although Mignani et al. (2019) dismissed the bow shock nebula scenario based
on HST observations, our analysis using MUSE data suggests that this scenario
cannot be excluded. Bow shock nebulae provide a valuable tool to probe the
shock physics and can be used to constrain the underlying pulsar (or neutron
star)’s spin-down energy $\dot{E}$ calculation. From X-rays observations of
the Vela Jr. SNR, an upper limit on the ambient density for thermally emitting
plasma with $kT$$\sim$0.5 keV is established as $n_{0}<2.9\times 10^{-2}$
cm-3(Slane et al., 2001). With the electron density set as 4.8$n_{0}$, we
estimate the thermal pressure $2n_{e}kT$ to be $\sim$$2.096\times 10^{-10}$
erg cm-3. Considering the standoff distance from the CCO to the Hα bow shock
nebula determined through MUSE data to be approximately 0.015 pc, and assuming
a spherical shape, as well as that the ram pressure of the wind termination
shock is balanced by the surrounding thermal pressure, the estimated upper
limit for any hidden pulsar’s $\dot{E}$ is $1.69\times 10^{35}$ erg s-1.

Mignani et al. (2019) set an upper limit of 300 mas yr-1 for the proper motion
of the Vela Jr. CCO, based on a comparison of Chandra observations taken 9
years apart. However, a recent study of SRG/eROSITA and XMM-Newton data of the
remnant has obtained a new coordinate for the geometric centre of Vela Jr.,
showing that the current position of the CCO is consistent with zero proper
motion (Camilloni et al., 2023). This suggests that the CCO has likely
remained in its birthplace without significant displacement.

The morphology of the optical nebula observed in the other emission lines, as
well as the overabundance of nitrogen, indicates a WR circumstellar nebula,
which consist of material ejected by a massive star during its mass loss
phase, that is currently being ionized by the Vela Jr. CCO. The circumstellar
nebulae commonly found in young SNRs are typically characterized by slow
expansion velocities and elevated [N ii]/Hα ratios (e.g. Hartigan et al.,
1987; Chu, 2001). It has been proposed that Vela Jr. CCO is a product of the
collapse of a massive star (Slane et al., 2001; Allen et al., 2015), and the
recent eROSITA observation of the SNR supports a massive progenitor (Camilloni
et al., 2023), although no specific mass or mass range had been inferred in
these studies. Our test model, with a 60 M☉ WN star as the progenitor using
MAPPINGS V, has provided circumstantial evidence that an object with a
blackbody temperature similar to Vela Jr. CCO can be sufficient to create an
optical nebula that is N-rich, H-weak, O-poor, and still shows faint emission
from ionized heavy elements, much like what we observe. Hence, it is
reasonable to suggest that the progenitor is a WN star, that generated a
nitrogen-rich circumstellar nebula during its mass loss. This nitrogen-rich
nebula would have been subsequently shocked by the SN explosion and it is then
being ionized by the CCO, leading to the nitrogen-abundant optical nebula we
observe today. To validate our hypothesis, further modeling is needed, and we
intend to address this in our forthcoming paper. If confirmed, this would be
the first detection of a WR nebula around a CCO.

##  5 Summary and Future Work

In this article, we have presented and discussed the results of a deep optical
integral field spectroscopy observation, obtained with VLT’s MUSE, of the
optical nebula spatially coincident with the CCO in the Vela Jr. SNR. The sub-
arcsecond spatial resolution of MUSE, combined with its AO capabilities,
enables us to spatially-resolve, for the first time, the emission distribution
of the long-debated optical nebula and characterize its emission line spectrum
as being dominated by [N ii] emission. The MUSE observation enabled the
following key findings:

  1. 1.

The optical nebula exhibits a significant overabundance of [N ii]λλ6548,6583
emission, approximately 33 times the intensity of Hα emission.

  2. 2.

We identified a complex set of spectral lines associated to the optical
nebula, including the nitrogen auroral line, sulfur lines, and the emission
lines of the products of advanced stages of stellar evolution.

  3. 3.

The nebula displays two distinct morphologies: an arc-like shape that
resembles a bow-shock nebula (observed in Hα) and smooth blob structures
(observed in [N ii], [S ii], [Ar iii], [Fe ii], and [S iii]; and some with
kidney-bean shape, observed on [N i], and He i emissions), while there is no
nebulosity structure detected in the Hβ and [O iii] lines.

  4. 4.

We observed spatial structuring within the optical nebula, including a clear
stratification of the ionized nitrogen levels.

The MUSE clear detection of the distinct morphology of the optical nebula in
Hα, which was previously unobserved, has revived the bow shock nebula
interpretation. The bow shock nebula scenario was dismissed based on the non-
detection of the nebula in the HST WFPC2 F656N data (Mignani et al., 2009) and
the very minor contribution from Hα in the nebula’s spectrum (Mignani et al.,
2019). The overabundance of nitrogen in the optical nebula is consistent with
a WR star progenitor and with the notion of the Vela Jr. SNR being associated
with a massive progenitor. Our preliminary test model using MAPPINGS V (i) has
demonstrated the feasibility of a massive star in its WR stage producing an
optical nebula akin to the one observed in Vela Jr. CCO and (ii) supports the
suggestion that very massive progenitors not necessarily making black holes.

To gain a deeper understanding of the physical conditions within the nebula, a
comparison of our observations with models is needed, particularly to simulate
the heating and cooling processes that lead to the observed emissions.
Photoionization models such as MAPPINGS and CLOUDY (Ferland et al., 2013)
solve the ionization and energy equilibrium equations, calculate the radiation
transfer, and provide estimates of the physical quantities and the line
intensities expected from a specific gaseous nebula. Through modeling, we also
seek to comprehend the origins of the concentrated emissions of sulfur, argon,
and iron within the nebula, and how they fit to the optical nebula
interpretation. A more detailed modeling of the optical nebula’s spectra using
photoionization codes is deferred to a future publication.

## Acknowledgements

This work made use of NASA’s Astrophysics Data System and various Python
packages: Numpy (Harris et al., 2020), Scipy (Virtanen et al., 2020), Astropy
(Astropy Collaboration et al., 2013, 2018, 2022), brutifus (Vogt, 2019) (a
Python module to process datacubes from integral field spectrographs, that
relies on statsmodel (Seabold & Perktold, 2010), matplotlib, astropy, and
photutils, an affiliated package of astropy for photometry), Matplotlib
(Hunter, 2007) and specifically cubehelix colormap (Green, 2011), and APLpy
(Robitaille & Bressert, 2012). JS acknowledges funding from the University of
Manitoba Graduate Fellowhips (UMGF) and the Ernst and Ingrid Bock Graduate
Award. SSH acknowledges support from the Natural Sciences and Engineering
Research Council of Canada (NSERC) through the Discovery Grants and the Canada
Research Chairs programs, and by the Canadian Space Agency. CJL acknowledges
the support of the NSTC grants 110-2112-M-001-020 and 111-2112-M-001-063. AJR
was supported by the Australian Research Council through award number
FT170100243.

## Data Availability

The Multi Unit Spectroscopic Explorer (MUSE) data can be downloaded from the
ESO Science Archive Facility at http://archive.eso.org/cms.html, under
programme ID 0104.D-0092(B), P.I.: F.P.A. Vogt.

## References

  * Aharonian et al. (2005) Aharonian F., et al., 2005, A&A, 437, L7
  * Aharonian et al. (2007) Aharonian F., et al., 2007, ApJ, 661, 236
  * Allen et al. (2015) Allen G. E., Chow K., DeLaney T., Filipović M. D., Houck J. C., Pannuti T. G., Stage M. D., 2015, ApJ, 798, 82
  * Arnaboldi et al. (1998) Arnaboldi M., Capaccioli M., Mancini D., Rafanelli P., Scaramella R., Sedmak G., Vettolani G. P., 1998, The Messenger, 93, 30
  * Arsenault et al. (2008) Arsenault R., et al., 2008, in Hubin N., Max C. E., Wizinowich P. L., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7015, Adaptive Optics Systems. p. 701524, doi:10.1117/12.790359
  * Aschenbach (1998) Aschenbach B., 1998, Nature, 396, 141
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, ApJ, 935, 167
  * Bacon et al. (2010) Bacon R., et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 773508, doi:10.1117/12.856027
  * Banovetz et al. (2021) Banovetz J., et al., 2021, ApJ, 912, 33
  * Braun et al. (2019) Braun C., Safi-Harb S., Fryer C. L., 2019, MNRAS, 489, 4444
  * Braun et al. (2023) Braun C., Safi-Harb S., Fryer C., Zhou P., 2023, arXiv e-prints, p. arXiv:2308.13693
  * Camilloni et al. (2023) Camilloni F., Becker W., Predehl P., Dennerl K., Freyberg M., Mayer M. G. F., Sasaki M., 2023, A&A, 673, A45
  * Chu (2001) Chu Y.-H., 2001, in Holt S. S., Hwang U., eds, American Institute of Physics Conference Series Vol. 565, Young Supernova Remnants. pp 409–418 (arXiv:astro-ph/0012520), doi:10.1063/1.1377129
  * Cleveland (1979) Cleveland W. S., 1979, Journal of the American Statistical Association, 74, 829
  * De Luca (2017) De Luca A., 2017, in Journal of Physics Conference Series. p. 012006 (arXiv:1711.07210), doi:10.1088/1742-6596/932/1/012006
  * Dopita & Sutherland (2017) Dopita M. A., Sutherland R. S., 2017, ApJS, 229, 35
  * Drew et al. (2014) Drew J. E., et al., 2014, MNRAS, 440, 2036
  * Ferland et al. (2013) Ferland G. J., et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137
  * Ferrand & Safi-Harb (2012) Ferrand G., Safi-Harb S., 2012, Advances in Space Research, 49, 1313
  * Freudling et al. (2013) Freudling W., Romaniello M., Bramich D. M., Ballester P., Forchi V., García-Dabló C. E., Moehler S., Neeser M. J., 2013, A&A, 559, A96
  * Gaensler & Slane (2006) Gaensler B. M., Slane P. O., 2006, ARA&A, 44, 17
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Gotthelf et al. (2013) Gotthelf E. V., Halpern J. P., Allen B., Knispel B., 2013, ApJ, 773, 141
  * Green (2011) Green D. A., 2011, Bulletin of the Astronomical Society of India, 39, 289
  * H. E. S. S. Collaboration et al. (2018) H. E. S. S. Collaboration et al., 2018, A&A, 612, A7
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hartigan et al. (1987) Hartigan P., Raymond J., Hartmann L., 1987, ApJ, 316, 323
  * Hebbar et al. (2020) Hebbar P. R., Heinke C. O., Ho W. C. G., 2020, MNRAS, 491, 1585
  * Ho (2011) Ho W. C. G., 2011, MNRAS, 414, 2567
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Iyudin et al. (1998) Iyudin A. F., et al., 1998, Nature, 396, 142
  * Kuijken (2011) Kuijken K., 2011, The Messenger, 146, 8
  * Kuijken et al. (2002) Kuijken K., et al., 2002, The Messenger, 110, 15
  * Li et al. (2021) Li C.-J., Seitenzahl I. R., Ishioka R., Chu Y.-H., Ruiter A. J., Vogt F. P. A., 2021, ApJ, 915, 20
  * Long et al. (2020) Long X., Gaetz T. J., Plucinsky P. P., 2020, ApJ, 904, 70
  * Luridiana et al. (2015) Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573, A42
  * Maxted et al. (2018) Maxted N. I., et al., 2018, ApJ, 866, 76
  * Mereghetti (2001) Mereghetti S., 2001, ApJ, 548, L213
  * Mignani et al. (2007) Mignani R. P., de Luca A., Zaggia S., Sester D., Pellizzoni A., Mereghetti S., Caraveo P. A., 2007, A&A, 473, 883
  * Mignani et al. (2009) Mignani R. P., de Luca A., Pellizzoni A., 2009, A&A, 508, 779
  * Mignani et al. (2019) Mignani R. P., De Luca A., Zharikov S., Hummel W., Becker W., Pellizzoni A., 2019, MNRAS, 486, 5716
  * Osterbrock & De Robertis (1985) Osterbrock D. E., De Robertis M. M., 1985, PASP, 97, 1129
  * Parker & Phillipps (1998) Parker Q. A., Phillipps S., 1998, Publ. Astron. Soc. Australia, 15, 28
  * Pavlov et al. (2001) Pavlov G. G., Sanwal D., Kızıltan B., Garmire G. P., 2001, ApJ, 559, L131
  * Pavlov et al. (2004) Pavlov G. G., Sanwal D., Teter M. A., 2004, in Camilo F., Gaensler B. M., eds, Vol. 218, Young Neutron Stars and Their Environments. p. 239 (arXiv:astro-ph/0311526), doi:10.48550/arXiv.astro-ph/0311526
  * Pellizzoni et al. (2002) Pellizzoni A., Mereghetti S., De Luca A., 2002, A&A, 393, L65
  * Renaud et al. (2006) Renaud M., Vink J., Decourchelle A., Lebrun F., Terrier R., Ballet J., 2006, New Astron. Rev., 50, 540
  * Reynoso et al. (2006) Reynoso E. M., Dubner G., Giacani E., Johnston S., Green A. J., 2006, A&A, 449, 243
  * Robitaille & Bressert (2012) Robitaille T., Bressert E., 2012, APLpy: Astronomical Plotting Library in Python, Astrophysics Source Code Library, record ascl:1208.017 (ascl:1208.017)
  * Rogers & Safi-Harb (2016) Rogers A., Safi-Harb S., 2016, MNRAS, 457, 1180
  * Roy et al. (2020) Roy A., Sutherland R. S., Krumholz M. R., Heger A., Dopita M. A., 2020, MNRAS, 494, 3861
  * Safi-Harb (2017) Safi-Harb S., 2017, in Journal of Physics Conference Series. p. 012005 (arXiv:1712.06040), doi:10.1088/1742-6596/932/1/012005
  * Seabold & Perktold (2010) Seabold S., Perktold J., 2010, in 9th Python in Science Conference.
  * Slane et al. (2001) Slane P., Hughes J. P., Edgar R. J., Plucinsky P. P., Miyata E., Tsunemi H., Aschenbach B., 2001, ApJ, 548, 814
  * Ströbele et al. (2012) Ströbele S., et al., 2012, in Ellerbroek B. L., Marchetti E., Véran J.-P., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 8447, Adaptive Optics Systems III. p. 844737, doi:10.1117/12.926110
  * Sutherland & Dopita (2017) Sutherland R. S., Dopita M. A., 2017, ApJS, 229, 34
  * Sutherland et al. (2018) Sutherland R., Dopita M., Binette L., Groves B., 2018, MAPPINGS V: Astrophysical plasma modeling code, Astrophysics Source Code Library, record ascl:1807.005 (ascl:1807.005)
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogt (2019) Vogt F. P. A., 2019, brutifus: Python module to post-process datacubes from integral field spectrographs, Astrophysics Source Code Library, record ascl:1903.004 (ascl:1903.004)
  * Vogt et al. (2017) Vogt F. P. A., Seitenzahl I. R., Dopita M. A., Ghavamian P., 2017, A&A, 602, L4
  * Vogt et al. (2018) Vogt F. P. A., Bartlett E. S., Seitenzahl I. R., Dopita M. A., Ghavamian P., Ruiter A. J., Terry J. P., 2018, Nature Astronomy, 2, 465
  * Weilbacher et al. (2016) Weilbacher P. M., Streicher O., Palsa R., 2016, MUSE-DRP: MUSE Data Reduction Pipeline (ascl:1610.004)
  * Weilbacher et al. (2020) Weilbacher P. M., et al., 2020, A&A, 641, A28
  * Wray (1966) Wray J. D., 1966, PhD thesis, Northwestern University