An XMM-Newton EPIC X-ray view of the Symbiotic Star R~Aquarii

We present the analysis of archival XMM-Newton European Photon Imaging Camera (EPIC) X-ray observations of the symbiotic star R Aquarii. We used the Extended Source Analysis Software (ESAS) package to disclose diffuse soft X-ray emission extending up to 2.2 arcmin ($\approx$0.27 pc) from this binary system. The depth of these XMM-Newton EPIC observations reveal in unprecedented detail the spatial distribution of this diffuse emission, with a bipolar morphology spatially correlated with the optical nebula. The extended X-ray emission shares the same dominant soft X-ray-emitting temperature as the clumps in the jet-like feature resolved by Chandra in the vicinity of the binary system. The harder component in the jet might suggest that the gas cools down, however, the possible presence of non-thermal emission produced by the presence of a magnetic field collimating the mass ejection can not be discarded. We propose that the ongoing precessing jet creates bipolar cavities filled with X-ray-emitting hot gas that feeds the more extended X-ray bubble as they get disrupted. These EPIC observations demonstrate that the jet feedback mechanism produced by an accreting disk around an evolved, low-mass star can blow hot bubbles, similar to those produced by jets arising from the nuclei of active galaxies.

1. INTRODUCTION R Aquarii (R Aqr) is a symbiotic star (SS) consisting of a 385-day period pulsating Mira-like star and a hot white dwarf (WD) harboring an accretion disk. At its Gaia geometric distance of 385±60 pc, it is one of the closest and best studied SS, with detailed determinations of its orbital parameters (Hollis, Pedelty, & Lyon 1997;Gromadzki & Mikołajewska 2009;Bujarrabal et al. 2018), physical characteristics (Schmid et al. 2017) and dust content (see Mayer et al. 2013).
Its outstanding morphology has been also subject of exhaustive studies. R Aqr is surrounded by an exquisite filamentary and knotty hourglass nebula as shown in optical wavelengths (Fig. 1). The optical images show several bowl-shape cavities open toward the north and south directions (Solf & Ulrich 1985;Michalitsianos et al. 1988; Corresponding author: Jesús A. Toalá j.toala@irya.unam.mx Liimets et al. 2018). A close inspection of the VLT [O III] image presented by Liimets et al. (2018) reveals faint optical emission closing the southern lobe. In addition, images obtained with the 30 cm telescope at Terroux Observatory reveal a faint closed lobe in the northern region extending 2.8 from R Aqr. At the heart of this complex nebula, there is a long-studied precessing jet with an S-shape morphology (see, e.g., Wallerstein & Greenstein 1980;Paresce & Hack 1994;Melnikov, Stute, & Eislöffel 2018). R Aqr is a complex system where all components seem interconnected (mass loss of the primary, accretion onto the secondary, production of inner intricate morphologies, etc.). Hence, the dimming of the periodic giant Mira-type star and the occurrence of the jets have been linked to the periastron passage of the binary system (see for example Kafatos & Michalitsianos 1982).
The detection of high ionization features in the jets in the UV regime (Michalitsianos, Kafatos, & Hobbs 1980) prompted several early X-ray studies. X-rays from R Aqr were marginally detected with Einstein observations (Jura & Helfand 1984) and subsequently confirmed by EXOSAT and ROSAT observations (see Viotti et al. 1987 . 1998). The X-ray emission of R Aqr remained unresolved until the advent of Chandra. Kellogg et al. (2001) used Chandra observations to show that most of the X-ray emission is associated with the SS as well as with clumps of the S-shaped jet. The spectrum associated with the SS shows the contribution from the Fe Kα line at 6.4 keV, which corroborated the accretion phenomena with material falling onto the WD component. Kellogg et al. (2001) showed that the X-ray-emitting clumps in the S-shaped jet were shock-heated with plasma temperatures 10 6 K, similarly to what is found in extremely collimated systems such as HH objects (Rodríguez-Kamenetzky et al. 2019;Favata et al. 2002), proto-planetary nebulae (Sahai et al. 2003), nova-like systems (Toalá et al. 2020;Montez et al. 2021) and other SSs (Galloway & Sokoloski 2004;Stute & Sahai 2009;Stute et al. 2013). Subsequent X-ray studies performed by the same group demonstrated that the X-ray-emitting clumps detected by Chandra exhibit variability very likely due to the jet precession (e.g., Kellogg et al. 2007Kellogg et al. , 2009). Kellogg et al. (2007) showed that the spatial distribution of the X-ray-emitting gas changes within a few years with an apparent projected velocity of 580×(d/200 pc) km s −1 , where d is the distance to the object. The authors found that the clump located at the SW from the SS seemed to fade between 2000 and 2004 and suggested that it might be due to its adiabatic expansion and cooling, even though they predicted cooling times longer than 800 yr.
In this letter we present the analysis of archival XMM-Newton observations using reduction techniques for extended sources. We report the detection of extended X-ray emission beyond the S-shaped jets with a double-lobe morphology filling the space between the hourglass nebular structures. The observations and data preparation are presented in Section 2, while the results are presented in Section 3. A closing discussion addressing the origin of this extended emission and its spectral properties in light with contemporaneous archival Chandra observations is presented in Section 4.
2. OBSERVATIONS AND DATA PREPARATION R Aqr was observed by XMM-Newton on 2005 June 30 using the European Photon Imaging Camera (EPIC) (PI: E. Kellogg; Obs. ID. 0304050101). The three EPIC cameras, namely MOS1, MOS2, and pn, were used in the full frame mode with the medium optical blocking filter for exposure times of 72.2, 72.2, and 70.6 ks, respectively. The data were retrieved from the XMM-Newton Science Archive 1 and processed with the Science Analysis Software (SAS, version 17.0; Gabriel et al. 2004), using the most recent calibrations available by January 2022.
To take advantage of the good quality of the XMM-Newton observations of R Aqr, we analyzed the EPIC data in two different ways. First, we started our analysis making use of the Extended Source Analysis Software (ESAS) package, currently included as part of the SAS. The data were processed following the cookbook for analysis of extended objects and diffuse background 2 . The ESAS tasks apply restrictive event selection criteria, but their results leverage the presence of soft (E < 1.5 keV) extended emission (see Snowden et al. 2004Snowden et al. , 2008).
2 https://xmm-tools.cosmos.esa.int/external/sas/current/doc/esas/ After processing the data with the ESAS tasks, the final net exposure times of the MOS1, MOS2, and pn cameras are 33.8, 39.1, and 20.4 ks respectively. The processed files were used to create EPIC images in the 0.3-0.7 keV (soft), 0.7-1.2 keV (medium), and 1.2-5.0 keV (hard) energy bands. Individual pn, MOS1, and MOS2 images were extracted, corrected for exposure maps, and merged together. The soft, medium, and hard band images were adaptively smoothed using the ESAS task adapt requesting 30, 10, and 10 counts, respectively. The final images are presented in Figure  well as a color-composite X-ray image combining the three bands.
To study the physical properties of the extended X-ray emission detected with XMM-Newton, we need to extract and model its spectrum. Since the ESAS event selection criteria are too restrictive for this purpose, we reprocessed the EPIC data this time using the epproc and emproc tasks. To remove periods of high background levels in the data, we produced light curves in the 10-12 keV energy range and excised time periods with count rates higher than 0.2 and 0.45 counts s −1 for the MOS and pn cameras, respectively. After excising bad time intervals, the event files have net exposure times of 59.3, 60.0, and 35.8 ks for the MOS1, MOS2, and pn, respectively.
For comparison and discussion we also retrieved Chandra observations of R Aqr from the Chandra Data Archive 3 . Several sets of Chandra observations R Aqr are available, but, given its X-ray variability, we choose to retrieve those observed at the same epoch of the XMM-Newton EPIC data, i.e., those corresponding to the Obs. ID. 5438 (PI: E. Kellogg) obtained on 2005 October 9 with a total exposure time of 66.7 ks. The data were obtained with the Advanced CCD Imaging Spectrometer (ACIS)-S in the VFAINT mode and R Aqr was registered on the back-illuminated CCD S3. The Chandra data were processed with the Chandra Interactive 3 https://cda.harvard.edu/chaser/ Analysis of Observations (CIAO, version 4.7.4; Fruscione et al. 2006).

