© The Authors 2023

1 Introduction

The increase in angular resolution and sensitivity of γ-ray detectors like the High Energy Stereoscopic System (HESS; H.E.S.S. Collaboration 2018a) produced a considerable step forward in multiwavelength studies of supernova remnants (SNRs) in the last decade. The main mechanism responsible for particle acceleration in SNRs is diffusive shock acceleration (Blandford & Eichler 1987). This is particularly true for young remnants, which are often characterized by a featureless continuum spectrum from the radio to the X-ray and γ-ray bands (e.g., Bamba et al. 2005b; Miceli et al. 2016; Okuno et al. 2018), indicating that the detected radiation is nonthermal synchrotron emission. Moreover, multiband fits like those presented by Pannuti et al. (2010) and Kishishita et al. (2013) have allowed us to put constraints on the underlying electron distribution, relevant to characterizing the Galactic cosmic ray population. By the middle of 2022, the online catalog of TeV sources at the University of Chicago (Wakely & Horan 2008) lists 95 entries.¹ Among them, only eight are supernova remnants, of which Vela Jr (also known as RX J0852.0-4622 or G266.2-1.2) is the one with the largest extent. All these remnants are shell-type (H.E.S.S. Collaboration et al. 2018d,c), excluding contributions of pulsars or pulsar wind nebulae (PWN). Vela Jr was discovered by Aschenbach (1998) using data from the ROSAT all-sky survey, which was the first all-sky survey performed by an imaging X-ray telescope. With the launch of the Spectrum-Röntgen-Gamma (SRG) mission in June 2019 (Sunyaev et al. 2021), for which the extended ROentgen Survey with an Imaging Telescope Array (eROSITA) is the soft X-ray instrument on board (Predehl et al. 2021), another possibility for observing the two-degree wide Vela Jr remnant in its entirety by a spectro-imaging X-ray telescope became available. In this work we report on the results of an analysis of this data. In order to increase the photon statistics for our analysis, we also made use of archival XMM-Newton calibration data taken on the northwestern (NW) region of the remnant between 2001 and 2021.

Soon after the discovery of Vela Jr, Aschenbach et al. (1999) found evidence in favor of a remnant distance of less than 1 kpc and an estimated remnant age of 680 yr. This study relied mostly on data from the COMPTEL γ-ray observatory and an apparent ⁴⁴Ti signal from Vela Jr (Iyudin et al. 1998; Chen & Gehrels 1999). The isotope ⁴⁴Ti has a half-life of only 60 yr, so that its presence lends strong support for a young remnant age. However, the significance of that ⁴⁴Ti emission line was at the 3σ level only, (Renaud et al. 2006) and later observations indeed did not confirm its existence (Hiraga et al. 2009; Weinberger et al. 2020). Today, more than 30 yr after the discovery of Vela Jr in the ROSAT all-sky survey, there is still no agreement about its distance and age. One of the main observational problems with this remnant is its location along the line of sight of the much more extended Vela SNR. The presence of Vela SNR makes the spectral analysis of Vela Jr very complex. Especially in the soft X-ray band, spectral modeling is limited by the problem of disentangling the emission from the two remnants. Most of the spectral studies in X-rays conducted so far have ignored the energy range below 1 keV for that specific reason (Iyudin et al. 2005; Katsuda et al. 2008; Pannuti et al. 2010; Takeda et al. 2016).

Using ASCA data, Slane et al. (2001) estimated a remnant distance of 1–2 kpc based on the significantly higher column density measured for Vela Jr than for the Vela SNR. If correct, Vela Jr would be a background object relative to the Vela SNR, which a parallax measurement of the central pulsar placed at a distance of ${287}_{-17}^{+19}$ pc (Dodson et al. 2003). Conversely, from the estimate given in Aschenbach et al. (1999), it is possible that Vela Jr is a foreground object. Several later studies concluded that Vela Jr was much closer than 1–2 kpc. For example, Bamba et al. (2005a) studied the filamentary structure of the NW rim with Chandra and obtained a distance of around 0.33 kpc, which is almost compatible with the distance of the Vela SNR. Katsuda et al. (2008) followed a different approach by studying the expansion rate of the NW rim in order to test the scenarios presented by Aschenbach et al. (1999) and Slane et al. (2001): the rather low expansion rate of 0.84 ± 0.23″ yr⁻¹ made them conclude that Vela Jr was most likely a supernova remnant of age 1700–4300 yr, located at a distance of around 750 pc. However, also these authors were not able to completely rule out the hypothesis that Vela Jr is indeed a very young remnant that was rapidly decelerated, for example in the interaction with a dense interstellar medium.

At X-ray energies, the northwestern region of Vela Jr was found to show a filamentary trailing emission behind a shock front, tentatively associated with a reverse shock (Iyudin et al. 2005). This region also appears to be the brightest part of the remnant. To explore it in more detail, it was in the focus of several XMM-Newton and Chandra observations. XMM-Newton in particular monitored the NW region in regular observations between 2001 and 2021. Allen et al. (2015) measured the expansion rate of the NW rim, comparing two Chandra observation performed in 2003 and 2008. They obtained an expansion rate of 0.42 ± 0.10 arcsec yr⁻¹, which is half that measured using XMM-Newton by Katsuda et al. (2008). Inserting this expansion rate as a constraint on a numerical simulation based on the models of Truelove & McKee (1999), the authors concluded that Vela Jr has an age in the range 2.4–5.1 kyr. Using synchrotron cutoff frequency measurements from Pannuti et al. (2010), these authors set a lower limit of 0.5 kpc to the distance of Vela Jr.

As described above, there is no clear tendency for the distance estimate and age of Vela Jr. Summarizing the results from the literature sets the lower and upper limits for the remnant distance at about 0.3 and 2 kpc. while there is general consent that the remnant is young, whether it is a historical remnant of age ~700 yr or whether its age is more comparable with that of Puppis A (around 4000 yr; see Winkler et al. 1988; Mayer et al. 2020) is still under debate.

In this work, we use eROSITA all-sky survey data to conduct a detailed spectro-imaging analysis of the X-ray emission of Vela Jr. Our data set constitutes the most sensitive observation covering the entire remnant in the 0.2–10.0 keV energy band, with an angular resolution of 26 arcsec in survey mode and 15 arcsec in pointed mode (NW rim). It further offers a relatively uniform exposure over the extent of Vela Jr, eliminating the need to create mosaics from many individual observations. In addition, to better constrain the emission from the NW region of Vela Jr, we included archival XMM-Newton data in our analysis. These data were taken for calibration purposes² on a yearly basis between 2001 and 2021.

Our paper is organized as follows. Section 2 presents the basic characteristics of our data sets and describes initial data reduction steps taken to ensure its correct treatment. We describe our methods and results in imaging and spectroscopic analysis in Sect. 3 and Sect. 4; the core results of spatially resolved spectroscopy of Vela Jr are presented in Sect. 4. In Sect. 5, we present our distance analysis. Finally, in Sect. 6 we summarize our results and discuss their physical implications.

2 Data reduction

2.1 eROSITA data reduction

Vela Jr was in the eROSITA field of view during four consecutive all-sky surveys performed between December 2019 and December 2021. The unvignetted exposure time of all survey observations covering Vela Jr sum up to about 1270 s. The vignetting-corrected exposure averaged over the size of the remnant for the full 0.2–8.0 keV band is 658 s. In addition to the survey observations, eROSITA observed the NW rim of Vela Jr during its CalPV phase for an exposure time of 60 ks. This observation was focused on the same position as the XMM-Newton calibration observations taken between 2001 and 2019. This CalPV eROSITA observation was originally scheduled for the purpose of inter-calibrating its detector response with that of XMM-Newton, and was performed on October 29–30, 2019. A procedure very similar to that described in Sect. 2.2 was followed for the data reduction of the eROSITA survey and CalPV data. We used the eROSITA science analysis software (eSASS) version 211214³ (Brunner et al. 2022) with the calibration files CALDB 2021Q4 for our data analysis.

