© J.-M. Hameury and J.-P. Lasota 2020

1. Introduction

According to the disc instability model (DIM; see Lasota 2001; Hameury 2020, for reviews of the model), accretion discs around compact objects are subject to a thermal–viscous instability if the rate at which matter is brought to their outer edge is less than a critical value strongly increasing with radius. This basic tenet of the DIM has been convincingly confirmed by observations of dwarf novae (Dubus et al. 2018) and transient X-ray sources (Coriat et al. 2012).

Since the critical accretion rate increases with radius, for sufficiently large orbital periods even potentially very bright systems should have large unstable discs and exhibit some types of outbursts. This could be the case of ultraluminous X-ray sources (ULXs; see Kaaret et al. 2017, for a review), the majority of which, as has been well established, are apparently super-Eddington luminosity X-ray binaries containing either stellar-mass black holes or neutron-star accretors (as predicted a long time ago by King et al. 2001). Although it is very difficult to determine the binary parameters of these distant systems (distances up to 20 Mpc) it is clear that a large fraction of them contain giant or supergiant stars¹ and have orbital periods larger than two days, and thus could have large unstable discs. The recent observation of a very well-sampled transient ULX source in the galaxy M 51 by Brightman et al. (2020) has provided a confirmation of this hypothesis. Since some of the usual, sub-Eddington X-ray transients show long-lasting outbursts (tens of years), this could also be the case for some apparently steady ULXs.

Despite these potential applications, until recently modelling the accretion-disc instability has focused on systems with relatively short orbital periods, and disc sizes not exceeding 10¹¹ cm, with the notable exception of Bollimpalli et al. (2018) who applied the DIM to two symbiotic stars, Z And and RS Oph, which have orbital periods of 759 d and 454 d respectively, implying outer disc radii larger than 10¹² cm. Dubus et al. (2001) also considered discs in LMXBs with radii up to 10¹² cm, but considered relatively low accretion rates (sub-Eddington).

The reason for the paucity of models with long orbital systems is that modelling large discs is numerically challenging, mainly because of the large ratio between the inner and outer disc radii (up to six orders of magnitude). Lasota et al. (2015) attempted to describe outbursts of large discs using analytical formulae based on the properties of models for standard X-ray transients (XRTs) but the accuracy of such a method and its application to ULXs should be subject to caution and requires testing through numerical calculations.

The main aim of the present work is to extend and generalize the irradiated-DIM to the case of large (outer radius larger than 10¹² cm) accretion discs and apply it to the description of systems with super-Eddington accretion rates at outburst and, in some cases, super-Eddington mass transfer rates. Such models of ultraluminous X-ray transients can be compared to observations; this helps in understanding both the nature of ULXs and the mechanisms driving accretion in astrophysical discs.

We do not claim that all (or even most) ULXs are transient systems in which the mass transfer rate from the secondary is sub-Eddington or moderately super-Eddington. Some ULXs, however, are observed to be transient, such as HLX-1 in ESO 243-49 (Farrell et al. 2009) or M51 XT-1 (Brightman et al. 2020), and one needs therefore to assess whether the DIM can account for the observed outburst properties, and, if the answer is yes, under what conditions. On the other hand, many ULXs appear to be permanently in a bright state, but because of the long timescales involved in large discs, some of them could be caught during outbursts lasting for years or even decades, while other sources, presumably the vast majority if not all, are genuinely permanent² sources, in which case the DIM sets constraints on the disc size and the mass transfer rate.

As recalled in Sect. 3, the thermal–viscous stability of discs depends mainly on the mass transfer rate from the secondary, the disc size, the viscosity parameter, and disc irradiation. Our results are therefore, to a very good approximation, independent of the nature of the secondary³. It can be a low-mass or a massive star, it may or may not fill its Roche-lobe; if the mass transfer rate is lower than a critical value depending mainly on the size of the outer disc radius, the system will be transient.

In Sect. 2, we first briefly summarize the ingredients that enter the DIM as applied to XRTs; we give our results for a grid of models in Sect. 3, and we show that large accretion rates in outbursts exceeding the Eddington value can be attained when the mass transfer rate is large enough and the heating front propagates to large distances in the accretion disc, possibly reaching its outer edge. This, however, requires that the irradiation efficiency increases and compensates for the decrease in the accretion efficiency when the Eddington luminosity is exceeded. In Sect. 4 we then provide analytical approximations for important properties of the outbursts (including peak accretion rate, outburst decay time, and duration) and compare our results with selected observations in Sect. 5.

2. Model

We follow the thermal–viscous evolution of an accretion disc using our code described in Hameury et al. (1998) and Dubus et al. (2001). This code solves the equations for mass and angular momentum conservation, in which the viscosity is parametrized according to the Shakura & Sunyaev (1973) prescription, as well as a thermal-balance equation that includes viscous dissipation, heating by tidal torques, and energy advection. Because the disc is geometrically thin, the vertical structure of the disc can be decoupled from its radial variations; the heating and surface cooling terms Q⁺ and Q⁻ that enter the thermal equations are known functions of the radius r, surface density Σ, and mid-plane temperature T_c. We have extended our grid of pre-calculated Q⁺ and Q⁻ to radius values of 10¹³ cm; at these large distances, the disc temperatures can be very low (below 1000 K) in quiescence, implying that the opacities may have to be extrapolated, and are therefore uncertain. For low temperatures, we use the tables from Alexander (1975) that extend down to 700 K; more recent and more accurate tables do exist (see e.g. Alexander & Ferguson 1994), but it is unclear whether the assumptions used for estimating the abundance of grains and molecules are appropriate for accretion discs. However, we note that, as shown in Sect. 3, the structure of the distant regions is used only to describe the heating and cooling front propagation when the temperature is larger than typically 10⁴ K, so that the inaccuracy of the opacities at low temperatures is not really a problem.

The viscosity parameter α is taken to be bimodal, with a value α_c on the cold branch, α_h on the hot branch, and a smooth transition at temperature T_crit that is the average of the turning points and of the S-curve that describes the equilibrium in the Σ – T_c plane. We use the analytical fits given by Dubus et al. (2001), and we checked that these fits are also appropriate for much larger radii than considered there.

2.1. Accretion luminosity

When the accretion rate Ṁ reaches

(1)

the accretion luminosity corresponds to the Eddington value for an accretion efficiency of 0.1. The Eddington accretion rate is therefore Ṁ_Edd. For a 7 M_⊙ black hole, it is equal to 9.1 × 10¹⁸ g s⁻¹.

When the accretion rate is super-Eddington, Ṁ(t) cannot be considered to represent the X-ray light curve any longer. In calculating light curves of super-Eddington outbursts, we assume that the apparent luminosity is given by (King 2009)

(2)

In this relation, ṁ = Ṁ/Ṁ_Edd, and is a beaming term; the larger the beaming parameter , the larger b and hence the smaller the beaming effect. King (2009) found that the beaming factor has the form and determined . Here, we substituted the original beaming term , which is valid only for by in order to get a smooth transition with the case where beaming is negligible.

King (2009) determined using both theoretical arguments and observational constraints; this value is therefore uncertain, but, as shown below, the peak ṁ we obtain in our models is often moderate, typically less than about ten, so that the influence of the beaming term is also moderate. It is also worth noting that, although the comparison with observations requires a relation between L_x and Ṁ, our models are independent of beaming, but depend on the relation between the bolometric luminosity and Ṁ in an indirect and, as we see in the next section, hidden way, via the disc irradiation term.

The actual accretion efficiency η_t, defined as the ratio of the true luminosity (without the beaming factor) to c², is therefore

(3)

