Free Access
Issue
A&A
Volume 618, October 2018
Article Number A155
Number of page(s) 11
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201832879
Published online 26 October 2018

© ESO 2018

1. Introduction

Supernova remnants (SNRs) are widely considered to be the main candidates for the acceleration sites of Galactic cosmic rays (CRs), which places them among the most interesting and most studied astrophysical objects. Vela Jr. is one of just a few SNRs with detected nonthermal emission and a well resolved shell-like morphology across the whole electromagnetic spectrum from radio to very-high-energy (≿100 GeV; VHE) gamma-rays. Its proximity, large size, and strong nonthermal emission (specifically from the north-western rim of the remnant) allow for detailed spatially resolved observations of the source both in X-rays and at TeV energies, making it one of the prime sources for studies of particle acceleration and radiation mechanisms in SNRs. Despite a rather late discovery of the source by the ROSAT satellite at X-ray energies above ∼1.3 keV (Aschenbach 1998) owing to its location in a very crowded region with the bright and large Vela SNR in the foreground, it became one of the most observed remnants across the electromagnetic spectrum.

Vela Jr. was one of the first sources detected in the VHE range (Aharonian et al. 2005) and benefits from a large amount of data collected with the H.E.S.S. telescopes since 2004. It was also detected by Fermi-LAT in the GeV energy range (Tanaka et al. 2011). Both Fermi-LAT and H.E.S.S. maps resolve a clear shell-like morphology of the remnant. The spectrum in the GeV range obtained by Fermi-LAT and the spectrum in the TeV range obtained by H.E.S.S. connect remarkably well and can be nicely described by a power law with an exponential cutoff with a cut-off energy of (6.7 ± 1.2stat ± 1.2syst) TeV (H.E.S.S. Collaboration 2018).

The process responsible for the gamma-ray emission from the remnant remains to be identified. Both hadronic (proton interactions with a subsequent π0-decay) and leptonic (inverse Compton (IC) scattering of relativistic electrons on ambient radiation fields) scenarios can adequately reproduce the shape of the gamma-ray spectrum (H.E.S.S. Collaboration 2018). A fit of the combined GeV–TeV spectrum with the IC emission from a parent electron population of the form of a power law with exponential cutoff in a pure leptonic scenario yields an electron spectral index of 2.33 ± 0.03stat ± 0.33syst and an electron cutoff energy of (27 ± 1stat ± 12syst) TeV. In the hadronic scenario, the fit of the observed spectrum results in a parent proton spectrum with a spectral index of 1.83 ± 0.02stat ± 0.11syst and a cutoff energy of (55 ± 6stat ± 13syst) TeV (H.E.S.S. Collaboration 2018). Both scenarios, however, face some difficulties. The hadronic scenario is supported by a good spatial correlation of the TeV emission with the distribution of the HI gas (Fukui et al. 2017), but the lack of thermal X-ray emission from the remnant places the upper limit on the ambient density at ∼0.03 cm−3 (Slane et al. 2001), which makes the hadronic scenario rather implausible as it would require the total energy in protons to be comparable to the explosion energy. For RX J1713.7–3946, which shares a lot of observational properties with Vela Jr., it was suggested that the lack of thermal X-ray emission can still be explained within the hadronic scenario if the remnant is expanding in a medium with small dense clouds that survive the passage of the shock and provide target material for cosmic rays (Inoue et al. 2012; Gabici & Aharonian 2014). However, this scenario also has some difficulties (Federici et al. 2015). The large amount of mass in clouds required for a fit of the high-energy spectrum in this scenario requires pre-existing clouds unaffected by the stellar wind of the progenitor star, which is rather implausible as the stellar wind would form Kevin–Helmholtz and Rayleigh–Taylor instabilities destroying the outer layers of the clouds. This would be followed by the formation of ∼(1−10) cm−3 streams of gas, which would strongly emit in X-rays when interacting with the forward shock (Federici et al. 2015).

On the other hand, high-resolution X-ray observations reveal thin and bright filamentary structures in the north western (NW) part of the shell (Bamba et al. 2005). These filaments are usually explained by fast synchrotron cooling in a strongly amplified magnetic field (B ≿ 100 μG; see e.g., Berezhko et al. 2009). This strong magnetic field is in contradiction with the leptonic scenario of the gamma-ray emission, as in this case a relatively weak field of ∼10 μG is required to explain the overall X-ray emission from the remnant (H.E.S.S. Collaboration 2018; Lee et al. 2013). However, it is also possible that X-ray filaments are limited by magnetic field damping downstream of the forward shock (Pohl et al. 2005; Rettig & Pohl 2012). In this scenario the magnetic field downstream of the shock can be weaker as the width of the filament is determined by the length scale at which strong magnetic turbulence formed upstream is damped downstream.

Recently a detailed analysis of the archival XMM-Newton data revealed a gradual softening of the spectrum (from a photon index of 2.56 to a photon index of 2.96) in the north-western rim towards the interior of the remnant (Kishishita et al. 2013). This softening was interpreted as the decrease of the maximum electron energy with time due to synchrotron cooling. Using a simple spectral-evolution model assuming that the downstream time evolution of the maximum electron energy is determined only by synchrotron losses, Kishishita et al. (2013) found that the softening of the spectrum can be explained by a rather low magnetic field of ∼10 μG. However, this requires an electron maximum energy of about 50 TeV which is higher than that suggested by H.E.S.S. observations at TeV energies (H.E.S.S. Collaboration 2018). More importantly, these calculations were done neglecting a projection effect which can potentially strongly weaken the spatial variation of the spectrum depending on the three-dimensional (3D) structure of the remnant. It can also be shown that such a gradual softening of the spectrum can be similarly explained by magnetic field damping. Indeed, the cut-off frequency of the synchrotron spectrum decreases linearly with the decrease of the magnetic field strength, and therefore the damping of the magnetic field downstream of the shock can potentially result in similar spectral variation to that seen during synchrotron cooling. Therefore, in this work we test both scenarios, properly accounting for the projection effect.

