© A. Mastichiadis et al. 2022

1. Introduction

Black-hole X-ray binaries (BHXRBs) are systems in which mass is transferred from a donor star onto a stellar-mass black hole, converting its gravitational energy into radiation. These systems spend most of their lives in quiescence but occasionally undergo major outbursts that may reach the Eddington limit (see, e.g., Motta et al. 2017; Miller-Jones et al. 2019; Koljonen & Tomsick 2020, for the recent outburst of V404 Cyg). Studies of the outburst–quiescence cycle of BHXRBs have revealed dramatic changes in their X-ray spectral characteristics and temporal properties in addition to changes in their luminosity (Remillard & McClintock 2006).

It is well known that, during their outbursts, BHXRBs exhibit some characteristic spectral and temporal states, which, according to Belloni (2010), are: the quiescent state, the hard state, the hard-intermediate state, the soft-intermediate state, and the soft state. A different classification scheme was introduced by Remillard & McClintock (2006), and a comparison of the two schemes was presented by Motta et al. (2009). For the purposes of this paper, we will use the Belloni (2010) classification.

The standard picture of BHXRBs in the hard and hard-intermediate states involves an accretion disk and a hot corona around (or above) the black hole. The accretion disk emits soft X-ray photons as part of its multi-temperature black-body spectrum. A fraction of the disk soft radiation is directly seen by the observers, while the rest goes into the corona, gets up-scattered by the hot electrons residing therein, and results in a hard X-ray power-law spectrum. Part of the hard photon radiation can be observed directly, while the rest is either absorbed by the disk and re-emitted as soft secondary photons or reflected. In the hard state, as the name suggests, the X-ray emission is dominated by the hard component, with a very faint thermal disk component sometimes identified.

While the presence of a corona in the hard state is central in describing the phenomenology of BHXRBs, its nature remains largely uncertain. For example, in some models (for a review, see Done et al. 2007) the corona is taken to be the inner part of the accretion flow around the black hole, which in the hard and hard-intermediate states is geometrically thick, optically thin, and hot (Narayan & Yi 1994, 1995). Alternatively, the corona might be the base of an outflow that emanates from the hot inner flow (Markoff et al. 2005) or even the entire outflow, since photons entering at its base can travel through the length of the whole outflow (Reig et al. 2003; Reig & Kylafis 2021).

A unique feature of BHXRBs is their time variability. The power spectra of BHXRBs exhibit characteristic frequencies, called quasi-periodic oscillations (QPOs; for a recent review, see Ingram & Motta 2019). When BHXRBs are in the so-called hard or hard-intermediate state, where the energy spectrum is dominated by the hard power law, the observed QPOs are called type-C (Casella et al. 2005). In the soft state, where the energy spectrum is dominated by the soft multi-temperature black body emitted from the disk, QPOs are rarely detected (see, e.g., Motta et al. 2012). The frequency of type-C QPOs in BHXRBs is in the range of a few tens of millihertz to a few tens of hertz, and it is an increasing function of luminosity (Belloni et al. 2005). It also correlates with spectral parameters, such as the photon index of the power law (Vignarca et al. 2003). Several scenarios have been proposed for explaining type-C QPOs, including oscillations of boundary layers (Titarchuk & Fiorito 2004), some type of disk instabilities (Tagger & Pellat 1999), or Lense-Thirring precession (Stella & Vietri 1998; Ingram et al. 2009).

The idea that there is an interaction between the corona and the accretion disk was first put forward by Haardt & Maraschi (1991) to explain the X-ray spectra from radio-quiet active galactic nuclei (AGN). Usually, the implicit assumption is that there is a balance between heating and cooling, and the corona is thereby in some stationary state. Even when variations in one or more physical parameters are considered, the changes are assumed to occur in such a way that this balance is not violated. However, there is no a priori reason why this should always be the case. This issue could be clarified only with a study of the dynamical behavior of the corona-disk system, and such research, to the best of our knowledge, had not been performed until now. The aim of the present paper, therefore, is to make an exploratory study of the dynamical aspects of the radiative feedback between the corona and the accretion disk in the hard (and hard-intermediate) state of BHXRBs.

To treat the above nonlinear problem, we write two coupled time-dependent equations that describe the evolution of the energy density of the hot electrons in the corona and the energy density of the hard photons emanating from it. The electrons are cooled due to the Comptonization of soft photons from the disk, including those produced by the reprocessing of the impinging hard photons on the disk. While these equations are by necessity simplified – they consider only the integrated energy densities of each population – they are self-consistent, in the sense that they conserve energy, and they keep the physical essence of the interactions. Thus, for instance, the hot electrons in the corona are assumed to be heated by mechanisms that are left unspecified, but are connected with the accretion flow. On the other hand, since the electrons interact with the disk photons, they cool through inverse Compton scattering, and the electron energy losses are taken to be equal to the radiated hard photon luminosity.

The corona-disk system, as described above, has many similarities to the prey-predator dynamic models commonly used in the study of biological and ecological systems (e.g., Brauer & Castillo-Chávez 1996). The prey-predator relationship basically describes the interactions between two species and their mutual effects. For example, in the case of the lynx and the snowshoe hare, when the number of predators (lynx) is low, the preys (hares) thrive. The population of predators, however, will increase as they can reproduce and survive by feeding on the increasing number of preys. With more lynx hunting, the hare population will rapidly decline, eventually bringing down the population of lynx due to starvation, and the cycle repeats itself. In the corona-disk system under study, hot electrons are the prey and soft disk photons (from the disk and reprocessed) are the predator. When the energy density of the latter is low, the electron energy density builds up in the corona because the inverse Compton scattering cooling is negligible. When the electrons “thrive”, the energy density of hard photons impinging on the disk also increases, thus leading to an increase in the reprocessed soft radiation and electron cooling. This, in turn, suppresses the “growth” of the electron energy density, eventually leading to a decrease in the electron cooling agents, and the cycle repeats.

