### 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.12^{43} and CALDB 4.9.0 (https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/caldb_intro.html).

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 shock^{26,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 cloud^{30,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” procedure^{45}. 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.1f^{46} 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 (*k**T* = 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-Newton*^{29}. 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 *N*_{H} = 7 × 10^{20} cm^{−2}, in agreement with radio observations^{32}. We performed the F-test in all regions, finding that the quality of the spectral fittings does not improve significantly by letting *k**T*, or *N*_{H} 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 case^{47}, since this model is particularly well suited for SN 1006^{48} (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 works^{27,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 × 10^{15} cm at a distance of 2.2 kpc^{34}). 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, *α* = 15^{h}: 02^{m}: 55.74^{s}, (delta=-4{1}^{circ }:56^{prime} :56.603^{primeprime})). We verified the precision and reliability of our method by considering more regular regions, like those adopted in previous works^{29}, 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 *E**M*. 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 measure^{49}, 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=intnolimits_{{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 *E**M*) and *τ*, in region 0 and region 5. Since both *E**M* 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-2times 1{0}^{2}) yr).

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.

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, MOS2^{50} and pn^{51} 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 limb^{29} and also shown in Table 3 and Fig. 3.

### Modeling the shock modification

Efficient acceleration of CRs has always been associated with an increase in the shock compressibility^{16,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 radiative^{53,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 SNRs^{55,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 itself^{20,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 spectra^{57}. The relevance of the postcursor is controlled by the post-shock Alfvén velocity^{20} 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 theory^{21,58}.

It has been shown^{20} 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 pool^{20,22,59} and re-accelerated seeds^{35}, which both are expected to produce magnetic turbulence via the Bell instability for strong shocks^{35,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 simulations^{22,35,60} and successfully applied to the study of individual SNRs^{11,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 protons^{61}.

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 limbs^{24}) 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.