eROSITA consists of seven telescope modules (TMs), each equipped with its own detector in the focal plane. A detailed description of the instrument is given in Predehl et al. (2021). Due to a light leak from the backside of the focal plane, TM5 and TM7 can be less suitable for spectroscopic analysis for certain Sun angle constraints, making the calibration of these modules less solid (at that time of the mission). This problem primarily affects the low energy sensitivity (Predehl et al. 2021). For the spectral analysis of the Vela Jr data we therefore employed only the merged data from TMs 1, 2, 3, 4, and 6. The cleaned photon event files were extracted using the eSASS command evtool to filter the good time intervals, and selecting all the available detection patterns (PATTERN = 15). We then extracted spectra, background, redistribution matrix file (RMF), and ancillary response file (ARF) using the srctool tasks.

2.2 XMM-Newton data reduction

The XMM-Newton data used in our analysis of the NW region are summarized in Table 1. These data were reduced using the XMM Science Analysis System (SAS version 19.1.0) together with the latest calibration files. We first created filtered event files. To do so, we generated one light curve per observation to look for the presence of solar flaring events, and cleaned the data by accepting only the good time intervals when the sky background was low. The same process was applied for both EPIC-PN and EPIC-MOS data, whose spectra were employed in the final spectral analysis, which we restricted to the 0.3–10 keV energy band. After this cleaning procedure, we extracted the images and used them to spatially divide the NW shock region into seven sectors. The regions were tailored in accordance with the regions selected in the work by Pannuti et al. (2010). We finally produced spectra for each of the seven sectors with the highest possible statistics, selecting events with PATTERN<=4,12 for PN and MOS and FLAG==0 settings in the filtering expression option of the evselect command. For each region a source spectrum, background, RMF and ARF were extracted using the commands evselect, rmfgen, and arfgen, respectively. Since we modeled the background (see Appendix B), we did not apply any grouping to the spectra. The same procedure was applied for all observations listed in Table 1, which sum up to a total exposure time of more than 410 ks.

3 Spatial analysis

In survey mode eROSITA provides an almost unlimited field of view, and thus allowed us to map the whole Vela Jr remnant with an XMM-Newton-like sensitivity. Figure 1 depicts a multispectral image of Vela Jr as obtained from stacked eROSITA data of all four sky surveys (eRASS:4) taken between December 2019 and December 2021. Impressively, the RGB color-coding reveals for the first time the different emission contributions of the Vela SNR (in red and green) and of Vela Jr (in blue) over its entire extent.

Having data of that high quality of the entire remnant allowed us for the first time after the ROSAT all-sky survey to address the question of the remnant’s geometrical center and to compare it with the position of the central compact object (CCO) CXOU J085201.4-461753 (cf. Becker et al. 2006, and references therein). The first measurement of the geometrical center of Vela Jr was obtained by Aschenbach (1998). Our measurement is based on the eRASS:4 intensity image shown in Fig. 2. The image was created from all events recorded within the 1.1–8.0 keV energy band using all seven telescope modules. In this energy band Vela Jr stands out most clearly above the Vela SNR emission. To demonstrate this, we applied an energy-sliced image analysis within the 0.7–1.3 keV band, comparing the emission contributions of the Vela SNR and Vela Jr in steps of 0.1 keV. The result is depicted in Fig. 3. In order to enhance the visibility of diffuse emission in these images while leaving point sources unsmoothed to the greatest possible extent, we applied the adaptive smoothing algorithm of Ebeling et al. (2006) with a Gaussian kernel allowed to vary between σ_min = 3.0 pixels and σ_max = 4.0 pixels. As can be seen in Fig. 2, the north to northeastern part of the remnant appears almost perfectly circular. Assuming that only this part of the remnant tracks the position of the geometric remnant center, we fitted the upper half only with a circle so that it follows closely the outer remnant boundary. In the first instance we did this fitting by eye. To overcome possible subject bias, these fits were performed independently by several members of the team. In Sect. 6, we discuss the choice of using this part of the remnant alone to determine the geometric center. From all measurements, we computed the mean and the 1σ standard deviation and found the center of the circle at RA = 08h 51m 58s(±1s), Dec = −46° 17′52″(±22″). The radius of the fitted circle is 0.95° ± 0.02°. The best-fitting circle is indicated in Fig. 2. The two lines indicated in Fig. 2 connect the four corners of the square exactly surrounding that best-fitting circle. As can be seen, the lines cross exactly at the location of the CCO. In eROSITA we find the position of the CCO at RA = 08h52m01.42s, Dec = −46° 17′ 54.52″, which is 1.26″ different from the position RA = 08h52m01.4s, Dec = −46° 17′53.3″ obtained with a typical error of 0.6″ in a Chandra observation by Pavlov et al. (2001). Based on this CCO position, the separation between the center of the fitted circle and the CCO is only 35.2 ± 15.8″, which is within the 2σ uncertainty range in agreement with zero proper motion. This result was supported by applying an alternative technique based on the edge detection filter EDGE_DOG⁴ to Fig. 2 by using the Interactive Data Analysis Language IDL. This IDL function determined the outer boundary of Vela Jr on which we overlaid a circle region in ds9. The radius and geometrical center of this circle turned out to be in full agreement with the results described above. Although the very southern part of the remnant is not well described by the circle that fits the northern part of the remnant, the fact that the position of the CCO is in agreement with the center of the best-fitting circle makes it appealing to speculate that it is still at its birthplace and has not not moved significantly from there since its formation.

Mignani et al. (2019) measured an upper limit of 10 mas yr⁻¹ for the proper motion of an infrared candidate counterpart for the CCO in the center of Vela Jr. Employing the geometric center value given by ROSAT and the CCO position obtained with Chandra (Mignani et al. 2019), the separation is close to 4 arcmin, while, conversely, it should be comprised between 24″ and 51″ with an age of 2.4–5.1 kyr and proper motion of 10 mas yr⁻¹. From our measurements, if only the NW rim is considered when fitting a circle to the rim, we find that the geometric center of Vela Jr is very close to the CCO position, with a separation of (35.2 ± 15.8)″, so this is in accordance with the proper motion limit set by the infrared counterpart. Using the transverse velocity relation employed in Mignani et al. (2019) $$V_{T,100}~5\times {\mu }_1\times d_{100},$$(1)

where d₁₀₀ is the distance in units of 100 pc, V_T,100 the velocity in units of 100 km s⁻¹ and µ₁ the proper motion in units of 1 arcsec yr⁻¹, it is possible to infer information on the velocity of the CCO. The (35.2 ± 15.8)″ of separation measured in this paper implies an upper limit on the proper motion of (14.7 ± 6.6) mas yr⁻¹ , if an age of 2400 yr is assumed for the remnant. Assuming a distance of 750 pc (Katsuda et al. 2008), we infer a transverse velocity υ = (55 ± 25) km s⁻¹. If instead an age of 5100 yr is assumed (upper limit on the age range given by Allen et al. 2015), the proper motion is around (6.9 ± 3.0) mas yr⁻¹, which translates into a velocity of υ = (26 ± 12) km s⁻¹ at the same distance. In both cases the value is fully in accordance with the upper limit on the proper motion set by Mignani et al. (2019) for the possible IR counterpart. Nevertheless, we note the number of uncertainties present in our analysis: the broad range of values for the age of the remnant, the uncertainty in our determination of the geometric center, and, for the velocity, the distance.

Fig. 1 False color image of Vela Jr. The energy of the X-ray photons are color-coded: red (0.2–0.7 keV), green (0.7–1.2 keV), and blue (1.2–8.0 keV). In order to enhance the visibility of the diffuse X-ray emission, the images of the different bands were convolved with a Gaussian filter (expressed in pixels; red: radius = 7, σ = 3.5; green: radius = 8, σ = 4; blue: radius = 8, σ = 4).

Fig. 2 Intensity representation of the emission of Vela Jr and its CCO J085201.4-461753 in the 1.1–8 keV band as seen by eROSITA’s TM1–TM7. The colors are linearly distributed so that black corresponds to a pixel intensity value of 0 and white to 3 cts pixel−1. The point source at the southeastern rim of the Vela Jr remnant is the 65 ms pulsar PSR J0855-4644.