Indeed, the solution of the ordinary differential equations that describe the corona-disk coupling shows that the energy densities of the interacting particle species undergo quasi-periodic changes before they relax to a steady state, that is to say, the system undergoes damped oscillations. The frequency of the oscillations depends on the mass accretion rate, the radius of the corona, the fraction of hard radiation that is being reprocessed to soft, and the black hole mass. Interestingly, for nominal values of these parameters, the derived frequencies are in the range of type-C QPOs (i.e., from a millihertz to tens of hertz).

This paper is structured as follows. In Sect. 2 we outline the model and present the equations that describe the interaction between the electrons in the corona and the soft radiation from the disk. In Sect. 3 we present numerical solutions to the equations of the corona-disk system and investigate the effect of the model parameters on its variability properties. We conclude in Sect. 4 with a discussion of our results.

2. Model description

We first consider a corona with a thermal electron population of temperature Θ = k_BT/m_ec² immersed in a low-energy (soft) photon bath with energy density U_s. The cooling timescale of the electrons due to Comptonization is t_cool = m_ec²(4σ_TcU_s)⁻¹ (e.g., Rybicki & Lightman 1986) and depends solely on the energy density of soft photons. If the electron heating rate is independent of U_s, as assumed here, it is likely that cooling might not be able to balance heating at all times, and a steady state might not be reached for all conditions.

The irradiation of the disk by the hard coronal photons results in a reflected component and an absorbed one, which is reprocessed. The reprocessed radiation is assumed to be soft, and as such, it also contributes to the cooling of the electrons. Since the reprocessed radiation is originally produced by the electrons, it constitutes a nonlinear cooling term in analogy to the well-known synchrotron self-Compton process for relativistic electrons (Bloom & Marscher 1996; Mastichiadis & Kirk 1997; Zacharias & Schlickeiser 2010).

The study of the nonlinear cooling of the electrons in the corona and its impact on the emergent radiation (hard and reprocessed components) requires a time-dependent approach. In this work we construct a simplified set of equations that describe the interplay between: (i) hot electrons in the corona, which we denote by the subscript e, (ii) hard photons produced through Comptonization, which we denote by the subscript h, (iii) soft photons produced from the disk thermal emission, which we denote by the subscript ds, and (iv) soft photons produced via the reprocessing of the irradiating hard component in the disk, which we denote by the subscript rs.

Our simple model does not consider the evolution of the energy dependence of the particle densities, only the energy-integrated densities of all components. A similar approach, which has been adopted for the study of coupled nonthermal radiative processes in leptohadronic plasmas, was proven successful at capturing the main dynamical properties of the system (Petropoulou & Mastichiadis 2012). In the following, we write the equations and relations that describe the temporal evolution of each component.

Hot electrons. We assume that a fixed number of electrons N_e resides in the corona, which is modeled as a sphere of radius R_c. The electrons are heated by some mechanism, which is left unspecified (for heating models, see, e.g., Tout & Pringle 1992; Miller & Stone 2000; Beloborodov 2017; Yuan et al. 2019; Sridhar et al. 2021), but it is assumed to depend on the mass accretion rate Ṁ. More specifically, we assume that the energy injection rate to hot electrons is a fraction f ≤ 1 of the rate with which gravitational energy is being released as mass is being accreted from infinity to the radius of the innermost stable orbit (i.e., fGM_∗Ṁ/R_isco; e.g., Haardt & Maraschi 1991). Here, M_* is the mass of the compact object, R_isco = ξR_S is the radius of the inner stable circular orbit, R_S = 2GM_*/c² is the Schwarzschild radius of the compact object, and ξ = 3 or 1/2 for nonrotating or maximally (prograde) rotating compact objects, respectively (Kaplan 1949; Bardeen et al. 1972). We note that the specific choice for the heating rate would not alter our main findings as long as it is proportional to Ṁ.

The electron cooling rate (in the limit of small optical depths) can be written as −N_eΘm_ec²/t_cool (e.g., Rybicki & Lightman 1986). Then, the equation that describes the energy balance of electrons reads

$$\begin{array}{c}\hfill \frac{\mathrm{d}{E}_{\mathrm{e}}}{\mathrm{d}t}=f\frac{G{M}_{\ast }\dot{M}}{\xi R_{\mathrm{S}}}-4N_{\mathrm{e}}\mathrm{\Theta }{\sigma }_{\mathrm{T}}c{U}_{\mathrm{s}},\end{array}$$(1)

where E_e ≈ (3/2)N_eΘm_ec² is the total energy content of the electrons in the corona. Because the number of electrons is considered fixed inside the corona, there is no escape term (−E_e/t_e, esc) on the right-hand side of the above equation. We discuss the effects of particle escape in Sect. 4.

Dividing by the volume of the corona, Eq. (1) can be written in terms of the thermal electron energy density, U_e, as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{U}_{\mathrm{e}}}{\mathrm{d}t}=\frac{3f}{4\pi }\frac{G{M}_{\ast }\dot{M}}{\xi R_{\mathrm{S}}{R}_{\mathrm{c}}^3}-\frac{8}{3}\frac{{\sigma }_{\mathrm{T}}c}{m_{\mathrm{e}}{c}^2}{U}_{\mathrm{s}}{U}_{\mathrm{e}}.\end{array}$$(2)