Equation (2) also shows that, unless the value of the Eddington luminosity is increased by the presence of a magnetar-type magnetic field, beaming is necessary for the luminosity to significantly exceed the Eddington critical value as is the case in the pulsating ultraluminous X-ray sources (PULXs), which contain neutron stars that are very unlikely to be magnetars (King & Lasota 2019) and have apparent luminosities up to a few hundred times larger than the Eddington value (see Table 2 in King & Lasota 2020).

2.2. Outer and inner disc radii

The outer disc radius is allowed to vary, and is limited by the tidal torques exerted on the outer parts of the disc. Here, as in Viallet & Hameury (2008), we assume that these torques have a steep exponential dependence on radius, so that they are negligibly weak for a radius smaller than the tidal truncation radius r_tid and very large beyond that radius, effectively setting r_tid as the maximum disc radius.

We assume that the inner disc is truncated during quiescence because of the formation of a hot, optically thin, and radiatively inefficient flow (see e.g. Lasota et al. 1996; Narayan et al. 1997; Narayan & McClintock 2008). For simplicity, and because the precise dependence on the transition radius between the classical disc and the hot flow is of little importance as far as the time evolution of the disc is concerned – what matters is the actual inner disc radius in quiescence just prior to the onset of an outburst – we assumed that the inner disc is that of an accretion disc truncated by a fictitious magnetic field with effective magnetic moment μ. We considered values of μ of the order of 10³¹–10³² G cm³ so that the inner disc radius is a few 10⁹ cm during quiescence. The “magnetic field” we use is just a proxy for the unknown but necessary (at least at some phases of the disc’s evolution, see e.g., Basak & Zdziarski 2016; Zdziarski & De Marco 2020) “evaporation” mechanism.

In outburst, the inner disc radius should be close to the compact accretor (assuming the neutron star is not strongly magnetized), which is ∼10⁷ cm. In the case of very large discs (outer radius ≳10¹² cm) the resulting large radius ratio is numerically too challenging to be implemented. Therefore, as in Dubus et al. (2001) we limit the inner disc radius to 10⁹ cm. Since near the outburst peak the inner disc is hot, the viscous timescale at which density perturbations evolve is short enough there for the adjustments of its structure to be instantaneous compared to the outburst evolution timescale, as long as the inner radius is low enough. Dubus et al. (2001) have checked that light curves produced with inner radii of 10⁷ and 10⁹ cm are practically identical and we expect this to be also true in our case. Truncation at 10⁹ cm also has the advantage of allowing us to avoid problems with the infamous radiation-pressure instability that appears persistently in models (Jiang et al. 2013, see, however, Sądowski 2016; Lančová et al. 2019) but seems to be absent in the real Universe since neither steady sources, nor black hole transients decaying from outburst show any characteristic variability when observed in the range of luminosities where this instability is supposed to operate⁴ (see also Gierliński & Done 2004). In our 10⁹ cm–truncated models the contribution of radiation pressure never exceeds 25%.

If, during outbursts, the accretion rate exceeds Ṁ_Edd, it is very likely that large outflows are generated at small radii, typically below the spherization radius of r_sph = (27/4) ṁ × GM₁/c² (Shakura & Sunyaev 1973). For all the models we calculated, r_sph was found to be significantly smaller than 10⁹ cm; in any case, it would be impossible to calculate the vertical disc structures for r <  r_sph. However, this effect is taken into account when calculating the luminosity variations using Eq. (2).

In three of the models we calculated, the disc temperature in the innermost disc regions exceeded the limits of the pre-calculated grid providing Q⁺ and Q⁻; we found it convenient to truncate the disc at a radius of 3 × 10⁹ cm or 5 × 10⁹ cm for these two particular models when the accretion rate exceeded 2  ×  10¹⁸ erg s⁻¹, and we checked on other models that this approximation did not introduce any change in the properties of the light curves.

2.3. Disc self-irradiation

Observations show that the outer parts of accretion discs in low-mass X-ray binaries (LMXBs) are irradiated by the central X-ray source (van Paradijs & McClintock 1994). This irradiation determines the disc stability criteria (van Paradijs 1996; Coriat et al. 2012) and strongly influences the properties of outburst light curves (Dubus et al. 2001; King & Ritter 1998). Unfortunately, in spite of several observational and theoretical attempts to discover what the irradiation mechanism and geometry are (it cannot be direct irradiation, Dubus et al. 2001), we still do not know how and by what accretion discs in X-ray binaries are illuminated.

Faced with this situation, we had no choice but to reach for the prescription by Dubus et al. (1999) that we used twenty years ago (Dubus et al. 2001). This ansatz has the advantage of being simple, physically motivated, and, what is most important, providing a correct stability criterion (Coriat et al. 2012).

The irradiation formula in question is obtained considering the vertical structure of an accretion disc modified by adding an extra term to the standard boundary condition , where σ is the Stefan–Boltzmann constant, F is the vertical heat flux in the disc generated by viscous dissipation, T_s is the disc surface temperature, T_irr is the irradiation temperature given by

(4)

𝒞 contains all the physics of the irradiation process, and η_t is defined by Eq. (3) (except for low accretion rates, see below). Dubus et al. (2001) found that the light curves of low-mass X-ray transients are reasonably well reproduced, and Coriat et al. (2012) found that the corresponding stability criterion provides the observed division of sources into steady and outbursting if one uses a constant η_t𝒞 of the order of 10⁻³. We stress that our definition of 𝒞 is the same as in Dubus et al. (2001), but differs from the one used in Dubus et al. (1999) in that it does not include the efficiency η_t. In the following, we use

(5)

to quantify the effect of irradiation. As in Dubus et al. (2001), η_t is reduced at low accretion rates, because for those low rates, the accretion flow below the truncated disc becomes radiatively inefficient. We use here a cut-off term [1 + (Ṁ/10¹⁶ g s⁻¹)⁻⁴]⁻¹.

For large mass accretion rates, we consider two extreme cases: (a) C is constant below the Eddington luminosity, and is reduced by a factor (1 + ln ṁ) above it, and (b) η_t𝒞 is constant at all accretion rates. Case (a) corresponds to a constant irradiating flux above the Eddington value, while case (b) implies that the decrease in η_t is compensated by an equivalent increase in 𝒞. We also consider an intermediate case where C is constant (case a1).

The increase in 𝒞 could be due to strong thermal outflows driven by X-ray irradiation of the outer accretion disc that would scatter efficiently the X-rays emitted in the vicinity of the accreting compact object, as proposed by Dubus et al. (2019). They estimated in a self-consistent manner the efficiency of disc irradiation from a wind model. However, their approach can apply to sub-Eddington rates only and we consider systems with a larger orbital separation than theirs. In addition, when the model is applied to an observed system, the conclusion is that scattering in the wind is still not sufficient to produce the observed heating, even in combination with direct illumination (Tetarenko et al. 2020). Here, since discs are observed to be irradiated, we adopt a simple (even simplistic) “pragmatic” approach in order to produce models that can be tested by observations, as has been the case for X-ray transient models tested by Tetarenko et al. (2018a,b).

3. Results