4 Spectral analysis

All the spectra were fitted with the software PyXSPEC, the Python interface of XSPEC (Arnaud 1996). The given errors correspond to the 1σ confidence interval and the statistic employed is CSTAT (Cash 1979). Solar abundances were set to those of Wilms et al. (2000).

4.1 The eROSITA survey data eRASS:4 of Vela Jr

In this chapter we present a systematic spectral study of Vela Jr by using the stacked data of the first four eROSITA all-sky surveys (eRASS:4).

4.1.1 Voronoi binning approach

We gained a first impression of the global characteristics of non-thermal emission across Vela Jr via a Voronoi binning approach. This analysis was carried out in order to provide a useful look at the nonthermal shell of Vela Jr, which is agnostic to any morphological features in the emission. We followed a similar approach as in Mayer et al. (2022), which can be summarized as follows. We ran an adaptive Voronoi tessellation algorithm (Cappellari & Copin 2003), with a target S/N = 100 on the eRASS:4 count image of Vela Jr, which was extracted in the 0.2–2.3 keV band. The resulting Voronoi bins were fed into srctool to extract spectra from each region using only data from telescope modules which have an on-chip filter (i.e., TMs 1, 2, 3, 4, 6). These spectra were then fitted in the energy range 0.2–8.5 keV with a source model combining foreground absorption (Wilms et al. 2000) with a thermal contribution from shocked plasma with non-equilibrium ionization (Borkowski et al. 2001) to fit the contribution from the Vela SNR. This was combined with a power law component that accounts for the presence of nonthermal emission from the Vela Jr remnant itself. In order to avoid the challenging task of independently constraining the absorption toward Vela Jr (see Sect. 4.2), we fixed the absorption of this nonthermal component to that of the CCO, N_H = 3.5 × 10²¹ cm⁻² (Becker et al. 2006). In Xspec, the combined source model is thus expressed as TBabs*vpshock+TBabs*powerlaw. The background contribution was modeled using a combination of three templates, reflecting the instrumental background (Yeung et al. 2023), the extragalactic X-ray background due to unresolved AGN (De Luca & Molendi 2004), and a combination of thermal emission models. The last consists of components taking into account emission from the local hot bubble, the Galactic halo, and unresolved emission from the Galactic disk. The relative shape of the thermal background was fixed to that fitted to the spectrum of an elliptical background region located outside the Vela SNR shell, centered on (α, δ) = (08^h44^m30^s, −39°26′), and with semi-axes of 2.6° × 1.3°.

From the fit results, we created maps illustrating the spatial distribution of the best-fit physical parameters across Vela Jr, which are displayed in Fig. 4. One can clearly see the signature of the Vela Jr shell in the map of the nonthermal brightness Σ_Γ as it exceeds the local background outside its shell by around an order of magnitude. While the brightest emission is clearly detected along the northwestern, western, and southern rims, there is a significant spatially inhomogeneous excess of emission inside the shell. We note, however, that the map indicates a nonzero level of what appears to be nonthermal emission outside the shell of Vela Jr. We believe this is likely to be caused by the contribution of a hotter thermal emission component in the Vela SNR, and probably does not correspond to a truly nonthermal contribution. In fact, modeling the thermal background as a combination of two thermal components in equilibrium (vapec+vapec) is a valid alternative approach (see Lu & Aschenbach 2000), which we tested. We recovered a very similar morphology for the nonthermal emission of Vela Jr, while the detected level of emission in the region outside was reduced. Importantly, even without considering alternate models, the measured level of nonthermal emission on the inside of the shell visible in Fig. 4 is significantly above the background level outside, especially in the western half. Thus, the presence of diffuse nonthermal emission on the inside of the shell of Vela Jr appears likely, whereas we do not claim to observe diffuse nonthermal emission from the Vela SNR. The properties of the suspected diffuse nonthermal component in Vela Jr are investigated in more detail in Sect. 4.1.2.

The distribution of the photon index Γ of the power law component appears to indicate comparatively hard nonthermal emission along the northwest and especially the southeast portions of the rim of Vela Jr. However, we believe that some skepticism is warranted regarding the physical origin of suspiciously small photon indices (Γ ~ 1.4–1.8) in the southeast. On one hand, visual inspection of selected spectral fits in this region shows that the detected nonthermal emission likely does have an astrophysical origin within the shell of Vela Jr, as its level is above that expected for nonthermal or instrumental backgrounds by a factor of a few. On the other hand, a possible issue with our approach is that we assume the absorption toward the entire Vela Jr remnant to be compatible with that of the CCO (N_H = 3.5 × 10²¹ cm⁻²). If the true absorption were much higher than that in the southeast, the nonthermal component would experience significant downward curvature even above 1 keV, artificially reducing the inferred power law slope, in particular since the thermal component is relatively soft in this region. We tested this hypothesis by carrying out an array of simple XSPEC simulations using parameters for the thermal component typical for the southeast, but fixing the intrinsic power law photon index to Γ = 2.5 and varying the true absorption of the nonthermal component in the range 4−14 × 10²¹ cm⁻². This showed that, when fitted with a fixed absorption as described above, the measured photon index reaches values of 1.8 and 1.4 at intrinsic column densities of N_H = 9 × 10²¹ cm⁻² and N_H = 1.2 × 10²² cm⁻², respectively. The physical origin of this suspected additional foreground absorption may be in the Vela molecular ridge (VMR), which spatially overlaps the southeastern part of the shell, if we assume that the shell of Vela Jr is located within or behind it, requiring a distance to Vela Jr ≳700pc (Liseau et al. 1992). For completeness, we note that the pulsar PSR J0855–4644 and its faint X-ray nebula overlap the southeastern rim as well (Maitra et al. 2018). However, since the extent of this nebula is limited to less than an arcminute in the Chandra data (Maitra et al. 2017), and no extended emission from it is visible in the eRASS:4 intensity image (see Fig. 3), it seems unlikely to be responsible for the hard emission detected in the southeast rim (see Acero et al. 2013, for a discussion of the pulsar and of the association with the VMR).

Fig. 3 Energy-selected eROSITA (TM1-TM7) view of the Vela Jr region. The images are centered on the position of the CCO. Intensity vs. color are distributed according to the ASINH-function. The color vs. cts pixel−1 distribution in each image is depicted. For the 0.7–1.0 keV images a scale of 0–6 cts pixel−1, and for the 1.1–1.3 keV images a scale of 0–2 cts pixel−1 were used. Below ~ 1 keV emission from the Vela SNR dominates. Above ~ 1 keV the nonthermal emission from Vela Jr starts to stand out above the emission from Vela SNR.

Fig. 4 Properties of nonthermal emission across Vela Jr. The top panels illustrate the distribution of the photon index Γ (left) and the surface brightness ΣΓ (right), that is, the integrated flux per surface area of the nonthermal model components in spectral fits to Voronoi bins. The energy range for the computation of the integrated nonthermal flux is 1.0–5.0 keV. The bottom panels illustrate the corresponding 1σ uncertainties of photon index (left) and surface brightness (right) measurements in each bin. The gray circle gives the extent of Vela Jr, and is identical to the circle displayed in Fig. 2. The position of the CCO and the energetic pulsar PSR J0855–4644, located in the southeast region of Vela Jr, are respectively indicated with a gray and green plus sign (+).

4.1.2 Spectral modeling of selected regions

Due to the limited photon statistics in the eROSITA survey data of Vela Jr, we divided the entire remnant into seven distinct regions, which allowed us to obtain spectral parameters with higher accuracy than is possible in the Voronoi binning approach. As indicated in Fig. 5, the defined areas include three boundary shock regions and four sectors located well inside the remnant. From Figs. 1 and 3 it is clear how the contribution of the Vela SNR is dominant at energies below ~ 1 keV and how Vela Jr contributes to the emission above this energy. In a first attempt, we therefore tried to restrict our spectral analysis to the 1–5 keV energy band to exclude most of the contribution from the Vela SNR. However, after a detailed analysis, we concluded that this interval did not allow us to fully constrain all the physical spectral components. Finally, we used the full 0.3–10 keV bandwidth for the spectral analysis and fitted the contribution from the Vela SNR first in order to get a background model. Building up on that, we fitted the source spectra for the distinct regions of Vela Jr simultaneously with the background model. It is worth noting that the background was taken from an elliptical ring surrounding Vela Jr so that the background model obtained this way may be only an approximation of the real Vela SNR contribution for regions within Vela Jr. Details of the background modeling and spectral fitting are described in Appendix A.