Comptonized (hard) photons. These photons are produced by the hot electrons in the corona. Their energy production rate should be equal to the electron loss rate, and their escape from the source is given, in the limit of optically thin coronae, by the source crossing time t_c = R_c/c (see, e.g., Fabian et al. 1986). The equation that describes the hard photon energy density can be written as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{U}_{\mathrm{h}}}{\mathrm{d}t}=\frac{8}{3}\frac{{\sigma }_{\mathrm{T}}c}{m_{\mathrm{e}}{c}^2}{U}_{\mathrm{s}}{U}_{\mathrm{e}}-\frac{c}{R_{\mathrm{c}}}{U}_{\mathrm{h}}.\end{array}$$(3)

Disk (soft) photons. These photons are the soft photons produced at the inner part of the cold accretion disk. For simplicity (see also Karpouzas et al. 2021), they are assumed to be produced within a sphere of radius R_d by the dissipation of the remaining 1 − f fraction of the accretion power. In our formalism, disk photons do not depend on the other populations. Therefore, their energy density just outside this sphere can be written directly as

$$\begin{array}{c}\hfill U_{\mathrm{ds}}=\frac{3(1-f)}{4\pi }\frac{G{M}_{\ast }\dot{M}}{\xi R_{\mathrm{S}}{R}_{\mathrm{d}}^{2}c}.\end{array}$$(4)

Reprocessed (soft) photons. These photons are produced by the hard photons, which illuminate the accretion disk and get absorbed there. Reprocessing in the disk transforms the incident hard radiation to soft. For simplicity, we neglected the slight temperature increase in the disk and assumed that the reprocessed photons have the same average energy as the disk photons. Thus, as long as the production rate of the reprocessed soft photons is proportional to the incident rate of hard photons, one can relate their energy densities as

$$\begin{array}{c}\hfill U_{\mathrm{rs}}=a{U}_{\mathrm{h}},\end{array}$$(5)

where a < 1 is a factor that incorporates both geometrical and albedo effects. Since both the disk and the reprocessed photons are soft, they contribute to electron cooling¹ and one can write

$$\begin{array}{c}\hfill U_{\mathrm{s}}=U_{\mathrm{rs}}+f_{\mathrm{dil}}{U}_{\mathrm{ds}},\end{array}$$(6)

where f_dil = (R_d/R_c)², if R_d ≤ R_c, is a geometric factor that accounts for the dilution of the soft disk photon density in the corona.

It is useful to introduce normalized quantities at this point, namely $\stackrel{\sim }{t}=ct/R_{\mathrm{S}}$, r_c, d = R_c, d/R_S, and ${\stackrel{\sim }{u}}_{\mathrm{i}}=4\pi R_{\mathrm{S}}^{2}c{U}_{\mathrm{i}}/L_{\mathrm{Edd}}$, where L_Edd = 4πGMm_pc/σ_T and ṁ = Ṁ/Ṁ_Edd = Ṁc²/L_Edd. Then, Eqs. (2)–(5) can be rewritten in dimensionless form as

$$\begin{array}{cc}& \frac{\mathrm{d}{\stackrel{\sim }{u}}_{\mathrm{e}}}{\mathrm{d}\stackrel{\sim }{t}}=\frac{3f}{2\xi }\frac{\dot{m}}{r_{\mathrm{c}}^3}-\frac{4}{3}\frac{m_{\mathrm{p}}}{m_{\mathrm{e}}}{\stackrel{\sim }{u}}_{\mathrm{e}}({\stackrel{\sim }{u}}_{\mathrm{rs}}+{\stackrel{\sim }{u}}_{\mathrm{ds}}{f}_{\mathrm{dil}}),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill & \frac{\mathrm{d}{\stackrel{\sim }{u}}_{\mathrm{h}}}{\mathrm{d}\stackrel{\sim }{t}}=-\frac{{\stackrel{\sim }{u}}_{\mathrm{h}}}{r_{\mathrm{c}}}+\frac{4}{3}\frac{m_{\mathrm{p}}}{m_{\mathrm{e}}}{\stackrel{\sim }{u}}_{\mathrm{e}}({\stackrel{\sim }{u}}_{\mathrm{rs}}+{\stackrel{\sim }{u}}_{\mathrm{ds}}{f}_{\mathrm{dil}}),\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill & {\stackrel{\sim }{u}}_{\mathrm{ds}}=\frac{3(1-f)}{2\xi }\frac{\dot{m}}{r_{\mathrm{d}}^2},\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill & {\stackrel{\sim }{u}}_{\mathrm{rs}}=a{\stackrel{\sim }{u}}_{\mathrm{h}}.\hfill \end{array}$$(10)

The steady-state energy densities can easily be obtained from Eqs. (7) and (8) after setting the time derivatives equal to zero. They read

$$\begin{array}{cc}& {\stackrel{\sim }{u}}_{\mathrm{h}}^{(\mathrm{ss})}=\frac{3f}{2\xi }\frac{\dot{m}}{r_{\mathrm{c}}^2},\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill & {\stackrel{\sim }{u}}_{\mathrm{e}}^{(\mathrm{ss})}=\frac{3}{4r_{\mathrm{c}}}\frac{m_{\mathrm{e}}}{m_{\mathrm{p}}}\frac{f}{1-f+af}\hfill \end{array}$$(12)

or, in dimensional quantities,

$$\begin{array}{cc}& U_{\mathrm{h}}^{(\mathrm{ss})}=\frac{3f}{\xi }\frac{\dot{M}{c}^2}{8\pi c{R}_{\mathrm{c}}^2},\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill & U_{\mathrm{e}}^{(\mathrm{ss})}=\frac{3}{8}\frac{m_{\mathrm{e}}{c}^2}{{\sigma }_{\mathrm{T}}{R}_{\mathrm{c}}}\frac{f}{1-f+af}·\hfill \end{array}$$(14)

