A multi-wavelength view of distinct accretion regimes in the pulsating ultraluminous X-ray source NGC 1313 X-2

The XMM-Newton observations were carried out with the EPIC-pn and MOS detectors in full-frame imaging mode, and were reduced using the Science Analysis Software (SAS)v.16.¹¹1https://www.cosmos.esa.int/web/xmm-newton/xsa We summarise the essential features of these observations in Table 1. After examining the radial profile to confirm there were no extended components to the point spread function, we extracted source events within an aperture of 30 arcseconds centred on the Chandra source position (RA = 49.592, DEC = -66.601; Zampieri et al. 2004). The background region was chosen to lie adjacent to the source on the same chip. Care was taken to avoid areas with a large internal quiescent background in the EPIC-pn detector, which is primarily traced by the intensity of the Copper-K$\alpha$ emission line (at $\sim$ 8 keV). We excluded events associated with bad pixels by implementing the standard flagging criteria (FLAG==0) and grade selection (PATTERN $\leq 12$ for MOS and PATTERN $\leq 4$ for PN). We extracted light curves in the 10 - 12 keV energy band using the evselect task from chips devoid of a large number of point sources, in order to inspect and remove periods of intense soft proton flares. We ran the SAS task bkgoptrate to determine a safe count rate threshold below which the background level remains steady. As such, periods during which the count rates exceeded this threshold (in the 10 - 12 keV band) were excluded; this resulted in a loss of up to $\sim 46\%$ of the original exposure for the pn camera in the March 24 observation, with other losses being smaller. Finally, using the cleaned event files we extracted the source and background spectra in the 0.3 - 10.0 keV energy range for each instrument. The spectra were grouped into a minimum of 25 counts after background subtraction per energy bin, ensuring not to oversample the intrinsic detector resolution by a factor larger than 3. Hence, it was acceptable to use $\chi^{2}$ statistics to characterise the goodness-of-fit. Finally, in order to deconvolve the energy spectrum through the detector response, we generated the auxillary response files and redistribution matrix files using the arfgen and rmfgen tasks respectively.

