Thermal stability of winds driven by radiation pressure in super-Eddington accretion disks

TL;DR

This paper models the thermal stability of radiation-driven winds in super-Eddington accretion disks of ULXs, revealing stable wind conditions and explaining observed spectral features and correlations.

Contribution

It presents the first calculation of photoionization balance in winds from thick, super-Eddington accretion disks, highlighting their thermal stability and spectral behavior.

Findings

Winds are generally thermally stable in super-Eddington accretion.

Variations in accretion rate and precession affect wind appearance.

Model explains correlation between spectral residuals and ULX spectral state.

Abstract

Ultraluminous X-ray sources (ULXs) are mainly powered by accretion in neutron stars or stellar-mass black holes. Accreting at rates exceeding the Eddington limit by factors of a few up to hundreds, radiation pressure is expected to inflate the accretion disc, and drive fast winds that have in fact been observed at significant fractions of the speed of light. Given the super-Eddington luminosity, the accretion disc will be thicker than in sub-Eddington accretors such as common active galactic nuclei and X-ray binaries, leading to a different spectral energy distribution and, possibly, a different thermal status of the wind. Here we show the first attempt to calculate the photoionization balance of the winds driven by strong radiation pressure in thick discs with a focus on ULXs hosting black holes or non-magnetic neutron stars. We find that the winds are generally in thermally stable equilibrium, but long-term variations in the accretion rate and the inclination due to precession may have significant effects on the wind appearance and stability. Our model trends can explain the observed correlation between the spectral residuals around 1 keV and the ULX spectral state. We also find a possible correlation between the spectral hardness of the ULX, the wind velocity and the ionization parameter in support of the general scenario.

Tables

Table 1. Table 1: Characteristics of the brightest ULXs and XMM- Newton exposure times

Source	d (Mpc) (a)	$L_X^{40}$	Flux	$HR$	Wind	$v_w$	$\log \xi$	t${}_{}{}^{\mathrm{tot}}{}_{\exp }{}^{}$ (b)	Similar sources
NGC 300 ULX-1	2.0	$0.5$	$4$	$0.7$	D	$0.24$	$3.9$	212	NGC 7793, NGC 5907, M 82 (PULXs)
Holmberg IX X-1	3.8	$4$	$15$	$0.6$	T	$0.23$	$4.5$	161	IC 342 X-1, NGC 1313 X-2
NGC 1313 X-1	4.2	$2$	$5$	$0.5$	D	$0.19$	$2.3$	750	NGC 5204 X-1
NGC 5204 X-1	4.8	$1$	$2$	$0.4$	D	$0.34$	$3.0$	160	NGC 1313 X-1
NGC 5408 X-1	5.3	$1$	$3$	$0.2$	D	$0.20$	$2.1$	655	NGC 6946 X-1, Holmberg II X-1
NGC 55 X-1	2.0	$0.2$	$3$	$0.1$	D	$0.19$	$3.3$	120	NGC 247 X-1, M81 ULS

Table 2. Table 2: List of non-ULX sources used for comparison.

Source	Source type	$L_{\mathrm{Edd}}$ peak (a)	$M_{\odot }$
NGC 5548	Seyfert 1 / AGN	0.05	$4\times {10}^7{}^{\left(b\right)}$
IRAS 13224	NLS1 / AGN	5	$1\times {10}^6{}^{\left(c\right)}$
V2491 Cygni	Nova / WD	1	${1.3}^{\left(d\right)}$
Blackbody	50 eV model	–	–

Table 3. Table 3: XMM- Newton observations used to build ULX SEDs.

Source	State	$HR$	OBS_ID	t${}_{}{}^{\mathrm{tot}}{}_{\exp }{}^{}$ (ks)
NGC 300	Hard	0.7	0791010101	201.0
NGC 300	Hard		0791010301
NGC 1313	Int-Hard	0.5	0405090101	345.6
NGC 1313	Int-Hard		0693850501
NGC 1313	Int-Hard		0693851201
NGC 1313	Bright-BB	0.35	0301860101	21.3
NGC 1313	Low-Soft	0.04	0205230401	14.3
NGC 5408	Soft	0.2	0302900101	644.9
NGC 5408	Soft		0500750101
NGC 5408	Soft		0653380201
NGC 5408	Soft		0653380301
NGC 5408	Soft		0653380401
NGC 5408	Soft		0653380501

Table 4. Table 4: Wind properties and spectral hardness in ULXs

Source	$HR$	$v_w$	$\log \xi$	REF
NGC 300 ULX-1	$0.68\pm 0.03$	$0.24\pm 0.01$	$3.9\pm 0.2$	[5]
Holmberg IX X-1	$0.61\pm 0.03$	$0.23\pm 0.01$	$4.5\pm 0.2$	[4]
NGC 1313 X-1 (a)	$0.46\pm 0.02$	$0.23\pm 0.02$	$4.5\pm 0.2$	[1,3]
NGC 1313 X-1 (b)	$0.46\pm 0.02$	$0.19\pm 0.01$	$2.3\pm 0.1$	[1,3]
NGC 5204 X-1	$0.40\pm 0.02$	$0.34\pm 0.01$	$3.0\pm 0.2$	[4]
Holmberg II X-1	$0.24\pm 0.01$	$0.20\pm 0.01$	$3.3\pm 0.2$	[4]
NGC 6946 X-1	$0.23\pm 0.01$	$0.08\pm 0.01$	$2.0\pm 0.2$	[1,4]
NGC 5408 X-1 (a)	$0.19\pm 0.01$	$0.20\pm 0.01$	$2.1\pm 0.2$	[1]
NGC 5408 X-1 (b)	$0.19\pm 0.01$	$0.10\pm 0.01$	$1.7\pm 0.2$	[1]
NGC 55 X-1 (a)	$0.11\pm 0.01$	$0.19\pm 0.01$	$3.3\pm 0.2$	[2]
NGC 55 X-1 (b)	$0.11\pm 0.01$	$0.06\pm 0.01$	$0.5\pm 0.3$	[2]
NGC 247 X-1	$0.03\pm 0.01$	$0.14\pm 0.01$	$3.0\pm 0.3$	[2]

Table 5. Table 5: Wind properties and spectral hardness correlations

Parameters	Spearman	Pearson
$HR-v_w$	$0.71$	$0.61$
$HR-\log \xi$	$0.58$	$0.62$
$v_w-\log \xi$	$0.76$	$0.70$

Equations

$$\xi =\frac{L_{\mathrm{ion}}}{n_{\mathrm{H}}{R}^2}$$

$$\begin{array}{rlrl}{L}_w& =0.5{\dot{M}}_w{v}_w^2& & \\ & =2\pi m_p\mu \Omega C\frac{v_w^3}{\xi }{L}_{\mathrm{ion}}\sim 10^{39-41}erg/s& & \end{array}$$

$$v_w=\sqrt{2L_w/{\dot{M}}_w}\approx 0.19c\sqrt{L_{39}{\tau }_6/\delta m}$$

$$R_{\mathrm{bubble}}=0.76{\left(\frac{{\dot{M}}_w{v}_w^2}{2{\rho }_{\mathrm{ISM}}}\right)}^{1/5}{\tau }^{3/5}$$

$$\begin{array}{rlrl}{v}_w& =(0.06\pm 0.04)+(0.46\pm 0.15)HR& & \\ \log \xi & =(1.5\pm 0.5)+(3.9\pm 1.5)HR& & \end{array}$$

Full text

††pagerange: Thermal stability of winds driven by radiation pressure in super-Eddington accretion discs–Thermal stability of winds driven by radiation pressure in super-Eddington accretion discs††pubyear: 2018

1 Introduction

It is thought that when the Universe was young, black holes were accreting matter at impressive rates in order to build up “fully grown” supermassive black holes with masses of $\gtrsim 10^{9}\,M_{\odot}$ in a few hundred million years (see, e.g., Fan et al. 2003, Volonteri et al. 2003). At accretion rates approaching the Eddington limit, $\dot{M}_{\rm Edd}$, radiation pressure inflates the thin disc producing a geometry similar to that of a funnel (see, e.g., Shakura & Sunyaev 1973) and a wind is launched (see, e.g., Poutanen et al. 2007). These systems significantly differ from most common active galactic nuclei (AGN) and X-ray binaries (XRBs) which accrete at sub-Eddington rates (see, e.g., Grimm et al. 2002; Aird et al. 2018). Winds driven by radiation pressure in a super-Eddington (or supercritical) regime can reach mildly-relativistic $\sim 0.1c$ velocities (see, e.g., Takeuchi et al. 2013) and therefore may carry out a huge amount of matter and kinetic power, potentially altering both the accretion onto the black hole and the star-formation in the host galaxy (feedback, see, e.g., Fabian 2012 and references therein).

1.1 Highly-accreting supermassive black holes