Recent observations with Suzaku (Takeda et al. 2016) also reveal a faint hard X-ray emission from the NW rim of the remnant. The spectrum in the energy range from 12 to 22 keV can be described with a power law with a spectral index of . This result challenges the existence of a roll-off in the X-ray spectrum or at least suggests that the cut off of the electron spectrum is slower than exponential. However, the uncertainties on the spectral parameters are too large to draw strong conclusions.

It should also be mentioned that the southern part of the shell is coincident with a pulsar wind nebula (PWN) detected in X-rays with XMM-Newton (Acero et al. 2013). H.E.S.S. observations provide a hint that a fraction of up to 8% of the total VHE emission from the Vela Jr. region might be associated with this PWN (H.E.S.S. Collaboration 2018). If confirmed in future observations, for example with CTA, this should be taken into account for the broadband modeling of the emission.

In this paper we provide modeling of the nonthermal emission from the Vela Jr. SNR, focusing on an explanation of the morphology of the remnant and the constraints the observed morphology imposes. In this regard we consider two different profiles of the magnetic field downstream of the shock, damped magnetic field and transported magnetic field, to test which of the two effects described above determines X-ray filaments and spectral softening observed in these filaments. For our simulation we use the Radiation Acceleration Transport Parallel Code (RATPaC) described in a number of previous papers (Telezhinsky et al. 2012, 2013; Brose et al. 2016). The code solves the time-dependent transport equation for cosmic rays in one dimension (1D) in the test particle regime and subsequently simulates the nonthermal radiation from accelerated particles. The shock evolution is simulated using 1D hydrodynamic simulations performed by the Pluto code (Mignone et al. 2007) which is the most recent addition to RATPaC. The calculation of the particle acceleration can also be coupled to the evolution of isotropic Alfvénic turbulence upstream of the shock (Brose et al. 2016).

2. Modeling

2.1. Physical parameters of the SNR and its environment

We assume that Vela Jr. was created in a core collapse supernova (SN) explosion, which is supported by the detection of the central compact object (CCO) AX J0851.9−4617.41 close to the center of the remnant (Aschenbach 1998; Aschenbach et al. 1999; Slane et al. 2001). However, it is not clear whether the association with the SNR is correct because no pulsations were detected and later it was also suggested that a potential CCO might be an unrelated planetary nebula (Reynoso et al. 2006). Additionally, the absence of broad Ca II absorption lines in spectra of background stars together with the constraints set by the Ti44 gamma-ray line also suggests a core collapse origin with SNe types of Ic or Ibc (Iyudin et al. 2010). In our hydrodynamic simulations (see below) we assume an explosion energy of 1051 erg and an ejecta mass of 3 M which roughly corresponds to the properties of a type Ibc SN.

The current density of the ambient medium is constrained by the lack of the X-ray thermal emission which places the upper limit at , where D is the distance to the remnant and f is the volume-filling factor of the X-ray-emitting gas (Slane et al. 2001). That limit applies to gas with a temperature above 1 keV, because at energies below 1 keV the observed X-ray emission is dominated by the strong thermal emission from the Vela SNR.

The distance to the remnant and its age can be constrained by its angular size and its angular expansion rate measured by Katsuda et al. (2008) and Allen et al. (2015) using XMM-Newton and Chandra data, respectively. XMM-Newton data taken between 2001 and 2007 imply an expansion rate of 0.84 ± 0.23 yr−1 (Katsuda et al. 2008) while Chandra data from 2003 to 2008 suggest an expansion rate of 0.42 ± 0.10 yr−1 (Allen et al. 2015). The difference between the two measurements is somewhat puzzling because the data used to estimate the expansion rate were taken from almost the same region and over almost the same time period. In our simulations we use the measurement provided by XMM-Newton to constrain the hydrodynamic model but in Sect. 4.3 we also discuss how uncertainties on the expansion rate measurement can influence our results. The allowed range for the distance is 500−1000 pc. The lower limit is determined based on the requirement that the shock speed should be higher than 1000 km s−1 to effectively accelerate particles (Allen et al. 2015). The upper limit of 1 kpc is determined by the comparison of the detected 44Ti line emission with evolution models for different types of SN explosions (Iyudin et al. 2010). The comparison of the X-ray absorption column density to the distribution of 12CO and HI gas places Vela Jr. in the foreground or at the interface of the Vela molecular ridge suggesting a similar upper limit of 900 pc (Acero et al. 2013).

2.2. Hydrodynamics

The standard gas-dynamical equations,

(1)

(2)

are solved, where ρ is the density of the thermal gas, v the plasma velocity, m = the momentum density, P the thermal pressure of the gas, I the unit tensor, and E the total energy of the ideal gas with γ = 5/3. We assume that the magnetic field is dynamically unimportant due to its low strength and the remnant not being in the radiative phase yet (Petruk et al. 2016). The equations are solved in 1D for a spherical symmetry.

We initialize the ejecta profile by a plateau in density with the value ρc up to the radius rc followed by a power-law distribution up to the ejecta-radius Rej:

(3)

The exponent for the ejecta profile is set to n = 9 for the core-collapse explosion. The velocity of the ejecta is defined as

(4)

where TSN = 1 yr is the initial time set for hydrodynamic simulations. Then for the assumed mass, Mej, and energy, Eej, of the ejecta, and defining the radius of the ejecta as multiple of rc, Rej = xrc, the initial conditions for simulations can be written as

(5)

(6)

(7)

The initial temperature is set to 104 K everywhere and the initial pressure is calculated using the equation of state. We used 200 000 linearly distributed grid-points and 14 pc grid-size for the hydro simulations. The heaviest ejecta bin contains about 0.1M. For our simulations we set x = 3 resulting in rc = 0.006 pc and Rej = 0.018 pc. To use the hydro-output for the particle acceleration, we had to resharpen the shock. We further assume in Eq. (3) that the SNR is expanding into a wind zone, i.e., the density radial distribution follows nr−2, and require the current upstream density to be lower than 0.03 cm−3. Simulations are stopped at the moment when the expansion rate of the remnant, that is, the ratio of the shock speed to the distance, reaches the observed value (see the previous section).