In the following, the primary is either a black hole with mass M₁ = 7 M_⊙ or a neutron star with M₁ = 1.4 M_⊙; the secondary mass is in all but one cases M₂ = 0.4 M_⊙. We consider three possible values for the orbital period P_orb: 155 h, 400 h, and 1200 h, corresponding to an outer disc radius of approximately 10¹², 2 × 10¹², and 4 × 10¹² cm respectively in the black hole case, and to 5 × 10¹¹, 10¹², and 2 × 10¹² cm when the accretor is a neutron star. The value of M₂ = 0.4 M_⊙ is obviously not suitable for mass transfer rates up to 3 × 10⁻⁷ M_⊙ yr⁻¹ but we use this mass to make the models homogeneous. In any case, we are interested in considering the impact of the disc radius on the outcome of the model, and not on modelling a specific system, with the exception of V404 Cyg and M51 XT-1. When one modifies both the secondary mass and the orbital period in such a way that the outer disc radius remains the same, the only disc parameter that is different is the specific angular momentum of matter incorporated at the outer edge of the disc, which is usually parametrized by the so-called circularization radius r_circ. As we show below (models 34 and 34a in Sect. 3.1), the outburst properties are essentially independent of r_circ. The secondary mass we have chosen corresponds to the properties of V404 Cyg, the binary parameters of M51 X-1 are not known. We finally note that if the disc is fed via a wind from the secondary instead of via Roche lobe overflow, the outer disc radius is smaller than what we infer from the orbital separation; again, because of the very weak dependence on r_circ, our results are also valid in this case provided one considers the actual value of the disc radius.

In the following, the mass transfer rate is taken as a free parameter; in real systems, it will be determined by the orbital separation and the secondary properties in a complex way that we do not address here (see e.g. Wiktorowicz et al. 2017, for a discussion of the various formation channels of ULXs). But even in the case where the secondary is more massive than the primary and experiences mass loss on a thermal timescale, the disc can in principle still be unstable if it is large enough.

3.1. Outburst sequences

We have calculated the outburst properties for various sets of parameters. Tables 1 and 2 summarize our results in the case where the compact object is a 7 M_⊙ black hole (Table 1), or a 1.4 M_⊙ neutron star (Table 2). For each model, defined by the mass transfer rate from the secondary Ṁ₂, illumination factor f_irr, magnetic moment μ₃₀ in units of 10³⁰ G cm³, α_c (we kept α_h constant and equal to 0.2), we provide the duration of a cycle of outbursts t_c, defined as the periodicity of the light curve that may cover several outbursts of different intensities, which we classified as large (L; these are the largest outbursts of the sequence), intermediate (m, with peak luminosities ranging between 10 and 100% of the main outburst values), small (s, 1 to 10% of the main outburst peak luminosity), and very small (μ, less than 1% of the main outburst peak luminosity) and for each model, we provide the pattern of the outburst sequence. Tables 1 and 2 also provide the average recurrence time t_r between outbursts, which is equal to t_c divided by the number of outbursts in the pattern, and, for the largest outbursts, the outburst duration d in years, the maximum accretion rate Ṁ_max, the maximum distance reached by the heating front r_tr, 12 in units of 10¹² cm, the total mass ΔM accreted during the outburst, and ΔM/M_d, where M_d is the disc mass.

Models 34 and 34a differ only by the secondary mass, taken arbitrarily to be 5 M_⊙ in model 34a, and the orbital period chosen in such a way that the outer disc radius is the same in both cases. Although r_circ differs by 40% between both cases, the outburst properties are, at the percent level, identical. Whether systems with such secondary masses may exist in nature is irrelevant; our point here is merely to show that the outburst characteristics do not depend on M₂ for a given outer disc radius.

Figures 1–3 show examples of light curve patterns for model 1, a classical case of sub-Eddington outbursts, and model 12, a slightly super-Eddington case with ṁ = 5.5, at the outburst’s maximum; its light curve exhibits a complex outburst sequence. Figure 3 shows that for moderately super-Eddington accretion rates (ṁ < 8.5) the difference between the L_x and Ṁ “light curves” is negligible for all practical purposes. We note that the jump in the transition radius observed when r_tr reaches the outer disc radius is due to the fact that, because of tidal heating, the outer disc edge is already in the hot state shortly before the heating front reaches it.

Fig. 1. Time evolution for model 1, which corresponds to parameters appropriate for V404 Cyg. Upper panel: disc mass. Intermediate panel: mass accretion rate. Bottom panel: outer disc radius (red curve), inner disc radius (blue curve), and transition radius rtr (black curve). For such sub-Eddington outbursts, mass-accretion variations represent the X-ray light-curve.

Fig. 2. Same as Fig. 1 for model 12.

Fig. 3. Detailed view of a large outburst in Fig. 2. Top panel: bolometric light curve (red) and accretion rate times 0.1c2. Bottom panel: outer disc radius (red curve), inner disc radius (blue curve), and transition radius rtr (black curve).

Figure 4 explains why a complex sequence of outbursts can be found in large discs. After a major outburst, the disc density is low; after a time of the order of the diffusion time in the inner parts of the disc, Σ reaches the critical density at which the disc can no longer remain in a cold quasi-steady state, and an inside-out outburst is triggered. Because Σ in the outer disc is still low, the heating wave cannot propagate very far in the disc, and the outburst amplitude is small. The outburst amplitude, however, has been sufficient to significantly modify Σ in intermediate regions of the disc, so that, during the next outburst, the heating front is able to propagate to larger distances, which empties these intermediate regions; during the following outburst, the heating front stops at a smaller distance than for outburst number 2. The sequence continues, and the mass builds up in the outer disc as a result of mass transfer. Eventually, the heating front reaches the outer disc edge, triggering a major outburst. Such a complex sequence of outbursts is reminiscent of the sequences found when modelling discs around active galactic nuclei (AGNs; Hameury et al. 2009), in which heating fronts are never able to reach the outer edge of the disc.

Fig. 4. Radial surface-density profiles of the accretion disc just prior to outbursts for model 12 (solid lines). The curve labelled 1 shows Σ prior to the first weak outburst after a major one, curve 2 shows Σ before the second outburst, curve 3 is just before the last weak outburst of the sequence, and curve 4 represents Σ at the onset of the last outburst. Dashed lines show the critical surface densities Σmin and Σmax of the S-curves.

3.2. Impact of disc irradiation

In all models, except model 24, we have used a constant value of f_irr. Figure 5, corresponding to model 24, illustrates what happens when in the black-hole case, irradiation is Eddington-limited by setting f_irr to be

(6)

Fig. 5. Time evolution for model 24, in which the irradiation efficiency firr decreases when the accretion rate reaches the Eddington value.

As can be seen, after an initial rapid decay, the accretion rate remains approximately constant at a level slightly above the Eddington value. During this plateau, irradiation remains constant, which prevents the cooling front from propagating inwards; the whole disc, even in the cool outer regions, is close to steady state, and evolves slowly on a characteristic timescale of M_d/Ṁ until Ṁ falls below the Eddington value.

The long plateau observed for model 24 is somewhat reminiscent of the one observed in XTE J1550-564 (Sobczak et al. 2000), and is due to the constant irradiation temperature at a given disc radius when the accretion rate exceeds the Eddington value, which keeps the position of the cooling front constant. Model 24a corresponds to the same parameters as model 24, with the difference that 𝒞 is now assumed to be constant, so that f_irr varies as 1 + ln ṁ for ṁ > 1. Outbursts now have the usual shape; they are weaker and shorter than in the case where f_irr is constant for all luminosities.

For completeness, we have also calculated two models in which irradiation is not taken into account (models 4 and 5). As mentioned above, the recurrence time is short and the total outburst fluence is also much lower than in the irradiated case, even though the peak outburst luminosity is comparable to the irradiated case, leading to a short outburst duration. One should also note that outbursts of non-irradiated discs are terminated by a sequence of reflares, as shown in Dubus et al. (2001) for smaller accretion discs.

We note that in all our neutron-star models, the accretion rate at maximum is super-Eddington (Ṁ_Edd = 1.8 × 10¹⁸ g s⁻¹) and it is also the case of fourteen (out of twenty-four) black hole models. The mean mass transfer rate is only slightly super-Eddington in models 23, 27, and 35; it is sub-Eddington for all other models, so that none of them has a mass transfer rate usually envisaged for “steady” ULXs.