Interestingly, the steady-state energy density of hot electrons does not depend on the mass accretion rate but solely on the size of the corona and a (for f ≈ 1). Equation (14) can also be expressed in terms of the electron temperature and Thomson optical depth (τ_T = n_eσ_TR_c), leading to τ_TΘ = f/(4(1 − f + af)) ≈ 1/(4a) for f ≈ 1. Thus, for τ_T ∼ 1 and a ∼ 0.5, our simple model predicts kT_e ∼ 250 keV. It can also easily be shown that the hard photon luminosity at steady state, $L_{\mathrm{h}}=4\pi R_{\mathrm{c}}^{2}c{U}_{\mathrm{h}}^{(\mathrm{ss})}/3$, is equal to the energy dissipation rate at R_isco that goes into electron heating in the corona: fṀc²/(2ξ).

A time-dependent analysis of Eqs. (7) and (8) would reveal how different populations reach their steady-state values. This behavior can be qualitatively understood in two limiting cases. If ${\stackrel{\sim }{u}}_{\mathrm{ds}}\gg {\stackrel{\sim }{u}}_{\mathrm{rs}}$, then Eq. (7) reduces to a linear differential equation of the electron density, since the disk photon density is solely determined by the accretion rate. In this regime, the hot electron energy density reaches a steady state without exhibiting any intrinsic variability. In the case, however, where ${\stackrel{\sim }{u}}_{\mathrm{ds}}\ll {\stackrel{\sim }{u}}_{\mathrm{rs}}$, then Eqs. (7) and (8) constitute a nonlinear system of differential equations, which has, as we will show, some interesting temporal properties.

3. Results

3.1. Constant mass accretion rate

To better understand the properties of the dynamic corona-disk system under study, we first solved² Eqs. (7) and (8), assuming that both energy densities are initially very small, namely ${\stackrel{\sim }{u}}_{\mathrm{e}/\mathrm{h}}(\stackrel{\sim }{t}=0)={10}^{-10}$, and that the mass accretion rate remains constant. Motivated by observations of accreting systems in the hard state, where the emission is dominated by the hard component, we used f = 0.999. We verified that our results are not sensitive to the exact value of f as long as it is close to unity (for more details, see later in this section). In what follows, we set ξ = 3 and r_d = 10, unless otherwise stated.

Figure 1 shows the temporal evolution of the dimensionless energy density of the various components for different choices of the accretion rate ṁ, the corona radius r_c, and the parameter a. The solutions for all components (except the disk soft photons) are, in general, oscillatory for several cycles before they die out and the system settles into the steady state described by Eqs. (9), (12), and (11). The frequency and amplitude of these transient oscillations depend on the physical parameters, as exemplified in Fig. 1. For instance, an increase in the accretion rate leads to a higher disk photon density, which acts as a stabilizing factor and drives the electrons in the corona into a steady state faster (see the qualitative discussion at the end of the previous section). Similarly, the time window of the transient oscillations decreases for smaller coronae, while the frequency of the oscillations increases (compare panels a and b). As the parameter a decreases, the ratio of the reprocessed soft photon energy density to the hard incident photon energy density decreases. As a result, the cooling of electrons due to the reprocessed soft photon radiation becomes less important, thus weakening the nonlinear coupling between the electron and hard photon populations (see Eqs. (7) and (8)). Therefore, solutions for lower a values do not exhibit pronounced oscillations, and the system reaches a steady state (compare panels a and d).

Fig. 1. Temporal evolution of energy densities for different model parameters: ṁ = 5 × 10−4, rc = 30, a = 0.2 (panel a), ṁ = 5 × 10−4, rc = 15, a = 0.2 (panel b), ṁ = 5 × 10−3, rc = 30, a = 0.2 (panel c), and ṁ = 5 × 10−4, rc = 30, a = 0.05 (panel d). In all cases, f = 0.999, ξ = 3, rd = 10, and time is normalized to RS/c.

Figure 2 illustrates the prey-predator behavior of the dynamic system under study. Whenever the reprocessed soft photon density increases, the electron energy density decreases because of Compton cooling. Hence, hot electrons represent the prey, and soft photons are the predators. In contrast to the classical prey-predator system, where the oscillations (or limit cycles in phase-space plots) are constant in time, the disk-corona system is always led to a steady state, which is represented by a point in the phase-space plot. Nevertheless, even small perturbations from this steady state lead the system back into an oscillatory mode. We address this point more in the next section and provide analytical arguments in Appendix A.

Fig. 2. Phase-space plot illustrating the dynamic behavior of the corona-disk system. Here, the same parameters as those in panels a, b, and d of Fig. 1 are used. For a better display, the energy densities are normalized to a peak value of one.

Figure 3 illustrates the dependence of the system’s intrinsic variability on ṁ, r_c, a, and f. Here, we display the energy density of hard photons, but similar behavior is also found for the reprocessed soft radiation component and the hot electron population (as previously shown in Fig. 1). The corona-disk system exhibits oscillatory behavior for a wide range of parameter values, except for f. In this specific example, we find oscillations only for f ≳ 0.997. However, we typically find that the system reaches a steady state without signs of intrinsic periodicity if the energy transferred to the corona is ≲99 per cent. Therefore, our simple model predicts no periodicity in soft-state accreting systems where much less power is transferred to the corona. We also notice that, as expected, the steady-state value of the hard photon density does not depend on a, while it scales linearly with ṁ and $r_{\mathrm{c}}^{-2}$ (see also Eq. (11)).

