© The Authors 2025

1 Introduction

Pulsars (PSRs) are one of the most powerful particle accelerators of our Galaxy and give rise to highly energetic phenomena that can be observed up to the very-high-energy (VHE, 50 GeV ≲ E < 100 TeV) domain (The VERITAS Collaboration 2011; H.E.S.S. Collaboration 2023). Goldreich & Julian (1969) proved that pulsars are surrounded by the magnetosphere, a region filled with charged particles that can be accelerated with large Lorentz factors through the open magnetic field lines. Later on, Rees & Gunn (1974) showed that the accelerated particles escaping the magnetosphere interact with the surrounding medium, producing the pulsar wind nebula (PWN) via synchrotron emission.

During their evolution, pulsars, especially older ones, may travel distances greater than the radius of their original PWN, eventually escaping from it. When this occurs, the original PWN begins to evolve independently of the parent pulsar, giving rise to the so-called relic PWN, while the escaped pulsar generates a new PWN (van der Swaluw et al. 2004). In some cases, neutron stars may acquire a kick velocity during the supernova explosion and, when old enough, can also escape the parent supernova remnant (SNR), giving rise to the phenomenon of bow-shock PWNe (BSPWNe). These PWNe are produced by pulsars travelling with a supersonic speed inside the interstellar medium (ISM) after exiting the SNR (Kargaltsev et al. 2017; Olmi 2023), and they exhibit a variety of different morphologies.

Pulsar wind nebulae are observed from radio to X-rays through their synchrotron emission, and they can generate gamma rays through inverse Compton scattering (ICS) of the ambient photons up to the tera-electronvolt regime and even above 100 TeV, at ultra-high-energies (UHEs; de Oña Wilhelmi et al. 2022; Breuhaus et al. 2022). Relic PWNe are expected to be X-ray faint, because they host a population of cooling old electrons, but still observable at VHEs thanks to ICS emission (de Jager & Djannati-Ataï 2009). Many PWNe have been detected at VHEs and UHEs, hence they have been considered one of the possible candidates for the origin of leptonic ‘PeVa-trons’, i.e. the Galactic sources responsible for the production of electrons and positrons at energies above 1 PeV. On the other hand, BSPWNe are commonly detected as Hα nebulae or in the radio and X-ray bands (Pellizzoni et al. 2004; Kargaltsev et al. 2017).

A new class of sources known as pulsar halos (or TeV halos) has been recently established. They were discovered in 2017 (Abeysekara et al. 2017) as regions of degree-wide teraelectronvolt emission around the Geminga and Monogem PSRs, which could not be attributed to the associated PWNe due to their significantly larger extension of approximately tens of parsecs. They are thought to have originated via ICS between cosmic microwave background (CMB) photons and VHE electrons and positrons accelerated by the central PSR that have escaped into the ISM. Following the discovery, Giacinti et al. (2020) defined a pulsar halo as the third stage of the evolution of a PSR-PWN system. The first stage is the free expansion of the PWN, lasting until the age of the PSR is t_PSR ≲ 10 kyr. This is followed by the evolution of the system after the SNR’s reverse shock encounters the PWN, which usually lasts up to tp_SR ≲ 100 kyr. In both stages, the pulsar-originated particles are still confined in the parent SNR, or only a few of them can escape, and no pulsar halo can be observed. VHE particles are free to travel within the ISM and generate halos only when the PSR becomes middle-aged (i.e., t_PSR ≳ 100 kyr), either due to its escape from the parent SNR or to the fading of the SNR itself. Adopting this definition, halos are identified as systems where the PSR, which generated the relativistic particles, no longer dominates the ISM dynamics or composition anymore (Giacinti et al. 2020). In particular, the energy density can be used as a potential phenomenological indicator to distinguish between halos and tera-electronvolt PWNe by comparing it to that of the ISM. As a consequence, only the extended emission detected by HAWC around Geminga and Monogem (Abeysekara et al. 2017) and that detected by LHAASO around PSR J0622+3749 (Aharonian et al. 2021) can be classified as a proper pulsar halo. However, a second definition has been proposed by Sudoh et al. (2019) and describes a halo as a region hosting relativistic particles where diffusion dominates the particle transport. This definition is independent of the age of the parent PSR and implies that halos could also be found in younger systems. For this reason, many other candidate sources have been proposed in the latest VHE catalogues (Albert et al. 2020; Cao et al. 2024).

The first catalogue published by the Large High Altitude Air Shower Observatory, LHAASO (Cao et al. 2024), reported many new unidentified sources detected at both VHEs and UHEs. Of the listed 90 sources, 43 were detected above 4σ at energies larger than 100 TeV, 25 were discovered for the first time, and 35 are located in the vicinity of pulsars. 1LHAASO J1740+0948u belongs to all the previous categories. It has been detected only by the KM2A instrument, i.e. above 25 TeV, with a significance of ~12σ (TS = 156) that decreases to ~6σ (TS₁₀₀ = 37) above 100 TeV. Recent results by LHAASO Collaboration (2025) obtained with additional observation time on the source reported an improved detection at 17.1σ and 9.4σ above 25 and 100 TeV, respectively. They derived a 95% confidence level upper limit on the extension of 0.147° and a log-parabola best-fit spectrum showing a peak at -30 TeV and extending up to ~300 TeV. The closest pulsar to 1LHAASO J1740+0948u reported in the Australia Telescope National Facility (ATNF) Pulsar Catalogue¹ (Manchester et al. 2005) is PSR J1740+1000, with an angular separation of 0.22°. Note that the 95% confidence level positional uncertainty on 1LHAASO J1740+0948u is -0.02° on both right ascension and declination. Despite this offset, the absence of other sources near the LHAASO coordinates makes the association with the pulsar likely.

PSR J1740+1000 (J1740 hereafter) is a radio pulsar discovered during a 430 MHz survey of the Arecibo telescope (McLaughlin et al. 2000) and later detected in X-rays by Chandra and XMM-Newton (Kargaltsev et al. 2008a; Rigoselli et al. 2022), as well as in the gamma-rays by Fermi-LAT (Abdollahi et al. 2022). It has a period of P = 154 ms and a period derivative of P˙ = 2.1 · 10⁻¹⁴ s⁻¹, from which the characteristic age (τ_C = P/2P˙ = 114 kyr) and spin-down power (Ė = 2.3 · 10³⁵ erg s^-1) can be obtained. The distance estimated through the most updated dispersion measure data is d = 1.23 kpc (Yao et al. 2017). The source is located significantly above the Galactic plane (b = 20.4°) and McLaughlin et al. (2002) suggested a halo origin of the progenitor star, a hypothesis reinforced by the X-ray detection of a long ~5-7′ cometary tail extended in the south-western direction (Kargaltsev et al. 2008a).

