H.E.S.S. and Suzaku observations of the Vela X pulsar wind nebula

Pulsar Wind Nebulae (PWNe) represent the most prominent population of Galactic very-high-energy gamma-ray sources and are thought to be an efficient source of leptonic cosmic rays. Vela X is a nearby middle-aged PWN, which shows bright X-ray and TeV gamma-ray emission toward an elongated structure called the cocoon. Since TeV emission is likely inverse-Compton emission of electrons while X-ray emission is synchrotron radiation of the same electrons, we aim to derive the properties of the relativistic particles and of magnetic fields with minimal modelling. We use data from the Suzaku XIS to derive the spectra from three compact regions in Vela X covering distances from 0.3 pc to 4 pc from the pulsar along the cocoon. We obtain gamma-ray spectra of the same regions from H.E.S.S. observations and fit a radiative model to the multi-wavelength spectra. The TeV electron spectra and magnetic field strengths are consistent within the uncertainties for the three regions, with energy densities of the order $10^{-12}\rm\,erg\,cm^{-3}$. The data indicate the presence of a cutoff in the electron spectrum at energies of $\sim$100 TeV and a magnetic field strength of $\sim$$6\,\rm\mu G$. Constraints on the presence of turbulent magnetic fields are weak. The pressure of TeV electrons and magnetic fields in the cocoon is dynamically negligible, requiring the presence of another dominant pressure component. Sub-TeV electrons cannot account completely for the missing pressure, that may be provided either by relativistic ions or from mixing of the ejecta with the pulsar wind. The electron spectra are consistent with expectations from transport scenarios dominated either by advection via the reverse shock or by diffusion. Constraints on turbulent magnetic fields and the shape of the electron cutoff can be improved by spectral measurements in the energy range $\gtrsim 10\rm\,keV$. (abridged)


Introduction
Pulsars eject relativistic winds which are thought to be loaded primarily with electrons and positrons.The wind is highly supersonic leading to the formation of a termination shock at the distance where the wind ram pressure becomes comparable to the external pressure.Beyond the shock lies a bubble of magnetised relativistic plasma that originated in the pulsar magnetosphere.The formation of these so-called Pulsar Wind Nebulae (PWNe) is accompanied by efficient particle acceleration.Thus, PWNe are bright non-thermal emitters with spectra extending from radio to gamma rays, and represent the dominant class of identified Galactic sources observed at the highest-energy end of the electromagnetic spectrum (H.E.S.S. Collaboration et al. 2018a,b).However, the exact sites and mechanisms of particle acceleration up to PeV energies in PWNe, and how much PWNe contribute to the electron/positron component observed in cosmic rays still remain to be established (e.g., Amato 2014).
Most of the recent progress in the understanding of PWNe, including acceleration and propagation of high-energy particles, has come from theoretical studies of the properties of the magnetohydrodynamic (MHD) flow in the nebula (Kennel & Coroniti 1984b;Bogovalov & Khangoulian 2002;Lyubarsky 2002;Del Zanna et al. 2006;Volpi et al. 2008;Amato 2014;Porth et al. 2014), and observations in the X-ray and TeV gamma-ray domains (for a review see, e.g., Kargaltsev et al. 2015).Under the conditions inferred in PWNe, synchrotron radiation of relativistic electrons should dominate the emission at X-ray energies, and the TeV gamma rays are generated through inverse-Compton (IC) scattering.For typical magnetic fields in PWNe, X-ray and TeV gamma-ray emission is generated by particles with similar energies.The synchrotron emissivity is determined by the electron density and the strength of the local magnetic field.MHD simulations show that the magnetic field strength can vary strongly throughout the nebula.Thus, X-ray data alone do not enable us to determine the particle density or to perform detailed studies of the particle evolution in PWNe.
On the other hand, the IC emissivity is determined by the electron density and the energy density of target photon fields.The latter are, for the production of TeV emission, predominantly the cosmic microwave background (CMB), and farinfrared (FIR) Galactic dust emission, which are unrelated to the nebular processes, contrary to magnetic fields that determine the synchrotron emission.This means that gamma rays are a direct tracer of the particle densities in a nearly model-independent way (except for uncertainties related to the FIR photon field).Once the particle density is estimated, the X-ray emissivity constrains the strength of the magnetic field.Therefore, in principle, the combination of X-rays and gamma rays can provide invaluable information to study particle acceleration and transport in PWNe, and validate numerical MHD simulations.
However, combining information from the X-ray and TeV wavebands is often complicated, not only because of the unavoidable projection effects along the line of sight, but also due to the limited angular resolution of gamma-ray measurements.TeV signals are typically registered from large structures with varying magnetic field strength and particles accumulated over a long fraction of the pulsar's lifetime (e.g., de Jager & Djannati-Ataï 2009).
The latter limitation can be overcome observationally for a sufficiently extended and bright PWN, such as Vela X. Originally discovered as a 3 • × 2 • radio source in the Vela constellation (Rishbeth 1958), owing to its level of polarisation and its flat spectrum Vela X was later interpreted by Weiler & Panagia (1980) as the PWN formed by PSR B0833−45, the Vela pulsar (spin-down power of Ė 7 × 10 36 erg s −1 and characteristic age of τ = 1.1 × 10 4 yr, Manchester et al. 2005).This association places the PWN at a distance of 287 +19 −17 pc from the Earth (Dodson et al. 2003, from parallax measurement of the neutron star).ROSAT revealed an X-ray shell with a diameter of 8 • associated with the supernova remnant (SNR) connected to the Vela pulsar and enclosing Vela X (Aschenbach et al. 1995).ROSAT also revealed a 1 • elongated structure within Vela X that seemingly emanates from the pulsar and was dubbed the "cocoon" (Markwardt & Ögelman 1995) 1 .H.E.S.S., an array of gamma-ray imaging atmospheric Cherenkov telescopes located at an altitude of 1800 m above sea level in the Khomas highland of Namibia, unveiled the TeV emission associated with the cocoon (Aharonian et al. 2006a), making Vela X one of the archetypal TeV PWNe.
To explain multiwavelength observations of Vela X, de Jager et al. (2008) proposed the existence of two distinct electron populations, one corresponding to the extended radio nebula and another to the X-ray/TeV cocoon, both interacting with magnetic fields of similar strength of ∼ 5 µG.Hinton et al. (2011) proposed that the extended radio nebula is filled with a relic electron population, devoid of high-energy particles (> 10 GeV) owing to energy-dependent escape, while the cocoon is filled with electrons accelerated more recently.This scenario is in agreement with hydrodynamical simulations that suggest that the cocoon was formed in an evolved stage of the system, when the reverse shock from the SNR crushed the PWN (Blondin et al. 2001;Slane et al. 2018).These simple two-zone models, however, cannot reproduce the details of the multiwavelength morphology revealed by the latest observations, such as the larger-scale TeV emission overlapping the extended radio nebula (Abramowski et al. 2012), or hints of energy-dependent morphology at GeV energies (Grondin et al. 2013).
In this work we take advantage of archival data from the Suzaku X-ray Imaging Spectrometer (XIS) and of an extended H.E.S.S. dataset to improve the constraints on the properties of the highest-energy particles and magnetic fields in the Vela X cocoon.The Suzaku XIS (Koyama et al. 2007) has a large field of view with a diameter of 18 and provides spectral coverage in the range from 0.2 keV to 12 keV with low background contamination, which enables us to probe emission from the highest-energy electrons in the extended Vela X cocoon (for previous studies, see, e.g., Mattana et al. 2011).The same electrons are also responsible for the multi-TeV emission observed with H.E.S.S., which indicates a cutoff in the underlying particle spectrum (Aharonian et al. 2006a;Abramowski et al. 2012).
This energy range is very interesting because the shape of the cutoff encodes information on the particle acceleration and transport mechanism (see Zirakashvili & Aharonian 2007;Romoli et al. 2017, and references therein), and also because the highest X-ray energies may reveal a spectral hardening from nonuniform-strength magnetic fields (Kelner et al. 2013).Owing to the proximity of Vela X, and subsequently its large apparent size, and its extremely high flux at TeV energies, we can now extract spatially-resolved spectra extending to the highest energies for compact regions in X-rays and gamma rays, covering different distances from the pulsar wind termination shock, thus overcoming some limitations that affect multiwavelength studies of other PWNe.
This paper is structured as follows: in Section 2 we describe the observations and the definition of the analysis regions.Then we describe the X-ray and gamma-ray data analysis in Section 3 and Section 4, respectively.In Section 5 the multiwavelength spectral energy distributions (SEDs) are used to constrain the properties of particles and magnetic fields in the cocoon.Finally, we discuss the results in Section 6 and summarise our conclusions in Section 7.