Fig. 3. Temporal evolution of the dimensionless hard photon energy density for different values of ṁ (panel a), rc (panel b), a (panel c), and f (panel d). Unless stated otherwise, the parameters used are: ṁ = 5 × 10−3, rc = 30, a = 0.05, and f = 0.999. Time is normalized to RS/c.

Figure 4 shows the light curve of hard Comptonized photons (in linear scale and after rescaling to a peak value of 100) for the same parameters used in panel a of Fig. 1. For the time units, we assumed a black hole with mass M_* = 10 M_⊙, which corresponds to R_S/c ≃ 9 × 10⁻⁵ s. The raw power spectrum of the time series (multiplied by frequency) is displayed in the right panel of the same figure (in blue). The average power spectrum computed after splitting the time series into m = 5 segments is overplotted (in orange). The peak at ∼3.5 Hz (dashed line) corresponds to the natural frequency of the system (the first harmonic is also visible). This feature resembles a QPO, which is loosely defined as a resolved peak in the power spectrum (van der Klis 1989). In reality, the natural frequency may be hidden in the periodogram depending on the noise level and the light curve duration.

Fig. 4. Evolution of the hard photon density in time and frequency domains. Panel a: light curve of hard (Comptonized) photons (rescaled to a maximum value of 100) for the same parameters used in panel (a) of Fig. 1. Panel b: power spectrum of the light curve multiplied by frequency using the Leahy normalization. Overplotted in orange is the averaged power spectrum computed using five segments of the light curve.

Next we computed the power spectra of noisy light curves, the same way we did for the one displayed in Fig. 4, for different values of the accretion rate, corona radius, and parameter a. In all cases, where at least three peaks could be identified in the light curve of hard photons with a peak value of at least 5% of the asymptotic steady-state value, we determined the natural frequency (f₀) of the system. In Fig. 5, we plot f₀ as a function of the dimensionless accretion rate for different values of r_c and a (as indicated in the color bars of panels a and b, respectively). The oscillations die out, meaning that the system reaches a steady state after exhibiting one or two bursts above a certain value of the mass accretion rate (as indicated by the dotted black lines). We can see that for low enough values of the a parameter, the oscillations also disappear below a certain ṁ value. For this example, in particular, we find no oscillatory solutions for a ≲ 0.02. For a black hole with mass 10 M_⊙, as assumed here, we find that the natural frequency of the system lies in the range between a hertz and tens of hertz and depends on the parameters as $f_0\propto {\dot{m}}^{1/2}{r}_{\mathrm{c}}^{-3/2}{a}^{1/2}$. This result is robust – it does not depend on our method for numerically determining f₀. In fact, the above relation can be obtained analytically by performing a stability analysis of Eqs. (7) and (8) near the steady state (for more details, see Appendix A). This analysis also reveals the dependence on ξ and M_*, which were not varied in the numerical investigation. The full expression for the natural frequency then reads

$$\begin{array}{c}\hfill f_0\simeq 5\mathrm{Hz}{\left(\frac{a}{0.2}\right)}^{1/2}{\left(\frac{\dot{m}}{{10}^{-3}}\right)}^{1/2}{\left(\frac{r_{\mathrm{c}}}{30}\right)}^{-3/2}{\left(\frac{\xi f}{3}\right)}^{-1/2}{\left(\frac{M_{\ast }}{10M_{\odot }}\right)}^{-1}.\end{array}$$(15)

Fig. 5. Natural frequency of the dynamic corona-disk system for a nonrotating black hole with mass M* = 10 M⊙ as a function of the mass accretion rate (in units of the Eddington rate) for different values of rc (panel a) and a (panel b), as indicated in the color bar. Beyond (or below) a certain mass accretion rate, the system reaches a steady state without exhibiting oscillations (dotted black lines). Other parameters used are: a = 0.2 (panel a) and rc = 30 (panel b), with f = 0.999 in both panels.

Observations of BHXRBs in the hard and hard-intermediate states indicate that the QPO frequency increases with X-ray luminosity or, equivalently, the mass accretion rate (e.g., Belloni et al. 2005; Heil et al. 2011). In addition, the corona shrinks in the hard-intermediate state (Kylafis & Reig 2018; Kara et al. 2019), also resulting in an increase in the QPO frequency according to Eq. (15). Finally, as the source evolves in the hard-intermediate state, the inner radius of the accretion disk decreases (García et al. 2015), resulting in higher values of the parameter a and higher f₀ values. It is intriguing that the dependence of the predicted natural frequency on physical parameters in our simple model agrees qualitatively with observed trends.

3.2. Variable mass accretion rate

In the previous section, we demonstrated that the corona-disk system can exhibit quasi-periodicity in the case of a constant accretion rate and under specific initial conditions. To further examine whether these transient oscillations are an intrinsic property of the system, we performed runs in which we allowed the mass accretion rate to vary in time.

We generated Gaussian-distributed noise with a power-law spectrum (∝f⁻²) using the Python software colorednoise³, which is based on the algorithm described in Timmer & Koenig (1995). The generated red-noise time series, {y_i}, is then used to compute a time series for the accretion rate as ṁ>_i = (y_i + δ)^ϵ, for i = 1, …, N. Here, N = 2¹⁶, δ was chosen to ensure positive values of {ṁ_i}, and ϵ was chosen such that the maximum and minimum values of the input always differ by a factor of 2. This normalization allowed us to control the dynamic range of the variations in ṁ and gives us a basis to compare the output photon light curves. We verified that the power spectrum of {ṁ_i} is a power law with the same index as the generated time series {y_i}.