Due to its X-ray tail, J1740 can be considered an example of a BSPWN. Benbow et al. (2021) attempted to detect its tera-electronvolt counterpart, as well as for a few other X-ray BSPWNe, using the Very Energetic Radiation Imaging Telescope Array System (VERITAS) telescopes. In their work, they analysed XMM-Newton observations of J1740 in the range 0.310 keV and found that the emission is best fitted by a power law with spectral index Γ = 1.75 ± 0.04. For the VHE counterpart, instead, only upper limits were obtained. They then developed a model to describe the synchrotron losses of the electron distribution as the particles propagate along the tail, which was assumed to have a cylindrical geometry. They fitted the model to match the X-ray observations and used it to derive the spectral energy distribution (SED) of the corresponding VHE-emitting population. According to their modelling, with improved flux sensitivity of tera-electronvolt instruments, it could become possible to detect a counterpart of the X-ray tail of J1740.

Considering its offset with respect to the pulsar position, LHAASO Collaboration (2025) argued that 1LHAASO J1740+0948u may not be the exact counterpart of the X-ray tail, but it could represent the VHE detection of re-accelerated electrons coming from the bow-shock pulsar tail. In particular, they performed a phenomenological model of the synchrotron and ICS losses from a leptonic population and concluded that, for magnetic field values of B = 3, 5 μG, an X-ray counterpart at the location of 1LHAASO J1740+0948u could be detected with appropriate exposure time. They also explored different scenarios, including a hadronic hypothesis, which was excluded due to the absence of molecular clouds in the surroundings of the pulsar location, and a pulsar halo scenario, which was disfavoured due to the small size and offset of the LHAASO source and J1740. However, they claim that projection effects or very slow diffusion could possibly still be compatible with a pulsar halo origin of the source.

Another study on 1LHAASO J1740+0948u was carried out by Xie et al. (2025). They collected multi-wavelength data at the source location, including Planck, Swift-XRT, Fermi-LAT, CfA ¹²CO, and IceCube to complement the LHAASO spectrum reported in Cao et al. (2024). No counterpart has been found at other wavelengths. They then performed a leptonic SED modelling, assuming the source is powered by a Geminga-like PWN, and derived the main parameters of the relativistic particle population.

Due to their interplay with VHE emission and the better resolution of the instruments when compared to the gamma-ray ground-based facilities, X-rays are an important tool to help us better understand the nature of both PWNe and halos. In particular, since pulsar halos are produced via ICS of the ambient photon fields by the highly relativistic particle wind generated from the pulsar, we would expect to detect synchrotron radiation from the same population of particles and observe it in the X-ray band (Liu et al. 2019a; Li et al. 2022). No counterpart of a pulsar halo in the X-rays has been detected so far, but some attempts have been performed, see for example Khokhriakova et al. (2024); Manconi et al. (2024); Adams et al. (2025). In a recent work by Niu et al. (2025), the authors claimed the detection using eROSITA of extended X-ray emission around PSR B0656+14, i.e. the Monogem pulsar, but independent work by Khokhriakova et al. (2025), which adopted the most up-to-date eROSITA point spread function (PSF) models, proved that the soft diffuse emission detected by Niu et al. (2025) can be modelled as the Monogem pulsar emission, suggesting the diffuse emission is due to PSF leakage instead of an additional component.

In this work, we aim to characterise for the first time the diffuse X-ray emission around PSR J1740+1000 using deep XMM-Newton observations in order to investigate the presence of a putative counterpart of the VHE source 1LHAASO J1740+0948u. We also revisit the analysis of the tail-like PWN, performing the first study of the spectral index as a function of the distance from the pulsar. We then perform a multiwavelength SED fitting, using both the X-ray fluxes derived from the spectral analysis of the tail and the upper limits obtained from the diffuse emission. We also used the newly available teraelectronvolt spectrum from LHAASO Collaboration (2025) to test different possible acceleration scenarios that may explain the nature of the dark VHE source. The paper is divided as follows: in Sect. 2 we present the methodology followed for the X-ray analysis of both the diffuse emission and the tail PWN, in Sect. 3 we show the derived results and the SED fitting, and in Sect. 4 we discuss the results interpretation.

2 X-ray observations and data analysis

The European Photon Imaging Cameras (EPIC) on board XMM-Newton observed PSR J1740 between 2017 and 2018 for a total of five observations (ObsIDs: 0803080101–0803080501) and twice in 2006 (ObsIDs: 0403570101, 0403570201). We focused on the 2018 sample due to the large total exposure time compared to the 2006 sample. Out of the five observations, the first one (ObsID: 0803080101) was not considered due to both EPIC-MOS1/2 and EPIC-pn instruments being affected by high background from solar flares. The exposure time for the four observations together reaches a total of 532.5 ks. The EPIC-pn instrument was operated in ‘Small Window’ mode, i.e. with only a part of CCD4 used to collect data, while both EPIC-MOS1 and MOS2 were operated in ‘Full Frame’ mode, i.e. all CCDs in read-out mode (ESA: XMM-Newton SOC 2024). Moreover, CCD3 and CCD6 of the MOS1 have not been operational since 2005 and 2012, respectively, due to impact with micrometeorites (ESA: XMM-Newton SOC 2024). For these reasons, in the analysis of the diffuse emission around the pulsar, we used only data collected by the MOS2 camera that recorded the full field of view (FoV) available to the camera. The same holds for the analysis of the tail, which, for some observations, was not fully covered by the available MOS1 chips. We performed the data reduction using version v21 of the Science Analysis System (SAS)², but used two different approaches for the analysis of the diffuse emission and the PWN.

2.1 Analysis of the diffuse emission

To study the diffuse emission, we employed the XMM-Newton Extended Source Analysis Software (XMM-ESAS) pack-age³ available in SAS. The pipeline we followed for the analysis started with the raw data reduction with emchain. Then, we filtered for the soft proton flares using the histogram method for the task espfilt, setting a flare-free region at 3σ around the peak of the Gaussian distribution of the count rates. The next step was to run the task cheese for the source detection in the entire FoV of MOS2, for which we used a minimum maximum likelihood detection value of mlmin=5 (corresponding to a significance of ~4σ). This task produced a mask for the selected background sources that was later used in the spectral extraction phase.

The goal of this part of the analysis was to detect diffuse nonthermal X-ray emission potentially associated with 1LHAASO J1740+0948u. The main challenge we encountered when selecting the source region was that the full extension of 1LHAASO J1740+0948u is not covered by the XMM-Newton FoV. For this reason, we investigated two regions. The first, dubbed ‘Region 1’ hereafter, was selected as a r = 0.11° circle centred on the pulsar and excluding an ellipse comprising both the pulsar and the PWN. This region covers part of the extension of the LHAASO source and approximately the entire area of the central CCD, where the PSF of XMM-Newton has minimal distortions when compared to the on-axis PSF. The second, ‘Region 2’ from now on in the paper, was selected as a circular r = 0.06° region in the outer part of the XMM FoV, corresponding to the majority of the 68% source extension reported by LHAASO Collaboration (2025). Both regions are shown in panel (b) of Fig. 1, along with the extension of 1LHAASO J1740+0948u.