Table 1 summarizes the properties of the progenitor used as input parameters to our hydrodynamic model and lists the current properties of the SNR obtained from hydrodynamic simulations. The mass-loss rate of the progenitor star, , is a free parameter, which is tuned to get the right combination of the current wind density and expansion rate of the remnant. The choice of is, however, in agreement with that expected for red supergiants, and it leads to a total mass in the wind of about 13 M and a time of about 900 kyr to build the wind zone. The time evolution of the shock radius and shock velocity is shown in Fig. 1.

Table 1.

Parameters of the hydrodynamic model.

thumbnail Fig. 1.

Time evolution of the forward shock radius (top panel) and shock velocity (bottom panel) up to the current age as obtained from our hydrodynamic model.

Open with DEXTER

2.3. Magnetic field

We assume that the magnetic field in the wind zone of the progenitor star, Bwind, follows a 1/r dependence on the radial distance. It should be noted that the fit of the observational data for that radial dependence of the circumstellar magnetic field results in similar values of the current magnetic field strength as compared to a magnetic field of constant strength. However, stronger synchrotron cooling at earlier stages of the SNR evolution might significantly modify the resulting particle spectrum in the case of the 1/r profile.

We also assume that the circumstellar magnetic field, Bwind, is amplified upstream of the shock by streaming instabilities. The magnetic field strength in the immediate upstream region is then given by Bu = kBwind. Subsequently, for the shock compression ratio of 4, the immediate downstream magnetic field strength is . It is assumed here that the magnetic field upstream of the shock is turbulent and fully isotropised. In addition to the magnetic turbulence created by streaming instabilities, we expect pre-existing turbulence in the stellar wind (see Bruno & Carbone 2013, for the review of the turbulence in the solar wind).

We apply two different parametrizations to describe the distribution of the magnetic field downstream of the shock. In the first scenario we assume that turbulently amplified magnetic field is effectively damped downstream of the forward shock of the remnant. The magnetic turbulence created upstream and at the shock by various instabilities is transferred downstream where it is eventually damped due to the lack of the turbulence driving (Pohl et al. 2005). To describe the magnetic field profile downstream of the shock in this case we adopt a simple parametrization (Pohl et al. 2005),

(8)

where B0 is the magnetic field far downstream of the shock, and is set in our simulations to be equal to Bwind, and ld is the damping length scale. The value of the damping length scale is chosen to be ld = 0.05Rsh (∼0.6 pc) to fit the peak of the X-ray brightness profile. This large damping scale is consistent with the low magnetic field strength in Vela Jr. with correspondingly large eddy size. The damping length can be estimated based on the eddy turnover time at the largest scales:

(9)

where χ is a nonlinear efficiency factor of the order of a few, and vd is the downstream flow speed. We can identify the wavelength, λ, with the Larmor radius of cosmic rays at the highest energy attained in the remnant, λrL(Emax). For the velocity fluctuations we assume energy equipartition between magnetic and kinetic fluctuations,

(10)

Combining all expressions, we find for the damping scale

(11)

where E100 = Emax/(100 TeV), n−1 = n−d/(0.1cm−3), and B−5 = Bd/(10 μG). Noting that all quantities here are measured in the downstream region, except for the shock speed, vsh, we find with the best-fit parameter values of our model that ld ≃ 0.6 pc for a reasonable χ ≃ 7.

Alternatively, we assume that the immediate downstream magnetic field Bd is transported inside the SNR with the plasma flow and evolves following the induction equation for ideal magnetohydrodynamics.

Parameters which describe the magnetic field in both scenarios are collected in Table 2. Figure 2 illustrates magnetic field profiles in both cases. The parameters are chosen in order to fit the broadband photon spectrum (see Sect. 3).

Table 2.

Models description.

thumbnail Fig. 2.

Magnetic field strength as a function of the normalized radius between the contact discontinuity and the forward shock. Red solid line shows the profile for the magnetic damping scenario and the blue dashed line for the transported magnetic field. Both distributions are shown for the models with the Bohm diffusion coefficient.

Open with DEXTER

The strength of the magnetic field in the wind at the current stage of evolution adopted in our models is approximately 1 μG (see Table 2). This value is in agreement with the experimental measurements of the surface magnetic field of red supergiants of 1−10 G (see e.g., Tessore et al. 2017). The assumed 1/r profile yields a surface magnetic field of B ≃ (6 G) (R/(100 R))−1.

2.4. Particle acceleration

We simulate the particle density evolution solving the transport equation in the form

(12)

where N is the differential number density of cosmic rays, D is the spatial diffusion coefficient, v is the plasma velocity, represents the energy losses (in our case synchrotron losses) and Q is the source term.

We use the general form of the spatial diffusion coefficient

(13)

and consider two different cases: Bohm diffusion with ζ = 1 and α = 0 and Kolmogorov-type diffusion with α = −2/3 and ζ chosen in a way to fit the broadband photon spectrum. The type of the diffusion determines the shape of the cut-off in the resulting particle spectrum: Kolmogorov-type diffusion implies a slower cut-off than Bohm diffusion. We also calculate the diffusion coefficient by solving the transport equation for magnetic turbulence dominated by Alfvén waves (Brose et al. 2016). The equation is solved in 1D for spherical symmetry assuming that Alfvénic turbulence is isotropic and accounts for compression, advection, cascading, damping and growth due to resonant amplification of Alfvén waves. The shape of the cut off of the resulting spectrum of accelerated particles falls in between Kolmogorov-type and Bohm diffusion. In reality, Alfvén waves are not the only waves in the magnetic turbulence, which might also modify the resulting spectrum. Therefore, two parametrizations considered for the diffusion coefficient can be considered as extreme cases defining the range of possible scenarios.

Particle injection is determined by the source term which is given by

(14)

where η is the injection efficiency parameter, nu the plasma number density in the upstream region, Vsh the shock speed, Vwind the wind velocity upstream of the shock, Rsh the shock radius and pinj = ξpth is the injection momentum, defined as a multiple of the momentum at the thermal peak of the Maxwellian distribution in the downstream plasma with temperature Td, . We assume the thermal leakage injection model (Blasi et al. 2005), where pinj is the minimum momentum for which a thermal particle can cross the shock and enter the acceleration process and the injection efficiency for the compression ratio of 4 is determined as