The main differences between black hole and neutron star systems are that, for a given orbital period, the disc around a neutron-star binary is smaller, and the Eddington luminosity is lower, which implies that the disc might be more prone to mass loss via a wind for example. From Tables 1 and 2, one can also note that the characteristics of an outbursts do not depend on the disc size, as long as the heating front does not reach the outer disc edge. This is not surprising as the outer disc regions are not affected by the propagation of heating and cooling fronts and play no role in the outburst outcome.

3.3. Recurrence time

Since the condition for the occurrence of outside-in outbursts is (Lasota 2001)

(7)

where δ ≤ 2 and r_out is the outer disc radius in quiescence, all models considered here correspond to inside-out propagation of the heating front, meaning they correspond to outbursts starting near the edge of the inner disc.

For such outbursts, the recurrence time is approximately equal to the viscous time at the innermost parts of the accretion disc, and should thus depend only on α_c and on the inner disc radius in quiescence. Tables 1 and 2 show that the recurrence time is the same for all models that share the same α_c, with the exception of models 4, 5, and 24, and is almost independent of the mass transfer rate. In models 4 and 5, the recurrence time is short because neglecting irradiation causes reflares that alter the surface density profile as compared to the case where a cooling front is able to propagate directly to the disc inner edge (Menou et al. 2000; Dubus et al. 2001). In model 24, the inner disc is largely depleted after long outbursts and slow decay to quiescence (see below); the inner disc radius is larger than in other models, leading to a long recurrence time.

The recurrence time of outbursts in X-ray irradiated discs is longer than in dwarf novae (except for WZ Sge-type “super-outbursts”), and is not given by Eq. (53) in Lasota (2001), which applies to systems in which outer-disc irradiation is negligible (such as dwarf novae). In such a case the quiescent Σ is in the range [Σ_min, Σ_max], where Σ_min and Σ_max are the turning points of the S curve in the Σ – temperature diagram. As explained in Lasota (2001), this is no longer true when, due to irradiation, the disc, and in particular its inner portions, is strongly depleted during outbursts. This is illustrated by Fig. 6, which shows the Σ(r) profile at the end of an outburst; Σ is much below the non-irradiated (irradiation plays no role in quiescence) Σ_min, which explains the long recurrence time in these models: at the start of quiescence, quiescent discs can be almost empty. As seen in Tables 1 and 2, up to 90% of the disc mass can be accreted during outbursts.

Fig. 6. Quiescent radial surface-density profile 173 days after the end of an outburst for model 26. The dashed lines correspond to critical surface densities Σmax (upper line) and Σmin (lower line).

The recurrence time can also be estimated from the disc refilling time,

(8)

in agreement with the t_r values given in Tables 1 and 2.

In contrast with the dwarf nova case, the disc mass in our transients is never close to maximum (corresponding to Σ ∼ Σ_max) as clearly seen in Fig. 4. This said, one should keep in mind that our standard assumption that in quiescence the disc accretion is driven by a process describable by an α–prescription is far from being guaranteed. As studied in detail by Scepi et al. (2018), in accretion discs with T ≲ 3500 K and Σ ≲ 180 g cm⁻², the magnetic Reynolds number R_m = c_sH/η_R, where c_s is the sound speed, H the pressure scale-height of the disc, and η_R the Ohmic resistivity coefficient, is less than 10⁴ so that diffusion of the magnetic field becomes too important for the disc to sustain the magnetohydrodynamic turbulence (Hawley et al. 1996; Fleming et al. 2000, see also Gammie & Menou 1998). Observed quiescent X-ray fluxes prove that accretion is also occurring during this phase of the outburst cycle (see e.g. Lasota 2000) but since in quiescent X-ray transient discs T is less than 3500 K and Σ less than 180 g cm⁻² (see for example Fig. 6), it is not known what drives accretion in quiescent XBTs. The calculated values of recurrence times should be therefore subject to caution, even if they are in good agreement with observations (compare for example with Coriat et al. 2012).

4. Analytical relations

4.1. Peak accretion rate

As noted by Lasota et al. (2015), a tight correlation exists between the maximum accretion rate during an outburst Ṁ_max and the maximum distance reached by the transition front r_tr, max, because at the outburst peak, the portion of the disc that is in the hot state is close to being steady, and the mass accretion rate is thus equal to the minimum (critical) rate for which a hot stable disc can still exist. This critical rate is well approximated by

(9)

Here we use the formulae for critical quantities from Lasota (2001) that are slightly different from those in Lasota et al. (2015) who used fits obtained in Lasota et al. (2008). The version of the DIM code used in this paper differs in some minor aspects from that used in this reference.

Figure 7 shows all our models in the plane (r_tr,max, Ṁ_max) where Ṁ_max has been normalized by . We have omitted models 4 and 5, in which irradiation is not taken into account, as well as model 24 in which f_irr varies.

Fig. 7. Relation between the maximum mass accretion rate reached during an outburst multiplied by and the maximum distance travelled by the heating front. Stars correspond to black hole systems, circles to neutron star systems. Red symbols correspond to an orbital period of 155 h, green symbols to 400 h, and blue symbols to 1200 h. Symbols are filled when the heating front reached the outer edge of the disc. The dashed line represents the critical mass transfer rate as given by Eq. (9).

As can be seen, the agreement between the maximum accretion rate during an outburst and is remarkable as long as the heating front does not reach the outer-disc edge. Moreover, if we use the fits for appropriate for the unirradiated case given in Bollimpalli et al. (2018), we find that the peak accretion rate for model 4 is 1.36 times the critical rate at the transition radius, and this ratio is 1.07 for model 5, so that this agreement also holds for the two models that are not shown in Fig. 7.

If the heating front reaches the outer edge of the disc, the peak accretion rate is larger than what is predicted by Eq. (9). This corresponds to filled symbols in Fig. 7. The Ṁ_max can reach values higher than predicted by up to a factor 50. The higher Ṁ₂, the higher Ṁ_max, with a limit due to the fact that Ṁ₂ cannot be too high in transient systems, as we discuss in the next section. The upper red dot in Fig. 7 corresponds to model 27, for which Ṁ₂ = 2 × 10¹⁸ g s⁻¹ is rather close to the stability limit (4.9 × 10¹⁸ g s⁻¹). However, this does not prevent the disc from exhibiting fully-fledged outbursts.

It is worth noting that the r_tr, max–Ṁ_max relation also approximately holds to within a factor of two for the actual accretion rate during an outburst, with Ṁ slightly below the fit during the rise and slightly above during decay. In Sect. 4.2.2 we use this property to derive the propagation time of cooling fronts.

The heating front can reach the outer disc edge if the mass transfer rate is high enough. A close look at Table 1 shows that, to first order, Ṁ_max/Ṁ₂ depends mainly on the ratio α_h/α_c, provided that the heating front does not reach the outer disc edge. Ṁ_max/Ṁ₂ is of the order of 25–35 when α_h/α_c = 5, and of the order of 70–100 when α_h/α_c = 10. Deviations are observed when there is a long tail in the outburst, as in the case of Eddington limited irradiation; the low Ṁ_max/Ṁ₂ found in this case results from the depletion of the accretion disc during the plateau that prevents the heating front from propagating at large distances during the next outburst. It also appears that f_irr has a limited impact on Ṁ_max/Ṁ₂; a decrease in irradiation changes the relation between r_tr and Ṁ_max, but not Ṁ_max itself. The relation

(10)

also holds when the transition front reaches the outer disc edge, but with some deviation that increases with increasing Ṁ₂, in particular when the outburst duration is no longer much shorter than the recurrence time.