Finally, we performed the spectral extraction using mosspectra in the energy range 0.5-10 keV. This task also generated the corresponding response matrix file (RMF) and ancillary response file (ARF). Note that we did not run the task mosback because we wanted to fit the background, both instrumental and physical components, with a custom-made model. We adopted this approach as our goal was to investigate the presence of an extended component that may extend beyond the FoV of the instrument.

The whole procedure was repeated for each observation individually. Once the final spectrum with its RMF and ARF were obtained for all four observations, we stacked them to obtain the total spectrum. To do so, we used some of the tasks of the ftools package: mathpha, to stack the source spectra, and addrmf and addarf to produce an average RMF and ARF. Computing an average response is necessary to take into account the differences between the individual pointings. The final source spectrum was grouped using grppha with a minimum of 200 counts per bin. The spectral modelling and fitting of both the source and background were performed with version v12.13.1 of the Xspec package (Arnaud 1996).

Fig. 1 (a) ROSAT all-sky map showing the large-scale diffuse emission of the North Polar Spur and the location of the XMM-Newton pointings, i.e. the small green circle. The red, green and blue colours represent the 0.11-0.28 keV, 0.47-1.21 keV and 0.76-2.04 keV energy bands, respectively. (b) Image of the XMM-Newton FoV, showing both Region 1 and Region 2 (solid circles), and 1LHAASO J1740+0948u (dashed circles). The plotted radii of the source are taken from the 1 σ extension of the Gaussian (RG,LHAASO) and the diffuse (RD,LHAASO) models for photon energy above 25 TeV, shown in LHAASO Collaboration (2025). The 0.5-1 keV band is shown in red, the green emission corresponds to the range 1-2 keV, and the range 2-4.5 keV is represented in blue. Note that the observed point sources were later removed during the diffuse emission analysis.

2.1.1 Background modelling

The instrumental background of XMM-Newton consists of several components: the quiescent particle background (QPB), the contamination by soft proton flares (SPF) and the contamination by solar wind charge exchange (SWCX). The first was already treated during the data reduction, so we focused on the QPB component that can be described as a combination of a flat continuum emission and instrumental lines. We followed the prescriptions of Leccardi & Molendi (2008) and modelled the continuum as a broken power law (bknpower) with Γ₁ = 0.32, Γ₂ = 0.22 and E_b = 3 keV. Being this component intrinsic to the instrument, we did not convolve it with the ARF. We performed the fit by leaving all the parameters of the broken power law, except the normalisation, free. We then added the instrumental lines as Gaussian components (gauss) one by one to verify whether the new component was significant. Out of all the lines listed in Leccardi & Molendi (2008) we included: Al Kα, Si Kα, Cr Kα, Mn Kα, Cr Kβ, Fe Kα, Ni Kα and Au Lα. We left the energy and normalisation of all lines free to vary and fixed the width to zero for all except the Al Kα and Si Kα. We also verified the presence of SWCX contamination in the dataset. We found no significant flux enhancement in the soft band among the four observations, and the fluxes of all the background components are compatible, considering the statistical uncertainties.

The physical background, on the other hand, is not intrinsic to the instrument and can change depending on the observation and the pointing. The model we used comprised several components:

1. Unabsorbed thermal component representing the local hot bubble (LHB) with fixed kT = 0.097 keV and Z = Z₃ (Ponti et al. 2023) and free normalisation.

2. Absorbed thermal component representing the Galactic halo emission, with fixed Z = 0.1 Z₃ and free temperature and normalisation.

3. Absorbed thermal component representing the North Polar Spur (NPS), a region of enhanced soft X-ray emission that is part of the radio Loop I (Kataoka et al. 2018). The projected position of J1740 is on the edge of the NPS, as can be seen from Fig. 1, so its emission cannot be neglected for the background modelling. We fitted it with all the parameters free to vary because fixing a temperature and an abundance would be too restrictive, since the NPS is probably a multi-component plasma.

4. An absorbed double broken power law representing the cosmic X-ray background (CXB), following the method of Ponti et al. (2023), in turn based on the model of Gilli et al. (2007). We let the normalisation free to vary, but it was later fixed to avoid a degeneracy with the power law component later added to model the source (see Sect. 2.1.2).

For the absorption, we first employed the tbabs model described in Wilms et al. (2000) and fixed column density n_x to the value at the coordinates of the pulsar, n_x = 8 · 10²⁰ cm⁻², derived with the Xspec task nh. Then, to properly take into account the possible fluctuations of the column density inside the studied regions, we substituted tbabs with the disnht model developed by Locatelli et al. (2022). This model assumes the column density to be distributed as a log-normal function and can be considered as a generalisation of tbabs, better suited for extended and faint sources. In the end, the expression for the model used in the fit of the physical background is apee_LHB + disnht*(apee_Halo + apec_NPS + bkn2pow_CXB).

2.1.2 Source modelling

To search for diffuse non-thermal X-ray emission in the two regions, we added an absorbed power law (disnht*powerlaw) component to the model, leaving both normalisation and index free to vary. The parameters of disnht were set to be equal to those of the physical background model. To avoid a degeneracy with the CXB component, we fixed its normalisation to the expected value for the area considered in the source region. To do so, we used the best-fit value of the model by Gilli et al. (2007) and derived the normalisation for the areas of both Region 1 and Region 2. The results are n_CXB = 8.28 · 10^-5 ph s cm² and n_CXB = 2.57 · 10⁻⁵ ph s cm² for Region 1 and 2, respectively. We then performed the fit in the range 0.5-10 keV.

2.2 Analysis of the tail-like PWN

As for the diffuse emission analysis, we employed MOS2 data only and reduced them using the task emproe to obtain the event files. We lowered the contribution of the SPF by selecting the good time intervals after applying a maximum limit on the count rates to the lightcurve above 10 keV. For the source spectral extraction, we selected the elliptical region shown in Fig. 2 (solid ellipse), excluding the three point sources (small black circles) displayed in the same figure. The background was chosen as a circular region located in the same CCD but sufficiently far from the pulsar and its PWN to exclude possible contamination from both, and did not include any significant point-like source present in the FoV (see dashed circle in Fig. 2). Even though the background region is included in the extension of Region 1, the contamination by X-ray emission possibly associated with the LHAASO source can be considered negligible, since no excess of non-thermal X-rays was found by performing the analysis of the diffuse emission (see Sect. 3.1). The task baekseale was used on both the source and background spectra to properly take into account the different areas of the two regions and the possible chip gaps or bad pixels. The RMF and ARF were produced using the rmfgen and arfgen tasks, respectively. Also, in this case, we stacked the four observations by using mathpha on the spectra of both the source and background, and we produced an average RMF and ARF with addrmf and addarf. The final source spectrum was grouped using grppha with a minimum of 175 counts per bin to ensure at least 25 counts per bin in the background-subtracted spectrum.