Starting from these, we employed an approach based on the Markov chain Monte Carlo (MCMC) technique to have a better understanding of the correlation between the various spectral parameters of the tested models. This delivers a more robust error estimation (see, e.g., van Dyk et al. 2001; Borkowski et al. 2010, respectively for a review of the method and an applied case to SNR). We ran the Python library emcee-based code (Foreman-Mackey et al. 2013) with 40 000 steps to make sure the chains could converge. We then manually inspected the trace plot of the chains, selecting only the last 2000 steps. This choice ensured that most of them were converged. Even though MCMC-based approaches are not designed to give a quantitative estimate of the goodness of the fit (unlike the χ² estimator), we present in Table 2 the values of the median Cash statistic and degree of freedom values drawn from the converged MCMC chains. In this sense, we employ the ratio CSTAT/d.o.f. and also the reduced ${\chi }^2/\text{d.o.f.=}{\chi }_{\text{d.o.f.}}^2$.

The first model tested for each of the seven regions was VPSHOCK, which represents a single temperature plane-parallel shocked plasma with non-equilibrium collisional ionized material (Borkowski et al. 2001), allowing us to fit the single abundance parameters. We left free to vary C, N, O, Ne, Mg, Si, and Fe. We also tested the power law model in order to compare a featureless model spectrum against a thermal one. Both models included Galactic absorption using the TBabs model.

Starting from the NW rim region, we fitted a dominating nonthermal model to the featureless spectrum above 1 keV. Comparing the results of the different models, the reduced statistic values of VPSHOCK and power law are identical (χ²/d.o.f. = 1.1, see Table 2). To assess which was the most reliable model, we compared the corner plots given by the two models, concluding a power law model gives more robust constrained parameters. Our preferred conclusion is that the spectrum is featureless and a power law is the best-fitting model for it.

Regarding the west region (labeled “W” in Fig. 5), it seems to be dominated by nonthermal emission as well, showing no strong lines above 1 keV. However, comparing it to the NW rim region, the spectrum of the W rim appears flatter. Using a power law model, we observe Γ = 3.05 ± 0.17 at the NW rim against Γ = 2.5 ± 0.12 of the W rim obtained with the same model (Fig. 6). We further observed that the W rim has a lower value for the column density than those of other regions, while the normalization and photon indices are more similar to each other. Another featureless spectrum is obtained from the S rim. The best-fitting model is clearly a power law for this region (Fig. 7). The photon index and column density values for that region are indeed similar to those of the NW rim.

The four inner regions of the remnant are all characterized by the presence of several unresolved lines above 1 keV, while the outer regions, for example NW, W, and S, show featureless spectra at this energy range. Sectors 2 and 3 are basically dominated by background emission (which is emission from the Vela SNR). The VPSHOCK and power law model fits for sectors 1 and 4 have a comparable goodness (see Table 2), whereas for sectors 2 and 3 the spectral fits are definitely worse when a power law model is fitted.

Regarding sectors 2 and 3, we tested a two-component model made by a VPSHOCK plus a power law component (i.e., a thermal plus a non-thermal component). In both regions we obtain good fitting results; the extra thermal component consistently improves the fit, compared to the single power law model (see Table 2). If in sector 2 the power law photon index is set to a value close to those of the other regions (Γ ~ 2.5) then the thermal component is largely unconstrained. Conversely, the thermal component is stronger in sector 3 and the power law is not well constrained. However, a combination of the two is always needed to improve the fit statistic. We note again how these spectra suffer from degeneracy of several factors, mainly the superposition of the Vela SNR and Vela Jr along the line of sight. In conclusion, the spectrum of sector 2 is probably featureless, as are those of sector 1 and 4, while in sector 3 an extra thermal component (probably due to the incomplete background modeling of the Vela SNR contribution) is needed to have a good fit. In Fig. 6, it is shown how the boundary regions selected in NW, W, and S have strong nonthermal spectra. Therefore we conclude that the spectrum of Vela Jr is uniformly nonthermal across the whole remnant.

Fig. 5 Vela Jr and the surrounding region in the 1–8 keV band, taken from stacked eRASS:4 data. Source regions selected for spectral analysis are highlighted in red, whereas the region used to model the background spectrum is indicated in blue.

Fig. 6 Distribution of spectral parameters for the power law model when fitted to the spectra of the selected regions defined in Fig. 5. The photon index Γ plot clearly shows how sectors 1 and 4 have well constrained values, similarly to the NW, W, and S regions, while sectors 2 and 3 do not.

4.2 eROSITA CalPV and XMM-Newton data of the NW rim

In addition to the survey observations described in the previous section, eROSITA observed the NW rim of Vela Jr during its CalPV phase in a 60 ks deep pointed observation. XMM-Newton has also observed the NW rim of Vela Jr on a regular basis for the purpose of instrument calibration. In total, we found six data sets in the public archive (see Table 1) taken between 2001 and 2021, which sum up to a total exposure of more than 410 ks. We used this eROSITA CalPV and XMM-Newton archival data to search for temporal spectral changes in the emission properties of the NW rim and to further constrain the results obtained from eROSITA’s survey data. The much higher photon statistics of this data set allowed us to apply a spectral analysis for seven distinct regions of the NW rim (see Fig. 8). These regions were tailored according to Pannuti et al. (2010) in order to make the results comparable with the one obtained from a Chandra observation of that spot. The background spectra for XMM and the CalPV data were extracted from an area north of the NW rim, located at RA = 08h48m42s, Dec = −45°30′18″ and covering an elliptical area of approximately 33 × 11 arcmin.

The XMM-Newton observations were taken in calibration mode, except the one from 2001. This mode can create two disturbing effects in the data. The first is an unnatural rise in the continuum spectrum due to Bremsstrahlung produced by the production of Auger electrons from the onboard calibration source. Both the source and background spectra are affected by that, but the effect is not homogeneous over the detector. The second effect is an uneven illumination of the PN and MOS 1/2 detectors, making the background even more varying. After a careful background modeling (see Appendix B), we did not observe any significant temporal variation in the spectra extracted from each of the six archival XMM-Newton data sets taken between 2001 and 2021. This allowed us for each of the seven regions to fit the spectra from all XMM-Newton observations simultaneously. The energy band was set to 0.3–10 keV. From this analysis, we obtained the best-fit value for the photon index and the column density of each of the seven regions (Table 3), averaged over six different observations. We then applied the MCMC approach as in Sect. 4.1.2. For the XMM-Newton data set we used 5000 steps because of the very high photon statistics available, whereas for the eROSITA CalPV data we had to apply 10000 steps because of the lower photon statistics.

Based on the results obtained in the previous sections, we tested an absorbed power law model first. In order to correctly model the absorption, we employed the Tuebingen-Boulder model (TBABS in XSPEC; Wilms et al. 2000). The fit gave well-constrained parameters in all regions. Based on the XMM-Newton data, the fitted column density is found to be in the range ≃2−8 × 10²¹ cm⁻² (Table 3). With an uncertainty of typically 0.1 × 10²¹ cm⁻² there is just a small variation in the column absorption among the seven defined regions of the NW rim. The fitted photon indices show a similar soft variation and are fitted in the range 2.6–3.1. We ran the same fit on the CalPV data and found good agreement for all regions, except region 5. Comparing Table 3 and Table 4 it seems CalPV derived parameters are systematically higher by approximately 20% for regions 1, 2, and 3 compared to those derived with XMM-Newton. These differences are easily explained with the early calibration status of eROSITA.