Using this approximate relation together with enables us to determine r_tr, max. If the latter quantity is larger than the outer disc radius r_out, the heating front reaches the outer disc edge, and the full disc is brought into the hot state. In the opposite case, a cooling front is initiated at r_tr, max.

4.2. Outburst decay time and duration

During outburst maximum and decay, the hot disc region is close to being steady (Lasota 2001). The surface density is then

(11)

where α_0.2 = α/0.2, Ṁ₁₉ = Ṁ/10¹⁹ g s⁻¹, r₁₂ = r/10¹² cm, and ψ ∼ 1 describes the deviation of the opacities from the Kramers’ law used by Lasota (2016) to obtain Eq. (11). Our calculations of the effective temperature as a function of surface density yield values for ψ varying between 1.0 and 1.5; in the following, we use ψ = 1.3. If the full disc is hot, its mass is

(12)

The estimates provided by Eqs. (11) and (12) require that irradiation does not significantly alter the disc vertical structure, which is the case when the irradiation temperature is significantly lower than the disc central temperature. We have checked that the disc mass obtained using Eq. (12) corresponds by better than a factor of 1.5 to the actual calculated disc mass when the heat front does reach the outer disc edge.

4.2.1. Decay of a fully hot disc

We first consider the case where the heat front reaches r_out. During the initial decay, the disc mass decreases while remaining entirely in the hot state; Eq. (12) relates the disc mass to its time derivative; it can be solved and we obtain (Ritter & King 2001):

(13)

where we corrected a misprint in Ritter & King (2001) Eq. (10). The constant of integration t₀ is given by

(14)

where Ṁ_max,19 = Ṁ_max/10¹⁹ g s⁻¹. Although related to it, t₀ is not the characteristic timescale of the disc evolution. We can define a characteristic timescale from Eqs. (12) and (9) as

(15)

where . This ratio is given by

(16)

and is less than ϕ since the mass transfer rate has to be less than the critical rate for the system to be transient.

We note that contrary to the common opinion (e.g. King & Ritter 1998; Dubus et al. 2001; Lasota et al. 2015) the light curves of outbursts produced by the instability of irradiated discs in X-ray transients systems are not exponential (see also King 1998).

Figure 8 compares the time evolution of the mass accretion rate as found in model 14 from our numerical simulations with the analytical estimate given by Eq. (13). The little wiggle seen at time t ∼ 4.2 yr is due to the ad hoc truncation of the disc during the outburst maximum for this model (see Sect. 2). The approximation is quite good and relies on two free parameters: the time at which the decay from maximum starts and the peak mass accretion rate. In the case shown in Fig. 8, we use Ṁ_max = 1.94 × 10²⁰ g s⁻¹, and the corresponding t₀ is 5.0 yr. As mentioned above, t₀ is longer than the characteristic evolution time, which is 2.0 yr. We also stress again that, although the decay is not very different from exponential, a −10/3 power law is by far a better fit.

Fig. 8. Time evolution of model 14. Bottom panel: transition radius (black curve) as well as the outer disc radius (red curve). Top panel: actual mass accretion rate (black curve) and critical rate (green curve). The blue curve gives the analytic solutions (Eq. (13)) when the entire disc is in the hot state, and the red curve corresponds to a propagating cooling front (Eq. (22)).

We can relate the peak accretion rate to the characteristic decay time by combining Eqs. (9) and (15); we find

(17)

If, from observations, both Ṁ_max and τ are known, this relation determines f, and hence and the size of the accretion disc.

A slightly different approach was used by Lasota et al. (2015) who assumed that the decay from maximum of a fully hot, irradiated disc begins with an approximately exponential phase, during which matter is accreted at a viscous time at the constant (during this phase) outer disc radius (Dubus et al. 2001; King & Ritter 1998), the characteristic viscous time being given by t ∼ r²/(3ν).

Therefore the decay time is equal to the viscous time and can be written as

(18)

where c_s = kT_c/m_p, with being the critical minimal midplane temperature of the hot irradiated disc. Taking T_c ≈ 27200 (r_out/10¹² cm)^0.05 K (Lasota 2001), we get

(19)

This is close to the decay time found in Eq. (15), although the parameter dependence is rather different. This difference strongly affects the Ṁ_max(τ) relation and explains why our Eq. (17) and Eq. (6) of Lasota et al. (2015) are so dissimilar.

4.2.2. Propagation of a cooling front

When Ṁ falls below Ṁ_crit(r_out), a cooling front starts propagating from the outer disc edge. The situation is similar if the heating front does not reach the outer disc edge but in this latter case, the cooling front starts inside the disc. We can then use the same method for determining the time evolution of the disc, the main differences being that the disc mass in Eq. (12) now refers to the mass of the fraction of the disc that is in the hot state, and Ṁ is no longer the time derivative of the hot disc mass. Instead we can write

(20)

where the dots indicate time derivatives and Ṁ_tr is the mass flow at the transition radius⁵. Our simulations indicate that Ṁ_tr is positive (i.e. mass flows outwards) and comparable to Ṁ (see Menou et al. 1999; Dubus et al. 2001, for a detailed discussion of the structure of the transition fronts). Since , we can express r_tr as a function of Ṁ from Eq. (9), and after some algebra we get

(21)

If we make the additional assumption that ξ does not vary much with time, we can solve Eq. (12) using Eqs. (9) and (21); we find:

(22)

where is a constant that is determined by the condition that, when the cooling front starts, Ṁ is equal to at the maximum transition radius. We can then write as

(23)

Figure 8 compares the time evolution of the accretion rate as determined by simulations for model 14 with the analytical approximation above. Here, there is only one free parameter, ξ, and we took ξ = 6.3. As can be seen, Eq. (22) represents the results of numerical simulations quite well. Interestingly, the differences are mainly due to the fact that the accretion rate is not exactly equal to the critical rate at the transition radius; Eq. (22) fits the critical rate very well, much better than the actual accretion rate Ṁ.

It must also be stressed that is equal to the duration of this phase, and it can thus be directly compared to observable quantities. Figure 9 shows the ratio of the duration d of the outburst divided by , for ξ = 5, corresponding to Ṁ_tr = Ṁ. Because the rise from quiescence is much shorter than the decay, d is a good estimate of the duration of the final decay phase. When the heating front does not reach the outer disc edge, d is close to , to better than a factor of 1.5. This is not the case when the full disc is brought to a hot state, because a phase exists during which the disc decays while being entirely in the hot state.

Fig. 9. Comparison between the outburst duration and the duration of the propagation time of the cooling wave as estimated by Eq. (22) with ξ = 5. The symbols have the same meaning as in Fig. 7.

It is interesting to note that and τ, as given by Eqs. (23) and (15), have the same functional dependence and are of the same order. This is not a surprise, since both equations result from Eqs. (9) and (12). As a consequence, the relation between Ṁ_max and τ is essentially the same as in Eq. (17), with f = 1; it is written

(24)

As a note of caution, these results have been obtained under the assumption that ξ is constant. There is no reason to believe that this is the case, nor that ξ should not depend on parameters such as M₁, C, or α. The goodness of the fits seems to indicate, however, that ξ is not a sensitive function of these parameters. One reason for this is that, by definition, ξ >  2.47; and, from our experience, Ṁ_tr is never much bigger than Ṁ.

5. Comparison with observations

As we have known for some time, the irradiated-DIM accounts well for the observations of classical galactic transients (Coriat et al. 2012) and provides basic information about their observed light curves (Dubus et al. 2001) but requires refinement if it is to reproduce outbursts that have been really observed (Tetarenko et al. 2018a,b). For accretion rates less than the Eddington value, the luminosity is proportional to ṁ and the decay time of the luminosity τ_L is equal to the decay time of the accretion rate. From Eq. (17), in order to have ṁ_max ≤ 1 we need