Finally, we used χspee to model the tail emission as an absorbed power law tbabs*powerlaw, where we fixed the column density to the value at the pulsar’s coordinates n_χ = 8 · 10²⁰ cm⁻² (see Sect. 2.1.1). In this case, we employed the simpler tbabs model, since the source is not as extended as the diffuse emission. The fit was performed between 0.5 keV and 8 keV, after which the instrumental background of the CCDs becomes dominant over the physical signal, ignoring the 1.3-1.8 keV band. Due to the choice of the background region, located on one of the edges of the central CCD, some residuals around the Al Kα (1.49 keV)-Si Kα (1.74 keV) complex could still be observed after the background subtraction. This occurs because of the peculiar distribution of the lines’ strength in the detector plane, as shown in Kuntz & Snowden (2008). We produced the same plots with the most up-to-date merged filter wheel closed (FWC)⁴ data. The same line strength distribution of both the Al Kα and Si Kα reported in Kuntz & Snowden (2008) is observed.

The tail region previously defined was subsequently divided into five smaller and concentric elliptical regions, as shown in Fig. 2 (dashed ellipses). The procedure described above was repeated on each one of them to investigate the evolution of the spectral model as a function of the distance from the pulsar. We fitted the smaller tails in the range 0.5-6 keV, again excluding the 1.3-1.8 keV band, to remove the last energy bins where the minimum of 25 counts in the background-subtracted spectrum was not satisfied. The restriction does not affect the final result, since most of the tail emission is mostly concentrated where the instrumental background is not dominating.

Fig. 2 XMM-Newton image of the tail region in the energy range 0.58 keV. The solid ellipse represents the region where the overall tail emission has been extracted, while the dashed ellipses mark the borders of the sub-tail regions. The dashed circle is the background extraction region. We also display the Galactic coordinates grid.

3 Results of the X-ray analysis

3.1 Diffuse emission

When producing the model including only the background components, we first compared the best-fit results of the CXB normalisation to that expected from Gilli et al. (2007), n_CXB = 8.28 · 10⁻⁵ ph s cm² and n_CXB = 2.57 · 10⁻⁵ ph s cm² for Regions 1 and 2, respectively (see Sect. 2.1.2). In the case of Region 1, the best-fit result is n_CXB = (8.2 ± 0.4) · 10⁻⁵ ph s cm², with errors at 90% confidence level. The result is fully compatible with the theoretical prediction, suggesting that no additional components are probably needed in the fit. Instead, for Region 2, the best-fit normalisation is n_CXB = (3.9 ± 0.4) · 10⁻⁵ ph s cm², with 90% confidence level uncertainties. The value is higher than the expectations, although we note that it is reported without an associated error. This hints that some residual non-thermal emission could likely be present in the region.

As explained in Sect. 2.1.2, we subsequently fitted the diffuse emission with an absorbed power law with free index and normalisation. In the case of Region 1, the non-thermal component was not detected, and its addition did not improve the χ² value. This is compatible with what we observed when studying the normalisation of the CXB. In the case of Region 2, instead, the best-fit parameters are $\mathrm{\Gamma }={2.1}_{-1.1}^{+1.9}$ and $N_0={2.8}_{-1.9}^{+3.1}\cdot {10}^{-5}$ phkeV⁻¹ s⁻¹ cm⁻², with statistical uncertainties at 3σ level. The improvement in the χ² after adding the absorbed power law is from χ² = 240.6 (178 d.o.f.) to χ² = 235.8 (176 d.o.f.). Even though it could be considered a detection at 3σ, we treat the result as an upper limit for several reasons. First, we see that the improvement in the χ² after adding the source component is not very significant, hinting at a weak source detection. More importantly, some assumptions were made when building the background model, and fluctuations in the choice of the parameters for the adopted functions could easily lead to artefacts that can be observed as true detections. As shown in Fig. 1, there does not seem to be any strong flux in Region 1, while the emission from Region 2 might be affected by the presence of the NPS or fluctuations in the background.

Since we do not have a detection of the diffuse X-ray emission, the spectral index is not determined, so we derived the total upper limit on the X-ray emission by considering different values of the photon index in the spectrum. We obtained the 3σ upper limit on the normalisation at 1 keV for values of the photon index between Γ = 3.0 and Γ = 0.5. We represented all the derived power laws in the same plot, and we computed an overall upper limit in 0.5-10 keV. We computed it for both Region 1 and 2, and we subsequently used them together with the spectrum reported in the LHAASO catalogue to perform the SED fitting, see Sect. 3.3. For the SED, we assumed that the emission is uniform as a function of radius. Both upper limits for Region 1 (radius r ₁ = 0.11°) and Region 2 (radius r₂ = 0.06°) were rescaled to the area of the LHAASO region, assuming r_LHAASO = 0.147°, as provided by LHAASO Collaboration (2025).

3.2 Tail-like PWN analysis

We obtained a good spectral fit (χ² = 70.57 for 75 d.o.f.) of the tail emission for a spectral index of Γ = 1.76 ± 0.06 and normalisation at 1 keV of N₀ = (3.49 ± 0.14) · 10⁻⁵ phkeV⁻¹ s⁻¹ cm⁻² or N₀ = (33.1 ± 1.3) phkeV⁻¹ s⁻¹ cm⁻²sr⁻¹ when expressing it in surface brightness after setting the proper value of the AREASCAL keyword. The 0.5-8 keV unabsorbed flux is F_x = (1.86 ± 0.12) · 10⁻¹³ ergs⁻¹ cm⁻², corresponding to an unabsorbed luminosity of L_x = (3.4 ± 0.2) · 10³¹ erg s⁻¹ assuming a distance of d = 1.23 kpc (Yao et al. 2017). From the luminosity, we derived the X-ray efficiency as η = L_x/Ė = 1.5 · 10⁻⁴. The result for the flux is compatible within the statistical uncertainties with that of Benbow et al. (2021), F_X = (1.93 ± 0.06) · 10⁻¹³ erg s⁻¹ cm⁻², even though they reported the absorbed flux of the PWN in the energy range 0.3-10 keV. On the other hand, there is a discrepancy between our result for the luminosity and theirs, L_X = (5.20 ± 0.16) · 10³¹ erg s⁻¹, most likely due to different source and background extraction regions and energy bands employed in the two analyses and a different assumed value for the pulsar distance. Finally, we extracted the SED flux points for this dataset directly from Xspec and used them together with the LHAASO power law spectrum in the SED fitting, also comparing them to the upper limit by Benbow et al. (2021) (see Sect. 3.3).