Extended X-ray emission
The XMM-Newton EPIC images produced with the ESAS tasks unambiguosly reveal the presence of extended X-ray emission associated with R Aqr. Figure 2 reveals the presence of X-ray emission up to 2.2 from the SS. In addition, Figure 2 shows that the bulk of the emission is emitted in the soft band. There is some extended emission detected in the medium band distributed in the inner ∼30 very likely associated with the X-ray-bright clumps resolved in the Chandra data (see, e.g., Kellogg et al. 2007). The hard X-ray emission is only associated with the SS. This energy-dependent spatial distribution of the hot gas can be appreciated in the color-composite picture presented in the bottom-right panel of Figure 2. Figure 2 also suggests that the soft X-ray emission has a bipolar morphology. A pair of conical structures open up to the northern and southern directions protruding from the central region. This situation is further illustrated in Figure 1-left that compares the XMM-Newton soft X-ray emission with narrow band images obtained with the VLT, making evident the match between the distributions of the X-ray-emitting gas and the large-scale optical nebular structures. It is worth noting that the EPIC images does not resolve completely the emission in the central region of R Aqr, unlike the Chandra observations that trace in detail the hot gas associated with the bipolar jets of R Aqr (see Kellogg et al. 2007;Nichols et al. 2007). On the contrary, the Chandra data present negligible contribution of extended emission. This situation is illustrated in Figure 1, which shows composite optical and X-ray pictures using XMM-Newton EPIC (left) and Chandra ACIS-S (right) data. We confirm the different spatial-scales sampled by XMM-Newton and Chandra in Figure 3 by comparing the event files of the EPIC and ACIS-S cameras. The Chandra event files clearly resolve the contribution from the SS, the NE clump, and some extended emission towards the SW (see Fig. 3 right panel), which can only be hinted in the XMM-Newton event files. Contours of the soft, medium, and hard emission detected by Chandra and generated with the CIAO tasks aconvolve are also compared with the EPIC images in Figure 2. 3.2. Physical properties of the extended X-ray emission XMM-Newton spectra were extracted from the EPIC pn camera, given its superior effective area compared with the MOS cameras. The spectrum of the extended emission was extracted from a region encompassing the soft X-ray emission. We excluded the central region which is dominated by the emission from the bipolar jet as well as a pair of background point-like sources (see Fig. 3). The background was extracted from a circular region with a radius of 80 towards the SE from R Aqr (see Fig. 3 left panel). The spectrum and associated calibration matrices were created using the SAS tasks eveselect, rmfgen and arfgen. The resultant background-subtracted EPIC pn spectrum of the extended emission is presented in the left panel of Figure 4. The net count rate detected from the extended emission is 23.0±2.3 counts ks −1 .
For comparison we also extracted the spectrum from the brightest X-ray clump detected ∼10 NE of R Aqr and resolved in the Chandra ACIS-S data (see Fig. 3 right panel).
Its extraction region corresponds to a 4 -radius circular aperture, whilst the background was selected from a 15 -radius circular aperture towards the SW. The background-subtracted ACIS-S spectrum and calibration matrices were produced using the CIAO task specextract. The resultant background-subtracted Chandra ACIS-S spectrum of the NE clump is presented in the right panel of Figure 4. The NE clump has a net count rate of 6.9±0.3 counts ks −1 .
Both spectra are extremely soft and resemble those presented and analyzed by Kellogg et al. (2007). They exhibit the presence of the O VII triplet at 0.58 keV and prominent emission at energies below 0.5 keV. The latter includes the very likely emission from N VI at 0.43 keV and the Lyα C VI emission line at 0.37 keV. No significant emission is detected above 1.0 keV.
The spectral analysis was performed using XSPEC (Arnaud 1996).
We started by fitting single plasma emission models to the data, which we note it was enough for the EPIC spectrum, but according to Kellogg et al. (2007) the Chandra spectrum requires two-temperature apec optically-thin emission plasma models. The absorption was included with the tbabs model (Wilms et al. 2000) and all spectra were fit adopting solar abundances from Lodders (2003).
We started by fitting the Chandra ACIS-S spectrum of the NE clump. The best fit model (χ 2 =1.04) resulted in plasma temperatures of kT 1 = 0.09 +0.02 −0.03 keV and kT 2 = 0.24 +0.09 −0.05 keV with normalization parameters 4 of A 1 = 4.3 × 10 −4 cm −5 and A 2 = 1.2 × 10 −5 cm −5 , respectively. The absorption column density of the model was N H = (8 +11 −5 ) × 10 20 cm −2 , consistent with the extinction estimate of A V =1-2 mag by Bujarrabal et al. (2021) towards R Aqr. The total observed flux in the 0.3-1.5 keV energy range is f X = (4.8 ± 0.8) × 10 −14 erg cm −2 s −1 which corresponds to an intrinsic flux of F X = (1.5±0.2)×10 −13 erg cm −2 s −1 and an X-ray luminosity 5 of L X = (2.7±0.4)×10 30 erg s −1 . We note that the softer component contributes to 85% of the total X-ray emission. Finally, we estimated an electron density of n e ≈ 40 cm −3 by adopting an averaged radius of 4 for the NE clump. The model is compared with the observed Chandra ACIS-S spectrum in the right panel of Figure 4 and the details are listed in Table 1.
The spectral fit of the EPIC pn spectrum of the extended emission resulted in very similar properties as those of the NE clump. However, we note that we had to fix the absorption column density, because when left as a free parameter it converged to unphysical values. The N H in the extended regions away from R Aqr should have less extinction than that of the central regions. The HEASARC N H column density webpage 6 estimates that in the vicinity of R Aqr this is ∼ 2×10 20 cm −2 . Thus, we decided to fix the absorption column density to an intermediate value of N H = 4 The normalization parameter can be estimated as A ≈ 10 −14 n 2 e dV /4πd 2 , where ne is the electron number density, V is the volume of the X-ray-emitting region and d is the distance to the object. 5 Here we adopt the geometric Gaia distance, but we note that smaller distance estimates have been reported in the literature. The more recent is that of Liimets et al. (2018) which obtained a distance of ∼180 pc using the expansion parallax method. 6 https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl 4 × 10 20 cm −2 . With this, the best fit model (χ 2 =1.00) has a dominant plasma temperature of kT 1 = 0.09 +0.01 −0.02 keV with normalization parameter of A 1 = 1.6 × 10 −4 cm −5 . The total observed flux in the 0.3-1.5 keV energy range is f X = (4.4 ± 1.5) × 10 −14 erg cm −2 s −1 with an intrinsic flux of F X = (7.7 ± 2.6) × 10 −14 erg cm −2 s −1 , which corresponds to an X-ray luminosity of L X = (1.4 ± 0.5) × 10 30 erg s −1 . The estimated electron density for the most extended emission of n e ≈0.5 cm −3 adopting an averaged angular radius of 2 . The best fit is compared to the EPIC pn spectrum in presented in the left panel of Figure 4 and the model details are also listed in Table 1.
We note that the observed surface brightness of the extended emission detected by XMM-Newton is a few times 10 −17 erg s −1 cm −2 arcsec −2 . This value is just at the detection limit of the ACIS-S instrument, thus explaining why it was not detected by Chandra.
Other models to the X-ray spectra were attempted. For example, we also tried plasma emission models with variable N abundance in order to fit the spectral feature at 0.43 keV, but such models did not improve the fit. Models adopting a power law distribution for the hardest component were also attempted. Such models resulted in good-quality fits (χ 2 ≈ 1) too, suggesting that the physical origin of the second component can not be firmly established. Lastly, we note that other models taking into account non equlibrium thermal components or non equlibrium collisional model are difficult to constrain given the relative poor-quality of the spectra (see Kellogg et al. 2007).