Nowadays, it is still difficult to study the most distant ($z>6$) and rare quasars, but there is evidence that some nearby supermassive black holes are accreting at significant fractions of the Eddington limit or beyond such as those powering Narrow Line Seyfert 1 galaxies (NLS1, e.g. Boroson 2002). However, it is difficult to determine the bolometric luminosity in Eddington units of NLS1s as most disc emission is produced in the UV, significantly affected by interstellar absorption and several systematic effects by up to an order of magnitude, similarly to the highly uncertain mass estimates (see, e.g., Vasudevan & Fabian 2007; Jin et al. 2012; Castelló-Mor et al. 2016; Buisson et al. 2018).

Evidence of relativistic winds in NLS1 has indeed been found (see, e.g., Hagino et al. 2016; Kosec et al. 2018c) and in some cases the strength of the absorption lines correlates with the variability of the ionizing continuum as expected for radiatively-driven winds (see, e.g., Matzeu et al. 2017; Parker et al. 2017; Pinto et al. 2018). However, these “ultrafast outflows” are very similar to those found in several AGN accreting at about the Eddington limit or 1-2 orders of magnitude below (see, e.g., Pounds et al. 2003; Reeves et al. 2003; Tombesi et al. 2010), which could be driven by magnetic fields (Fukumura et al. 2017) rather than radiation pressure (King & Pounds 2015). It is not easy to find and compare unique observable predictions of the theoretical wind models with the observations of winds. The study of the wind thermal stability can be a step forward in this field and, therefore, it is proposed in this paper.

1.2 Ultraluminous X-ray sources

The detailed study of accretion physics at extreme rates has been boosted in the last decade by extraordinary discoveries in the class of objects known as Ultraluminous X-ray sources (ULXs), currently the best examples of compact objects surpassing the Eddington limit for long periods of time. ULXs are bright, point-like, off-nucleus, extragalactic sources with X-ray luminosities above $10^{39}$ erg/s, resulting from accretion onto a compact object (Kaaret et al., 2017). Recent studies have shown that some ULXs are powered by accretion onto neutron stars with strong magnetic fields ($10^{9-14}G$, e.g. Bachetti et al. 2014, Fürst et al. 2016, Israel et al. 2017a, Israel et al. 2017b, Tsygankov et al. 2016, Carpano et al. 2018, Brightman et al. 2018, Rodriguez Castillo et al. 2019, Sathyaprakash et al. 2019, Middleton et al. 2019), which confirmed earlier speculation that ULXs are in majority stellar-mass compact objects $(<100M_{\odot})$ at or in excess of the Eddington limit (King et al. 2001, Poutanen et al. 2007, Gladstone et al. 2009, Middleton et al. 2013, Liu et al. 2013, Motch et al. 2014). A few ULXs might still host intermediate mass black holes ($10^{3-5}M_{\odot}$) at more sedate Eddington ratios (see, e.g., Greene & Ho 2007, Farrell et al. 2009, Webb et al. 2012, Mezcua et al. 2016).

Before the discovery of ultraluminous pulsars (PULX) it was already thought that most ULXs were super-Eddington due to the presence of a strong turnover below $\sim 7$ keV in most ULX spectra (see, e.g., Stobbart et al. 2006, Bachetti et al. 2013). In the softest states of ULXs (e.g. in NGC 55 ULX) the spectral curvature may start as low as 0.7 keV (see, e.g., Gladstone et al. 2009). It is difficult to explain the presence of this turnover in combination with the high $10^{39-41}$ erg/s luminosities invoking sub-Eddington accretion models. ULXs represent a very complex category of objects with spectra of different shapes and slopes. In Fig. 1 we show a hardness sequence of ULX X-ray spectra uncorrected for interstellar absorption similar to Pintore et al. (2017) and Pinto et al. (2017). When modelled with a power-law, the spectra show the photon index ranging from $\sim 1$ for the hardest ULXs, to 2.5 for the soft ULXs and above 4 for the ultraluminous supersoft sources (ULS or SSUL, see, e.g., Feng et al. 2016, Urquhart & Soria 2016).

In general, ULX broadband X-ray spectra require 2-3 emission components (with the hardest component possibly related to the accretion column onto the poles of a neutron star since it dominates the energies with the largest pulsed fractions in the known ULX pulsars when pulsations are detected; Walton et al. 2018a). A broad blackbody-like emission with temperatures around 1 keV is likely produced by the supercritical thick disc itself, whose photon emission is highly distorted by hot electrons in the inner regions (up-scattering) and by down-scattering with cool electrons in the wind (see, e.g., Gladstone et al. 2009, Sutton et al. 2013 and Middleton et al. 2015). In the case of highly-magnetized neutron stars with moderate accretion rates a substantial fraction of the emission around 1 keV may come from near the neutron star surface. A soft component with temperatures around 0.1 keV most likely arises from upper regions of the outer disc, where a radiatively-driven wind is expected to be launched at accretion rates comparable to or higher than the Eddington limit (see, e.g., Ohsuga et al. 2005; Poutanen et al. 2007; Takeuchi et al. 2013). It is also possible that some portions of the wind, particularly phases of the outflow with much lower velocities, will be launched by thermal heating close to the Compton temperature (e.g., Begelman et al. 1983; Done et al. 2018). Outflows can also be launched by strong magnetic fields (see, e.g., Romanova & Owocki 2015; Romanova et al. 2018; Parfrey et al. 2017).

1.2.1 The winds discovery and the unification scenario

The first discovery of powerful winds possibly driven by radiation pressure in ULXs was achieved by the recent detection of blueshifted absorption lines in high-resolution soft X-ray spectra provided by the XMM-Newton Reflection Grating Spectrometers (RGS, Pinto et al. 2016; Pinto et al. 2017; Kosec et al. 2018a). Further confirmation was obtained in the Fe K hard X-ray energy band (Walton et al. 2016) and, in particular, for the ultraluminous X-ray pulsar NGC 300 ULX-1 (or NGC 300 PULX) both in the soft and the hard energy bands (Kosec et al., 2018b).

Several features were already seen, albeit spectrally unresolved, in the low-to-moderate resolution CCD spectra (see, e.g., Stobbart et al. 2006). The study of high-spatial-resolution Chandra data (Sutton et al., 2015) and timing variation (Middleton et al., 2015) suggested that the residuals are mainly produced within a small region and are intrinsic to the ULX itself, providing the first evidence of a wind (see also Middleton et al. 2014).

On the basis of this theoretical, computational and observational evidence, a unification scenario has been proposed according to which ULX spectral properties depend on their accretion rates and inclination angles (see, e.g., Gladstone et al. 2009, Middleton et al. 2011, Middleton et al. 2015, Sutton et al. 2013) and on the ionization state and complexity of the wind in the line of sight (Pinto et al., 2017).

Fig. 2 shows an artistic impression of a wind blowing from the supercritical region around the spherization radius, $R_{\rm sph}$, of an accretion disc of a black hole. In general, we should expect high ionization states of the wind at low inclinations from the disc rotation axis (on-axis) and progressively lower ionization at higher inclinations. However, the outflow is at least in part clumpy and optically thick, which implies that each line-of-sight (LOS) is dominated by a fraction of the wind exposed to a certain ionizing field and therefore characterized by a different thermal equilibrium.

A broadly similar scenario might apply for neutron star ULXs, provided that the magnetic field is not so strong that it truncates the disc before the thick inner regions form and launch the associated disc wind (see, e.g., Takahashi & Ohsuga 2017). Mushtukov et al. (2019) demonstrated that this should be the case for magnetic fields of at least $\sim 10^{12}$ G (and below). While there are cases which may have much stronger magnetic fields (Tsygankov et al. 2016; Brightman et al. 2018), in general the field strengths in neutron star ULXs are still hotly debated (see, e.g., Middleton et al. 2019 and King & Lasota 2019), and there are also cases where there are indications that the field could be much more modest (e.g. Walton et al. 2018b). Indeed Kosec et al. (2018b) report the detection of a powerful outflow in the ULX pulsar NGC 300 ULX-1, implying that the field strength and configuration permit these winds in at least some cases. Nonetheless, the scheme in Fig. 2 is optimal for black holes and non-magnetic neutron stars.

So far, no work has focused on the study of the thermal state of the winds in ULXs and its effect on the variability and detectability of the features imprinted on the ULX X-ray spectra. This is crucial in order to probe and understand whether the proposed unification scenario actually works. This is further complicated as three-dimensional radiation-hydrodynamic simulations of outflows from supercritical accretion have shown that the wind is clumpy and may vary on very short (seconds) time scales (Kobayashi et al., 2018). Regardless of the exact launching mechanism (e.g. radiation pressure or magnetic pressure), it is still important to understand the thermal stability and other physical properties of ULX winds and, more generally, of winds observed in super-Eddington accretors. ULX winds, for instance, might be relevant to understand the existence of superbubbles with sizes of about $\sim 100s$ pc around many ULXs (e.g., Pakull & Mirioni 2002; Pakull et al. 2010).

