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

TL;DR

This study combines X-ray and gamma-ray observations of Vela X to analyze relativistic particles and magnetic fields, revealing electron energy cutoffs, magnetic field strengths, and pressure components in the nebula.

Contribution

It provides a multi-wavelength analysis of Vela X, deriving properties of particles and magnetic fields with minimal modeling, and constrains the electron spectrum and magnetic turbulence.

Findings

Electron spectra show a cutoff at ~100 TeV.

Magnetic field strength is approximately 6 microGauss.

TeV electrons and magnetic fields have negligible pressure, implying other dominant pressure sources.

Abstract

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)

Tables

Table 1. Table 1: Parameters of the spectral fits to X-ray data

Region	Photon index	Fluxa𝑎aa𝑎a$a$The flux within the extraction region is reported for the energy range $2\mathrm{keV}$ to $10\mathrm{keV}$. (${10}^{-7}\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-2}{\mathrm{sr}}^{-1}$)
0	$1.920\pm 0.014$	$13.86\pm 0.15$
1	$2.251\pm 0.013$	$5.72\pm 0.05$
2	$2.31\pm 0.02$	$5.66\pm 0.10$

Table 2. 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 ​ σ 1 𝜎 1\sigma 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 $>1\mathrm{TeV}$	$\Gamma$	$E_{\mathrm{cut}}$	TS
		(${10}^{-8}{\mathrm{cm}}^{-2}{s}^{-1}{\mathrm{sr}}^{-1}$)		(TeV)
0	PL	$3.6\pm {0.4}_{\mathrm{stat}}\pm {0.7}_{\mathrm{sys}}$	$1.75\pm {0.10}_{\mathrm{stat}}\pm {0.14}_{\mathrm{sys}}$	–	–
	ECPL	$3.5\pm {0.5}_{\mathrm{stat}}\pm {0.7}_{\mathrm{sys}}$	$1.2\pm {0.3}_{\mathrm{stat}}\pm {0.1}_{\mathrm{sys}}$	$16\pm 7_{\mathrm{stat}}\pm 6_{\mathrm{sys}}$	6.4
1	PL	$4.1\pm {0.2}_{\mathrm{stat}}\pm {1.0}_{\mathrm{sys}}$	$1.92\pm {0.06}_{\mathrm{stat}}\pm {0.10}_{\mathrm{sys}}$	–	–
	ECPL	$4.1\pm {0.3}_{\mathrm{stat}}\pm {1.5}_{\mathrm{sys}}$	$1.47\pm {0.16}_{\mathrm{stat}}\pm {0.16}_{\mathrm{sys}}$	$14\pm 5_{\mathrm{stat}}\pm 6_{\mathrm{sys}}$	13.4
2	PL	$4.6\pm {0.3}_{\mathrm{stat}}\pm {0.9}_{\mathrm{sys}}$	$1.84\pm {0.06}_{\mathrm{stat}}\pm {0.11}_{\mathrm{sys}}$	–	–
	ECPL	$4.6\pm {0.3}_{\mathrm{stat}}\pm {1.4}_{\mathrm{sys}}$	$1.27\pm {0.14}_{\mathrm{stat}}\pm {0.10}_{\mathrm{sys}}$	$12\pm 3_{\mathrm{stat}}\pm 6_{\mathrm{sys}}$	25.7

Table 3. 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 .

Pointing	$W_e$ ($>1\mathrm{TeV}$)a𝑎aa𝑎a$a$Total electron energy for particle energies $>1\mathrm{TeV}$. We note that the solid angle subtended by region 0 is 41% of that in regions 1 and 2.	$\alpha$b𝑏bb𝑏b$b$Parameters of the electron spectrum as defined in Equation 3.	$E_{\mathrm{cut}}$b𝑏bb𝑏b$b$Parameters of the electron spectrum as defined in Equation 3.	$\beta$b𝑏bb𝑏b$b$Parameters of the electron spectrum as defined in Equation 3.	$B$c𝑐cc𝑐c$c$Strength of the magnetic field.	BICd𝑑dd𝑑d$d$Bayesian information criterion, i.e. $k\cdot \ln n+{\chi }^2$, where $k$ is the number of parameters estimated from the model and $n$ is the number of data points.
	(${10}^{44}$ erg)		(PeV)		($\mu$G)