Observations and Analysis Region Definition
We use three archival Suzaku observations of the Vela X cocoon conducted in 2006.From North (closer to the pulsar) to South, they have observation ID 501109010 (hereafter Pointing 0), 501107010 (hereafter Pointing 1), 501110010 (hereafter Pointing 2), for an exposure of 60 ks, 61 ks, and 18 ks, respectively.
We define three regions for spectral analysis corresponding to the Suzaku pointings as illustrated in Figure 1.Later we use these very same regions to derive a SED from the X-ray data and the gamma-ray data.For Pointings 1 and 2 we use circular spectral extraction regions with a radius of 7.5 centred in the middle of the Suzaku XIS field of view (R.A. = 128 • .7666,Dec. = −45 • .4581,and R.A. = 128 • .6368,Dec. = −45 • .8007,respectively).For Pointing 0 we define a spectral extraction region such that we avoid the region immediately adjacent to the pul-sar.The pulsar itself emits X-rays up to a few keV due to thermal emission from its surface, and to higher energies from magnetospheric particle acceleration and radiation (Pavlov et al. 2001).Furthermore, within 1.5 (0.13 pc) from the star, the pulsar wind creates a complex jet-torus structure very bright in X-rays (Kargaltsev et al. 2003;Manzali et al. 2007), not resolved in gamma rays.We exclude from the analysis a circular region centred on the pulsar with a radius of 3.6 (0.We apply the following selection criteria to the dataset: observations where at least three of the 12 m telescopes were operational to improve high-energy performance, since we are interested in gamma rays at energies > 1 TeV, produced by the same electrons responsible for the emission measured by Suzaku in the keV range; observations free from occasional hardware failures (i.e., malfunctioning camera pixels) to closely match the nominal description of the instruments in the Monte Carlo simulations used to evaluate the array performance; observations with good atmospheric transparency, determined by a cut on the "Cherenkov transparency coefficient" (Hahn et al. 2014) so that the observing conditions are also well-matched to the description of the atmosphere in the Monte Carlo simulations.
The selection yields 100 h of livetime spent observing within 3 • .5 from the centre of Vela X, at R.A. = 128 • .75,Dec. = −45 • .6 (all coordinates are given in J2000 equinox in this paper).

X-ray Analysis
We use the HEASoft package 2 version 6.19 for data reduction and analysis.The data are processed through the Suzaku standard pipeline in which standard event screening criteria are applied to data obtained with all the XIS Charge-Coupled Devices (CCDs, XIS 0-3).We then extract X-ray spectra from the analysis regions in Pointings 0, 1, and 2 (see Figure 1).The instrument response is calculated using the xisrmfgen and xissimarfgen tools.The instrument response is generated assuming uniform emission that extends to a circular region with a radius of 20 , which is sufficiently larger than the field of view and a standard setting of xissimarfgen for data analysis of diffuse emission.We evaluate instrumental background based on night Earth's observations by using xisnxbgen.We use only photons with energies above 2.25 keV to avoid thermal emissions from the Vela SNR, and set an upper energy limit at 10 keV for front-illuminated CCDs (XIS 0, 2, 3) and to 7 keV for the back-illuminated CCD (XIS 1) to minimise effects from particleinduced instrumental background.Each X-ray spectrum is well described by a single power law modified by a fixed Galactic interstellar absorption N H = 2.59×10 20 cm −2 (Manzali et al. 2007) plus underlying cosmic Xray background (CXB), parametrised as in Miyaji et al. (1998), i.e., a power law with a photon index of 1.42 and a normalisation of 10.0 s −1 cm −2 sr −1 keV −1 at 1 keV.We note that, owing to the proximity of the system and the high-energy threshold, the effect of interstellar absorption on the Vela X parameters is negligible.A potential contamination from the Galactic ridge X-ray emission is negligible in this region as Vela X is sufficiently far from the Galactic centre.The spectra and best-fit parameters obtained from the fits are shown in Figure 2 and Table 1, in which all power-law fluxes are calculated over an energy range of 2-10 keV.
The fit results do not show significant differences between Pointings 1 and 2 in the flux level or spectral slope (∼ 2.2-2.3),consistent with XMM-Newton results in Slane et al. (2018), who find that the X-ray spectrum softens at angular distances exceeding 60 from the pulsar, i.e., beyond the bright TeV cocoon.Conversely, in Pointing 0 we find a significantly harder spectrum (photon index 1.92 ± 0.014) and the emission is a factor > 2 brighter than in the other two regions (the flux within the extraction region is similar, but the solid angle subtended is 41% of that in regions 1 and 2).This is consistent with previous studies of X-ray emission in the region with XMM-Newton and BeppoSAX (Mangano et al. 2005).
Based on the spectral analysis, we derived X-ray SED points in 10 energy bins between 2.25 keV and 9.5 keV.We include 2 https://heasarc.gsfc.nasa.gov/docs/software/heasoft/.more, we include additional uncertainties to the SED points of Pointing 0 due to the leakage into the spectral extraction region from emission close to the pulsar.We estimate this uncertainty contribution as the spill-over due to the XIS PSF of the flux from a point source at the position of the pulsar that accounts for the total flux measured in the circular region with radius of 3.6 centred on the neutron star.It amounts to a fraction of the flux that varies from 35% to 50% going from low to high energies and it constitutes the dominant source of uncertainty for Pointing 0.  (a) The flux within the extraction region is reported for the energy range 2 keV to 10 keV.

Gamma-ray Analysis
Since we are mainly interested in the energy range > 1 TeV, to minimise systematics that may be induced by combining different instrumental configurations we use only data from the four 12 m telescopes for the whole time period (the contribution from the 28-m telescope to the effective area in the energy range of interest is unimportant).The gamma-ray energy and direction reconstruction, and separation from the background of cosmic-ray hadrons are based on a multivariate technique applied to parameters of the Cherenkov images of the atmospheric showers (Ohm et al. 2009).Given the large extension of Vela X and low surface brightness, we use tight selection criteria to make the residual contamination from hadrons misclassified as gamma rays as low as possible, i.e., we apply the "hard cuts" from Ohm et al. (2009).Additionally, to ensure uniform energy thresholds among different observing periods, we include in the subsequent analysis only candidate gamma-ray events with reconstructed energy > 0.6 TeV, thus reducing systematics in spectral reconstruction at the lower end of the energy range.
All the results in the paper are cross-checked using a second independent software/analysis chain for calibration, event reconstruction, and selection based on an air-shower model template approach (de Naurois & Rolland 2009).Simulated image templates are fitted to the measured image, in order to obtain the physical properties of the shower progenitor gamma ray.The goodness of fit is used as a discriminant variable between gamma-ray and background events.As for the main analysis pipeline, we use "hard cuts" as defined in Aharonian et al. (2006b).Results are given for the main analysis, and the alternative analysis is used to derive systematic uncertainties when noted.
Figure 1 (right) shows the gamma-ray detection significance for the whole region of the Vela X cocoon.The map is derived using the formalism by Li & Ma (1983) after evaluating the residual hadronic background using the ring method (e.g., Berge et al. 2007).In the latter, the residual background is estimated from the data within the same observation, excluding regions of the field of view where significant gamma-ray emission is potentially present.Specifically, we exclude from background estimation regions with gamma-ray emission detected at significance > 5σ in the H.E.S.S. Galactic plane survey (H.E.S.S. Collaboration et al. 2018b, 3.2.2),or where the brightness temperature at 44 GHz measured with Planck (Planck Collaboration et al. 2014) is > 1.5 mK, which indicates the presence of relativistic electron populations that may radiate in gamma rays 3 .Figure 1 shows that there is very significant gamma-ray emission at energies > 0.6 TeV in all three analysis regions.
We evaluate the residual background for spectral estimation applying the reflected-region method (e.g., Berge et al. 2007).The requirement to have at least two reflected regions for background estimation outside the exclusion region (the same described above for the derivation of the significance map) reduces the livetime to 70 h, 75 h, and 80 h for Pointings 0, 1, and 2, respectively.Taking into account the background count spectra thus estimated, we fit analytic functions to the count spectra in the spectral extraction regions using a maximum likelihood method based on Poisson statistics for both counts and background counts.In Pointing 0 counting statistics above a few tens of TeV are low, with no events recorded at reconstructed energies > 60 TeV, except for one single event at a reconstructed energy of ∼ 90 TeV in the spectral extraction region (no events with reconstructed energies > 60 TeV are found in the background regions).We have verified that all the results presented in the paper are not affected, within statistical uncertainties, by including the ∼ 90 TeV event in the spectral analysis or not.We report in the following the results obtained by selecting only events with reconstructed energy < 60 TeV in Pointing 0.
We consider two analytical functions to represent the source photon spectrum.The first is a simple power law (PL), for which the number of photons, N, as a function of energy, E, varies as where A is the differential flux that normalises the spectrum at E 0 , and Γ is the spectral index of the distribution.Alternatively we consider an exponentially-cutoff power law (ECPL) with the characteristic energy of the cutoff, E cut , as additional parameter.We note that we use Eqs.( 1) and ( 2) just to describe the shape of the gamma-ray spectrum.To infer more detailed information about the physical processes behind the gamma-ray emission in Section 5 we determine the energy distribution of high-energy electrons in the system by fitting their non-thermal radiation to Suzaku and H.E.S.S. data.Let L PL and L ECPL be the maximum likelihood values for the PL and ECPL models, respectively, given the source region and background regions counts.The test statistic TS = 2 log(L ECPL /L PL ) is used to assess whether the additional degree of freedom associated to the cutoff significantly improves the fit.Although the formal criteria to apply the likelihood ratio test (e.g., Protassov et al. 2002) are not applicable 4 , a larger TS Table 2. Best-fit parameters for the spectral models of the gamma-ray spectrum of the Vela X cocoon in the three regions shown in Figure 1.Errors correspond to 1σ uncertainties.The functional forms used to model the spectra are given in Equations 1 and 2. TS is the test statistic for the improvement of the fit when the ECPL model is used, see text for definition.

Pointing Model
Flux  value can still be taken to represent a significant improvement in the fit when using the ECPL model.For Pointings 1 and 2 we obtain TS equal to 13.4 and 25.7, respectively, therefore we select the ECPL as best-fit model.For Pointing 0 TS is 6.4 formally corresponding to a significance of the cutoff at 2.5σ statistical level.The presence of a cutoff in Pointing 0 is not significant.However, we note that the value of the cutoff energy, 16 ± 7 TeV, is consistent with those obtained for the other two Pointings.
The resulting spectral parameters and gamma-ray fluxes for the best-fit functions are reported in Table 2. Systematic uncertainties on the fit parameters are evaluated combining the differences between results from the two different calibration, recon-struction, and event selection schemes used in this work with the other sources of systematic uncertainties as evaluated in Aharonian et al. (2006b), namely 20% on flux, 0.1 on the spectral index, and 30% on the cutoff energy.These include, for the flux, uncertainties due to the Monte Carlo models of the atmosphere and particle interactions used to evaluate the instrument performance, effect of missing camera pixels, uncertainty in the livetime estimate, and, for all parameters, the effect of the choice of event selection criteria and background estimation method, and systematic fluctuations observed either within single runs or on a run-by-run basis.The results from the two different calibration, reconstruction, and event selection schemes are consistent with Article number, page 7 of 19 each other, with the exception of the cutoff energy in Pointing 2, which differs by ∼ 3σ.This deviation is taken into account in the systematic uncertainties.The spectral parameters are consistent with previous analyses of H.E.S.S. data covering larger regions of the Vela X cocoon (Aharonian et al. 2006a;Abramowski et al. 2012).Additional information on the spectral fitting is provided in Appendix A.
The results from the spectral fitting are used to derive a binned SED for each pointing.The spectra are rebinned such that for each point the significance of gamma-ray signal detection is > 2σ, and we require that each bin has at least two excess5 counts.The binned SEDs are shown along with the full-range spectral functions in Figure 3.For Pointings 1 and 2 we consider only the ECPL hypothesis that is strongly preferred based on the data, while for Pointing 0 we show both the results based on the PL and ECPL hypotheses.Systematic uncertainties on the SED points are shown as the sum in quadrature of the systematic error on flux from the sources discussed above, amounting to 20%, with the average difference between the points obtained from the two analyses.The average is performed over energy independently for every pointing, and it amounts to 35%, 24%, and 29% for Pointing 0, 1, and 2, respectively.Figure 3 shows that within statistical uncertainties the binned SEDs agree with the best-fit function determined from the whole energy range.For Pointing 0 we note that the SEDs derived based on the PL and ECPL hypotheses are consistent within the level of the statistical uncertainties.Since the cutoff model is not statistically favoured, for this pointing in the following section 5 we adopt the SED derived from the PL fit.We have verified that all the results are consistent with those obtained from the ECPL SED within the quoted uncertainties.

Radiative Modelling of the Multiwavelength Spectral Energy Distribution
To study the obtained X-and gamma-ray spectra in a consistent way, we compute the synchrotron and IC emission from a population of high energy electrons.The energy distribution of particles is assumed to be a power-law with sub/super-exponential cutoff: where N e is the electron number, E e is the electron energy, A is the normalisation factor for the distribution, α is the spectral index of a power-law distribution with reference energy E 0 , E cut is the cutoff energy, and β is the cutoff index.Leptons produce X-ray emission through synchrotron radiation in a magnetic field, which is assumed to have a random orientation but uniform strength B. The gamma-ray emission is generated through IC scattering on CMB and FIR photon fields.For the latter we use a recent model by Popescu et al. (2017) at the position of Vela X.The relevance of FIR radiation, which was often overlooked in past studies, and the possible uncertainties imposed by the model for the FIR field as seed for the IC scattering in Vela X are discussed in Appendix B (for a general discussion of the contribution of different photon fields to IC scattering in extended gamma-ray sources see also Aharonian et al. 1997).
The computation of the SED models and subsequent fit to the multiwavelength SEDs are performed using the naima Python package (Zabalza 2015).Specifically, synchrotron emission is computed based on the formalism in Aharonian et al. (2010), and IC emission on the formalism in Khangulyan et al. (2014); Aharonian & Atoyan (1981).
naima enables us to derive the best-fit values and posterior probability distributions of the model parameters given the SED points from the χ 2 , calculated assuming that the SED point uncertainties are Gaussian and uncorrelated.Owing to the presence of systematic uncertainties, that, to this end, we combine in quadrature with statistical ones, and to the instruments' energy dispersion, this is only an approximate assumption, that should be overcome in future works through multi-mission fullforward analyses (e.g., Vianello et al. 2015).The model parameters are scanned using the Markov Chain Monte Carlo (MCMC) method, as implemented in the emcee software package (Foreman-Mackey et al. 2013).For all model parameters we assume a flat prior probability distribution, within physical constraints on the parameters values (e.g., particle densities are positive).We scan the cutoff energy E cut in logarithmic space, so that the fit parameter is actually log 10 (E cut /1 TeV).
The fit results are shown in Figure 4 and Table 3. Figure 4 shows the best-fit model SED and model uncertainties compared to the measured SEDs.Table 3 gives median and upper and lower uncertainties on the parameter values.In addition, appendix C contains the probability density distributions of the model parameters.
Figure 4 shows that leptonic models can reproduce the Xray and gamma-ray SEDs for plausible model parameter values (see Table 3).The inferred properties of the electron spectra and strength of the magnetic field are consistent within the uncertainties over distances from 0.3 pc up to 4 pc from the pulsar wind termination shock.As expected for IC radiation in the Thomson regime, the values of α in Table 3 and those of Γ in Table 2 are consistent with the equation Γ = (1 + α)/2.The cutoff exponent β is only weakly constrained by the data.Furthermore, the constraints on the cutoff for Pointing 0 are weaker than for the other regions studied, consistent with the harder X-ray spectrum (Section 3) and the lack of a significant cutoff detection in gamma rays (Section 4).For Pointing 0 the weaker constraints also stem from the larger systematic uncertainties on the X-ray SED, increasing with increasing energy due to the spill-over from the region immediately around the pulsar.
Although we used the X-and gamma-ray spectra extracted from the compact regions, the emission is not necessary produced in completely homogeneous zones.The spectrum extraction regions span distances of ∼ 1 pc, which significantly exceeds the gyro radius of TeV particles.Thus, it cannot be excluded that some MHD inhomogeneity affects the magnetic field strength across the production region.It is important to note that the magnetic field strength can also vary along the line of sight, therefore future improvement of the X-ray and gamma-ray instruments' sensitivity will not remove this uncertainty completely.
If the magnetic field varies in the production area, the synchrotron spectra can be deformed significantly (Kelner et al. 2013).The key parameter that determines the radiation regime in an inhomogeneous magnetic field is the ratio of the magnetic field fluctuation length, λ B , to the characteristic photon formation length, λ ph = m e c 2 /eB (with m e electron mass, e electron charge, c speed of light, and B strength of the magnetic field, see Rybicki & Lightman 1979).Landau damping imposes a rather fundamental constraint on the minimum length of inho- where n e and are electron number density and mean energy per electron, respectively.Since /(m e c 2 ) ∼ Γ PW 1 (Γ PW is the bulk Lorentz factor of the pulsar wind) and B 2 /(4π n e ) ∼ 1 (for energy equipartition between magnetic fields and particles expected in PWNe, Kennel & Coroniti 1984b), one obtains that λ B /λ ph 1.Thus, in PWNe the impact of magnetic field inhomogeneities on the X-ray spectra can be modelled as a superposition of synchrotron spectra produced in magnetic fields of different strength (Kelner et al. 2013).
Therefore, we gauge magnetic turbulence in Vela X against the multiwavelength SEDs by assuming that the probability density function (PDF) corresponding to a magnetic field strength B is: The first term of the summation represents a magnetic field with a fixed strength, i.e., B 0 , with δ being the Dirac distribution.The second term in the summation represents a magnetic field with power-law distribution of index ζ between minimum and maximum strengths B min and B max , respectively, with H being the Heaviside step distribution.The parameter a represents the mixing between the constant and power-law magnetic field components.For illustration we take6 ζ = 3/2, and B max = 100×B 0 , i.e., sufficiently exceeding the mean value to modify the synchrotron spectrum.
The parameters C and B min in Eq. ( 5) are set respectively by requiring that the integral of the PDF is 1, and that the root mean square expectation value for the power-law magnetic field (and, therefore, the total magnetic field as well) is B 0 (i.e., B RMS = B 0 ).The only free parameters are then the mixing a and the root mean square magnetic field expectation value, B 0 , that can be fit to the multiwavelength SEDs.
Table 4 reports median, upper and lower uncertainties on the parameter values.Results for the electron spectra parameters and magnetic field expectation value are generally consistent with those derived for a fixed magnetic field (Table 3).The mixing parameter a is not constrained by the observations in Pointing 0, while for the other two regions small contributions from the power-law magnetic field distribution are favoured: the 99th percentiles of the posterior probability distribution of a lie below the values 0.57 and 0.56 (for Pointings 1 and 2, respectively).The Bayesian Information Criterion (BIC) values in Tables 3 and 4 are smaller for the fixed magnetic field case.Thus, this study does not allow us to claim a detection of the emission generated in a turbulent magnetic field.To a large extent this is not surprising, since the magnetic field turbulence predominantly affects the high energy part, i.e., a region above SED peak which is not well covered by current X-ray measurements.This results also in large uncertainties on the cutoff index β, which is the electron spectrum parameter most coupled with the effect of turbulent magnetic fields.This prevents us from studying in Fig. 4. Multiwavelength spectral-energy distributions (SEDs) from the three regions studied in the Vela X cocoon.Blue points are derived from X-ray measurements with the Suzaku XIS (Section 3), and green points from the gamma-ray measurements with H.E.S.S. (Section 4).Error bars combine statistical and systematic uncertainties as described in text.Top, middle, and bottom panels correspond to the extraction regions 0, 1, and 2, respectively.In the upper subplot of each panel, radiative models are overlaid to the SED points.In the lower subplot, the residuals for the best fit parameters are shown.The models are based on a single population of leptons producing synchrotron emission in a uniform-strength magnetic field and inverse-Compton emission on the cosmic microwave background and diffuse infrared radiation.The black line shows the best-fit (minimum χ 2 ) model, while the shaded bands represent the 95% confidence level bands from the Markov Chain Monte Carlo scan.
Article number, page 9 of 19 Table 3. Best-fit parameters for the radiative model of the X-ray and gamma-ray spectrum of the Vela X cocoon in the three regions shown in Figure 1.

17.3
Notes.Median values from the MCMC scan, with lower and upper uncertainties based on the 16th and 84th percentiles of the posterior distribution. (a) Total electron energy for particle energies > 1 TeV.Note that the solid angle subtended by region 0 is 41% of that in regions 1 and 2.
(b) Parameters of the electron spectrum as defined in Equation 3. (c) Strength of the magnetic field. (d) Bayesian information criterion, i.e., k where k is the number of parameters estimated from the model and n is the number of data points.
Table 4. Best-fit parameters for the radiative model of the X-ray and gamma-ray spectrum of the Vela X cocoon for the case of turbulent magnetic field in the three regions shown in Figure 1. Pointing 0.09 ± 0.04 2.6 +1.9 −1.3 7.6 ± 0.9 0.39 +0.11   −0.14 19.9 Notes.Median values from the MCMC scan, with lower and upper uncertainties based on the 16th and 84th percentiles of the posterior distribution.
(a) Total electron energy for particle energies > 1 TeV.Note that the solid angle subtended by region 0 is 41% of that in regions 1 and 2.
(b) Parameters of the electron spectrum as defined in Equation 3.
(c) Parameters of the magnetic field distribution as defined in Equation 5. (d) Bayesian information criterion, i.e., k • ln n + χ 2 , where k is the number of parameters estimated from the model and n is the number of data points.
detail the shape of the particle spectrum cutoff, that encodes information about the particle acceleration and transport mechanisms (see Zirakashvili & Aharonian 2007;Romoli et al. 2017, and references therein).

Pressure balance in the cocoon
As shown in Table 3, the electron distribution and magnetic field strength are consistent within the uncertainties between the three regions considered.For a magnetic field strength of 6 µG the energy density is On the other hand, for a total energy of electrons with energies > 1 TeV of W e ∼ 10 44 erg (see Table 3), the particle energy density is where we have taken r = 0.63 pc (which corresponds to the spectrum extraction radius of 7.5 at the source distance of 290 pc), and is the rather uncertain size of the emitting region along the line of sight.Another rather uncertain parameter, F , is the filling factor that determines the fraction of the volume filled by relativistic electrons.For large filling factors, F 1, the cocoon is exclusively occupied by relativistic electrons, and smaller values imply a significant mixing of relativistic plasma and SN ejecta.The projected size of the cocoon is about 2 pc × 7 pc (which corresponds to angular size of 0.4 • × 1.4 • at the source distance of 290 pc), thus it is feasible that 1 pc.Unless is large, 1 pc, the TeV particle pressure coincides within a factor of a few with the magnetic field pressure in the cocoon for F 1.
Observations of the jet-torus structure in the inner nebula by Chandra point to a size of the wind termination shock R s 1.3 × 10 17 cm (i.e., 30 for a distance of 290 pc, where 30 is the projected radius of the X-ray arc in Fig. 2 of Helfand et al. 2001).For the pulsar spin-down luminosity of Ė 7 × 10 36 erg s −1 the Rankine-Hugoniot relations for a weakly magnetised ultrarelativistic pulsar wind (Kennel & Coroniti 1984a) where R snr and t snr are SNR radius and age, and n ism is the nucleon number density in the interstellar medium (ISM) surrounding the shell, which is inferred to be 1 − 2 cm −3 (Dubner et al. 1998).
Equations 8 and 9 rely on completely different observational constraints and therefore provide a robust order-of-magnitude estimate of the pressure of 10 −9 erg cm −3 .The pressure estimates by Eqs. 8 and 9 significantly exceed the pressure of the TeV particles and magnetic field in the cocoon (Eqs.6 and 7).This is somewhat surprising, since the MHD flow in PWNe is expected to be subsonic (see, e.g., Kennel & Coroniti 1984b), i.e., to be nearly isobaric.Therefore, there should be another dominant contribution to the pressure.The energy required to support the pressure in the cocoon can be estimated as where we estimated the cocoon volume as V c 14× pc 2 (which corresponds to angular size of 0.4 • × 1.4 • at the source distance of 290 pc).Slane et al. (2018) obtained that the size of the Vela X PWN and SNR can be reproduced if one adopts an initial spin-down power of Ė0 = 10 39 erg s −1 and braking index n b = 3.For these parameters of the pulsar, the total energy injected to the PWN can be estimated as 3 × 10 49 erg.Although the size of the entire nebula is much bigger than the cocoon and particles may undergo significant energy losses due to radiative and adiabatic cooling, it is possible that relativistic particles could provide the missing pressure.
The pressure deficit cannot be explained by sub-TeV electrons in the cocoon, since Tibaldo et al. (2018) found that their spectrum is very hard.Alternatively, the missing pressure could be provided by relativistic electrons in the extended radio nebula.According to Grondin et al. (2013) their total energy is about 10 48 erg.These electrons are distributed in a region significantly larger than the cocoon (probably, a factor of ∼ 10 in volume), thus their contribution is not sufficient to explain completely the pressure deficit either.Finally, even mildly relativistic electrons with such energy density could violate the upper limits derived at 400 MHz on the radio brightness from the Vela X direction (Haslam et al. 1982;Grondin et al. 2013).
Since the radio data do not allow us to constrain a possible contribution to pressure from relativistic ions, they may be considered a candidate to explain the pressure deficit.Models of pulsar winds that include relativistic ions were presented, e.g., by Arons & Tavani (1994).Horns et al. (2006) proposed a model of gamma-ray emission from Vela X in which all the emission is produced by relativistic ions with a total energy of 10 49 erg for protons or 10 48 erg for iron nuclei.We note, however, that this scenario may be problematic when considering the energetics of relativistic ions and leptons over the whole PWN, including the extended radio nebula, as discussed in de Jager & Djannati-Ataï (2009).We also note that in this scenario the interpretation of our gamma and X-ray data would need to be significantly modified.
Alternatively, the pressure deficit could be considered as an indication of a significant contribution from non-relativistic ejecta material to the pressure in the cocoon region, i.e., by taking a small filling factor F ≤ 10 −2 in Equation 7.This is qualitatively consistent with the results from the hydrodynamic simulations of Slane et al. (2018) that show how the impact of the reverse shock also causes a pronounced mixing of the ejecta with young PWN material in the cocoon area (Fig. 13 of Slane et al. 2018).We note, however, that the study by Slane et al. (2018) was performed in the hydrodynamic limit, thus it does not allow robust conclusions regarding the properties of the magnetic field.In the future, this scenario should be revisited in the MHD framework.

Origin of the cocoon and particle transport
We consider two scenarios to explain the origin of the relativistic particles in the Vela X cocoon.In the first scenario particles diffuse from the acceleration site, which is conventionally associated with the pulsar wind termination shock, to the cocoon (e.g., Hinton et al. 2011).In this case the cocoon represents just a preferred particle leakage path due to some peculiar configuration of the magnetic field (see, e.g., the discussion in de Jager & Djannati-Ataï 2009).In the second scenario, particles are quickly advected from the acceleration site to the cocoon by the highly asymmetric reverse shock (Blondin et al. 2001;Slane et al. 2018).We note that proper modelling of these two processes requires detailed MHD/particle transport simulations, that are beyond the scope of this paper.For the following we just compare some of the results obtained from our analysis to general expectations for the two scenarios.
As it was shown above, the TeV particle and magnetic field energy densities are comparable in the cocoon, therefore particle diffusion must proceed in the non-linear regime.Nevertheless it is worthwhile to compare the required diffusion coefficient with the Bohm limit: where κ is a dimensionless scaling constant.The characteristic diffusion distance is then which is comparable to the size of the cocoon for E ∼ 1 TeV if κ ≥ 100.Although the diffusion should proceed much faster than in the Bohm limit, this value of the diffusion coefficient is still significantly smaller than typical ISM values (e.g., Strong et al. 2007), which is expected, to some extent, in media perturbed by an SN explosion or due to non-linear particle transport (Evoli et al. 2018).
Even for a weak magnetic field in the cocoon, synchrotron emission should be the dominant cooling channel.Thus, the cooling time can be estimated as In the diffusion scenario, the particles accelerated in the inner nebula may experience a strong magnetic field.The pressure near the termination shock estimated above corresponds to an equipartition magnetic field of B s 100 µG, which for 100 TeV electrons yields a cooling time of t syn ∼ 10 yr, which is very short Article number, page 11 of 19 as compared to the age of the source.However, particles may escape quickly from the high-magnetic field strength region (for ballistic propagation ultrarelativstic particles would take ∼1 yr to reach the lower-magnetic field strength region in pointing 0), and it cannot be excluded that the magnetic field strength is in fact weaker than the equipartition value.To assess if the diffusion scenario is viable the transport of the particles near the termination shock must be studied in detail, which is, however, beyond the scope of this paper (for some studies in this direction, see, e.g., Van Etten & Romani 2011;Porth et al. 2016;Ishizaki et al. 2018).Alternatively, if the cocoon was created by the reverse shock, then it would correspond to a region that was originally close to the termination shock and which was then swept away.Before the arrival of the reverse shock, the radius of the termination shock was significantly larger than at present, and consequently the magnetic field at the termination shock could be smaller.In this case since the formation of the cocoon proceeds by a hydrodynamic process, the particle transport occurs in an energy independent manner.For the observed magnetic field strength of B 6 µG and a time of 6000 yr, the age of the cocoon suggested by Slane et al. (2018), a synchrotron cooling feature appears at E cut 100 TeV, consistent with our results (Table 3).In this case, the reconstructed spectrum below 100 TeV still corresponds to the acceleration spectrum.Table 3 shows that the values of the power-law index are broadly consistent with the canonical value of α 2 − 2.2 expected at late times for Fermi acceleration at relativistic shocks (e.g., Achterberg et al. 2001;Sironi & Cerutti 2017), but uncertainties remain large.

Summary and Conclusions
We have combined X-ray data from three pointings of Suzaku (total observation time of 139 ks) with ∼ 100 hours of gammaray observations with H.E.S.S. to extract the SEDs in three compact (∼1 pc) regions of the Vela X cocoon covering distances from the pulsar that range from ∼ 0.3 pc to ∼ 4 pc (for a distance of 290 pc).The gamma-ray spectra are best modelled by power laws with exponential cutoffs in pointings 1 and 2 (farther from the pulsar).For pointing 0 (closer to the pulsar) the best-fit exponential-cutoff power law has a cutoff energy consistent within statistical uncertainties with the other two regions, although the presence of a cutoff is not statistically significant (< 3σ).The X-ray spectra are all well described by power laws modified by Galactic interstellar absorption.The X-ray spectral properties of regions 1 and 2 are consistent within uncertainties, while region 0, closer to the pulsar, shows a hint of harder emission, although uncertainties from the spill-over of emission outside our spectral extraction region are larger.
We have fit the X-ray and TeV SEDs with a simple radiative model with an electron population producing synchrotron and IC emission.This enabled us to reconstruct the electron spectrum and magnetic field properties with minimal modelling assumptions.For the electron spectral distribution we adopted a four-parameter family of power laws with sub/super-exponential cutoff.For the magnetic field we considered two different scenarios: a fixed strength, and a case in which the magnetic field also has a component with power-law strength distribution that aims to account for the contribution from turbulent fluctuations.The results did not favour the presence of a significant turbulence level: the constant magnetic field scenario with fewer parameters can satisfactorily reproduce the data.The parameters of the electron spectral distribution are consistent in the two cases within the uncertainties, which are particularly large for the exponential cutoff index.
TeV electrons and magnetic field are in a state close to energy equipartition in all three regions.The pressure of both components (∼10 −12 erg cm −3 ) is small compared to what is inferred from the properties of the wind termination shock and the SNR shell (∼10 −9 erg cm −3 ).This points to the existence of another dominant source of pressure.Sub-TeV electrons can hardly explain the pressure deficit completely, but a contribution from relativistic ions cannot be excluded based on the current data.An alternative explanation may be found in a significant contribution from non-relativistic matter due to mixing of PWN and ejecta (small filling factor for relativistic particles), consistent with recent hydrodynamical simulations of the system (Slane et al. 2018).
We compared the electron spectra with general expectation from particle-transport scenarios dominated by either diffusion or advection via the reverse shock.In the first case the diffusion coefficient needs to be 100 times larger than the Bohm limit, but still significantly smaller than typical ISM values.Significant radiative losses expected near the termination shock require to carefully study particle transport in this region in order to assess if the scenario is viable in the light of the measured cutoff energies of ∼100 TeV.Conversely, in the reverse-shock scenario the cutoff energies are consistent with radiative cooling in the cocoon, and the spectra below the cutoff, corresponding to the acceleration spectra, are broadly compatible with expectations for Fermi acceleration at relativistic shocks.
The constraints on magnetic field turbulence and the shape of the particle cutoff are very sensitive to the spectrum in the hard X-ray band.Thus, they can be improved with future observations in this domain, e.g., with NuSTAR (Harrison et al. 2013).They can also benefit from better gamma-ray measurements with CTA (Actis et al. 2011).

Fig. 1 .
Fig. 1.Regions used for spectral analysis (blue box and circles) overlaid onto two maps of the Vela X region.Left panel: X-ray count map from the ROSAT survey at energies > 0.5 keV.Right panel: significance map from H.E.S.S. at energies > 0.6 TeV (see Sec. 4 for details on how the map is derived, the map is oversampled with a correlation radius of 0 • .2 for display).The red star indicates the position of the Vela pulsar and its size approximately corresponds to that of the jet-torus structure measured with Chandra (Manzali et al. 2007).The dashed red circle represents the 95% containment radius of the Suzaku point spread function around the pulsar.The dashed green squares indicate the borders of the Suzaku field of view for the three pointings used in this paper.
3 pc), which corresponds to the 95% containment radius of the Suzaku XIS point spread function (PSF).The final spectral extraction region in Pointing 0 is therefore defined as a rectangle centred at R.A. = 128 • .81,Dec. = −45 • .286,with major edge of 17.4 , minor edge of 4.2 , and tilted by 157 • with respect to the R.A. axis.From 2004 to 2012 H.E.S.S. was operated as an array of four telescopes with 12 m diameter.In 2012 a fifth telescope with 28 m diameter was added, improving the performance at low energies.We use a dataset spanning from 2004 to 2016, comprising observations taken as part of the survey of the Galactic plane (H.E.S.S. Collaboration et al. 2018b), from studies of nearby sources such as the Vela pulsar (H.E.S.S. Collaboration et al. 2018c), Puppis A (H.E.S.S. Collaboration et al. 2015), and Vela Junior (H.E.S.S. Collaboration et al. 2016), as well as dedicated observations.