This paper is therefore dedicated to the study of the stability of ULX winds and their comparison to winds found in better known X-ray sources such as Narrow Line Seyfert 1, classical sub-Eddington Seyfert 1 nuclei and classical supersoft sources such as white dwarfs burning hydrogen (some of which are novae). This paper is structured as follows. We present the X-ray objects of this study in Sect. 2 and the computation of stability curves in Sect. 3. We discuss the results and provide insights on future X-ray missions in Sect. 4 and give our conclusions in Sect. 5. Some technical detail on our analysis is reported in Appendix A.

2 The sources

In this section we report the objects of our study. We include archetypal ULXs focusing on three cornerstone states (soft, intermediate, hard) represented by NGC 5408 X-1, NGC 1313 X-1 and the pulsar NGC 300 ULX-1. All three show evidence of outflows albeit at different statistical significance owing to the different exposure times (Pinto et al. 2016, Kosec et al. 2018b). Some characteristics of these sources are shown in Table 1 along with other ULXs among the brightest ones and their highest significance wind components. These three sources have been observed with several satellites yielding broadband spectra covering the energy domain from IR to hard X-rays, which is necessary to compute photoionization balance. We also give examples of some other bright ULXs with similar spectral and timing behaviour.

On a broader context, we also want to compare ULX outflows to those found in other sources, particularly those where radiation pressure may play an important role (Table 2). Therefore, we include in our study the supersoft source nova V2491 Cygni, which exhibited a thick 3000 km/s outflowing envelope 40 days after its outburst (see, e.g.,   Ness et al. 2011 and Pinto et al. 2012) and the highly-accreting NLS1 IRAS 13224-3809 that shows an ultrafast $0.2c$ outflow responding to continuum changes in agreement with theoretical models of radiation pressure driven winds (Parker et al., 2017; Pinto et al., 2018).

We also test a blackbody spectrum with a temperature of 50 eV to mimic the spectra of supersoft ultraluminous X-ray sources (see, e.g., Urquhart & Soria 2016) and tidal disruption events (TDEs, see, e.g., Leloudas et al. 2016). A TDE is the disruption of a star that occurs near the event horizon of a ‘small’ supermassive black hole, $M_{\rm BH}\sim 10^{5-7}\,M_{\odot}$. They are believed to reach fallback luminosities $100\times$ Eddington for a short period (e.g. Wu et al. 2018) and show evidence of radiation-pressure driven winds (see, e.g., Miller et al. 2015).

For useful diagnostics, we also incorporate the sub-Eddington Seyfert type 1 AGN NGC 5548 because the spectral code that we use (SPEX, Kaastra et al. 1996) has the photoionization balance calculation for this source as default (Steenbrugge et al., 2005) and was used in previous work on ULX winds (Pinto et al. 2016, Pinto et al. 2017, Kosec et al. 2018a, Kosec et al. 2018b).

We will first compare the ULX NGC 1313 X-1 (where we found the first evidence of winds) to different types of non-ULX sources. Then we will focus on the thermal stability of ULX winds according to their ionization field.

Here we do not consider different cases of obscured / unobscured AGN as a detailed analysis was already performed by Mehdipour et al. (2016) including a description of the systematics introduced by using various photoionization codes. Testing different codes, they found variations of 10–30 % in the value of the ionization parameter of the winds and the optical depth of the absorption lines, which are difficult to detect with the current X-ray telescopes.

3 Thermal stability

The SPEX code111https://www.sron.nl/astrophysics-spex is a powerful package to compute ionization balance of high-energy astrophysical plasmas. In particular, the newly implemented pion code is optimized to perform instantaneous calculation of ionization balance in the regime of photoionization, modelling simultaneously the input continuum and the emission/absorption lines produced by the ionized gas.

Photoionization equilibrium is parametrized with the well known relationship

[TABLE]

where $\xi$ is the ionization parameter (a measure of the number of photoionizing photons per particle), $L_{\rm ion}$ the ionizing source luminosity (usually taken between 1 and 1000 Rydberg, i.e. 13.6 eV and 13.6 keV), $n_{\rm H}$ the hydrogen density and $R$ the distance between the plasma and the ionizing source (see, e.g., Krolik et al. 1981). A crucial input to the photoionization calculation is therefore the ionizing continuum or the broadband spectral energy distribution (SED). The ionization balance also slightly depends on element abundances; here we adopt the recommended proto-Solar abundances of Lodders & Palme (2009), which are also the default in SPEX.

3.1 Spectral energy distribution

In order to build the SEDs of our sources we extract fluxes through all the available data from optical to hard X-ray energies. In Fig. 3 (top) we compare the SED of ULX NGC 1313 X-1 with the SEDs of the other non-ULXs sources detailed in Table 2. All SEDs have been corrected for interstellar absorption and E(B-V) reddening.

Briefly, we take the time-averaged XMM-Newton and NuSTAR X-ray spectra of IRAS 13224-3809 from Parker et al. (2017) and the optical/UV data from the optical monitor (OM) time-average fluxes (Buisson et al. 2018). For NGC 5548 we use the default unobscured SED used by SPEX (Steenbrugge et al., 2005), while for nova V2491 Cygni we adopt the SED produced in Pinto et al. (2012).

For NGC 1313 X-1 we adopt the time-average XMM-Newton/EPIC spectrum (0.3–10 keV) where the wind was discovered (Pinto et al., 2016) along with the hard X-ray fluxes provided by NuSTAR (3–30 keV, Bachetti et al. 2013). Optical and UV fluxes are taken from Yang et al. (2011), whilst FUV fluxes are extracted from the XMM-Newton/OM using a 4 arc seconds circle onto its X-ray centroid and cross-checking with the source list produced by the XMM-SAS task OM CHAIN. The effects due to the time variability of the SED and the optical/UV fluxes of NGC 1313 X-1 are explored in Sect. 3.4–3.5.

For the other two ULXs, namely NGC 5408 X-1 and NGC 300 ULX-1 we use the X-ray spectra from Pinto et al. (2016) and Kosec et al. (2018b), respectively (see also Table 3). Optical and UV fluxes are obtained from HST measurements published in the literature (Grisé et al. 2012, Lau et al. 2016, Villar et al. 2016) and from the optical monitor on board XMM-Newton where some literature data points are missing. We prefer data from HST observations which are deeper and have higher spatial resolution than XMM-Newton/OM. Overall, the SEDs of the three ULXs look rather smooth and do not show sharp jumps of flux when compared to the other sources with more accurate measurements outside the X-ray band like AGN and novae, see Fig. 3 and 4. They resemble the broadened disc spectra in Sutton et al. (2013), which suggests that the fluxes should be accurate with uncertainties less than an order of magnitude, small enough to avoid dramatic uncertainties in the ionization balance calculations.

The X-ray portions of the SEDs for the three ULXs under investigation (NGC 1313 X-1, NGC 5408 X-1 and NGC 300 ULX-1) are determined by fitting the time-averaged XMM-Newton/pn spectrum with a multi-component model consisting of blackbody, multicolour blackbody and powerlaw as discussed in Sect. 1 (see also Kosec et al. 2018b). The blackbody component ($T\sim 0.15-0.25$ keV) fits the disc upper photosphere, the multicolour blackbody ($T\sim 1-3$ keV) reproduces the overall disc emission (and/or any contribution from the boundary layer near a neutron star) and the powerlaw ($\Gamma\sim 1.5-2.5$) accounts for emission from the inner disc (and eventual contribution from the accretion column onto a neutron star, see e.g. Walton et al. 2018a). Neutral absorption from the Galactic interstellar medium and circumstellar medium near the ULXs is reproduced by the $hot$ model in SPEX adopting a low temperature $T=0.5$ eV (see, e.g., Pinto et al. 2013 and references therein).

In summary, we have obtained the UV-to-IR fluxes at 2120, 2310, 2910, 3600, 4300 and 5500 and 9000 Å from the literature or the XMM-Newton/OM, the X-ray fluxes between 0.3–10 keV ($\sim 1.2-40$ Å) from the XMM-Newton RGS+pn spectral fits. The SED portions above 10 keV are either provided by NuSTAR or extrapolated from our XMM-Newton spectral model. The SED range covering the FUV domain is missing due to interstellar absorption and is therefore interpolated. This adds uncertainty in the photoionization calculation but we stress that our SED shapes agree with those calculated for irradiated discs of black holes in the regime of super-Eddington accretion (Ambrosi & Zampieri, 2018). We investigate some systematic effects due to uncertainties in the low-energy flux in Sect. 3.5.

3.2 NGC 1313 ULX and non-ULX objects