(15)

In our simulations we inject all the particles at the same temperature at the position of the shock with the momentum pinj assuming that the injection parameter ξ is the same for electrons and for protons. In this case the resulting electron to proton ratio at high energies is determined by the mass ratio, (Pohl 1993). The injection parameter ξ is a free parameter which is determined by the fit of the expected emission to the observed gamma-ray data.

We solve the transport equation in 1D in the test particle regime separately for electrons and protons using the RATPaC code as described in Telezhinsky et al. (2012, 2013) taking into account only a forward shock and ignoring a reverse shock as its contribution to the overall particle spectrum is negligible due to the low density in the upstream of the reverse shock at the current age. The resulting particle spectrum at the current age of the remnant is then used to simulate nonthermal emission. The particle spectrum at low energies follows a power law with a spectral index of 2 because of the test-particle approximation. The choice of the test-particle approximation is well justified as the ratio of the cosmic ray pressure to the ram pressure of the incoming flow, Pcr/Pram, does not exceed 2% in any of our models (see Tables 2 and 4).

2.5. Nonthermal emission

We calculate nonthermal emission from the remnant considering three emission processes: synchrotron radiation, IC scattering of accelerated electrons on the cosmic microwave background (CMB) and proton interactions with subsequent pion decay. For IC scattering we do not consider other possible radiation fields as the Vela Jr. SNR is located far away from the Galactic center, and local optical and infrared radiation fields are weak and should not play a strong role. Synchrotron emission is calculated accounting for the turbulent magnetic field (Pohl et al. 2015).

3. Results

In this section we present the results of our modeling of the broadband emission from the Vela Jr. SNR. The emission from the remnant is calculated assuming spherical symmetry. We construct four models based on our assumptions of the form of the spatial diffusion coefficient and the shape of the magnetic field distribution downstream of the shock. Main parameters used in these models are collected in the Table 2.

3.1. Broadband photon spectrum

Figure 3 shows the observed broadband spectral energy distribution (SED) overlaid with the curves representing simulated emission in four different models. In all four models the simulated SED strongly underpredicts the observed radio emission. This discrepancy is further discussed in Sect. 4.1.

thumbnail Fig. 3.

SED of the Vela Jr. SNR. Lines represent simulated emission from the source produced in different processes as specified in the legend in the magnetic field damping model (left panels) and synchrotron cooling model (right panels). Green triangles indicate the observed radio emission as detected by Parkes (Duncan & Green 2000), an orange bow-tie shows the spectral fit of the X-ray XMM-Newton data (Aharonian et al. 2007), red squares represent the Fermi-LAT data points (Tanaka et al. 2011) and blue circles show the H.E.S.S. data points (H.E.S.S. Collaboration 2018).

Open with DEXTER

The choice of the spatial diffusion coefficient determines the shape of the cut-off in the electron spectrum and therefore the shape of the X-ray and gamma-ray spectra. In the case of the Bohm diffusion we fail to reproduce the slope of the X-ray spectrum as the synchrotron emission exhibits a faster-than-observed cutoff of the spectrum (Fig. 4, top panel). Our BOHM models also slightly underpredict GeV emission but at the same time provide a very good fit of the TeV spectrum (Fig. 4, bottom panel). On the other hand, the models with the Kolmogorov-type diffusion coefficient much better reproduce the observed X-ray slope and GeV data, but yield a slightly worse fit of the cut-off at VHEs, although still compatible with data. The observational data suggest a faster cut-off in the gamma-ray spectrum with respect to the X-ray spectrum. In addition, the detection of hard X-rays from the remnant with Suzaku (Takeda et al. 2016) also suggests a much slower cutoff in the X-ray spectrum which is much better reproduced by the Kolmogorov-type diffusion coefficient. The observed GeV data suggest a softer electron spectrum than the one obtained in DSA in the test-particle approximation. However, it is still possible to reproduce the observed GeV spectrum with the Kolmogorov-type diffusion coefficient due to the lower cutoff energy and slower cutoff.

thumbnail Fig. 4.

Zoom into X-ray (top panel) and GeV–TeV (bottom panel) spectra for all four models.

Open with DEXTER

The IC emission in the BOHM_BTRAN model shows a shoulder around 1 GeV due to synchrotron cooling. This feature is created because of the 1/r profile of the ambient magnetic field. At the early stages of SNR evolution the ambient magnetic field is very strong, cooling electrons efficiently. This feature is not particularly prominent in the synchrotron spectrum because the cooled particles responsible for this break are accumulated at the contact discontinuity where the magnetic field is lower and the synchrotron emission is therefore weaker; it is also not visible in KOLM_BTRAN model due to a weaker ambient magnetic field. The impact of the ambient magnetic field on the broadband radiation spectrum will be further discussed in a follow-up paper.

The contribution of the pion decay process to the overall gamma-ray flux is negligible for the type of injection we implement. If we relax the requirement that the injection parameter is the same for electrons and protons, and reduce the electron-to-proton ratio, the observed SED could be reproduced with a higher contribution from the pion decay in the gamma-ray energy range. However, it is easy to see that matching the level of the GeV–TeV emission exclusively by the pion decay process would require a total energy in protons comparable to or higher than the supernova explosion energy, as this would require ∼100 times more protons.

The difference between the models with damped and transported magnetic fields is not particularly pronounced in the total volume integrated spectral distribution of the SNR. However, as is shown in the following subsection, significant differences arise in the morphology of the remnant and spatial spectral variations. It should be noted though, that the fit of the observed SED requires different values of the magnetic field strength downstream of the forward shock in these two scenarios (see Table 2). Damping models require a stronger field in the immediate downstream region.

Our results are sensitive to variations of the model parameters. Although small changes do not strongly modify the models, the parameter space is well constrained by the observational data. The X-ray to gamma-ray flux ratio constrains the magnetic field strength which in the immediate downstream cannot be higher than ∼15 μG, even in the damped field scenario. Magnetic field strength in turn determines the total energy in accelerated electrons. We also notice that the maximum energy to which particles can be accelerated is sensitive to the magnetic field strength in the far upstream of the shock as it determines the diffusion of particles away from the shock. The sensitivity of our results on the parameters of the hydrodynamic model is discussed in Sect. 4.3.

