VLT MUSE observations of the bubble nebula around NGC 1313 X-2 and evidence for additional photoionization

The bubble nebula surrounding NGC 1313 X-2 is believed to be powered by high velocity winds from the central ultraluminous X-ray source (ULX) as a result of supercritical accretion. With the Multi-Unit Spectroscopic Explorer (MUSE) observation of the nebula, we find enhanced OIII emission at locations spatially coincident with clusters of stars and the central X-ray source, suggesting that photoionization in addition to shock-ionization plays an important role in powering the nebula. The X-ray luminosity of the ULX and the number of massive stars in the nebula region can account for the required ionizing luminosity derived with MAPPINGS V, which also confirms that pure shocks cannot explain the observed emission line ratios.


INTRODUCTION
Ultraluminous X-ray sources (ULXs) are extremely accreting compact objects, and occupy the high end of the luminosity function of high mass X-ray binaries (Mineo et al. 2012). They display X-ray luminosities exceeding the Eddington limit of a typical stellar mass black hole, with spectral and timing behaviors distinct from Galactic X-ray binaries (for a review see Kaaret et al. 2017). Identification of neutron stars in ULXs (Bachetti et al. 2014;Fürst et al. 2016;Israel et al. 2017a,b;Carpano et al. 2018;Sathyaprakash et al. 2019;Rodríguez Castillo et al. 2020) and their ordinary spectral properties (Pintore et al. 2017;Walton et al. 2018) suggest that supercritical accretion occurs in the majority of these systems.
Interaction of the wind with the interstellar medium (ISM) may produce shock-ionized bubble nebulae, which have been found around some ULXs (Pakull & Mirioni 2002Ramsey et al. 2006;Abolmasov et al. 2007;Abolmasov 2008;Russell et al. 2011;Soria et al. 2021). These bubbles have a size of some 100 pc, expanding at a velocity of about 100 km s −1 . Based on the emission line luminosity, bubble size, and expansion velocity, the mechanical power of the wind or jet that drives the bubble is estimated to be about 10 39 -10 40 erg s −1 , suggesting that the accretion is supercritical (Pakull & Mirioni 2003;Abolmasov et al. 2007), in agreement with the detection of diffuse X-ray emission in one of them .
NGC 1313 X-2 is a prototype ULX among the first discovered (Colbert et al. 1995). A bubble nebula is found to be coincident with the X-ray source position and has been extensively studied (Pakull & Mirioni 2002Zampieri et al. 2004;Mucciarelli et al. 2005;Pakull et al. 2006;Ramsey et al. 2006;Ripamonti et al. 2011). The nebula is elongated with a geometry of 590 pc × 410 pc (26 × 18 ). The presence of strong [S II] and [O I] emission lines and the supersonic expansion suggests that the nebula is shock-ionized (Pakull & Mirioni 2002, while the most likely powering source is the wind/jet from the central ULX ) with a mechanical power estimated to be about 1.5 × 10 39 erg s −1 (Pakull & Mirioni 2002).
In this paper, we present a spatially resolved spectroscopic study of the bubble nebula around NGC 1313 X-2 using the MUSE instrument (Bacon et al. 2010) on the Very Large Telescope (VLT). The observations are described in § 2, the data analysis and results are presented in § 3, and the physical nature is discussed in § 4. We assume a distance of 4.6 Mpc to NGC 1313 based on Cepheids measurements (Qing et al. 2015).

