{'fraction_non_alphanumeric': 0.03748738778361317, 'fraction_numerical': 0.028018523788580448, 'mean_word_length': 4.08739141879442, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The gas density around SN 1006

The density of the ambient medium where the supernova remnant evolves is a relevant parameter for its hydrodynamical evolution, for the mechanism of particle acceleration, and for the emission at TeV energies. Using XMM-Newton X-ray observations, we present a study of the ambient medium density of the historical supernova remnant SN 1006. We modelled the post-shock thermal emission to constrain the ambient medium density. Our study is focused on the North-West and the South-East rims of the remnant, where the thermal emission dominates. We used a plane-parallel shock plasma model plus another component for the ejecta that are not negligible in the regions of our study. The importance of the synchrotron component is also studied. In order to improve statistics, we combined several observations of the remnant. The density found in the South-East rim is low, roughly 0.05 cm-3, and seems to be representative of the rest of the remnant. However, in the North-West rim (close to the bright optical filament), the density is significantly higher (about 0.15-0.25 cm-3). This confirms a picture of SN 1006 evolving in a tenuous ambient medium, except in the North-West where the remnant has recently encountered a denser region. A density this low is compatible with the non-detection of the remnant by the HESS gamma-ray observatory. The lower density in the South-East implies a higher shock speed of 4900 km/s, higher than that of 2890 km/s measured in the North-West. This new estimate of the velocity could increase the maximum energy that accelerated particles can reach to energies of about 1 PeV.


Introduction
The thousand-year-old Supernova Remnant (SNR) SN 1006 is a type Ia SNR located at a distance of 2.2 kpc (Winkler et al. 2003). The high Mach number collisionless shock of the remnant heats and ionizes the interstellar medium (ISM). Furthermore, the blast wave is thought to be an efficient accelerator of cosmicrays (CR).
Evidence that electrons are accelerated at the shock is found at several wavelengths. The radio synchrotron emission (Reynolds & Gilmore 1986) arises from electrons accelerated to energies of at least 1 GeV. The X-ray synchrotron emission in SN 1006 (Koyama et al. 1995) comes from electrons of ∼ 20 TeV. Evidence concerning proton acceleration is more scarce. The best hope is to detect the gamma-ray emission of pion decay resulting from the interaction of accelerated protons (and heavier ions) with the thermal gas. However, the TeV emission can also result from Compton scattering of electrons at ≃ 20 TeV on the CMB background. Currently there are only a few shell SNRs detected at those energies (Cas A, Vela Junior, RX J1713.7-3946; Aharonian et al. 2001Aharonian et al. , 2005Aharonian et al. , 2006 and SN 1006 is not yet one of them. To interpret the TeV emission, it is important to estimate the gamma-ray flux from the π 0 mechanism. One of the most important parameters for predicting this flux is the gas density in the γ-ray emitting region. Ksenofontov et al. (2005) have shown that, in order to be consistent with the non-detection of SN 1006 by the HESS experiment, the ambient medium density must be ≤ 0.1 cm −3 . The ambient medium density is also a key parameter for the hydrodynamical evolution of the remnant. As SN 1006 results from a type Ia SN, the surrounding medium is expected to be unaltered by the progenitor. Moreover, the remnant is relatively isolated, as it is about 500 pc above the galactic plane. Thus it is a good case study for comparing observations with theoretical models of CR acceleration.
In the literature, the gas density around SN 1006 ranges from 0.05 ≤ n H ≤ 1.0 cm −3 . The different estimates were based on radio HI study (Dubner et al. 2002), UV spectra (Laming et al. 1996;Korreck et al. 2004;Raymond et al. 2007), and X-ray emission (Hamilton et al. 1986;Winkler & Long 1997;Long et al. 2003). In SN 1006, while the limb-brightened synchrotron is dominant in the North-East and in the South-West, the thermal X-ray emission from the shocked gas dominates in the South-East (SE) and in the North-West (NW). It is possible to deduce the characteristics of this shocked gas, and in particular its density, with X-ray observations. In the NW, such a study has already been carried out by Long et al. (2003) with Chandra.
Using XMM-Newton, we present a study of the density in the South-East and in the North-West by modeling the thermal emission.

Data processing
The first observations of the remnant's South-East region by the XMM-Newton satellite (Rothenflug et al. 2004) were not very deep. The effective exposure time in this region was 6 ks (ObsID 0111090601) in comparison with 33 ks for the North-West (ObsID 0077340101). A second observation was carried out to provide better coverage of this region, on 2005 August 21 and 22 with the medium filter for a duration of 35 ks. Unfortunately, the proton flare contamination was high and the remaining exposure time is only 11 ks. Note that for this new observation, the CCD number 6 (no longer working) and the CCD number 4 (very noisy) of the camera MOS1 were not used. All observations used for our study are listed in Table 1 with their remaining exposure time after screening. We used only MOS data due to its better spectral resolution. All the spectra used are the average of MOS1 and MOS2 detectors. The data were processed using the Science Analysis System (SAS version 6.5). All fitting was carried out using unbinned spectra with the C-statistics implemented in Xspec (v12.2.1). Binning was used for graphical purposes only and fixed at 3σ for all spectra. We have fitted and plotted all the data from 0.3 keV to 10 keV.
To clean proton flare contamination from the event files, we used the MOS-FILTER process available in the Extended Source Analysis Software (XMM-ESAS). The following method was used for background subtraction. To take into account the astrophysical background, we extract a spectrum from a region outside of the SNR (regions A and B on Fig. 1); for statistical reasons, the area selected is as large as possible. The instrumental background spectrum is derived from a 960 ks blank sky observations compilation (Carter & Read 2007) in the same detector area as the source and renormalized in the 10-12 keV band (where astrophysical photons are negligible) over the full field of view.

Spectral analysis
To estimate the density of the ambient medium surrounding SN 1006, we model the thermal X-ray emission at the shock. This gives us the post-shock density which is related to the gas density by the compression ratio, which is 4 if particle acceleration is negligible.
Our study is focused on the SE and NW regions where thermal emission dominates. Figure 1 shows the extraction regions used for the spectra. In the NW, we used two regions: one is the region noted NW-1 in the study of the remnant by Long et al. (2003) with Chandra (to relate our study to what has been done previously) and the second one is the faint emission region and is noted NWf. Because the emission is faint in this region, we used a box as large as possible (1.2'x7.3'). The NW-1 region, in comparison, is much smaller (17"x3'). In the SE (region ESE and SSE), the width and length were fixed to 1' and 6.5' respectively, as a compromise between good statistics and being close to the shock.

Method
Our method is based on fitting the spectral data with thermal models and infering the post-shock density from the emission measure (EM) of the shocked ISM. Using Xspec, we get the parameter norm=10 −14 Ω/4π n e n H dl which is proportional to the EM. In this expression, Ω is the solid angle of emitting gas contained in our extraction region and dl the integration element along the line of sight. We assume that the shocked ISM fills all the volume inside the SNR. This implies a filling factor of 1.
For a solar abundance plasma, the electron and hydrogen density are related by n e = 1.21×n H . Therefore, for a given value of norm (obtained from the fit), we can solve the previous equation if we know the density profile to compute the integral. We assumed a Sedov radial profile for that purpose. An important advantage of this method is that the hydrogen density is proportional to the square root of the ratio norm/l. Thus the density is weakly sensitive to variations of EM when using different thermal models, and to errors due to line of sight estimate. The error due to a filling factor of 1 is also minimized for the same reasons.

What kind of model can we use ?
The pertinence of our estimate of the density will depend on the quality of the spectral model for the shocked ISM. As the ionization of the shocked ISM is out of equilibrium, we need to use non-equilibrium ionization (NEI) models. Moreover, the variation inside the SNR of the ionization timescale, defined as τ = n e dt, is important. Starting at zero at the shock, the ionization timescale increases toward the center of the remnant down to ∼ 0.9 R shock for a Sedov model and then decreases. Therefore we integrate cells of gas with different τ along the line of sight.
To model the emission, we used a simple model of planeparallel shock (PSHOCK) with constant temperature, which allows the use of a distribution of τ instead of a single ionization timescale like in a standard NEI model. A plane-parallel shock model is a good local approximation for radially small regions close to the shock, where τ is still increasing inwards. The projection of the radial profiles onto the sky is not modelled by PSHOCK.
To check the consistency of our model, we also determine the density by comparing the maximum ionization timescale τ resulting from the fit with the ionization timescale computed using Sedov solutions (refered to later as the τ method). Sedov solutions have the advantage of taking into account the temporal evolution of the density and the shock velocity. At a given relative angular distance θ/θ shock , the only parameter required is the age of the remnant, which is well known. The method using the norm parameter, which will be refered to later as the EM method, is much more stable and will be used for the results. The τ method is useful for checking the consistency between the models we use. In the following subsections, we use the region ESE as a test region to try different spectral models.

Shocked ISM model
As a first step, we assumed that synchrotron emission is negligible and that all of the emission comes from thermal emission of the shocked ISM. We therefore used a PSHOCK model with solar abundances (Lodders 2003) and we fixed the interstellar absorption to N H = 7 × 10 20 cm −2 (Dubner et al. 2002). In the NEI model (v1.1 in Xspec) the oxygen He-K δ and ǫ lines are not included but they seem to be necessary to reproduce the large "shoulder" of the He-like oxygen line complex. To take them into account, we added one gaussian at 0.723 keV to represent the oxygen series He-K δ and higher. The spectral fit of region ESE is shown in Fig. 2 (top panel) and the resulting best-fit parameters are given in Table 2 (see "one-component"). The fit is poor and we clearly see that some emission lines like magnesium (at 1.35 keV) and silicon (at 1.86 keV) cannot be reproduced by this model. Such emission lines require oversolar abundances and must come from the ejecta.

Shocked ISM with ejecta model
To reproduce the spectrum correctly, we added another component corresponding to the ejecta emission (the same PSHOCK model as above but with non-solar abundances).
To reduce the number of free parameters in this twocomponent model, we fixed most of the abundances of the ejecta. As there is strong evidence that SN 1006 is a Type Ia SN, we can obtain an estimate of the synthesized mass of elements and thus of the abundances in the ejecta through numerical simulations (Iwamoto et al. 1999). Since there is neither hydrogen nor helium in type Ia explosions, the abundances are given relative to oxygen with respect to cosmic abundances from Lodders (2003). We fixed the abundances of elements not constrained by emission lines in the spectrum to the following values: O=1.0 (by definition); C=1.0; N=0.1; Ar=2.0; Ca=3.0. Abundances from elements with prominent emission lines (Ne, Mg, Si) were let free.
SN Ia are expected to produce an iron mass of 0.6 M ⊙ , which is equivalent to a relative abundance of 20. This value is not compatible with our data because the iron emission lines are very faint in the observed spectrum. The value preferred by the fit is compatible with zero. In practice, we fixed it to Fe=0.1. We can only see the emission from the shocked ejecta, and the spatial distribution of iron is complex (Wang & Chevalier 2001). In our extraction region there appears to be little iron.
The fit is much improved in comparison with the onecomponent model with ∆C-stat=495 (see Table 2 for C-stat and best-fit parameters). We note that the ejecta component is required for a good fit at Mg and Si emission lines, but that it also contributes to the emission at low energy and in particular to the oxygen lines. In the first model, all the emission was attributed to the shocked ISM, whereas the emission of the ejecta is important, as seen in Fig.2 (middle panel). Using the EM method for the shocked ISM component, we obtained a density n H =0.053 (0.048-0.057) cm −3 . The value derived from the τ method is consistent with the EM method : n H =0.059 (0.05-0.07) cm −3 . This model reproduces the spectrum well and the two methods give a consistent estimate of the density.

Impact of a potential synchrotron component
SN 1006 is well known to be a non-thermal X-rays emitter. Most of the emission in the bright SW and NE limbs is due to synchrotron emission. In the regions of our study, the thermal emission dominates at low energy, but can we completely neglect the non-thermal emission?
In order to answer to that question, we built a threecomponent model by adding a SRCUT model. We fixed the radio spectral index to 0.6 (Dyer et al. 2001). The radio flux was extracted from the same regions as the X-rays using the radio map from Rothenflug et al. (2004). This radio image is a combination of VLA data at the average frequency of 1517.5 MHz and single-dish Parkes radiotelescope observations.
The spectral fit is improved with a ∆C-stat=43. The resulting cut-off frequency of 4.5 × 10 16 Hz (about 200 eV) is compatible with the frequency found by Rothenflug et al. (2004) in their Fig. 7 (at a position angle 130 degrees). As we can see in Fig.2 (bottom panel), the synchrotron component dominates the spectrum at high energy, at low energy, the shocked ejecta component dominates.
In the two-component model, the shocked ISM component was constrained by the spectrum at high energy, while the ejecta were constrained mainly by the Mg and Si emission lines. Here, as the synchrotron reproduces the high energy, the shocked ISM is no longer well constrained. In order to correctly reproduce the data, we need a continuum and a thermal model for the emission lines. The continuum can be handled either by a high temperature thermal model or by the synchrotron model. When we include both in the same fit, it is hard to disentangle them. The emission lines can be more easily reproduced by the shocked ejecta as several abundances are left unconstrained. Thus the shocked ISM component is almost not needed. The best-fit for a model with only synchrotron and shocked ejecta has a C-stat that is higher by only ∆C-stat=2. . Bottom: same model as above but with a synchrotron component (solid line). In the bottom panel, we note that the shape of the shocked ISM component is similar to the synchrotron component. The shocked ISM is almost not needed, and the power law shape is due to a very low ionization timescale τ ≃ 10 8 cm −3 s (no emission lines). For the middle and the bottom panel, the global fit (not shown) is a very good fit to the data.
Including the synchrotron in the model could lower the estimate of the density. However, using hydrodynamical arguments we can set a lower limit on the density of the ISM, as explained later in Section 5.4. Therefore we think that the statistically preferred model is not the best physical model.
In conclusion, there are good physical reasons to include the synchrotron, but the data are not able to estimate its contribu-tion. We decided to use the simpler two-component model. The resulting densities should be viewed as upper limits.

Using different sizes of extraction region
The width of our extraction regions ranges from 17" for region NW-1 up to 1' for the regions ESE and SSE. We checked the impact of varying the size of the extraction region on the density estimate. In order to do so, we reduced the width of the region ESE to 17", keeping it immediately behind the shock. Using the two-component model described in Sect. 3.4, we derive densities from the small ESE region with the EM and τ methods of 0.065 (0.055-0.075) cm −3 and 0.06 (0.048-0.09) cm −3 respectively. The density estimate is fully consistent with that determined with a larger ESE region (see Sect 3.4). Therefore the density estimate is not sensitive to the size of the extraction region.

General results
Following the spectral analysis of region ESE (Sect. 3), we now use the same model for all the regions: shocked ISM plus shocked ejecta (as in Sect. 3.4). The results for the four regions (ESE, SSE, NWf, NW-1) studied here are listed in Table 3.
For all the regions, the ejecta component has greatly improved the quality of the spectral fit in comparison with a onecomponent model. The emission lines are then well reproduced with our model as shown in Fig. 3.
The density obtained with the EM method in the SE rim (ESE and SSE region) and in the NW faint region (NWf) is globally low (n H =0.049-0.053 cm −3 ). Moreover the plasma parameters of the shocked ISM obtained from NWf are similar to those found in the SE. This is consistent with the smooth appearance of the Hα image (Winkler et al. 2003). We also note that the temperature of the shocked ISM is comparable in those three regions: kT ISM =1.75-2.1 keV. This ISM temperature is much higher than the temperature of the shocked ejecta (0.5 keV). The shocked ISM component dominates the spectrum at high energy.
As we can see on the spectra in Fig. 3, the ejecta component is important below 2 keV. Its contribution to the Mg, Si and Ne emission lines is obvious, and we can also see that part of the emission from the oxygen lines is produced by shocked ejecta. In the NW regions, an overabundance of Ne is required. This has also been noted by Long et al. (2003) in their spectral analysis of the bright filament.

The bright filament
The bright filament seen in X-rays in the NW is coincident with a bright Hα filament. A possible explanation is a "recent" interaction of the expanding SNR with a denser ISM. The sudden lower shock velocity flattens the shape of the blast wave in this region as it is seen in optical and X-rays. The thermal emission at the shock increases because the shock propagates a denser region and thus creates a brighter filament than in the SE. The fact that the shock is entering in a denser region in the NW is particularly visible in Hα. The radio study by Dubner et al. (2002) also detected a higher density in the NW. In Hα and in radio, we also see the faint structure of the region NWf. A possible interpretation of this faint border is that the remnant is expanding in front of (or behind) the bright filament in a region of lower density similar to that of the SE. Table 3. Best-fit parameters obtained with the shocked ISM plus shocked ejecta model (two-component) for the different regions. The gas density derived with the EM and the τ methods is also given (see Sect. 3.1). Due to the faint emission in region NWf, the temperature and the abundances of the shocked ejecta component (parameters in square brackets) were fixed to the values of region NW-1. The lower and upper bounds are given in parentheses at a 90 % confidence level. * : The density estimate for region NW-1 assumes constant density and velocity rather than Sedov profiles (see Sect. 4.2).  To derive the density from the plasma parameters, we used a simple Sedov model to compute τ and EM and then compare them with the best-fit parameters (see Sect. 3.1). However, in the NW-1 region it seems that part of the shock wave has recently encountered a denser region. For this reason, the geometry of the shell is too complex and probably not well represented by a Sedov density profile. The density derived with this method is probably not reliable in this particular region. Nevertheless, with a simpler model (assuming constant density and constant velocity) we can derive the density from the ionization timescale ). Indeed, we know the shock speed in this region (2890 km/s, Winkler et al. 2003) and know that region NW-1 has a radial extent of 17" (corresponding to an elapsed time of 240 years since the gas was shocked, at a distance of 2.2 kpc). Thus for τ=4.5 (3.5-5.5) ×10 9 cm −3 s, the density is n H ≃ 0.15 (0.11 − 0.18) cm −3 . We note that kT ISM is lower in region NW-1 (1.45 keV) than in the SE and NWf (1.7-2.1 keV). This is consistent with the fact that the density in region NW-1 is higher.
In comparison, the study of the bright filament region by Long et al. (2003) with Chandra found a density of 0.25 cm −3 . To see if there are differences between the data observed by XMM-Newton and Chandra, we tried the same simple solarabundance PSHOCK model (described in Sect. 3.2 of Long et al. 2003) on our data. The absorption column density they use is higher than our value: 8.9 vs 7.0 ×10 20 cm −2 . We modelled the data with both values. The best-fit parameters are given in Table  4. The model with low N H is statistically preferred with a ∆C-stat of 35. Table 4. Comparison of the best-fit parameters obtained with XMM-Newton and Chandra data in the NW-1 region. kT e (keV) τ (10 9 cm −3 s) XMM Chandra XMM Chandra N H =8.9 ×10 20 cm −2 0.6 0.6 6.8 6.9 N H =7.0 ×10 20 cm −2 0.7 ... 6.5 ...
Both data give similar parameters and the fact that we are not using the same value for N H does not have an important impact on the plasma parameters. Thus, the difference in the estimate of the density is not due to the data; it is mostly due to the difference of models. When taking into account the ejecta using a two-component model, the density derived is slightly lower than when using a simple solar abundance PSHOCK model.
The difference of density between the NW-1 and other regions was already pointed out by Korreck et al. (2004) and Long et al. (2003). Korreck et al. (2004) estimated the density in the NE to be ≤ 0.06 cm −3 , which is compatible with what we found in the SE.
We conclude that a low value of ∼ 0.05 cm −3 is probably representative of the ambient medium around SN 1006, except in the bright NW filament where the density is higher (0.15-0.25 cm −3 , depending on the model).

Comparison with other observations
By modelling the post-shock thermal emission of the remnant in the SE and in the NW, we found a density of ∼ 0.05 cm −3 . Many other studies (local and global) have been carried out to measure this density and provide estimates in the range of 0.05 cm −3 to 1.0 cm −3 .
Most of the local ones were focused on the bright Balmer line filament in the NW. Using HST observations of this filament, Raymond et al. (2007) derived a density in the range of 0.25-0.4 cm −3 by modelling the Hα brightness behind the shock. By using the high spatial resolution X-ray observatories, some studies also estimated the density in the NW region. Winkler & Long (1997), with ROSAT, used a model for the X-ray surface brightness profile behind the shock; they found a density of 1.0 cm −3 . However, using the Chandra X-ray Observatory, Long et al. (2003) studied the same region and found a much lower value. They derived a density of 0.25 cm −3 from the ionization timescale in a region right behind the shock. The spectral analysis of the same region using XMM-Newton observations gave an identical value when using the same model (see Sect. 4.2). A model taking into account the ejecta component gave a slightly lower value (0.15 cm −3 ). We note that the density derived from the X-rays is lower than that derived from the optical. By using X-ray observations, we average the density on the size of the extraction region, whereas in optical we measure the density right at the shock. If the shock wave is entering into a cloud, it is plausible that the difference of density between optical and X-rays is related to the gradient of density when entering the cloud.
There is a general agreement that the remnant is encountering a denser region in the NW, so the density there is not directly comparable with what we found in the SE.
The global studies average the density all around the remnant. A surrounding density of ∼ 0.3 cm −3 has been inferred by Dubner et al. (2002) with a survey of the HI emission. However, this result is the average density along a line of sight of several hundreds of pc, which is very large in comparison with the typical size of the remnant of 20 pc. Hamilton et al. (1986) derived a density of 0.05 cm −3 by modelling the global X-ray spectrum under the assumption that most of the emission is thermal. Spatially resolved spectral analysis showed that this assumption is only valid locally in the NW and in the SE. The global spectrum is in fact dominated by synchrotron emission. Therefore, even though their result agrees with ours, we view this as coincidental. Numerical simulations have also been carried out to study the evolution of SNRs in constant density ISM. In order to be consistent with the observed angular size and expansion rate of SN 1006, Dwarkadas & Chevalier (1998) derived from their simulation a density between 0.05 cm −3 and 0.1 cm −3 . Those results are compatible with what we found. We note that our density is lower than values of most other studies, but most of them were focused on the NW, where the ambient medium seems denser. The value found in the SE region is more likely to be representative of the overall ambient medium density, including in the bright synchrotron limbs.

Limitations of our model
Our method is based on estimating the density of the shocked ambient medium. To infer the density of the unshocked ambient medium, we have to assume a value for the compression ratio. In the efficient acceleration sites (NE and SW rims) the compression ratio can increase to 6 as calculated in Berezhko et al. (2002). However, in the SE and NW rims, the acceleration process is not as efficient, as we see very little X-ray or radio synchrotron emission. Thus a compression ratio of 4 seems reasonable.
To derive the density from the EM, we assumed that the shocked ISM filled all the volume i.e. a filling factor ε=1. As seen in Sect. 3.1, we need ejecta to reproduce the spectra and it must occupy some fraction of the volume. In the ejectadominated phase (Chevalier 1982), the fraction of the volume occupied by the shocked ISM in the global shocked region is about 0.8. SN 1006 is in a more advanced phase (towards Sedov) where the ejecta are more diluted in the shocked ISM. Moreover, closer to the shock, as in our regions, ε is supposed to be closer to unity. So in our model, we can reasonably say that ε = 0.8 is a lower limit to the filling factor. In that case, the density estimate for region ESE would increase from 0.053 cm −3 to 0.059 cm −3 .
To derive the density from the norm parameter, we assumed a Sedov density profile to estimate the value of the integral n 2 e dl. As SN 1006 is switching from the ejecta-dominated stage to the Sedov phase, the density profile that we assumed is steeper than it is in reality. We therefore overestimate the density corresponding to a given value of the EM. In the extreme case of a flat density profile, we obtain for the ESE region a density of 0.042 cm −3 in comparison with 0.053 cm −3 for the Sedov profile.
The distance of 2.2 kpc used throughout this paper was derived by Winkler et al. (2003) from the measurement of the shock speed of the bright filament of 2890±100 km/s and the apparent expansion of 0".28 yr −1 . However, Heng & McCray (2007), using the same data as in Winkler et al. (2003) but with a new model that incorporates additional physics, derived a lower shock speed of 2509±111 km/s. The corresponding distance to the remnant with this shock speed is 1.9 kpc. For our study this will change the length of the line of sight in the EM method, but it will have a low impact on the density as n H ∝ √ norm/dl. Concerning the contribution of the synchrotron emission, our three-component model was not able to disentangle the thermal emission of the shocked ISM from the synchrotron emission. We cannot exclude that part of the emission at high energy is due to synchrotron emission, so our density estimate is an upper limit in that respect. A lower limit to the density is discussed later in Sect. 5.4.

Velocity of the blast wave
With the density of 0.05 cm −3 inferred above, we can reestimate several parameters in order to account for the observed radius of 9.6 pc (at a distance of 2.2 kpc). The age is t=1000 years, and as SN 1006 seems to be a type Ia supernova, the ejecta mass is M ej = 1.4 M ⊙ . The impact of the distance uncertainty was low in the density estimate. However, parameters like energy or shock speed are really dependent on this distance. The following results are given for both d=2.2 kpc and 1.9 kpc. As the remnant is just switching from the ejecta-dominated stage to the Sedov stage, we decided to use the equation from Truelove & McKee (1999) for an homogeneous ISM with s=0 and a type Ia ejecta profile with n=7.
The kinetic SN energy consistent with the previous parameters is E SN =2.2/1.3×10 51 ergs. For example, for a density higher than 0.2 cm −3 , the energy required would be higher than E SN =6.0/3.3 ×10 51 ergs. Type Ia SN are expected to release an amount of energy of ∼ 10 51 ergs. Energies above 2×10 51 ergs seem hard to reach. Thus a uniform gas density ≥ 0.2 cm −3 surrounding all the remnant is excluded for energetic reasons.
The global velocity of the blast wave can also be reestimated. As the density is lower in the SE, the shock speed there must be higher than the shock speed in the NW which is 2900 km/s. Indeed, the shock speed derived from the Truelove & McKee equations with an ambient medium density of 0.05 cm −3 at t=1000 years is 4900/4600 km/s. If the density found in the SE region is representative of the rest of the remnant, this velocity is what we can expect for the whole remnant except in the NW. We note that there is also a difference between the global shock speed measured in radio and the velocity in the NW deduced from Hα observations. Indeed, even though the error bars are big, the mean proper motion measured by Moffett et al. (1993) in radio is 0".44±0".13 yr −1 , corresponding to 4500 km/s for a distance of 2.2 kpc. This is much higher than the velocity of 2890 km/s obtained in the NW, and it is compatible with our estimate.

Lower limit on the density
When a shock wave encounters a cloud, a transmitted shock is driven into the cloud. The pressure behind the transmitted shock is greater, at most, by a factor β than the pressure behind the blast wave not in interaction with the cloud. For a given density contrast n cloud /n ISM , we can derive this β factor using the equations from Sgro (1975). The pressure is proportional to n H .V 2 shock . What we are interested in here are the density and velocity right at the shock. We know those values in the NW (transmitted shock) from the optical studies : n H = 0.25 − 0.4 cm −3 and V 2 shock =2890 km/s (Raymond et al. 2007;Winkler et al. 2003, respectively). For any ambient medium density, the velocity (in a region where the shock is not encountering any cloud) is known using the Truelove & McKee equations. Therefore we can derive a "model" ratio of pressures for any ambient medium density. By comparing the theoretical ratio β with the "model" ratio, we can set a lower limit to the ambient medium density. The change of slope in the "model" ratio corresponds to the transition between the ejecta-dominated phase and the Sedov phase. For n cloud =0.25 cm −3 , the "model" ratio becomes higher than the maximum theoretical value of β(0.25/n ISM ) for n ISM lower than 0.03 cm −3 (see Fig. 4). For n cloud =0.4 cm −3 , the lower limit on the ambient medium density is 0.07 cm −3 . We conclude that the density of the ambient medium where the SNR evolves cannot be lower than 0.03 cm −3 . For our preferred value of the density of 0.05 cm −3 , we have a "model" ratio of 1.7 which is compatible with the theoretical value of β(0.25/0.05)=2.1.

Implication for CR acceleration
The efficiency of the acceleration of the particles is strongly dependent on the shock speed. In SN 1006, this velocity has been measured only in the NW region, by Ghavamian et al. (2002) and Heng & McCray (2007). However, the velocity there does not seem to be representative of the rest of the remnant, and in particular, of the efficient acceleration sites in the NE and SW rims. If we use for those regions the velocity estimated in the SE in Sect. 5.3, this has an impact on the highest energy the accelerator can reach. In fact, in a Bohm regime, at an age of 1000 yr we have E max = 32.7 B d 100 µG ( V shock 1000 km/s ) 2 TeV (adapted from Parizot et al. 2006) where B d is the downstream magnetic field. Thus with the higher velocity (4900 km/s) and B d =150 µG as in Ksenofontov et al. (2005), the maximum energy that particles could reach in the bright limbs would be E max ≃ 1 PeV. The TeV emission of those accelerated particles depends on the density of the ambient medium. The total γ-ray flux (IC+π 0decay) from the NE half of the remnant for different values of n H has been given in Ksenofontov et al. (2005). For an upstream magnetic field of 30 µG, a density of 0.05 cm −3 and a distance of d=2.2 kpc, this flux is a factor of 3 below the current HESS upper limit, and only a factor ∼ 2 for d=1.9 kpc. This would, however, need to be revisited in the light of our revised estimate of the shock velocity.

Electron temperature equilibration
A question that remains open in the physics of strong shocks is the mechanism by which the electrons are heated just behind the shock. In the Rankine-Hugoniot relations for high Mach number shocks, the particles are heated in relation to their mass. The proton temperature is higher than the electron temperature by a factor 1836. However, it has also been proposed by Cargill & Papadopoulos (1988) and Ghavamian et al. (2007) that collective plasma interaction can give rise to prompt (limited) electron heating at the shock. Such mechanisms produce a higher electron temperature at the shock than in the Rankine-Hugoniot case. Our model requires an electron temperature kT e ∼ 2 keV to reproduce the spectrum in the region ESE. The Rankine-Hugoniot relations (without particle acceleration) imply a proton temperature just at the shock of ∼ 40 keV for a shock speed of 4900 km/s, if the degree of temperature equilibration is low. Using the Coulomb equilibration mechanism with an initial electron-to-proton temperature ratio of 1/1836 (the electron-to-proton mass ratio), the predicted electronic temperature is lower than what is observed. Indeed, with an ionization timescale of 5×10 9 cm −3 s, kT e rises to 0.8 keV only. An initial ratio of at least T e /T i =0.04 (corresponding to kT e =1.7 keV just at the shock) is needed for the temperature to reach 2 keV. According to Ghavamian et al. (2007), the electronic temperature immediately behind the shock is constant (∼ 0.3 keV) for any shock speed above 400 km/s. In comparison, the electronic temperature just at the shock required in our model (kT e =1.7 keV) is much higher.
However, we note that this high temperature depends very much on the amount of synchrotron emission in the SE, as discussed in Sect. 3.5. To check if the three-component (synchrotron, shocked ejecta and shocked ISM) model is compatible with a low electronic temperature just at the shock, we fixed several parameters of the shocked ISM component for the lowest possible value of the ambient medium density i.e. 0.03 cm −3 . The norm parameter and the upper limit of the ionization timescale can be fixed from Sedov models and from the geometry of the emitting region. They are fixed respectively to 0.8×10 −5 cm −5 and 3 × 10 9 cm −3 s. From the upper limit of τ, we can derive an upper limit to the temperature but not the average (for the region) electron temperature that is needed for the PSHOCK model. In order to avoid this difficulty, we use a plane-parallel shock model with separate ion and electron temperatures (NPSHOCK) to reproduce the shocked ISM. In this model the two temperatures can be set precisely. The mean postshock temperature is derived from the shock speed (22 keV at 4900 km/s), and the electron temperature immediately behind the shock front is set to 0.3 keV. The global fit is good and the C-stat is higher only by ∆C-stat=5 in comparison to the three-component model (where the parameters for the shocked ISM component were left free). Thus the data are compatible with a low electron temperature just behind the shock if a model including synchrotron emission is assumed.

Summary
The study of the post-shock thermal emission in the SE and in the NW of SN 1006 with XMM-Newton observations leads to the following conclusions : 1. Most previous studies were focused on the NW region where the remnant seems to be encountering a denser ISM. By modelling the post-shock thermal X-ray emission, we were able to estimate the density in the SE rim. The low value of ∼ 0.05 cm −3 we found seems representative of the rest of the remnant. 2. The shock speed in the SE region, and maybe in the rest of the remnant, is ∼ 4900 km/s. The velocity of 2890 km/s was deduced from UV observations in the NW where the density is higher, causing the expansion to slow down. 3. The highest energy particles can reach scales as V 2 shock . The new estimate of the velocity could increase E max up to energies of about 1 PeV. 4. If we assume that the ambient medium density found in the SE region extends to the bright limbs where the acceleration is efficient, then according to Ksenofontov et al. (2005) the TeV gamma-ray flux would be a factor of ∼ 3 below the current HESS upper limit.