We compute the photoionization balance running the SPEX pion code version 3.05.00 onto all SEDs. In Fig. 3 (middle panel) we show the detail of the photoionization balance computation (i.e. the temperature $T$ – ionization parameter $\xi$ curves). Some differences can be distinguished between these curves such as the variable slopes at intermediate values of the ionization parameter ($\log\,\xi\sim 2-3$). The detail of how heating and cooling rates are computed for ULX NGC 1313 X-1 and, in particular, their trends with the ionization parameter can be found in Appendix A.2.

However, a more informative diagram is provided by the stability curves. Those can be produced by computing the ratio between the radiation pressure ($F/c$) and the thermal pressure ($n_{\rm H}kT$), which is given by $\Xi=F/n_{\rm H}ckT=19222\,\xi/T$, with $F=L_{\rm ion}/4\pi r^{2}$ (Krolik et al., 1981). The stability (or $S$) curves computed for our sample of sources are shown in Fig. 3 (bottom panel). The $S$ curves are powerful diagnosis tools because along them heating equals cooling and, therefore, the gas is in thermal balance.

On the left side of the $S$ curve cooling dominates over heating, whilst on the right side, heating dominates over cooling. Most importantly, where the $S$ curve has a positive gradient, the photoionized gas is thermally stable or, in other words, small perturbations upwards in temperature will be balanced by an increase in cooling. Similarly, small perturbations downwards will be balanced by increase in heating. Instead, in the points on the $S$ curve with a negative gradient, the plasma is thermally unstable and therefore any perturbation upwards (downwards) will cause further increase (decrease) of temperature with a dramatic change, if not loss, of the equilibrium. In principle, any two components on the $S$-curve with the same $\Xi$ value are expected to be in pressure equilibrium.

The stability curves computed for the different sources (Fig. 3) show notable differences, with the harder sources exhibiting several unstable branches. Here with hard sources we actually mean ‘X-ray hard sources’ i.e. those sources which have a high ratio between the fluxes calculated in the 2–10 keV and 0.3–2 keV energy bands, respectively. High-energy photons from the Fe K energy band and beyond significantly affect the ionization balance by increasing the size of the unstable region where, in principle, no plasma should be found (see also Krolik et al. 1981). On the other hand, the softest sources like novae or SSUL have rather stable curves. Moreover, the high-temperature turnovers differ by up to two orders of magnitude (Fig. 3, middle panel), which indicates great differences in the Compton temperature. Such variations are expected in objects of different accretion rates and are relevant for thermal and thermal-radiative winds since they can be suppressed if the Compton temperatures are too low (Done et al. 2018).

$S$ curves of AGN (and classical X-ray binaries) have been studied in great detail and their instability branches are ubiquitously found (see, e.g., Krolik et al. 1981; Reynolds & Fabian 1995; Netzer et al. 2003; Steenbrugge et al. 2005; Chakravorty et al. 2016; Mehdipour et al. 2016; Higginbottom et al. 2017). Of course, owing to its SED, the stability curve of archetypal ULX NGC 1313 X-1 fits in between those of AGN and soft sources, showing only some small branches of instability.

3.3 NGC 1313 ULX and other ULXs

A comparison between the $S$ curves of the intermediate-state ULX NGC 1313 X-1 with harder (Seyfert 1 AGN) and softer (supersoft source, SSS, nova) spectra provides a crude idea of how thermal instability increases for the wind parts which are subjected to harder SEDs. There is currently a strong likelihood that ULXs are characterized by geometrically-thick discs and optically-thick outflows, at least for ULXs powered by black holes or neutron stars with limited magnetic fields (e.g. $\lesssim 10^{12}G$). It is therefore thought that the spectral hardness of a ULX is due to a combination of accretion rate and inclination. In fact, at lower inclinations (face on) we should see harder spectra (e.g., Middleton et al. 2015; Feng et al. 2016; Urquhart & Soria 2016; Pinto et al. 2017).

It is very likely that the inner portions of the wind are exposed to a harder photoionizing field compared to the outer regions or, in other words, the inner regions may see beamed hard X-rays, while the radiation field in the outer regions may be dominated by UV and soft X-rays (see Fig. 2). If this scenario applies then the gas phases in these regions will have radically different stability curves.

In Fig. 4 (top) we compare the SEDs of three ULXs with different X-ray spectral shapes: NGC 300 ULX-1 (hard), NGC 1313 X-1 (intermediate) and NGC 5408 X-1 (soft, see also Table 1). We have computed the ionization balance for the other two ULXs with the same code used for NGC 1313 X-1 in Sect. 3.2 and found significant differences in their stability curves (see Fig. 4, bottom). As expected, the harder ULX has longer branches of instability. This means that, according to the ULX unification scenario, there should be a systematic increase of instability for wind portions towards the inner regions of the accretion disc.

In Fig. 4 (bottom) we also indicate the ionization parameters and the corresponding pressure ratios of the winds detected in these ULXs on each $S$ curve. For NGC 1313 X-1, there are two regions due to two alternative solutions for the Fe L and Fe K energy bands (Pinto et al. 2016; Walton et al. 2016, respectively). It appears that the wind ionization state increases with the spectral hardness of the source (and likely decreases with the viewing angle) and that the measurements are formally acceptable as none falls along instability branches. For NGC 5408 X-1, instead, the outflow has a more complex structure with two components of comparable ionization state ($\log\xi\sim 1.5-1.9$ and $1.7-2.1$) but different velocities (Pinto et al. 2016; see also Table 4).

3.4 NGC 1313 ULX: Spectral Variability

We normally consider as classical ULXs those sources which surpass the Eddington limit of a 10 Solar mass black hole for long periods of time. However, these persistent beacons can show remarkable X-ray variability on time scales from hours to several days. PULXs for instance show long-term variability associated to either precession (see, e.g., Middleton et al. 2018) or switching between two different regimes that may be interpreted as a propeller phase (see, e.g., Tsygankov et al. 2016). Luminosities of ULXs can vary by an order of magnitude, followed by changes in the spectral slope by 20% or more (see, e.g., Kajava & Poutanen 2009). It is therefore a useful exercise to compute the ionization balance for ULXs where the spectral shape significantly changes with time in order to probe any effects on the wind stability (or on the portion of the wind in our line-of-sight).

Of course, NGC 1313 X-1 provides an excellent workbench thanks to the dozens of dedicated observations taken with XMM-Newton since early 2000s and owing to its high spectral variability (see, e.g., Middleton et al. 2015). We therefore focus on three characteristic spectra of this ULX which show that the source alternates between soft-ultraluminous, hard-ultraluminous and bright-broadened-disc states, respectively (according to the classification of Sutton et al. 2013; for more detail see Table 3).

NGC 1313 X-1 normally spends most of the time in the hard-ultraluminous state with intermediate brightness compared to the soft (X-ray faint) and the broadened-disc (X-ray bright) states. The corresponding SEDs are shown in Fig. 5 (top). For simplicity, we adopt the same values of optical and UV flux as given by the time-averaged SED (see Sect. 3.5 for a discussion on effects due to variable optical/UV fluxes). These broadband spectra represent a small sample, but they provide a good indication of how much the ionization field may vary and, therefore, how different will be the SED seen by the outflowing gas depending on its location (and our line of sight).

The standard continuum model consisting of absorbed powerlaw, blackbody and modified blackbody overestimates the flux below 1 keV due to the powerlaw steep index ($\Gamma\sim 2$). This effect is stronger for the X-ray bright and faint spectra (BBD and SUL, see Fig. 5 top) due to their strong spectral curvature and requirement of high column density (an order of magnitude above the interstellar value, $N_{\rm H,ISM}\sim 5\times 10^{20}$ cm*-2*). For this reason, we cut their extrapolation below 1 keV forcing a flat plateau (see the blue dashed and red dashed-dotted lines). Alternatively, we also test another model where the powerlaw is substituted by another modified blackbody, i.e. $hot*(bb+mbb+mbb)$. This model reproduces the spectral curvature without the need of high neutral absorption (see the blue dashed-dotted and red dotted lines) and indicates the order of magnitude of uncertainty in the (unabsorbed) SEDs of these states. The intermediate, hard-ultraluminous state is not strongly affected by this issue and, therefore, we just keep the results obtained with the standard model.

We compute the ionization balance for these SEDs of NGC 1313 X-1 as done for the time-average spectrum and show the stability curves in Fig. 5 (bottom). As expected, the $S$ curves significantly differ; in particular the broadened-disc SED exhibits a large instability branch for ionization parameters $\log\xi\sim 1.8-2.3$, which is smaller in the intermediate-hard state. The soft state shows instead a generally stable curve very similar to soft ULXs, novae and other SSS (see Fig. 3 and 4). The uncertainties in the SEDs due to the soft X-ray extrapolation do not seem to have drastic effects on the thermal stability of the wind.

3.5 NGC 1313 ULX and optical/UV screening