OBSERVATIONS
The VLT/MUSE observations of NGC1313 X-2 were performed on the nights of 2019 October 16 and November 12 (UT). The seeing during the observations was from 0. 8 to 1. 0, and sky transparency conditions were clear on 2019 October 16 and partially cloudy on 2019 November 12. The MUSE WFM-NOAO-E mode was used to cover the wavelength between 4600Å and 9350Å. The observations include four 1400-second exposures: two of them with a position angle of 0 • and the rest two with an angle of 90 • .
The data were reduced by the EsoRex pipeline version 3.13.2 (Weilbacher et al. 2020). We used the muse_scibasic recipe to correct bias frames, lamp flats, arc lamps, twilight flats, geometry, and illumination exposures for individual exposures. Then the muse_scipost recipe was used for telluric and flux calibrations with the standard star taken in the same night. At last the individual exposures were aligned using the muse_exp_align recipe and combined to generate the final datacube using the muse_exp_combine recipe. The field of view of the final datacube is 1 × 1 and the spaxel scale is 0. 2 × 0. 2. Figure 1 shows the nebula images in the Hα band, one with line plus continuum emission and the other with line emission only. A bright object appears near the west end of the nebula, and is identified to be a foreground M1 dwarf at a distance of 750 pc (Gaia Collaboration et al. 2021) along the line of sight of the nebula (also see Ramsey et al. 2006).

DATA ANALYSIS AND RESULTS
We find that the nebula spectra are contaminated by unresolved stellar emission in the host galaxy projected on the nebula region (see the Hubble images in Grisé et al. 2008). The stellar emission contributes a continuum component with absorption lines, while the nebula spectrum mainly consists of emission lines. In order to remove the stellar contamination, we extract a template spectrum from a background region defined in Figure 1. For the nebula, we are only interested in the emission line properties in two bands, a blue band that contains Hβ and [O III], and a red band with [O I], Hα, [N II], and [S II]. We extracted the nebula spectrum from an elliptical region defined in Figure 1, and fit it with the template in the two bands, respectively, at wavelengths excluding emission lines. In the fitting, we employ a scale factor and a shift on the flux onto the template spectrum. Then the best-fit stellar template is removed from the nebula spectrum. The blue and red bands are not fitted jointly because different scales are needed. We have to assume that the stellar contamination is insensitive to the sky position, because the statistics from a small number of pixels does not allow us to perform the template fitting. The spectrum extracted from a single pixel after subtracting the stellar template exhibits fluctuations consistent with the statistical error, justifying the validity of the technique. The final spectra in the blue and red bands are shown in Figure 2. Each emission line is fitted locally with a Gaussian to a continuum subtracted spectrum. The continuum is determined with a linear function from nearby wavelengths free of the emission component.  (Storey & Zeippen 2000), which is in good agreement with the ratio derived from independent measurements. The line flux, velocity relative to the host galaxy (z = 0.001568 or v = 470 km s −1 ; Koribalski et al. 2004), and the full-width-half-maximum (FWHM) velocity dispersion corrected for instrument line broadening 1 are derived for bright lines and listed in Table 1  Rest wavelength ( ) 32 . nebula. The line fluxes derived from the observation on October 16 are systematically higher than those on November 12 because of different sky conditions, and have been corrected as in the clear sky condition. The average radial velocity of the whole nebula is around −80 km s −1 (negative means a blue shift) with respect to the bulk velocity of 470 km s −1 as mentioned above. They are consistent with the line velocities measured from three nearby, possibly H II regions. For example, one of them marked by a dashed green circle in Figure 1 has a velocity of −81 ± 2 km s −1 measured with Hα or −80 ± 2 km s −1 with [S II]. Thus, we adopt −80 km s −1 as the systematic velocity. This confirms the conclusion in Pakull & Mirioni (2002) that the bubble nebula is indeed associated with the host galaxy.
To investigate the spatial distribution of the line properties, we produce the flux, radial velocity, and FWHM velocity dispersion (instrument broadening removed) maps in Figure 3. Hα is the only emission line that has a signal-to-noise ratio high enough for this purpose. We employ a 2 × 2 pixel binning for the flux and velocity maps, and a 3 × 3 binning for the FWHM map to improve the statistics.
We mark three particular regions on the map (see Figure 4). Region 1 is a 2 wide elliptical annulus along the edge of the bubble, excluding the eastern end where it overlaps with region 3. This is the outermost region of the bubble, and represents the part with the lowest FWHM. Region 2, i.e., the ULX region, is a 2 -radius circle around the X-ray position, which is blue shifted compared with other regions. Region 3 is a low flux cavity in the eastern part of the bubble. The velocity and and velocity dispersion measured with Hα, [S II], and [O III] in the three regions are listed in Table 2.