(25)

This (taking ψ = 1.3) is in reasonable agreement with Eq. (10) in Lasota et al. (2015) who used f = 3, and found a critical τ_L smaller by a factor of about 1.3. However, Eq. (25) is more general since it explicitly contains f, which is implicitly related to the disc size at maximum.

If, on the other hand, ṁ ≫ 1 so that the apparent luminosity scales approximately as ṁ², the decay time of the luminosity τ_L is half the decay time of the accretion rate, and Eq. (17) can be written as

(26)

The comparison with Eq. (17) of Lasota et al. (2015) is a bit pointless because in both formulae the parameters are raised to rather high but somewhat different powers.

If the entire disc is not brought into the hot state, the total outburst duration is related to the maximum mass accretion rate via Eq. (24), and we must have

(27)

where now refers to the outburst duration, longer than the decay time that enters Eq. (26). Equation (27) does not apply when the entire disc is brought into the hot state, because in this case the outburst duration is longer than the propagation time of a cooling wave and because ṁ_max is higher than the accretion rate for which a cooling wave starts by a factor f, possibly much larger than unity.

These equations confirm the main argument of the present paper, that to study super-Eddington outbursts one has to consider large discs. A detailed comparison of Eqs. (26) and (27) with reality could be attempted and is encouraged, keeping in mind that this has been done already for sub-Eddington outbursts, suggesting that α_h >  0.2 (Tetarenko et al. 2018b) and that f_irr is constant neither in time, nor in space, with possible values larger than 1 (Tetarenko et al. 2018a). These papers, however, used the King & Ritter (1998) exponential-decay formalism so that checking the influence of this assumption on the conclusions about α_h and f_irr is probably worth trying. Also Eq. (17) could be used to check the consistency of the “observed” α_h and f_irr with the disc size. In principle, for transient ULXs, as for LMXBs, it should be possible to determine from observations if a phase during which the entire disc is brought to a hot state exists, because there should be a change of slope in the light curve; the existence of such a break might, however, be difficult to assess if the signal-to-noise ratio is not high enough.

This said, it is worth using the numerical models directly to compare them with the observations of X-ray outbursts with well sampled and not too extravagant light curves that allow to test the credibility of the values of the assumed parameters. We first discuss the case of V404 Cyg, a galactic transient source with an orbital period of 155 hr, before considering more extreme cases, in terms of luminosity and possibly of disc size.

5.1. V404 Cyg

The parameters of model 1 are adequate for a system such as V404 Cyg: M₁ = 7 M_⊙ and P_orb = 155 h. The orbital period of V404 Cyg is about 6.5 d (Casares et al. 1992), and the mass function is f(M) = 6.1 M_⊙, with a mass ratio (Casares & Charles 1994). The secondary is a K2–4 giant (Khargharia et al. 2010), located at a distance of 2.39 ± 0.14 kpc (Miller-Jones et al. 2009), in agreement with the Gaia parallax of 0.44 ± 0.10 mas (Gandhi et al. 2019). Integrating the observed X-ray light curve, Życki et al. (1999) estimated that the total accreted mass during the 1989 outburst was 6 × 10²⁵ g, assuming a radiative efficiency of 0.1 and a distance of 3.5 kpc; rescaling to the currently accepted distance of 2.39 kpc leads to an estimate of the mass transferred during this outburst of 2.6 × 10²⁵ g (Ziółkowski & Zdziarski 2018), implying an average accretion rate of the order of 2.5 × 10¹⁶ g s⁻¹ if the amount of mass accreted in quiescence is negligible.

Figure 1 shows the time evolution of the accretion disc for parameters that are compatible with those of V404 Cyg: the mass transfer rate from the secondary is 5 × 10¹⁶ g s⁻¹; we used α_h = 0.2 and α_c = 0.04, and took η_t𝒞 = 5 × 10⁻⁴, as in Dubus et al. (2001). We used μ = 10³¹ G cm³. We find a sequence of regular outbursts, lasting 11 months, with a peak accretion rate of 1.4 × 10¹⁸ g s⁻¹, and which recur every 12 years. These characteristics compare reasonably well with the observed properties of V404 Cyg, which showed three outbursts in 1938, 1989, and 2015, implying a recurrence time of the order of 30–50 yr. The recurrence time we obtain is shorter than that observed in V404 Cyg; longer recurrence times can easily be obtained by reducing the viscosity in the cold state α_c, as we have seen in Sect. 3.3. The peak luminosity in 1989 was about 8 × 10³⁸ erg s⁻¹, and the outburst lasted for about 400 d. The 2015 outburst was very different from the 1938 and 1989 events (see e.g. Casares et al. 2019). These authors suggest that the outburst was terminated by a very strong outflow from the outer disc regions and invoke the possibility of an enhanced mass transfer during the outburst. As discussed above, winds are not properly taken into account in the standard DIM. Enhanced mass transfer is routinely added to the DIM description of dwarf nova outbursts (see e.g. Hameury 2020), but only rarely invoked in the case of X-ray transient sources (see Esin et al. 2000, for an example).

5.2. Ultraluminous X-ray sources

ULXs are a more extreme case, both in terms of luminosity and of disc size, to which our model can be applied. As stated earlier, presumably only a fraction of these sources are transient; we first consider in detail M51 XT-1 and HLX-1.

5.2.1. M51 XT-1

Figure 10 shows how observational data of M51 XT-1 (Brightman et al. 2020) compare with our model when using Eqs. (13) and (22), combined with Eq. (2) with the beaming parameter to obtain the light curve. We have considered a 1.4 M_⊙ accreting neutron star, and the disc size was taken to be 4.8 × 10¹¹ cm. As can be seen, the agreement is reasonably good and as acceptable as the original fit with a power law with index −5/3. These data can also be fitted with a 10 M_⊙ accreting black hole, with a disc extending to 4.9 × 10¹¹ cm.

Fig. 10. Observed flux from M51 XT-1 as compared with model predictions for a 1.4 M⊙ neutron star (red curve) and a 10 M⊙ black hole (green curve). Black points correspond to Swift data and blue points to Chandra or XMM-Newton observations (courtesy of M. Brigthman, see also Brightman et al. 2020).

In the neutron star case, the maximum accretion rate is 6 × 10¹⁹ g s⁻¹; in the black hole case, it is 1.5 × 10²⁰ g s⁻¹. Hence the emission from a black hole system would not be beamed, but in the case of a neutron star accretor the beaming factor would be b = 0.06, as expected for ṁ = 33. The mass transfer rate can be very roughly estimated to be of the order of 1–2% of the maximum accretion rate, making it 1 − 3 × 10¹⁸ g s⁻¹. We checked that making instead of in Eq. (2) also provides an acceptable fit, but is excluded by the Chandra and XMM-Newton points.

5.2.2. HLX-1 in ESO 243-49

The X-ray source HLX-1 in the galaxy ESO 243-49 (Farrell et al. 2009) is the brightest of the few hyperluminous X-ray sources, which are defined by L_X >  10⁴¹ erg s⁻¹; one of them, NGC 5907 ULX1, with L_X ≳ 10⁴¹ erg s⁻¹ is a pulsating ultraluminous X-ray source (PULXs, see King & Lasota 2020, for their properties and references). The luminosity of HLX-1 is variable and can exceed 10⁴² erg s⁻¹ at maximum. Its association with the galaxy ESO 243-49 at a distance of 95 Mpc is rather well established (Wiersema et al. 2010; Soria et al. 2013) so there is not much doubt about the value of its luminosity (see, however, Lasota et al. 2015). The mass of the presumed accretor is widely believed to be close to or higher than 10⁴ M_⊙, which would make it one of the few intermediate-mass black holes (IMBH) claimed to reside in an ULX. The once almost universal belief that all ULXs host IMBHs (i.e. compact accretors with masses between 10² and 10⁴ M_⊙; Colbert & Mushotzky 1999) has been shattered by the discovery of an X-ray pulsar in M82 ULX-2 by Bachetti et al. (2014), followed by several other discoveries of such pulsing ULXs. The evidence that the compact object in HLX-1 is an IMBH is based mainly on the spectral (Servillat et al. 2011; Godet et al. 2012; Titarchuk & Seifina 2016; Soria et al. 2017) and radio properties (Webb et al. 2012) of this source.