It is difficult to build a proper spectral energy distribution without simultaneous observations in all energy domains going from the optical to the X-rays. In most cases quasi-simultaneity with optical/UV data taken at slightly different time from the X-ray can provide reasonable measurements as we might expect the optical/UV not to vary on time scales as short as X-rays owing to the larger emitting region. In ULXs, of course, the scenario is complex as the bulk of the UV (and possibly the optical) emission may come from the outer disc itself but some contamination can be expected from the companion star (e.g., Grisé et al. 2012; Ambrosi & Zampieri 2018). Sutton et al. (2014) have shown that most objects at $\sim 10^{39}$ erg/s exhibit disc reprocessing fractions similar to sub-Eddington objects, but at higher luminosities X-ray reprocessing seems to increase, possibly due to scattering of wind X-rays back to the outer disc.

It is important to notice that the full optical-UV-X-ray SED that we observe might not actually be representative of the SED that ionizes all the different regions of the wind. The thick inner disc and even the wind itself may self-screen the inner regions of the wind from some of the optical/UV emission arising from the outer disc and the companion star.

Fig. 2 shows a simplified, qualitative, description of a radiation-driven wind from a super-Eddington accretion disc around a stellar-mass black hole or a neutron star. The temperature at the spherization radius can be roughly expressed as $T_{\rm sph}=1.5\,m^{-1/4}\,\dot{m}_{0}^{-1/2}\,(1+0.3\,\dot{m}_{0}^{-3/4})\,\rm keV$ (see, e.g., Poutanen et al. 2007), which is of the order of 0.2 keV for a black hole with a mass $M=m\,M_{\odot}=10\,M_{\odot}$ and an accretion rate $\dot{M}=15\,\dot{M}_{\rm Edd}=\dot{m}_{0}\,\dot{M}_{\rm Edd}$ (where $\dot{m}_{0}=\dot{m}$ at $R>R_{\rm sph}$ and we adopt the common definition $\dot{M}_{\rm Edd}=L_{\rm Edd}/\epsilon c^{2}=1.86\cdot 10^{18}(M/M_{\odot}){\rm\,g\,s}^{-1}$ with a radiative efficiency $\epsilon=8$%). A comparable temperature is found in the soft X-ray component of almost all ULXs. Most X-rays will be produced within $T_{\rm sph}$, while at larger radii the flux will peak in the UV and optical energy bands.

The uncertainties in the SED shape and, particularly, in the UV/optical irradiation of the winds can be tested by computing the ionization balance for SEDs where the optical/UV fluxes are 1-2 orders of magnitude below the threshold measured with the XMM-Newton/OM or other facilities. We therefore construct ad-hoc SEDs with low-energy attenuation where only 1% or 10% of the optical/UV fluxes are transmitted and therefore received by the plasma in our LOS (see Fig. 6, top). This is a very simplified way to simulate partial covering of outer UV/optical photons by inner optically thick clouds. We use as template the SED of NGC 1313 X-1 during the most common hard-ultraluminous state where a wind was detected.

The stability curves show a progressive extension of the unstable branch for optical/UV transmission decreasing from 100% to 10% and then 1% (see Fig. 6, bottom). The effects of optical/UV screening is however more pronounced at low ionization parameters ($\log\xi\lesssim 1.5$) which might not affect the current wind detection ($\log\xi\gtrsim 2.2$).

4 Discussion

The most common AGN rarely reach the Eddington limit. Other extreme objects such as the tidal disruption events may reach fallback luminosities $100\times$ Eddington for a short period (e.g. Wu et al. 2018), but they are rare and difficult to observe due to their transient nature. It is also challenging to construct robust broadband SEDs for many AGN and TDEs because of extinction from neutral hydrogen, which mainly affects the energy range containing the SED peaks.

ULXs instead are X-ray bright and persistent, surpassing the Eddington limit for long periods of time. The search and study of radiatively driven winds in these sources may therefore provide clues on the way primordial black holes grew up when the Universe was a few hundred million years old. In this paper we provide a first attempt to understand the thermal state of ULX winds, their relationship with the LOS (that can change due to precession, e.g., Middleton et al. 2018) and the accretion rate (e.g., Poutanen et al. 2007; Middleton et al. 2015).

4.1 Thermal stability of winds in ULXs

In Sect. 3.1 we have shown broadband SED spectra of different types of sources (AGN, ULX, nova, etc). Then we have computed the ionization balance of plasma around ULXs for different spectral states and shown the stability curves, indicating combinations of $\xi-\Xi$ where the plasma is stable (Fig. 4). The winds detected in the three ULXs under investigation (i.e. NGC 300 ULX-1, NGC 1313 X-1 and NGC 5408 X-1) show $\xi-\Xi$ pairs corresponding to positive slopes of the $S$ curve, which suggests that they are in stable equilibrium (as it would be required for them to be detectable).

It is important to notice that our results predict the existence of an evacuated funnel, similar to that one previously invoked in ULXs (e.g. Middleton et al. 2015). In the innermost regions near the compact object, the wind likely becomes highly-ionized and/or thermally unstable, possibly providing a better view on the hard X-ray emitting region (see also Sec. 4.3 for evidence from observations). Moreover, density variations might cause thermal instabilities and multiphase gas (see, e.g., Owocki & Cohen 2006) as suggested by the relatively mild ionization parameters and complex wind structure found in some ULXs.

4.2 ULX winds in a broader context

4.2.1 Comparison with winds from AGN and XRBs

There have been several studies dedicated to outflows in AGN and classical X-ray binaries, as well as their thermal stability in the past decades. Winds are better characterized in these objects owing to their high brightness due to either proximity (Galactic X-ray binaries) or huge luminosity (AGN). Below we report a brief summary on such winds.

AGN show a variety of winds detected in the UV such as the $0-5000$ km/s warm absorbers (see, e.g., Reynolds & Fabian 1995; Steenbrugge et al. 2005; Laha et al. 2014), the UV/X-ray obscurers that create temporary occultations and have comparable velocities (see, e.g., Kaastra et al. 2014; Mehdipour et al. 2017; Turner et al. 2018), the $\sim 10\,000$ km/s broad-absorption-line quasars (BALQSO, e.g., Proga et al. 2000; Hewett & Foltz 2003; Dai et al. 2008; Leighly et al. 2014) and the most extreme $0.1-0.3c$ ultrafast outflows detected in the X-ray band (see, e.g., Pounds et al. 2003; Reeves et al. 2003; Tombesi et al. 2010). The highly ionized ultrafast outflows tend to occupy large branches of stability through the $S$ curves (see, e.g., Danehkar et al. 2018; Kraemer et al. 2018), while the warm absorbers and the other winds with intermediate ionization states can fall near the regions of thermal instability ($\log\xi\sim 1-2$, see e.g. Krolik et al. 1981; Krongold et al. 2003; Netzer et al. 2003; Detmers et al. 2011). The ultrafast outflows found in AGN are likely the most relevant comparison for the winds seen in ULXs as their similar properties ($\xi$, $v_{\rm outflow},N_{\rm H}$) suggest some analogies in their launching mechanisms.

X-ray binaries often show highly ionized winds detected primarily in the Fe K band (see, e.g., Miller et al. 2004, 2006; Neilsen & Lee 2009; Ponti et al. 2012). These winds are typically slow (a few 100s km/s, therefore associated with the disc); they are often detected in the soft states and disappear when the source is in a hard state, most likely due to the instability branch that stretches further at harder ionizing SEDs (see, e.g., Chakravorty et al. 2016; Higginbottom et al. 2017). However, an acceleration of the flow or the exhaustion of the plasma during the soft state are possible, viable solutions (see, e.g., Gatuzz et al. 2019).

SEDs of ULXs can significantly differ the one from the other, but they often show shapes resembling that of a broadened disc blackbody or of a complex, super-Eddington, accretion disc (see, e.g., Gladstone et al. 2009; Sutton et al. 2013; Middleton et al. 2015). Their stability curves are halfway between those of high-Eddington NLS1 and supersoft novae, having only small branches of instabilities (see Fig. 3). The hard state ULXs such as the ultraluminous pulsar NGC 300 ULX-1 have in proportion a larger fraction of hard X-ray photons which further extend the unstable section (see Fig. 4). For PULXs with stronger magnetic fields, we can naturally expect more significant truncation of the inner accretion disc and strong polar inflow (at a given $\dot{M}$). This not only decreases radiation pressure in the inner regions of the disc, but can also produce a harder SED, followed by increasing instabilities in the inner portions of the wind. In other words, strongly-magnetized PULX will likely have weaker winds (i.e. with lower outflow rates) than weakly-magnetized PULXs for the same values of $\dot{M}$.

As mentioned before, outflows can also be launched by mild magnetic fields (see, e.g., Romanova & Owocki 2015; Romanova et al. 2018; Parfrey et al. 2017). A study of the launching mechanism is however beyond the scope of our paper.