Line decomposition
We suppose that the line kinematics seen on each pixel of the nebula is a result of both approaching and receding motions, with respect to the systematic velocity (also see Soria et al. 2021). Region 1, which is along the edge of the bubble, shows a radial velocity well consistent with the systematic velocity, and a FWHM among the lowest in the bubble, indicating that the gases along the edge may be moving transversally. If this is the case, the FWHM in region 1 reflects the intrinsic line width of the shock, which is around ∼70 km s −1 (see Table 2).
Due to the presence of strong low ionization forbidden lines, such as [O I], [S II] and [N II], the shock is believed to be radiative instead of adiabatic (Heng 2010;Soria et al. 2021). For radiative shocks, it is not straightforward to infer the local shock velocity from the gas velocity dispersion (see Appendix A in Soria et al. 2021). Following their recipe, which is based on the assumption of a uniformly expanding   2) and of the whole nebula (Table 1), respectively. Previous studies using long slit spectroscopy argued that the shock velocity of the nebula could be as high as 80 km s −1 (Pakull & Mirioni 2002) or 100 km s −1 (Ramsey et al. 2006). Thus, in this study, we assume a shock velocity of 80 km s −1 . Then, we try to model the observed Hα line profile in the ULX region (region 2), where the bulk motion of shocked gas is radial. Here we assume that the shock velocity equals the bulk velocity of shocked gas, which is the observed velocity, in the case of fully radiative shocks (Dopita & Sutherland 2003). We fix the shock velocities to be ±80 km s −1 for the approaching and receding components, respectively, and assume that both have a FWHM of 70 km s −1 . The instrument line broadening and measurement noise are taken into account. However, the model line width is significantly larger than the observed line width ( Figure 5). We obtain the same conclusions if we fix the FWHM at any value in the range of 50-100 km s −1 .
This may suggest that the emission line has a more complicated velocity distribution (also see the 2D spectrum in Ramsey et al. 2006). Therefore, we add a zero-velocity component with an intrinsic FWHM = 0. The three-component model can fit the line profile adequately, with a flux ratio of 3:3:1 ( Figure 5). The 2D spectrum in Ramsey et al. (2006) also shows a hint of a stationary component. The presence of a low-velocity narrow component could be interpreted as due to photoionization.