3.2. X-ray morphology and spatial spectral variations

To study the morphology of the remnant we produce 2D energy-dependent intensity maps with a cell size of 0.67, equal to the width of regions 1–4 as defined for the X-ray spectral analysis in Kishishita et al. (2013). Then for the cells falling into these regions we extract the flux and the spectrum in the energy range between 2 and 10 keV. Regions 5–9 are twice as wide, and for these regions we combine the emission from two cells for the spectrum estimation.

Figure 5 compares the simulated surface brightness profile to the observed one. For a direct comparison, the observed flux (Kishishita et al. 2013) was transformed into the surface brightness accounting for the different size of the regions. All profiles are normalized in a way that the integral along the radius of the remnant above 0.85 Rsh (above 51 for the radius of 1) is equal to unity. The difference between Bohm and Kolmogorov-type models is almost negligible, but the shape of the brightness profile strongly depends on the magnetic field model. The models with the magnetic damping scenario (red markers) better reproduce the shape of the profile with the emission peaking around 1 away from the shock. In the case of the transported magnetic field (blue markers) the peak is located closer to the center of the remnant at the angular distance of around (3−5) from the shock. Neither case, however, can reproduce the magnitude of the flux decrease. While in the case of the transported field scenario the magnetic field stays almost constant downstream of the shock, resulting in substantial synchrotron emission from the interior, in the damped field scenario the magnetic field strength decreases very quickly with the distance from the shock, suppressing synchrotron emission; the plateau in the profile is therefore mainly determined by the projection effect and inevitable under spherical symmetry. Therefore, the observed decrease of the X-ray intensity towards the center of the remnant could indicate a strong deviation from spherical symmetry. This point is further discussed in the Sect. 4.2.

thumbnail Fig. 5.

Simulated projected X-ray radial flux profiles for different models assuming spherical symmetry compared to the experimental measurements (green squares Kishishita et al. 2013). Normalized surface brightness in the energy range 2−10 keV is plotted as a function of the angular distance from the forward shock of the SNR with negative values corresponding to the interior of the remnant. The radial profile obtained within the asymmetric model discussed in Sect. 4.2 is represented by the black filled circles connected with the dashed line.

Open with DEXTER

To study the spectral variations with the distance from the shock we extract the spectrum for each region and fit it in the same way as Kishishita et al. (2013). Simulated spatial distributions of the spectral index and of the cut-off energy in all four models follow the same trend as the observed radial distributions. However the observed variation of the spectral parameters is stronger. Close to the shock the spectral variation is stronger in the case of the damped magnetic field profile, but at the distance ≿3 the spectral shape stays roughly constant while in the case of the transported field the spectrum gets continuously softer with distance from the shock. This happens because in the case of the damped profile the magnetic field strength decreases much faster with distance from the shock and at ∼3the spectrum is already fully determined by the projection effect. In the case of the transported field, the magnetic field strength drops more slowly and the interior of the remnant still significantly contributes to the integrated intensity along the line of sight.

Distributions resulting from models using Bohm and Kolmogorov-type diffusion coefficients appear to mimic each other, but differ in the parameter values. Kolmogorov-type diffusion much better reproduces the spectral shape close to the forward shock where most of the X-ray emission is coming from. This is not surprising as the Kolmogorov-type models better fit the spectrum of the total X-ray emission, and the total emission is dominated by the emission from the shock region.

thumbnail Fig. 6.

Simulated projected radial profiles of the X-ray spectral parameters for four different models assuming spherical symmetry. Top panel: cut-off energy and the bottom one shows the spectral index in the 2−10 keV band as a function of the angular distance from the forward shock of the SNR. The distributions are obtained by extracting the spectrum for each region and fitting it with an exponential cut-off power law assuming a photon spectral index of Γ = 1.6(top panel) and with a simple power law (bottom panel). In both panels, green squares indicate the observational data (Kishishita et al. 2013).

Open with DEXTER

3.3. Gamma-ray morphology

In Fig. 7 we compare our simulated gamma-ray profiles to the TeV profiles measured by H.E.S.S. Simulated profiles in the energy range from 300 GeV to 20 TeV are extracted from the simulated surface brightness map smoothed with a Gaussian kernel with σ = 0.06, which roughly resembles the H.E.S.S. PSF in the Aharonian et al. (2007) analysis. The H.E.S.S. profile was normalized so that the sum of the histogram contents between 0.3 and 1.2 is equal to unity (Aharonian et al. 2007). Simulated profiles were normalized in the same way for the direct comparison. We show the profiles only for the models with Bohm diffusion coefficient (red solid and blue dashed lines; green dash-dotted line represents the asymmetric model discussed in Sect. 4.2) as the ones for the Kolmogorov-type diffusion are very similar.

thumbnail Fig. 7.

Simulated radial surface brightness profiles of the gamma-ray emission in the 300 GeV−20 TeV energy range compared to the H.E.S.S. measurements for the northern part of the remnant (Aharonian et al. 2007, black filled circles). Red solid and blue dashed lines show simulated profiles for the spherically symmetric models BOHM_DAMP and BOHM_KOLM, respectively. The radial profile obtained within the asymmetric model discussed in Sect. 4.2 is represented by the green dash-dotted line. This profile is averaged over the azimuth angle for the hemisphere containing the cone of enhanced emission. The thick gray line shows the azimuthally averaged radial profile for the northern part of the SNR extracted from the recently published updated H.E.S.S. excess skymap (H.E.S.S. Collaboration 2018). This skymap is smoothed with a Gaussian function of width 0.08. All simulated profiles are normalized in the same way as H.E.S.S. measurements in Aharonian et al. (2007).

Open with DEXTER

It is clear that the model with the transported magnetic field resembles the data much more closely than the magnetic damping model. In the magnetic damping scenario, the peak appears to be significantly broader and the decrease of the emission towards the center of the remnant is slower. This is due to a much weaker synchrotron cooling of the particles downstream of the shock at distances r/RFS ≾ 0.9. This results in a larger number of high-energy particles inside the remnant and in turn a broader radial distribution.