0	${0.7}_{-0.4}^{+1.1}$	${2.2}_{-0.6}^{+0.3}$	${0.3}_{-0.2}^{+4.0}$	${2.0}_{-1.4}^{+2.6}$	${8.6}_{-1.7}^{+2.8}$	14.8
1	${1.7}_{-0.6}^{+0.8}$	$1.8\pm 0.5$	${0.05}_{-0.03}^{+0.07}$	${0.9}_{-0.2}^{+0.5}$	${6.7}_{-0.7}^{+0.9}$	18.2
2	${1.6}_{-0.5}^{+0.7}$	${1.9}_{-0.6}^{+0.3}$	${0.11}_{-0.06}^{+0.03}$	${2.0}_{-0.8}^{+2.3}$	${5.4}_{-0.6}^{+0.8}$	17.3

Table 4. 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	$W_e$ ($>1\mathrm{TeV}$)a𝑎aa𝑎a$a$Total electron energy for particle energies $>1\mathrm{TeV}$. We note that the solid angle subtended by region 0 is 41% of that in regions 1 and 2.	$\alpha$b𝑏bb𝑏b$b$Parameters of the electron spectrum as defined in Equation 3.	$E_{\mathrm{cut}}$b𝑏bb𝑏b$b$Parameters of the electron spectrum as defined in Equation 3.	$\beta$b𝑏bb𝑏b$b$Parameters of the electron spectrum as defined in Equation 3.	$B_0$c𝑐cc𝑐c$c$Parameters of the magnetic field distribution as defined in Equation 5.	$a$c𝑐cc𝑐c$c$Parameters of the magnetic field distribution as defined in Equation 5.	BICd𝑑dd𝑑d$d$Bayesian information criterion, i.e. $k\cdot \ln n+{\chi }^2$, where $k$ is the number of parameters estimated from the model and $n$ is the number of data points.
	(${10}^{44}$ erg)		(PeV)		($\mu$G)
0	${0.7}_{-0.3}^{+0.4}$	${2.3}_{-0.6}^{+0.4}$	${0.24}_{-0.19}^{+0.30}$	${2.3}_{-1.7}^{+2.4}$	$11\pm 3$	${0.4}_{-0.3}^{+0.4}$	17.3
1	${1.3}_{-0.4}^{+0.6}$	${2.0}_{-0.6}^{+0.3}$	$0.09\pm 0.04$	${2.6}_{-1.3}^{+1.9}$	$7.6\pm 0.9$	${0.39}_{-0.14}^{+0.11}$	20.5
2	${1.4}_{-0.3}^{+0.7}$	$1.7\pm 0.4$	$0.08\pm 0.03$	${2.4}_{-1.1}^{+1.5}$	${6.7}_{-1.2}^{+1.0}$	${0.32}_{-0.15}^{+0.14}$	19.9

Table 5. Table 5: Counting statistics in the spectral extraction region (ON) and background evaluation region (OFF) integrated over the whole energy range > 0.6 ​ TeV absent 0.6 TeV >0.6\rm\,TeV , and significance of the gamma-ray excess for the three analysis regions.

	Pointing 0	Pointing 1	Pointing 2
$N_{\mathrm{ON}}$	262	724	768
$N_{\mathrm{OFF}}$	2825	3777	4560
$\alpha$	0.048	0.092	0.078
$N_{\mathrm{background}}$	135.9	348.0	355.1
$N_{\mathrm{excess}}$	126.1	376.0	412.9
Significance ($\sigma$)	9.3	16.6	18.0

Table 6. Table 6: Best-fit parameters for the radiative model of the X-ray and gamma-ray spectrum of the Vela X cocoon for the case of increased FIR intensities in the three regions shown in Figure  1 .

Pointing	$W_e$ ($>1\mathrm{TeV}$)a𝑎aa𝑎a$a$Total electron energy for particle energies $>1\mathrm{TeV}$. We note that the solid angle subtended by region 0 is 41% of that in regions 1 and 2.	$\alpha$b𝑏bb𝑏b$b$Parameters of the electron spectrum as defined in Equation 3.	$E_{\mathrm{cut}}$b𝑏bb𝑏b$b$Parameters of the electron spectrum as defined in Equation 3.	$\beta$b𝑏bb𝑏b$b$Parameters of the electron spectrum as defined in Equation 3.	$B$c𝑐cc𝑐c$c$Strength of the magnetic field.	BICd𝑑dd𝑑d$d$Bayesian information criterion, i.e. $k\cdot \ln n+{\chi }^2$, where $k$ is the number of parameters estimated from the model and $n$ is the number of data points.
	(${10}^{44}$ erg)		(PeV)		($\mu$G)
