Quasi-isotropic UV Emission in the ULX NGC~1313~X--1

A major prediction of most super-Eddington accretion theories is the presence of anisotropic emission from supercritical disks, but the degree of anisotropy and its dependency with energy remain poorly constrained observationally. A key breakthrough allowing to test such predictions was the discovery of high-excitation photoionized nebulae around Ultraluminous X-ray sources (ULXs). We present efforts to tackle the degree of anisotropy of the UV/EUV emission in super-Eddington accretion flows by studying the emission-line nebula around the archetypical ULX NGC~1313~X--1. We first take advantage of the extensive wealth of optical/near-UV and X-ray data from \textit{Hubble Space Telescope}, \textit{XMM-Newton}, \textit{Swift}-XRT and \textit{NuSTAR} observatories to perform multi-band, state-resolved spectroscopy of the source to constrain the spectral energy distribution (SED) along the line of sight. We then compare spatially-resolved \texttt{Cloudy} predictions using the observed line-of-sight SED with the nebular line ratios to assess whether the nebula `sees' the same SED as observed along the line of sight. We show that to reproduce the line ratios in the surrounding nebula, the photo-ionizing SED must be a factor $\approx 4$ dimmer in ultraviolet emission than along the line-of-sight. Such nearly-iosotropic UV emission may be attributed to the quasi-spherical emission from the wind photosphere. We also discuss the apparent dichotomy in the observational properties of emission-line nebulae around soft and hard ULXs, and suggest only differences in mass-transfer rates can account for the EUV/X-ray spectral differences, as opposed to inclination effects. Finally, our multi-band spectroscopy suggest the optical/near-UV emission is not dominated by the companion star.


INTRODUCTION
Ultraluminous X-ray sources are defined as extragalactic off-nuclear point-like sources with an X-ray luminosity exceeding the Eddington limit of a 10 M ⊙ black hole (BH) (e.g.Kaaret et al. 2017;King et al. 2023).It is now established that the vast majority of these systems are powered by super-Eddington accretion onto a stellar-mass compact objects in binary configurations with a donor star.While it is speculated that the population of ULXs might be dominated by transient systems briefly reaching the ULX threshold (Brightman et al. 2023), most of the well-known systems shine persistently at such extreme luminosities, acting as laboratories for the study of sustained super-Eddington accretion.However, despite decades of studies, how such extreme luminosities are produced remains a matter of debate.
One major prediction of super-Eddington accretion theory is the presence of highly anisotropic emission (Shakura & Sunyaev 1973;Poutanen et al. 2007).As the mass-transfer rates reaches or exceeds the Eddington limit, powerful radiation-driven optically-thick outflows are launched from the accretion disc, creating an evacuated cone or funnel around the rotational axis of the compact object.This causes observers at high-inclination to see the reprocessed emission of the outflow photosphere, whereas observer peering down the funnel will see the hot emission from the inner parts of the accretion flow (Poutanen et al. 2007;Abolmasov et al. 2009).Even if ★ E-mail: a.gurpide-lasheras@soton.ac.uk the quantitative details may differ, there is now a body of numerical simulations which indeed reproduce the accretion flow geometry and anisotropic emission pattern envisioned by Shakura & Sunyaev (1973) (e.g.Kawashima et al. 2012;Narayan et al. 2017;Mills et al. 2023), .
The discovery of neutron stars (NSs) in ULXs through X-ray pulsations (Bachetti et al. 2014;Israel et al. 2017;Castillo et al. 2020) opened up new avenues in which super-Eddington accretion may proceed.For instance, it is known that magnetic fields reduce the cross-section for electron scattering, thereby increasing the allowed Eddington luminosities (Basko & Sunyaev 1976;Mushtukov et al. 2015).NSs are also expected to be more radiatively efficient compared to BHs, as the latter swallow the excess radiation which is emitted otherwise at the NS surface (Takahashi et al. 2018).Additionally, at high-mass transfer rates, the NS may be engulfed in an optically-thick magnetosphere, whose spectrum is predicted to emit instead as a multi-color blackbody with a dipolar temperature dependency (Mushtukov et al. 2017).
Constraining the degree of anisotropy in ULXs is thus not only key to understand the accretion flow geometry powering them and test existing theories, but also imperative to understand the effect of the ULX on its environment.For instance, the exact radiative output will inform about the role of ULXs/X-ray binaries on the epoch of re-ionization (Madau & Fragos 2017) as well as help explain how or whether ULXs can shape some galaxy properties.In this regard, explaining the presence of bright nebular Heii 4686 emission line in the integrated spectra of metal-poor galaxies remains a long-standing issue, as regular stars do not produce enough photons above its ionisation potential (IP = 54 eV) to explain it (Schaerer et al. 2019, and references therein).The discovery of HeIII regions around a few ULXs (Pakull & Mirioni 2002;Kaaret et al. 2004;Abolmasov et al. 2008) together with the fact that this line seems more prevalent in low-metallicity galaxies, where hard-ionising sources such as X-ray binaries and ULXs are more common (e.g.Shirazi & Brinchmann 2012;Kovlakas et al. 2020;Lehmer et al. 2021), has made ULXs receive attention as a potential explanation for this so-called 'Heii problem' (e.g.Simmonds et al. 2021;Kovlakas et al. 2022).Whether that is the case remains uncertain, mainly due the poor understanding of the ULX UV emission and anisotropy effects (Simmonds et al. 2021;Kovlakas et al. 2022).
A key discovery allowing to constrain the degree of anisotropy was the observation of extended (25-80 pc) EUV/X-ray photoionized gas around a handful of ULXs (Pakull & Mirioni 2002).Such nebulaewhen spatially-resolved -effectively provide a 2D map of the ionising SED not directed onto the line of sight, allowing observers to compare the nebular emission lines with those expected from the line-of-sight SED.If the expected emission lines and the observed ones from the nebula are comparable, then the degree of anisotropy must be relatively small.
Most works to date have focused on the HeIII region around the ULX Holmberg II X-1 (e.g.Pakull & Mirioni 2002;Kaaret et al. 2004;Berghea et al. 2010;Berghea & Dudik 2012).These works have found the nebula sees similar SED as that observed along the line of sight, arguing for isotropic emission.It must be noted however that most of these works were based on optical observations, with particular focus on the Heii 4686 line, which is most effective in probing the extreme-UV (EUV) emission.Instead, theoretical/numerical works suggest anisotropy must be strongest in the X-ray band (e.g.Poutanen et al. 2007;Narayan et al. 2017), although observations of the Galactic edge-on ULX-like system SS433 suggest collimation along the polar funnel may well take place in the EUV too (Waisberg et al. 2019).A crucial aspect is therefore the exact extrapolation of the X-ray spectrum to the unaccesible EUV.Evidence suggest that in general, direct extrapolation of X-ray derived models to the UV band is not accurate for ULXs (Dudik et al. 2016;Abolmasov et al. 2008).A way forward is therefore to combine broadband spectroscopy with nebular observations, of which only a few studies exist (Berghea et al. 2010;Berghea & Dudik 2012).
In this work, we attempt to constrain the degree of anisotropy in the ULX NGC 1313 X-1 using the ∼200 pc photoionized nebula we discovered in our previous work (Gúrpide et al. 2022).Here we combine state-resolved multi-band spectroscopy -which allows to us to constrain the SED along the line of sight and reduce uncertainties related to the extrapolation to the UV -with spatially-resolved line maps from IFU spectroscopy (Gúrpide et al. 2022) -which allow us to constrain the SED seen by the nebula along two different sight lines.
We will show that the degree of collimation of the UV/EUV emission is small (a factor ∼4) in agreement with previous works (e.g.Kaaret et al. 2004).However, we will show that unlike the ULXs Holmberg II X-1 or NGC 6946 X-1, NGC 1313 X-1 does not produce a strong HeIII region.We will argue that the reason most likely lies in the differences in mass-accretion rate between these three sources.
This paper is structured as follows: in Section 2 we present our multi-band data reduction and its classification into the spectral states of NGC 1313 X-1 based on the long-term behaviour of the source.In Section 3 we present spectral modelling of the multi-band spectral states of NGC 1313 X-1.Section 4 presents the Cloudy photoionization modelling of the nebular emission along two different sight lines and finally, Sections 5 and 6 present our Discussion and Conclusions.

DATA REDUCTION
In order to characterize the broadband SED of NGC 1313 X-1, we have investigated the available archival data to characterize the source temporal and spectral properties.In particular, we have considered the long-term monitoring provided by the Swift-XRT (Burrows et al. 2005), optical photometry from the Hubble Space Telescope and Xray spectroscopy from XMM-Newton (Jansen et al. 2001) and the NuSTAR (Harrison et al. 2013) observatories.
The long-term lightcurve of NGC 1313 X-1 along with the hardness-intensity diagram (HID) is shown in Fig. 1 and was extracted using the standard online tools (Evans et al. 2007(Evans et al. , 2009)).From the Figure, it is clear how the source transits through two main states (a behaviour also supported by XMM-Newton observations; Gúrpide et al. 2021a;Pintore & Zampieri 2012).We refer to them as the 'high' (≳ 0.25 Swift-XRT ct/s) and 'low' (< 0.15 ct/s) states, respectively, hereafter.Such bi-modality is commonly observed in other ULXs (Luangtip et al. 2016;Amato et al. 2023) but it appears less strong in the case of NGC 1313 X-1 (Weng & Feng 2018).As we discuss below in Section 4, our nebular modelling is insensitive to short-term changes in flux.Therefore we focused on building the broadband SED of the source for these two broadly defined spectral states.

X-rays
X-ray spectral products from XMM-Newton and NuSTAR were taken from Gúrpide et al. (2021a).In particular, based on the recurrent behaviour of NGC 1313 X-1 from the Swift-XRT lightcurve (Fig. 1) and the analysis from Gúrpide et al. (2021a), we selected XMM-Newton and NuSTAR observations taken in 2012 (corresponding to 10XN and 11XN in the notation used by Gúrpide et al. (2021a) or XN1 in Walton et al. (2020a)), March of 2017 (14XN or XN3 respectively) and August/September 2017 (17XXN or XN5 respectively) to characterise the low state.For the characterisation of the high state, we extracted the first 10 ks of XMM-Newton obsid 0803990101 together with the first 20 ks (owing to the lower number of counts) of NuSTAR obsid 30302016002 taken in June of 2017 (15XN or XN4).More details on the characterization of the high state are provided in the Appendix A.
We noticed recent important calibration updates for NuSTAR, so we reduced the data using nuproducts version 2.1.2with the most recent calibration files as of February of 2023.Source and background regions were extracted following Gúrpide et al. (2021a).All XMM-Newton (EPIC-PN and the MOS cameras) and the NuSTAR spectra were rebinned using the scheme proposed by Kaastra & Bleeker (2016) and fitted over the band where the source dominated above the background (this was the ∼ 0.3-10 keV band for the XMM-Newton data and typically up to 20-25 keV for the NuSTAR data).All spectra had sufficient number of counts per bin to use the  2 statistic.
In order to extract fluxes in the different filters from the optical counterpart identified by Yang et al. (2011), we performed aperture photometry using a 0.2"-radius circular aperture in the ACS/WFC and WFC3/UVIS detectors (corresponding to 4 and 5 pixels respectively) and 0.15" (6 pixels) for the ACS/HRC.The aperture centroid was determined using 2D Gaussian fitting.We then corrected the counts for the finite aperture, following a two step process, as recommended2 : first, because the aperture correction at small scales is known to vary with time and location on the detector, we estimated the aperture correction to 10 pixels using isolated bright stars in the field, typically selecting 1 to 4 stars (depending on the availability of isolated stars in the field) and averaging the results when possible.Next, we corrected the 10-pixel fluxes to infinity using the tabulated values for each combination of detector and filter3 .In instances where it was not possible to find any suitable star to estimate the aperture correction, we relied on the tabulated values for the full correction.
The background level was determined by taking the 3--clipped median count rate in a concentric annulus around the source region.For the UVIS filters these regions contained several bright stars, so we instead used a nearby ∼0.7" circular region relatively free of stars but still containing some of the galaxy diffuse emission.The uncertainties on the final background-subtracted count rates where derived assuming Poisson statistics for the source and background regions and accounting for the uncertainty on the aperture correction.
In order to derive extinction-corrected fluxes, we first defined a likelihood function as the sum of the uncertainty-weighted residuals between the background-subtracted count rates and the predicted count rates by a model when convolved with the corresponding HST filter.Along with the source model parameters (described below), we modelled extinction using two components: a Galactic component using Cardelli et al. (1989) reddening law, with  V = 3.1 and an additional, extragalactic component with  V = 4.05 using the extinction curve from Calzetti et al. (2000), which may be appropriate for a star-forming galaxy such as NGC 1313.The galactic component was fixed to the value along the line of sight ( ( − ) G = 0.11; Schlafly & Finkbeiner 2011) while the extragalactic component was included as a Gaussian ( E(B−V) = 0.15,  E(B−V) = 0.03) prior on  ( − ).The constraints on the extragalactic extinction come from our MUSE spectrum and are justified later in Section 3.1, where we show there is additional extinction ( v = 0.59±0.11mag) towards NGC 1313 X-1.By including  ( − ) in this manner we were able to potentially constrain it further based on the HST filters or in the worst case scenario, propagate its uncertainties to the final deabsorbed-fluxes.However, in all instances we could not constrained further and obtained our prior back in the posteriors of  ( − ).
For the source model, we assumed the same absorbed powerlaw (  =   where   are the fluxes in erg/cm 2 /s/Å and  are the pivot wavelengths of the filters) for all filters in a given epoch.Exceptions to this were epoch 2004-02-22 where the fluxes clearly deviated from a single powerlaw ( 2 > 20 for 1 degree of freedom; see also figure 4 in Yang et al. 2011) and epochs for which only one measurement was available.In these cases we assumed a flat spectrum in erg/cm 2 /s/Å, parametrised by its amplitude .
To estimate the best-fit model parameters and their uncertainties, we drew parameter samples (either , ,  ( −) or , ( −)) to sample the posteriors using 32 Markov-Chain Monte Carlo (MCMC) chains using the emcee package (Foreman-Mackey et al. 2013).The chains were run until a) the number of steps reached 100 times the integrated autocorrelation time (), which was estimated on the fly every 800 samples, and b)  changed less than 1% compared to the previous estimate.We then burned the first 30 samples and thinned the chains by /2.The final intrinsic fluxes and its uncertainties were estimated by drawing 2,000 realizations of , , (or  in cases of a flat spectrum) from the posteriors to estimate the distribution of best-fit intrinsic fluxes in each filter.The mean and the 1 confidence interval of the distribution were taken as our final estimates (these posteriors were symmetric and Gaussian-like in all instances).The final extinction-corrected ST magnitudes and fluxes are reported in Table 1.
While the exact values obviously differ from those reported by Yang et al. (2011) due to the different treatment of the extinction, our analysis is in agreement with the variability reported in the F555W filter by Yang et al. (2011).The rest of the filters for which multiepoch data exist are approximately consistent within uncertainties.We verified that correcting only for foreground extinction as in Yang et al. (2011) yielded fluxes in good agreement with their results.
We also attempted to supplement the multi-band data with UV fluxes from the optical monitor (OM) onboard XMM-Newton.We ran omichain on obsid 0803990101 and found NGC 1313 X-1 is detected with a significance above the 8 level in the UV filters, while it is undetected in the optical filters.We used the default aperture of 6" in radius and converted the instrumental-corrected count-rates to fluxes using the average tabulated values4 .However, upon comparison of the OM fluxes with those from HST in similar bands we found the OM fluxes (≳10 −16 erg/cm 2 /Å/s) to be 2 order of magnitude overestimated, owing to the large aperture which likely contains significant stellar contribution.We concluded that the contribution of NGC 1313 X-1 to the detection must be minimal and discarded the OM data from further analysis.

Source spectral states
In order to place the HST observations in context with respect to the X-ray spectral states, we looked in the archives for X-ray data taken simultaneously with the HST observations.We found three Chandra observations simultaneous with the first two sets of HST observations (taken in November 2003 andFebruary 2004), which were presented also in Yang et al. (2011).The details of the Chandra observations are given in Table 2.The 2014 WFC3/UVIS observations were instead covered by the Swift-XRT long-term monitoring (Fig. 1).
The Swift-XRT places the HST 2014 observations when NGC 1313 X-1 was in the low state.The Chandra data simultaneous with the other HST observations requires more detailed modelling, particularly given the presence of pile-up in some observations.We therefore relegate the full analysis to the Appendix (Section A) but briefly, we found that despite the presence of pile up, we can confidently place NGC 1313 X-1 in the high state during the November 2003 observations, and in the low state during the February 2004 observations.The unabsorbed fluxes and luminosities for the best fits for both Chandra observations are reported in Table 2 along with the determined spectral state.Fig. 1 shows the two equivalent Swift-XRT count rates determined from the Chandra observations on the Swift-XRT lightcurve.
As a summary, Table 3 presents all the data gathered to characterized the broadband SED for the low and high states, respectively.We proceed to characterise the state-resolved broadband SED of the source in the following Section.

MULTI-BAND SED MODELLING: WHAT WE SEE
Our aim here is mostly focused on constraining the EUV emission along the line of sight by testing different spectral models.We will then test these models using Cloudy photo-ionization modelling of the emission-line nebula.Additionally, we wish to determine whether the optical emission is dominated by the companion star or reprocessing in the outer disc, which remains unclear (e.g.Grisé et al. 2012) and can only be tackled with strictly simultaneous broadband data, of which only a few studies exist (e.g.Soria et al. 2012;Sathyaprakash et al. 2022).

Extinction
Modelling of the optical data requires knowledge of the level of extinction towards NGC 1313 X-1.We used the MUSE data presented Notes.Magnitudes have been corrected for Galactic extinction along the line of sight and for additional extinction towards NGC 1313 X-1.Uncertainties represent the 1 confidence level and take into account the uncertainties on the aperture correction, the best-fit model parameters and the extinction used to derive the intrinsic fluxes.'−' Indicates there is no X-ray simultaneous information to categorise the source in one of its two states. Aperture correction.
Total extinction used for the correction estimated from the Balmer decrement.In all fits  H was frozen to the value derived by Gúrpide et al. (2021a).Uncertainties are given at the 90% confidence level for one parameter of interest.
Photon index of the powerlaw or temperature of the diskbb component.
Unabsorbed luminosity over the same band. This observation was not used as it was simultaneous with obsid 4750, which instead was free of pile up.
in Gúrpide et al. (2022) to estimate the level of extinction from the ratio of the Balmer lines (H and H).To this end, we extracted an average spectrum from cube 1 in that work from a circular region of 1" (corresponding roughly to the PSF FWHM) around the optical counterpart.We then corrected the spectrum from Galactic extinction using the Cardelli et al. (1989) extinction curve with  V = 3.1.From this foreground-extinction-corrected average spectrum we measured (H) = 4.2±0.1 and (H) = 14.2±0.1,both in units of 10 −18 erg/s/cm 2 , using a simple constant model for the local continuum around each line and a Gaussian for the line itself.The Balmer decrement (H)/(H) = 3.39±0.11suggests there is additional extinction towards NGC 1313 X-1.The extinction correction requires knowledge of the intrinsic Balmer decrement, which depends on the electron density ( e ) and temperature of the gas () (Osterbrock & Ferland 2006).While we do not have access to temperature-sensitive line diagnostics, the typical values we found in these regions for the electron-density sensitive line ratio were 3. These values correspond to the low-density regime which is nearly temperature-insensitive, indicating  e < 100 cm −2 .Assuming case B recombination and  = 10000 K (broadly consistent with our nebular modelling; Section 4) the intrinsic ratio is H/H= 2.863 (noting the likely absence of shocks in this region also constrains  < 20, 000 K; Osterbrock & Ferland 2006).Therefore we found  ( − ) = 0.15 ± 0.03 adopting the Calzetti et al. (2000) extinction curve with  V = 4.05 as stated in Section 2. The total absorption is thus  v = 0.93 ± 0.11 mag.
Another argument supporting the additional level of extinction can be made considering the relationship between the neutral hydrogen absorption column ( H ) and reddening found by Güver & Özel (2009) using simultaneous X-ray and optical observations of Galactic supernova remnants: Using the Galactic value along the line of sight  ( − ) = 0.11 (Schlafly & Finkbeiner 2011) with  v = 3.1 would suggest  H ∼ 7.5 × 10 20 cm −2 , which is an order of magnitude below the value ).We will additionally show that all models used below also require  H > 2 × 10 21 cm −2 .Instead, using the total  V derived above suggests  H = (2.1 ± 0.3) × 10 21 cm −2 , which is consistent within the order of magnitude with the values derived from X-ray spectral fitting and the values derived below.This suggest our estimate for the amount of reddening is reasonable.

Spectral fits
To model the spectra we decided to use two physically motivated models along with a third phenomenological model typically used to describe the 0.3-25 keV emission in ULXs, which we summarize below.The low-state presented the well-known strong residuals at around 1 keV (Middleton et al. 2015b) which have been associated with radiatively-driven relativistic outflows (Pinto et al. 2016(Pinto et al. , 2020)).These residuals have a minor impact on the continuum estimation but they can affect the estimation of the level of neutral absorption.Moreover, because we wanted our  2 to reflect closely the goodness of fit, we modelled the residuals with a Gaussian at  ∼ 1keV with  ∼0.09 keV, which resulted in Δ 2 improvements upwards of 100 (depending on the exact model) for 3 degrees of freedom.Neutral X-ray absorption was modelled with a tbabs model component as explained in Appendix A. Whether  H is variable in ULXs remains a contested issue virtually unstudied owing to the lack of physically motivated models.In NGC 1313 X-1, Gúrpide et al. (2021a), using phenomenological models, showed  H varies only in an unusual state of the source -termed 'obscured state' in that work (see also Middleton et al. 2015b).Since this state is not considered here, we decided to tie  H between the high and low states, while leaving the rest of the parameters free to vary unless stated otherwise.
Below, we summarise the models employed and the resulting fits to the data.Results from the spectral fitting are reported in Table 4 and Fig. 2 shows the best-fit spectral models and residuals.
• diskir: The diskir model is an extension of the multicolour disc blackbody (Mitsuda et al. 1984) which takes into account selfirradiation effects by a Compton tail and by emission from the disc inner regions (Gierliński et al. 2009).The most relevant parameters for the UV/optical data are the reprocessed flux fraction which is thermalized in the disc (  out ) and the outer disc radius ( out ), whereas the rest of the parameters are constrained by the X-ray data.This model has been used to fit the broadband spectrum of ULXs by analogy with XRBs (Grisé et al. 2012;Tao et al. 2012;Vinokurov et al. 2013;Sutton et al. 2014;Dudik et al. 2016) finding satisfactory agreement with the data, although always with some degeneracy with the companion star (Grisé et al. 2012;Tao et al. 2012).It must be noted that most of these works did not use simultaneous data, which has been shown to be crucial to discriminate between models (e.g.Dudik et al. 2016).Compared to existing works (Grisé et al. 2012;Vinokurov et al. 2013;Dudik et al. 2016) where some parameters were frozen (namely those concerning the Compton tail: the electron temperature k e and the ratio of luminosity in the Compton tail to that of the unilluminated disc,  C / D ) we were able to constrain all parameters simultaneously thanks to the high-energy coverage provided by the NuSTAR data.We found good consistency between the  C / D parameter between the high and low states, which is unsurprising given the lack of variability of NGC 1313 X-1 above ∼10 keV (Gúrpide et al. 2021a;Walton et al. 2020a).We therefore tied this parameter between the high and the low states.The fit provides a good description of the continuum with  2 /dof = 1616/1264 as can be seen from Fig. 2. We note however, that the reprocessing fraction for the low state (  out ∼ 4 × 10 −2 ) is about an order of magnitude higher than found in XRBs (see Urquhart et al. 2018, and reference therein).In using this model, we do not necessarily have a physical interpretation in mind, particularly because we expect the accretion flow geometry to deviate substantially from the sub-Eddington accretion flow geometry assumed by the model.We rather consider the diskir as a proxy to obtain a realistic extrapolation to the inaccessible EUV.
• sirf: The self-irradiated multi-color funnel can be considered an extension of the diskir to the supercritical regime, in which the disc/wind acquires a cone-shaped geometry altering the selfirradiation pattern (Abolmasov et al. 2009).As a caveat, we note that this model does not include Comptonization or the irradiation of the disc wind onto the outer disc itself (e.g.Middleton et al. 2022).
For simplicity, we fixed some parameters which we considered as nuisance for our analysis to their fiducial values.These were the velocity law exponent of the wind, which we fixed to −0.5 (i.e.parabolic velocity law) and the adiabatic index () which we fixed to 4/3 as the gas is radiation-pressure supported.Leaving the outflow photosphere  out free for both states yielded a loosely constrained and highly degenerate fit.We further noted that this parameter did not strongly affect the fits.We managed to obtain successful fits by leaving  out free for the high state, while for the low state, we fixed  out to 100 R sph .We also included self-irradiation effects and found 4-5 iterations were sufficient to reach convergence.
The fit was insensitive to the inclination  so as long as  was smaller than the half-opening angle of the funnel  f so we fixed it to 0.5 • (i.e.nearly face-on) as expected for a ULX such as NGC 1313 X-1 (Middleton et al. 2015a;Gúrpide et al. 2021a).Solutions with  >  f were statistically excluded by the data.
The model fitted reasonably the continuum (Fig. 2), although substantially worse than the diskir ( 2 /dof = 1714.2/1268).This model cannot reproduce the fluxes redward of ∼ 5555 Å particularly in the low state (where the residuals  ≳7).Perhaps unsurprisingly the parameters are similar to those found for Holmberg II X-1 and NGC 5204 X-1 (Gúrpide et al. 2021b) but with a narrower opening angle of the funnel ( f ∼ 37 • ) and a lower Eddington mass ejection rate  eje compared to Holmberg II X-1.
As can be seen from Fig. 2 the main difference with respect to the diskir model is the predictions in the EUV: the sirf predicts a UV flux more than an order of magnitude higher than the diskir, with the UV flux peaking around the He + ionization potential of 54 eV (Fig. 2).Contrary to the diskir model, this model predicts little difference in the UV/optical emission between the high and low states.
• phenomenological: Finally, we tested a phenomenological model based on an absorbed dual thermal diskbb including upscattering of the hard diskbb through the empirical simpl model (Steiner et al. 2009) (see e.g.Gúrpide et al. 2021a;Walton et al. 2020a).The complete model in XSPEC was tbabs⊗(diskbb + simpl⊗diskbb).Here too we tied the photon index Γ between the high and low states owing to the aforementioned lack of variability above ∼10 keV.Because the model make no account for the optical/near UV emission, we fitted this model to the X-ray data only.Fig. 2 shows this model severely underestimates the UV/optical fluxes compared to the other models.Similar underpredictions of the optical/UV fluxes from direct extrapolation of X-ray spectral models were reported for other ULXs such as Holmberg IX X-1 (Dudik et al. 2016) and NGC 6946 X-1 (Abolmasov et al. 2008).

Donor star or irradiated disc?
As we have seen the only model capable of explaining the broadband data is the diskir, in agreement with earlier works (Berghea & Dudik 2012;Grisé et al. 2012).Instead, the phenomenological and sirf struggle to describe the optical data.We thus considered an alternative description of the data, in which the optical/near UV emission may have additional contribution from the donor star (in Section 5.4 we discuss whether such interpretation holds physically).We approximated its spectrum using a blackbody (bbody in XSPEC) alongside the emission from the accretion flow.We further assumed the putative star parameters remain constant between the high and low states.
Table 4 shows the resulting best-fit parameter, including the constraints on the star temperature  * and radius  * , while Fig. 3 shows the best-fit models and residuals.Both the phenomenological and the diskir would favour  * ∼7 ⊙ and  * ∼ 30, 000K.This would imply an O-type star with a  ∼ 78, 000 ⊙ and a mass of >20 ⊙ (Ekström et al. 2012).On the other hand, the sirf requires a ∼26 ⊙ and  * ∼ 8, 000K star, which would correspond to an F0-A star with a mass of <9 ⊙ (Ekström et al. 2012).We discuss the implication of these results in more detail in Section 5.4.
The blackbody provides a significant fit improvement for both the phenomenological and the sirf (for instance a Δ 2 = 84 for 2 degrees of freedom for the sirf model).However, the diskir again provides the best overall fit ( 2 /dof = 1613.6/1262),without much improvement with respect to the starless model, as the optical fluxes were already well described by the disc emission alone.Similar result was found by Berghea & Dudik (2012) in NGC 6946 X-1.The sirf provides only slightly worse fit in terms of  2 (1630.2/1266)but as can be seen from Figure 3 this model offers a poor description of the high-energy (> 10 keV) tail.This may be expected as the model does not include Comptonization.However, if we consider the trade off between the likelihood and the number of parameters, in terms of Bayesian Information Criterion (BIC; Schwarz 1978) it could be argued that the sirf provides a better description of the data, as it offers the lowest BIC (ΔBIC ≃ -12 compared with the diskir).
The phenomenological provides a marginal description of the data ( 2  ≈ 1.4), offering a poor description of the high-energy (> 10 keV) tail and the optical data.We have verified that the high-energy residuals persists even if we allow Γ to vary between the high and low states ( 2  =1.42), indicating that the phenomenological model cannot account simultaneously for both the high-energy and optical emission.In terms of overall fit, therefore the best representation of the line-of-sight SED is given by the diskir because it can account simultaneously for both the optical and the high-energy tail.In Section 4.2 we put further constraints on the line line-of-sight SED, and provide supporting evidence for the UV extrapolation provided by the diskir.
Crucially, these models make different predictions about the UV flux.This is clearly reflected in Figures 2 and 3.The UV luminosity (defined as the luminosity in the 3 eV-0.1 keV band) is also reported in Table 4.The phenomenological predicts the lowest amount, differing in more than an order of magnitude from the prediction made by the sirf.We are now ready to test these predictions against the emission-line nebula.

Cloudy MODELLING: WHAT THE NEBULA 'SEES'
Having constrained the broadband SED of NGC 1313 X-1, we now turn to examine the effects on its environment through two different sight-lines: sideways, by studying the surrounding, extended nebula reported in Gúrpide et al. (2022) (Section 4.1) and along the line of sight, by studying the line-of-sight integrated nebular spectrum (Section 4.2).In order to select an appropriate ionising SED for our photo-ionization modelling, we noted the numerical calculations presented by Chiang & Rappaport (1996), who studied nebulae photoionized by variable supersoft X-ray sources.Their calculations considered SEDs and densities to a good degree physically relevant for ULX nebulae.These authors showed that so long as the duty the cycle of the source is much shorter than the recombination timescales, the nebula will effectively 'see' or react to the equivalent of a time-averaged SED of the source.In diffuse nebulae the longest recombination timescale is that of hydrogen (Osterbrock & Ferland 2006), which based on typical values from our nebular modelling we find of the order of  (H + ) ∼2×10 5 yr.Shorter timescales are typically (HeII ++ ) ∼2,000 yr and (O +++ ) ∼ 900 yr, therefore much longer than the variability timescales of the source.Assuming the Table 4. Best-fit parameters to the multi-band SED of NGC 1313 X-1.All models were fitted to the multi-band data except for the phenomenological without the blackbody, which was fitted to the X-ray data only.

Notes.
Uncertainties are given at the 1 level.All luminosities are corrected for absorption.
Note the change in number of degrees of freedom here is due to the addition of the HST data in the fit which includes the bbody.Swift-XRT variability is representative of the overall variability of NGC 1313 X-1 (also supported by XMM-Newton observations; Gúrpide et al. 2021a), the temporal average (black dashed line in Fig. 1) suggest the SEDs from the low state are a good approximation for the ionising SED, as the excursions to the high state are too rapid for any meaningful impact on the overall ionising radiation.We return to this point in Section 5.5.We therefore created a set of Cloudy (version C22.02; Ferland et al. 2017) models using our best-fit (extinction-corrected) low-state SEDs from Table 4, totalling 6 different SEDs (3 of them including contribution from the putative companion star).We tested two sets of metallicities for the gas,  = 0.15 ⊙ , 0.3 ⊙5 , corresponding to 12 + log(O/H)= 7.86 and 12 + log(O/H)= 8.17, representative of the nebula around NGC 1313 X-1 (Gúrpide et al. 2022), assumed a filling factor of 1 and open geometry.The center of the nebula was taken as the ULX position in the data cube presented in Gúrpide et al. (2022).
We carried out two set of calculations: one where we assumed the ULX was the only source of ionisation and another where we included an ionising stellar-background from the stars in the field.In order to include a realistic stellar ionising background, we obtained a composite stellar spectrum from the Binary Population and Spectral Synthesis (BPASS) v2.1 (Eldridge et al. 2017) including masses up to 100 M ⊙ .We adopted a metallicity matching that of the gas and an age of the stellar population matching that of the nearby stars ( = 10 7.5 yr; Yang et al. 2011).We followed Simmonds et al. (2021) and rescaled the spectrum based on the star-forming rate (SFR) of NGC 1313.To do so, we noted Suzuki et al. (2013) measured the SFR surface density in NGC 1313 to be ∼0.01M ⊙ /yr/kpc 2 .Considering the region of interest here (∼0.04 kpc 2 ), we found the local SFR = 0.0004 ⊙ yr −1 .Using the scaling between UV luminosity and SFR from Kennicutt (1998), a bolometric luminosity for the stellar background of ∼9.5×10 39 erg s −1 is suggested.This value is only accurate to the order of magnitude owing to uncertainties related to the spatial distribution of the stars in the field, potential contribution from stars outside the photo-ionized region and their exact distance(s) to the nebula.Therefore we have ran our calculations for three different luminosity values for the stellar-background,  = 2.75 × 10 40 , 9.5×10 39 and 2.75×10 39 erg/s, in order to consider the effects of varying this parameter.As we show below, independent on the exact luminosity, all calculations strongly support the fact that stars contribute to the line of ratios of the nebula.The best results were found for  = 2.75 × 10 40 erg/s and we focus on the results obtained for this luminosity throughout, although we also discuss the less luminous Photon Energy (keV)  As stated in Gúrpide et al. (2022), the Heii 4686 was not detected in cube 2 taken in extended mode.The Cloudy predictions we make below prompted us to examine carefully the presence of this line or at least derive an upper limit on its flux in order to constrain further our models.Because we found the line too faint for a pixel-by-pixel fit, we constructed a flux map by integrating the spaxels around the expected position of the Heii 4686 line based on the systemic redshift of NGC 1313 ( = 0.001568) -from 4688.5 Å to 4698.6 Å -and subtracting the mean value of the nearby continuum.The left panel of Fig. 5 shows the resulting map resampled by a factor 3 to highlight a tenuous feature close to the ULX position.We extracted an average spectrum from this region (shown in blue in the right panel of the  2022), we measured an average  (−) value in the same region of  (−) = 0.179±0.005mag.With this value and using the Calzetti et al. (2000) extinction curve with  v = 4.05 we arrived at a total Heii 4686 luminosity of (7.5±1.4)×10 35erg/s.While we do not use this value directly in our model-data comparison, we discuss it further below.

The Side View
In order to study the extended surrounding nebula discovered in Gúrpide et al. ( 2022), we converted the 1D Cloudy predictions to 2D images and resampled them to the MUSE pixel scale of 0.2".We then blurred the images applying a 2D Gaussian kernel with FWHM matching the datacube's PSF (∼1"; Gúrpide et al. 2022).Next we compared the line ratios from these set of Cloudy-generated images to the real MUSE-derived line ratio images.We restricted the datamodel comparison to the regions where the gas was identified as being EUV/X-ray photoionized, where the ULX contribution dominates (see Figure 7 below) and where the influence of shocks is minimized (Gúrpide et al. 2022).The choice to work on line ratios rather than in fluxes was in order to reduce uncertainties related to distance, extinction and geometry.
The geometry above implicitly assumes that the cloud is approximately co-planar with the ULX in the plane of the sky and that the dimension along the line of sight is much smaller than the dimensions on the plane, such that the line ratios along a given line of sight are approximately constant.We are working on extending our modelling to more complex 3D structures and we will present a more refined treatment of the cloud geometry in a future publication, but note our preliminary results assuming a projected, spherical sector give results qualitatively consistent with those presented here.
For each line ratio ([O iii]5007/H, [O i]6300/H, [N ii]6583/H, [S ii]6716/H and [S ii]6716/[S ii]6730, [Ar iii]7135/Hand [S iii]9069/H) we computed a  2 using the model-and data-pixel values and uncertainties on the flux ratio propagated from the measurement error of each line in each spaxel (estimated from the pixel-by-pixel Gaussian-line fitting presented in Gúrpide et al. 2022).By adding each individual line ratio  2 , we were able to select the model with the overall lowest  2 .We did not use [O iii]4959 nor [N ii]6548 as their fluxes were tied by theoretical constraints (Storey & Zeippen 2000).Given other sources of uncertainty such as the exact density profile, geometry, clumpiness of the cloud, projection effects, other sources in the field, etc. and the fact that some of the pixels in our models will have no values in them (see below), our  2 should be regarded as an heuristic to select the best model, rather than the usual goodness of fit.
We produced models for a range of constant hydrogen number densities log( H ), varying from 0.0 to 1.0 in steps of Δ log( H ) =0.05 (noting that the [S ii] lines indicate  e < 100cm −3 ) 7 and inner radius of the cloud  in = 1, 12.5, 22.3, 40 pc, setting values inside the cavity or outside of the Cloudy calculation to 0. The outer radius was set to 200 pc, roughly matching the extent of the high [O i]6300/H ratio and to avoid the emission from a stellar cluster further south (see the bright blob to the far south of the ULX in Figures B1).Fig. 6 shows an example- 2 contours derived for the diskir_star model (circular markers) for the two metallicity values (green and blue for  = 0.15 ⊙ and  = 0.3 ⊙ respectively), which also illustrates that the effects of varying  in are minimal.The Figure also 7 We ran models with higher densities but they clearly failed to reproduce the extent of the nebula.
shows the resulting  2 contours for the same ULX model when the stellar background is included (star-shaped markers).It is clear that the inclusion of the stellar background improves significantly the fit.Fig. 7 shows a 2D data-model comparison for the best-fit disk_star model for  = 0.15 ⊙ for the ULX and ULX + stellar background runs, to illustrate the region we probed and the type of comparisons we carried out.
Tables C2 and 5 show the resulting overall  2 alongside the maximum observed line ratios for the best-fit models and the data for  = 0.3 ⊙ and  = 0.15 ⊙ , respectively.The models with  = 0.3 ⊙ overpredict all line ratios by a factor of ∼2, regardless of whether the stellar background is included or not (Table C2).We therefore were able to reject this metallicity based on the data.This may be surprising as this value matches more closely the metallicity inferred around NGC 1313 X-1.However, as we show below we clearly find that the models with  = 0.15 ⊙ match more reasonably the observed line ratios (see also Fig. 6).We note that in Gúrpide et al. (2022) we were unable to measure the metallicity in the photoionized region itself as metallicity estimators are mostly calibrated for standard HII regions.Hence it may be possible that the nebula has a lower metallicity content than the neighbouring gas.In any case, although such metallicity is low, it is not unrealistic as it is at the lower end of the values measured in Gúrpide et al. (2022).We therefore focused on the results for  = 0.15 ⊙ .
In the case of no stellar background, both diskir and phenomenological offer comparable level of agreement with the data, and superior to that offered by the sirf.The main difference between the phenomenological and the diskir is the peak [O iii]5007/H ratios, due to the lower UV level of the former model.We note that it is actually possible to reproduce peak [O iii]5007/H ratios in the 5-6 range with the phenomenological model (cf.5-7.5 for the diskir and 6-7.5 for the sirf), therefore compatible with the data.However, due to the lower overall values coupled with the smaller region of enhanced [O iii]5007/H produced by the phenomenological due to its lower UV flux, the spatial smoothing smears these values below 3, which highlights the importance of considering the spatial scale.
As shown in Fig. 7 (upper panels), the extent of the [O iii]5007/H area is roughly matched although it is too compact with respect to the data, while the peak values are in fair agreement with those observed.For the [O i]6300/H ratio, although the extent and morphology are well matched, the peak values are overpredicted by a factor 2. This overprediction occurs for most low-ionisation lines (e.g.[S ii]6716/H) as can be seen in Table 5 and for all models.This can be understood due to a lack of strong soft optical spectra from the background stars, which will excite more readily the Balmer lines and increase them compared to e.g.[O i]6300, which is produced in the outer neutral parts of the nebula via highlypenetrating soft X-rays.The inclusion of the putative companion star does not change this basic conclusion and cannot account for these differences.This prompted us to take into account the contribution from the stellar background.
The inclusion of the stellar background not only significantly improved all  2 (Fig. 6; Table 5 for  = 2.75 × 10 40 erg/s and Tables C1 for the less luminous stellar-background cases), but now the peak line ratios in low ionisation lines are much closer to the observed values, particularly for the less-UV bright diskir and phenomenological models.To show the effects of adding the stellar background more clearly, Fig. 8  Green and blue colors show the results for metallicities Z= 0.15Z ⊙ and Z= 0.30Z ⊙ , respectively.We can see that for fixed density, the effects of varying  in are negligible and that the inclusion of the stellar background significantly improves the results.also widens the area over which [O iii]5007 is excited.This effect is similar to that observed by Berghea et al. (2010); Berghea & Dudik (2012) in Holmberg II X-1 and NGC 6946 X-1, where it was found that the low-energy photons from the companion star create low-ionisation states which are further ionised by the high-energy ULX photons (Figure 8 in Berghea et al. 2010).The effect on the low-ionisation lines is instead the opposite.Both these effects support there is an additional source of ionisation along with the ULX.Despite the widening of the enhanced [O iii]5007/H region introduced by the stellar background, all models fail to match the width of the observed radial profile.In Fig. 9 we show a 1D radial profile along the peak of the excitation region comparing the models and the data to illustrate this.The width of the profile was set to 10 pixels to smooth out variations and avoid gaps due low signal-to-noise ratio pixels.Although strictly speaking we compared the 2D generated images with the data, we find these 1D visualization provide an accurate summary of our results.We can see that [O iii]5007/H is overpredicted in all models while the extent over which [O iii]5007 is produced is underpredicted.We can confirm that part of the reason is due to projection effects.As alluded earlier, here we have assumed that the nebula is thin enough in the line of sight direction such that the line ratios along a given line of sight are approximately constant.Instead, if the nebula has considerable structure in the dimension along the line of sight, for instance in the case of a spherical sector, the line ratios along a given sightline will be averaged over regions with different degree of ionization.This will smooth out the line ratios over a wider region compared to the planar geometry assumed here, hence lowering their peak values and broadening their profiles compared to the line ratio profiles shown in Figure 9.We can preliminary confirm these effects from ongoing work, but leave the treatment of more complicated cloud geometries for future work.
Another effect that could be affecting our results is the treatment of the stellar background emission as a point-like source, instead of an extended component.In reality the optical stellar emission will be approximately uniformly distributed over the nebula, both diluting and extending the [O iii]5007/H region compared to the pointlike treatment afforded by Cloudy.This may also explain the wider region of high [N ii]6583/H in the data compared to the models.
To inspect whether this is the case, we have run another Cloudy simulation with the stellar background alone with a typical log( H [cm −3 ]) = 0.60 to inspect the line ratios it would produce.may be partially affected by this and this has to be borne in mind when interpreting the results.
While the fact that the stellar background can produce rather high [O iii]5007/H ratios may call into question whether the nebula is actually produced by the ULX, we see that the stellar background by itself cannot account for the morphology of the nebula (Fig. 8).In particular, the [O iii]5007/H emission is concentrated around the point-like source, while in the models including the ULX the gas close to the source is too ionized to produce [O iii]5007.Sim-ilarly, the stellar background alone does not produce the extended enhanced [O i]6300/H region, which can only be produced by the large mean free path of the X-rays.Therefore, our modelling shows undoubtedly that a high-energy source is needed to explain the nebular morphology.
In terms of  2 , the preferred model is the phenomenological, as it provides the best-overall  2 .From Fig. 9 we can see that both the diskir and phenomenological offer comparable level of agreement with the data, while the sirf overpredicts more severely the [O i]6300/H and [S ii]6716/H ratios.Despite the differences and the clear more-complex profiles in the data due to the effects mentioned above, the overall structure in most lines is well reproduced.The exception to this is the [S iii]9069/H (gray line in the plots).While all models peak at around 120 pc, the observed peak is seen at just 50 pc from the ULX.We suspect this line might be more strongly affected by the sky subtraction and/or shocks, which may explain this discrepancy.As discussed above, the remaining differences may be attributed to projection effects, the treatment of the stellar background as a point-like source and slight differences in abundances and density profile of the gas.Nevertheless, the peak values in the first two models match those observed in the data (see also Table 5).
As stated earlier, we have run an additionally set of simulations by lowering the stellar contribution to study the sensitivity of our results to this component.The results are shown in Table C1 and show that the fits worsen in all instances.In particular, we can see that the models now again overpredict most low-ionisation lines by a factor 1.5-2.We conclude that the bright ( = 2.75×10 40 erg) is a better match to the data, with the caveats outlined above.Nevertheless, the phenomenological and diskir continue to be the preferred models regardless of the exact treatment of the background stars.On this basis, we consider the UV flux predicted by the sirf to be too high to match the nebular emission.
To test our final conclusion, based on our best-fit models, we calculated the expected Heii 4686 flux by multiplying each of the models Heii 4686/H ratios by the observed H extracted from the datacube (Gúrpide et al. 2022).We then extracted the fluxes per spaxel expected from the same blue region shown in Fig. 5.The lower panel of Fig. 9 shows the range of values expected for each model (averaged over the region), the average  (Heii 4686) value measured in the data (Figure 5) and the estimated 3 detection limit.To estimate the latter, we inspected the individual spaxels from the blue region in Fig. 5 and found the lowest flux at which the 3 error on the flux was consistent with 0. This value was ∼10 −18 erg s −1 cm −2 and is shown in the Figure as a black dashed line.We further managed to detect the line in 13 individual spaxels although the relative 3 errors are quite high (70-95%).These detections are also shown in the same Fig. 9 in the lower panel (orange colored histogram).
From the histograms in Fig. 9, we can see that the sirf model produces  (Heii 4686) values that are too high with respect the observed values.Instead, the best agreement is again provided by the phenomenological model, because most of the spaxels have values below the detection threshold, as expected based on the lack of strong detections, and the overall histogram shows the best agreement with the data.Therefore, the marginal detection of the Heii 4686 again reinforces the idea that the predicted UV flux by both the sirf and diskir is overestimated.

The Front View
In Gúrpide et al. (2022) we showed there is a lack EUV/X-ray photoionization signatures around NGC 1313 X-1 in other directions, most likely due to a lower-ISM density in those areas.However, as can be appreciated from the BPT diagrams presented in Gúrpide et al. (2022) (their Figure 11) the [O iii]5007/H ratio is slightly enhanced (∼3.2) at the position of NGC 1313 X-1.The enhanced [O iii]5007/H ratio can also be observed in Figure 10 (left panel), where we show the extinction-corrected spectrum extracted in Section 3.1 in order to measure the extinction towards NGC 1313 X-1.More specifically, from this spectrum we measured [O iii]5007/H = 2.93±1.2,([S ii]6716 + [S ii]6731)/H = 0.65±0.01 and [O i] 6300/H = 0.176±0.006,which again are rather unusual for Hii regions.We recall this is the average spectrum per pixel from a circular region of 1" in radius around the source and therefore can be considered to be from a spatial scale of a single MUSE pixel (0.2" or ∼4 pc at 4.25 Mpc).
While these ratios are not as extreme as in its vicinity (Section 4.1), this is to be expected if we are observing the integrated emission of the photo-ionized nebula along the line of sight.Particularly, should the supercritical funnel in NGC 1313 X-1 be orientated towards us, then it may be reasonable to consider it as a potential source of ionisation.The lack of clear diagnostics, such as the extent of the [O iii] 5007/H and [O i] 6300H regions, makes attributing these enhanced ratios to photo-ionization by the ULX more uncertain (for instance the high ([S ii]6716 + [S ii]6731)/H ratio could be due to shocks, although we have also seen these ratios are also produced by EUV/X-ray photo-ionization; Figure 9).Therefore as a first step we considered whether we could reproduce these ratios as photoionization by stellar continua.
To this end, we attempted to reproduce the fluxes from the lines in the spectrum above -indicated in the Figure with black tick underneath -using photo-ionization models.As in Section 3.1, in order to extract line fluxes we fitted the lines using a Gaussian and a constant for the local continuum.Extinction-corrected fluxes (averaged over the extraction region) for all lines of interest are reported in Table 6.The flux of Heii 4686 was derived using the same region but from cube 2. The 3 negative error was consistent with zero, therefore we considered this measurement an upper limit.Its extinction-corrected flux and the 1 uncertainty are also reported in Table 6.These results were consistent with those obtained using instead a smaller aperture of 0.5 ′′ in radius.
We then compared the observed line fluxes from this spectrum to the fluxes obtained by integrating the Cloudy nebulae along the line of sight using again  2 statistics.The upper limit on Heii 4686 was  Notes. To calculate the reduced  2 , the number of degrees of freedom is defined as the number of pixels covering the photo-ionized region (11474) minus 2 variables ( H and  in ).These  2 need to be understood as an heuristic to rank the models, rather than a goodness-of-fit.taken into account in the modelling following Hoof (1997): where  obs and  model refer to the observed and predicted (Heii 4686) fluxes and   obs is the (1) uncertainty on  obs .
We assumed we only see nebular emission due to gas located be-tween the ULX and ourselves (i.e.no contribution from an equallyextended nebula behind the ULX).Assuming we only observe the gas in front of the ULX may be reasonable considering much of the emitting nebular gas behind the ULX will be absorbed by the system itself.However, such fraction may be small considering we are observing gas averaged (or integrated) over an angular (physical) diameter of 2 ′′ (∼40 pc).Assuming a symmetric case, where the gas is equally distributed in front and behind the ULX, is not expected to affect the results significantly as the effect would be to simply rescale the line fluxes for a given density/inner radius by a factor 2, without altering the predicted line ratios.We have rerun the calculation under this assumption and verified that while there are small numerical differences between the two approaches, our overall conclusions below do not change.We present the results for the former assumption (termed asymmetric case) and discuss the calculations for the latter assumption (symmetric case hereafter) when relevant.
We performed the data-model comparison following two approaches and found very good agreement between them (Δ 2 < 6 for all models), so we only present the results from the latter approach.In the first method, we considered the averaged fluxes over the spatial region and assumed we are observing a column of  ×  × out , where  is the physical size of a MUSE pixel at  = 4.25 Mpc and  out is the (unknown) outer edge of the nebula along the line of sight.We found best results when considering all models to be ionization-bounded.Therefore the nebula outer radius  out was set to the maximum of the Cloudy calculation ( out = 200 pc) or lower in cases of higher density where the hydrogen ionization front was reached before.In the second approach, instead of considering the average spectrum over the 1" aperture, we have instead used the integrated one.We then have compared the integrated line fluxes to those obtained in Cloudy by using the aperture command to integrate the simulated nebula over the same spatial region.
We noticed some of the models  2 minima were at the limit of our density calculation (log( H [cm −3 ]) = 1), so we extended the calculations to densities up to log( H [cm −3 ]) = 1.6 and  in = 60 pc to ensure we found the absolute minima for each model.
As stated above, we first verified whether we could explain the line fluxes as photo-ionized by a population of stars by running models with the stellar background alone.In particular, we have run models for stellar backgrounds with  = 2.75 × 10 38,39,40 erg/s for the range of densities and radii quoted above.We have found  2 /dof upwards of 470 in all instances (also for the symmetric case), with unusually lowly extended (≲ 30 pc) nebulae due to the hydrogen ionization front being reached very close to the source as a result of the gas being optically thick to the soft radiation.Moreover, sulfur or [O i] 6300 lines were predicted to be ≳5 and ≳ 15 lower than observed, respectively.This agrees with the fact that these pixels were classified as 'AGN' in the BPT diagram presented in Gúrpide et al. (2022).Therefore, we considered whether we could put further constraints on the SED of the ULX along the line of sight by modelling the nebular emission in this direction as photo-ionization by the ULX instead.
Table 6 shows the obtained data/model line ratios and resulting  2 for all models, including those for which the stellar contribution discussed in Section 4.1 is added to the ULX SEDs. Figure 10 (right panel) shows the  2 contours obtained for the three models.As opposed to the extended nebula (Section 4.1), where the inclusion of the stellar background clearly improves the results, here most models worsen or show little improvement when the background stars are added.It is also clear that best results are obtained for the non-background or moderate ( = 2.75 × 10 39 erg/s) background cases (Table 6 and Figure 10, right panel).Similar trend was found for the symmetric case.This is reasonable considering here we are modelling a single sight line, where the major (or only) contributor is likely to be the ULX.The best match is obtained for the diskir model alone (or with moderate background stars) and suggest the nebula is denser ( H ∼ 9 cm −3 ) and ∼112 pc long along the line of sight.Expectedly, the additional gas implicitly present in the symmetric case instead lowers the required cloud density, to  H ∼ 5.6 cm −3 .Note that although the reduced  2 is high, almost all lines are predicted within a factor 2 or less with the diskir model.Thus the ULX clearly provides a much better match to the line fluxes than the stellar continua above.The worst agreement is found for the sulfur lines.As alluded in the previous Section, we suspect [S iii]9068 is likely affected by the sky subtraction.The discrepancy between the model predictions and the other sulfur lines is less clear cut, but may be attributed either to some contribution from shocks and/or to a higher abundance of sulfur.
Along with the predicted line fluxes and  2 , in Table 6 we also report the integrated hydrogen absorption column along the line of sight, together with that derived from spectral fitting (Section 3).Most derived  H values from the photo-ionization modelling, although not unrealistic, are higher than those derived from X-ray spectral fitting after subtraction of the galactic contribution to the total  H .In principle, we would expect  H derived from spectral fitting to be higher than that of the nebula due to additional contribution from the system itself.One possibility is that indeed there is some contribution to the nebular lines from behind the ULX.For the symmetric cases,  H is reduced approximately by ∼ (0.3 − 0.4) × 10 21 cm −3 , which would bring them to a level below that derived from spectral fitting, arguably in more reasonable agreement with the X-ray data.It is also likely that there are additional sources of uncertainty due to the exact underlying model and the fact that we have not included the exact (unknown) abundances into account when deriving  H from spectral fitting.Therefore these values should be taken with caution.Nevertheless, our analysis suggest most contribution to the absorption column may be due to the photo-ionized nebula itself.
While our photo-ionization modelling here is more uncertain due to additional contribution from possible shocks and the unknown extent of the nebula, we can confidently rule out the sirf model as a good line-of-sight SED, as it cannot account for the low-ionisation lines (as they are all underestimated by a factor 3 or more) while producing too strong Heii 4686.At the same time, it suggest our best-fit line of sight SED (the diskir; Section 3) is not unrealistic based on the nebular line fluxes observed along the line of sight.Therefore, while we cannot rule some contribution from shocks, we can at least ascertain that a high-energy source (namely the ULX) is a good match to the line fluxes if these are to be explained by a photo-ionization and that these are best explained by the diskir model.

DISCUSSION
Through multi-band spectroscopy and modelling of the nebular emission we have been able to confirm our previous assertion (Gúrpide et al. 2022) that NGC 1313 X-1 powers an extended ∼200 pc EUV/Xray photo-ionized region, with additional contribution from the stars in the field.Using state-resolved multi-band spectroscopic data, we have constrained the SED of the ULX along the line of sight.The best-fit model, capable of explaining both the high-energy tail (>10 keV) and the UV/optical data simultaneously is the diskir model (regardless of whether we add any contribution from the putative stellar counterpart; Section 3).This result agrees with earlier works that have found that such model, despite not describing the accretion flow geometry envisioned for a supercritically accreting compact object (e.g.Lipunova 1999;Poutanen et al. 2007;Abolmasov et al. 2009), fits ULX broadband data satisfactorily (e.g.Kaaret & Corbel 2009;Grisé et al. 2012;Berghea & Dudik 2012;Tao et al. 2012).We do not however take this as evidence for NGC 1313 X-1 being powered by a standard self-irradiated disc -an unlikely scenario owing to the unusual X-ray spectral shape (Bachetti et al. 2013), the Notes.Uncertainties at the 1 level. 3 upper limit along with the 1 uncertainty.
Degrees of freedom are defined as 10 (lines) -2 variables (log( H ) and  in ).
Extra-galactic neutral absorption column derived from spectral fitting (quoted from the diskir model; Table 4) i.e. the Galactic contribution along the line of sight ( G H = 7.07×10 20 cm −2 ; HI4PI Collaboration et al. 2016) has been subtracted.
presence of relativistic outflows (Pinto et al. 2016) and the ∼400 pc shock-ionized bubble surrounding the source (Gúrpide et al. 2022).Instead, we consider the diskir as a physically reasonable extrapolation to the unaccessible EUV.Such extrapolation to the UValthough uncertain -is supported by our modelling of the nebula along the line of sight (Section 4.2).We have further shown that direct extrapolation of X-ray only models (the phenomenological model; Section 3) fail to account for the full SED when the optical data is considered, an issue highlighted also in previous works (Abolmasov et al. 2008;Dudik et al. 2016).
Through photo-ionization modelling of the extended emission-line nebula (Section 4.1) we have instead attempted to constrain the SED seen by the nebula through a different sight line.The fact that the extended [O iii]5007 and [O i]6300 emission do not overlap (Fig- ures 7 and 8) is a strong indication that the nebula sees the accretion flow in NGC 1313 X-1 sideways.Therefore NGC 1313 X-1 and its EUV/X-ray excited nebula offers us an opportunity to study super-Eddington accretion flows effectively from two different sightlines.
Here we have shown that the nebular lines are best described with a model with a lower UV flux that than constrained along the line of sight.In particular, the best match to the nebular lines is provided by the phenomenological model, whose UV flux is about a factor 4 lower compared to the diskbb (Table 4).This model not only provides the best match to the nebular lines (Table 5) but also accounts for the lack of strong nebular Heii 4686 detection (Fig. 9).One may argue that similar results were found by Berghea et al. (2010); Berghea & Dudik (2012) in Holmberg II X-1 and NGC 6946 X-1, wherein it was found that the models with the lowest UV levels were a better match to the nebular lines (their Figures 4 and 3 respectively or their PLMCD and MBC models respectively)-although we caution that these works did not use spatially-resolved data, which is needed if the degree of beaming is to be determined.
Below we elaborate on how these findings, that is, the discrepancy between the best-fit line-of-sight SED (with  UV ∼ 2 × 10 39 erg/s) and the best-fit nebular model ( UV ∼ 0.5 × 10 39 erg/s) allows us to put constraints on the degree of anisotropy of the UV emission in the ULX NGC 1313 X-1.

Quasi-isotropic UV emission
Our results may suggest the extended nebula is not seeing the same SED as that derived from the line-of-sight.That is, our observations may be interpreted as evidence for a mild degree of anisotropy in the emission of NGC 1313 X-1.Because of the nature of the emission lines used in this work, with the highest ionisation potentials falling in the UV band (Figures 2 and 3), our observations constrain the degree of anisotropy mostly in the EUV, and any extrapolation to the soft/hard X-rays remains more speculative.This can be observed, for instance, by the fact that the Heii 4686 line is insensitive to the X-rays (≳0.1 keV), as the three models predict vastly different Heii 4686 fluxes (Fig. 9) despite all having similar X-ray luminosities constrained by the data (Table 4).Observations probing higher excitation lines found in the IR (Berghea et al. 2010;Berghea & Dudik 2012) will allow to put tighter constrains and extend our measurements to higher energies.
Although below we provide quantitative calculations for the degree of beaming based on our results, we would like to highlight that there are inevitably additional sources of uncertainty we cannot account for, potential contribution of shocks, other SED extrapolations we have not tested for and the treatment of the stellar background (we discuss these in more detail below in Section 5.5).Nevertheless, while the quantitative details may be uncertain, we believe our results can be confidently interpreted as a lack of strong anisotropy in the UV emission.
The differences in UV luminosity between the best-fit line-of-sight SED (diskir) and the nebular one (phenomenological) may be used to constrain the beaming factor  proposed by King et al. (2001).King et al. (2001) defines the beaming factor as: where  is the true emitted radiative luminosity and  sph is the observed luminosity under the assumption of isotropic emission.
Crucially, this factor must also depend on the inclination of the system and energy.The dependence of  with the inclination complicates estimating this value quantitatively from an observational point of view, particularly due to the uncertain inclination of the ULX with respect to our line of sight and the nebula.Nevertheless, we may approximate its value as: where for  sph we assume that NGC 1313 X-1 is observed close to face on ( = 0 • ) and where   nebula is the luminosity observed by the nebula at an unknown inclination angle .We may get a crude estimate of  from the diskir normalization, which is ∝ cos().We have reran our analysis in Section 4.1 by 'inclining' the diskir_star to inclination angles  = 45 • , 60 • , 80 • , finding the best-fit ,  in in each case.All calculations were carried out with the stellar background stars and assuming the flux of the putative companion is isotropic.We have found the best match to the nebular lines is given for  = 45 • with log( H [cm −3 ]) = 0.45,  in = 1 pc and  2 /dof = 89 (cf. 2 /dof = 92 for the diskir_star).However, this inclined version not only provides a worse fit than the phenomenological model, but also still produces Heii 4686 in comparable numbers to the nominal diskir_star.This suggests the differences between line-of-sight and the ionising SED are more complex than a simple scaling factor.
Alternatively, we may find  by considering what value of  is needed to reduce the UV luminosity of the diskir to a level comparable to that of the phenomenological (i.e. a factor 4 dimmer).This would imply  = 80 • , although this model gives worse fits to the nebular emission than the nominal diskir.These values are obviously uncertain due to the uncertainties related in using the diskir to describe a super-Eddington accretion flow, but may suggest the nebula sees an inclined  = 45 • − 80 • version of the line-of-sight SED.Nevertheless, we have already argued that the non-overlapping [O iii]5007/H and [O i]6300/H regions strongly suggest the nebula sees the emission sideways.We additionally note that the sirf instead predicts a nearly isotropic UV flux and no amount of inclination can yield the necessary reduction in UV flux.
Therefore considering the differences in UV luminosity between the line-of-sight SED (diskir) and that inferred from the extended nebula (phenomenological), we constrain the beaming factor in the EUV of NGC 1313 X-1 to about 0.3-0.15 and suggest the nebula sees an inclined  = 45 • − 80 • version of the line-of-sight SED.Such estimates are in agreement with measurements of the photoionized nebula around Holmberg II X-1 (Kaaret et al. 2004) who found beaming factors  >> 0.1 (their Table 1) -with the caveat of the extrapolation of the X-ray spectrum (see the Introduction).Such estimates are also consistent with detailed analytical calculations by Abolmasov et al. (2009).If such results can be extrapolated to the Xray band, then we may be able to rule out the strong beaming invoked to explain the emission from some PULXs (King & Lasota 2020) and would support the constraints derived on  from modelling of the observed PULX pulse-fractions (Mushtukov & Portegies Zwart 2023;Mushtukov et al. 2021).
Alternatively, the beaming factor may indeed be more pronounced in the X-ray band as alluded in the Introduction, as these photons emanate from within the wind funnel from the supercritical disk.Instead, the wind photosphere, which is expected to dominate the UV emission, is expected to emit quasi-isotropically (see e.g.Shakura & Sunyaev 1973;Weng & Feng 2018).Hence our results would be consistent with this picture.In this regard, our results also seem broadly consistent with the general relativistic radiation-magnetohydrodynamic simulations presented by Narayan et al. (2017), although it is hard to derive a quantitative comparison.Broadly speaking, their postprocessed spectra (e.g.their Figure 13) show the UV is reduced by about a factor ∼6 when going from  = 10 • to 60 • .This would be broadly consistent with our reasoning based on the diskir model and the increase in inclination needed to match a reduction in flux of about 4.

ULXs in the UV
Regardless of whether we consider the intrinsic UV flux in NGC 1313 X-1 as the line-of-sight value or that inferred from the extended nebula, both estimates are significantly lower than those measured in NGC 6946 X-1 (Abolmasov et al. 2008;Kaaret et al. 2010).The measurement in the F140LP reported by Kaaret et al. (2010) is shown in Fig. 4 for comparison and is about a factor 5 than predicted by any of our models.Similarly, we have integrated the Heii 4686 line predicted by the best-nebular model (phenomenological with the background stars) over the whole nebula (assuming isotropy), and found  (Heii4686) ∼9×10 35 erg/s, slightly above our observed value of (3.4±0.6)×10 35erg/s (Section 4).The value derived from the modelling may be considered an upper limit as we have considered the same Heii 4686H ratios everywhere around the source.Still, such value is significantly lower than the measured  (Heii4686) = (2±0.2)×10 37erg/s in the MF16 nebula surrounding NGC 6946 X-1 reported by Abolmasov et al. (2008).If these UV differences were due to an inclination effect, then we should have seen evidence for strong UV (comparable to that of NGC 6946 X-1) either along the line of sight or in the nebular lines.Such differences suggest the UV emission in NGC 6946 X-1 and NGC 1313 X-1 is intrinsically different, suggesting we are probing differences in the mass-transfer rate.
Because the UV emission is thought to be linked to the wind photosphere (∝  0 −3/4 ; Poutanen et al. 2007, where  0 is the mass-transfer rate at the companion in Eddington units) our analysis suggest a lower mass-transfer rate in NGC 1313 X-1 compared to NGC 6946 X-1, despite the stronger X-ray luminosity of the former (Gúrpide et al. 2021a).We suggest the inclination of these two systems might be comparable, but NGC 6946 X-1 might possess a narrower funnel due to higher accretion-rate, creating a strong soft X-ray/EUV source and mimicking a highly inclined source.In NGC 1313 X-1 instead we peer down the funnel most of the time owing to the wider opening angle of the funnel, except in the extremely soft and unusual 'obscured state' (Gúrpide et al. 2021a), where we showed the source becomes even softer than NGC 6946 X-1.A brighter UV and softer X-ray spectrum in NGC 6946 X-1 due to a higher mass-transfer rates would be fully consistent with predictions from R(M)HD simulations (Kawashima et al. 2012;Narayan et al. 2017) and arguments made by Abolmasov et al. (2007) based on the observed nebular lines in a sample of ULXs.
In Table 7 we have collated literature results regarding the properties of high-excitation nebulae surrounding ULXs for which the X-ray spectral regime is known, which we have taken from Sutton et al. (2013); Urquhart & Soria (2016); Gúrpide et al. (2021a).Notice the apparent differences in nebular emission between hard ULXs such as NGC 1313 X-1 and M81 X-6 and soft ULXs such as NGC 5408 X-1 (Kaaret & Corbel 2009), Holmberg II X-1 (Kaaret et al. 2004) or NGC 6946 X-1 (Abolmasov et al. 2008).These differences are not trivial as we have already stressed the fact that the Heii 4686 is insensitive to the X-rays (Fig. 9 lower panel) as its ionization cross-section falls roughly as  −3 (Osterbrock & Ferland 2006).The Heii 4686 is obviously sensitive to the density or distribution of material around the ULX.For this reason in the Table we also report the Heii 4686/H ratio, which should be less sensitive to density differences.As can be seen, the ratios are generally lower in hard ULXs, indicating that the differences in Heii 4686 luminosity are not due to a density difference.Thus the differences in nebular Heii 4686 around hard and soft ULXs are a telltale that hard and soft ULXs are not only distinct in their X-ray spectral properties, but also in their UV/EUV.Because strong EUV is linked to the mass-transfer rate according to RMHD simulations (Narayan et al. 2017) and analytical estimates (Poutanen et al. 2007), bright Heii 4686 around soft ULXs signals these systems must possess higher mass-transfer rates compared to their hard counterparts.Instead, hard ULXs, due to their faint EUV, can produce strong [O iii]5007/H ratios but dim or no detectable Heii 4686 (e.g. as it is the case in NGC 1313 X-1, Holmberg IX X-1 and NGC 1313 X-2).Our findings are consistent with earlier assertions made by Abolmasov et al. (2007) and suggest the mass-transfer rate may be more relevant parameter distinguishing hard and soft ULXs, as opposed to their inclination (Sutton et al. 2013).
In Table 7 we have also collated mechanical powers  mec inferred from observations of optical/radio bubbles surrounding ULXs.Because we expect the outflow power to increase with  0 (e.g.Kitaki et al. 2021), we should expect soft ULXs to possess higher  mec .From the current limited sample, it does not appear that soft ULXs show higher  mec , but certainly a systematic study is needed here, which is beyond the scope of this paper.It is also unclear whether some of these nebulae are comparable: for instance, the bubble around Holmberg IX X-1 shows a nearly spherical morphology (Abolmasov & Moiseev 2008), while Holmberg II X-1 or M51 ULX-1 show instead bipolar bubbles associated with collimated jets (Cseh et al. 2014;Urquhart et al. 2018).Soft ULXs seem also to be less likely associated with shock signatures (e.g. the case of Holmberg II X-1 or NGC 5408 X-1 Kaaret et al. 2004;Cseh et al. 2012), which may suggest hard ULXs are more likely to clear out the material around them, hindering the detectability of the Heii 4686 line.Other factors such as the mass, or even the nature of the accretor are also likely to play a role in explaining such differences.In particular, hard ULXs have been systematically shown to be more likely to host NSs (Pintore et al. 2017;Walton et al. 2018;Gúrpide et al. 2021a;Amato et al. 2023).Whether the nature of the accretor could also explain such differences remains to be seen, but from the limited sample it would seem that differences in the X-ray spectral are translating into differences in the interaction with the environment.

The nebular Heii 𝜆4686 problem
Whether ULXs can produce Heii emission in enough numbers to account for the Heii 4686 line observed in the integrated spectra of metal-poor galaxies (the Heii 4686 problem) has been recently examined in Simmonds et al. (2021) and Kovlakas et al. (2022), reaching contradictory results.Simmonds et al. (2021) found that a the multi-band diskir model presented by Berghea & Dudik (2012) could produce Heii 4686 in enough numbers to explain it.Kovlakas et al. (2022) on the other hand, built empirical models based analytical descriptions of super-Eddington accretion discs (Shakura & Sunyaev 1973;Lipunova 1999;Poutanen et al. 2007) and found ULXs do not produce enough ionising UV photons to explain the Heii 4686 line.
The main uncertainty on these works was the lack of knowledge about the UV emission (see Figure 6 in Kovlakas et al. (2022)), to which the Heii 4686 line is most sensitive.Based on our discussion above and our analysis, we suggest that there should be a dichotomy between hard and soft ULXs in terms of their capacity to excite Heii 4686, with hard ULXs ruled out as potential candidates to explain the nebular Heii 4686 in metal-poor galaxies.Finally, we also note that the relatively isotropic EUV found here would render uncertainties related to beaming and its dependency with  in the study of the nebular Heii 4686 problem nearly unimportant (e.g.Kovlakas et al. 2022).We aim to pursue similar studies in other ULXs with different X-ray spectral hardness to reliable confirm these results.

The origin of the optical/near-UV light
In Section 3 we have included a blackbody component to model the optical/near-UV fluxes from the HST.While this component is not required by the diskir, it is strongly required by the sirf and phenomenological models to fit the HST data.We initially presented this component as a proxy for the contribution from the companion star.However, based on its luminosity and temperature, we now investigate whether a stellar origin is physically plausible and consider alternative explanations.
We have seen both the phenomenological and the diskir give similar parameters for this blackbody.The temperature and luminosity of this component imply an O-type star, requiring a star >20 M ⊙ .This would be consistent with the donor OB-types inferred from direct optical modelling of the ULX spectra (Tao et al. 2011).However, here we cannot rule out this type of star as it would be at odds with constraints from the population and age of the nearby stars presented in Yang et al. (2011), which suggest masses under 10 M ⊙ instead.Therefore for NGC 1313 X-1, we can rule out the optical/near-UV fluxes being dominated by the donor star, which cast doubts on the OB-type stars inferred in other ULXs (e.g.Grisé et al. 2012).This is consistent with the optical short-term variability observed by Yang et al. (2011) (and confirmed here Section 2) and the strong optical variability observed during the nascent ULX in M83 (Soria et al. 2012).Yang et al. (2011) and Tao et al. (2011) presented additional diagnostics that can help shed light on the nature of the emission.Based on our reanalysis of the HST data and the state-resolved optical and X-ray data we updated the Johnson-Cousin Vega extinction-corrected magnitude  0 = 23.28 ± 0.07 and the optical to X-ray ratio diagnostic from van Paradĳs (1981)  =  0 + 2.5 log( X ) = 21.8 ± 0.7, reported in Yang et al. (2011).As stated in Yang et al. (2011), such values are typical for low-mass X-ray binaries and suggest the emission is dominated by the disc itself (or reprocessed emission from it) rather than the companion star (Kaaret et al. 2017).Indeed, Tao et al. (2011) found through careful analysis of the HST data of a handful of ULXs counterparts that the optical emission is not consistent with any stellar type.
An alternative explanation is that this component represents the emission from the wind photosphere.The temperature of this component ( ∼ 3 × 10 4 K) is comparable to that measured in SS433 ( = (7 ± 2) × 10 4 ; Dolan et al. 1997) and in NGC 6946 X-1 ( = 3.1 × 10 4 K; Kaaret et al. 2010) both using similar photometric filers and a single blackbody model.We note however that the value on SS433 is highly uncertain as relies on a likely overestimated level of extinction, as noted by the detailed X-SHOOTER analysis on SS433 (Waisberg et al. 2019).
Nevertheless, the radii inferred for both SS433 and NGC6946 X-1 are of the order of 10 12 cm, while for NGC 1313 X-1 we measured about half that value.Whether we can associated any predictive power to this component is questionable owing to the uncertainty on the exact modelling (and in fact most likely additional components are present/needed; Kaaret et al. 2010;Soria et al. 2012), however given the photosphere radius is expected to scale with  3/2 0 , this is indeed in line with our previous suggestion of a lower mass-transfer rate in NGC 1313 X-1.Such interpretation is consistent by the lower overall UV luminosity in NGC 1313 X-1 (≲ 10 39 erg/s) compared to NGC 6946 X-1 (Abolmasov et al. 2008) or SS433 (Dolan et al. 1997;Waisberg et al. 2019), which both have a UV luminosities in excess of 10 40 erg/s.
In this regard, while we have shown the level of UV predicted by the sirf is too high to explain the nebular emission (Section 4), it may offer a more accurate description of the optical/near-UV data.This is because when employing the sirf model, it is the (supercritical) disc that dominates the optical/near-UV fluxes, with the putative companion adding a negligible contribution (Fig. 3).Such description of the SED would fit more accurately our expectation of the optical/near-UV fluxes based on the reasoning above.Given the putative star parameters, we infer an F0-A-type star with a mass of ≲ 9M ⊙ (Ekström et al. 2012) which would be consistent with the population of stars around NGC 1313 X-1 (Yang et al. 2011).We further note that this is the stellar-type inferred for the companion star of SS433 (see Goranskĳ 2011, and references therein), which may add supporting evidence for the link between SS433 and ULXs.However, we must stress that this ignores binary evolutionary effects and irradiation by the X-rays from the disc/wind, which will distort the spectrum of the companion star (see discussion in Ambrosi et al. 2022;Sathyaprakash et al. 2022, and references therein).In any case, our analysis suggest that the optical/UV emission in NGC 1313 X-1 is not dominated by the companion star.Similar situation was found in Holmberg II X-1 (Tao et al. 2012) and may suggest the same is true in NGC 5408 X-1 (Grisé et al. 2012).
Finally, we consider the results presented by Vinokurov et al. (2018), who pointed out that the brightest ULXs in the optical (approximately  V < −5.5 where  V is the absolute magnitude in the Johnson  band) show a powerlaw-like spectra, whereas dimmer ULXs instead have blackbody-like optical spectra.Vinokurov et al. (2018); Fabrika et al. (2021) argue that the dimmer appearance of some ULXs is due to a lower contribution of the wind photosphere, linked to the mass-accretion rate.Such interpretation would reinforce our conclusion that indeed the  0 in NGC 1313 X-1 is lower compared to softer ULXs such as e.g.NGC 6946 X-1 or NGC 5408 X-1, which instead have  V <-6 (Tao et al. 2011).Vinokurov et al. (2018) argues that the blackbody-like spectra in the dimmer systems may represent emission from the donor.While the low optical luminosity of NGC 1313 X-1 ( V ∼ −4.9) may suggest the latter scenario applies here as well, we have already seen that the companion star inferred (O-type) from blackbody fits is at odds with the stellar population around NGC 1313 X-1.Ignoring geometrical and irradiation effects (Abolmasov et al. 2009), assuming the wind has the virial velocity at the spherization radius (Poutanen et al. 2007), and that the wind photosphere radiates as a spherical blackbody, one can show that its luminosity is roughly equal to a third of the Eddington luminosity and independent of  0 (Poutanen et al. 2007): Thus the luminosity of the optical blackbody (∼ 2×10 38 erg/s) can be easily explained a by a ∼5 M ⊙ BH.However we note that we could not match its temperature ( ∼3×10 −3 keV) for a reasonable value of  0 .
To illustrate this, consider the temperature of the wind photosphere may be expressed as Poutanen et al. (2007): where   =  wind /( wind +  rad ), with  wind the kinetic luminosity of the wind and  rad the observed radiative luminosity (Poutanen et al. 2007).Assuming   = 0.5 (Pinto et al. 2016;Gúrpide et al. 2022), in order to match the observed blackbody temperature we would need an abnormally high  0 ∼ 1000 for the aforementioned BH mass.Thus, as alluded above, it is unlikely that we can associate this component solely to the wind photosphere and the picture may be more complex due to geometrical, irradiation effects (Abolmasov et al. 2008) and deviations from a blackbody due to scattering (Lipunova 1999).

Caveats
There are several limitations of our study that are worth highlight and need to be borne in mind when interpreting the results.The first is the timescales involved in the nebular emission with respect to the variability of the source, which will unavoidably affect any study of this kind.As alluded in Section 4, the recombination timescales of the nebula are of the order of thousands of years, while the observation baseline of NGC 1313 X-1 is only of dozens of years, if we consider XMM-Newton observations prior to the Swift-XRT.Hence the time-averaged spectrum could be higher or lower depending on the (unaccessible) activity history of NGC 1313 X-1.However, the fact that the diskir provides a reasonable description of the line of sight nebular lines (all lines accurately predicted within a factor 2) -which are obviously subject to similar recombination timescales -suggests that the time-average SED cannot deviate substantially from the present-day estimate.In fact, we can rule out a brighter time-average spectra as the sirf substantially overpredicts the line fluxes both along the line of sight and in the extended nebula (Tables 5  and 6).We have explicitly confirmed that the high-state SEDs did not offer better match to the nebular lines than the low state SEDs.There remains the possibility that the time-average spectrum is dimmer than the diskir but still brighter than the phenomenological, because the level of UV provided by the phenomenological already provides a worse match to the line-of-sight line fluxes than the diskir (Table 6).It may thus be possible to explain both the line of sight and extended nebular line fluxes better by a time-average spectrum whose UV brightness sits between these two spectra, which may suggest the UV is closer to being isotropic than we have estimated.
Another possibility is that some of the line ratios (or fluxes) are slightly affected by shocks.In our modelling of the extended emission (Section 4.1) we do not expect this to be an issue affecting significantly our results.The first reason is that we already showed in Gúrpide et al. (2022) that shocks are mainly concentrated in the edges or rim of the bubble, whereas the inner parts where instead dominated by EUV/X-ray excitation.Secondly, the morphology of the nebula, that is the observed ionization gradient, most remarkable in the [O iii] 5007 and [O i] 6300 lines, can only be produced by EUV/X-ray photo-ionization.Shocks would instead produce rather co-spatial [O i] 6300 and [O iii] 5007 regions.Moreover, we already showed in Gúrpide et al. (2022) that to produce the observed levels of [O iii] 5007 would require shock velocities >300 km/s (see also Berghea et al. 2010), clearly ruled out by the kinematic data.Therefore, if anything, we may expect a slight contribution from shocks in low-ionization lines in the outermost parts of the photoionized region.However, winds can still alter the distribution of the gas, by depleting and compressing it onto a thin shell (e.g.Siwek et al. 2017;Garofali et al. 2023;Gúrpide et al. 2024).To some extent this has been taken into account, since we have marginalized our results over a range of radii and density.In practice however, it is often hard to derive an accurate description of the geometry of the nebulae, due to projection effects.We are also working on extending our modelling to 3D structures to provide potentially a more accurate description of the cloud geometry.However, preliminary results suggest that extending to 3D, at least for the present work, is not likely to have a strong impact on the results.
Shocks may be important however in our photo-ionization modelling of the nebula along the line of sight.As stressed in Section 4.2, the lack of clear diagnostics such as the extent of the regions with enhanced oxygen to Balmer line ratios makes the attribution of the observed line ratios to photo-ionization by the ULX less certain.A more refined modelling, self-consistently accounting for shocks and photo-ionization may provide a more accurate picture.Nevertheless, once again we do not expect this to strongly affect lines such as [O iii] 5007 and Heii 4686 and therefore we can confidently ascertain that the line of sight SED is best described by the diskir or the phenomenological models.Thus we consider our assertion that beaming in the UV must be small to be robust to these caveats.We leave for future work a more detailed treatment self-consistently accounting for photo-ionization and shocks.
Lastly, moving forward it may be advantageous to include the stars in the field in a spatially resolved manner, here and for other ULXs in crowded fields such as Holmberg II X-1 (Pakull & Mirioni 2002).
At present, this cannot be done in Cloudy, but may be explored in the future with three-dimensional photo-ionization codes such as MOCASSIN (Ercolano et al. 2003).

CONCLUSIONS
Coupling multi-band spectroscopy with detailed modelling of the EUV/X-ray excited emission-line nebula surrounding NGC 1313 X-1, we have attempted to constrain the degree of anisotropy of the UV emission in this archetypal ULX.Our results suggest that the UV emission is mildly beamed, by a about a factor ∼4 at most.We have also shown that the optical emission in NGC 1313 X-1 is unlikely to be dominated by the companion star.
We have also discussed the weak detection of Heii 4686 in the nebula surrounding NGC 1313 X-1, which seems to be a common finding around other hard ULXs.Instead, bright nebular Heii 4686 seems to be a common finding around soft ULXs.We suggest differences in mass-transfer rate may explain such dichotomy, since according to analytical calculations and numerical simulations only at high-mass transfer rates a ULX will become an extreme EUV source.This implies that only a subset of the whole ULX population may excite Heii 4686 in high enough numbers to account for the observed Heii 4686 line in metal-poor galaxies.
Moving forward, a better understanding of the ULX SED or observations targeted at reducing the uncertainty on the line-of-sight SED by probing the FUV emission would be of great interest to reduce the uncertainties in our work.While the lines probed here do not allow us to constrain the degree of anisotropy of the soft X-ray emission, probing the high-excitation lines found in the IR with James Webb Space Telescope would enable to constrain the degree of anisotropy of the ULX emission at higher energies, allowing us to test existing super-Eddington accretion theories and improving our understanding of ULXs and the feedback on their environments.

Figure 1 .
Figure 1.Swift-XRT observations showing the variability and spectral states of NGC 1313 X-1.(Left) Swift-XRT lightcurve.The dashed blue line shows the time when the HST/WFC3/UVIS observations were performed (too close in time to be distinguished here).The red and green shaded areas show the equivalent Swift-XRT countrates of the 2003 and 2004 Chandra observations respectively.The uncertainties on the former includes any uncertainties associated with the pile up modelling (see text for details).The black solid and dashed line shows the mean Swift-XRT count-rate and its standard error.(Right) hardness ratio given as the count rate in the 1.5 -10 keV band over the 0.3 -1.5 keV band.The red and green stars mark the high and low state of NGC 1313 X-1 respectively, which were derived taking the mean and standard deviation of the snapshots above and below 0.225 ct/s respectively.

Figure 4 .
Figure 4. BPASS stellar template alongside the ULX models used for the low state of NGC 1313 X-1.Both stellar templates have an age of 10 7.5 yr and were approximately rescaled to the local SFR around NGC 1313 X-1.The shaded bands are as per Fig. 2 and correspond to the bands used to calculate the fluxes in Table4.
Figure) along with a spectrum from the nearby stellar cluster (shown in orange) to compare the presence of the Heii 4686 line.The resulting spectra are shown in Fig. 5 (right panel).The line is clearly detected 6 in the patch next to the ULX at the expected position based on the redshift of NGC 1313 (shown in the Figure by a black dashed vertical tick).Fitting a Gaussian and a constant for the local continuum (shown as a magenta dashed line) we measured an average flux (Heii 4686) = 3.6±0.7×10−19 erg/s/cm 2 .From the  ( − ) map derived in Gúrpide et al. (

Figure 5 .Figure 6 .
Figure 5. (Left) Heii 4686 flux map, constructed by integrating the spaxels around the Heii 4686 line (4688.5-4698.6Å) -accounting for the systemic redshift of NGC1313 -and subtracting the nearby local continuum (using Cube 2 from Gúrpide et al. 2022).The map has been resampled by a factor 3 compared to the original resolution to enhance a feature (blue rectangle) close to the ULX position (white circle).(Right) Average spectra from the blue and orange regions.The orange spectrum has been divided by 1.4 for visual clarity.The spectrum extracted from the patch close to the ULX shows a clear emission line at the expected position of the Heii 4686 based on the redshift of NGC 1313 (indicated by a vertical black dashed line).The magenta dashed line and shaded area shows the best-fit with a Gaussian for the Heii 4686 line and a constant for the local continuum, and its 1 uncertainties.

Figure 7 .
Figure 7. Cloudy-generated 2D best-fit models (contours) and data (background image) comparison for the diskir model with the putative companion star (top panels) and same model but now including the stellar background (bottom panels), for  = 0.15 ⊙ .[OI]6300/H and [OIII]5007/H maps are shown on the left and right panels, respectively.The contours show the Cloudy model predictions, with the numbers for each colored annuli given in the legend.The contours also show the region used in the comparison, where the EUV/X-ray nebula is located.All images show a 22"×22" region centered around the ULX.

Figure 8 .
Figure 8. Effects of including the stellar-background on the [O iii]5007/H (left axis) and [O i]6300/H (right axis) ratios (shown for the phenomenological_star model).The stellar background increases the peak [O iii]5007/H ratio and 'widens' the area over which [O iii]5007 is produced, whereas it has the opposite effect on [O i]6300/H.These 1D plots are for illustration purposes only, as they are not resampled and smoothed to match MUSE spatial resolution.

Figure 9 .
Figure9.(Top) Comparison of 1D radial profiles extracted from the data (solid and shaded areas) and the best-fit 2D Cloudy models (dashed lines), after resampling and smoothing to match MUSE spatial resolution.The profiles were extracted with a width of 10 pixels.For the data, the weighted average and standard deviation are shown as solid lines and dashed regions, respectively.Gaps in the data are due to low signal to noise ratio or bad pixels in those areas.The left axes show the line ratios in low ionisation lines, whereas the right axes (black lines) shows the [OIII]5007/H ratio.The sirf model overpredicts the low ionization lines in the outer part of the nebula.Note the increase in [O iii]5007/H towards the end in the data is due to another nearby group of stars.(Bottom) Range of Heii 4686 flux predicted for each best-fit model (green histograms) averaged over the blue region shown in Figure5, compared to the average value observed in the data in the same region (blue dashed line showing the ±1 uncertainty) and the 13 individual detections (orange histogram).Both diskir and sirf predict too strong Heii 4686 compared to the data, particularly considering the 3 upper limit on Heii 4686.

Figure 10 .
Figure10.Cloudy analysis of the nebular spectrum along the line of sight.(Left) Spectrum extracted around NGC 1313 X-1, averaged over a circular region of 1" in radius roughly matching the cube PSF FWHM.Lines used in the analysis are labelled and their expected position based on the redshift of NGC 1313 are indicated by a black vertical tick.The spectrum has been resampled by a factor 2 compared to MUSE 1.25 Å spectral sampling.(Right) Reduced  2 obtained for the three photo-ionization models with and without the stellar background ( = 2.75 × 10 39 erg/s; Z = 0.15Z ⊙ ).Lines and symbols as per Figure6, including now the results for  in = 60 pc with loosely dotted lines.

Figure B1 .
Figure B1.maps for [Ar iii]7135 (left) and [S iii]9069 (right) derived as per Gúrpide et al. (2022) (see their Figure5).Pixels below a S/N of 5 have been masked.The brighter blob to the far south of the ULX is due to a nearby stellar cluster.

Table 1 .
Derived HST fluxes for the counterpart of NGC 1313 X-1.

Table 2 .
Fits to the Chandra data used to determine the spectral state of NGC 1313 X-1 during the simultaneous HST observations.

Table 3 .
Datasets used to characterised the state-resolved broadband SED of NGC 1313 X-1.

Table 5 .
Cloudy modelling results for the extended nebula for  = 0.15 ⊙ .

Table 6 .
Cloudy modelling results for the emission-line nebula along the line of sight for  = 0.15Z ⊙ .Line fluxes are extinction-corrected and averaged over an extraction region of 1 ′′ .

Table 7 .
Comparison of ULX spectral properties and their high-excitation emission-line nebulae.Notes. Nebular Heii 4686 luminosity.Different values correspond to different measurements reported in the references. Whether unusually high [O iii]5007/H ratios indicative of high-excitation have been reported.
UV luminosity inferred from the nebula when available.

Table C1 .
As per Table5but with different contribution from the background stars. = 0.15 ⊙ .