This result is in conflict with the result for the X-ray morphology which is better reproduced by the magnetic damping scenario. For the synchrotron emission in the magnetic damping scenario the key parameter is the magnetic field strength, which drops exponentially with the distance from the shock and thus determines the shape of the radial profile.

4. Discussion

4.1. Radio emission

All four models strongly underpredict the radio emission from Vela Jr. This could be an indication that our models do not account for some other possible processes taking place in SNRs. For example, stochastic re-acceleration of particles at the turbulence in the immediate downstream region of the shock (Pohl et al. 2015) can soften the resulting synchrotron spectrum, which means that it would be possible to explain the radio emission within the models presented here by simply adding another process.

Another possible approach to reproducing the radio emission from Vela Jr. would be to account for the nonlinear feedback of the cosmic-ray pressure on the hydrodynamic evolution of the remnant. In this case, cosmic rays would modify the shock and in turn the particle spectrum, resulting in a softer spectrum. For example, Lee et al. (2013) succeeded in reproducing the radio emission together with the X-ray emission and gamma-rays using the nonlinear DSA calculations with a semianalytic solution (e.g., Blasi et al. 2005; Caprioli et al. 2009). A drawback of their model is that they do not simulate the shape of the cut-off of the particle spectrum but simply assume it to follow ∝ exp(−(p/pmax)β), where pmax is the maximum momentum of the particle and is obtained from simulations, and β is simply chosen in such a way that the model explains the X-ray spectrum. Therefore, in this approach the model parameters can be tuned in a way that makes it possible to obtain the pmax required to fit the radio points even for the particle spectral index of 2.

Another effect, which has not been considered here and can cause a softening of the particle spectrum, is Alfvénic drift (see e.g., Morlino & Caprioli 2012). However, the magnetic field suggested by the simultaneous fit of the X-ray and gamma-ray emission is too low to cause this effect, because the Alfvén velocity is VA ∼ 100 km s−1. For a higher magnetic field strength, radio and X-ray spectra could possibly be fit even without Alfvénic drift, as in this case strong synchrotron cooling would create a break in the synchrotron spectrum. However, stronger magnetic field would require fewer electrons to produce the X-ray emission and therefore the model would strongly underpredict the gamma-ray emission.

It should also be noted that the reported value of the detected radio flux can be misleading as the remnant is located in a highly crowded region with a very complicated background and there are indications that a part of the radio continuum emission coincident with Vela Jr. might be associated with the Vela SNR (Maxted et al. 2018).

High-resolution radio observations and detailed morphological studies could be of great importance for constraining the magnetic field downstream of the shock. Damped and transported magnetic field profiles predict different radio morphology of the remnant. In the case of the damped field, similarly to the X-ray emission, the radio emission from the SNR would be determined by the strength of the magnetic field and therefore the radio morphology of the remnant should be similar to the X-ray morphology. In the transported field scenario the magnetic field strength is almost constant downstream (see Fig. 2), meaning that the morphology is determined mainly by the number of particles emitting synchrotron radiation. The X-ray emission would still peak close to the shock because of the lack of high-energy electrons far downstream due to synchrotron cooling, but the radio emission is expected to have its maximum closer to the contact discontinuity due to a high abundance of low-energy electrons combined with the projection effect. Radio images obtained with the 64 m Parkes radio telescope (Duncan & Green 2000) show a good correlation of the radio emission with the X-ray emission, suggesting that the magnetic field should quickly decrease downstream. A more quantitative comparison is difficult, however, due to the poor angular resolution of radio observations and possible contamination by the emission from the Vela SNR.

4.2. Deviation from spherical symmetry

The observed X-ray flux profile exhibits a very fast decrease of the flux towards the interior of the remnant. This is a strong indication for a significant deviation from spherical symmetry. Indeed, even if we assume that the magnetic field strength falls to zero and that electrons capable of radiating X-ray synchrotron can be found only in the immediate vicinity of the shock, it would still be impossible to obtain such a decrease within the assumption of spherical symmetry simply due to the projection effect. The projection effect also causes a slower softening of the spectrum towards the interior as the projected emission from regions close to the shock dominates the total projected emission. Therefore, the observed softening of the spectrum which is faster than the one simulated within the spherical symmetry assumption also suggests the asymmetry of the remnant. Finally, the observed X-ray morphology of the remnant also deviates from spherical symmetry showing a number of regions with local brightening of the X-ray emission, one of which was used by Kishishita et al. (2013) to detect the spectral softening.

To test the asymmetry hypothesis and check how well we can reproduce the morphology, we assume that in one radial direction within the solid angle Ω the synchrotron emission is brighter than in the rest of the remnant. This could be achieved, for example, if in this cone the injection of particles was more efficient due to higher ambient density in this region. In this case, the radial profile taken in this direction would much more closely resemble the data as the contribution of the shock region emission to the projected emission inside the remnant would be lower. We solve the transport equation in 1D assuming there is no tangential transport of particles and therefore we can simply modify one of our models calculated for spherical symmetry by assuming that synchrotron emission is stronger in the region limited by the cone with the solid angle Ω. For this exercise we use model KOLM_BDAMP as it reproduces both the right position of the peak in the X-ray flux profile and the X-ray spectral index value in the immediate downstream of the shock. We assume that the synchrotron emission within the solid angle Ω = 2π(1−cos θ) is stronger than in the rest of the remnant by a factor κ, but at the same time keep the total synchrotron emission from the remnant constant, so that our modified model still reproduces the SED in the same way as KOLM_BDAMP. The opening angle θ is set by the size of the bright filament, which roughly corresponds to the size of the field of view studied in Kishishita et al. (2013).

The simulated radial profile of the filament for θ = 15 and κ = 12 is shown with black filled circles connected with the dashed line on Fig. 5 compared to the observed radial profile. The regions are defined in the same way as for the spherically symmetric models. It can be seen that the modified asymmetric model can reproduce the observed radial profile fairly well, reproducing both the shape of the profile and the magnitude of the flux decrease.