However, since 2008, HLX-1 has been observed to exhibit outbursts that first appeared at about one year intervals; then the recurrence became less frequent and the last outburst was observed in April 2017 (see Lin et al. 2020). The outbursts have the typical fast-rise slow-decay (FRSD) shape of X-ray transient outbursts produced in accretion discs around stellar-mass compact accretors. Although they look like disc-instability X-ray transient events and have the same timescales (decay times of about 30 days, duration 180 days, recurrence times longer than a year or more), they cannot be their analog for an IMBH precisely because of these timescales (Lasota et al. 2011): the outbursts are much too short for a mass of 10⁴ M_⊙. Several models that often invoke episodic accretion from a companion star in a very eccentric orbit have been proposed to explain HLX-1 outbursts assuming that the accretor is an IMBH (see e.g. Lasota et al. 2011; Godet et al. 2014; Lin et al. 2020), but all face severe difficulties both for observational and theoretical reasons (see e.g. Soria et al. 2017; Lin et al. 2020; van der Helm et al. 2016).

Another model, according to which the HLX-1 outbursts would have their origin in a Compton-heated wind instability appearing at wind mass-losses much bigger than the central accretion rate (Shields et al. 1986) has been proposed by Soria et al. (2017). This is an interesting possibility, however, the light curves generated by such an instability do not, and cannot, even remotely resemble what is observed in HLX-1 (see e.g. Shields et al. 1986; Ganguly & Proga 2020).

It is not an accident since, as mentioned above, the outbursts of this hyperluminous source have a shape naturally explained by the thermal–viscous disc instability model. Since no other model, as yet, is able to reproduce such light curves, but the DIM applied to a disc around a 10⁴ M_⊙ black hole produces incorrect timescales (Lasota et al. 2011), it is legitimate to test whether at least the outburst shape, peak luminosities, and characteristic times can be reproduced assuming that, as is probably the case for most ULXs, HLX-1 possesses a stellar-mass accretor.

In 2012, Webb et al. (2014) observed the beginning of an HLX-1 outburst quasi-simultaneously in X-rays and optical, and found that there is a possibility that the flux in the optical V band began to rise about two days before the X-rays. Such a few-days “X-ray delay” is typical of X-ray transient outbursts in LMXBs (Russell et al. 2019) and is explained by the viscous-time filling of the quiescent disc’s inner hole at the outburst start. Because of the short two-day timescale, this cannot be the explanation if the accretor is an IMBH (Webb et al. 2014). In principle this X-ray delay could be explained by propagation at the sound speed of a disturbance in the disc produced by a black-hole-orbiting star, but this would just add another speculative element to models that, as mentioned above, have also to face other difficulties.

Therefore, if the X-ray delay is real and produced in the disc, it would nicely fit with the other properties of the outbursts (shapes, timescales) that are typical of stellar-mass XBTs and there would be no escape from the conclusion that also here we are dealing with a low-mass black hole or a neutron star. Timescales are more basic than interpretation of spectra⁶.

This has motivated Lasota et al. (2015) to propose a model according to which HLX-1 contains a 3 M_⊙ accretor surrounded by an unstable disc whose maximum apparent luminosity 10⁴² erg s⁻¹ corresponds to ṁ = 170 (see Eq. (2)). This solution was obtained in the framework described at the end of Sect. 4.2.1, and we revisit it here.

Assuming that ṁ_max ≫ 1 and there is a peak luminosity of 10⁴² erg s⁻¹, Eq. (2) gives, in the case ,

(28)

which can easily be solved when M₁ is given. If the decay is due to the viscous decay of a disc fully in the hot state, then Eq. (26) must also be satisfied; the decay time varies between 18 and 62 d (Yan et al. 2015). We assume here that τ_L = 40 d; this leads to

(29)

We also require that f does not exceed 50–100 for the disc to be unstable. This is equivalent to setting an upper limit on M₁, as it can easily be seen that increases with increasing M₁. For a neutron star primary with M₁ = 1.4, ṁ_max = 248 and . This is possible if, as in sub-Eddington transients, α_h is significantly larger than 0.2. If, for example, α_h = 0.6 and f_irr = 1, one would have f = 18. This corresponds to a peak accretion rate of 3.8 × 10²⁰ g s⁻¹, Ṁ₂ ≃ 10¹⁹ g s⁻¹, and an outer disc radius of 7.9 × 10¹¹ cm. This is much smaller that the determination by Soria et al. (2017) who found that the outer radius is about 10¹³ cm by fitting the optical emission with an irradiated disc model plus a red, stellar component; however, we do have deep concerns about the accuracy of such models of the optical emission from the system in particular in the presence of strong outflows that may extend to large distances and contribute in a significant way to the optical emission. A solution with f <  100 with a primary black hole would imply unrealistic parameters for the viscosity and the irradiation efficiency. The difference between our conclusion and that of Lasota et al. (2015) is essentially due to the fact that Lasota et al. (2015) used the duration of the outburst instead of the decay time in the ṁ_max − τ relation; given that τ_L enters Eq. (26) to the power 3.39, the impact is quite significant.

If, on the other hand, the observed decay is due to the propagation of a cooling front that has not (quite) reached the outer disc edge, then Eq. (26) does not apply, but one should use instead the relation between the maximum accretion rate and outburst duration (Eq. (27)), with d. As for the case of the decay of a fully hot disc, the higher M₁, the higher is f. Assuming again a 1.4 M_⊙ accreting neutron star, we find that the 180 d duration can be explained with α_0.2 = 3 and f_irr = 0.5 if ξ = 7.7.

The amplitude of HLX-1 outbursts is about 50, too low for a standard X-ray transient; it is not clear that the X-rays observed in quiescence are all produced by the disc, but it could also be “fitted” by increasing α_c, which at the same time would account for the relatively short recurrence time. And there would be no problem explaining and reproducing a two-day X-ray delay. We conclude that a stellar-mass accretor is still a viable (but admittedly not compelling) option for HLX-1 in ESO 243-49 (see also King & Lasota 2014), which should be seriously considered since no model assuming the presence of an IMBH is able to explain the basic properties of the observed outbursts. More work is obviously needed to consolidate this option and to account for all of the observed spectral and optical properties of HLX-1, but this is a formidable task that is clearly outside the scope of this paper.

5.2.3. Other transient sources

It is also interesting to note that the DIM predicts that even if the mass transfer rate from the secondary is constant, sequences of outbursts with very different outburst peak luminosities and fluences can be produced. NGC 925 ULX-3 is a good example of such a situation (Earnshaw et al. 2020): this source was detected with a luminosity of (7.8 ± 0.8)×10³⁹ erg s⁻¹ in November 2017 by Chandra; it had been detected in 2011 at a similar luminosity by Swift, and also in January 2017 by XMM-Newton, albeit at a level twenty times lower. Although the January 2017 observation does not necessarily coincide with the peak of the outburst, such a possibility should be kept in mind.

5.2.4. Steady sources