4.2.2 Effects of ULX winds on their surrounding

This study of ULX winds at small scales might be relevant to understand their effects at larger scales such as their role in the existence of ULX superbubbles with sizes of about $\sim 100s$ pc that are thought to be similar to the (radio) W 50 nebula around SS 433 (e.g., Pakull & Mirioni 2002; Pakull et al. 2010). For instance, Pakull et al. (2006) have shown that the bubbles have supersonic expansion speeds of 80–250 km/s, derived from the width of H $\alpha$, which suggests that (at least some parts of) the bubbles are shock-excited by winds or jets rather than photoionized by X-rays from the central ULX (see also Siwek et al. 2017). They estimated that about $10^{39-40}$ erg/s is required to ionize the optical bubbles. This is comparable to the estimates of wind kinetic power obtained in ULXs so far (Pinto et al., 2016; Pinto et al., 2017; Walton et al., 2016; Kosec et al., 2018a; Kosec et al., 2018b):

[TABLE]

where $\dot{M}_{w}=4\,\pi\,R^{2}\,\rho\,v_{w}^{2}\,\Omega\,C$ is the outflow rate, $\Omega$ and $C$ are the solid angle and the volume filling factor (or clumpiness), respectively, $\rho$ is the density and $R$ is the distance from the ionizing source. Here we have used Eq. (1) to get rid of the $R^{2}\rho$ factor where $\rho=n_{\rm H}\,m_{p}\,\mu$ with $m_{p}$ the proton mass and $\mu=0.6$ the average particle weight of a highly ionized plasma. A solid angle $\Omega/4\pi=0.3$ and volume filling factor $C=0.3$ were chosen as fiducial average values from MHD simulations of winds driven by radiation pressure in super-Eddington winds (Takeuchi et al., 2013).

The filling factor of the wind might actually be much smaller, at least for the cool and outer phases. Using Eq. (23) in Kobayashi et al. (2018) and assuming that the outflow is comparable to the accretion rate, we obtain $C\sim 5\times 10^{-3}$, in which case the wind power would be much lower, although still significantly affecting the surrounding medium.

Following the approach of Pakull et al. (2006), we know that the maximum mass ejected in the wind, $\dot{M}_{w}\,\tau$, throughout the ULX life time $(\tau)$ cannot be larger than the total mass transferred from the companion star to the Roche lobe of the compact object, which can be expressed as $\delta m\,M_{\odot}$, where $\delta m$ is a factor of a few (for a standard high-mass companion star). From the wind kinetic luminosity $L_{w}$ in Eq. (2) we can derive the wind velocity as

[TABLE]

where $L_{39}$ is the wind kinetic luminosity in units of $10^{39}$ erg/s, $\tau_{6}$ is the ULX life time in $10^{6}$ yr and $\delta m$ is the mass in solar mass units. This broadly agrees with the typical ULX wind velocities detected so far.

Alternatively, following the approach of Castor et al. (1975); Weaver et al. (1977); Soker (2004); Begelman et al. (2006), we can determine the expected radius of a bubble that has been inflated by a ULX wind. Assuming energy conservation, the bubble radius can be approximated as

[TABLE]

where $\rho_{\rm\footnotesize{\,ISM}}\sim 10^{-24}$ g cm*-3* is the interstellar density. If we adopt a wind kinetic power of a few $10^{39}$ erg/s, we obtain bubble with sizes of about 100 pc, which is consistent with the measurements quoted above.

These results have to be taken with some caution due to several uncertainties in the properties of both the winds and the bubbles. For instance, we adopt a lifetime of $10^{6}$ yr from the size $\sim 100$ pc, the current expansion velocities of $\sim 100$ km/s and the expansion low $R\sim t^{3/5}$, but it is possible that in the past the bubbles were expanding faster. This would correspond to shorter lifetimes. We have adopted solid angles of $0.3$, according to MHD simulations, but Begelman et al. (2006) suggest slightly lower values ($\Omega/4\pi=0.1$). Finally, the total mass transferred from the companion to the compact object might be smaller than a solar mass ($\delta m\,M_{\odot}<1\,M_{\odot}$). However, these uncertainties would tend to cancel each other out in the equations above. Our predictions of radii for the wind-inflated bubbles and of wind velocities are likely correct within the order of magnitude, which makes the comparisons meaningful and therefore still likely provide useful points of comparison. This would provide further evidence on the detection of winds in ULXs and their role on the surrounding environment.

Currently, the ULX wind sample is very small due to instrumental limitations and short exposure times. Once a larger sample would be available it will be of high interest to compare the presence of winds in ULXs with that of superbubbles around them – and their characteristics – in order to better understand their formation.

4.3 Effects from inclination and accretion rate on the appearance of ULX winds

Fig. 2 shows an artistic impression of a wind launched from an accretion disc due to high radiation pressure. As mentioned before, this simplified geometry would apply to the cases where the compact object is a black hole or a non-magnetic neutron star (with the maximum magnetic field depending on the accretion rate). The wind gives the system a funnel geometry and scatters, i.e. geometrically beams, the X-ray photons coming from the inner regions (see, e.g., Poutanen et al. 2007; Middleton & King 2016). Therefore, we expect to see a hard X-ray spectrum followed by high ionization in the wind when the source is viewed at low inclinations (i.e. on-axis) and the more the inclination increases, the softer the X-ray spectrum and the lower the ionization of the wind should be, for a given accretion rate (e.g., Kawashima et al. 2012; Pinto et al. 2017).

We show the trend between the velocity and the ionization parameter of the winds with the spectral hardness in Fig. 7. We notice that a third of these detections were marginal ($\lesssim 3\sigma$, see those shown as diamonds) and that the results have to be taken with caution. We retrieve all the necessary information on the wind properties from previously reported detections and double-check the results with the newer SPEX version used in this paper. More detail on the wind data is reported in Sect. A.3 and, particularly, in Table 4. We find a possible correlation between the ULX spectral hardness, $HR$, the wind ionization parameter, $\xi$, and velocity, $v_{w}$. We fit a straight line to the data using the PYTHON routine odr222https://docs.scipy.org/doc/scipy/reference/odr.html that performs orthogonal distance regressions and obtained the following best fits:

[TABLE]

where $HR$ is the hardness ratio. The Spearman and Pearson correlation coefficients range between 0.6 and 0.8 (see Table 5), which indicates a substantial positive correlation among the three parameters. Although there is significant scatter, the apparent trend is qualitatively consistent with the scenario for which we see hotter and faster wind component at lower inclination. The velocity trend may suggest that the optically-thin component of the wind (which is responsible for the blueshifted absorption lines) is likely not moving along the equatorial plane or, more simply, that at larger inclinations we see material expelled from an outer region where the escape velocity is naturally lower. The ULX wind sample is currently too small to obtain conclusive results, but our measurements agree with the general picture.

Some of the scatter is likely due to systematics. Here we show the $(HR,\xi,v_{w})$ values for the observations with statistics good enough to enable wind detections, but ULXs vary in hardness ratio. Therefore, we also show the $HR$ fractional variability in Fig. 7 (dotted horizontal lines) as retrieved from the literature (Middleton et al. 2015; Kosec et al. 2018a; Kosec et al. 2018b; Pinto et al. 2017; Carpano et al. 2018 and references therein). Although, it might not be physically correct to fit a model to these points, the $HR$ variability could still explain why they deviate from a tight correlation.

For the several reasons discussed above, the hard NGC 300 ULX-1, the intermediate NGC 1313 X-1 and the soft NGC 5408 X-1 ULXs are potentially being seen from increasing viewing angles. The comparison of their thermal stability curves in Fig. 4 would then confirm that the inner portions of the wind have to be necessarily more ionized and located along the $\xi-\Xi$ Compton cooling arm or they would disappear due to strong instabilities as often observed in the hard states of Galactic X-ray binaries (Chakravorty et al., 2016).

Super-orbital periods have been reported in several ULXs (see, e.g., Kaaret et al. 2006; Motch et al. 2014; Walton et al. 2016). These may be evidence of precession which suggests that different portions of the accretion disc and the wind would be visible to us at different times (see, e.g., Luangtip et al. 2016; Middleton et al. 2018). Variability in the accretion rate can cause similar behaviour because the opening angle of the wind is expected to be related to the accretion rate (Poutanen et al. 2007). ULXs that have hard spectra at moderately-high accretion rates (dashed line in Fig. 5, top) could potentially undergo substantial obscuration of the inner region and appear softer after a large increase of accretion rate (dotted line). In principle, this should also strengthen the soft emission due to the outer disc emission (dotted line). At even higher rates or inclinations the ULX would look like a supersoft source (solid line).

