Introduction

Cosmic rays (CRs) are extremely energetic particles, mainly composed by protons. It is widely accepted that the bulk of cosmic rays (below the knee at approximately 3 × 1015 eV) stems within our Galaxy, playing a significant role in its energy budget1. Several convincing indications point towards supernova remnants (SNRs) as their most likely source2.

Shock waves in SNRs can accelerate particles through the first-order Fermi mechanism (or diffusive shock acceleration, DSA). The capability of SNRs of accelerating electrons is testified by the ubiquitous synchrotron radio emission detected at their shock fronts3, associated with approximately 1–10 GeV electrons. X-ray synchrotron emission from ultrarelativistic (about 10 TeV) electrons has been first detected in the remnant of the supernova observed in 1006 AD4 (hereafter SN 1006) and, afterward, in other young SNRs5. Accelerated protons in SNRs might provide γ − ray emission through interaction with protons in their environment leading to π0 production and subsequent decay. TeV γ − ray emission from a handful of shell-like SNRs has been observed6, while GeV emission has been detected in about 30 remnants7. In SN 1006, γ − ray emission has been detected in both the TeV and GeV bands, with the HESS observatory8 and Fermi telescope9, respectively, but its nature is uncertain, since γ − ray emission can also be produced by inverse Compton from ultrarelativistic electrons (leptonic scenario). Indications for a hadronic origin of the γ − rays have been obtained in a few cases10,11,12,13, confirming that SNR shocks can indeed accelerate hadrons.

The identification of SNRs as the main Galactic factory of CRs is still based on plausibility arguments and many important open issues need to be addressed14. To prove that SNRs are CR factories, it is necessary to show that they indeed supply the power needed to sustain the Galactic CRs (of the order of 1041 erg s−1). Considering the rate of supernovae in our Galaxy (about 2 per century), SNRs should transfer about 10–20% of their characteristic kinetic energy (approximately 1051 erg) into CRs15. The loss of such a large fraction of the ram energy is expected to alter the shock dynamics with respect to the adiabatic case. Nonlinear DSA predicts the formation of a shock precursor (travelling "ahead" of the main shock) that modifies the shock structure. This shock modification should result in an increase of the total shock compression ratio, rt, and a decrease of the post-shock temperature with respect to the Rankine-Hugoniot (adiabatic) values16,17,18,19. Recent self-consistent hybrid (kinetic ions and fluid electrons) simulations show that efficient acceleration of CRs leads also to the formation of a shock postcursor where non-linear magnetic fluctuations and CRs drift away from the shock front, moving downstream20,21. This postcursor is quantitatively more important for shock modification than the precursor; it acts as an additional energy sink, providing an increase of rt, even with a moderate acceleration efficiency: when the CR pressure is about 5–10% of the bulk ram pressure, the total compression ratio ranges approximately between rt = 5–7.

SN 1006 is an ideal target to reveal shock modification. Thanks to its height above the Galactic plane (approximately 600 pc), the remnant evolves in a fairly uniform environment (in terms of density and magnetic field). Moreover, SN 1006 shows regions with prominent particle acceleration, where shock modification might be present, together with regions where there are no indications of particle acceleration and where we do not expect shock modification. In particular, the bilateral radio, X-ray and γ − ray morphology of SN 1006 clearly reveals nonthermal emission in the northeastern and southwestern limbs. X-ray synchrotron emission in nonthermal limbs, seen notably above 2.5 keV, identifies sites of electron acceleration to TeV energies, whereas the lack of nonthermal emission in the southeastern and northwestern limbs marks regions without considerable particle acceleration, where the X-ray emission is mainly thermal (see Fig. 1a). It has been shown that the efficiency of diffusive shock acceleration increases by decreasing the angle θ between the ambient magnetic field and the shock velocity22. There is strong evidence that the bilateral morphology of SN 1006 can be explained in the framework of this quasi-parallel scenario, with the ambient magnetic field, B0, oriented approximately in the southwest-northeast direction23,24,25. If these nonthermal limbs are also sites of efficient hadron acceleration, the signature of shock modification is expected to be stored in the X-ray emission of the shocked interstellar medium17 (ISM). The amount of shock modification, namely an increase in the post-shock density, should raise near the nonthermal limbs (i.e., in quasi-parallel conditions), being smaller in thermal regions, where the shock velocity and B0 are almost perpendicular. Shock modification is expected to reshape the remnant structure, by reducing the distance between the forward shock and the outer border of the expanding ejecta (i. e., the "contact discontinuity"). However, this effect (observed in SN 100626,27) can also be explained as a natural result of ejecta clumpiness28, without the need of invoking shock modification.

Fig. 1: SN 1006 in different energy bands.
figure 1