Since the timescale of the mass accretion variations, Δτ = t_i + 1 − t_i, is not known a priori, we studied the system using Δτ = 2^k R_S/c with k = 0, 1, 2, .... For any value of k > 0, we performed a linear interpolation between successive values (i.e., ṁ_i+1 and ṁ_i) and solved Eqs. (7) and (8) with the accompanying relations for the disk and reprocessed soft photon densities. While we started with the same initial conditions as before, here we discarded the first few thousand time steps of the computed light curves in order to remove any dependence on them. When k is small, the variations in ṁ will be faster than the light-crossing time of the corona, therefore effectively setting a lower limit on k. On the other hand, for large k values changes in ṁ happen on a much longer timescale, during which the system will have time to approach a quasi-steady state. For intermediate values of k, however, variations in the input will be faster than the time required to reach steady state, and the system should show its characteristic natural frequency, if it exists.

Our results for different choices of k are presented in Fig. 6. On the left-hand side we show the hard photon light curves (arbitrarily normalized to a peak value of 100), and on the right-hand side we display the corresponding average power spectra. For the averaging of the spectra, we split each light curve into nine segments. As can be seen, a broad feature appears between ∼3 and 4 Hz for k < 11. For k = 11 the feature starts diminishing and disappears altogether for larger values of k, resulting in featureless power-law spectra with the same slope as the input time series. Our findings imply that the natural frequency of the system is evident if small-amplitude variations in ṁ occur on a timescale Δτ ≃ (32 − 512)R_S/c = (3 − 50)(M_*/10 M_⊙) ms. These values, however, should be taken with some caution. If the changes in ṁ are assumed to be discontinuous between successive time steps (i.e., in the form of small-amplitude step functions), then the QPOs persist to much higher values of Δτ, yet another manifestation of the strong nonlinearity of the system.

Fig. 6. Hard photon light curves (panel a) and power spectra (panel b) obtained for red-noise mass accretion rates with a minimum to maximum ratio of 2. Other parameters used are: rc = 30, a = 0.2, and M* = 10 M⊙. The initial accretion rate is ṁ0 = 10−3.

4. Summary and discussion

In this work we have presented a simple model to describe the dynamic coupling of a hot electron population in the corona of accreting objects with the soft radiation emerging from the disk after the reprocessing of the hard corona radiation. Our approach is simple because we address only the temporal evolution of the system through a set of coupled differential equations, without treating any spectral features of the radiation. Nonetheless, it is self-consistent, as we couple the hard photon luminosity to the thermal electron losses in the corona. We find that the temporal properties of the corona-disk system are similar to those of the classical prey-predator dynamical system (Lotka 1920; Volterra 1926), with electrons acting as the prey and soft photons acting as the predator.

The main physical parameters of the problem are the mass of the central black hole M_*, the mass accretion rate (in units of the Eddington rate) ṁ, the radius of the corona r_c (in units of the Schwarzschild radius), and the ratio of hard-to-soft (reprocessed) radiation a. Even if all parameters are constant in time, the dynamic corona-disk system exhibits damped oscillations with a characteristic frequency $f_0\propto {\dot{m}}^{1/2}{r}_{\mathrm{c}}^{-3/2}{a}^{1/2}{M}_{\ast }^{-1}$. The frequency persists even in cases where ṁ varies stochastically and is independent of the initial conditions, strongly indicating that it is an inherent property of the system. The natural frequency of the corona-disk system ranges between ∼1 Hz and tens of hertz for typical parameter values for BHXRBs, with a tendency to increase with increasing ṁ. Because of the derived scaling of f₀ with mass, our model predicts quasi-periodicity on timescales of a few hours (10⁻⁵ − 10⁻⁴ Hz) to several weeks (10⁻⁸ − 10⁻⁷ Hz) for supermassive black holes in AGN with masses between 10⁶ and 10⁹ M_⊙ (for candidate QPOs in AGN, see, e.g., Middleton & Done 2010; Carpano & Jin 2018; Ashton & Middleton 2021). On the other end, it predicts kilohertz-scale QPOs for neutron stars with masses of ≃1 M_⊙, which indeed seems to be the case for the Rossi X-ray Timing Explorer observations of 4U 1728-34 and Sco X-1 made more than 25 years ago (Strohmayer et al. 1996; van der Klis et al. 1996).

Disk soft photons tend to stabilize the system by washing out the oscillations and bringing it to a steady state. We find that the nonlinear coupling of hot electrons and hard photons appears only when the amount of available accretion power goes overwhelmingly to heating up the corona and not the disk (see panel d in Fig. 3). Therefore, the recorded oscillations would appear only in the so-called hard and hard-intermediate states of the sources (for a review, see Ingram & Motta 2019), or, in other words, the system can switch from linear to nonlinear behavior depending on its state.

We note that, by construction, our model cannot calculate the (dimensionless) temperature Θ and the Thomson optical depth of the corona τ_T separately, but rather the product of the two quantities, namely τ_TΘ ≃ 0.25α⁻¹, which is of order unity. The inferred electron temperature generally falls within the range of the one derived from the spectral modeling of BHXRBs in the hard state for 1 ≲ τ_T ≲ 5 (see, e.g., Malzac 2012, and references therein). A value of τ_T > 1, however, can delay the escape of photons due to their spatial diffusion, as they scatter off the electrons, and thus our choice of the hard photon escape timescale, t_c = R_c/c (see Eq. (8)), can be considered as a lower limit. In such cases, one should usually introduce a photon diffusion timescale t_c, diff = R_c(1 + τ_T/3)/c that takes the above effect into account. Since in most corona models 1 ≲ τ_T ≲ 5, introduction of this correction would have modified f₀ by a factor of ∼2.