0	${0.6}_{-0.2}^{+0.5}$	${2.2}_{-0.6}^{+0.4}$	${0.3}_{-0.2}^{+2.0}$	${2.1}_{-1.6}^{+2.6}$	$9_{-2}^{+4}$	14.8
1	${1.4}_{-0.4}^{+0.6}$	${1.8}_{-0.5}^{+0.4}$	${0.05}_{-0.03}^{+0.06}$	${0.9}_{-0.3}^{+0.6}$	${7.3}_{-0.9}^{+1.1}$	18.1
2	${1.7}_{-0.5}^{+0.6}$	${2.0}_{-0.4}^{+0.3}$	${0.12}_{-0.05}^{+0.03}$	${2.3}_{-1.0}^{+2.2}$	${5.8}_{-0.8}^{+0.9}$	17.3

Equations

$$\frac{\mathrm{d}N}{\mathrm{d}E}=A{\left(\frac{E}{E_0}\right)}^{-\Gamma },$$

$$\frac{\mathrm{d}N}{\mathrm{d}E}=A{\left(\frac{E}{E_0}\right)}^{-\Gamma }\exp \left(-\frac{E}{E_{\mathrm{cut}}}\right),$$

$$\frac{\mathrm{d}{N}_e}{\mathrm{d}{E}_e}=A{\left(\frac{E_e}{E_0}\right)}^{-\alpha }\exp \left[-{\left(\frac{E_e}{E_{\mathrm{cut}}}\right)}^{\beta }\right],$$

$$\frac{{\lambda }_B}{{\lambda }_{\mathrm{ph}}}>\sqrt{\frac{ϵ}{4\pi e^2{n}_e}}\frac{eB}{m_e{c}^2}=\frac{ϵ}{m_e{c}^2}\sqrt{\frac{B^2}{4\pi ϵn_e}},$$

$$\mathrm{PDF}(B)=(1-a)\delta (B-B_0)+aC{B}^{-\zeta }\mathcal{H}(B-B_{\min })\mathcal{H}(B_{\max }-B).$$

$$w_{\text{textsc}b}=\frac{B^2}{8\pi }\simeq 1.4\times 10^{-12}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}.$$

$$w_{\mathrm{e}}=\frac{W_{\mathrm{e}}}{\pi r^2\ell \mathcal{F}}\simeq 2.7\times 10^{-12}{\left(\frac{\ell }{1\mathrm{pc}}\right)}^{-1}{\mathcal{F}}^{-1}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3},$$

$$p_{\mathrm{s}}\simeq \frac{\dot{E}}{6\pi R_{\mathrm{s}}^{2}c}\simeq 7\times 10^{-10}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}.$$

$$p_{\mathrm{shell}}$$

$$E_{\mathrm{c}}=3V_{\mathrm{c}}p\simeq 10^{48}\left(\frac{\ell }{1\mathrm{pc}}\right)\mathrm{erg},$$

$$D=\kappa D_{\mathrm{B}}\sim 10^{25}\kappa \left(\frac{E}{1\mathrm{TeV}}\right){\left(\frac{B}{6\text{upmu}G}\right)}^{-1}\mathrm{c}{\mathrm{m}}^2{\mathrm{s}}^{-1},$$

$$\Delta r\sim \sqrt{D\Delta t}\sim 0.5{\kappa }^{1/2}{\left(\frac{E}{1\mathrm{TeV}}\right)}^{1/2}{\left(\frac{B}{6\text{upmu}G}\right)}^{-1/2}{\left(\frac{\Delta t}{10^4\mathrm{yr}}\right)}^{1/2}\mathrm{pc},$$

$$t_{\mathrm{syn}}=3.5\times 10^5{\left(\frac{E}{1\mathrm{TeV}}\right)}^{-1}{\left(\frac{B}{6\text{upmu}\mathrm{G}}\right)}^{-2}\mathrm{yr}.$$

Full text