After testing a power law, we tried a more complex model in order to characterize the underlying electron population. The SRCUT model (see Reynolds 1996, 1998, for a comprehensive discussion), which is implemented in XSPEC (Arnaud 1996), convolves the synchrotron volume emissivity of a single electron with an exponentially cut off power law electron distribution. In other words, this model sets an upper limit to the cutoff energy of the electrons, without further physical assumptions (for a detailed description of this model, see also Reynolds & Keohane 1999). The SRCUT model takes as parameter the radio spectral index α (f(v) ∝ v^−α, where f is the radio flux) and the normalized radio flux at 1 Jy for the normalization. Therefore, in order to exclude a degeneracy between the parameters, one has to rely on radio measurements (either flux or radio spectral index) to constrain the break frequency. We thus fixed the radio spectral index to -0.4 (see, e.g., the discussion in Maxted et al. 2018) and left the normalization free to vary. As with the power law model, for each region we loaded all the spectra from the different observations and fitted them together. The background was handled as described in Appendix B. In all regions we find a cutoff frequency consistent with 10¹⁶ Hz. Following Pannuti et al. (2010), we can convert that frequency (parameter “break” in SRCUT) to the actual cutoff energy of the electrons according to $${\nu }_{\text{cutoff}}\approx 1.6\times {10}^{16}\left(\frac{\text{B}}{10\text{\hspace{0.17em}μG}}\right){\left(\frac{E_{\text{cutoff}}}{10\text{\hspace{0.17em}TeV}}\right)}^2\text{Hz,}$$(2)

where B is the magnetic field and E_cutoff is the cutoff energy of the electron population. If on the one hand the TeV emission detected from the NW rim of Vela Jr. suggests a small value for the magnetic field, on the other hand a much higher value is needed to explain the small width of X-ray-emitting synchrotron filaments (Bamba et al. 2005a; Kishishita et al. 2013). If we assume a primarily leptonic origin of the TeV emission and the electron cutoff energy of 27 TeV, which is the value inferred by H.E.S.S. Collaboration (2018b) for the leptonic case, we can apply Eq. (2) to the values in Table 4, giving an estimate on the magnetic field. We obtain that the magnetic field can range from 2 µG in region 5 up to 16 µG in region 7. It is also possible to use the inverse argument and calculate the maximum energy reached by the electrons, in this case assuming a value for the magnetic field.

Finally, we implemented in XSPEC the following model (Zirakashvili & Aharonian 2007, 2010), $$F\left(E\right)\propto E^{-2}{\left[1+0.38{\left(\frac{E}{E_0}\right)}^{0.5}\right]}^{11/4}\exp \left(-\sqrt{\frac{E}{E_0}}\right),$$(3)

which we called 1cutoff. This model represents the spectrum of an electron distribution accelerated by a shock wave, limited by radiative losses (see, e.g., Sapienza et al. 2022, for a recent application). The relation provides a direct estimation of the expansion velocity V_f, independently of other factors than the parameter E₀: $$E_0=\frac{2.2\text{\hspace{0.17em}keV}}{{\left(1+{\kappa }^{1/2}\right)}^2{\eta }_{\text{B}}}{\left(\frac{V_f}{3\times {10}^3\text{km\hspace{0.17em}}{\text{s}}^{-1}}\right)}^2.$$(4)

The parameter E₀ depends on κ, the ratio of the upstream to the downstream magnetic field, and η_B, which accounts for the deviations from the Bohm diffusion regime. We assume $k=1/\sqrt{11}$, which is, according to Zirakashvili & Aharonian (2007), realistic for young and nonthermal SNRs. Rearranging Eq. (4), we can write a formulation with an explicit dependence on η_B (Table 5). Since the expansion rate µ measured by Allen et al. (2015) is mostly valid for our regions 3 and 4, assuming as expansion velocity 2140 km s⁻¹ (mean value $V_f{\eta }_{\text{B}}^{-1/2}$ for these two regions), we can derive the following approximate estimate of the distance (D = V_f/µ): $$D\approx \frac{2140{\text{km s}}^{-1}}{420{\text{mas yr}}^{-1}}{\eta }_{\text{B}}^{1/2}\approx 1100{\eta }_{\text{B}}^{1/2}\text{pc.}$$(5)

Assuming that our assumptions and modeling results are valid, the expression obtained in this way directly relates the distance to Vela Jr to the Bohm factor η_B. In particular, it implies that distances significantly below 1 kpc can be considered unlikely, as they would require η_B < 1.

Fig. 7 eRASS:4 spectra (from top left to bottom: NW, S, W, sector 1, sector 2, sector 3, and sector 4) fitted with non-equilibrium parallel shock thermal model (VPSHOCK). Red indicates the source model, green and blue represent the components of the X-ray background model, and magenta indicates the most significant component of the particle background. The spectra were rebinned for graphical purposes.

Fig. 8 Exposure- and vignetting-corrected eROSITA image of the NW rim of Vela Jr, restricted to the 1 1–8 0 keV energy band The labeled subregions indicate the sectors whose spectra were extracted The image pixel size is 4″ and a Gaussian smoothing with σ = 1.3 pixels has been applied. The color scale ranges from zero (dark blue) to a maximum of 0.0011 counts s−1 (yellow).

5 Distance estimate and the role of the Vela SNR

Looking at our survey and pointed observation results, we obtain column densities ranging from 10²⁰ and 10²² cm⁻², which is in agreement with the values presented in Sect. 1. Such a large variation can be addressed observing the data in other bands; for example, in Fig. 6 of Maxted et al. (2018) there seems to be some material overlapping the NW rim, along our line of sight, which could justify a measured column density of 10²¹ cm⁻² or higher. In addition, the synchrotron emission detected by Fukui et al. (2017) strongly correlates with the TeV emission observed by Aharonian et al. (2007) in the same region. This would justify the production of TeV photons via the hadronic channel, which requires a dense material to scatter the high energy protons coming from the SNR shock. There are two possible scenarios giving rise to the γ-ray emission (Lee et al. 2013; H.E.S.S. Collaboration 2018b). The first possibility is the so-called hadronic scenario, as the TeV photons would be produced in the decay of neutral pions, which are created during the violent collision between two protons. The second scenario, the lep-tonic model, involves ambient relativistic electrons upscattering seed photons up to TeV energies via the inverse-Compton (IC) mechanism.

In conclusion, the presence of material close to the acceleration region would justify the signal observed in the TeV band, assuming its origin is mainly hadronic. Nevertheless, some other parts of Vela Jr appear completely free of the covering HI emission. From all these different observations, we conclude that the column density is a parameter that is strongly variable across the remnant, probably due to the presence of a large amount of material along the line of sight.

Motivated by this evidence, we looked at the Gaia-based extinction database given in Lallement et al. (2019); the observations provide parallax together with extinction measurements for each single star, giving an indication for the distance of the absorbing material. Becker et al. (2006) observed the CCO with XMM-Newton and found N_H = (3.8 + 0.3) 10²¹ cm⁻². By using the relation of Predehl & Schmitt (1995) $$N_{\text{H}}\left[{\text{cm}}^{-2}/A_{\upsilon }\right]=1.79\times {10}^{21},$$(6)

we obtain the expected extinction A_υ ~ 1.92. Comparing this to the Gaia data set of Lallement et al. (2019), this extinction corresponds to a distance of 2.4 kpc (Fig. 9). Considering the radius is θ ~ 0.95° and the distance d = 2.4 kpc, we can estimate the physical radius as r = θ × d ~ 40 pc. This would imply Vela Jr is a very large remnant if it turns out to be young: this is around half the physical radius of the Vela SNR, which is at least ~ 11 kyr old. In conclusion, it is unlikely this is the true age of Vela Jr.