Fig. 2 .
Fig. 2. Suzaku X-ray count spectra measured in the regions as illustrated in Figure 1.In each panel, we show spectra obtained with XIS 0 (black), XIS 1 (blue), XIS 2 (red), and XIS 3 (green), superposing the best-fit model as solid lines.The residuals with respect to the best-fit model are shown in the lower panel.Error bars on the flux are 1σ statistical uncertainties.

Fig. 3 .
Fig. 3. SEDs of gamma-ray emission measured by H.E.S.S. for the three extraction regions/pointings shown in Figure 1.In the top panel the lines show the best-fit spectral models, and the bands display the statistical uncertainties based on quadratic error propagation from the inversion of the likelihood Hessian matrix.Points show the binned SEDs with their statistical uncertainties (coloured error bars with end caps) and the sum in quadrature of statistical and systematic uncertainties (grey error bars without end caps).The bottom panel shows the deviation of the SED points from the best-fit function as number of statistical σ.For Pointing 0 orange corresponds to the PL fit, and blue to the ECPL fit.

Table 1 .
Parameters of the spectral fits to X-ray data Region Photon index Flux a (10 −7 erg s −1 cm −2 sr −1 ) Notes.For fixed N H and CXB parameters, see text for details.The effect of interstellar absorption is negligible.All errors are reported as 1σ confidence intervals.