The same fit procedure has been applied to the single subtails, and the results are summarised in Table 1, along with the best-fit values of the model for the overall tail emission. The reported normalisation values are computed by including the region area in the AREASCAL parameter to obtain a surface brightness value at 1 keV. The spectra of the sub-tails are also depicted in Fig. 3. No evidence for an evolution of the spectral index is found, i.e. no synchrotron cooling is detected across the length of the tail, as it can also be observed in Fig. 4 (left). The trend of the surface brightness at 1 keV, illustrated on the right plot of Fig. 4, shows that the innermost part of the tail, closer to the pulsar, is brighter than the outer part. The horizontal bars associated with the distances correspond to the inner and outer boundaries of the ellipses.

3.3 SED modelling

We performed an SED modelling over the X-ray and teraelectronvolt range using the package naima (Zabalza 2015), which allows the user to study the radiative losses of a relativistic population of particles, considering a purely leptonic scenario. We did not treat a hadronic scenario because it is beyond the aim of this paper. Moreover, LHAASO Collaboration (2025) already discussed the low likelihood of this kind of model for 1LHAASO J1740+0948u due to the absence, in the vicinity of the source, of molecular clouds that could give rise to hadronic emission.

The goal of this SED modelling was to try to put constraints on the origin of the emission observed by LHAASO by exploring different scenarios.

Fig. 3 Spectra of the single tails in the energy range 0.5-6 keV, excluding the 1.3-1.8 keV band due to the contamination of instrumental lines. Tail 1 is the innermost with respect to the pulsar, Tail 5 is the outermost.

3.3.1 Scenario 1: potential VHE detection of the X-ray tail

We first investigated the hypothesis that 1LHAASO J1740+0948u is the counterpart of the X-ray tail. To this aim, we performed the SED modelling using the previously produced flux points in the range 0.5-8 keV and the KM2A spectrum of LHAASO Collaboration (2025) as well as their extrapolation for the WCDA instrument. We assumed the electron distribution at injection to be an exponential cut-off power law (ECPL) $$f(E)=A(\frac{E}{E_0}{)}^{-p}\exp [-(\frac{E}{E_{cut}}{)}^{\beta }],$$(1)

where E₀ is the reference energy, fixed to E₀ = 1 TeV, p is the spectral index, E_cut is the cut-off energy, β the cut-off parameter, fixed to β = 2 as in Zirakashvili & Aharonian (2007), and N₀ is the normalisation, which is related to the total energy content W_e. The minimum and maximum energy of the particle distribution were set to E_min = 1 GeV and E_max = 960 TeV, respectively. The maximum energy was obtained from Eq. (3) of de Oña Wilhelmi et al. (2022), setting η_e = η_B = 1, where η_e represents the ratio between the electric field strength and the magnetic one, and η_B represents the fraction of pulsar wind energy flux that goes into magnetic energy density (de Oña Wilhelmi et al. 2022). The minimum, instead, was chosen taking into account that no points below the X-ray band are included in the fit. We considered synchrotron and ICS as the emission mechanisms for the particle population. For the latter we included as photon field seeds the CMB, far-infrared (FIR) radiation with T_FIR = 33 K and u_FIR = 0.2 eV/cm³, near-infrared (NIR) radiation with T_NIR = 350 K and u_NIR = 0.02 eV/cm³, and the starlight contribution with T_star = 4800 K and u_star = 0.3 eV/cm³, where all the values were obtained from Popescu et al. (2017) at the pulsar location.

We first performed the fit, fixing the spectral index to the expectation value from the X-ray spectrum p = 2Γ_x - 1 (Kargaltsev & Pavlov 2010) that corresponds to p = 2.52 for Γ_x = 1.76. The derived best fit model, with parameters E_cut = 329-27 TeV, B = 6.8 ± 1.9 μG and W_e = 1.3 · 10⁴⁵ erg, only matched the X-ray points and was not able to reproduce the ICS component observed by LHAASO. The magnetic field obtained through this fit, however, is compatible with previous literature results (Kargaltsev et al. 2008a; Benbow et al. 2021). To better fit both X-ray and gamma-ray data, we performed several tests changing the spectral index of the electron distribution. We compared the results of each test by using the Bayesian information criterion, BIC (Schwarz 1978), defined as BIC = k ln(n) - 2 ln(L), where k is the number of free parameters, n is the number of points in the sample and L represent the likelihood of the model. The BIC is computed by naima during the fitting procedure, along with the maximum logarithmic likelihood value. Lower values of the BIC usually indicate a statistical preference for the associated model with respect to the other. The model that returned the lowest value of BIC and also had a match with the LHAASO spectrum was the one with p = 1.6, with best-fit results $E_{cut}={304}_{-21}^{+32}$ TeV, B = ${1.31}_{-0.07}^{+0.10}$ μG and W_e = 1.34 · 10⁴⁴ erg. The best-fit SED and its residuals are depicted in Fig. 5. The electron energy distribution from the best-fit model shows a harder spectrum (p = 1.6) compared to the expected value p = 2Γ_x - 1 = 2.52, and the expected flux at tens of tera-electronvolts is below the observed spectrum by a factor of two or three. Moreover, the best-fit magnetic field value is significantly smaller than that obtained for p = 2.52 and previous literature results. Even though with naima we are estimating only the perpendicular component of the magnetic field, it has been proven that the total magnetic field is dominated by that component (Liu et al. 2019b), so the real value of the magnetic field will still be lower than the expected ~6 μG. This hints that the ICS component generated by the electrons responsible for the X-ray emission is insufficient to explain the tera-electronvolt flux detected by LHAASO. Therefore, it suggests that LHAASO probably did not detect the counterpart of the X-ray tail and that the two sources are related to different populations of electrons, most likely a young electron population generating the X-ray emission and an older population emitting at tera-electronvolt energies.

3.3.2 Scenario 2: diffuse X-ray emission from the LHAASO source

The second hypothesis we investigated was whether 1LHAASO J1740+0948u could be associated with the diffuse X-ray emission in the surroundings of J1740. To evaluate this scenario, we assumed that both the LHAASO source and the diffuse X-ray emission originated from the same relativistic population of electrons and adopted the same ECPL distribution with fixed E₀ = 1 TeV, β = 2, E_min = 1 GeV and E_max = 960 TeV as for the previous case. Since there is no constraint on the spectral index from the X-ray results, we first froze it and tested values between p = 1.2 and p = 2.4, then we let it free to vary. We also included in the fit the ICS emission from the X-ray tail, assuming the whole emission is due to the synchrotron-emitting particles, i.e. p = 2.52, to take into account its possible contribution inside the LHAASO energy range. To model the synchrotron emission, we tested different values from 0.5 to 1.6 μG to determine the maximum value of the magnetic field of the diffuse emission before overshooting the X-ray limit. In Table 2 we report the best-fit values for E_cut (in both electron and photon energy) and We obtained with naima when fixing p of the ECPL in the range 1.2-2.4, as well as the maximum logarithmic likelihood and the BIC value derived during the fit. We depicted in Figure 6 the models for p = 1.2 (left) and p = 2.4 (right). The derived values of both the BIC and maximum logarithmic likelihood slightly favour the models with the highest values of p. This is also confirmed by the best-fit result when using a free spectral index, for which we obtained p = 2.2 ± 0.2. Leaving p free, however, did not improve the quality of the fit. In all the tested models, the values of the magnetic field required not to overshoot the X-ray upper limit are ~0.6 μG for Region 1 and ~ 1.1-1.2 μG for Region 2, respectively. The different values in the two regions are likely due to the excess of non-thermal emission detected in Region 2 that causes the upper limits to be higher. Another possibility is related to the fact that Region 2 is located on the edge of the XMM-Newton FoV, where the sensitivity is lower than at the centre, so the estimation could be affected by additional uncertainties. Such low values of the magnetic field are compatible with the hypothesis of them originating from an old electron population, too faint to be X-ray-emitting but only visible in the tera-electronvolt domain.