*Aims. *Since TeV emission is likely inverse-Compton emission of electrons, predominantly from interactions with the cosmic microwave background, 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.

*Methods. *We used data from the Suzaku XIS to derive the spectra from three compact regions in Vela X covering distances from $0.3\rm\,pc$ to $4\rm\,pc$ from the pulsar along the cocoon. We obtained gamma-ray spectra of the same regions from H.E.S.S. observations and fitted a radiative model to the multi-wavelength spectra.

*Results. *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\rm\,TeV and a magnetic field strength of \sim$$6\,\rm\upmu G. Constraints on the presence of turbulent magnetic fields are weak.

*Conclusions. *The pressure of TeV electrons and magnetic fields in the cocoon is dynamically negligible, requiring the presence of another dominant pressure component to balance the pulsar wind at the termination shock. Sub-TeV electrons cannot completely account for the missing pressure, which 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, but for the latter the role of radiative losses near the termination shock needs to be further investigated in the light of the measured cutoff energies. 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$.

††offprints: H.E.S.S. collaboration,

;

11footnotemark: 1 Spectra are only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/

22footnotemark: 2 Corresponding authors

33footnotemark: 3 Deceased

1 Introduction

Pulsars eject relativistic winds that 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 to what extent the PWNe contribute to the electron and positron component observed in cosmic rays still remain to be established (e.g. Amato, 2014).

Most of the recent progress in understanding 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 emissions are 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. For the production of TeV emission, the latter are predominantly the cosmic microwave background (CMB), and far-infrared (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 that have varying magnetic field strengths 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 with observations for a sufficiently extended and bright PWN, such as Vela X. Originally discovered as a $3^{\circ}\times 2^{\circ}$ 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 $\dot{E}\simeq 7\times 10^{36}\rm\,erg\,s^{-1}$ and characteristic age of $\tau=1.1\times 10^{4}\rm\,yr$, Manchester et al. 2005). This association places the PWN at a distance of $287^{+19}_{-17}\rm\,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)111The choice of name cocoon was motivated by the interpretation, now superseded, that this structure confines particles injected by a jet from the pulsar. H.E.S.S. is 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. In 2006 it unveiled the TeV emission associated with the cocoon (Aharonian et al., 2006a), making Vela X one of the archetypal TeV PWNe.