In summary, either precession or variable accretion could expose our LOS to a wind portion whose ionizing field may be dominated by a region of the accretion disc with different average temperature. In Sect. 3.4, we have used the spectral variations of NGC 1313 X-1 to understand the effects of dramatic changes in the SED. When the source brightens from an intermediate-hard (HUL) state to the bright-broadened-disc state (BBD) and reaches its X-ray luminosity peak, the $S$ curve changes, extending the region of instability and, therefore, affects the wind “warm” component with low ionization parameter ($\log\xi\sim 1.8-2.3$, dotted line in Fig. 5, bottom). This means that this wind component is either not in equilibrium or simply out of our LOS, in each case basically undetectable. This can explain the weakening of the 1 keV residuals previously seen in Middleton et al. (2015) and seems to be confirmed by new, deeper, XMM-Newton observations of NGC 1313 X-1 (Pinto et al., 2019).

4.4 Caveats

It is worth noticing that the black holes powering AGN - which accrete at much lower Eddington ratio than those in ULXs - exhibit ultrafast outflows within a broad range of inclination angles and their wind parameters do not show strong correlations with disc characteristics such as the inclination (see, e.g., Tombesi et al. 2014). This would suggest that AGN winds might differ from those seen in highly super-Eddington sources like ULXs, possibly due to a different launching mechanism. In fact, low-Eddington AGN winds are likely driven by magnetic pressure (see, e.g., Fukumura et al. 2017 and references therein).

We also address briefly the possibility that some portions of the wind can be launched by thermal heating close to the Compton temperature (see, e.g., Begelman et al. 1983; Done et al. 2018). In our case, the Compton radius will be $R_{IC}=GM\mu m_{\rm p}/(kT_{IC})\sim 3\times 10^{5}R_{S}$ for a $10\,M_{\odot}$ black hole and $T_{IC}\sim 10^{7}$ K (i.e. the temperature where the $T-\xi$ curve flattens for the ULX in Fig. 3, middle panel). This is much larger than the spherization radius $R_{\rm\,sph}=5/3\dot{M}k/(8\pi c)\sim 50R_{S}$ and the photospheric radius $R_{\rm\,phot}=\dot{M}_{w}k/(4\pi v_{w}\cos\theta)\sim 500R_{S}$, the latter indicating where the wind becomes optically thin. Here we adopt $\dot{M}=10\dot{M}_{\rm Edd}$, $\dot{M}_{w}=\dot{M}$, $v_{w}=0.2c$, $\theta=45^{\circ}$ and $\dot{m}=10$ (see, e.g., Poutanen et al. 2007; Fiacconi et al. 2017). The hard X-rays from the inner region will likely be screened by the thick disc atmosphere or the wind itself, and it will be difficult for them to reach the putative Compton radius (just like qualitatively shown in Fig. 2). Therefore, we do not expect to see a classical Compton wind with slow velocities ($\sim 500$ km/s) from the outer disc.

A word of caution is required regarding the velocity of the wind for magnetic accretors. The current estimates suggest some winds as fast as $0.2c$ (see e.g. Fig. 7). If this corresponds to the escape velocity, for a neutron star with $B\sim 10^{12}$ G, it would mean that the wind has to be launched from inside the magnetosphere ($R_{M}\sim 100R_{S}$), with possible additional contributions from the magnetic fields themselves. This would not be the case for much lower magnetic fields, e.g. $B\lesssim 10^{10}$ G, or if the detected high velocities are due to later vertical acceleration (e.g. Takeuchi et al. 2013). Further observations and comparisons of winds in pulsating vs. non-pulsating ULXs might help to understand the phenomenology.

Another caveat for the applicability of our model and calculations of thermal stability concerns the recombination timescales of the plasma. If such timescales are much longer than the timescale of the expansion necessary to drastically change the ionizing field then the gas might be out of equilibrium. The recombination timescales depend on the specific ion, they tend to be longer for higher ionization states. Through the pion model in SPEX it is possible to calculate the product between the electron density $n_{\rm e}$ and the recombination time $t_{\rm rec}$ for a given SED. Using the time-averaged SED of NGC 1313 X-1 (see Fig. 3, top panel) and adopting $\log\xi=2.3$ from the best fits (see Table 4 and Pinto et al. 2016), we estimate $n_{\rm e}\cdot t_{\rm rec}=n_{i}/(n_{i+1}\alpha_{i+1}-n_{i}\alpha_{i})\sim 10^{9-10}{\rm cm}^{-3}{\rm s}$, where $n_{i}$ is the density of the ion in the state $i$ and $\alpha_{i+1}$ is the recombination rate coefficient (in ${\rm cm}^{3}{\rm s}^{-1}$) from state $i+1$ to $i$. The range indicates the values calculated for the most relevant transitions of the following ions: O vii-viii and Ne ix-x. Pinto et al. (2019) placed a constraint on the density of the gas responsible for the emission through the O vii triplet: $n_{\rm H}\sim 10^{10-12}{\rm cm}^{-3}$ (with $10^{11}{\rm cm}^{-3}$ as a median value) similar to disc coronae in Galactic X-ray binaries. The outflowing gas coming from the inner regions might be at even higher densities. This would suggest $t_{\rm rec}\sim 0.001-1{\rm s}$ or shorter (with a median value around $0.05{\rm s}$). On average we would expect the gas to have travelled around 2 $R_{\rm\,sph}$ or 0.2 $R_{\rm\,phot}$ during the recombination time scale (for $v=0.2c$, $\dot{M}=10\dot{M}_{\rm Edd}$ and $M_{\rm BH}=10\,M_{\odot}$). In the most pessimistic case it would have travelled a few $R_{\rm\,phot}$, where the photospheric radius is the region with optical depth close to unity. We expect a similar condition for a neutron star with an accretion rate of $100\dot{M}_{\rm Edd}$. This means that there can be enough time for the outflowing gas to recombine before dramatical changes would happen in the ionizing field.

In this paper we have shown a first step towards the characterization and the understanding of winds in super-Eddington accretion discs of stellar-mass black holes. In the future it will be important to further investigate systematic effects due to relativistic corrections due to the velocity of the wind and to time scales of its expansion as compared to the recombination and cooling time scales.

4.5 Future missions: XRISM and ATHENA

XMM-Newton and Chandra have provided strong insights in the knowledge of ULXs and accretion physics in general since their launches 20 years ago. However, newly planned missions will significantly improve our ability to detect winds in these objects. XRISM333X-Ray Imaging and Spectroscopy Mission, a.k.a XRISM, to be launched in 2022, see https://heasarc.gsfc.nasa.gov/docs/xrism., for instance, will enable for the first time observations of ULXs with a good collecting area and an outstanding spectral resolution of about 5 eV in the hard X-ray band ($\sim 2-10$ keV), a combination missing in current telescopes (Guainazzi & Tashiro, 2018).

Here we perform a simulation of NGC 1313 X-1 as it may be observed by XRISM using its predecessor Hitomi/SXS 5 eV response matrix (see Fig. 8, bottom panel). We adopt the best fit model from Pinto et al. (2016) consisting of a spectral continuum (blackbody $T=0.3$ keV and powerlaw $\Gamma=1.9$ components), an isothermal rest-frame collisionally-ionized plasma with a temperature of 0.8 keV, an outflowing photoionized plasma (column density $N_{\rm H}=5\times 10^{21}{\rm cm}^{-2}$, $\log\xi=2.3$, $v_{\rm LOS}=0.2c$, $v_{\rm turb}=500$ km s*-1*), and Galactic absorption through the $hot$ model in SPEX with $N_{\rm H}=2\times 10^{21}{\rm cm}^{-2}$ (see, e.g., Pinto et al. 2013). The simulation performed with the XRISM/Resolve microcalorimeter provides a clear view of the Fe L region around 1 keV and boosts the detection of lines above 2 keV with the resolution necessary to break degeneracy such as between the ionization parameter and the outflow velocity. However, we will still need observations of 100 ks or above in order to get $5\,\sigma$ significance detections for the strongest lines. Below $\sim 0.7$ keV the XMM-Newton/RGS grating spectrometer will still be the most sensitive instrument.

X-ray astronomy will radically change in the 2030s after the launch of ATHENA (Nandra et al. 2013; Guainazzi & Tashiro 2018). The mission will have both a wide field imager (WFI) with collecting area an order of magnitude above any current X-ray mission over the $0.3-10$ keV energy band and a field of view of 40’$\times$40’, along with the X-ray Integral Field Unit (X-IFU). The X-IFU will be the first X-ray microcalorimeter to provide both high spatial (5”) and spectral (2.5 eV) resolution with an effective area an order of magnitude higher than any current high-spectral resolution X-ray detector. This corresponds to an improvement of almost 2 orders of magnitude in the way we detect and resolve narrow features in the canonical 0.3–10 keV X-ray energy band.