Fig. 4 Evolution of the spectral parameters as a function of the angular distance from the pulsar position. The uncertainties on the index and surface brightness values are shown at 1σ level. Left panel : the trend of the spectral index. Right panel : the trend of the normalisation of the power law at 1 keV considering the AREASCAL parameter.

Fig. 5 SED of the tail-like PWN of J1740, including the X-ray points derived in this work and the LHAASO-KM2A spectrum along with the LHAASO-WCDA extrapolation of LHAASO Collaboration (2025). The electron distribution is modelled as an ECPL with fixed p = 1.6, β = 2, E0 = 1 TeV and free Ecut and normalisation. We included ICS, considering all the photon fields obtained from Popescu et al. (2017), and Synchrotron, where the magnetic field value was left free to vary. We also report the VERITAS upper limit of Benbow et al. (2021) obtained for an extraction region of 0.1°. The model residuals, obtained as (data - model)∕error, are shown in the bottom panel.

3.4 Cooling times and distances

To further support the hypothesis that the two populations, i.e. the one of the X-ray tail and that of the LHAASO source, are different, we estimated their cooling times. We followed the prescriptions of Khangulyan et al. (2014) to compute the cooling times, and we used Eq. (1) of Van Etten & Romani (2011) to estimate the particle energy from the kilo-electronvoltemitting photons and Eq. (1) of LHAASO Collaboration (2021) for the tera-electronvolt photons. We obtained a total cooling time for the particles observed in the energy range 0.5-8 keV of τ_tot ~ 1.2-4.3 kyr, assuming the best-fit value of the PWN magnetic field obtained with naima. On the other hand, we derived a total cooling time in the range τ_tot ~ 10-26 kyr for the tera-electronvolt-emitting particles, assuming a magnetic field of 1 μG. If we focus on the synchrotron cooling time only, we obtained τ_syn ~ 1.2-4.7 kyr for the X-ray population and τ_syn ~ 13-100 kyr for the VHE population.

In Fig. 4, we see no spectral evolution with the distance from the pulsar, meaning that no cooling is detected along the entire visible length of the tail. To further support this finding, we computed the cooling distances of the electrons of the tail assuming different regimes for the particle transport, ballistic transport (v = c), advection (v = c/√3), and diffusion. For the latter, we used the following relation $$D(E)=D_0(\frac{E}{10\text{ GeV}}{)}^{\frac{1}{3}},$$(2)

where D₀ = (3 - 8) · 10²⁸ cm² s⁻¹, taken from Giacinti et al. (2018), and we chose D₀ = 5 · 10²⁸ cm² s⁻¹ as the reference value at 10 GeV. We found distances in the range 0.4-1.3 kpc, 200750 pc, and 88-140 pc for the ballistic transport, advection and diffusion, respectively. We also tested a slow-diffusion scenario, assuming D₀ = 1.2 · 10²⁶ cm² s^-1 at 1 GeV, which is the average value of the best-fit results reported in Albert et al. (2024) for the Geminga and Monogem halos. In this case, we obtained distances in the range of 6.1-9.6 pc. All the derived distances are larger than the physical offset observed for the LHAASO source, corresponding to d = 4.7 pc for a distance of J1740 of 1.23 kpc (d = 5.4 pc for 1.4 kpc used in the LHAASO Collaboration (2025) paper), suggesting particles could be able to travel to the LHAASO location without invoking any re-acceleration process.

Fig. 6 SED model of the LHAASO spectrum along with the X-ray ULs for Region 1 (solid UL) and Region 2 (dashed UL) using an ECPL distribution for the relativistic electrons generating the emission with fixed spectral index p. The contribution to the ICS from the X-ray tail is shown as the dash-dotted line. Values of the magnetic field between 0.5 μG and 1.6 μG are tested. Left panel : plot for p = 1.2. Right panel : plot for p = 2.4. In both panels, the upper limits were computed by rescaling the results to the area of the LHAASO region, assuming a radius of 0.147° as LHAASO Collaboration (2025).

4 Discussion and conclusions

In this work, we conducted a deep analysis of the X-ray emission around PSR J1740+1000 (J1740), a faint radio and gammaray pulsar located outside the Galactic plane, using ~500 ks of XMM-Newton observations. We re-analysed the long taillike BSPWN and investigated for the first time the evolution of the spectral index as a function of the distance from J1740. We also performed the first study of the diffuse emission in the surroundings of J1740, using a custom model for the background emission. Our goal was to constrain a possible X-ray association with the recently detected UHE source 1LHAASO J1740+0948u, located at a 0.22° offset from the pulsar position.

The analysis of the cometary tail proved the best-fit model in the energy range 0.5-8 keV to be an absorbed power law with photon index Γ = 1.76 ± 0.06, compatible with the result obtained by Benbow et al. (2021) on the same dataset but with an independent analysis. The efficiency of the PWN, derived from the X-ray luminosity and assuming a dispersion measure distance of 1.23 kpc, is η_X = 1.5 · 10⁻⁴, a typical value observed for middle-aged pulsars (Kargaltsev & Pavlov 2008). For the first time, we explored the evolution of the spectral index with distance and found no evidence for spectral cooling along the length of the tail. A decrease in the absolute value of the spectral index is found in the outermost region of the tail, at an angular distance of ~0.1° from the pulsar position, but the poor statistic implies large error bars on the result and makes the result still compatible with that of the other regions. It is still not well understood why some BSPWNe, like the Lighthouse PWN (Pavan et al. 2016) or the Mouse PWN (Klingler et al. 2018), show a spectral evolution, while other systems, such as the one of J1740 studied in this work or that of PSR B0355+54 (Klingler et al. 2016), do not show this feature. A possible explanation could be related to the different environments where the PWNe are located or different interactions between the particles and the surrounding medium, possibly connected to the age of the population.