We suggest that some of the apparently steady ULXs could well be transient sources caught in long-lasting outbursts. As shown earlier, outbursts can last for years. Several ULXs have been observed for a long time, and have not shown any sign of long term variability, but this option cannot be completely excluded since the DIM can easily produce outbursts lasting for ten years (see e.g. models 22, 23, 24, 24a) or even more if the disc is larger than what we have considered here and the mass transfer rate is sufficient for the heating front to reach the outer disc edge, or if the viscosity is small (see Eqs. (15) and (23)). It is unclear whether any of the observed steady sources would fall into this category, but this option should be kept in mind when trying to account for the observed luminosity of ULXs, which might not directly reflect the mass transfer rate from the secondary.

6. Conclusions

We successfully extended and generalized the irradiated-DIM to the case of large accretion discs and super-Eddington accretion rates. Assuming that, during decay, the inner disc extending between the inner edge and the position of the cooling front is close to being steady with an accretion rate equal to the critical rate, we have been able to derive a relation between the peak accretion rate during an outburst and the maximum distance reached by the heating front, that closely matches the results of numerical simulations as long as the heating front does not reach the outer edge of the disc. We have also been able to solve explicitly the time evolution of the outburst decay, which consists of one or two phases: the propagation of the cooling front throughout the disc, possibly preceded by a phase during which the full disc is in the hot state. In both cases, we have been able to determine the characteristic timescales. Again, these analytic solutions are in very good agreement with the results of numerical simulations. These results are important, because they limit the need for numerical simulations.

We have shown that, in most cases, the peak accretion rate obtained during outbursts is proportional to the mass transfer rate, with a proportionality coefficient that depends mainly on the ratio α_h/α_c, and is of the order of a few tens, possibly reaching 100 when α_h/α_c = 10. We found that sub-Eddington outbursts of systems with large orbital periods, such as V404 Cyg, are well accounted for by the DIM. We have shown that, provided that the fraction of the X-ray flux that irradiates the accretion disc increases when the X-ray luminosity exceeds the Eddington luminosity and hence the accretion efficiency decreases, super-Eddington accretion rates can be obtained in large accretion discs during outbursts. We have shown that ULX transients (ULXTs), such as M51 XT-1 are well described by the DIM with parameters consistent with the general properties of such sources. Some of these outbursts might last several years, and could thus account for some other types of variability of observed ULXs. Outbursts of HLX-1 in ESO 243-49 can be described by our model, on the condition that the accretor is a neutron star.

In our models, a number of transient ULXs are assumed to be beamed. This certainly must be the case for ULXs containing neutron stars, since luminosities exceeding up to one hundred times the Eddington luminosity cannot be attained if L_x increases as ln ṁ. Other options have been suggested, such as the presence of a magnetar-type magnetic field, but they are not all viable (see King & Lasota 2020). Beaming implies that a large number of sources exist that would be seen off axis, and therefore would not appear as ULXs, although they would be surrounded by a nebula; an example of such a source is SS 433, which is most probably an ULX seen from the side.

We also note that not all ULXs need to be beamed since according to Eq. (2) this should only be the case for . Indeed, using simulations of stellar populations, Wiktorowicz et al. (2019) have shown that observed ULXs harbouring black hole accretors are typically emitting isotropically, whereas systems containing neutron stars would be beamed. This is in line with observational arguments showing that, at least in some systems, the luminosity averaged over solid angle is not very different from the apparent luminosity. For example, a comparison of the luminosities from the soft X-rays and from the helium line emitted by the surrounding nebula suggests that the beaming effect is at best moderate in a few sources (Pakull & Mirioni 2002; Moon et al. 2011; Russell et al. 2011); NCG 300 ULX-1, which contains a neutron star, is in a similar situation (Binder et al. 2018), but since its apparent luminosity is 4.4 × 10³⁹ erg s⁻¹ it corresponds according to b ≈ 0.2 (Binder et al. 2018; King & Lasota 2019). However, bright PULXs with > 10⁴⁰ erg s⁻¹ must be strongly beamed since no physically motivated alternative exists.

We finally note that in our model a potentially crucial ingredient, namely winds, is included only implicitly and in a very crude way. Winds are required here to avoid the decrease in the illumination factor f_irr when the accretion rate exceeds the Eddington value, but they also significantly modify the outburst duration and could have a major impact on angular momentum transport if the wind is magnetized. Unfortunately, although they have an important effect on the whole outburst cycle, the physics of outflows from accretion discs are still not fully understood (Dubus et al. 2019; Tetarenko et al. 2020).

Thus the DIM applied to ULXTs suffers from the same weaknesses as the “standard” DIM used to describe LMXB transient sources (some are even enhanced). However, the latter model has the advantage of being based on sound physical assumptions and could, in principle, be easily tested by observations, although in practice the large distances to ULXs are an obvious drawback.

Acknowledgments

We thank Murray Brightman for having provided us with the observational data points of M51 XT-1. We are grateful to Guillaume Dubus, Andrew King, and Matt Middleton for helpful comments and suggestions. J.P.L. was supported by a grant from the French Space Agency CNES.

References

All Tables

All Figures

	Fig. 1. Time evolution for model 1, which corresponds to parameters appropriate for V404 Cyg. Upper panel: disc mass. Intermediate panel: mass accretion rate. Bottom panel: outer disc radius (red curve), inner disc radius (blue curve), and transition radius rtr (black curve). For such sub-Eddington outbursts, mass-accretion variations represent the X-ray light-curve.
In the text

	Fig. 2. Same as Fig. 1 for model 12.
In the text

	Fig. 3. Detailed view of a large outburst in Fig. 2. Top panel: bolometric light curve (red) and accretion rate times 0.1c2. Bottom panel: outer disc radius (red curve), inner disc radius (blue curve), and transition radius rtr (black curve).
In the text

Fig. 4. Radial surface-density profiles of the accretion disc just prior to outbursts for model 12 (solid lines). The curve labelled 1 shows Σ prior to the first weak outburst after a major one, curve 2 shows Σ before the second outburst, curve 3 is just before the last weak outburst of the sequence, and curve 4 represents Σ at the onset of the last outburst. Dashed lines show the critical surface densities Σmin and Σmax of the S-curves.

In the text

	Fig. 5. Time evolution for model 24, in which the irradiation efficiency firr decreases when the accretion rate reaches the Eddington value.
In the text

	Fig. 6. Quiescent radial surface-density profile 173 days after the end of an outburst for model 26. The dashed lines correspond to critical surface densities Σmax (upper line) and Σmin (lower line).
In the text

Fig. 7. Relation between the maximum mass accretion rate reached during an outburst multiplied by and the maximum distance travelled by the heating front. Stars correspond to black hole systems, circles to neutron star systems. Red symbols correspond to an orbital period of 155 h, green symbols to 400 h, and blue symbols to 1200 h. Symbols are filled when the heating front reached the outer edge of the disc. The dashed line represents the critical mass transfer rate as given by Eq. (9).

In the text

	Fig. 8. Time evolution of model 14. Bottom panel: transition radius (black curve) as well as the outer disc radius (red curve). Top panel: actual mass accretion rate (black curve) and critical rate (green curve). The blue curve gives the analytic solutions (Eq. (13)) when the entire disc is in the hot state, and the red curve corresponds to a propagating cooling front (Eq. (22)).
In the text

	Fig. 9. Comparison between the outburst duration and the duration of the propagation time of the cooling wave as estimated by Eq. (22) with ξ = 5. The symbols have the same meaning as in Fig. 7.
In the text

	Fig. 10. Observed flux from M51 XT-1 as compared with model predictions for a 1.4 M⊙ neutron star (red curve) and a 10 M⊙ black hole (green curve). Black points correspond to Swift data and blue points to Chandra or XMM-Newton observations (courtesy of M. Brigthman, see also Brightman et al. 2020).
In the text