The field containing the counterpart to the ULX was observed during two separate HST visits with the WFC3 and ACS detectors in seven independent filters. We summarise the details of these observations in Table 2. We first ran the calwf3 pipeline (version 3.3²²2https://hst-crds.stsci.edu/) on the raw event files of all WFC3 observations using the latest calibration files available in the STScI website³³3The equivalent pipeline calacs was used to reduce the ACS data. The pipeline produces event files that have been bias subtracted, flat fielded and corrected for the effect of charge transfer inefficiency (CTI). We ensured that all sub-exposures were aligned to within 0.1" accuracy (using the drizzlepac tool TweakRegv1.4). This is especially important for the near-infrared (NIR) image, in which the ULX counterpart resides in a crowded environment, such that any misalignments between individual sub-exposures could compromise the photometric accuracy. We then implemented the PSF-fitting technique to derive instrumental magnitudes for the ULX counterpart in all six WFC3 filters using the DOLPHOTv1.2 software (Dolphin 2000). This method identifies point sources above a preset detection threshold from a stack of flat-fielded and cosmic-ray cleaned images. It then uses the TinyTim PSF models (Krist 1995; Hook et al. 2008) to simultaneously fit the PSF of each source and the local background. This approach is preferred over aperture photometry in scenarios where the crowding is severe (e.g. in the near infra-red waveband), which causes the PSFs of closely spaced point sources to overlap. The software performs all the necessary pre-processing steps prior to running the photometry routine (including the masking of bad pixels affected by cosmic rays and image combination). The instrumental magnitudes were initially measured using a circular aperture of 3 pixel radius, and converted to an equivalent 10 pixel (or 0.4") magnitude by applying aperture corrections. To do this, we selected relatively bright and well isolated stars in the field-of-view, with signal-to-noise ratios larger than that of the ULX counterpart. The vertical axis intercept of a linear fit to the 3 pixel versus 10 pixel magnitude of the selected objects yields the aperture correction $\delta$. Here, we emphasise that a fitting function accounting for uncertainties in both the independent and dependent variables (equivalent to the fit_exy routine in IDL; see also Press et al. 1992) was used to ensure that the errors in $\delta$ are not underestimated.

The initial reduction steps were similar for the ACS/SBC data. We used the few bright and well resolved stars in the field of view to align the drizzled SBC image with the F225W and F336W filters, to confirm that the brightest object in the F140LP image is the ULX counterpart. However, we were unable to perform PSF fitting on the ACS/SBC data due to the lack of both a model ACS/SBC PSF in DOLPHOT, and of a sufficient number of well resolved stars in the field of view. Hence aperture photometry was used within an aperture of radius $0.4$ arcseconds (equivalent to the 80% encircled energy for ACS/SBC). We were also hampered by the lack of bright resolved stars in calculating an aperture correction, hence we resorted to using that provided by the instrument manual ⁴⁴4http://www.stsci.edu/hst/acs/documents/isrs/isr1605.pdf. The results of the photometry are quoted in Table 2, and we show a cleaned, 3-colour image of the field that highlights the position of the ULX counterpart in Figure 1.

We began by fitting the XMM-Newton spectra of both observations separately. In order to obtain an understanding of the spectral shape in each observation, and any changes between observations, we fitted the data using a simple model used previously to characterise ULX spectra, i.e. an absorbed multicolour disc (diskbb) plus power-law model in xspec. X-ray absorption was modelled using two instances of the Tuebingen-Boulder ISM component (tbabs; Wilms et al. 2000). We fixed the first component to the Galactic extinction value in the direction of the source (3.97 $\times 10^{20}$ cm²) using the heasarc tool N${}_{\text{H}}$⁶⁶6https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl, and allowed the second to vary. We added cross normalisation factors to account for any differences in calibration between the EPIC and MOS detectors, but these were found to be consistent to within 5 percent. The model provides a good fit to the data as indicated in Table 3. As can be seen in Figure 4 the ratios of the power-law and disc fluxes, in addition to the disc temperatures, suggest that the source spectrum changes from being in the hard ultraluminous (HUL) regime (in the first observation) to a broadened disc (BD) spectrum (see Sutton, Roberts & Middleton, 2013). As discussed in that work, one caveat regarding this spectral classification is that the HUL spectra of some ULXs with pronounced curvature can be difficult to distinguish from BD spectra. While the latter were predominantly found near the Eddington threshold in Sutton et al. (2013), their X-ray luminosity was never part of the identification for the regime, and indeed the realisation that some ULXs have spectra that evolve from HUL to BD-like with increasing luminosity occurred subsequent to that work (e.g. see Luangtip et al. 2016). This appears to also be the case for NGC 1313 X-2, so in what follows we use the term broadened disc to describe the shape of the X-ray spectrum while noting that this does not directly infer any physical meaning.

Critically, the observed change in spectrum supports the notion that there are real physical changes in the ULX between the low, stable 0.3-10 keV flux behaviour in December 2015 and the higher flux, more variable periods in March 2016.

In Section 3.1 we demonstrated that the X-ray emission varies between the two observational epochs (i.e. over a three month period). The HST data were obtained simultaneously to the X-ray data in the two different epochs specifically to test how the UV/optical/IR colours change with the X-ray emission. We present the observed fluxes in each filter at these epochs (converted from the magnitudes given in Table 2) in Table 4, alongside the ratio of the fluxes between the two observations. We also present the X-ray flux separated into three bands; this again demonstrates the variability of the X-ray emission, with the largest variations occurring at medium X-ray energies (1-4 keV).

There appears to be no statistically significant variation exceeding the 3$\sigma$ error range between the two epochs in the UV/optical bands. However, there is some evidence for variability at longer wavelengths, particularly in the NIR (F125W) band, where the source dims by about $25\%$ ($0.33\pm 0.04$ mags) as the X-ray emission rises from epoch 1 to epoch 2. To validate this result, we carefully cross-matched relatively bright and well-isolated sources in both observations in the WFC/IR filter images, and in three other filters. We rejected objects that were too saturated, or those with too many bad pixels (due to cosmic ray hits), by applying a flag threshold within dolphot. Aside from the ULX, our final selection included a number of stars belonging to the OB association, which should be intrinsically invariant on monthly timescales⁷⁷7Although we note that B stars are sometimes surrounded by circumstellar excretion discs due to their fast rotation, which can enhance variability in the infrared.. In any case, the bottom right panel of Figure 5 shows that the ULX counterpart has the largest change in the NIR magnitude (defined as $\Delta$ mag) between the two epochs, when compared to the reference sources; its $\Delta$ mag value is clearly an outlier, and deviates from the mean NIR variability of the sample by about 3 $\sigma$ (see Figure 6). This hints towards a possible first detection of NIR variability from this object, which had only been observed at UV/optical wavelengths prior to this work. We will elaborate on the possible origin of the NIR emission in the discussion.

Given that the optical emission appears to be relatively invariant between the two epochs, we examine whether the UV/optical emission of the source is consistent with arising solely from the donor star. The ULX counterpart is thought to be associated with a young OB star cluster in its vicinity, since it has a similarly large UV/optical luminosity and a distinctly blue colour, as per other members of the cluster (cf. Figure 1). By fitting stellar isochrones to the colour magnitude diagrams (CMDs) of its member stars, Grisé et al. (2008) inferred the cluster age to be $20\pm 5$ Myrs, implying that the mass of the ULX donor is $<16M_{\odot}$. Given this constraint on the counterpart, in addition to its absolute visual magnitude ($M_{\text{V}}\sim-5$) and colour index ($B-V\sim-0.25$), a range of possible OB spectral types were suggested in the work of Grisé et al. (2008). Here, we attempt to refine these constraints by re-fitting the HST magnitudes (see Table 2) with a grid of O/B type stellar templates from Castelli & Kurucz (2003). The temperature of individual spectral types were adopted from Martins, Schaerer & Hillier (2005)⁸⁸8for O-type stars and Straizys & Kuriliene (1981)⁹⁹9for all other spectral types. The normalisation of each template depends on the stellar radius and the source distance, where we allowed the former to vary and fixed the latter to a value of 4.0 Mpc (based on a weighted average of the distances to NGC 1313 in the SIMBAD database ¹⁰¹⁰10http://simbad.u-strasbg.fr/simbad/). We included two additional free parameters, namely the metallic abundance and interstellar reddening (using the Cardelli et al. 1989 extinction curve). We fit the data from both epochs simultaneously given that we have seen little or no change in the flux between epochs in most bands.

The best fitting stellar templates are shown in Figure 7, and the fits are summarised in Table 5. A B0III giant provides the most reasonable fit (i.e. $\chi^{2}$ = 14.5 for 9 degrees of freedom), but only if the NIR data points are excluded, with the fit otherwise worsening by $\Delta\chi^{2}=20$ for 2 degrees of freedom. A B0Ia supergiant may have a better reduced $\chi^{2}$, but the best-fit stellar radius is at least a factor of 2.5 smaller than expected for a star of this type (Straizys & Kuriliene 1981). A main-sequence O6V star can provide a reasonable fit, but it is a worse fit to the data than B-star models with $\chi^{2}=24.9$ for 9 degrees of freedom.

We also attempted a reddened power-law model as a fit to the data. It strongly over-predicts the far UV flux, but provides an acceptable fit at larger wavelengths. This is intriguing for interpreting the optical emission, as reprocessed emission will have a power-law-like spectrum at optical wavelengths. Hence it is clear from the optical data alone we cannot unambiguously identify the physical origin of the counterpart’s emission.

Despite invariance of the optical emission in the current data set, Grisé et al. (2008) used archival VLT observations to report that the optical counterpart to the ULX shows significant variability in the B-band (by $\sim 0.2$ mag) on timescales ranging from several hours to days, with respect to other non-variable stars in the field. These timescales are much faster than that at which the mass accretion rate varies at the outer disc, which occurs on the viscous timescale:

where $\alpha$ is the viscosity coefficient, $H$ is the disc scale-height at radius $R$ and $t_{\text{dyn}}$ is the Keplerian dynamical timescale. We can crudely estimate that $\tau_{\text{visc}}$ is of order a year assuming $\alpha\sim H/R\sim 0.1$ for a geometrically thin outer disc, whose size is determined by the orbital period of the system containing a neutron star accretor (which we assume is $\sim 2$ days; Sathyaprakash et al. 2019) and a B-type companion (with a mass $\sim$ 10 M_⊙). The latter is not tightly constrained, but even if it is taken to be two orders of magnitude smaller, $\tau_{\text{visc}}$ remains much larger than the B-band flickering timescale. Variation of the outer disc itself cannot therefore explain the previous reports of optical variability.

A possible alternative to explain this is via a disc reprocessing model: variable X-rays from the inner accretion flow can illuminate the outermost regions of the disc (Gierliński, Done & Page 2009), where a fraction of the incident X-rays are absorbed and heat the underlying disc material. The highest energy photons penetrate deeper into the disc atmosphere and can thermalise with the surrounding material, re-emerging with a blackbody continuum peaking in the UV/optical band (depending on the local temperature). In this scenario, the longer wavelength emission should vary coherently with the inner X-ray emission across a broad range of timescales, although the variability will be suppressed on the shortest timescales due to a finite propagation time between the two regions. In particular, it will be limited by the light travel time in the simple picture of X-rays directly illuminating a geometrically thin outer disc. If the disc is tidally truncated with an outer radius $R_{\text{out}}$ of $\sim$ 0.7 times the radius of the neutron star’s Roche lobe, the light travel time is $\lesssim$ a minute (for the orbital period estimate used above). Thus, one might expect changes in the UV/optical and X-ray fluxes to be correlated over the much longer three month (inter-observational) timescale. However, a comparison of the fluxes in the two separate bands demonstrates no evidence for such a correlation. The broadband 0.3-10.0 keV flux increases by $\sim 30\%$ between the two observations, while the UV/optical fluxes are mostly consistent between the two epochs and only change significantly (and in the opposite sense) for the longest wavelength filter (see Table 4 and Figure  5). Therefore, we do not observe any correlated X-ray-to-optical long-term variability in the system, at least between the two epochs at which our new observations were taken. Nevertheless given the previous evidence we cannot rule out X-ray reprocessing as a contributor to the optical emission, and so seek to constrain its presence using further techniques.

We have seen that the X-ray emission changes between epochs but that the UV/optical/IR largely doesn’t; that the optical emission alone could be described by either stellar light or reprocessed X-ray emission; but that a state-of-the-art sub-Eddington irradiation model does not provide a meaningful constrain on the possible reprocessing. We therefore attempt to construct a new model that takes into account the super-Eddington nature of the central regions and that can also account for the optical light.

Our prime means of achieving this is to replace the multi-colour disk blackbody component in the X-ray spectrum with the physically-motivated slim disc model of Kubota & Done 2019, which is referred to as agnslim. In previous studies of the two-component ULX X-ray spectra, this component has been interpreted as the photospheric region at the base of the wind (Kaaret et al. 2017 and references therein) although more recent models based on magnetic neutron stars have more classically attributed it to the outer accretion disc (Koliopanos et al., 2017)¹¹¹¹11This attribution is very problematic, given that the outer discs are assumed to be ’normal’ optically-thick geometrically-thin sub-Eddington discs (which means their emission cannot be substantially beamed), yet still apparently radiate well in excess of the Eddington limit for a neutron star.; in broadened disc spectra this component has frequently been attributed to an advection-dominated slim disc (Soria et al., 2015). This model suppresses the luminosity released by the inner accretion flow, relative to that predicted by the accretion rate $\dot{M_{\text{in}}}$ through the outer disc. Hence, beyond some critical value of $\dot{M_{\text{in}}}$ the intrinsic X-ray luminosity of a slim disc saturates to $\sim 10$ $L_{\text{Edd}}$. The model retains the same three-radial zone structure as optxrplir, with two distinct (warm and hot) coronal regions. However, a critical assumption is that the puffed-up disc at the radius $r_{\text{c}}$ self-shields the hard X-ray emission from the inner accretion flow, such that any reprocessed emission from the outer disc is negligible (see Figure 12 of Kubota & Done 2019 for a visualisation of this effect).

We made slight adjustments to some of the model parameters before it was appropriate for fitting the multi-wavelength spectrum of a neutron star ULX. These sources can have strong magnetic fields that impart an opposing pressure to the inward ram pressure of the infalling material, thus halting the accretion flow at a radius $R_{\text{M}}$ much larger than the innermost stable orbit $R_{\text{isco}}$. The value of $R_{\text{M}}$ depends on the dipolar field strength and the mass accretion rate at that radius (see King et al. 2017). So, we enforced the inner disc radius of agnslim to take on values much larger than $R_{\text{isco}}$ (while still allowing it to vary). Moreover, we modelled the emission from the accretion envelope surrounding the NS surface (reference) using the thermal comptonisation model (nthcomp) following Middleton et al. (2019). Further, we adopted a canonical NS mass of 1.5 M_⊙, and retained the same neutral gas absorption and dust reddening components as in the previous section.

In Table 7, we present the best-fit parameters for the two cases where (i) we fit the multi-wavelength SED of X-2 with a composite model including agnslim, nthcomp and a single-temperature blackbody (bbodyrad) to account for the donor star’s emission, and (ii) when the latter is ignored. Discarding any stellar contribution in the UV/optical yields an upper limit to $\dot{m}$ that we annotate as $\dot{M}_{0}$. But, we emphasise that due to uncertainties in the donor spectral type and any other contributions to the HST spectrum, the latter will invariably have some unknown systematic uncertainty. We set the temperature and normalisation of bbodyrad to match a B0III spectral type, since this gives the most reasonable (although not an ideal) fit to the UV/optical magnitudes in section 3.3. We illustrate the fits to the SED of NGC 1313 X-2 with this new model both with and without a stellar contribution in Figures 9 and 10 respectively.

Although the models generally provide good fits across the SED, it is noticeable that the IR datapoint presents an outlier to both models. This is because neither the stellar blackbody produces sufficient flux at these longer wavelengths to match the data, nor does the accretion disc become large enough (given the suspected orbital constraints) for its outer radii to radiate with a characteristic temperature in this band. We must therefore seek a different explanation for the flux in this band, and do so in the Discussion section.

Table 7 shows that statistically acceptable fits are obtained for both cases (i) and (ii), with the major difference being the requirement of a substantially larger $\dot{M}_{0}$ when the stellar emission is not considered (i.e. $\gtrsim 6000$ times the Eddington value). One reason to disfavour this solution is the following: for a neutron star accretor, any advected energy must be released near their surface (unlike the case of a black hole, where it falls beyond the event horizon). In the absence of strong disc winds (as the Lasota et al. 2016 prediction implies), we can estimate the luminosity radiated near the surface to be $\sim 5\times 10^{41}$ erg s^-1, using the accretion efficiency applicable to a neutron star ($\eta=\frac{GM_{\text{NS}}}{R_{\text{NS}}c^{2}}\sim 0.2$) and the lower best-fit $\dot{M}_{0}\sim 6000\dot{M}_{\text{Edd}}$. This is clearly at odds with the much smaller observed luminosity of the thermal comptonisation component ($3\times 10^{39}$ erg s^-1), implying that the outflow rate must be more than than 95% of the accretion rate through the outer disc, and that implausibly strong winds are required¹²¹²12Estimates from ULX bubbles, including for NGC 1313 X-2, are that winds carry $\lesssim 10^{40}~\rm erg~s^{-1}$ of kinetic energy (e.g. Roberts et al. 2003), over an order of magnitude smaller than required here.. Admittedly, we could be underestimating the luminosity from the accretion column above 10 keV, with analyses of the pulsed emission of other NS ULXs indicating that it peaks above the bandpass of XMM-Newton. However, this discrepancy is likely to be well within a factor of 2 (see e.g. Walton et al. 2018b). Therefore, the best-fit mass inflow rate inferred by only considering the slim disc contribution at UV/optical wavelengths leads to inconsistencies on the following front: the prediction of an unreasonably large luminosity from the NS surface and/or an implausibly powerful wind. We note that the best-fit $\dot{m}$ is also influenced by the inclination $i$, since the projected emission from the accretion flow scales as $\cos i$ (i.e. larger viewing angles will require an increase in $\dot{m}$ to fit a given spectrum). Setting the inclination to an extremal 0 degrees (fully face on view) does reduce the best-fit $\dot{m}$, but not to the extent that the above problems are alleviated. Hence, we strongly favour the fits including a stellar component on a physical basis.

The derived mass accretion rate $\dot{m}$ hardly changes between the two observations. This is expected given a lack of observed variability in the UV/optical bands, which correspond to emission from the outermost regions of the disc that are much more sensitive to changes in $\dot{m}$ (since photon trapping effects are unimportant here). Despite an unchanged mass accretion rate, the inner disc radius takes on a much smaller value in the second epoch (decreasing by a factor $\sim 5$), which allows agnslim to account for an increase in the 1.0 - 6.0 keV flux. This is difficult to understand if $R_{\text{in}}$ is interpreted as the magnetospheric radius, which scales primarily with the mass accretion rate and magnetic field strength as $\dot{M}^{-\frac{2}{7}}B^{\frac{4}{7}}$ (Ghosh & Lamb 1979). Unless the field strength decreases considerably (by a factor $\sim 10$) on the three month timescale, there is no reason to expect any variation in $R_{\text{M}}$.

NGC 1313 X-2 is remarkable in that it is one of the half-dozen well-established PULXs, albeit the one with the faintest pulsations (Sathyaprakash et al., 2019). The presence of a large bubble nebula around this object strongly argues for the occurrence of outflows and/or jets, an expectation of supercritical accretion models. In this work, we explored whether the two distinct modes of X-ray spectral behaviour it shows (as described in section 3.1) are also consistent with such models, derived originally from Shakura & Sunyaev (1973) and adapted for the presence of a magnetic neutron star by Mushtukov et al. (2019). We note that to fully understand this, exploration of the luminosity-temperature plane (e.g. Robba et al. 2021) via a continual monitoring in X-rays and at other wavelengths is necessary, rather than the two snapshots we have obtained; however we will proceed with this caveat in mind.

We notice that the change in X-ray spectral shape between the two observations occurs primarily through an increase in the 1.0-6.0 keV flux, while the lowest ($<1$ keV) and highest energy ($>6$ keV) components remain stable. Such an occurrence is qualitatively similar to at least two other ULXs (Holmberg IX X-1 and NGC 1313 X-1; Luangtip et al. 2016; Walton et al. 2020), and is difficult to reconcile with an intrinsic change in the mass accretion rate alone (see Table 7). This is supported by the seemingly invariant UV/optical and soft X-ray fluxes, which are (in part) emitted by the outermost regions of the disc, implying that the mass accretion rate through the outer disc remains steady over the three month timescale. As noted by Luangtip et al. (2016), the observed X-ray spectral transitions may instead be caused by a bulk precession of the thick inner disc and wind, plausibly due to the Lense-Thirring (or frame dragging) effect (Middleton et al. 2018). This naturally changes the observer’s viewing angle to different parts of the accretion flow; a line-of-sight corresponding to the inner-most regions presents a view of the highest energy photons, which are geometrically collimated by the ultra-fast wind. Conversely, a view of outer wind photosphere reduces the hard X-ray flux, as the highest energy photons are Compton down-scattered by the outflowing material and/or scattered out of the line-of-sight. This picture is supported both by observational studies of how the soft X-ray spectral residuals correlate with luminosity (Middleton et al. 2015b), in addition to numerical simulations of super-critical discs (Ohsuga & Mineshige 2011; Sądowski & Narayan 2016; Narayan, Sądowski & Soria 2017).

Clearly, the precession scenario also predicts a stronger modulation in the hard X-ray band over long timescales, since the soft X-rays are less geometrically beamed and distributed over a larger solid angle, and should thus be less variable as the wind moves in/out of the observer’s sight line. Indeed, recent studies have reported quasi-periodic hard X-ray variability in several (P)ULXs on weekly-monthly timescales (e.g. Weng & Feng 2018; Brightman et al. 2019; Vasilopoulos et al. 2020), invoking the above precession model as a possible explanation. In the particular case of X-2, Weng & Feng (2018) conclude that the variability is strongly periodic, although they do not consider red noise into their significance estimation, and the periodicity is not fully obvious in their folded light curve. In any case, Ogilvie & Dubus (2001) discuss that a radiatively warped disc may undergo unstable precession, which need not induce strictly periodic modulations on the X-ray flux (although it is currently unclear whether this mechanism is preferred over the general relativistic Lense-Thirring effect for ULXs).

The key question is why, if the X-ray emission changes, the optical does not? Our analysis suggests that this may be due to the optical emission being dominated by the companion star, with a weaker contribution from an irradiated outer disc (although other factors may also play a role, as we discuss below). Indeed, this is perhaps not unexpected if the hard X-rays are collimated away from the accretion disc and (presumably) the companion star, and so cannot be reprocessed. Hence, precession cannot affect the optical appearance of the PULX in the same way that it can affect the X-ray emission.

Our analysis of the optical light from the ULX counterpart shows it to be largely invariant between the two observational epochs, despite a change in the 0.3-10 keV X-ray flux. Using the magnitudes obtained in each of the different HST filters we can obtain a coarse spectrum of the counterpart, finding that its shape is broadly consistent with a B0III giant (if the NIR data point is ignored), but is also well explained by a power-law model (if the far UV data point is ignored), which is reminiscent of reprocessed X-rays. We therefore attempted to fit the broadband X-ray - UV/optical SED with a disc reprocessing model, but this provided an unreasonably large reprocessed fraction ($>80\%$ of the X-ray flux), although this is at least partly attributable to the unsuitability of the model to fitting super-Eddington accretion flows. While a ‘slim’ disc model is more justified in this respect, it does not consider X-ray reprocessing, and implicitly assumes that the outer disc is largely shielded from hard X-rays by the inflated inner disc. However, usage of this model implied that either unreasonably high accretion rates are required to account for the observed levels of UV/optical emission, or that an additional stellar component (dominating the UV/optical flux) is necessary. We emphasise that, while the light from the star dominating the counterpart appears probable from our results, we cannot fully rule out reprocessing since it is not accounted for by the slim disc model.

Even if the optical light were to arise in reprocessed X-ray emission, changes in the reprocessed fraction may be difficult to observe in our case, given that the broadband X-ray flux increases only by $\sim 25\%$ between the observations probed by XMM-Newton. In the simple disc geometry applicable to sub-Eddington low mass X-ray binaries with a geometrically thin disc, $f_{\text{out}}$ can be predicted analytically. It is a function of the X-ray flux irradiating the outer disc $f_{\text{out}}\propto\big{(}L_{\text{X}}\big{)}^{1/7}$ (Cunningham 1976), for a given compact object mass and disc size. This relation implies that a 25% increase in $L_{\text{X}}$ should lead only to a 3% increase in $f_{\text{out}}$, which cannot be demonstrated given the current quality of our data (see Table 6). While this is consistent with the lack of variability we see, our accretion disc should not be in a geometrically thin configuration (as emphasised before). Indeed, it is noted above that in supercritical discs one might expect the hard X-rays to be funnelled away from the disc, making large reprocessing fractions unlikely.

A final point to consider on reprocessing is that the longer wavelength emission may respond to changes in the X-ray flux with a time delay that is much larger than the duration of the snapshot HST observations ($15-20$ mins). This would clearly mask any correlated variability between these bands, although a crude estimate of the light-travel time between the X-ray and optically emitting regions of the disc may suggest that it is unlikely (see section 3.4). Indeed, our estimate does not account for the emergence of an ultrafast wind at accretion rates above Eddington (Pinto et al. 2016, Kosec et al. 2018b). The lag between the X-ray and UV/optical flux could then be influenced by the photon diffusion timescales across the outflow, which will partly depend on the density of material within it. Some evidence for such lags in ULXs is already emerging within the X-ray regime (Kara et al., 2020).

Our modelling does however favour the optical emission arising in the companion star. It is worth noting that the best-fitting stellar model reported here may not reflect the true donor spectral type, partly due to reasons outlined in Patruno & Zampieri (2008). They demonstrate that the luminosity-temperature evolution of an isolated star can substantially differ from one in a binary. For example, an isolated star on the main sequence will tend to become more luminous as its radius expands. By contrast, a binary main sequence star (of the same initial mass) cannot expand above its Roche Lobe radius, and tends to decrease both in luminosity and temperature at the onset of Roche Lobe overflow. It will thus appear cooler and fainter than an isolated star at the same stage of its evolution, posing problems for a direct comparison of the observed magnitudes and colours with standard stellar templates. In addition, the donor star can also be irradiated by X-ray photons from the accretion flow, further affecting its brightness and colour (e.g. Motch et al. 2014, although see comments above). Hence, our attempts at identifying the stellar type may not be very accurate.

Finally, we note that our work presents an intriguing link to other PULXs. The first ULX for which the counterpart was unambiguously demonstrated to be stellar was NGC 7793 P13, also a PULX (Motch et al., 2014). A red supergiant companion has now been revealed for a second PULX, NGC 300 X-1 (Heida et al., 2019b), although attempts to uncover the counterparts for two more PULXs have proven inconclusive due to their distance, being heavily affected by obscuration and/or residing in crowded locations (Heida et al., 2019a). Thus at least half of the PULXs have plausible stellar-dominated counterparts. This contrasts strongly with the remainder of the ULX population where reprocessed X-ray emission remains a strong explanation for the optical counterpart (Tao et al., 2011). This could perhaps support the hard X-ray beaming in NS ULXs predicted by King (2008), and seems consistent with the notion that NS ULXs would have smaller accretion discs than in their black hole equivalents (Heida et al., 2019a).

In almost all of our modelling, the NIR data point from the HST F125W filter appears to sit above the expected model emission, demonstrating a likely NIR excess in this ULX. It is also notable that this filter shows the strongest evidence for variability in any of the HST data between the two observational epochs, with a $\sim 35\%$ decrease in NIR flux observed. This is in stark contrast to an increase in X-ray flux between the same two epochs.

The fact that this variability occurs only in the NIR (whilst being absent in the UV/optical) immediately rules out an origin due to X-ray reprocessing, or arbitrary changes in the line-of-sight interstellar extinction. It is also unlikely to come from the counterpart being a red supergiant star, which appears to be excluded by the observed SED¹³¹³13We note that López et al. (2020) catalogue a possible red supergiant companion to X-2. However, this object (their object B) is only just within their formal error region. They also detect the object we regard as the counterpart (their object A) much closer to the centre of their error region, and derive magnitudes from the HST data consistent with ours.. We speculate that it may instead be attributed to transient jet emission, although we require additional evidence to confirm this. Classical Galactic black hole binaries (GBHBs) sometimes show NIR flaring episodes during accretion state transitions, from the low-hard state (featuring a persistent radio jet) to the high-soft/thermal dominant state (when jet emission is suppressed). These episodes are hypothesised to be caused by ejections of material from the hot inner flow. Similar potentially explosive jet formation has also been seen in more than one ULX, for example in Ho II X-1 (Cseh et al., 2015) and in a transient ULX in M31 (Middleton et al., 2013).