As expected, by improving the fit of the surface brightness profile, we can also better reproduce the observed softening of the spectrum. Figure 8 shows simulated radial distributions of the spectral index and the cut-off energy compared to the observational data. Simulated distributions start to deviate from the measured values at ≿5 from the forward shock. This happens because at this distance the contribution from the cone of the brighter emission becomes negligible in comparison to the rest of the remnant along the line of sight and therefore simulated spectral parameters start to resemble the spherically symmetric model. This results in an upturn for the cut-off energy and downturn for the spectral index at ∼5 which is not seen in the spherically symmetric model because the emission is evenly distributed.

thumbnail Fig. 8.

Simulated cut-off energy (top panel) and spectral index (bottom panel) of the X-ray spectrum in the filament as a function of the angular distance from the shock compared to the experimental measurements. Simulated distributions are obtained using the asymmetric model described in Sect. 4.2.

Open with DEXTER

If the brightening of the synchrotron emission is due to the higher number of electrons in the region limited by the cone, this would also naturally enhance the inverse Compton emission from that region by the same factor. Observationally, this would result in a local brightening of the gamma-ray emission from the location of the X-ray filament which is in agreement with the observed correlation between the TeV and X-ray morphology. Remarkably, the choice of the cone parameters (the opening angle θ and the brightening factor κ), motivated by the size and the properties of the X-ray filament, also results in a much better fit of the observed azimuthally averaged surface brightness profile of the gamma-ray emission above 300 GeV as compared to the model with the damped magnetic field under the spherical symmetry assumption (Fig. 7). This means that for the asymmetric case the damped field model can fit both X-ray and gamma-ray profiles as opposed to the spherically symmetric case where the gamma-ray intensity profile prefers the transported field model.

Our cone model for the X-ray filament appears to also be in agreement with the azimuthal intensity distribution of the X-ray and TeV emission. H.E.S.S. Collaboration (2018) provide azimuthal intensity profiles for an annulus with an inner radius of 0.6 and outer radius of 1 around the center of the SNR split into forty equally sized bins. Defining the bins in the same way we extract the flux from the bin containing the bright filament and compare it to the flux from the bin outside the cone of the brighter emission. The azimuthal intensity contrast estimated as the ratio of these two fluxes is about five (for both X-ray and TeV emission) which is in agreement with the observed ratio.

As mentioned above, Fukui et al. (2017) detected a good spatial correlation of the TeV emission with the distribution of the HI gas, indirectly supporting the hadronic scenario for the gamma-ray emission because of the higher density of the ambient medium. This, however, also implies more efficient particle injection and thus might explain local brightening of X-ray and TeV emission within the leptonic scenario simply because of the higher number of accelerated electrons. At the same time, gamma-ray emission produced in hadronic interactions might still be too low to significantly contribute to the observed gamma-ray emission. Therefore the HI gas does not necessarily indicate the hadronic scenario, but it might be the reason for the existence of regions with local brightening of the emission. This idea is also supported by the coincidence of the HI emission with a radio filament as noticed by Maxted et al. (2018).

4.3. Uncertainties surrounding the hydrodynamic model

To test how sensitive our results are to the parameters of our hydrodynamic model we built an alternative model based on the expansion rate measured using the Chandra data (Allen et al. 2015). We kept the same ejecta mass, Mej = 3 M, and wind speed, vwind = 15 km s−1, and varied all other parameters to obtain the right angular size of the remnant and to reproduce the observed expansion rate. As before, we constrain the ambient density by the lack of X-ray thermal emission and set the upper limit on the distance at 1 kpc. The parameters of the alternative model are listed in Table 3. A lower expansion rate naturally yields a lower shock speed and higher age of the remnant. To keep the angular size unchanged we also need to slightly increase the density, the distance to the remnant, and the mass-loss rate of the progenitors wind. A higher mass-loss rate, however, implies a rather high total mass in the wind of about 33 M.

Table 3.

Parameters of the alternative hydrodynamic model.

Based on the new hydrodynamic model, we create a modified model of the broadband emission for the Bohm diffusion coefficient and damped magnetic field. The model parameters are listed in Table 4. We tune the parameters so that the model provides a fit to the data that is similar to that of the primary model. It appears that the fit requires only small changes of the main parameters (Table 4).

Table 4.

Description of the emission model for the alternative hydrodynamic solution.

5. Conclusions

In this work we studied the Vela Jr. SNR through the modeling of its broadband nonthermal emission aiming to explain the observed SED together with spatial variations of intensity and spectral parameters. One of the main questions that we try to address is the nature of observed X-ray filaments which feature a fast decrease in intensity and strong spectral softening towards the interior of the remnant. These filaments are usually explained by fast synchrotron cooling in a strongly amplified magnetic field of ∼100 μG, but can also be caused by the fast decrease of the magnetic field strength downstream of the shock due to the magnetic field damping. The latter effect does not require as strong a magnetic field downstream of the shock and therefore might be a more plausible explanation for the X-ray filaments observed in the Vela Jr. SNR. Indeed, our simulations favor the leptonic scenario, suggesting that the observed gamma-ray emission is produced predominantly via IC scattering of relativistic electrons accelerated at the shock of the remnant, implying that the downstream magnetic field is 8−13 μG. The contribution of the gamma-ray emission produced in hadronic interactions appears to be negligible. We consider two different descriptions for the magnetic field profile downstream of the shock, for the damped and transported magnetic fields. In both cases the softening of the X-ray spectrum as well as the X-ray intensity profile is determined by the decrease of the magnetic field strength with the distance from the shock. The synchrotron cooling plays a less important role due to the low value of the magnetic field strength. The projection effect in the spherically symmetric case smears the spatial variation of both intensity and spectral parameters resulting in a similar trend but much weaker variation. The observed depth of the decrease of the X-ray intensity cannot be reproduced even for the ideal case when all the emission is coming from the location of the shock, strongly indicating a significant deviation from spherical symmetry. We constructed an asymmetric model for the damped field scenario by assuming that the X-ray emission is enhanced within a cone with an opening angle of θ = 15, roughly corresponding to the size of the observed X-ray filament. It is interesting that in this scenario both the intensity profile and the spectral variation can be explained quite well. The reason for the enhanced emission inside this cone could be more efficient particle injection due to the higher ambient density. This idea seems to be in agreement with the recently detected correlation of the HI gas distribution with bright TeV emission regions which also correspond to local brightenings at X-ray energies. Moreover, if the enhanced X-ray emission is really caused by the more efficient electron injection, this would also cause the enhancement of the TeV emission in a similar fashion. Remarkably, the choice of the cone parameters tuned for the fit of the X-ray properties also results in a much better fit of the azimuthally averaged TeV intensity profile as compared to the spherically symmetric model for the damped field. As a result, for the asymmetric case, both X-ray and TeV intensity profiles can be reproduced within the damped field scenario. The choice of the cone parameters is also in perfect agreement with the observed azimuthal profiles at X-ray and TeV energies.