Even though the correlation between optical extinction and column density is very well established on large scales up to more than 10 kpc, there may be local deviations. Dust destruction by the shock wave of the supernova itself could be a reason (Micelotta et al. 2016; Zhu et al. 2019) for deviations from Eq. (6). However, we note that the Lallement et al. (2019) data show a steep increase in A_υ at around 750 pc (Fig. 9), which implies the distance is greater than or equal to this value unless the association of A_V and N_H is incorrect by a factor ~ 10. Given the large amount of uncertainties, this is in accordance with the lower limit set by Allen et al. (2015) using Chandra expansion measurements at 500 pc. This may imply a hypothetical dust cloud at the same distance as Vela Jr, especially since massive progenitors are very likely to be found in dust-rich star-forming regions. Considering the result of Allen et al. (2015) and the uncertainties associated with the Eq. (6), the new evidence provided by the optical data of Lallement et al. (2019) is in agreement with the findings of previous papers, leading us to conclude that a distance of Vela Jr of 750 pc or higher is very likely.

Fig. 9 One-dimensional extinction map obtained from Gaia extinction cube and extracted using the tool (https://astro.acri-st.fr/gaia_dev/) in the direction of CXOU J085201.4-461753 (l = 266.2°, b = −1.2° in galactic coordinates).

6 Summary and discussion

In Sect. 6 we concluded that Vela Jr. has an overall featureless spectrum. This confirms the early indications that the NW rim of Vela Jr. is likely to be an acceleration site for high energy particles. Starting from a single power law, we tested different increasingly more complex models, specifically SRCUT and the so called lcutoff, to characterize in the detail each region. For the latter, we considered the electron cutoff energy derived with HESS (27 TeV for the leptonic case, considered more suitable than the hadronic scenario given the observation of X-ray synchrotron emission), finally obtaining the magnetic field can vary between 2 µG and 16 µG. Employing lcutoff, we manage to estimate the distance as 1100 pc, but we note how this value is a function of the Bohm factor η_B, which can significantly vary across the remnant (Tsuji et al. 2021). In conclusion, our findings confirm this region as probable acceleration site for particles, in accordance with the observations reporting TeV emission from the same direction. This scenario collocates Vela Jr. firmly in the class of the so called shell-type supernova remnants, with the eROSITA observation of the whole remnant confirming the results of the NW rim pointed observations.

Putting these findings in the broader context of TeV observations, as mentioned in Sect. 1 only 8 TeV emitting sources are firmly identified with SNRs, and all of them are shell-type. This is quite a small number if the observed number of Galactic cosmic rays is supposed to be due mainly to acceleration in shock fronts. Motivated by this observation, Aschenbach (2016) discusses how three well-known TeV emitter SNRs, SN1006, RX 1713.7-3946 and Vela Jr may be good objects to study this discrepancy. He proposed a modified version of the Sedov-Taylor equations solution in order to account for a loss of energy in cosmic rays, instead of being dispersed in the expansion of the remnant. Specifically, the shock front is reaccelerated after the explosion. The assumed fraction of energy transferred to cosmic rays for the modified ST version is A = 0.9075 (see Fig. 10). This would place the explosion at around 1200–1300 AD, making it compatible with the hypothesis of a much closer distance proposed by Burgess & Zuber (2000). This discussion is based on the ⁴⁴Ti detection reported in Sect. 1. However, such a small distance is ruled out by past observations (see Sect. 1) and our new ones (Sects. 4, 5). In conclusion, the model proposed by Aschenbach (2016) provides an explanation for the small number of observed objects of this kind because the energy lost in particle acceleration accounts for the low surface brightness of these remnants.

However, cosmic ray acceleration as the supplier of additional momentum in the late stages of SNRs was treated with a different approach by Diesing & Caprioli (2018). The cosmic ray impact on the expansion of the remnant is boosted by the presence of dense ISM or even by molecular clouds. Our detection of multiple velocities and different energies in the NW rim, gives solid support to a picture with the remnant expanding in a nonuniform medium. Moreover, similar cosmic ray acceleration efficiencies were theoretically derived by Haggerty & Caprioli (2020) and observationally tested by Giuffrida et al. (2022). Instead of starting from Sedov-Taylor solutions, all these papers describe microscopically the transfer of energy from the shock to cosmic rays. We find that these models need a much lower transfer of energy in cosmic rays, compared to what was proposed by Aschenbach (2016). Additional theoretical studies connecting the Sedov-Taylor solutions with a microscopic particle acceleration treatment would be very useful to ultimately have a unified picture, eventually determining the cosmic ray acceleration efficiency.

Comparing the POWERLAW and VPSHOCK models, we conclude Vela Jr. has an overall featureless spectrum with region NW, S, and W presenting flatter photon index values than inner sectors 1, 2, 3, and 4. Having steeper spectra in the inner regions is an effect likely associated with synchrotron cooling, already discussed for the case of Vela Jr. by Kishishita et al. (2013). This is not surprising since the shock regions should be more energetic due to the interaction with the ISM, while the inner early shocked regions should progressively cool toward the center, unless a reverse shock is in action. However, we note how the uncertainties for the inner sector measurements are large, as shown in Fig. 6. In conclusion, we have significant indications that the spectra of Vela Jr. are mainly nonthermal.

We also observe how sector 2 and sector 3 are poorly described by a power law, and only a two-component model definitely improves the fit. The additional thermal component is probably needed due to the presence in these regions of bright thermal filaments. Further studies, for instance optical spectroscopy of the Ha line in order to detect two possible velocity components, might help to associate them finally with Vela SNR or Vela Jr.

Even though the most straightforward explanation for such different spectra is an excess of emission by Vela in that region, we cannot exclude that this change in the plasma can be associated with an extra thermal component belonging to Vela Jr. However, we observe with the two-component modeling super solar abundances in the thermal component. This is similar to abundances measured in other regions of Vela SNR by Mayer et al. (in prep.), which makes us speculate that the enhanced abundances belong to this remnant and not to Vela Jr. Further studies would be needed to clarify the origin of this additional thermal component.

However, in the Bayesian framework we can try to assess which model fits the data better by employing the Bayesian information criterion (BIC). It is defined as $$\text{BIC}=-2\log L_{\max }+k\log N,$$(7)

where L_max is the maximum likelihood of the model, k is the number of free parameters, and N the number of data points. We assumed the median likelihood of the model represents the maximum likelihood. According to Wall & Jenkins (2012), the model that retrieves the smallest value of the BIC is to be considered the best to fit the data set. From Table 6, the power law model has to be preferred for the boundary shock regions (i.e., NW, W, and S, and sector 1). On the other hand, sectors 2 and 3 clearly favor a power law with an added thermal component, from a purely statistical point of view. For sector 4 the improvement of the fit is only marginal. Therefore, our overall conclusion is that the spectrum of Vela Jr is uniformly featureless across all the remnant, even though (especially in sectors 2 and 3) it is unclear whether the additional thermal component can be associated with Vela Jr or Vela SNR. To test this scenario, as a consistency check we ran again the VPSHOCK fit, this time linking all the abundances to the O parameter. The result is given in Table 7 and clearly shows how leaving all the abundance parameters (C, N, O, Ne, Mg, Si, Fe) free to vary independently is needed in sectors 2 and 3, indicating a strong variation in the plasma condition. This supports our idea that the background varies in those regions, given also the evidence provided in Sect. 4.1.1.

As presented in Sect. 3, we addressed the question of the geometric center of the remnant. Assuming the center found by fitting a circle only through the northern remnant rim represents the true explosion site, we showed that the CCO is practically still at its birthplace. Although there are two Chandra data sets in the public archives that were taken several years apart from each other, missing X-ray bright calibration stars in the Chandra field of view prevented us from taking detailed proper motion measurements on the Chandra CCO data itself. Mignani et al. (2019) proposed the existence of an IR counterpart of the CCO, which they put into question due to its missing proper motion. They set an upper limit of less than 10 mas yr⁻¹ which would not be in agreement if the CCO had traveled from the ROSAT data-based geometrical remnant center which is 4 arcmin apart (Aschenbach 1998). However, given the new geometrical remnant center determined in this work, the missing IR candidate proper motion is in full agreement with our finding that the CCO is practically still at its birthplace, invalidating the argument of a missing proper motion of the IR counterpart.

There are various reasons why the shape of the remnant is not exactly spherical in all places. One possible explanation for its small asymmetry is a massive progenitor, specifically a luminous blue variable (LBV). Such a progenitor would account for low density and high explosion velocities, predicted in order to have a low surface brightness remnant (Aschenbach 2016). Remarkably, our expansion velocities measured in the NW rim (Table 5) are all quite high, which fits this picture. Recent papers (Ustamujic et al. 2021; Chiotellis et al. 2021; Das et al. 2022) clearly demonstrate how the wind emitted from a massive progenitor can also shape the circumstellar medium (CSM), and this affects the final shape of the remnant after the supernova explosion. However, a preexisting difference of density in the ISM among the different regions, especially if the remnant is large (up to 40 pc), can alternatively explain the small asymmetry. In conclusion, the primordial differences of density in the ISM and/or the action of the wind from a massive progenitor on the CSM can be a very likely explanation for the small asymmetric shape of Vela Jr. An even simpler and probably the most straightforward explanation would be that of an intrinsically asymmetric explosion that then leads to a deformation from a perfect spherical shape.

On the other hand, the small asymmetry that manifests itself especially in the southern rim of Vela Jr can be described by a second circle itself, posing the question of whether Vela Jr could have been formed by two SN explosions. In this scenario, one formed the CCO and the northern part of the remnant and the second SN formed PSR J0855-4644 and the southern part of the remnant. This hypothesis justifies our initial fitting of the NW rim only. The presence of the other young and powerful rotation-powered pulsar PSR J0855-4644 within the supernova remnant boundary makes this an interesting hypothesis, which Fig. 11 tries to illustrate. Nonetheless, the association of PSR J0855-4644 with Vela Jr was considered unlikely by various authors in the past, mostly because of the distance mismatch between Vela Jr and the pulsar. Its radio dispersion measure-based distance was found to be 5.4 kpc, which greatly exceeds the distance estimates for Vela Jr, regardless of all their uncertainties. However, Acero et al. (2013) concluded from a comparison of column densities a new upper limit to the distance of PSR J0855-4644, which could be as close as 900 pc, bringing the pulsar distance in a range that seems to be in agreement with the distances discussed above for Vela Jr. The question arises of whether it would be possible for a pulsar formed at the center of this putative southern circle to travel to its location in the remnant shell in a time span that is in agreement with the possible age of Vela Jr. The angular distance between the center of the southern circle shown in Fig. 11 and the position of PSR J0855-4644 is on the order of 0.5 deg. Assuming a distance of 1 kpc and a pulsar proper motion velocity of V_PSR = 1000 km s⁻¹ (implying a proper motion of 0.2 arcsec yr⁻¹), it would require ~8500 yr to reach its location in the remnant shell. Although this assumed pulsar velocity is at the upper end of the range of possible values, considering the uncertainties in these parameters this scenario cannot be completely ruled out on the grounds of the required velocity alone. It does however seem to be in contradiction with the expansion-based age estimate of 2.4–5.1 kyr for Vela Jr by Allen et al. (2015).

Maitra et al. (2017) detected a compact but elongated X-ray PWN around PSR J0855-4644. Assuming that the elongated diffuse emission indicates the pulsar's proper motion direction, it would not be in agreement with a birthplace in a second southern remnant. Using uGMRT continuum observations of the pulsar surroundings, taken at 1.35 GHz, Maitra et al. (2018) found a bow shock nebula in front of the pulsar, indicating that the pulsar’s ram pressure p_ram = ρ_ambV² exceeds the ambient gas pressure, probably due to a high pulsar velocity. Although no proper motion measurement of that pulsar have been published to date, it is possible to get at least a rough estimate for its proper motion direction from the geometry of its bow shock nebula. Both estimates for the pulsar proper motion directions are indicated in Fig. 11. As can be seen, they are more in agreement with a direction pointing toward the position of the CCO rather than toward the center of a putative hypothetical second SNR which formed the southern part of Vela Jr. If we suppose a common origin of the CCO and PSR J0855-4644, perhaps in a binary disruption event, an important question to ask is whether the pulsar could have reached its present-day location within the age of Vela Jr. Considering the present-day distance from the CCO to the position of PSR J0855-4644 of ~0.74°, and the measured SNR expansion rate of (13.6% ± 4.2%) kyr⁻¹ (Allen et al. 2015), even in the limit of free expansion, the required proper motion of the pulsar would be at least (0.36 ± 0.11) arcsec yr⁻¹. This however implies that the pulsar needs to have moved by around ~1.5″ within the 4 yr separating the two existing Chandra observations of the object (observation IDs 13780 and 18640), which can be excluded even by a brief inspection of the two data sets. Thus, overall, a physical association of PSR J0855-4644 with the explosion that formed the Vela Jr SNR has to be considered quite unlikely. It is however not excluded for the progenitor stars of the pulsar and Vela Jr to have an origin in the same star-forming region, as long as one allows for the pulsar to be significantly older than the CCO. This scenario is in agreement with the actual formation theory of CCOs (see, e.g., Enoto et al. 2019), which predicts that they are young neutron stars still embedded in the fall-back material of the supernova explosion.

Fig. 10 Expansion velocity from Sedov-Taylor (ST) equations and from the modified version of Aschenbach (2016). A peak appears in the latter formulation, implying a reacceleration of the front shock. The velocity range predicted for 0.75−1.1 kpc distance is also plotted.

Fig. 11 Vela Jr. in the 1.1–8 keV band. The left image indicates a possible circular geometry of the southern remnant region. The circle has a radius of 41 arcmin. The right image indicates possible proper motion directions of PSR J0855-4644 as indicated from the geometry of the elongated PWN seen in the Chandra data (blue) by Maitra et al. (2017) and from the geometry of the radio bow shock nebular (orange) according to Maitra et al. (2018).

Acknowledgements

We would like to thank all the eROSITA team for the helpful discussions and suggestions provided during the realization of the paper. We thank the anonymous referee for the useful and valuable comments. MS acknowledges support from the Deutsche Forschungsgemeinschaft through the grants SA 2131/13-1, SA 2131/14-1, and SA 2131/15-1. MF acknowledges support from the Deutsche Forschungsgemeinschaft through the grant FR 1691/2-1. WB acknowledges support from the Deutsche Forschungsgemeinschaft through the grant BE 1649/11-1 and thanks Bernd Aschenbach for inspiring discussions and Luciano Nicastro for his support in data processing. FC acknowledges support from the Deutsche Forschungsgemeinschaft through the grant BE 1649/11-1 and from the International Max-Planck Research School on Astrophysics at the Ludwig-Maximilians University (IMPRS). MGFM acknowledges support from the International Max-Planck Research School on Astrophysics at the LudwigMaximilians University (IMPRS). This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors, and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP) and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS software system developed by the German eROSITA consortium. This work makes use of the Astropy Python package⁵ (Astropy Collaboration 2013, 2018). A particular mention goes to the in-development coordinated package of Astropy for region handling called Regions⁶. We acknowledge also the use of Python packages Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020), PyLaTex⁷ and NumPy (Harris et al. 2020).