To explain multi-wavelength 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 $\sim 5\rm\,\upmu G$. Hinton et al. (2011) proposed that the extended radio nebula is filled with a relic electron population, devoid of high-energy particles ($>10\rm\,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 cannot, however, reproduce the details of the multi-wavelength 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 have taken 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^{\prime}$ and provides spectral coverage in the range from $0.2\rm\,keV$ to $12\rm\,keV$ with low background contamination, which enables us to probe emission from the highest-energy electrons in Vela X (for previous studies, see, e.g. Mattana et al., 2011; Katsuda et al., 2011). The same electrons are also responsible for the multi-TeV emission observed with H.E.S.S. indicating 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 non-uniform-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 multi-wavelength 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 multi-wavelength 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.

2 Observations and analysis region definition

We used three archival Suzaku observations of the Vela X cocoon conducted in 2006. From north (closer to the pulsar) to south, they have observation IDs 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 used circular spectral extraction regions with a radius of 7.5*′* centred in the middle of the Suzaku XIS field of view (R.A. = 128$\aas@@fstack{\circ}$7666, Dec. = $-45\aas@@fstack{\circ}4581$, and R.A. = 128$\aas@@fstack{\circ}$6368, Dec. = $-45\aas@@fstack{\circ}8007$, respectively). For Pointing 0 we defined a spectral extraction region such that we were able to avoid the region immediately adjacent to the pulsar. 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\rm\,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 excluded from the analysis a circular region centred on the pulsar with a radius of 3.6*′* ($0.3\rm\,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$\aas@@fstack{\circ}$81, Dec. = $-45\aas@@fstack{\circ}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 used 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.

We applied the following selection criteria to the dataset:

• •

observations for which at least three of the $12\rm\,m$ telescopes were operational to improve high-energy performance, since we are interested in gamma rays at energies $>1\rm\,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\aas@@fstack{\circ}5$ from the centre of Vela X, at R.A. = 128$\aas@@fstack{\circ}$75, Dec. = $-45\aas@@fstack{\circ}6$ (all coordinates are given in J2000 equinox in this paper).

3 X-ray analysis

We used the HEASoft package222https://heasarc.gsfc.nasa.gov/docs/software/heasoft/. 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^{\prime}$, which is sufficiently larger than the field of view and a standard setting of xissimarfgen for data analysis of diffuse emission. We evaluated instrumental background based on night Earth’s observations by using xisnxbgen. We used only photons with energies above $2.25\rm\,keV$ to avoid thermal emissions from the Vela SNR, and set an upper energy limit at $10\rm\,keV$ for front-illuminated CCDs (XIS 0, 2, 3) and to $7\rm\,keV$ for the back-illuminated CCD (XIS 1) to minimise effects from particle-induced instrumental background.

Each X-ray spectrum is well described by a single power law modified by a fixed Galactic interstellar absorption $N_{\mathrm{H}}=2.59\times 10^{20}\rm\,cm^{-2}$ (Manzali et al., 2007) plus underlying cosmic X-ray background (CXB), parametrised as in Miyaji et al. (1998), that is, a power law with a photon index of 1.42 and a normalisation of $10.0\rm\,s^{-1}\,cm^{-2}\,sr^{-1}\,keV^{-1}$ at $1\rm\,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\rm\,keV$.

The fit results show no significant differences between Pointings 1 and 2 in the flux level or spectral slope ($\sim 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^{\prime}$ from the pulsar, that is, beyond the bright TeV cocoon. Conversely, in Pointing 0 we find a significantly harder spectrum (photon index $1.92\pm 0.014$) and the emission is a factor of greater than two 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 ten energy bins between $2.25\rm\,keV$ and $9.5\rm\,keV$. We included in the subsequent multi-wavelength analysis a 10% systematic uncertainty on flux measurements (Sekiya et al., 2016). Furthermore, 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 estimated 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.

4 Gamma-ray analysis

Since we are mainly interested in the energy range $>1\rm\,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 used tight selection criteria to make the residual contamination from hadrons misclassified as gamma rays as low as possible, in other words, we applied the ‘hard cuts’ from Ohm et al. (2009). Additionally, to ensure uniform energy thresholds between different observing periods, we included in the subsequent analysis only candidate gamma-ray events with reconstructed energy $>0.6\rm\,TeV$, thus reducing systematics in spectral reconstruction at the lower end of the energy range.

All the results in the paper were cross-checked using a second independent software and 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. Also in this case we used ‘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.

The right-hand panel of Figure 1 shows the gamma-ray detection significance for the whole region of the Vela X cocoon. The map was 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 was estimated from the data within the same observation, excluding regions of the field of view in which significant gamma-ray emission is potentially present. We specifically excluded from the background estimation regions that had gamma-ray emission detected at significance $>5\sigma$ in the H.E.S.S. Galactic plane survey (H.E.S.S. Collaboration et al., 2018b, 3.2.2), or for which 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 rays444The second criterion increases the area of the excluded region by $<5\%$.. Figure 1 shows that there is very significant gamma-ray emission at energies $>0.6\rm\,TeV$ in all three analysis regions.

We evaluated 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 fitted 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\rm\,TeV$, except for one single event at a reconstructed energy of $\sim 90\rm\,TeV$ in the spectral extraction region (no events with reconstructed energies $>60\rm\,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 $\sim 90\rm\,TeV$ event in the spectral analysis or not. We report in the following the results obtained by selecting only events with reconstructed energy $<60\rm\,TeV$ in Pointing 0.

We considered 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

[TABLE]

where $A$ is the differential flux that normalises the spectrum at $E_{0}$, and $\Gamma$ is the spectral index of the distribution. Alternatively we considered an exponentially-cutoff power law (ECPL)

[TABLE]

with the characteristic energy of the cutoff, $E_{\mathrm{cut}}$, as additional parameter. We used Eqs. (1) and (2) simply 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 $\mathcal{L}_{\mathrm{PL}}$ and $\mathcal{L}_{\mathrm{ECPL}}$ be the maximum likelihood values for the PL and ECPL models, respectively, given the source region and background regions counts. The test statistic $\mathrm{TS}=2\log(\mathcal{L}_{\mathrm{ECPL}}/\mathcal{L}_{\mathrm{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 applicable555Eq. 2 reduces to Eq. 1 in the limit $E_{\mathrm{cut}}\rightarrow+\infty$, therefore the null hypothesis (PL) lies on the border of the parameter space for the test hypothesis (ECPL)., a larger TS 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 selected the ECPL as best-fit model. For Pointing 0 TS is 6.4 formally corresponding to a significance of the cutoff at $2.5\sigma$ statistical level. The presence of a cutoff in Pointing 0 is not significant. However, we note that the value of the cutoff energy, $16\pm 7\rm\,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, reconstruction, 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 each other, with the exception of the cutoff energy in Pointing 2, which differs by $\sim 3\sigma$. 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\sigma$, and we require that each bin has at least two excess666The excess is defined as the number of events in the spectral extraction region minus the background estimate derived from the events in the background regions. 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, which is strongly supported by 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 Section 5 we have adopted 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.

5 Radiative modelling of the multi-wavelength spectral energy distribution

To study the obtained X- and gamma-ray spectra in a consistent way, we computed 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 or super-exponential cutoff:

[TABLE]

where $N_{e}$ is the electron number, $E_{e}$ is the electron energy, $A$ is the normalisation factor for the distribution, $\alpha$ is the spectral index of a power-law distribution with reference energy $E_{0}$, $E_{\mathrm{cut}}$ is the cutoff energy, and $\beta$ 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 multi-wavelength 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).

The naima package enables us to derive the best-fit values and posterior probability distributions of the model parameters given the SED points from the $\chi^{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 full-forward 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 scanned the cutoff energy $E_{\mathrm{cut}}$ in logarithmic space, so that the fit parameter is actually $\log_{10}(E_{\mathrm{cut}}/1\,\mathrm{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 X-ray 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\rm\,pc$ up to $4\rm\,pc$ from the pulsar wind termination shock. As expected for IC radiation in the Thomson regime, the values of $\alpha$ in Table 3 and those of $\Gamma$ in Table 2 are consistent with the equation $\Gamma=(1+\alpha)/2$. The cutoff exponent $\beta$ 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 necessarily produced in completely homogeneous zones. The spectrum extraction regions span distances of $\sim 1\rm\,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, and that 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, $\lambda_{B}$, to the characteristic photon formation length, $\lambda_{\mathrm{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 inhomogeneities:

[TABLE]

where $n_{e}$ and $\epsilon$ are electron number density and mean energy per electron, respectively. Since $\epsilon/(m_{e}c^{2})\sim\Gamma_{\textsc{PW}}\gg 1$ ($\Gamma_{\textsc{PW}}$ is the bulk Lorentz factor of the pulsar wind) and ${B^{2}/(4\pi\epsilon n_{e}})\sim 1$ (for energy equipartition between magnetic fields and particles expected in PWNe, Kennel & Coroniti 1984b), one obtains that $\lambda_{B}/\lambda_{\mathrm{ph}}\gg 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 gauged magnetic turbulence in Vela X against the multi-wavelength SEDs by assuming that the probability density function (PDF) corresponding to a magnetic field strength $B$ is:

[TABLE]

The first term of the summation represents a magnetic field with a fixed strength, i.e. $B_{0}$, with $\delta$ being the Dirac distribution. The second term in the summation represents a magnetic field with power-law distribution of index $\zeta$ between minimum and maximum strengths $B_{\mathrm{min}}$ and $B_{\mathrm{max}}$, respectively, with $\mathcal{H}$ being the Heaviside step distribution. The parameter $a$ represents the mixing between the constant and power-law magnetic field components.

For illustration we take $B_{\mathrm{max}}=100\times B_{0}$, that is, sufficiently exceeding the mean value to modify the synchrotron spectrum, and $\zeta=3/2$. Although this is a typical power-law index of interstellar MHD turbulence (Kraichnan-type turbulence), we note that the parameter $\zeta$ in our formulation has a different physical meaning, related to the distribution of the magnetic-field strength. The effect of choosing a different power-law index, for example $\zeta=5/3$, has an impact on the results that is negligible compared to other sources of uncertainty.

The parameters $C$ and $B_{\mathrm{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_{\mathrm{RMS}}=B_{0}$). The only free parameters are then the mixing $a$ and the root mean square magnetic field expectation value, $B_{0}$, which can be fitted to the multi-wavelength 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. This is largely unsurprising, since the magnetic field turbulence predominantly affects the high energy part, in other words, a region above SED peak which is not well covered by current X-ray measurements. This also results in large uncertainties on the cutoff index $\beta$, which is the electron spectrum parameter most coupled with the effect of turbulent magnetic fields. This prevents us from studying in 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).

6 Discussion

6.1 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 $\simeq 6\,\rm\upmu G$ the energy density is

[TABLE]

On the other hand, for a total energy of electrons with energies $>1\rm\,TeV$ of $W_{\rm e}\sim 10^{44}\rm\,erg$ (see Table 3), the particle energy density is

[TABLE]

where we have taken $r=0.63\rm\,pc$ (which corresponds to the spectrum extraction radius of $7.5^{\prime}$ at the source distance of $290\rm\,pc$), and $\ell$ is the rather uncertain size of the emitting region along the line of sight. Another rather uncertain parameter, ${\cal F}$, is the filling factor that determines the fraction of the volume filled by relativistic electrons. For large filling factors, ${\cal F}\simeq 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{\rm\,pc}\times 7{\rm\,pc}$ (which corresponds to angular size of $0.4^{\circ}\times 1.4^{\circ}$ at the source distance of $290\rm\,pc$), thus it is feasible that $\ell\gtrsim 1{\rm\,pc}$. Unless $\ell$ is large, $\ell\gg 1\rm\,pc$, the TeV particle pressure coincides within a factor of a few with the magnetic field pressure in the cocoon for ${\cal F}\simeq 1$.

Observations of the jet-torus structure in the inner nebula by Chandra point to a size of the wind termination shock $R_{\rm s}\simeq 1.3\times 10^{17}\rm\,cm$ (i.e. $30^{\prime\prime}$ for a distance of $290\rm\,pc$, where $30^{\prime\prime}$ is the projected radius of the X-ray arc in Fig. 2 of Helfand et al., 2001). For the pulsar spin-down luminosity of $\dot{E}\simeq 7\times 10^{36}\rm\,erg\,s^{-1}$ the Rankine-Hugoniot relations for a weakly magnetised ultrarelativistic pulsar wind (Kennel & Coroniti, 1984a) yield a total pressure at the termination shock of

[TABLE]

On the other hand, the dynamic pressure inside the SNR shell can be obtained under the assumption that it is in the Sedov-Taylor phase of its evolution by combining Eqs. 27 and 28 from Khangulyan et al. (2018):

[TABLE]

where $R_{\textsc{snr}}$ and $t_{\textsc{snr}}$ are the SNR radius and age, and $n_{\textsc{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), in other words, 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

[TABLE]

where we estimated the cocoon volume as $V_{\rm c}\simeq 14\times\ell{\rm\,pc^{2}}$ (which corresponds to angular size of $0.4^{\circ}\times 1.4^{\circ}$ at the source distance of $290\rm\,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 $\dot{E}_{0}=10^{39}\rm\,erg\,s^{-1}$ and braking index $n_{\rm b}=3$. For these parameters of the pulsar, the total energy injected to the PWN can be estimated as $3\times 10^{49}\rm\,erg$. Although 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 $\simeq 10^{48}\rm\,erg$. These electrons are distributed in a region significantly larger than the cocoon (probably, a factor of approximately ten 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\rm\,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, for example, 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 require significant modification.

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, that is, by taking a small filling factor ${\cal F}\leq 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.

6.2 Origin of the cocoon and particle transport

We have considered 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 and particle transport simulations, that are beyond the scope of this paper. In the following we compare some of the results obtained from our analysis to general expectations for the two scenarios.

As 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:

[TABLE]

where $\kappa$ is a dimensionless scaling constant. The characteristic diffusion distance is then

[TABLE]

which is comparable to the size of the cocoon for $E\sim 1\rm\,TeV$ if $\kappa\geq 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

[TABLE]

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_{\rm s}\simeq 100\rm\,\upmu G$, which for $100\rm\,TeV$ electrons yields a cooling time of $t_{\rm syn}\sim 10\rm\,yr$, which is very short as compared to the age of the source. However, particles may escape quickly from the high-magnetic field strength region (for ballistic propagation ultra-relativstic particles would take $\sim$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\simeq 6\rm\,\upmu 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_{\rm cut}\simeq 100\rm\,TeV$, consistent with our results (Table 3). In this case, the reconstructed spectrum below $100\rm\,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 $\alpha\simeq 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.

7 Summary and conclusions

We have combined X-ray data from three pointings of Suzaku (total observation time of 139 ks) with $\sim 100$ hours of gamma-ray observations with H.E.S.S. to extract the SEDs in three compact ($\sim$1 pc) regions of the Vela X cocoon covering distances from the pulsar that range from $\sim 0.3\rm\,pc$ to $\sim 4\rm\,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\sigma$). 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 suggestion of harder emission, although uncertainties from the spill-over of emission outside our spectral extraction region are larger.

We have fitted 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 or super-exponential cutoffs. 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.

Magnetic field and TeV electrons are in a state close to energy equipartition in all three regions. The pressure of both components ($\sim$10*-12* erg cm*-3*) is small compared to what is inferred from the properties of the wind termination shock and the SNR shell ($\sim$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 $\sim$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 could be improved with future observations in this domain, for example, with NuSTAR (Harrison et al., 2013). They may also benefit from better gamma-ray measurements with CTA (Actis et al., 2011).

Acknowledgements.

The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Alexander von Humboldt Foundation, the Deutsche Forschungsgemeinschaft, the French Ministry for Research, the CNRS-IN2P3 and the Astroparticle Interdisciplinary Programme of the CNRS, the U.K. Science and Technology Facilities Council (STFC), the IPNP of the Charles University, the Czech Science Foundation, the Polish National Science Centre, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Innsbruck University, the Austrian Science Fund (FWF), and the Austrian Federal Ministry for Science, Research and Economy, the University of Adelaide and the Australian Research Council, the Japan Society for the Promotion of Science and by the University of Amsterdam. We appreciate the excellent work of the technical support staff in Berlin, Durham, Hamburg, Heidelberg, Palaiseau, Paris, Saclay, and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation. This research has made use of data obtained from the Suzaku satellite, a collaborative mission between the space agencies of Japan (JAXA) and the USA (NASA). We acknowledge the usage of the following open-source software: aplpy (Robitaille & Bressert, 2012), astropy (Astropy Collaboration et al., 2018), corner (Foreman-Mackey, 2016), emcee (Foreman-Mackey et al., 2013), matplotlib (Hunter, 2007), naima (Zabalza, 2015).

Appendix A Additional Information from the H.E.S.S. spectral fitting

Appendix B Contributions from infrared fields to inverse-Compton emission from Vela X

The interpretation of the SED of the Vela X cocoon has often taken into account the CMB alone as target radiation field for IC scattering of the accelerated electrons (e.g. Aharonian et al. 2006a). In Figure 5 we show the contributions from different target radiation fields that approximate the interstellar radiation field observed near the Sun to the IC emission from a source with spectrum analogous to Vela X.

Due to the Klein-Nishina suppression IC emission from higher-energy target photons is less important. The NIR component in Figure 5 contributes $<2\%$ to the total IC emission at any energies, that is, direct starlight is unimportant. On the other hand, the FIR component in Figure 5 accounts for $>10\%$ of the total IC emission at any energies. Moreover, if a cutoff of the electron spectrum at few tens of TeV is present, as in Vela X, the competition between the Klein-Nishina suppression and the energy of IC photons upscattered from higher-energy target photons being higher can make IC emission from the FIR field prevail over that from the CMB. This shows that accounting for thermal emission from dust can be important, depending on the exact spectrum of the electron population (see also Aharonian et al. 1997).

In our analysis we adopted the FIR model by Popescu et al. (2017). This should be considered as a large-scale average, and therefore close to a lower limit on the intensities of the FIR fields at the position of our source. We verified that individual IR sources in the IRAS Catalogue (Helou & Walker 1988) contribute negligible fluxes even under the hypothesis that they all lie at the same distance from the Earth as Vela X. However, the diffuse fluxes may be larger due to local interstellar structures not accounted for in the model by Popescu et al. (2017). We have assessed the impact on the results of an increase of a factor of two of the FIR energy density on the fit parameters. The results are reported in Table 6.

The comparison with Table 3 shows that the effect is smaller than other uncertainties in the parameters (in agreement with the general conclusion by Aharonian et al. 1997).

Appendix C Probability density distributions of radiative model parameters fit to the multi-wavelength data