In our simulations we use two different forms of the spatial diffusion coefficient (Bohm diffusion and Kolmogorov-type diffusion) which determines the shape of the cut-off of the particle spectrum. These two parametrizations can be considered as extreme cases defining the range of possible scenarios. Independent simulations including the solution of the transport equation for magnetic turbulence driven by Alfvén waves result in a particle spectrum with a cut-off that is slower than for the Bohm diffusion but faster than for the Kolmogorov-type diffusion. The detected TeV and X-ray spectra favor different types of diffusion. The X-ray spectrum is better reproduced by the Kolmogorov-type diffusion. The Bohm diffusion implies a much faster cut-off which is also in contradiction with the recent detection of hard X-rays from the northern rim of the remnant. However, the shape of the cut-off at TeV energies slightly favors the Bohm diffusion.

All our models strongly underpredict radio emission from the remnant. This might indicate that we do not account for some other processes taking place in the SNR and responsible for the radio emission. One of the possibilities could be stochastic reacceleration of electrons at the turbulence in the immediate downstream region of the shock. However, it is also possible that a part of the radio emission coincident with Vela Jr. might be associated with the Vela SNR located in the foreground. If that were the case, the real radio flux from the remnant would be lower.

Further theoretical studies of the Vela Jr. SNR would strongly benefit from future observations with CTA which will be able to offer higher sensitivity and better spatial resolution at VHEs. A spatially resolved comparison of the TeV and X-ray data is therefore crucial for the understanding of the physical processes taking place in this SNR.


1

Also known as CXOU J085201.4−461753.

References

All Tables

Table 1.

Parameters of the hydrodynamic model.

Table 2.

Models description.

Table 3.

Parameters of the alternative hydrodynamic model.

Table 4.

Description of the emission model for the alternative hydrodynamic solution.

All Figures

thumbnail Fig. 1.

Time evolution of the forward shock radius (top panel) and shock velocity (bottom panel) up to the current age as obtained from our hydrodynamic model.

Open with DEXTER
In the text
thumbnail Fig. 2.

Magnetic field strength as a function of the normalized radius between the contact discontinuity and the forward shock. Red solid line shows the profile for the magnetic damping scenario and the blue dashed line for the transported magnetic field. Both distributions are shown for the models with the Bohm diffusion coefficient.

Open with DEXTER
In the text
thumbnail Fig. 3.

SED of the Vela Jr. SNR. Lines represent simulated emission from the source produced in different processes as specified in the legend in the magnetic field damping model (left panels) and synchrotron cooling model (right panels). Green triangles indicate the observed radio emission as detected by Parkes (Duncan & Green 2000), an orange bow-tie shows the spectral fit of the X-ray XMM-Newton data (Aharonian et al. 2007), red squares represent the Fermi-LAT data points (Tanaka et al. 2011) and blue circles show the H.E.S.S. data points (H.E.S.S. Collaboration 2018).

Open with DEXTER
In the text
thumbnail Fig. 4.

Zoom into X-ray (top panel) and GeV–TeV (bottom panel) spectra for all four models.

Open with DEXTER
In the text
thumbnail Fig. 5.

Simulated projected X-ray radial flux profiles for different models assuming spherical symmetry compared to the experimental measurements (green squares Kishishita et al. 2013). Normalized surface brightness in the energy range 2−10 keV is plotted as a function of the angular distance from the forward shock of the SNR with negative values corresponding to the interior of the remnant. The radial profile obtained within the asymmetric model discussed in Sect. 4.2 is represented by the black filled circles connected with the dashed line.

Open with DEXTER
In the text
thumbnail Fig. 6.

Simulated projected radial profiles of the X-ray spectral parameters for four different models assuming spherical symmetry. Top panel: cut-off energy and the bottom one shows the spectral index in the 2−10 keV band as a function of the angular distance from the forward shock of the SNR. The distributions are obtained by extracting the spectrum for each region and fitting it with an exponential cut-off power law assuming a photon spectral index of Γ = 1.6(top panel) and with a simple power law (bottom panel). In both panels, green squares indicate the observational data (Kishishita et al. 2013).

Open with DEXTER
In the text
thumbnail Fig. 7.

Simulated radial surface brightness profiles of the gamma-ray emission in the 300 GeV−20 TeV energy range compared to the H.E.S.S. measurements for the northern part of the remnant (Aharonian et al. 2007, black filled circles). Red solid and blue dashed lines show simulated profiles for the spherically symmetric models BOHM_DAMP and BOHM_KOLM, respectively. The radial profile obtained within the asymmetric model discussed in Sect. 4.2 is represented by the green dash-dotted line. This profile is averaged over the azimuth angle for the hemisphere containing the cone of enhanced emission. The thick gray line shows the azimuthally averaged radial profile for the northern part of the SNR extracted from the recently published updated H.E.S.S. excess skymap (H.E.S.S. Collaboration 2018). This skymap is smoothed with a Gaussian function of width 0.08. All simulated profiles are normalized in the same way as H.E.S.S. measurements in Aharonian et al. (2007).

Open with DEXTER
In the text
thumbnail Fig. 8.

Simulated cut-off energy (top panel) and spectral index (bottom panel) of the X-ray spectrum in the filament as a function of the angular distance from the shock compared to the experimental measurements. Simulated distributions are obtained using the asymmetric model described in Sect. 4.2.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.