DISCUSSION AND CONCLUDING REMARKS
The analysis of the XMM-Newton EPIC data with the ESAS tasks has revealed the presence of extended emission surrounding R Aqr in the soft (0.3-0.7 keV) X-ray band. The diffuse emission detected by EPIC seems to fill the gaps in between the nebular hourglass structures, extending to distances more than 2 towards the northern and southern regions from R Aqr (see, e.g., Fig. 1 left panel). The extension of the diffuse X-ray emission is consistent with the results discussed in Bujarrabal et al. (2021) where the orbital plane of the SS is seen almost edge-on and most of the dense gas is distributed in the equatorial plane, with lower density along the N-S directions that facilitates the expansion of the hot gas.
We found that the spectrum of the extended X-ray emission detected by XMM-Newton EPIC shares a very cool component with the spectrum of the NE clump associated with the S-shaped jet resolved by Chandra ACIS-S. This suggests that the jet is an ongoing structure that keeps feeding the extended, X-ray-emitting hot bubble. The cooling time of such hot but low density gas is long (see Kellogg et al. 2007) and, accordingly, no significant reduction of the temperature is observed between the jet and the most extended X-ray emission. The estimated n e of the most extended gas is almost 80 times lower than that associated with the NE clump, implying a significantly smaller differential emission measure for the former. Thus, the gas expands adiabatically and cools down.
The nature of the harder component required in the spectral analysis of the X-ray-emitting jet is still unknown. This seems to be only spatially correlated with the collimated jet and it is not surprising that magnetic phenomena are responsible for the shaping (Burgarella et al. 1992;Melnikov, Stute, & Eislöffel 2018). However, the search for non thermal emission in radio bands has been challenging (Hollis et al. 1985;Bujarrabal et al. 2018). Further analysis combining multi-epoch Chandra data will help shed light on the physical properties of the X-ray-emitting gas associated with jets (R. Montez Jr private communication).
The jet appears to be collimated, its S-shape is clearly indicative of precession. The precession angle is wide, ≈ 50 • , as suggested by high-resolution Hubble Space Telescope images (Melnikov, Stute, & Eislöffel 2018), and actually theoretical predictions state that opening angles larger than 40 • produce nebulae with pairs of inflated bubbles (Soker 2004). We therefore suggest that blister-like structures at the tip of the precessing jet get disrupted and feed hot gas into the most extended hot bubble. The creation-destruction of the blister-like structures seems to be relatively fast given the evident expansion of the X-ray-emitting clumps resolved by Chandra (e.g., Kellogg et al. 2007). It is possible that the high pressure of the hot gas helped inflating the most extended lobe detected in optical wavelengths (see figure 2 of Liimets et al. 2018).
We also suggest that the continuous production and disruption of blister-like structures has created the bowl-like structures opening towards the northern and southern regions of R Aqr. Disk formation and jet launching has been studied through numerical simulations of binary systems in great detail (see, e.g., López-Cámara et al. 2020;Sheikhnezami & Fendt 2018;Waters & Proga 2018, and references therein), but the XMM-Newton EPIC observations analyzed here leverage the role of accretion process onto a compact object, a WD in this case, on the formation of hot bubbles and large-scale nebulae. This process might be similar, but in an evidently smaller scale, to the creation of hot bubbles produced by active galaxies hosting supermassive black holes (see, e.g., the case of Cen A; Król et al. 2020;Hardcastle et al. 2003, and references therein).
Several theoretical works have proposed the jet feedback mechanism (JFM) as the physics behind the production of inflated bipolar cavities in a wide variety of astronomical phenomena. The JFM successfully explains a diverse variety of objects regardless of differences of several orders of magnitude in size, energy injection, mass and timescales (Soker 2016), from young stellar objects through planetary nebulae to galaxy clusters (e.g., Cliffe et al. 1995;Soker 2002Soker , 2004Soker , 2016.
Deeper and higher-quality observations, as those provided by eROSITA and future missions such as Athena, might be able to disclose more extended X-ray emission with lower surface brightness of X-ray-emitting gas in R Aqr. A complete mapping of the diffuse X-ray hot bubble will help us study the feedback from SSs to create fast evolving nebulae around them.