In our model we have only considered inverse Compton scattering as the process coupling photons and electrons. However, there are other competing processes that contribute to the spectral formation and electron cooling, such as bremsstrahlung and double Compton scattering. We have examined the effects of the former by adding an energy-loss term in Eq. (2) and a similar source term in Eq. (3) (for rates, see, e.g., Rybicki & Lightman 1986), assuming an electron-proton plasma. We have found no significant changes in our results; that is to say, bremsstrahlung, despite its linear nature, cannot stabilize the system. This is to be expected to an extent, since inverse Compton losses are typically much more important as an energy loss mechanism in BHXRBs than bremsstrahlung.

In our analysis, electrons were confined to the corona. However, inspection of Eq. (1) reveals that electron escape can only play a role in cases when the associated escape timescale becomes shorter than the electron cooling timescale. We typically find that the system quickly reaches a steady state without undergoing oscillations if electrons escape on timescales shorter than ∼(100 − 500)R_c/c. Therefore, a way of quenching type-C QPO activity is by making electrons escape quickly from the corona. This could occur as the source goes from the hard state to the hard-intermediate one, then to the soft-intermediate state, and eventually to the soft state, because the corona shrinks and therefore the outflow becomes continuously narrower. In the soft-intermediate state, the outflow is quite narrow and, for the first time, we see type-B QPOs, which are explained quantitatively as precession of the outflow (Kylafis et al. 2020). When the outflow is wide, its precession is not visible; thus, no type-B QPOs are seen in the hard state, only type-C QPOs, which, as we explain, arise from the interaction of the hot corona with the cold disk. Of course, another factor that quenches the type-C QPO activity in the soft state is the increased soft disk radiation, as explained in Sect. 3.

We should emphasize that, due to the damped oscillatory behavior of the disk-corona system that we have found, one needs some type of flickering in ṁ in order to sustain the QPOs. Interestingly enough, as we have shown in Sect. 3, even small-amplitude perturbations are adequate for that, provided that each perturbation occurs before the oscillations caused by the previous one die out. Therefore, perturbations around a constant value of ṁ could produce QPOs of a constant frequency for the whole time the system remains in the hard or hard-intermediate state. We note that such rapid variations (flickering) are evident in the power spectra of BHXRBs.

To explain observed changes in the QPO frequency on longer timescales (e.g., days to weeks), a longer variability timescale needs to be introduced in our model, and this cannot be calculated self-consistently, because it depends on external causes. For instance, if the average accretion rate slowly increases on a timescale of weeks, the QPO frequency will also increase. However, the model still requires the aforementioned flickering in ṁ to sustain the QPOs. We note that since the required flickering is many orders of magnitude faster than the timescale for the hypothesized secular change in ṁ, it can safely be assumed that this occurs around the (slowly) changing value of ṁ. Clearly, no jumps are expected in the QPO frequency evolution unless the flickering in the mass accretion rate changes abruptly by a large factor (see Eq. (A.11)).

Our findings motivate detailed spectro-temporal calculations of the disk-corona system with sophisticated codes that follow the spectral evolution of the interacting particle populations in each time step (e.g., Belmont et al. 2008). Such calculations could also possibly self-consistently address issues such as the possible dependence of r_c on ṁ and the soft-hard lags. Numerical investigation of feedback processes in nonthermal leptonic and hadronic plasmas (Petropoulou & Mastichiadis 2012; Mastichiadis et al. 2020) has confirmed the presence of the QPOs found in simplified models. Numerical confirmation of the natural frequency computed here would make the corona-disk radiative feedback a novel scenario for QPOs in BHXRBs.

Acknowledgments

We would like to thank the anonymous referee for useful comments that helped us improve the manuscript and Dr. Georgios Vasilopoulos for brainstorming discussions. M.P. acknowledges support from the MERAC Foundation through the project THRILL. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), and Stingray, a Python library for spectral timing (Huppenkothen et al. 2019a,b).

References

Appendix A: Analytical derivation of the natural frequency

Even though finding explicit formal solutions to Eqs. (7) and (8) is not trivial (e.g., for the Lotka-Volterra system, see Evans & Findley 1999), one can gain insight into the dynamical properties of the system by performing a stability analysis of its fixed (stationary) points (see also Petropoulou & Mastichiadis 2012).

Equations (7) and (8) can be recast into the following form,

$$\begin{array}{cc}\hfill \dot{x}& =\alpha -\beta x(\gamma y+\delta )\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill \dot{y}& =-yϵ+\beta x(\gamma y+\delta )\hfill \end{array}$$(A.2)

where $\dot{x}\equiv \mathrm{d}x/\mathrm{d}\stackrel{\sim }{t}$, $\alpha =3f\dot{m}/(2\xi r_{\mathrm{c}}^3)$, β = 4m_p/(3m_e), γ = a, $\delta =3(1-f)\dot{m}/(2\xi r_{\mathrm{c}}^2)$, and ϵ = 1/r_c. These equations have one equilibrium point (x₀, y₀) – commonly referred to as the steady state – which is given by

$$\begin{array}{cc}\hfill x_0& =\frac{\alpha }{\beta }\frac{1}{\gamma \frac{\alpha }{ϵ}+\delta }\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill y_0& =\frac{\alpha }{ϵ}.\hfill \end{array}$$(A.4)

We examined the behavior of the system when it is slightly perturbed from its steady state (i.e., x = x₀ + x′, y = y₀ + y′). We linearized Eqs. (A.1)-(A.2) with respect to the perturbed quantities x′,y′ and find

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\dot{x}}^{\prime }\\ {\dot{y}}^{\prime }\end{array}\right)=\mathbf{J}(x_0,y_0)\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right),\end{array}$$(A.5)

where