Appendix A Background modeling and fitting procedure in eRASS:4 data

Here we describe the background modeling carried out for the eRASS:4 data. In Figure 5 we highlight background and source regions in blue and red, respectively. We extracted the background from an elliptical region around Vela Jr, representing mainly emission from Vela SNR. The spectrum extracted from this region is called X-ray background to distinguish it from the instrumental background. The latter is the contribution given by the sum of the internal noise and the signal created by incoming particles. More specifically, the camera is continuously hit from any direction by cosmic rays creating an X-ray continuum, produced mainly via the Bremsstrahlung effect. The components described so far are included in the instrumental background, which is not multiplied by the ARF. The components employed for the modeling were the following:

1. Source model: TBABS × CONSTANT × MODEL

2. X-ray background model: TBABS × CONSTANT × (VPSHOCK+POWERLAW)

3. instrumental model: CONSTANT × FWC_MODEL

The X-ray background is modeled using a VPSHOCK model plus a power law. We fixed the photon index of the power law, representing the contribution of the unresolved extragalactic sources (mostly AGN), to Γ = 1.46. This is the standard value given in De Luca & Molendi (2004) and from the same paper we fixed also the normalization. This assumption was made considering that the background due to extragalactic sources should be almost spatially constant and not particularly dominant in this region, especially considering the amount of absorbing material along the line of sight.