We perform another simulation of NGC 1313 X-1 using the same baseline model but only assuming an exposure time of 10 ks and the response matrix of ATHENA/X-IFU (dated late 2018). The expected results are extremely promising (Fig. 8, top). Even a short exposure with X-IFU will provide a high signal-to-noise spectrum with a forest of emission and absorption lines all over the soft X-ray band, a spectrum quality comparable to the XMM-Newton/RGS spectra of very bright Galactic novae (see, e.g., Ness et al. 2011), with the remarkable difference that the ULX is located 4 Mpc away, while novae are within a few kpc from us.

The X-IFU spectrum will detect dozens of emission and absorption lines well above $5\,\sigma$ each and provide strong constraints on emission lines from triplets, such as Ne IX, and other lines sensitive to plasma density. More importantly, we will be able to detect winds in ULXs at much larger distance, therefore boosting the explored volume (dozens of ULXs with winds detections) and the wind parameter space including its energetics and variability. This will be crucial to fully address the relationships between the spectral hardness of the source and the wind properties already pointed in Eq. (5) and Fig. 7. It is also worth mentioning that the X-IFU will enable, for the first time, short-term variability studies of ULX winds and their connections with X-ray lags.

5 Conclusions

In this work we have provided the first insights on the thermal stability of winds in ultraluminous X-ray sources (ULXs), currently the best candidates of astronomical persistent sources to accrete beyond the Eddington limit. We focus on the case where the compact object is either a black hole or a non-magnetic neutron star for which a thick disc would still form and launch powerful winds through radiation pressure. Using three archetypal ULXs with wind detections, we have probed the thermal stability for three strategic spectral states that broadly describe the ULX phenomenology: soft, hard and broadened-disc state. We have found that winds in ULXs are likely in stable thermal equilibrium with the harder states progressively extending the unstable branch. The stability curves resemble those of Narrow Line Seyfert 1 and supersoft sources as expected at high-Eddington rates. Our results predict the existence of the evacuated funnel in ULXs as the wind is unstable and/or highly-ionized in innermost regions, likely providing a better view of the hard X-ray region.

We have also investigated the consequence of X-ray spectral variability and, in general, dramatic variations in the ionizing field by studying some ad-hoc spectral states of ULX NGC 1313 X-1. The results have confirmed that over the periods of time where strong variability occurs (typically weeks-to-months in ULXs) the thermal stability also significantly varies, which can explain the disappearance of strong residuals (and therefore winds) during bright states. Precession and accretion rate changes could both contribute to this. We also find a possible correlation between the spectral hardness of the ULX, the wind velocity and the ionization parameter in support of the overall scenario and the geometry of super-Eddington accretion discs.

Appendix A Technical detail

Here we show some technical information that we avoid in the main body of the paper to facilitate the reading.

A.1 Observation log

In Table 3, we report the XMM-Newton archival observations used to build the SEDs. In particular, the three longest observations taken during the intermediate-hard state where NGC 1313 X-1 spends most time and two specific observations yielding the softest and brightest state recorded for this source. The latter two have been used in Sect. 3.4 to compute the ionization balance of NGC 1313 for different spectral shapes (see Fig. 5). It is remarkable how NGC 1313 X-1 reaches extremely faint, soft states (e.g. for observation 0205230401) where its hardness ratio agrees (within error bars) with the brightest state of the supersoft ULX NGC 247 (see Table 1). This is of course expected from the unification scenario owing to a combination of lower accretion rates and system precession (see Middleton et al. 2015 for more detail on spectral variability of ULXs).

The other ULX observations and RGS spectra used to produce the results shown in Table 1 and 4 and in Fig. 7 are the following: 0200980101 (Holmberg IX), 0200470101 (Holmberg II), 0691570101 (NGC 6946 X-1), 0142770101-301, 0405690101-201-501, 0693850701-1401 and 0741960101 (NGC 5204, e.g. Kosec et al. 2018a), 0655050101 (NGC 55) and 0728190101 (NGC 247, e.g. Pinto et al. 2017).

A.2 Heating and cooling rates

In Fig. 9 we show the heating and cooling rates computed for the broadband SED of NGC 1313 X-1 using the intermediate-hard time-average spectrum (see Sect. 3.1, Fig. 3 and Table 3). These rates are computed using the SPEX pion code. These plots help to understand the behaviour of different physical processes at different ionization parameters, i.e. in different wind conditions. For our calculations, we set the hydrogen density to $n_{\rm H}=10^{8}$ cm*-3* and the total column density to $N_{\rm H}=10^{22}$ cm*-2* similarly to Mehdipour et al. (2016). We notice that in the regime of non-optically-thick gas the main effect of the density is a systematic scaling of all heating and cooling rates.

At first, we notice that the total heating and cooling rates are broadly consistent, which means that a converging, equilibrium solution is found at nearly any value of the ionization parameter, $\xi$. At low $\xi$ the dominant processes are heating by photo-electrons and cooling by collisional excitation. At high $\xi$ the largest contribution to heating is given by Compton scattering and to cooling by inverse Compton scattering. There is therefore a radical change in the dominant physical process with the ionization parameter and, therefore, the wind temperature. This is consistent with what Mehdipour et al. (2016) found comparing blackbody-like SEDs to harder spectra of AGN where other processes, such as auger electrons and recombination, provide more significant contribution to heating and cooling, respectively. They also compared the results obtained by SPEX pion with those using different spectral codes such as XSTAR and CLOUDY, finding some interesting differences of 10–30% in the corresponding estimates of $\xi$ and the concentration of several ionic species. However, they argue that these differences are unlikely to impact the scientific interpretation of current (AGN) outflow observations.

A.3 Wind properties in ULXs

In Table 4 we quote the hardness ratios of each source calculated as $HR=L_{2-10\rm keV}/L_{0.3-10\rm keV}$ using a $bb+mbb+pow$ spectral continuum (soft blackbody + broad modified blackbody + powerlaw, see e.g. Kosec et al. 2018b) for the XMM/EPIC-pn spectra from observations with evidence of winds (see Fig. 1 and Sect.1.2). We also report the velocities and the ionization parameters of the wind along with the references used. These data are used to produce Fig. 7.

For NGC 5204 X-1 we fit the data with the new $pion$ model in emission adopting $\Omega=2\pi$ in order to obtain a measurement of the ionization parameter (see, e.g., Psaradaki et al. 2018). This is done because Kosec et al. (2018a) show spectral fits using only line-emission from gas in collisional equilibrium, which reproduces the lines better than photoionization equilibrium (but deeper observations are necessary to distinguish between a collisionally-ionized and a photo-ionized nature of the outflow).

For NGC 1313 X-1 we use the fastest and the slowest solutions found by Pinto et al. (2016) and Walton et al. (2016) because they yield the best description of the soft and hard X-ray components of the wind and also provide an idea of the systematic uncertainties.

We double-check the values for Holmberg II X-1, Holmberg IX X-1 and NGC 6946 X-1 by fitting the RGS+EPIC/pn spectra used in Kosec et al. (2018a) and Pinto et al. (2016). Here, we report the results obtained for the solution that correspond to the highest significance. We still quote the wind detection in these sources as tentative because the Monte Carlo simulations yield a significance below $3\,\sigma$, likely owing to the limited statistics. The previous results have all been double-checked by testing the new version of SPEX used in this paper (3.05.00).

In Table 5 we quote the Spearman and Pearson correlation coefficients for the relationships between the ULX spectral hardness, $HR$, the wind ionization parameter, $\xi$, and velocity, $v_{w}$. They are calculated using the scipy.stats.spearmanr and scipy.stats.pearsonr PYTHON routines for the dataset shown in Table 4. Despite the large scatter, the results indicate a substantial positive correlation among the three parameters, which agrees with the predictions for super-Eddington disc geometry and winds.

Acknowledgments

This work is based on observations obtained with XMM-Newton, an ESA science mission funded by ESA Member States and USA (NASA). CP is supported by European Space Agency (ESA) Research Fellowships. MM is supported by the Netherlands Organisation for Scientific Research (NWO) through the Innovational Research Incentives Scheme Vidi grant 639.042.525. ACF and PK acknowledge support from ERC Advanced Grant Feedback 340442, the European Union’s Horizon 2020 Programme under the AHEAD AO5 project (grant agreement n. 654215) and the Faculty of the European Space Astronomy Centre (ESAC) under the project number 459. DJW and MJM thank STFC for support in the form of an Ernest Rutherford Fellowship. TPR acknowledges funding from the STFC via consolidated grant ST/P000541/1. We thank Jelle de Plaa for support in optimising SPEX, Douglas Buisson and William Alston for support with the XMM-Newton/OM. We also thank Manfred Pakull for useful discussions on the effects of ULX winds onto the surrounding bubbles and Ioanna Psaradaki for useful support with the PYTHON and the SPEX pion codes. We finally acknowledge the anonymous referee for useful comments that significantly improved the paper.