As mentioned above, the strong [S II] and [O I] emission
lines in the nebula are suggestive of predominant shockionization (Pakull & Mirioni 2002. To produce high ionization emission lines like [O III], one requires either high velocity shocks or photoionization. We display the [O III] λ5007 flux map on top of an Hubble Space Telescope (HST) ACS F435W image in Figure 6. As one can see, regions with high [O III] fluxes are spatially coincident with dense stellar populations shown on the HST image. No similar coincidence is seen with the Hα map. This implies that some of the [O III] emission may be attributed to photoionization by massive stars. The [O III] flux in the ULX region (region 2) is also enhanced, possibly caused by additional X-ray irradiation. As one can see in Table 2, the FWHM of [O III] in region 2 is significantly smaller than that of Hα and [S II], implying that the [O III] emission is perhaps due to a different ionizing mechanism.

Comparison with MAPPINGS V
To check if photoionization is needed in addition to shock-ionization, we perform simulations with MAPPINGS V (Sutherland & Dopita 2017). We assume an ISM density of 1 cm −3 , a magnetic field of 3 µG (Han 2017); these two parameters are not sensitive to the results. We also assume Z = 0.5Z as suggested by previous studies (Ripamonti et al. 2011;Pintore & Zampieri 2012) and will discuss its influence. The models include dust calculations and allow grain destruction; these two options also have a small impact on our conclusions. The pre-ionization state is calculated in a self-consistent manner, taking into account inter- nal photoionization due to post-shock radiation and external photoionization from X-ray sources and O stars. First, we compare observations with simulations of pure shocks (Model A), and set the shock velocity v s = 80 km s −1 . However, the simulated [O III] λ5007 to Hβ flux ratio (on the order of 10 −6 ) is significantly lower than that observed (see Table 3). If we set a solar abundance, [O III] λ5007/Hβ is found to be 0.016, still about two orders of magnitude lower than that observed. To match the observed [O III] λ5007/Hβ, one requires a shock velocity over 100 km s −1 at 0.5Z , or ≈ 85 km s −1 at solar abundance. In these cases, however, the [O I] λ6300 to Hα flux ratio is found to be significantly lower than that observed, because most of the oxygens stay at a higher ionization state. To conclude, pure shock-ionization cannot account for both [O III] λ5007/Hβ and [O I] λ6300/Hα.
We thereby add photoionization in addition to shockionization (Model B). We add an X-ray source and a certain number of O stars for photoionization illuminating at a distance of 300 pc (length of the semi-major axis of the nebula). The X-ray power-law spectral index is fixed at −1.5, the typical value found from X-ray fitting (Qiu & Feng 2021). Each O star has an effective temperature of 40000 K. We vary the X-ray luminosity and O star quantity to fit the observation, and find that an X-ray luminosity of 1.6 × 10 40 erg s −1 plus 60 O stars with a total luminosity of 6.4 × 10 40 erg s −1 can reasonably match the observations (Table 3). The simulated [N II] λ6583/Hα is higher than that observed; this is consistent with previous studies by Ripamonti et al. (2011) who propose that the nitrogen abundance in the nebula could be further lower. This ULX also displays a soft blackbody component, which could be the ionizing source. We replace the powerlaw component with a blackbody component and set the temperature at two extremes from observations, 0.1 keV and 0.25 keV (Qiu & Feng 2021). However, it cannot fit the observed flux ratios in the luminosity range of (0.5 − 50) × 10 39 erg s −1 . If we further increase the blackbody luminosity to match the observed [O III] λ5007/Hβ, the predicted [S II] λ6716/Hα is way below that observed.
We also perform a test with a flat ionizing spectrum (F ν ∝ ν 0 ) in the EUV to soft X-ray band (40 eV to 1 keV) to mimic emission from a multicolor disk with an innermost temperature of about 0.1 keV. In this case, an X-ray luminosity of 5 × 10 40 erg s −1 plus 250 O stars can provide a good fit with the observations except [O I] λ6300/Hα, which is overpredicted by a factor of 1.6. Thus, we refer to this case as a marginal fit.
The observed Hα to Hβ flux ratio is 2.87 ± 0.14 after correction with the Galactic extinction E(B −V ) = 0.075 along the line of sight (Schlafly & Finkbeiner 2011). This sets the upper bound on the intrinsic Hα/Hβ. As one can see, the pure shock model (Model A) predicts a higher Hα/Hβ, inconsistent with the observation, while the model with additional photoionization (Model B) produces an acceptable ratio, if there is no extra reddening in the host galaxy. This is consistent with previous studies ) that there is no extragalactic extinction for the NGC 1313 X-2 nebula.
In Figure 7 The top panel shows the continuum removed [O III] λ5007 (±5Å around the line centroid) flux images around NGC 1313 X-2, with a binning of 9 × 9 pixels. Those having a signal-to-noise ratio less than 2 are not shown. The middle and bottom panels display an HST ACS F435W image with [O III] λ5007 or Hα flux contours, respectively.
[O III] λ5007/Hβ need more contribution from photoionization, including the ULX region and the southeastern and northwestern parts of the nebula. For comparison, the [O I] λ6300/Hα map is also shown in Figure 7. The distribution of [O I] λ6300/Hα is anti-correlated with that of [O III] λ5007/Hβ, as they require oxygens at different ionization states. Particularly, the southwestern region shows high [O I] λ6300/Hα and low [O III] λ5007/Hβ. We try to model the emission line ratios in this region with Model B, and find that less X-ray luminosity (7 × 10 39 erg s −1 ) and less number (20) of O stars are needed, compared with that needed for the whole nebula (see Table 3), indicative of a lower pre-ionization state in this region.