$$\begin{array}{c}\hfill \mathbf{J}(x_0,y_0)=\left(\begin{array}{cc}-\beta (\gamma y_0+\delta )& -\beta \gamma x_0\\ +\beta (\gamma y_0+\delta )& -ϵ+\beta \gamma x_0\end{array}\right)\end{array}$$(A.6)

is the Jacobian matrix of the linearized system of equations at the stationary point.

The evolution of the perturbed quantities near its steady state can be described by the eigenvalues of J(x₀,y₀), namely $x^{\prime },y^{\prime }\propto e^{\lambda \stackrel{\sim }{t}}$ (e.g., Strogatz 2000).

To simplify the calculations, we dropped the constant term δ from the Jacobian matrix since δ ≪ γy₀ (or equivalently, 1 − f ≪ af). The eigenvalues of J(x₀,y₀) are the roots of its characteristic polynomial,

$$\begin{array}{c}\hfill P(\lambda )={\lambda }^2+\frac{\alpha \beta \gamma }{ϵ}\lambda +\alpha \beta \gamma .\end{array}$$(A.7)

The system has two complex conjugate eigenvalues,

$$\begin{array}{c}\hfill {\lambda }_{1,2}=\frac{-\frac{\alpha \beta \gamma }{ϵ}\pm i\sqrt{|D|}}{2},\end{array}$$(A.8)

if the discriminant D of P(λ) = 0 is negative, which translates to the following constraint

$$\begin{array}{c}\hfill \frac{\alpha \beta \gamma }{{ϵ}^2}<4\Rightarrow \left(\frac{a}{0.1}\right)\left(\frac{\dot{m}}{{10}^{-3}}\right){\left(\frac{\xi }{1}\right)}^{-1}{\left(\frac{r_{\mathrm{c}}}{10}\right)}^{-1}<{10}^2.\end{array}$$(A.9)

The real part of the eigenvalues is always negative, which means that the steady state is a stable point. In other words, if the system is perturbed from its steady state (e.g., due to a small change in the accretion rate), it will return to it after exhibiting damped oscillations at a frequency given by the imaginary part of the eigenvalues,

$$\begin{array}{c}\hfill \stackrel{\sim }{\omega }=\sqrt{\alpha \beta \gamma }\sqrt{1-\frac{\alpha \beta \gamma }{4{ϵ}^2}}\approx {\left(\frac{2m_{\mathrm{p}}}{m_{\mathrm{e}}}\frac{\mathit{af}}{\xi }\frac{\dot{m}}{r_{\mathrm{c}}^3}\right)}^{1/2}.\end{array}$$(A.10)

The right-hand side of the equation was obtained using the approximation αβγ ≪ 4ϵ², which is satisfied for most parameter values (see also Eq. A.9). The natural frequency of the system (in Hz) can then be written as $f_0=\stackrel{\sim }{\omega }c/(2\pi R_{\mathrm{S}})$

$$\begin{array}{c}\hfill f_0\simeq 30\mathrm{Hz}{\left(\frac{a}{0.1}\right)}^{1/2}{\left(\frac{\dot{m}}{{10}^{-3}}\right)}^{1/2}{\left(\frac{r_{\mathrm{c}}}{10}\right)}^{-3/2}{\left(\frac{\xi f}{1}\right)}^{-1/2}{\left(\frac{M_{\ast }}{10M_{\odot }}\right)}^{-1}.\end{array}$$(A.11)

For parameter values that do not satisfy condition (A.9), the eigenvalues are both real and negative. In this case, the system will return to its steady state (i.e., a stable fixed point) without showing oscillatory behavior.

All Figures

	Fig. 1. Temporal evolution of energy densities for different model parameters: ṁ = 5 × 10−4, rc = 30, a = 0.2 (panel a), ṁ = 5 × 10−4, rc = 15, a = 0.2 (panel b), ṁ = 5 × 10−3, rc = 30, a = 0.2 (panel c), and ṁ = 5 × 10−4, rc = 30, a = 0.05 (panel d). In all cases, f = 0.999, ξ = 3, rd = 10, and time is normalized to RS/c.
In the text

	Fig. 2. Phase-space plot illustrating the dynamic behavior of the corona-disk system. Here, the same parameters as those in panels a, b, and d of Fig. 1 are used. For a better display, the energy densities are normalized to a peak value of one.
In the text

	Fig. 3. Temporal evolution of the dimensionless hard photon energy density for different values of ṁ (panel a), rc (panel b), a (panel c), and f (panel d). Unless stated otherwise, the parameters used are: ṁ = 5 × 10−3, rc = 30, a = 0.05, and f = 0.999. Time is normalized to RS/c.
In the text

Fig. 4. Evolution of the hard photon density in time and frequency domains. Panel a: light curve of hard (Comptonized) photons (rescaled to a maximum value of 100) for the same parameters used in panel (a) of Fig. 1. Panel b: power spectrum of the light curve multiplied by frequency using the Leahy normalization. Overplotted in orange is the averaged power spectrum computed using five segments of the light curve.

In the text

Fig. 5. Natural frequency of the dynamic corona-disk system for a nonrotating black hole with mass M* = 10 M⊙ as a function of the mass accretion rate (in units of the Eddington rate) for different values of rc (panel a) and a (panel b), as indicated in the color bar. Beyond (or below) a certain mass accretion rate, the system reaches a steady state without exhibiting oscillations (dotted black lines). Other parameters used are: a = 0.2 (panel a) and rc = 30 (panel b), with f = 0.999 in both panels.

In the text

	Fig. 6. Hard photon light curves (panel a) and power spectra (panel b) obtained for red-noise mass accretion rates with a minimum to maximum ratio of 2. Other parameters used are: rc = 30, a = 0.2, and M* = 10 M⊙. The initial accretion rate is ṁ0 = 10−3.
In the text