For the study of the diffuse emission, instead, we focused on two regions of the XMM-Newton FoV, due to the limited coverage of the extension of 1LHAASO J1740+0948u provided by the instrument. The first (Region 1) is a r = 0.11° circular region centred on the pulsar and excluding both the tail and the pulsar to avoid their contamination. The second (Region 2) is a r ~ 0.06° circle comprising the majority of the 68% extension of the VHE source. We fitted custom models for the instrumental and the physical backgrounds, where the latter comprised several components: the cosmic X-ray background, the Galactic halo, the local hot bubble and the North Polar Spur, a region of enhanced X-ray emission whose position coincides with J1740's location. To reproduce the non-thermal diffuse component we aimed to study, we added an absorbed power law corresponding to the source term. No detection was found for Region 1, while a weak component was detected for Region 2. This component, however, was only marginally detected, and deeper observations or coverage of a more extended area would be required to firmly exclude that such a feature is not due to a background fluctuation. Therefore, we treated it as a non-detection. We then derived the 3σ upper limits on the emission in the range 0.5-10 keV for both regions by combining the results obtained for fixed values of the photon index between Γ = 0.5 and Γ = 3.0. To constrain the nature of 1LHAASO J1740+0948u, we investigated two main hypotheses: (1) LHAASO detected the X-ray tail in the VHE band, and (2) the VHE source may be related to the diffuse non-thermal X-ray emission.

To explore hypothesis (1), we performed an SED fitting of the tail emission together with the spectrum displayed in LHAASO Collaboration (2025), assuming the particles to be distributed as an exponential cut-off power law. The approach is similar to that presented in Benbow et al. (2021), published before the discovery of 1LHAASO J1740+0948u, but we were able to further improve it. We performed a multi-wavelength MCMC fitting, including both a spectrum in the X-rays and in the VHEs, allowing for a better constraint of the parameters of the electron distribution when compared to a fit of the X-ray data alone. The model obtained when fixing the spectral index to the expected value from the X-ray emission, i.e. p = 2.52, is not able to fit the VHE part of the SED, but constrains the tail magnetic field to B = 6.8 ± 1.9 μG, compatible with the predictions by models of BSPWNe (Olmi & Bucciantini 2019) and with previous literature results on J1740 (Benbow et al. 2021 ; Kargaltsev et al. 2008a). The projected location of J1740 in the sky is coincident with the NPS (see Fig. 1), and what we obtained for the tail magnetic field is of the same order as the average value measured for the NPS (Iwashita et al. 2023). However, the distance to the pulsar as measured by Yao et al. (2017) is d ~ 1.23 kpc, whereas the root of NPS is measured to have a distance of several kilo-parsecs (Zhang et al. 2024). Therefore, the magnetic field measured from SED fitting in this work corresponds to the vicinity of the PWN, not the NPS. In order to have a better fit of the LHAASO spectrum, we conducted several tests and found that the model with p = 1.6 returned the best-fit results, both in the match with the VHE points and in the results of the BIC and likelihood. This index is considerably harder than the expected p = 2.52 and underestimates both the magnetic field of the tail and the VHE flux at ~10 TeV by a factor of two or three. These results show that LHAASO likely did not detect the counterpart of the X-ray tail and that the two populations, the X-ray bright one and that observed in the VHEs, are distinct: the VHE emission traces older electrons that are not emitting in the X-rays any more due to cooling, while the tail is populated by young and fresh electrons injected from the pulsar.

We then explored hypothesis (2) by carrying out an SED fitting using the spectrum of 1LHAASO J1740+0948u LHAASO Collaboration (2025) and the X-ray upper limits on the diffuse emission obtained for Regions 1 and 2. With this approach, it was possible to constrain the magnetic field to B ≤ 0.6 μG for Region 1 and to B ≤ 1.1-1.2 μG for Region 2, where the difference between the second and the first value could be due to the non-thermal excess marginally detected for Region 2, to physical fluctuations of the magnetic field in the medium or to the result being affected by the worse performance of XMM-Newton when observing off-axis. Similar methods for the SED modelling were adopted in LHAASO Collaboration (2025) and Xie et al. (2025). Both works attempted to put constraints on the electron population generating the LHAASO emission and its possible X-ray counterparts. In our study, we included in the picture the very first upper limits on the diffuse emission as obtained by XMM-Newton and we constrained the magnetic field of the diffuse emission based on the real observations. In addition, we relaxed some of the assumptions that were used in both previous works for the SED modelling parameters.

Considering these results, we argue that the X-ray tail is produced via synchrotron radiation of freshly injected electrons in a magnetic field of ~6.8 μG, while the LHAASO emission traces older electrons. This population is not emitting in the X-rays anymore due to cooling and the decrease of the magnetic field to ≤0.6-1.2 μG, but it is still shining via ICS of the ambient photon fields. One possible interpretation of the nature of this system is that the LHAASO emission represents the relic PWN of J1740. For these sources, it is common to observe offsets between the centroid of the VHE source and the pulsar position or the X-ray nebula, if detected. This is compatible with the 0.22° displacement observed between 1LHAASO J1740+0948u and the radio and X-ray coordinates of J1740. There are several examples of relic PWNe systems of different ages showing similar offsets: HESS J1825-137, showing an extended asymmetric morphology 0.2° offset from PSR J1826-1334 and its X-ray nebula (Aharonian et al. 2006; H.E.S.S. Collaboration 2019); HESS J1616-508, without a detected counterpart in the X-rays and located at ~10′ from PSR J1617-5055 and its faint PWN (Kargaltsev et al. 2008b); HESS J1356-645, associated with PSR J1357-6429, which hosts an extended nebula and is located at a distance of ~7′ from the VHE source (H.E.S.S. Collaboration 2011); or HESS J1303-631, which centroid is at ~3.1′ from PSR J1301-6305 and was one of the first discovered ‘dark’ sources without an X-ray counterpart (H.E.S.S. Collaboration 2012). There are two main explanations for the observed offsets between X-rays and VHEs. They could be either caused by asymmetries in the medium and, consequently, in the SNR reverse shock that interacts with the expanding PWN (Blondin et al. 2001), or due to the supersonic motion of the pulsar related to the kick velocity that it has acquired during its birth (Fiori et al. 2022). Considering J1740 is classified as a BSPWN, the second scenario seems to be the most likely, even though an asymmetric expansion of the parent SNR, not detected either in radio or X-rays so far, in principle could have occurred in the past. Considering that the ~0.22° offset between the LHAASO centroid and the position of J1740 corresponds to only ~5 pc, for both d = 1.23 kpc used in this paper and d = 1.4 kpc adopted by LHAASO Collaboration (2025), this kind of interpretation does not need to invoke any re-acceleration process to explain the connection between the X-ray upper limit and the VHE emission.