a Chandra flux images (photons cm−2 s−1) in the 0.5–1 keV (green) and 2.5–7 keV (blue) bands, and column density (in 1020 cm−2) of HI in the [+5.8, +10.7] km s−1 range30 (red). Lines mark angles (θ) relative to the ambient magnetic field, B. The white segment indicates an angular distance of \(5^{\prime}\). b radio continuum map62 (Jy beam−1) at 1.4 GHz (red), with Balmer Hα emission31 (10−17 erg s−1 cm−2, in light blue). c, d 0.5–1 keV and 2.5–7 keV maps, respectively, with the 9 regions selected for spectral analysis superposed in blue. e close-up view of the Hα map. The dashed circle marks the position of the forward shock. Minimum and maximum values for each panel are indicated in the corresponding color bar. All scales are linear, except for the light blue in b, which is in square root scale.

Full size image

A first attempt to observe azimuthal variations in shock compressibility was performed by analyzing the XMM-Newton X-ray observations of the southeastern limb of SN 100629. In this region, the X-ray emission is mainly thermal but contaminated by the contribution of shocked ejecta. Nevertheless, a spatially resolved spectral analysis showed a faint X-ray emitting component associated with the shocked ISM. The density of this component can be inferred from its emission measure (EM, defined as the integral of the square of the thermal particle density over the volume of the X-ray emitting plasma). The post-shock density shows hints of azimuthal modulation, with a minimum around θ = 90∘ (i.e., quasi-perpendicular conditions) and a rapid increase toward the nonthermal limbs, suggestive of a higher shock compressibility. This analysis was restricted to a small portion of the shell (approximately θ = 90∘ ± 35∘), without including regions with quasi-parallel conditions, and limited by the spatial resolution of the XMM-Newton telescope (comparable to the distance between the shock and ejecta). It was thus not possible to isolate regions with ISM emission only, and all the spectra were contaminated by the ejecta contribution, requiring additional assumptions (namely, pressure equilibrium between ejecta and ISM) to estimate the volume occupied by the ISM and then its density.