Total power
Following Pakull & Mirioni (2002), here we estimate the mechanical power of the wind/jet that inflates the bubble. The Hβ surface brightness (S Hβ ) is a function of the shock velocity (v s ) and the pre-shock ISM number density (n ISM ; Dopita & Sutherland 1996), (1) Observationally, the surface brightness can be calculated via the line luminosity (L Hβ ) or observed line flux (f Hβ ), where R B and θ B are the physical radius and angular radius, respectively, of the nebula. Here we take the size of the semi-major axis of the nebula into calculation. We take a shock velocity v s = 80 km s −1 and assume Galactic extinction only. Thus, for the whole nebula, we obtain an average n ISM = (0.45 ± 0.02) cm −3 . Assuming a pressure driven nebula (Weaver et al. 1977), the wind/jet power (P ) and the age of the nebula (t) are related with other properties of the nebula as follows, where ρ 0 = µm p n ISM is the pre-shock density of the ISM, m p is the proton mass, and µ ≈ 1.2 is the average atomic weight assuming neutral gas with Z = 0.5Z . Plugging in R B and n ISM and considering the measurement uncertainties, we get t ≈ 2.1 Myr, and P = (6.5 ± 0.3) × 10 39 erg s −1 .
4. DISCUSSION With the MUSE data, we obtain results consistent with previous studies for the whole nebula. The flux, velocity, velocity dispersion maps help reveal interesting structures inside the bubble. The surface brightness is scaled with shock velocity and ISM density. Therefore, the low surface brightness in region 3 can be explained as a result of a low ISM density.
In region 2 (ULX region), where the gas bulk motion is radial, we find that a simple two-component (approaching and receding) model cannot fit the line profile, and a third low-velocity narrow component is needed. Such a component could be a hint of photoionization. We find spatial correlation between the [O III] emission and clusters of stars, as well as enhanced [O III] emission in region 2 around the ULX, suggesting that some of the [O III] emission may be due to photoionization. Simulations with MAPPINGS V also suggest that pure shocks are unable to account for the observed line ratios; shocks with ≤80 km s −1 cannot produce [O III] λ5007/Hβ as high as that observed, while shocks with higher velocities that are able to reproduce the observed [O III] λ5007/Hβ cannot explain the observed high value of [O I]/Hα due to over-ionization. Adding photoionization to shocks of 80 km s −1 can reproduce the observations. In that case, the X-ray luminosity is derived to be close to 10 40 erg s −1 , which is comparable to but slightly higher than the typical isotropic luminosity (several times 10 39 erg s −1 ) of NGC 1313 X-2 inferred from X-ray observations (Qiu & Feng 2021). We note that both luminosities are derived assuming an illuminating distance of 300 pc. If the emission line clouds are located at a smaller distance to the source, a lower luminosity is required. For X-rays, a reduced distance by a factor less than 2 can match the luminosity seen from X-ray observations. For O stars, as the brightest [O III] region is spatially coincident with star clusters, they could be close in distance and only a few of them can account for the needed luminosity (e.g., one or two O stars at a distance of 50 pc). Therefore, the presence of these ionizing sources can account for the power needed for photoionization. However, the power-law component extrapolated from the X-ray band may cut off at energies above EUV, while the blackbody component observed in the soft X-ray band cannot provide sufficient ionization. Alternatively, a flat ionizing spectrum with a higher X-ray luminosity (5 × 10 40 erg s −1 ) plus more O stars (∼250), both at 300 pc, can provide a marginal fit. If this is the case, it may suggest that the X-ray ionizing source is a multicolor accretion flow and the ionizing distance is smaller than that assumed, e.g., 1.4 × 10 39 erg s −1 at 50 pc. We note that, although the total luminosity seems to fit the observation, one should keep in mind that the origin of the X-ray source for ionization is still uncertain.
Additional photoionization may affect the the estimation of the wind/jet power that is based on a shock model (Weaver et al. 1977). However, the current data do not allow us to accurately quantify the luminosity and spectrum of the photoionization source. Also, the extra photoionization changes the pre-ionization state and has a nonlinear contribution to the production of emission lines. Therefore, we may treat the inferred power as an upper limit of the wind/jet power for the ULX.