The other possible scenario is that 1LHAASO J1740+0948u represents the pulsar halo of J1740. This could be a realistic hypothesis considering the age (τ_c = 114 kyr) and its classification as a BSPWN that has likely already escaped the parent SNR and is interacting directly with the ISM. It is interesting to notice that the High Altitude Water Cherenkov (HAWC) experiment detected a weak (TS=28) VHE source near J1740, dubbed 3HWC J1739+099, but it is a point-like source and its centre is off-site compared to the LHAASO source (Albert et al. 2020). The source, which has been classified as a candidate pulsar halo, has not been included in our work due to the uncertainty in the connection with 1LHAASO J1740+0948u and the lack of spectral points of the source. To further investigate the interpretation of 1LHAASO J1740+0948u as a pulsar halo, we followed a phenomenological approach. A criterion to understand whether a VHE source can be classified as a pulsar halo is to compare its energy density to that of the ISM: if the former is lower than the latter, it is evidence for outflowing electrons towards a region not influenced by the pulsar (López-Coto et al. 2022). We computed the energy density of 1LHAASO J1740+0948u by following the method explained in Giacinti et al. (2020), using the W_e value computed by naima instead of performing the integral of the spectrum. We assumed r = 0.147°, which is the 95% confidence level upper limit on the source size under the pointlike source assumption (LHAASO Collaboration 2025), as the radius of the region. The results for the spectral indices between p = 1.2 and p = 2.4 are in the range 0.03-3.6 eV/cm³. Excluding the extreme value obtained for p = 2.4, which is likely a too soft index for an injection spectrum, the energy densities are in the range 0.03-0.67 eV/cm³. When leaving the spectral index free to vary in the fit, we obtained p = 2.2 ± 0.4 with a corresponding E_e = 3.45 · 10⁴⁵ erg and an energy density of ε = 0.57 eV/cm³. The conventional value assumed for the ISM energy density is 0.1 eV/cm³, but sources with energy densities lower than ~1 eV/cm³ may already display halo-like features. When studying the cooling of the tail electrons, we found they can reach very large distances from the pulsar. Assuming slow diffusion, instead, we derived a range of ~6-9 pc, which is larger than the offset than the physical offset observed between the pulsar and the LHAASO centroid (d ~ 5 pc). Even in a slow diffusion context, electrons may be able to reach the LHAASO location and generate a pulsar halo. Moreover, projection effects could explain both the offset and the small size, as already discussed in LHAASO Collaboration (2025). Therefore, we argue 1LHAASO J1740+0948u may likely be the pulsar halo of J1740 and, if this hypothesis is confirmed, we constrained the magnetic field around the pulsar to be as low as B ≤ 0.6-1.2 μG, lower than the typical value of the ISM magnetic field and compatible with other results obtained for the search of X-ray pulsar halos (Khokhriakova et al. 2024; Manconi et al. 2024; Adams et al. 2025). Future on-site X-ray observations should help unravel the nature of this UHE source.

Due to the lack of multi-wavelength coverage on the source and the absence of dedicated X-ray observations covering the full extension of 1LHAASO J1740+0948u, it was beyond the aims of this paper to develop a detailed treatment of diffusion to further support the pulsar halo scenario. A more complicated modelling, combining a spectral and spatial analysis of the multi-wavelength emission, is likely needed to investigate the source more in detail and have a better understanding of the relation between the X-ray bow-shock tail and the LHAASO emission. Future observations by particle detectors, such as LHAASO and HAWC, or by the upcoming Cherenkov Telescope Array Observatory, could help us understand more about this dark VHE source and whether the two detections, i.e. 1LHAASO J1740+0948u and 3HWC J1739+099, represent the same source. Future X-ray facilities, such as NewAthena (Cruise et al. 2025), could allow us to distinguish better the diffuse emission around J1740 from the background and improve the upper limits obtained in this work.

Acknowledgements

This research is based on observations performed by the satellite XMM-Newton, a European Space Agency science mission with instruments and contributions directly funded by ESA Member States and NASA. This work has also made use of data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. GP and HZ acknowledge financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program HotMilk (grant agreement No. 865637). GP also acknowledges support from Bando per il Finanziamento della Ricerca Fondamentale 2022 dell’Istituto Nazionale di Astrofisica (INAF): GO Large program and from the Framework per l’Attrazione e il Rafforzamento delle Eccellenze (FARE) per la ricerca in Italia (R20L5S39T9).

References

All Tables

All Figures

Fig. 1 (a) ROSAT all-sky map showing the large-scale diffuse emission of the North Polar Spur and the location of the XMM-Newton pointings, i.e. the small green circle. The red, green and blue colours represent the 0.11-0.28 keV, 0.47-1.21 keV and 0.76-2.04 keV energy bands, respectively. (b) Image of the XMM-Newton FoV, showing both Region 1 and Region 2 (solid circles), and 1LHAASO J1740+0948u (dashed circles). The plotted radii of the source are taken from the 1 σ extension of the Gaussian (RG,LHAASO) and the diffuse (RD,LHAASO) models for photon energy above 25 TeV, shown in LHAASO Collaboration (2025). The 0.5-1 keV band is shown in red, the green emission corresponds to the range 1-2 keV, and the range 2-4.5 keV is represented in blue. Note that the observed point sources were later removed during the diffuse emission analysis.

In the text

	Fig. 2 XMM-Newton image of the tail region in the energy range 0.58 keV. The solid ellipse represents the region where the overall tail emission has been extracted, while the dashed ellipses mark the borders of the sub-tail regions. The dashed circle is the background extraction region. We also display the Galactic coordinates grid.
In the text

	Fig. 3 Spectra of the single tails in the energy range 0.5-6 keV, excluding the 1.3-1.8 keV band due to the contamination of instrumental lines. Tail 1 is the innermost with respect to the pulsar, Tail 5 is the outermost.
In the text

	Fig. 4 Evolution of the spectral parameters as a function of the angular distance from the pulsar position. The uncertainties on the index and surface brightness values are shown at 1σ level. Left panel : the trend of the spectral index. Right panel : the trend of the normalisation of the power law at 1 keV considering the AREASCAL parameter.
In the text

Fig. 5 SED of the tail-like PWN of J1740, including the X-ray points derived in this work and the LHAASO-KM2A spectrum along with the LHAASO-WCDA extrapolation of LHAASO Collaboration (2025). The electron distribution is modelled as an ECPL with fixed p = 1.6, β = 2, E0 = 1 TeV and free Ecut and normalisation. We included ICS, considering all the photon fields obtained from Popescu et al. (2017), and Synchrotron, where the magnetic field value was left free to vary. We also report the VERITAS upper limit of Benbow et al. (2021) obtained for an extraction region of 0.1°. The model residuals, obtained as (data - model)∕error, are shown in the bottom panel.

In the text

Fig. 6 SED model of the LHAASO spectrum along with the X-ray ULs for Region 1 (solid UL) and Region 2 (dashed UL) using an ECPL distribution for the relativistic electrons generating the emission with fixed spectral index p. The contribution to the ICS from the X-ray tail is shown as the dash-dotted line. Values of the magnetic field between 0.5 μG and 1.6 μG are tested. Left panel : plot for p = 1.2. Right panel : plot for p = 2.4. In both panels, the upper limits were computed by rescaling the results to the area of the LHAASO region, assuming a radius of 0.147° as LHAASO Collaboration (2025).

In the text