Here, we identify shock modification in SN 1006 by studying the azimuthal profile of the post-shock density, and overcome the aforementioned limitations by combining deep X-ray observations performed with XMM-Newton (https://www.cosmos.esa.int/web/xmm-newton) and Chandra telescopes (https://chandra.harvard.edu/), to benefit from spectroscopy sensitivity and high-spatial resolution of the post-shock region.

Results

Spatially resolved spectral analysis

Figure 1a shows the soft (0.5–1 keV, mainly thermal emission) and hard (2.5–7 keV, nonthermal emission) X-ray maps of SN 1006, together with HI distribution in the ambient medium30. Panel b shows the Balmer Hα image31 and the radio continuum map. The ambient magnetic-field orientation is superimposed on the maps, with angle θ = 0∘ assumed at the center of the northeastern limb. The position of the forward shock is indicated by its Balmer optical emission (Fig. 1b, e) in thermal limbs, and by the hard X-ray synchrotron filaments in nonthermal limbs (Fig. 1a).

The spatial distribution of the ambient atomic hydrogen around SN 1006 is quite homogeneous30,32 and the only structures detected are localized in the western and northwestern rims (Fig. 1a). None are observed in the regions considered for our analysis, where the ambient medium appears uniform. While we cannot exclude small fluctuations in the ambient density, they are undetected with current radio data. We consider two additional pieces of evidence supporting a uniform upstream medium in this sector of the remnant. The first one is the fairly circular shape of the shock front from the southern to northeastern rim: the shock velocity (and then its radius) depends on the ambient density, and the almost constant shock radius shown in Fig. 1 indicates a uniform environment. The second piece of evidence is the faint and fairly uniform surface brightness of the Hα emission (whose intensity depends on the local density) in this sector of the shell (Fig. 1e). These two conditions support the scenario of a tenuous and homogeneous ambient medium. Therefore, we interpret azimuthal modulations in the post-shock density as ascribed only to variations in rt.

We define nine narrow spatial regions immediately behind the forward shock for X-ray spectral studies, excluding the regions contaminated by the ejecta (see lower panels of Fig. 1). The superior spatial resolution of Chandra allows us to identify the outermost ejecta, whose projected position in the plane of the sky is marked by ripples of thermal emission. Abrupt variations in the X-ray surface brightness of thermal emission show the position of the contact discontinuity. In particular, the X-ray surface brightness of the outermost ejecta is more than 10 times larger than that of the background (Sout, i.e., the surface brightness outside the shell), while the X-ray surface brightness within regions in the southeastern limb (i.e., regions −3, −2,−1, 0, +1) is only 2 − 5 × Sout. The inner border of these regions corresponds to a surface brightness contour level at 6 × Sout, thus marking the sharp separation between ejecta and ISM emission. In regions 2, 3, 4, 5, the contribution of synchrotron emission to the X-ray surface brightness dominates. Therefore, in this part of the remnant we selected very narrow regions, immediately behind the shock front, by carefully excluding visible ejecta clumps. We investigate a large portion of the shell (θ = 0∘–122∘), including regions with quasi-parallel conditions where shock modification is expected to be at its maximum.

Figure 2a shows the X-ray spectrum extracted from region 0, revealing the shocked ambient medium at approximately θ = 90∘, where we do not observe prominent particle acceleration and the shock is adiabatic (rt = 4). The spectrum can be modelled by an isothermal optically thin plasma in non-equilibrium of ionization (parametrized by the ionization parameter τ, defined as the time integral of the density reckoned from the impact with the shock). The post-shock density of the plasma is \(n=0.16{4}_{-0.016}^{+0.014}\) cm−3, in good agreement with previous estimates in this part of the shell29, (as well as the electron temperature, approximately kT = 1.35 keV, and \(\tau=4.{8}_{-4.7}^{+0.9}\times 1{0}^{8}\) s cm−3), but our values are obtained without any ad-hoc assumptions. The interstellar absorption is expected to be uniform in the portion of the shell analyzed, and we fix the absorbing column density to NH = 7 × 1020 cm−2, in agreement with radio observations32. We detect the shocked ISM emission in all regions, with a statistical significance >99%. Figure 2b shows that the ISM contribution is clearly visible at low energies even in the spectrum of region 5 (at approximately θ = 0∘) where the X-ray emission is dominated by synchrotron radiation. We found that the electron temperature in the shocked ISM is consistent with being constant over all the explored regions (though with large uncertainties), and letting it free to vary does not improve the quality of the fits. So we fixed it to the best-fit value obtained in region 0 (kT = 1.35 keV), where we obtain the most precise estimate (see, methods subsection X-ray data analysis). In case of shock modification, we expect the ion temperature to be lower in regions with efficient hadron acceleration. This should also reduce the electron temperature, but this effect is predicted to be quite small33. Moreover, this reduction may be compensated by the fact that strong magnetic turbulence in quasi-parallel regions should favor electron heating (by heat exchange with ions). We expect much larger variations of the shock compression ratio, so we trace the shock modification by focusing on the post-shock density. We estimate the density from the best fit value of the EM in each spectral region (by numerically computing the volume of the emitting plasma, see, methods subsection X-ray data analysis).

Fig. 2: Spatially resolved spectral analysis.
figure 2

Chandra X-ray spectra (black crosses) of regions 0 (a) and 5 (b) of Fig. 1, with the corresponding best fit models (solid lines) and residuals (lower insets in each panel). Error bars are at the 68% confidence level. Data and models are folded through the instrumental response (ACIS-I and ACIS-S in region 0 and region 5, respectively). Thermal (ISM) and nonthermal (synchrotron) contributions are highlighted with dotted lines.

Full size image

Azimuthal profile of the post-shock density

Figure 3a shows the azimuthal modulation of the post-shock density of the ISM. We verified that the estimates of the density are not affected by contamination from the ejecta. Ejecta can be easily identified because of their higher surface brightness with respect to the ISM. If we extract the spectra by slightly changing the shape of the extraction regions, so as that their inner boundaries include ejecta knots, the plasma density increases artificially (and the quality of the fit decreases). For example, by moving inwards the inner boundary of region 0 (region 3), and enhancing its projected area by less than 10% (<40% for region 3), we find that the plasma density changes from \(n=0.16{4}_{-0.016}^{+0.014}\) cm−3 up to n = 0.189 ± 0.011 cm−3 in region 0 (and from \(n=0.2{1}_{-0.04}^{+0.05}\) cm−3 to n = 0.25 ± 0.03 cm−3 in region 3). Conversely, by reducing the size of the extraction regions, the plasma density stays constant, thus indicating that the medium within each region is fairly uniform (e.g., by reducing the projected area of region 0 and region 3 by about 25%, we find \(n=0.15{8}_{-0.021}^{+0.015}\) cm−3, in region 0, and \(n=0.2{2}_{-0.05}^{+0.06}\) cm−3 in region 3). We conclude that Fig. 3 traces the azimuthal density modulation of the shocked ISM.

Fig. 3: Modulation of the shock compressibility.
figure 3

Azimuthal profile of post-shock density (a) and compression ratio (b) derived from Chandra (blue crosses) and XMM-Newton (grey boxes) spectra. Errors are at the 68% confidence level. Angles are measured counterclockwise from the direction of ambient magnetic field. Compression ratios were obtained by assuming a compression ratio of 4 in Chandra region 0 and in XMM-Newton region e (see Table 1, respectively. The solid curve marks the profile expected for parallel efficiency ξp = 12%, normalized magnetic pressure ξB = 5%, and efficiency of CR re-acceleration ξs = 6%, dashed curve is for ξp = 18% and ξs = 0%, dotted curve is for ξp = 12%, ξs = 6% and ξB = 0% (i. e., without including the effects of the postcursor). Source data are provided as a Source Data file.

Full size image

As explained above, it is thus hard to explain this density modulation as a result of inhomogeneities in the upstream medium. These inhomogeneities, if present, would affect the shock velocity (vs, which is proportional to the inverse of the square root of the ambient density), inducing a Δvs of approximately 1200 km s−1 between region 0 and region 5 (for a shock velocity in region 5 of31 5000 km s−1). This would produce a difference in the shock radius of about 1018 cm (corresponding to \(0.5^{\prime}\) at 2.2 kpc34) between region 0 and region 5 in only 250 yr. This is at odds with observations, which show a very circular shape of the shock front (see Fig. 1e), whose radius vary less than \(0.15^{\prime}\) between region 0 and region 5 (see, methods subsection X-ray data analysis for details). We then consider the density modulation as a tracer of azimuthal variations in the shock compression ratio. Assuming rt = 4 in region 0 (i. e., at θ = 90∘, where the acceleration is inefficient), Fig. 3b shows a higher compressibility in quasi-parallel conditions, where the shock compression ratio raises up to approximately rt = 7.

To further constrain the observed increase of the post-shock density towards quasi-parallel conditions, we added the data from the XMM-Newton Large Program of observations of SN 1006 (approximately texp = 750 ks). We updated the previous results obtained for 8 regions in the southeastern limb29 (from around θ = 55∘ to approximately θ = 120∘) to correct for the effects of the telescope point spread function (see, methods subsection X-ray data analysis). In addition, we extended the study to quasi-parallel regions by analyzing the XMM-Newton spectra including region 3 and region 4–5 (together). We adopted the same model as that adopted for Chandra data. The agreement between results obtained with the two telescopes is remarkable (see Fig. 3) and the combination of the reliability provided by the high Chandra spatial resolution (which excludes contamination from ejecta) and the sensitivity of XMM-Newton (which provides more precise estimates) confirms the azimuthal modulation of the post-shock density.

In the XMM-Newton spectra, the higher post-shock density in region 4–5 with respect to region 0 is confirmed even by letting the electron temperature and interstellar absorption free to vary in the fitting process. In particular, we found that the quality of the fit of the spectrum of region 4–5 worsens significantly by imposing the plasma density to be the same as that observed in region 0, even by letting both NH and kT free to vary (χ2 = 182.3 with 181 d.o.f., with \(kT=1.{0}_{-0.3}^{+0.5}\) keV and NH = 6.3 ± 0.6 × 1020 cm−2, to be compared with χ2 = 179.6 with 182 d.o.f., see Table 1). Moreover, we notice that by imposing a low post-shock density in region 4-5, the best-fit value of the ionization parameter (τ = 1.3 ± 0.2 × 109 s cm−3) increases by a factor of about 2 with respect to that reported in Table 1, thus still indicating the need for a higher post-shock density in quasi-parallel regions.

Table 1 Best fit values of emission measure (EM), ionization parameter (τ), cutoff energy of the synchrotron emission (Ecut) and post-shock density (nISM) derived from Chandra and XMM-Newton spectra extracted from the regions shown in Fig. 1, with the corresponding values of χ2 and degrees of freedom (d.o.f.)
Full size table

Discussions

From observations, the azimuthal profile of the post-shock density shown in Fig. 3 can be explained by a higher shock compression ratio in quasi-parallel regions than in quasi-perpendicular regions.

From theory, self-consistent kinetic simulations22 unveiled that the injection of thermal protons into DSA is maximum at quasi-parallel regions, where efficiency can reach 10–15%, and suppressed at quasi-perpendicular shocks. Prominent shock modification (rt ≲ 7) is therefore expected only in quasi-parallel regions20,21, while quasi-perpendicular shocks should show approximately rt = 4. While injection of thermal ions is typically suppressed for22θ ≳ 45∘, re-acceleration of pre-existing Galactic CRs (seeds) can provide a modest acceleration efficiency (approximately 2–6%), thus playing a role in shock modification up to35 θ ≳ 60∘, in agreement with the observed trend in Fig. 3. CR acceleration also leads to the amplification of the pre-shock magnetic field in the quasi-parallel regions, where synchrotron emission is more prominent and magnetic fields are more turbulent, as attested by radio polarization maps25. Also, quasi-perpendicular regions exhibit synchrotron emission in the radio but not in the X-rays, which points to the presence of GeV but not TeV electrons, again consistent with acceleration in the unperturbed Galactic magnetic field2.

The solid curve in Fig. 3b shows the expected trend of rt vs. θ, obtained by assuming a quasi-parallel acceleration efficiency ξp = 12%, with a cut-off at approximately 45∘. Such values are chosen based on self-consistent hybrid simulations36, which attest to thermal ions being spontaneously injected into DSA for quasi-parallel shocks37. Always guided by simulations, which show that more oblique shocks are able to efficiently accelerate pre-existing CR seeds38, we have also a component of reaccelerated CRs with efficiency ξs = 6% and cutoff at 70∘. Since both accelerated and re-accelerated particles can effectively amplify the initial magnetic field, we pose the normalized magnetic pressure ξB = 5%. We also note that magnetic-field amplification and high rt are consistent with the non-detection of X-ray emission from the precursor of SN 100639.

The profile shown by the solid curve in Fig. 3b is in line with the observed data-points. As a comparison, we also include the expected profiles obtained without CR re-acceleration (with ξp = 18%, dashed curve), and without postcursor effects20 (ξp = 12%, ξs = 6%, ξB = 0, dotted curve).

The theoretical curves in Fig. 3, while capturing the zero-th order azimuthal dependence of rt on θ, do not account for the possible refined geometry of the magnetic field in the SN 1006 field. A comparison between radio maps and MHD simulations (not including shock modification)24 suggests that the local magnetic field may be tilted by ϕB ≈ 38∘ ± 4∘ with respect to the line of sight, with a gradient of the field strength of the order of ∇ ∣B∣ = 1.5B over 10 pc, roughly laying in the plane of the sky (parallel to the limbs) and pointing toward the Galactic plane. In general, a finite tilt ϕB would stack regions with different shock inclinations along each line of sight, thereby reducing the contrast between the regions with maximum and minimum rt; nevertheless, since the expected rt is maximum for θ ≲ 40∘ (see Fig. 3), any tilt ϕB ≲ 40∘ would not induce any major modification to the curves in Fig. 3. The presence of a gradient, instead, has been shown24 to reduce the angular distance between the two polar caps, producing a narrow minimum in the rt versus θ profile, remarkably similar to that observed (see, methods subsection Modeling the shock modification). Since the simple geometry assumed in this paper captures well the bilateral morphology of SN 1006 and its azimuthal variations, we defer the study of the corrections induced by more elaborate field geometries to a forthcoming publication.

Finally, we consider the effects that the shock modification should induce on the spectrum of the accelerated particles. In the context of the classical theory of non-linear DSA40,41, a shock compression ratio larger than 4 leads to CR spectra harder than E−2, while radio observations suggest for SN 1006 a radio spectral index of 0.642, corresponding to a CR spectrum ∝ E−2.2. Nevertheless, when the postcursor physics is taken into account in the calculation of the CR spectral index21, for the parameters chosen in the paper for the q-parallel regions (ξtot = ξp + ξs = 18%, ξB = 5%), one obtains rt = 6.34 and a CR spectrum ∝ E−2.19, in remarkable agreement with the observed radio index. In conclusion, our findings show an azimuthal modulation of the post-shock density in SN 1006, which is consistent with a substantial deviation of the shock compression ratio from the value of rt = 4 (expected for strong shocks) in regions of prominent particle acceleration, where electrons are accelerated to TeV energies. The inferred values of compression ratios and CR slopes are compatible with those expected in CR-modified shocks when the effects of the postcursor are also accounted for20,21. Moreover, the azimuthal variation of rt (Fig. 3) attests to the prominence of parallel acceleration and to the important role played by the re-acceleration of pre-existing Galactic CRs for oblique shocks.

Methods

X-ray data analysis

We analyzed Chandra observations 13737, 13738, 13739, 13740, 13741, 13742, 13743, 14423, 14424, 14435 (PI F. Winkler) performed between April and June 2012, with a total exposure time of 669.85 ks, and observation 9107 (PI R. Petre) performed on June 2008 for a total exposure time of 68.87 ks (see Table 2). All observations were reprocessed with CIAO 4.1243 and CALDB 4.9.0 (https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/caldb_intro.html).

Table 2 List of observations analyzed in this work
Full size table

Mosaic images of SN 1006 were obtained by combining the different pointings with the CIAO task merge_obs (https://cxc.cfa.harvard.edu/ciao/ahelp/merge_obs.html). In particular, we produced vignetting-corrected mosaic images of the flux (in counts s−1 cm−2) in the 0.5 − 1 keV band (shown in green in Fig. 1) and in the 2.5–7 keV band (light blue in Fig. 1).

The contact discontinuity in SN 1006 is very close to the forward shock26,27. To measure the ISM post-shock density, we then extract X-ray spectra by selecting narrow regions between the contact discontinuity and the shock front. Regions selected for spatially resolved spectral analysis are shown in Fig. 1. By assuming θ = 0∘ at the center of the northeastern radio limb, the azimuthal range explored is θ = 0∘ − 122∘. In this azimuthal range, the spherical shape of the shock front, combined with the extremely faint and uniform HI emission, clearly point toward a uniform ambient environment. We do not consider regions with negative values of θ because of the lack of spherical shape in the remnant therein, combined with the superposition of several shock fronts (which make it difficult to correctly estimate the volume of the X-ray emitting plasma). We do not consider regions with θ > 122∘ because it is not possible to select regions not contaminated by the ejecta emission, given that several ejecta knots reaching the shock front (and even protruding beyond it) can be observed in the soft X-ray image (Fig. 1c) for approximately θ = 122∘–150∘. Beyond approximately θ = 150∘ the shell loses its spherical shape and interacts with an atomic cloud30,44 (Fig. 1a).

Spectra, together with the corresponding Auxiliary Response File, ARF, and Redistribution Matrix File, RMF, were extracted via the CIAO tool specextract (https://cxc.cfa.harvard.edu/ciao/ahelp/specextract.html). Background spectra were extracted from regions selected out of the remnant, without point-like sources and, when possible, in the same chip as the source regions. We verified that the results of our spectral analysis are unaffected by the selection of other background regions. Spectra were rebinned by adopting the "optimal binning" procedure45. As a cross-check, we also rebinned the spectra so as to get at least 25 counts per spectral bin, obtaining the same results, though with slightly larger error bars. Spectral analysis was performed with XSPEC version12.10.1f46 in the 0.5–5 keV band, by adopting χ2 statistics. Spectra extracted from the same region of the sky in different observations were fitted simultaneously. We found out that all our results do not change significantly by modeling the spectrum of the background, instead of subtracting it, and by using Cash statistics instead of χ2-minimization in the fitting process.

Thermal emission from the shocked ISM was described by an isothermal plasma in non equilibrium of ionization with a single ionization parameter, τ (model NEI in XSPEC, https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node195.html, based on the database AtomDB version 3.09, see http://www.atomdb.org/Webguide/webguide.php). Though we adopt a state-of-the-art spectral model, we acknowledge that there may be limitations in the description of the emission stemming from an under-ionized plasma with a very low ionization parameter, such as the one studied here. However, we expect this effect to introduce almost the same biases (if any) in all regions, and not be responsible for the density profile shown in Fig. 3. We found that the electron temperature in the shocked ISM is consistent with being constant over all the explored regions and fixed it to the best-fit value obtained in region 0 (kT = 1.35 keV), where we get the most precise estimate (error bars approximately 0.4 keV at the 68% confidence level). This value is in remarkable agreement with that measured in a similar region with XMM-Newton29. We included the effects of interstellar absorption by adopting the model TBABS (https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node268.html) within XSPEC. The interstellar absorption is expected to be uniform in the portion of the shell analyzed and we fixed the absorbing column density to NH = 7 × 1020 cm−2, in agreement with radio observations32. We performed the F-test in all regions, finding that the quality of the spectral fittings does not improve significantly by letting kT, or NH free to vary. The ISM emission measure and ionization parameter, τ, were left free to vary in the fitting procedure.

We verified that this model provides an accurate description of spectra extracted from regions in the thermal southeastern limb (namely regions 0, −1, −2, −3, +1) and an additional nonthermal component does not improve significantly the quality of the fits, its normalization being consistent with 0 at less than the 99% confidence level. However, in regions +2, +3, +4, +5 there is a significant synchrotron emission. We then added a synchrotron component when fitting the spectra from these regions and modeled the synchrotron emission by considering the electron spectrum in the loss-dominated case47, since this model is particularly well suited for SN 100648 (our results and conclusions do not change by adopting an exponentially cut-off power-law distribution of electrons (XSPEC/SRCUT model, https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node228.html) to describe synchrotron emission, as done in previous works27,29). Normalization and break energy of the synchrotron emission were left free to vary in the fittings. The normalization of the thermal component is significantly larger than 0 at the 99% confidence level in all regions.

Table 1 shows the best fit results for all the regions, with errors at the 68% confidence level. We derive the average plasma density, \(\overline{n}\), in each spectral region from the corresponding best-fit value of the emission measure of the plasma (\(EM=\int {n}^{2}dV=\overline{{n}^{2}}V\), where n is the plasma density and V is its volume). The volume is calculated with the following method (see Supplementary Software 1 for further details). We project the regions shown in Fig. 1 on a uniform grid with pixel size 0.2″ × 0.2″ (0.2″ correspond to about 6.3 × 1015 cm at a distance of 2.2 kpc34). For each pixel, we calculate the corresponding depth as the length of the chord along the line of sight intercepted by the sphere that maps the shock front, and compute the volume accordingly, We then sum over all the pixels within a given region. The radius of the sphere marking the shock front slightly depends on the region considered, ranging from \({R}_{min}=14.4^{\prime}\) in region +5 to \({R}_{max}=14.55^{\prime}\) in region 0, but we use the same center for all the regions (namely, α = 15h: 02m: 55.74s, \(\delta=-4{1}^{\circ }:56^{\prime} :56.603^{\prime\prime}\)). We verified the precision and reliability of our method by considering more regular regions, like those adopted in previous works29, where the volume can be calculated analytically. We found differences < 0.4% between the numerical and analytical values. The volumes of the emitting plasma in the regions adopted for spectral analysis are listed in Table 1 and were used to derive \(\overline{n}\) from EM. We verified that our results and conclusions do not change by adopting the PSHOCK model within XSPEC to model the thermal emission. The PSHOCK model assumes a linear distribution of the ionization parameter versus emission measure49, ranging from zero (at the shock front) up to a maximum value (τmax which is a free parameter in the fit), instead of a single, "mean", ionization parameter as the NEI model. The best-fit ISM density obtained in region 0 and region 5 with the PSHOCK model is \({n}_{0}^{P}=0.16{3}_{-0.017}^{+0.014}\) cm−3 and \({n}_{5}^{P}=0.3{2}_{-0.08}^{+0.11}\) cm−3, respectively (to be compared with \({n}_{0}=0.16{4}_{-0.016}^{+0.014}\) cm−3 and \({n}_{5}=0.2{9}_{-0.07}^{+0.10}\) cm−3 obtained with the NEI model). As expected, the maximum ionization parameter is approximately a factor of 2 higher than the mean τ obtained with the NEI model (\({\tau }_{0}^{max}=8.{9}_{-1.5}^{+2.1}\times 1{0}^{8}\) s cm−3 and \({\tau }_{5}^{max}=1.{2}_{-0.4}^{+1.1}\times 1{0}^{9}\) s cm−3, to be compared with \({\tau }_{0}=4.{8}_{-0.7}^{+0.9}\times 1{0}^{8}\) s cm−3 and \({\tau }_{5}={7}_{-2}^{+5}\times 1{0}^{8}\) s cm−3).

Table 1 shows the best-fit values of the ionization parameter \(\tau=\int\nolimits_{{t}_{s}}^{{t}_{f}}ndt=\overline{n}\overline{{{\Delta }}t}\) (where \(\overline{n}\) is the time-averaged plasma density, and \(\overline{{{\Delta }}t}\) is the mean time elapsed since the shock impact within the region) in all regions. Figure 4 shows the confidence contours of the ISM density (as derived from the EM) and τ, in region 0 and region 5. Since both EM and τ depend on the plasma density, it is possible to estimate \(\overline{{{\Delta }}t}\). The figure includes isochrones in the (n, τ) space, showing that we obtain very reasonable estimates of the mean time elapsed after the shock impact (\(\overline{{{\Delta }}t}=1-2\times 1{0}^{2}\) yr).

Fig. 4: Ionization parameter and post-shock density.
figure 4

68%, 90%, and 99% confidence contour levels of the ISM density and ionization parameter derived from the Chandra spectra of region 0 (a) and region 5 (b). Red and blue lines mark isochrones after the interaction with the shock front.

Full size image

The radial size of the extraction regions changes from case to case, and so does their inner boundary, which is closer to the shock front in some regions (e.g., region 3, 4, 5) than in others (e.g., region −1 and region −2). Therefore, \(\overline{{{\Delta }}t}\) is not strictly the same for all regions, and the ionization parameter does not depend only on the plasma density (we expect lower \(\overline{{{\Delta }}t}\) in the narrow regions around the northeastern polar cap). However, we find that, especially for regions with similar radial size, higher values of the post-shock density are associated with higher values of τ, as shown in Table 1 and Fig. 5. The azimuthal profile of the ionization parameter shown in Fig. 5 is then consistent with the density profile shown in Fig. 3.

Fig. 5: Modulation of the ionization parameter.
figure 5

Azimuthal profile of the ISM ionization parameter derived from Chandra spectra (blue crosses) (see Table 1). Errors are at the 68% confidence level. Angles are measured counterclockwise from the direction of ambient magnetic field. The red curve marks the profile expected for parallel efficiency (ξp = 12%, with ξB = 5% and ξs = 6%), assuming the same \(\overline{{{\Delta }}T}\) for all regions. Source data are provided as a Source Data file.

Full size image

In the framework of the XMM-Newton Large Program of observations of SN 1006 (PI A. Decourchelle), we analyze the EPIC observation 0555630201 (see Table 2). XMM-Newton EPIC data were processed with the Science Analysis System software, V18.0.0 (see the Users Guide to the XMM-Newton Science Analysis System", Issue 17.0, 2022 https://xmm-tools.cosmos.esa.int/external/xmm_user_support/documentation/sas_usg/USG/). Event files were filtered for soft protons contamination by adopting the ESPFILT task (https://xmm-tools.cosmos.esa.int/external/sas/current/doc/espfilt/espfilt.html), thus obtaining a screened exposure time of 89 ks, 94 ks and 51 ks for MOS1, MOS250 and pn51 data, respectively. We selected events with PATTERN≤12 for the MOS cameras, PATTERN≤4 for the pn camera, and FLAG = 0 for both.

We extracted EPIC spectra from region 3 and the union of region 4 and region 5 (hereafter region 4–5, extraction regions are shown in Fig. 1). Spectra were rebinned by adopting the optimal binning precedure and spectral analysis was performed with XSPEC in the 0.5–5 keV band by adopting the model described above for the analysis of Chandra spectra. MOS and pn spectra were fitted simultaneously. We point out that Chandra and XMM-Newton spectra were fitted independently.

Best fit values are shown in Table 1, with errors at the 68% confidence level. Regions selected for spectral analysis are located at the rim of the shell and part of the ISM X-ray emission is spread outside the outer border of the regions, because of the relatively large point spread function of the telescope mirrors (corresponding to about 6″ full width half maximum). We quantified this effect by assuming that the ISM emission is uniformly distributed in each spectral region and found that approximately 7% of the ISM X-ray emission leaks out of each region. We address this issue by correcting the measured plasma emission measure accordingly. This has a small effect on the density estimate, considering that the density is proportional to the square root of the emission measure. However, we applied this correction to derive the density estimates shown in Table 1 and in Fig. 3 as well as to revise the previous values obtained in the southeastern limb29 and also shown in Table 3 and Fig. 3.

Table 3 Updated values of density of the shocked interstellar medium from previous XMM-Newton data analysis29
Full size table

Modeling the shock modification

Efficient acceleration of CRs has always been associated with an increase in the shock compressibility16,41,52 as due to the softer equation of state of relativistic CRs, whose adiabatic index is 4/3 (rather than 5/3) and the escape of particles from upstream, which effectively makes the shock behave as partially radiative53,54. In this case, though, CR spectra would become significantly harder than E−2 above a few GeV, at odds with γ − ray observations of individual SNRs55,56.

Unprecedentedly-long hybrid simulations of non-relativistic shocks have recently revealed an effect that was not accounted for in the classical DSA theory, namely that the CR-amplified magnetic turbulence may have a sizable speed with respect to the shocked plasma, resulting in a postcursor, i.e., a region behind the shock where both CRs and magnetic fields drift away from the shock faster than the fluid itself20,21. The postcursor-induced shock modification has two main implications: on one hand, it acts as a sink of energy, which leads to an enhanced compression, and on the other hand it advects CRs away from the shock at a faster rate, which leads to steeper spectra57. The relevance of the postcursor is controlled by the post-shock Alfvén velocity20 relative to the downstream fluid velocity, and can therefore be inferred from observations in which shock velocity and downstream density and magnetic field are constrained; simple estimates for both radio SNe and historical SNRs return a remarkable agreement between observations and theory21,58.

It has been shown20 that it is possible to calculate the shock compression ratio given the post-shock pressures in CRs and magnetic fields (ξc and ξB), normalized to the upstream bulk pressure. We then consider the contribution of CRs injected from the thermal pool20,22,59 and re-accelerated seeds35, which both are expected to produce magnetic turbulence via the Bell instability for strong shocks35,60. The dependence on the shock obliquity θ is modelled after hybrid simulations and modulated with

$$\xi (\theta )\equiv \frac{{\xi }_{i}}{2}\left[2+\tanh \left(\frac{{\bar{\theta }}_{i}-\theta }{\Delta \theta }\right)-\tanh \left(\frac{\pi -{\bar{\theta }}_{i}-\theta }{\Delta \theta }\right)\right].$$
(1)

We consider i = [p, s, B], corresponding to the pressure in CR protons injected from the thermal pool, reaccelerated seeds, and B fields, respectively; with such a prescription, each pressure is maximum at ξ(0) = ξ(π) and drops over an interval Δθ = 20∘ centered at \({\bar{\theta }}_{i}=[4{5}^{\circ },7{0}^{\circ },7{0}^{\circ }]\), respectively. The solid line in Fig. 3 shows the prediction for efficiencies ξp = 12%, ξs = 6%, and ξB = 5%; all these values are not the result of a best fitting, but rather motivated by simulations22,35,60 and successfully applied to the study of individual SNRs11,58.

Note that, with the current parametrization, the total efficiency at parallel regions is ξp + ξs ≈ 18%, which is reasonable when acceleration of He nuclei is added on top of protons61.

We also explored a different configuration of the ambient magnetic field, by including the effects of a gradient of the field strength (as suggested by the slantness of the radio limbs24) on the shock obliquity. By adopting the formalism described above, with ξp = 12%, ξs = 6%, and ξB = 5%, we obtain the profile shown in Fig. 6.

Fig. 6: Modulation of the shock compressibility in a non-uniform magnetic field.
figure 6

Same as Fig. 3, with the solid curve marking the profile expected for ξp = 12%, ξB = 5% and ξs = 6%, but including a gradient of the magnetic field strength, lying in the plane of the sky at θ = 90∘. Source data are provided as a Source Data file.

Full size image