The instrumental background was modeled employing the Filter Wheel Closed model (FWC) developed in Yeung et al. (2023). This model is a combination of power laws and Gaussian lines, adapted to the processing version c020. This model was multiplied by a CONSTANT initially set to the BACKSCAL keyword (since this component is modeled without ARF). The CONSTANT parameter was left free to vary during the fit.

Once the best fit of the background was achieved, we loaded it as background model for the source and fitted them together. The source model was multiplied by a CONSTANT set to REGAREA⁸ keyword expressed in sr, which represents the geometric area of the region, so it is kept frozen. Therefore, the model is normalized by the area of the source.

Without fitting the data again, we launched an emcee run leaving all the parameters of the source free to vary, while only column density (nH), temperature (kT), ionization timescale (Tau), and normalization (norm) were left free in the X-ray background. This choice was made to allow the background to slightly vary from region to region, in order to match the high variability of Vela SNR emission. For the instrumental background, only the constant in front of the FWC model is free to vary, while all the other parameters were frozen.

Appendix B Background modeling and fitting procedure of the XMM-Newton data

In this Appendix we describe in detail the background modeling of XMM-Newton data. As shown in Table 1, most of the observations are in calibration mode, which means bright emission lines are artificially created in the spectrum (Figure B.1). In the following we distinguish between X-ray background and instrumental background. The components employed for the modeling were the following:

1. Source model: TBABS × CONSTANT × MODEL

2. X-ray background model: TBABS × CONSTANT × (VPSHOCK+VAPEC+POWERLAW)

3. Instrumental background model: CONSTANT × (POWER-LAW+GAUSSIAN lines)

Fig. B.1 XMM-Newton spectra (PN, MOS1, MOS2 for each observation, for a total of 18) from region 4 loaded and modeled. Bright emission lines are present.

Briefly, the instrumental background is modeled employing a superimposition of eight Gaussian lines (energies fixed at 1.5 keV, 4.5 keV, 5.9 keV, 6.5 keV, 7.49 keV, 7.11 keV, 8.05 keV, and 8.82 keV; line width (σ) fixed to 0.001 keV) plus a power law with photon index fixed to 0.01.

The X-ray background is mainly due to Vela SNR and is modeled as superimposition of a VAPEC, VPSHOCK, and power law models. VPSHOCK (Borkowski et al. 2001) is a component accounting for thermal emission from material not in equilibrium, while VAPEC (Smith et al. 2001) accounts for emission from material in thermal equilibrium. Photon index and normalization were set as described in Appendix A, this time leaving the normalization free to vary. In addition to the physical meaning, we looked for an empirical background model that allowed us to handle the decaying intensity and the changing shape of the spectrum.

As the first step we loaded the background data: we ran a fit with all the norms of the model components of the instrumental and X-ray background left free to vary, with energies and width of the Gaussian lines frozen to the values listed above. This gave us the best-fit model for the background. Thus, we loaded this best-fit model of the background with the source data, freezing all the parameters except the normalization of the lines and a constant in front of the X-ray background component. This choice was made since the normalization of the lines and the continuum are variable, mainly for the action of the radioactive source, whose intensity decays with time. In the second run, we loaded only the source data that we fitted together with the background model, accounting for the variation in the instrumental background. After these steps, we froze all the parameters in the instrumental background and we launched an emcee run on the source model plus the global normalization of the X-ray background.

References

All Tables

All Figures

	Fig. 1 False color image of Vela Jr. The energy of the X-ray photons are color-coded: red (0.2–0.7 keV), green (0.7–1.2 keV), and blue (1.2–8.0 keV). In order to enhance the visibility of the diffuse X-ray emission, the images of the different bands were convolved with a Gaussian filter (expressed in pixels; red: radius = 7, σ = 3.5; green: radius = 8, σ = 4; blue: radius = 8, σ = 4).
In the text

	Fig. 2 Intensity representation of the emission of Vela Jr and its CCO J085201.4-461753 in the 1.1–8 keV band as seen by eROSITA’s TM1–TM7. The colors are linearly distributed so that black corresponds to a pixel intensity value of 0 and white to 3 cts pixel−1. The point source at the southeastern rim of the Vela Jr remnant is the 65 ms pulsar PSR J0855-4644.
In the text

Fig. 3 Energy-selected eROSITA (TM1-TM7) view of the Vela Jr region. The images are centered on the position of the CCO. Intensity vs. color are distributed according to the ASINH-function. The color vs. cts pixel−1 distribution in each image is depicted. For the 0.7–1.0 keV images a scale of 0–6 cts pixel−1, and for the 1.1–1.3 keV images a scale of 0–2 cts pixel−1 were used. Below ~ 1 keV emission from the Vela SNR dominates. Above ~ 1 keV the nonthermal emission from Vela Jr starts to stand out above the emission from Vela SNR.

In the text

Fig. 4 Properties of nonthermal emission across Vela Jr. The top panels illustrate the distribution of the photon index Γ (left) and the surface brightness ΣΓ (right), that is, the integrated flux per surface area of the nonthermal model components in spectral fits to Voronoi bins. The energy range for the computation of the integrated nonthermal flux is 1.0–5.0 keV. The bottom panels illustrate the corresponding 1σ uncertainties of photon index (left) and surface brightness (right) measurements in each bin. The gray circle gives the extent of Vela Jr, and is identical to the circle displayed in Fig. 2. The position of the CCO and the energetic pulsar PSR J0855–4644, located in the southeast region of Vela Jr, are respectively indicated with a gray and green plus sign (+).

In the text

	Fig. 5 Vela Jr and the surrounding region in the 1–8 keV band, taken from stacked eRASS:4 data. Source regions selected for spectral analysis are highlighted in red, whereas the region used to model the background spectrum is indicated in blue.
In the text

	Fig. 6 Distribution of spectral parameters for the power law model when fitted to the spectra of the selected regions defined in Fig. 5. The photon index Γ plot clearly shows how sectors 1 and 4 have well constrained values, similarly to the NW, W, and S regions, while sectors 2 and 3 do not.
In the text

Fig. 7 eRASS:4 spectra (from top left to bottom: NW, S, W, sector 1, sector 2, sector 3, and sector 4) fitted with non-equilibrium parallel shock thermal model (VPSHOCK). Red indicates the source model, green and blue represent the components of the X-ray background model, and magenta indicates the most significant component of the particle background. The spectra were rebinned for graphical purposes.

In the text

	Fig. 8 Exposure- and vignetting-corrected eROSITA image of the NW rim of Vela Jr, restricted to the 1 1–8 0 keV energy band The labeled subregions indicate the sectors whose spectra were extracted The image pixel size is 4″ and a Gaussian smoothing with σ = 1.3 pixels has been applied. The color scale ranges from zero (dark blue) to a maximum of 0.0011 counts s−1 (yellow).
In the text

	Fig. 9 One-dimensional extinction map obtained from Gaia extinction cube and extracted using the tool (https://astro.acri-st.fr/gaia_dev/) in the direction of CXOU J085201.4-461753 (l = 266.2°, b = −1.2° in galactic coordinates).
In the text

	Fig. 10 Expansion velocity from Sedov-Taylor (ST) equations and from the modified version of Aschenbach (2016). A peak appears in the latter formulation, implying a reacceleration of the front shock. The velocity range predicted for 0.75−1.1 kpc distance is also plotted.
In the text

Fig. 11 Vela Jr. in the 1.1–8 keV band. The left image indicates a possible circular geometry of the southern remnant region. The circle has a radius of 41 arcmin. The right image indicates possible proper motion directions of PSR J0855-4644 as indicated from the geometry of the elongated PWN seen in the Chandra data (blue) by Maitra et al. (2017) and from the geometry of the radio bow shock nebular (orange) according to Maitra et al. (2018).

In the text

	Fig. B.1 XMM-Newton spectra (PN, MOS1, MOS2 for each observation, for a total of 18) from region 4 loaded and modeled. Bright emission lines are present.
In the text