© The Authors 2026

1. Introduction

A black-hole X-ray binary (BHXRB) consists of a black hole accreting matter from a low-mass companion star. The X-ray spectra of these objects in their hard state (for a review of the BHXRB spectral states see Belloni 2010) consist of two components, a blackbody at low (< 1 keV) energies and a power law that extends up to hundreds of keV energies and, in some cases, up to MeV. While the former is attributed to accretion disk emission, the latter is thought to be produced by inverse Compton scattering (ICS) of energetic electrons on the soft disk component. The electrons themselves are assumed to be energized in a region close to the central black hole, known as the corona; for a review, see Done et al. (2007).

Timing analysis, especially during outbursts, reveals the presence of QPOs (van der Klis 1989; Psaltis et al. 1999; Giannios & Spruit 2004; Ingram & Motta 2019). Among these, the so-called type-C QPOs are characterized by a strong peak in their power spectra and appear when the source is in the hard, or hard-intermediate, state. The origin of this type of QPOs is still debated and models include, among others, instabilities in the accretion flow (Tagger & Pellat 1999), oscillations of boundary layers (Titarchuk & Fiorito 2004), or Lense-Thirring precession (Stella & Vietri 1998; Ingram et al. 2009).

In a recent work, Mastichiadis et al. (2022) (henceforth MPK) showed that interactions between the corona electrons and the accretion disk soft photons can be inherently non-linear, giving rise to an oscillatory pattern in the X-ray flux, reminiscent of limit cycles found in non-linear dynamical systems. Their model was a dynamical treatment of Haardt & Maraschi (1991): the soft photons of the accretion disk cool by ICS the corona electrons → part of the electron ICS hard radiation is reprocessed on the accretion disk, producing extra soft photons there → the extra soft photons further cool the electrons. This succession of steps forms a positive feedback loop between the electrons and the soft photons as illustrated in Fig. 1. If the energetic electrons are continuously replenished, the modeling of the electron-photon coupling leads to a type of Lotka-Volterra, or predator-prey, system of equations with electrons being the prey and the soft photons the predators. The solution of this type of equation leads to limit cycles. MPK’s treatment was simplified in the sense that they used energy integrated densities for the two species. Therefore, their results were limited to the temporal properties of the system; yet, their approach was self-consistent and kept the physical essence of the interactions.

Fig. 1. Schematic representation of the feedback loop between electrons and photons in the coronal environment. For artistic reasons the corona is depicted as a cloud, while in the paper its shape is assumed to be spherical. (The figure is AI generated).

In the present paper, we extend the above approach by considering the full spectro-temporal evolution of the system. This essentially means that the set of two coupled differential equations that described the system in the case of MPK has to be replaced by a large set of coupled partially differential equations that describe the simultaneous evolution in time and energy of the electron and photon distribution functions. In order to solve these equations we have used a numerical code that self-consistently treats the relevant processes between electrons and photons (Mastichiadis & Kirk 1995; Dimitrakoudis et al. 2012). Electrons are assumed to be energized inside the corona by some (unspecified) mechanism and, at the same time, lose energy on soft photons. These are not independent of the electrons but come from reprocessing of a part of the electron ICS radiation on the accretion disk. This coupling of electron energization and cooling can give QPOs, as a result of the limit cycle behavior, to the system, exactly as in MPK. However, in this case one can get more information because the numerical code gives the electron and photon spectrum at every single time instant.

The aims of the present paper are (i) to show that limit cycles are inherent in a corona-disk system and (ii) to make a broad study of the physical parameters that can give results resembling the spectral and temporal behavior of the BHXRBs in their hard X-ray state. It 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 briefly present the specifics of the code used. In Sect. 4 we give a typical example of the numerical solutions to the equations of the corona-disk system and define some temporal and spectral parameters that characterize these solutions. In Sect. 5 we systematize our results by investigating the effect of each of the main parameters of the problem. In Sect. 6 we construct and solve analytically a simplified set of equations showing explicitly why the non-linearity arises. In Sect. 7 we give an application on the system GRS 1915+105 and we conclude in Sect. 8 with a summary and a discussion.

2. Model description

We start by adopting the standard picture of an accretion disk - corona system. The exact geometry is not of relevance here; all we require is that the disk provides soft photons for the cooling of the electrons of the corona. Electrons enter the corona, which is a spherical region of radius R_c, where they are energized by some unspecified mechanism. For the purpose of the present work it is sufficient to attribute to this mechanism only a characteristic energization timescale t_en that is assumed to be independent of the electron energy for simplicity. Here we will adopt the approach of Kirk et al. (1998) – see also Petropoulou et al. (2024) for a recent application to non-thermal flares from SgrA^*.

In the absence of losses, the electron Lorentz factor obeys the characteristic equation

$$\begin{array}{c}\hfill \frac{d\gamma }{\mathit{dt}}=\frac{\gamma }{t_{\mathrm{en}}},\end{array}$$(1)

which, for initial condition γ(t₀) = γ₀, gives

$$\begin{array}{c}\hfill \gamma (t)={\gamma }_0{e}^{(t-t_0)/t_{\mathrm{en}}}.\end{array}$$(2)

As the electrons are energized, they start losing energy through various processes. If we concentrate on ICS, then we can write the characteristic equation for the electron energy as

$$\begin{array}{c}\hfill \frac{d\gamma }{\mathit{dt}}=\frac{\gamma }{t_{\mathrm{en}}}-\frac{4}{3}{\sigma }_{T}c\left(\frac{U_{\mathrm{s}}}{m_e{c}^2}\right){\gamma }^2,\end{array}$$(3)

where σ_T is the Thomson cross-section, m_e is the electron rest mass, and U_s is the soft photon energy density. Note that while for U_s= constant, the above equation leads to an equilibrium at

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{sat}}=\frac{3m_{e}c}{4{\sigma }_T{U}_s{t}_{\mathrm{en}}},\end{array}$$(4)

a steady state might not be established in the case U_s = U_s(t).

To describe, therefore, the behavior of electrons in time and energy inside the corona we assume that these enter homogeneously the region with initial Lorentz factors γ_inj at a rate q_inj per unit volume. Once inside the corona the electrons can be energized with the timescale t_en, lose their energy by ICS, or physically escape the source on a timescale t_esc. Thus, in the case where ICS is the sole energy loss mechanism, the temporal evolution of the differential density electron distribution, n_e(γ, t) (units volume⁻¹ energy⁻¹), is described by a partial differential equation (Kirk et al. 1998)

$$\begin{array}{c}\hfill \frac{\partial n_e(\gamma ,t)}{\partial t}+\frac{\partial }{\partial \gamma }\left[(\frac{\gamma }{t_{\mathrm{en}}}-b_c{\gamma }^2)n_e(\gamma ,t)\right]+\frac{n_e(\gamma ,t)}{t_{\mathrm{esc}}}=q_{\mathrm{inj}}\delta (\gamma -{\gamma }_{\mathrm{inj}}),\end{array}$$(5)

where b_c = 4σ_TU_s/(3m_ec) and δ(x) is the Dirac delta function. In the case where t_en and t_esc are independent of the electron energy, q_inj does not vary in time and γ is far from the upper cutoff of the distribution, the solution is described by a power law, i.e. n_e(γ)∝γ^{−1 − (t_en/t_esc)} which is asymptotically compared to −1 as t_esc → ∞. In such cases, most of the total electron energy content is carried by the high energy end of the electron distribution, and this will be one of our assumptions in this paper.

Eq. (2) implies that the electron cooling term produces radiation and should be coupled by an analogous equation for photons that describes their spectral and temporal properties inside the corona. Assuming again that ICS is the main mechanism for electron radiation, the equation for photons can be written as follows

$$\begin{array}{c}\hfill \frac{\partial n_{\gamma }(ϵ,t)}{\partial t}+\frac{n_{\gamma }(ϵ,t)}{t_{\mathrm{cr}}}={\mathcal{Q}}_{\mathit{ICS}}(ϵ,t),\end{array}$$(6)

where n_γ(ϵ, t) is the differential photon density (units volume⁻¹ energy⁻¹), ϵ is the photon energy, and t_cr = R_c/c is the crossing time of the corona. The term n_γ/t_cr denotes the photon escape while 𝒬_ICS is the IC emissivity; see, for example, Blumenthal & Gould (1970).

Eqs. (5) and (6) are coupled because 𝒬_ICS, which appears in the photon equation, is a function of n_e(γ, t). Furthermore, the solution of the electron equation depends on the photon energy density U_s, which for the corona-disk case might depend on n_e, thus making the problem nonlinear. Therefore, in order to calculate U_s, we have to make some further assumptions about its origin. Here we follow the approach of MPK and assume that there are two components contributing to U_s: (i) the disk thermal radiation of temperature T_d, which is assumed to be constant and (ii) the reprocessed hard radiation from the corona impinging on the disk, which depends on the ICS component calculated from Eq. (6). Therefore, we can write

$$\begin{array}{c}\hfill U_s=U_{\mathrm{ds}}+U_{\mathrm{rs}},\end{array}$$(7)

where U_ds is the energy density of the thermal radiation of temperatude T_d coming from the disk and U_rs is the energy density of the reprocessed radiation. For the latter, we assume that it also emits as a gray body of temperature T_eff and, to avoid the use of many parameters, we set T_eff = T_d. The energy density of the reprocessed component is calculated at each instant by the relation

$$\begin{array}{c}\hfill U_{\mathrm{rs}}(t)=\alpha {\displaystyle \int d}ϵϵn_{\gamma }(ϵ,t),\end{array}$$(8)

with the integration taking into account only the hard photons of energies > k_BT_eff while α is the fraction of the hard corona radiation that is reprocessed by the disk.

The above describe both the parameters and the assumptions of the problem and, therefore, one can solve the electron and photon equations in a straightforward manner to obtain their tempo-spectral behavior. In the next section, we give some details about the numerical code used and also relate the parameters to the physics of BHXRBs.

3. Code specifics

The numerical code used is based on the leptonic part of the time-dependent code ATHEνA (Mastichiadis & Kirk 1995; Dimitrakoudis et al. 2012), which takes into account mostly nonthermal processes. For an updated description and comparisons with other similar codes, see Cerruti et al. (2026). The code has been modified in two ways: (a) it uses an explicit electron energization scheme as used in Kirk et al. (1998) and Petropoulou et al. (2024), rather than an electron injection term in high energies as is the norm with this type of code, and (b) as explained in the previous section, it uses reprocessed photons for ICS cooling. This means that at each time step the ICS emissivity is integrated, turned into radiation density, and a fraction α of it is assumed to be radiated as a gray body of temperature T_eff. The energy snapshots have been recorded in time steps Δt = 0.1t_cr. The physical processes used are (i) inverse Compton scattering, (ii) photon-photon pair production, and (iii) electron-proton bremsstrahlung. It turns out that for the majority of the parameters used, the impact of the two latter processes is minimal; however, they have been used throughout for completeness reasons. In contrast, synchrotron radiation that can potentially have an effect on our results is not taken into account; the implications that this has on the value of the B-field will be discussed in Sect. 6 and Appendix C.

The free parameters of the problem are the following.

• The corona radius R_c.

• The electron energization timescale t_en.

• The electron injection rate q_inj or, equivalently, the mass accretion rate Ṁ_inj (see below).

• The electron escape timescale t_esc.

• The accretion disk thermal luminosity.

• The temperature of the disk and the reprocessed component T_eff, assumed to be equal.

• The fraction α of the hard radiation that is reprocessed by the disk and is radiated back as soft photons.

We note that most of the above parameters have estimates in the BHXRB literature. Only t_en and t_esc are practically unknown, therefore we use generic values for them that are simple multiples of the crossing timescale t_cr = R_c/c and we use ${\stackrel{\sim }{t}}_{\mathrm{en}}=t_{\mathrm{en}}/t_{\mathrm{cr}}$. For simplicity we set ${\stackrel{\sim }{t}}_{\mathrm{esc}}={10}^3{\stackrel{\sim }{t}}_{\mathrm{en}}$, i.e., we require that the electron distribution function be flat, carrying most of the electron energy close to its upper cutoff. Other values of the ratio t_esc/t_en are discussed at the end of Sect. 5.

For the electron injection term at low energies (see Eq. (5)) we use γ_inj = 10^0.05. Electrons with γ < γ_inj are considered cold. Note that the injection rate of electrons per unit volume q_inj, used in the same equation, does not have a straightforward connection to some observable quantity as it gives the rate of low energy electrons entering the energization process. However, if one assumes an electron-proton plasma, they can write the total injection rate

$$\begin{array}{c}\hfill Q_{e,\mathrm{inj}}\simeq \frac{4\pi }{3}{R}_c^3{q}_{\mathrm{inj}}{m}_e{c}^2=\xi \frac{{\dot{M}}_{\mathrm{inj}}}{m_p},\end{array}$$(9)

where Ṁ_inj is the mass accretion rate, m_p is the proton mass, and ξ is the fraction of electrons that enter the energization process and which, for simplicity, we assume to be equal to 1. Note that Eq. (9) does not imply a relation of Ṁ_inj to the luminosity of the system. We return to this point in the next section.

Finally, since the system of Eqs. (5) and (6) is time-dependent, we use as initial conditions n_e(γ, 0) = n_γ(ϵ, 0) = 0 everywhere with the exception of the electrons first bin at γ_inj = 10^0.05, which takes the boundary value n_e(γ_inj, 0) = q_injt_en; see Kirk et al. (1998).

In the next section, we give a typical example of how the system moves from linear to nonlinear behavior and use this to define some useful quantities that we will use extensively for the rest of the paper. However, a detailed description of the mechanisms behind the limit cycle behavior can be found in Appendix A.

4. The effect of disk photons on the temporal evolution of the system

We show below some results that represent the typical features of our model. In the next section we use normalized quantities. So, the hard luminosity, the soft energy density of the disk, and the energy density reprocessed on the accretion disk are given in terms of compactness

$$\begin{array}{cc}\hfill {\ell }_h& =\frac{L_h{\sigma }_t}{4\pi R_c{m}_e{c}^3},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\ell }_{\mathrm{ds}}& =\frac{U_{\mathrm{ds}}{\sigma }_T{R}_c}{m_e{c}^2},\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {\ell }_{\mathrm{rs}}& =\frac{U_{\mathrm{rs}}{\sigma }_T{R}_c}{m_e{c}^2},\hfill \end{array}$$(12)

respectively, where L_h is the hard luminosity produced from ICS. Furthermore, time is measured in light crossing time of the corona. The photon energy is given in keV.

As a first result, we show the effects of soft disk radiation, expressed through ℓ_ds, on the temporal behavior of ℓ_h. Figure 2 shows a plot of ℓ_h versus normalized time for various values of ℓ_ds with all other parameters fixed. As ℓ_ds decreases, the system becomes nonlinear and limit cycles start appearing; see the succession of light curves from violet to red. Generally speaking, as long as ℓ_ds ≪ ℓ_rs, the cooling occurs in the non-linear regime. From the light curves obtained various quantities can be extracted, such as: (i) The frequency of QPO ν_QPO. This is calculated directly from the light curves using standard Python routines. Here and for the rest of the paper, for computational resource reasons, we have used the first ten cycles of the light curves, excluding the first one in order to minimize the effects of the initial conditions. Despite the fact that ν_QPO can be calculated for the whole hard spectrum, here we have chosen the energy range 2–200 keV to maintain an analogy to the observations that are usually performed in this energy range. (ii) The fractional rms of QPO f_rms. This was calculated for the same section of the light curves and for the same energy range as above. (iii) The timescale ${\stackrel{\sim }{T}}_{10}$, expressed in units of t_cr, required for the X-ray light curve to drop to 10% of its value at the second peak. As the damping of the obtained light curves does not fit any simple analytical law, we resorted to the calculation of ${\stackrel{\sim }{T}}_{10}$. For this we used the expression for the difference norm

$$\begin{array}{c}\hfill D(\stackrel{\sim }{t})=|n_{\gamma }({ϵ}_i,\stackrel{\sim }{t}+\mathrm{\Delta }\stackrel{\sim }{t})-n_{\gamma }({ϵ}_i,\stackrel{\sim }{t})|,\end{array}$$(13)

Fig. 2. Light curves of ICS compactness in the 2–200 keV band for various values of the parameter ℓds. All other parameters have been kept constant and have the values $R_c={10}^9\mathrm{cm},{\dot{M}}_{\mathrm{inj}}=3.2\times {10}^{-7}{M}_{\odot }/\mathrm{yr},{\stackrel{\sim }{t}}_{\mathrm{en}}=2.5,t_{\mathrm{esc}}/t_{\mathrm{en}}={10}^3,T_{\mathrm{eff}}=4\times {10}^6\mathrm{K},\alpha =0.1$. The colors of the lines correspond to ℓds = 10−6,  10−5,  10−4,  10−3,  10−2, and 10−1 from red to violet. For the sake of clarity each curve is displaced from the previous one by unity.

and numerically determined ${\stackrel{\sim }{T}}_{10}$ from the equation

$$\begin{array}{c}\hfill D({\stackrel{\sim }{T}}_{10})=0.1D\left({\stackrel{\sim }{t}}_2^{\mathrm{pk}}\right),\end{array}$$(14)

where $D\left({\stackrel{\sim }{t}}_2^{\mathrm{pk}}\right)$ is the value of the norm at the second peak of the light curve. Here $\mathrm{\Delta }\stackrel{\sim }{t}=0.1$ and by $\stackrel{\sim }{t}$ we mean that time is normalized to t_cr. (iv) The phase lag Δϕ between two characteristic energy bands. Here we chose the 2–6 keV as the “soft” energy band and the 6–15 keV as the “hard” energy band in line with most of the QPO literature, for example, Zhang et al. (2020). We calculated the time lag ΔT_sh between the two bands and, consequently, calculated the corresponding phase lag by using the relation Δϕ = 2πν_QPOΔT_sh. In this work, a positive lag means that the hard photons lag the soft ones. A detailed description of the behavior of ΔT_sh can be found in Appendix B.

The thus-calculated quantities are shown in Table 1 as a function of the compactness of the disk ℓ_ds; here, in addition to Δϕ, we have tabulated ΔT_sh (in units of t_cr) as well, as this is the quantity that is directly computed from the numerical code. Several conclusions can be drawn that, despite covering only one set of parameters, have more general consequences: (a) As already mentioned above, as ℓ_ds is decreasing, the system becomes nonlinear and limit cycles start appearing. Generally speaking, as long as ℓ_ds ≪ ℓ_rs, the cooling occurs in the non-linear regime. We also note that once established, the period of the oscillations remains largely unaffected by the value of ℓ_ds, i.e., there is a limit that the system tends to reach, which is its natural frequency. The lower the value of ℓ_ds, the stronger the limit cycles. This result verifies the findings of MPK that QPOs appear when the disk is not radiating, but instead acts as a reprocessor of hard ICS radiation. On the other hand, a luminous disk damps very effectively the limit cycles as, in this case, the plentiful disk photons provide efficient cooling for the corona electrons linearizing, at the same time, the whole process. (b) In contrast to MPK, who found only damped oscillations, here we find cases that are, practically speaking, undamped. We explain this difference in the next section. (c) The quantity ${\stackrel{\sim }{T}}_{10}$ is a good alias for the damping of the system. Generally speaking the runs that show ${\stackrel{\sim }{T}}_{10}<50$ are strongly damped. On the other hand, runs with ${\stackrel{\sim }{T}}_{10}>500$ can practically be considered as undamped. (d) In contrast to ν_QPO, which is rather insensitive to damping, fractional rms is a sensitive function of it. Therefore, the values shown in the table for ${\stackrel{\sim }{T}}_{10}<100$, where damping starts to become stronger, can only be considered as upper limits. (e) The quantity ΔT_sh that measures soft/hard lags soft photons lag the hard ones.

Figure 3 shows the averaged, over the first ten cycles, photon spectra for the same cases as above. In these plots, one can see the superposition of the disk and reprocessed soft photons, which have been assumed to be gray bodies of temperature T_eff = 4 × 10⁶ K, and the inverse Compton component that reaches, in this particular example, MeV energies. Here, the y−axis is in units $x\frac{d\ell }{\mathit{dx}}=x^2{\dot{n}}_{\gamma }(x,t)$, where x = ϵ/m_ec². This is equivalent to $ϵ\frac{d{L}_{\gamma }}{dϵ}$ where L_γ is the spectral luminosity.

Fig. 3. Average photon spectra for the same cases shown in the previous figure. The color codes have been kept the same.

As one can observe, low values of ℓ_ds produce harder ICS spectra; however, as ℓ_ds increases, the disk photons hinder electron acceleration, and as a result both electrons are lower and the ICS photons reach lower energies (see blue line curve). Ultimately, in the case of strong disk emission, electrons cannot, practically speaking, be energized, and, therefore, the photon spectrum consists only of the input disk gray-body photon distribution (violet line curve).

From the spectra shown in Fig. 3 various quantities can be extracted. These are: (i) The bolometric inverse Compton luminosity, L_bol. (ii) The X-ray luminosity between 2 and 200 keV, L_X. (iii) The peak of the ICS component, E_pk, in Fig. 3 is in units such that it immediately shows this quantity. (iv) The slope of the photon spectrum between 2–200 keV. (v) The efficiency of the acceleration process η = Ṁ_eff/Ṁ_inj, where Ṁ_eff = L_bol/c². This quantity signifies the fraction of the injected rest mass that ultimately goes into radiation. (vi) The optical depth of the cooled electrons τ_T. This is the optical depth of the electrons that have cooled by ICS at Lorentz factors γ < γ_inj including the unenergized ones of the flow.

Table 2 shows the above quantities, with the exception of τ_T, which remains in most cases below 0.1 signifying an optically thin medium¹. All remain practically constant for low values of ℓ_ds; however, as ℓ_ds and the degree of damping increases, the luminosity drops and the spectrum steepens.

5. General trends

According to the findings of the previous section, by assuming a cold disk, we are maximizing the chance of a limit cycle appearance. Therefore, by fixing the disk photon energy density at some low value, i.e., ℓ_ds = 10⁻⁶, we search for the appearance of QPO by changing three key parameters. For this we first obtain light curves for various values of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ while keeping all other parameters fixed. Once we detect the limit cycle behavior, we calculate the quantities defined in the previous section and given in Tables 1 and 2. As a next step, we repeat the same procedure, this time changing Ṁ_inj, and, finally, we consider variations in R_c. At the end of the section we discuss briefly variations in auxiliary parameters such as the reprocessing factor α, the reprocessing temperature T_eff, and the ratio t_esc/t_en.

5.1. Energization timescale ${\stackrel{\sim }{t}}_{\mathrm{en}}$ variations

As we have not adopted any specific theory behind the electron energization process in the corona, we follow a generic approach and vary the energization timescale ${\stackrel{\sim }{t}}_{\mathrm{en}}$ from fast (${\stackrel{\sim }{t}}_{\mathrm{en}}=1.25$) to slow (${\stackrel{\sim }{t}}_{\mathrm{en}}=80$) using a multiplication step of 2 between successive runs. We note that in the context of the one-zone model we require ${\stackrel{\sim }{t}}_{\mathrm{en}}>1$.

Figures 4 and 5 show, respectively, the light curves and averaged photon spectra produced for the values described above of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ when Ṁ_inj = 16 × 10⁻⁸ M_⊙ yr and R_c = 10⁹ cm, while Table 3 tabulates the parameters defined in the previous section. The parameters for these runs are given in the captions of Fig. 4 and Table 3.

Fig. 4. Set of light curves all sharing the same value of Qe, inj (Minj = 16 × 10−8 M⊙/yr) but having different energization timescales. Here ${\stackrel{\sim }{t}}_{\mathrm{en}}=1.25$ (violet), 2.5 (dark blue), 5 (light blue), 10 (green), 20 (yellow), 40 (orange), and 80 (red), all given in units of tcr. The other parameters are Rc = 109 cm,  tesc/ten = 103,  α = 0.1,  Teff = 4 × 106 K,  ℓds = 10−6. For the sake of clarity each curve is displaced from the previous one by 1.

Fig. 5. Average multiwavelength spectrum for the case shown in the previous figure, keeping the color codes the same.

There are several conclusions that can be deduced from these runs: (i) QPOs appear for all values of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ with strong damping that becomes evident only for large values of the energization timescale, in the present example for ${\stackrel{\sim }{t}}_{\mathrm{en}}>40$. Therefore, according to the discussion in Sect. 4, these cases are inherently damped, i.e., their damping does not depend on the level of disk radiation. However, as we will show in the next subsection, strong damping can occur at different values of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ for other values of Ṁ_inj. (ii) The QPO frequency ν_QPO decreases with increasing ${\stackrel{\sim }{t}}_{\mathrm{en}}$ up to ${\stackrel{\sim }{t}}_{\mathrm{en}}=20$, then it starts to increase again. However, as mentioned above, the runs with ${\stackrel{\sim }{t}}_{\mathrm{en}}>40$ produce heavily damped oscillations, so these cases can be considered outliers. (iii) Faster acceleration timescales tend to produce spectra, which have higher bolometric luminosities, peak at higher photon energies, and have harder spectral indices Γ. Note that L_bol, and therefore η, increases roughly linearly with ${\stackrel{\sim }{t}}_{\mathrm{en}}^{-1}$. On the other hand, L_X changes much less than L_bol because as the photon spectra are hard, most of the luminosity is carried around the peak of the photon distribution, E_pk, which lies above the band 2–200 keV used to calculate L_x. However, for slow ${\stackrel{\sim }{t}}_{\mathrm{en}}$, E_pk is decreasing and L_X eventually becomes equal to L_bol. (iv) The fractional rms f_rms is on the order of tens of percent. However, these values must be considered strictly as upper limits, especially in those cases that are characterized by strong to intermediate damping. Generally, f_rms and ${\stackrel{\sim }{T}}_{10}$ have an erratic behavior that will be discussed later. (v) The quantity ΔT_sh used to measure soft and/or hard lags decreases with increasing ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and ultimately changes sign; therefore, the present model can exhibit both soft and hard lags depending on the initial parameters. An explanation of this is given in Appendix B. (vi) Both Fig. 5 and Table 3 show that the X-ray spectra are hard (Γ < 2). However, this does not apply for all cases. For example, cases with higher values of the reprocessing parameter α tend to produce steeper spectra.

5.2. Mass accretion rate ${\dot{M}}_{\mathrm{inj}}$ variations

We can repeat the procedure described above for different values of Ṁ_inj, while keeping the other parameters constant. Indicative results are shown in Table 4. We note that, contrary to the previous case, ν_QPO remains largely independent of the value of Ṁ_inj, while both L_bol and L_X tend to increase with Ṁ_inj. However, in the case of L_bol the increase is very slow, i.e., despite the fact that Ṁ_inj changes by more than two orders of magnitude (a factor of 2⁷ = 128 to be exact), L_bol changes only by about a factor of 3. The result is that the efficiency η is a decreasing function of Ṁ_inj. In the case of L_X the change is greater, by a factor of about 30, because E_pk increases as Ṁ_inj decreases.

Our numerical study reveals that the only non-monotonous quantities are the degree of damping, which we express in terms of the quantity ${\stackrel{\sim }{T}}_{10}$, and the fractional rms f_rms, which, however, is related to the former, in the sense that heavily damped oscillations have small values of both ${\stackrel{\sim }{T}}_{10}$ and f_rms and vice versa. Figure 6 shows a heat map for the ${\stackrel{\sim }{T}}_{10}$ parameter as a function of both ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁ_inj. Dark regions indicate light curves showing strong damping. We can deduce that only certain combinations of (i) slow ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and high Ṁ_inj (the dark patch on the upper right side) and (ii) fast ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and intermediate Ṁ_inj produce damped oscillations. In addition, all runs with Ṁ_inj > 3.2 × 10⁻⁷ M_⊙/yr and ${\stackrel{\sim }{t}}_{\mathrm{en}}<1.25$ produce heavily damped oscillations. At any rate, judging from Fig. 6 we deduce that the transitions from heavily damped to lightly damped limit cycles do not follow some easily explainable pattern.

Fig. 6. Heat map of the parameter ${\stackrel{\sim }{T}}_{10}$ as a function of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁinj for Rc = 109 cm. Other parameters are α = 0.1, ℓds = 10−6 and Teff = 4 × 106 K. Dark regions correspond to light curves with strong decay.

Figure 7 gives a plot of ν_QPO as a function of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ for various values of Ṁ_inj. The dashed lines connect the points that are characterized by strong damping (${\stackrel{\sim }{T}}_{10}<50)$. We note the uniformity of ν_QPO if we exclude the dashed lines. Even if Ṁ_inj varies by almost two orders of magnitude, ν_QPO varies in most cases by less than a factor of 1.5. Therefore, the QPO frequency depends very weakly on ${\dot{M}}_{\mathrm{inj}}$; note, however, that ${\dot{M}}_{\mathrm{inj}}$, as defined in the present paper, is not related directly to the X-ray luminosity, but instead to the electron injection at low energies; see Eq. (9). It is also interesting to note that roughly ${\nu }_{\mathrm{QPO}}\propto {\stackrel{\sim }{t}}_{\mathrm{en}}^{-1/2}$ and compare the slope of the black line to that of the ν_QPO vs. ${\stackrel{\sim }{t}}_{\mathrm{en}}$ curves.

Fig. 7. QPO frequency νQPO versus the energization timescale ${\stackrel{\sim }{t}}_{\mathrm{en}}$ for various values of Ṁinj. Values of Ṁinj, in units of 10−8 M⊙/yr, are 64 (violet), 32 (dark blue), 16 (light blue), 8 (green), 4 (yellow), 2 (orange), and 1 (red). Other parameters are Rc = 109 cm, α = 0.1, and Teff = 4 × 106 K. The black line indicates the relation ${\nu }_{\mathrm{QPO}}\propto t_{\mathrm{en}}^{-1/2}$ given by Eq. (21).

Figure 8 shows the relation of two observables, L_X vs. ν_QPO for various values of Ṁ_inj and R = 10⁹ cm; here the points characterized by strong damping (connected with dashed lines in the previous figure) have been omitted for the sake of clarity. The tendency is that the cases with higher values of Ṁ_inj reach higher luminosities; however, the QPO frequencies cannot exceed values of around three irrespective of Ṁ_inj. Note that high Ṁ_inj cases cannot reach low ν_QPO values due to damping.

Fig. 8. QPO frequency νQPO versus the 2 − 200 keV luminosity LX for various values of Ṁinj. Values of Ṁinj, in units of 10−8 M⊙/yr, are 32 (violet), 16 (dark blue), 8 (light blue), 4 (green), 2 (yellow), 1 (orange), and 0.5. Other parameters are Rc = 109 cm, α = 0.1, and Teff = 4 × 106 K.

As an example of the variations of an output parameter in the ${\stackrel{\sim }{t}}_{\mathrm{en}}-{\dot{M}}_{\mathrm{inj}}$ plane, we show in Figure 9 a heat map of E_pk. There is a tendency for this parameter to increase as ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁ_inj decrease. There are cases where E_pk can enter the MeV regime (light green and yellow patches); however, this is a sensitive function of the reprocessing parameter α. Higher values of α tend to reduce drastically E_pk; see Sect. 5.4.

Fig. 9. Heat map of the average peak energy of the photon distribution Epk as a function of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁinj for Rc = 109 cm. Other parameters are α = 0.1 and Teff = 4 × 106 K. Dark blue regions indicate regimes with low values of Epk.

Figure 10 gives another heat map, this time of the time lag ΔT_sh as a function of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁ_inj for R_c = 10⁹ cm. In Appendix B we give a physical interpretation of these time lags; however, as a rule of thumb, we note that negative time lags are produced from electron cooling while positive ones arise due to contributions from the gradual dominance of the spectrum from photons of the reprocessed component. As we show in Sect. 5.5, ΔT_sh is a sensitive function of T_eff.

Fig. 10. Heat map of the time lag ΔTsh (in units of tcr) as a function of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁinj for Rc = 109 cm. Other parameters are α = 0.1 and Teff = 4 × 106 K. Violet color indicates negative lags, i.e., hard photons come before the soft ones in the prescribed energy bands.

Bolometric luminosities, on the other hand, tend to be larger at high Ṁ_inj and low ${\stackrel{\sim }{t}}_{\mathrm{en}}$, reaching a value of 4 × 10³⁸ erg/s, which is ≃30% of the Eddington luminosity for a 10 M_⊙ black hole, for Ṁ_inj = 10⁻⁸ M_⊙/yr and ${\stackrel{\sim }{t}}_{\mathrm{en}}=1.25$; they are minimum (on the order of 6 × 10³⁶ erg/s) at the opposite corner of the heat map diagrams, i.e., for low Ṁ_inj and high ${\stackrel{\sim }{t}}_{\mathrm{en}}$. The X-ray luminosities follow more or less the same trend, but their ratio to the bolometric luminosities does not remain constant, as already explained.

5.3. Corona radius R_c variations

The present model shows a very simple scaling of the various parameters with R_c, as shown in a simplified analytical way in Sect. 6. For this to work, we need to scale R_c and Ṁ_inj by the same factor. Thus, if both R_c and Ṁ_inj increase by some factor, then ν_QPO decreases while L_bol and L_X increase, all by the same factor. The rest of the output parameters remain unchanged.

This helps us to use the results at R_c = 10⁹ cm as a template and easily derive the output parameters for other corona radii. Thus, following the same procedure as before, we show in Table 5 the results obtained if we keep ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁ_inj constant and vary R_c. Here, we notice that ν_QPO varies roughly as R_c⁻¹, there is a gradient of ΔT_sh and Δϕ from negative to positive values, while both L_bol and L_X are increasing functions of R_c (L_bol varies almost linearly with R_c, while L_X as R_C^1/2). In addition, E_pk increases with R_c and as a result the spectrum becomes flatter.

5.4. Reprocessing fraction α variations

Another parameter that is unknown, but plays a role in QPO formation, is the reprocessing parameter α. Table 6 shows the output parameters of the problem in the case where α is varying.

We notice that changes in α affect most output parameters, but only slightly; that is, despite the fact that α varies by a factor of eight, most parameters change by a factor of two to three, while some, such as ΔT_sh, remain practically unaffected. It is also worth mentioning that even a small value of α is enough to set the system through limit cycles. In this case, however, the electrons have to reach relatively high energies before they trigger the feedback from the reprocessed photons. As a result, the spectra become flatter as α decreases. In general, low values of α tend to produce flat spectra (Γ < 2) and high E_pk, while the opposite holds when α increases. A final comment is that high values of α cause strong damping; however, this also depends on the other parameters of the problem.

5.5. Reprocessing temperature T_eff variations

Another input parameter that has an impact on the spectro-temporal properties of the oscillations is the temperature T_eff that characterizes the reprocessed emission of the accretion disk. Table 7 tabulates the usual parameters for four values of T_eff. Here once again, we notice that most of the output parameters show much slower variations than the input; however, we should mention that ΔT_sh, and by extension Δϕ, show a rather radical change from negative to positive as T_eff increases. This means that phase lags are a sensitive function of the choice of T_eff. This is a result of the exponential cutoff in the Wien part of the reprocessed component; see Appendix B.

5.6. Ratio of escape to energization timescale t_esc/t_en variations

The last input parameter that has an impact on the spectro-temporal properties of the oscillations is the ratio of the escape to energization timescales. As has already been mentioned in Sect. 2, as this ratio increases, the slope of the electron distribution becomes flatter, tending asymptotically to −1. In terms of energy, this translates to more electron energy being transferred to the high end, making the feedback stronger. Table 8 shows that for t_esc = t_en the oscillations are damped, however, for t_esc/t_en > 10 the damping becomes negligible.

6. A simplified analysis of the electron-photon system

In this section we will show in a simplified manner the basic principles behind the natural frequency occurrence discussed in the main body of the paper. As in MPK, we will assume that the electrons and photons do not depend explicitly on energy, i.e., we will treat energy densities instead of differential densities. Then the electron and photon equations can be written

$$\begin{array}{c}\hfill \frac{\partial u_e}{\partial t}=\frac{u_e}{t_{\mathrm{en}}}-\frac{u_e}{t_{\mathrm{loss}}}\end{array}$$(15)

and

$$\begin{array}{c}\hfill \frac{\partial u_{\gamma }}{\partial t}+\frac{u_{\gamma }}{t_{\mathrm{cr}}}=Q_{\mathrm{sources}}.\end{array}$$(16)

Here, t_loss is the characteristic time scale for electron losses while Q_sources denotes the photon energy gains.

We next assume that there are two types of processes: (i) linear, like synchrotron radiation, ICS on disk photons, and so on, and (ii) non-linear like ICS on reprocessed photons; see Eqs. (7) and (8). Thus we write for the electron equation

$$\begin{array}{c}\hfill \frac{\partial u_e}{\partial t}=\frac{u_e}{t_{\mathrm{en}}}-\frac{u_e}{t_{\mathrm{lin}}}-\frac{u_e}{t_{\mathrm{non},\mathrm{lin}}},\end{array}$$(17)

with t_lin and t_non, lin the two loss timescales corresponding to the linear and non-linear processes, respectively. Similarly, we can write for the photon equation

$$\begin{array}{c}\hfill \frac{\partial u_{\gamma }}{\partial t}+\frac{u_{\gamma }}{t_{\mathrm{cr}}}=Q_{\mathrm{lin}}+Q_{\mathrm{non},\mathrm{lin}},\end{array}$$(18)

where, due to energy conservation, each photon rate (Q_lin and Q_non, lin) should be equal to the corresponding loss term entering the electron equation. We proceed next to examine this system of equations according to the combinations of interest in the present paper.

6.1. Purely non-linear case

We begin by ignoring, for the moment, the linear terms. We set $t_{\mathrm{non},\mathrm{lin}}^{-1}\simeq {\sigma }_{T}c(u_s/m_e{c}^2)$ where u_s = αu_γ is the reprocessed soft energy density and α the reprocessing parameter. We normalize the energy densities using ${\stackrel{\sim }{u}}_{e,\gamma }={\sigma }_T{R}_c(u_{e,\gamma }/m_e{c}^2)$ and time using $\stackrel{\sim }{t}=t/t_{\mathrm{cr}}$. Then the set of Eqs. (17) and (18) becomes

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{u}}_e}{\partial \stackrel{\sim }{t}}=\frac{{\stackrel{\sim }{u}}_e}{{\stackrel{\sim }{t}}_{\mathrm{en}}}-\alpha {\stackrel{\sim }{u}}_{\gamma }{\stackrel{\sim }{u}}_e\end{array}$$(19)

and

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{u}}_{\gamma }}{\partial \stackrel{\sim }{t}}=-{\stackrel{\sim }{u}}_{\gamma }+\alpha {\stackrel{\sim }{u}}_{\gamma }{\stackrel{\sim }{u}}_e.\end{array}$$(20)

This set of equations is the classic predator-prey (Lotka-Volterra) system. As is known from theory, see, for example, Strogatz (1994), this system shows undamped oscillations with a natural frequency that is given, for the present case, by ${\stackrel{\sim }{\nu }}_0={\stackrel{\sim }{t}}_{\mathrm{en}}^{-1/2}/2\pi$ or, introducing units,

$$\begin{array}{c}\hfill {\nu }_0=(c/2\pi R_c){\stackrel{\sim }{t}}_{\mathrm{en}}^{-1/2}.\end{array}$$(21)

This is depicted with a black line in Fig. 7. While it matches very well the numerically derived slopes, it overestimates the calculated ν_QPO by ∼50%. Furthermore, since the equilibrium points are ${\stackrel{\sim }{u}}_{e,\mathrm{eq}}={\alpha }^{-1}$ and ${\stackrel{\sim }{u}}_{\gamma ,\mathrm{eq}}={(\alpha {\stackrel{\sim }{t}}_{\mathrm{en}})}^{-1}$, and since, by definition, ${\ell }_{\gamma }\equiv {\stackrel{\sim }{u}}_{\gamma }$, it becomes evident that the compactness ℓ is also scale-free. From Eq. (21) and the arguments given above, it becomes evident that the natural frequency scales as R_c⁻¹ and the luminosity as R_c; see Section 5.3. Furthermore, the dependence found of ${\stackrel{\sim }{u}}_{\gamma ,\mathrm{eq}}$ on the reprocessing parameter explains, at least qualitatively, the calculated behavior of L_bol with α.

6.2. Introducing damping: Synchrotron radiation

The system of Eqs. (17) and (18), in addition to the non-linear term examined above, contains a linear term of the form u_e/t_syn where t_syn = (σ_Tcu_B/m_ec²)⁻¹ with u_B the magnetic field energy density. This term transfers energy from electrons to photons at a constant rate and causes damping in the system, which now has the form

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{u}}_e}{\partial \stackrel{\sim }{t}}=\frac{{\stackrel{\sim }{u}}_e}{{\stackrel{\sim }{t}}_{\mathrm{en}}}-\alpha {\stackrel{\sim }{u}}_{\gamma }{\stackrel{\sim }{u}}_e-d_B{\stackrel{\sim }{u}}_e\end{array}$$(22)

and

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{u}}_{\gamma }}{\partial \stackrel{\sim }{t}}=-{\stackrel{\sim }{u}}_{\gamma }+\alpha {\stackrel{\sim }{u}}_{\gamma }{\stackrel{\sim }{u}}_e+d_B{\stackrel{\sim }{u}}_e,\end{array}$$(23)

where d_B = σ_TR_c(u_B/m_ec²) is the magnetic compactness. From standard non-linear theory we obtain that the natural frequency becomes in this case

$$\begin{array}{c}\hfill {\stackrel{\sim }{\nu }}_0=\frac{1}{2\pi }{[(\frac{1}{{\stackrel{\sim }{t}}_{\mathrm{en}}}-d_B)-\frac{1}{4}{({\stackrel{\sim }{t}}_{\mathrm{en}}{d}_B)}^2]}^{1/2},\end{array}$$(24)

provided that the expression in the parenthesis is greater than 0. The latter gives the condition required for oscillations to occur

$$\begin{array}{c}\hfill d_B<d_{B,\mathrm{crit}}\equiv \frac{2}{{\stackrel{\sim }{t}}_{\mathrm{en}}^2}[{(1+{\stackrel{\sim }{t}}_{\mathrm{en}})}^{1/2}-1].\end{array}$$(25)

From here we can obtain a limiting value of B

$$\begin{array}{c}\hfill B_{\mathrm{cr}}=1.7\times {10}^5{R}_9^{-1/2}{\stackrel{\sim }{t}}_{\mathrm{en}}^{-1}{[{(1+{\stackrel{\sim }{t}}_{\mathrm{en}})}^{1/2}-1]}^{1/2}G,\end{array}$$(26)

above which the system becomes overdamped. We also note that since the damping rate is $\lambda =\frac{1}{2}{d}_B{\stackrel{\sim }{t}}_{\mathrm{en}}$, small values of d_B (for given ${\stackrel{\sim }{t}}_{\mathrm{en}}$) produce negligible damping, and this is the underlying assumption of the current paper where the synchrotron has been neglected. The damping increases with d_B, until it reaches the limiting value d_B, crit above which the oscillations are quenched. A numerical treatment of the effects of synchrotron radiation is given in Appendix C.

6.3. Introducing damping: ICS on disk photons

Assuming that in addition to the reprocessed photons, the disk radiates its own photons, which have energy density U_ds (see Sect. 3), the system becomes identical to Eqs. (22) and (23) with ${\stackrel{\sim }{u}}_{\mathrm{ds}}$ replacing d_B. For t_en = 2.5 we get ${\stackrel{\sim }{u}}_{\mathrm{ds},\mathrm{mx}}=0.3$, which can be easily compared to the case shown in Fig. 2 where the damping becomes significant for ${\ell }_{\mathrm{ds}}\equiv {\stackrel{\sim }{u}}_{\mathrm{ds}}\simeq 0.1$.

7. The case of GRS 1915+105

Although it would have been tempting to use the method outlined in Sect. 5 and try to produce detailed χ² fits to the observational data, the number of free parameters would not have offered a physical insight into the underlying mechanism of the source. In contrast, it would be more interesting to see whether a systematic change in some fundamental quantity of our model would produce results that could be compared favorably to observational data, because this could give us a hint of the physics of the source. Therefore, we have tried to compare the data with trends produced in our model by systematically changing one of the free parameters. As an example, we take the case of GRS 1915+105, which is a well-studied source; see, for example, Morgan et al. (1997), Trudolyubov et al. (1999), Reig et al. (2000), Qu et al. (2010). Recently, Zhang et al. (2020) analyzed 620 RXTE observations of the source and found that the QPO phase lag decreases with QPO frequency and changes sign from positive to negative at a frequency around 2 Hz.

The simplest approach to compare our model to the above results is to consider changes in only one of our parameters. Since the observations indicate that ν_QPO changes by about one order of magnitude, this could mean that it could be achieved either with variations in ${\stackrel{\sim }{t}}_{\mathrm{en}}$ or in R_c. The phase lag, however, changes from negative to positive by almost 1 rad and a brief inspection on Table 3 indicates that this cannot be achieved by changes in ${\stackrel{\sim }{t}}_{\mathrm{en}}$ alone. Therefore, we consider changes in R_c and the results are shown in Fig. 11. This shows the data from Zhang et al. (2020) (black triangular points) along with various lines produced with the present model.

Fig. 11. Phase lag Δϕ versus the QPO frequency νQPO. The violet line corresponds to Ṁinj = 8 ×10−8 M⊙/yr for Rc between 109.3 cm and Rc = 108.55 cm, the dark blue line corresponds to Ṁinj = 16 ×10−8 M⊙/yr for Rc between 109.6 cm and 108.55 cm, while the light blue line corresponds to Ṁinj = 32 ×10−8 M⊙/yr for Rc between 109.9 cm to 108.7 cm. In all cases the radius Rc is decreasing from left to right. The other parameters are ${\stackrel{\sim }{t}}_{\mathrm{en}}=2.5$, α = 0.4, and Teff = 4 × 106 K. The red lines show the results of varying ${\stackrel{\sim }{t}}_{\mathrm{en}}$ by a factor of 2 (i.e., between .25 and 5) when Ṁinj = 16 ×10−8 M⊙/yr and the corona radius is Rc = 109.3 cm,  109 cm, and 108.7 cm (l to r). As ${\stackrel{\sim }{t}}_{\mathrm{en}}$ decreases, the curves move up and to the right. Orange lines repeat the same procedure with Teff varying by a factor of 2 (i.e., between 2 and 8 × 106 K). As Teff increases, the curves move upwards. The data points are taken from Zhang et al. (2020).

The violet, dark and light blue lines represent the results obtained for Ṁ = 8, 16 and 32 × 10⁻⁸ M_⊙/yr, respectively, keeping all other parameters constant but varying R_c. For example, the dark blue line curve shows the results obtained by changing the corona radius from R_c = 10^9.6 cm (upper left side) to R_c = 10^8.55 cm (lower right side) using a logarithmic step of 0.15 between consecutive runs. Having obtained this particular curve, the others can be scaled according to the findings of Sect. 5.3. Therefore, the pair (ν_QPO, Δϕ) at R_c = 10^9.6 cm for Ṁ = 16 × 10⁻⁸ M_⊙ yr will be mapped at (ν_QPO/2, Δϕ) for R_c = 10^9.9 cm when Ṁ = 32 × 10⁻⁸ M_⊙ yr. Similarly, it will be mapped at (2ν_QPO, Δϕ) for R_c = 10^9.6 cm when Ṁ = 16 × 10⁻⁸ M_⊙ yr.

Having obtained a general trend that shows an overall similarity to the observations, the next step is to investigate how the system responds to changes in the other parameters. For this we have chosen three radii (R_c = 10^9.3 cm,  10⁹ cm, and 10^8.7 cm) for case Ṁ = 16 × 10⁻⁸ M_⊙ yr and increased and/or decreased ${\stackrel{\sim }{t}}_{\mathrm{en}}$ by a factor of two while keeping all other parameters fixed. The results are depicted with red lines in Fig. 11. The red lines moving up and to the right indicate a decrease in ${\stackrel{\sim }{t}}_{\mathrm{en}}$ – see Table 3.

We can repeat the same procedure for showing changes in the Δϕ − ν_QPO plane caused by variations in T_eff. So we performed runs increasing and/or decreasing T_eff by a factor of two while keeping all other parameters fixed. The results are shown with an orange line in Fig. 11. In all three cases shown, Δϕ increases with increasing T_eff. Furthermore, we note that the orange lines are almost vertical, i.e., temperature changes, at least in the range specified here, do not produce significant changes in ν_QPO; however, they have a significant effect on Δϕ. This means that the present model is rather sensitive to the reprocessed disk temperature, a fact already mentioned in Sect. 5.5.

Most of the observation points lie between the lines of Ṁ = 16 and 32 × 10⁻⁸ M_⊙/yr. Alternatively, one could assume that the spread of the data is caused by a constant accretion rate and fluctuations of a factor of two in either the energization timescale or the reprocessed disk temperature.

Figure 12 shows the calculated fractional RMS versus ν_QPO (colored lines) plotted over the GRS 1915+105 data (black symbols). The colors of the lines correspond to the cases examined in Fig. 11. Once again, most points lie between the dark blue and light blue lines, i.e., Ṁ = 16 and 32 × 10⁻⁸ M_⊙/yr. However, we caution that the fractional RMS was calculated from the light curves using a Python routine; also, the shown value should be considered as an upper limit because the RMS was calculated using only a few cycles of the light curves. Therefore, the results shown here are only indicative.

Fig. 12. Fractional RMS versus the QPO frequency νQPO for the cases depicted in Fig. 11. The fractional RMS increases with decreasing ${\stackrel{\sim }{t}}_{\mathrm{en}}$ (red lines) and increasing Teff (orange lines). The data points are taken from Zhang et al. (2020).

Clearly, we could expand the approach taken above by varying two or more of the input parameters. For example, we could argue that by increasing Ṁ_inj, R_c decreases, and then combine the results shown in the corresponding tables to obtain relations between the various observables. This, however, requires some further physical justification and it lies outside the aims of this paper.

Finally, it is interesting to compare our results with those of Karpouzas et al. (2021) who, however, used a different approach to the time-dependent Comptonization model with feedback. Both models require large corona radii (on the order of 10³R_g), assuming a black hole mass of 12.4 M_⊙ (Reid et al. 2014) especially when ν_QPO is small. Furthermore, both find that ν_QPO increases as R_c decreases. However, in the Karpouzas et al. (2021) model, R_c starts to increase again as ν_QPO increases past 2 Hz, while in our case R_c continues to decrease monotonically. Nevertheless, it should be mentioned that in Karpouzas et al. (2021) the emphasis was placed on simultaneously fitting the phase lag and/or RMS data, so all parameters changed to obtain the best χ²; however, in our case the emphasis was placed on the trends obtained by only changing R_c while keeping the rest of the parameters fixed.

8. Summary and discussion

The aim of this paper was to verify and extend the results of MPK, i.e., to show that QPOs arise naturally in a disk-corona system. For this we applied a one-zone model for the evolution of electrons and photons inside the corona and solved the coupled kinetic equations for the two species distribution functions. We assumed that the electrons were energized by some unspecified mechanism and that they were losing energy by ICS on soft photons. The origin of these photons was assumed to come from the reprocessing of the electrons inverse Compton radiation on the disk. This simple prescription is adequate to make the system go through limit-cycle behavior. As shown in Appendix A, as electrons gain energy, their ICS radiation increases, and so does the reprocessed radiation. At some point, the latter has increased beyond a critical point and is able to overtake energization, thus cooling the electrons efficiently. This leads, in turn, to a reduction of both the ICS radiation and its reprocessed radiation; therefore, the electrons are able to gain energy once again and the cycle repeats itself. A necessary condition for the occurrence of these cycles is that the disk has low luminosity, so the soft radiation is dominated by reprocessing and not by the disk thermal emission. A luminous disk produces plentiful photons by itself, which linearizes the system as they cause steady electron ICS losses that effectively balance the electron energization. This qualitatively agrees with the observations because QPOs do not appear when BHXRBs are in the high soft state; see, for example, Remillard & McClintock (2006).

Our approach aimed to examine the feasibility of the limit cycle occurrence when the full energy dependent electron and photon distributions are employed. As key parameters, we considered the electron energization timescale t_en, the mass injection rate Ṁ_inj, and the radius of the corona R_c where the heating and cooling of the electrons is supposed to take place. Other auxiliary parameters include the reprocessing factor α and the temperature of the reprocessed photons T_eff. Using a wide range of generally accepted values for the input parameters, with the exception of t_en which we treated phenomenologically, we found that the vast majority of the cases show some type of oscillatory behavior; only very fast t_en (on the order of t_cr) and/or very high values of Ṁ_inj do not favor oscillatory solutions. Furthermore, we found that the frequency of oscillations ν_QPO ∝ R_C⁻¹ and ν_QPO ∝ (t_en/t_cr)^−1/2 while it has a weak dependence on the other input parameters; see Sects. 5.3 and 5.1, respectively, as well as the analytical estimates presented in Sect. 6.

As emphasized earlier, in contrast to MPK, the present approach can provide spectral information on the photon distribution. Taking the average of the output quantities during a few cycles, we found that the X-ray spectra showed either hard indices coupled with high energy spectral peaks or softer indices coupled with lower peaks², the corona optical depth in cold electrons was consistently small (τ_T < 0.1), while the spectral peaks varied between a few tens to a few hundreds of keV, which corresponds to electrons reaching maximum Lorentz factors of ≃5 − 50. Generally, faster acceleration showed spectra with higher peaks and higher bolometric luminosities, as intuitively expected; see Table 3 and Fig. 5.

Another characteristic of the present model is that it can produce both soft (negative) and hard (positive) lags. Appendix B gives a detailed explanation for this, but generally speaking, one could associate soft lags with electron cooling, while hard lags come from the contributions of the reprocessed component on the tail of the photons produced from uncooled electrons. This gives the prediction that time lags measured between high energy regimes, where contributions from the disk are not expected, have to be always soft. Generally speaking, the soft and/or hard trend is that fast acceleration and low Ṁ_inj tend to produce soft lags, while the situation reverses for slow acceleration and high Ṁ_inj; see Fig. 10.

The results of the model were compared to the well-studied source GRS 1915+105. We found that the trend observed in the QPO phase lag – QPO frequency can be satisfactorily reproduced by varying the corona radius R_c while keeping the other input parameters fixed; see Sect. 7. However, as the model is phenomenological, it offers no physical explanation why R_c should change by almost a factor of ten, while Ṁ_inj varies only by a factor of two.

Apart from ICS, we have augmented our calculations by using two other relevant physical processes, namely photon-photon pair production and electron-proton bremsstrahlung. These remain marginal for the vast majority of the cases considered. The former can play some role for those cases when the energization pushes the electrons well beyond the GeV range but these require very fast ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and low Ṁ_inj; see bottom left corner of Fig. 9. Generally, since photon compactnesses, as defined in Eq. (10), remain well below unity and cooling is dominated by collisions in the Thomson regime, this process is considered to contribute minimally in the vast majority of the cases. The same can be said for electron-proton bremsstrahlung. This process can cause changes in the behavior of the system because it produces radiation that acts as a background photon field, so, practically speaking, it has the same effect as the disk radiation; see Sect. 6. However, its inclusion does not change the essence of our results, at least for the set of parameters examined here.

This study has neglected the effects of synchrotron radiation. However, as shown in Sect. 6 and Appendix C, synchrotron radiation introduces a damping effect and for values of B-field B ≃ B_cr (see Eq. (26)) it can impede the oscillatory behavior of the system. Therefore, all results presented here imply that the B-field of the corona must be less than B_cr, which depends on ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and R_c.

Our approach placed emphasis on non-thermal processes because the compactnesses and the associated optical depths of low energy electrons are lower than 0.1 for most of the parameters studied. For higher compactnesses, one should include thermal processes and treat a hybrid thermal and/or non-thermal plasma (Poutanen & Coppi 1998). However, the non-linearity of the system does not depend so much on the nature of the physical mechanisms used, but it is embedded in the equations to be solved. Once we treat the system of electrons-photons self-consistently, allow the electron cooling to occur on reprocessed photons, and include explicitly an energization method for electrons that is able to replenish the energy lost by the ones which have cooled, then, inevitably, we come to the Lotka-Volterra type of equations which show the limit cycle behavior, as found both in MPK and in this paper³. The only way to effectively damp the oscillations is to allow some physical mechanism, which can have a linearizing effect, to operate on a comparable timescale to ICS.

Several papers deal with time-dependent Comptonization to probe the origins of spectral variability and to reproduce the trends of the rms and lags with QPO frequency either in accreting neutron stars (Lee & Miller 1998; Kumar & Misra 2014; Karpouzas et al. 2020) or BHXRBs (Karpouzas et al. 2021; García et al. 2021; Bellavita et al. 2022). Although these works use the same assumptions as ours, they differ in the sense that they first assume the QPO frequency and consequently calculate the quantities of interest. Furthermore, the electron heating rate, when used, is given in each case separately as an external parameter, i.e., the papers above do not use an equation for the electrons that evolves simultaneously with the one for the photons. In our case, we simply write and solve the system of the electron and/or photon equations, and there is no guarantee a priori whether the solution will go to limit cycles or to a steady state; therefore, there is a conceptual difference in the two approaches in both methods and aims.

The one-zone model employed in the present analysis is the standard framework that is used in modeling the spectra of both active galactic nuclei and/or blazars and gamma ray bursts, the latter both in the prompt and afterglow phase. It has the advantage that it can be fully time-dependent and it is able to treat self-consistently the gains and losses of electrons as well as the sources and sinks of photons. The only change that our approach brings to this standard model is to allow for the soft photon feedback. Furthermore, we have assumed that electron gains come through a continuous process (see Eq. (5)), which we have simply characterized by a timescale t_en. This allows the electrons to replenish their energy, a feature that is central to our approach. In our case, we have chosen a systematic acceleration of electrons; other approaches that are also possible, but have not been treated here, include stochastic acceleration (Li et al. 1996; Stawarz & Petrosian 2008) or an injection term of readily accelerated particles as is usually the approach in blazars (Mastichiadis & Kirk 1997)⁴. Finally, it should be mentioned that the present setup of the problem prohibits the use of t_en ≤ t_cr.

The disadvantage of the one-zone model is that it assumes that both species are uniformly distributed inside the corona, which, in addition, is assumed to be homogeneous and static. Furthermore, it can treat only a spherical geometry, and it cannot take into account any details of the disk-corona configuration. Here we have simply assumed that the corona is spherical and it either enshrouds the disk or it is above it, i.e., the lamp-post model.

This paper can be used to model various black hole systems, both galactic and extragalactic⁵. However, it needs to be implemented with some physically motivated mechanism for particle energization and to include more radiation mechanisms such as synchrotron radiation and, depending on the assumptions, thermal mechanisms. Other simplifications, like the uniform corona assumed here, can be more difficult to address as multizone models have not been satisfactorily developed yet. However, the model can accommodate time-delay effects, i.e., it can allow for a finite time-interval for reprocessing to occur. Similarly, the effects of a corona moving outward with a mild relativistic velocity should not, in principle, change the basic concept because there is still the possibility of a disk-corona feedback (Reig & Kylafis 2021). At any rate, the results of this analysis show that disk-corona feedbacks are robust and, thus, promising to explain some of the vast BHXRBs phenomenology.

Acknowledgments

The author would like to thank Dr. S. Boula for many discussions, substantial support on technical issues and a thorough reading of the final manuscript; Dr. N. Kylafis for many discussions and hospitality at the University of Crete where this project was first conceived; Drs. M. Mendez and M. Petropoulou for discussions and many helpful comments on the manuscript; Dr. G. Vasilopoulos for discussions and technical support throughout; and, last but not least, an anonymous referee for their comments that helped improve the present work. The following Python libraries were used: Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020). For Fig. 1 we used Google (2025).

References

Appendix A: A physical interpretation of the limit cycle behavior

In order to gain a physical understanding of the limit cycle behavior, one needs to follow in detail the spectro-temporal evolution of the particle and photon distributions. The exact values of the input parameters are not important for this particular example because, qualitatively, all cases showing oscillatory behavior go through the same phases. However, we have chosen a case that shows strong limit cycles in order to make the figures clearer.

Figure A.1 shows a close up of the first two limit cycles of a generic case where ℓ_ds = 10⁻⁶, a value that we used throughout in the main text. The red and green lines separate the time interval according to whether electrons lose or gain energy as this will be evident further down. Note, however, that these do not coincide with the maxima and the minima of the bolometric photon light curve. We proceed now to analyze this behavior of the light curve by examining in detail the electron and photon distribution functions at each of these time intervals.

Fig. A.1. Generic light-curve showing the first two peaks of the bolometric ICS compactness. The red and green lines indicate the start of the phases where electron cooling and re-energization prevails respectively; for details see text.

1. Initial electron energization phase: This phase is physically artificial, because it depends strongly on the initial conditions which have been set equal to zero for both electrons and produced photons. Nevertheless, it has to be included in our analysis because it demonstrates some key features of the model.

Electrons enter the energization process at low Lorentz factors (γ = γ_inj) and, as time passes, they increase their energy according to Eq. 2. This is shown in Fig. A.2 (left) that draws snapshots of the successive electron distribution functions, which is drawn in γ²n_e units⁶, at every 2t_cr with the red line connecting their peaks, which occurs at γ_pk(t). The distributions gradually form a power law which is close to n_e ∝ γ⁻¹ since we have assumed, as in the main text, that t_esc ≪ t_en; see text below Eq. 5. The corresponding photon spectra, depicted in Fig. A.2 (right) for the same time intervals as the electrons, originally consist only of the disk photons which have a gray-body distribution of a very low compactness (ℓ_ds = 10⁻⁶ by assumption). However, as the energy γ_pk increases, they start producing increasingly harder ICS photons, according to the Thomson relation x_pk(t) = γ_pk(t)²x₀, where x₀ is the peak of the reprocessed component, the reprocessing of which increases the soft photon density according to Eq. 7. This positive feedback causes an exponential outgrowth of both hard ICS and reprocessed soft photons. The red line that connects the consecutive peaks of the ICS component moves, therefore, up to higher energies and to higher luminosities, with the former being linear and the latter exponential. This will last up to the instant when the soft photons have reached such densities that can cause substantial losses which overcome the electron energization, i.e. the loss term in the bracket of Eq. 3 becomes larger than the gain term. From that instant the electrons enter a phase where cooling becomes dominant. This transitional time is shown as a vertical red line in Fig. A.1 and it is a universal feature of the present model.

Fig. A.2. Snapshots of a generic case showing the electron distribution ne(γ) multiplied by γ2 (left panel) and the emitted photon spectra in x.dℓ/dx units (right panel) taken every 2tcr during the initial particle energization phase. Here x = ϵ/mec2. The photon light curve of this case is given in Fig. A.1 from t = 0 up to the first red vertical line. Red lines in the present figure indicate the peak of the electron instantaneous distribution (left) and of the ICS photons (right). As time passes both peaks move from left to right.

2. Electron cooling phase: Figures A.3 (left) and (right) depict the electron and photon spectral evolution, respectively, during the phase where ICS losses exceed the electron energization gains. Snapshots here are shown at every t_cr because the system goes rapidly through this phase. At early stages, the electrons close to the peak suffer severe losses and release their energy producing even more ICS photons, further intensifying their cooling. Thus, electrons cool abruptly (i.e. within a few t_cr), γ_pk drops, and, as a result, the rate of ICS and reprocessed photons is reduced. The red line in Fig. A.3 is now moving from right to left and as the cooling is gradually reduced electrons become almost stationary at a new, lower γ_pk = γ_cool. Note that the electron cooling is incomplete, i.e. the soft photon burst cannot force them to cool all way down to their rest mass, so usually 1 ≪ γ_cool. It is also worth mentioning that the electron distribution becomes slightly harder due to the cooled high energy electrons that have moved to lower energies causing a pile-up effect.

Fig. A.3. Same quantities as in Fig. A.2 for the electron cooling phase, i.e. for the time interval between the red and green lines of Fig. A.1. Snapshots are taken at every tcr. The violet lines indicate the distribution functions at the peak of the photon light curve. Light blue lines depict the distribution functions for times before the peak, while the dark blue lines for times after the peak. As time passes, the peaks of both distributions move from right to left.

The behavior of the photons during this phase is more complex. At its early stages, the photons keep increasing in luminosity (but not in energy) as electron cooling reaches a maximum. The peak of the photon light curve is typically achieved a few crossing times after the saturation of electrons. The light blue lines in Fig. A.3 (right) correspond to the outgrowth of photons at these early stages of the electron cooling phase. These still have a spectral shape similar to the ones during the earlier energization phase. After the photon bolometric peak (depicted here with violet color in both panels) the photon spectra decrease in energy density and are dominated by (i) photon escape at higher energies as the electrons producing them have cooled and (ii) fresh ICS production at lower energies produced by the cooled electron distribution at the particular time instance; these spectra are shown in Fig. A.3 (right) with dark blue lines. The red line showing the successive peaks of the photon distribution (in x.dℓ/dx units) moves in a characteristic counter-clockwise trajectory. The continuous photon escape from their earlier high stages in combination with the low current ICS rate causes the cooling rate to diminish and eventually to become smaller than their energization rate. From then onward, the system moves back to the phase where electrons gain energy. This time instant is shown as the first green vertical line in Fig. 2.

3. Electron re-energization phase: Physically, this is very similar to the initial energization phase as the electrons once again move from low to high energies. There are, however, two key differences. (i) Electrons do not start from γ_inj but from γ_cool, i.e. the value that photons forced them to cool in the previous phase and which depends on the choice of the initial parameters; see Fig. A.4 (left). Furthermore, their slope becomes slightly steeper tending again to p = −1 as the cooled electrons are re-energized. (b) Early in this phase, higher energy photons continue to drop in luminosity because their spectrum is still dominated by escape; their spectra are shown in Fig. A.4 (right) with dark blue lines. However, gradually a new component of ICS photons produced from the freshly re-energized electrons starts appearing; this is shown here with light blue lines. This population dominates at later times and the superposition of these two (old and fresh) components produces the spectra shown in Fig. A.4 (right). The red line depicting the position of the peak photon luminosity moves in a counter-clockwise convex trajectory. As before, this phase will end once ICS photons are built in high enough densities so that the reprocessed photons are able to overcome electron energization and cool them to lower energies. The cycle therefore will repeat all over again.

Fig. A.4. Same quantities as in Fig. A.2 for the electron re-energization phase, i.e. for the time interval defined between the green and the second red line of Fig. A.1. Snapshots are taken at every 2tcr. Light blue lines depict the photon distribution function for times before the photon minimum, while the dark blue lines for times after the minimum. As time passes, the peaks of both distributions move from left to right.

A different way to see the above is to plot again Figs. A.2 - A.4, but now as a function of time, for various characteristic electron and photon energies. Figure A.5 plots the quantity γ²n_e vs. normalized time for various electron energies for the same parameters that produced the previous figures. The left panel shows the evolution of the system up to a time as to include the first two photon peaks. As the initial conditions are n_e(γ, 0) = n_γ(x, 0) = 0 one sees electrons initially to increase and gradually to obtain a steady state with the lower energies getting there before the higher ones. Once the soft photons increase and cause electrons to cool, then it is the higher energies that cool faster; see right panel that shows a zoom of the same electron bins around the first photon peak. Note that it is possible that the higher energies cool before photons reach their peak. Note also that the electrons do not cool completely but only down to a certain energy γ_cool. As already mentioned, this occurs because the ICS cooling timescale depends on energy and photons decrease in density before there is enough time to cool the electrons completely. This can be seen in Fig. A.5: only electron energies with time curves that show dips suffer cooling as their energy is lost to ICS radiation. Therefore, only the cooled electrons play an active role in the photon spectrum and this explains the behavior that one observes in the next figure.

Fig. A.5. Electron density function ne multiplied by γ2 for various electron energies as a function of normalized time for the run with the parameters that produced the light curve shown in Fig. A.1. Red line corresponds to γ = 100.3, orange to 100.6, yellow to 100.9, green to 101.2, light blue to 101.5, dark blue to 101.8 and violet to 102.1. The bolometric photon light curve of Fig. A.1 is shown with black line. The right panel shows a zoom of the first peak for clarity.

Figure A.6 depicts the ICS photon light curves for seven characteristic energies x_i that correspond to electron energies shown in Fig. A.5 through the Thomson relation x_i = γ_i²x₀. The light curves are plotted around the first peak of the bolometric ICS luminosity, depicted here with black line. Several characteristic features become readily apparent: Not all photon energies peak at the same time but the lower ones tend to lag from the higher ones. This is to be expected because of the differential cooling of the electrons discussed above. The lag becomes zero at lower photon energies that correspond to the uncooled electrons. Therefore, in this example, lags start appearing at the green curve that corresponds to the minimum electron energy that shows some cooling; see Fig. A.5. Furthermore, it is a small width of energies (less than one order of magnitude) that contribute to the bolometric ICS light curve. Low energies do not contribute because the electrons producing them have not cooled and thus do not radiate substantially.

Appendix B: A physical interpretation of time lags

The present model shows a rich behavior of soft and hard time lags because the photon spectrum is dominated at high energies by ICS produced from electron cooling, at lower energies by ICS of uncooled electrons, while at even lower energies by the reprocessed soft photons. The two components (ICS and reprocessed soft photons) dominate the spectrum at different energy regimes while there is a rather narrow intermediate energy regime that both contribute partially. As an illustrative example we show in Fig. B.1 the time lags that have been derived from the first light curve peak for the case examined in Appendix A. Here the time lags ΔT_i have been calculated for the first cycle using the relation

$$\begin{array}{c}\hfill \mathrm{\Delta }{T}_i=T_{pk,i}-T_{pk,rpr}\end{array}$$(B.1)

where T_pk, i corresponds to the instant at which the light curve of energy E_i peaks and T_pk, rpr corresponds to the peak of the reprocessed distribution. The time resolution is 0.1t_cr as in the rest of the paper. For a better understanding of the spectral behavior at the various energy regimes, we plot in the same figure the average photon spectrum during that cycle.

• At high energies (in the example here approximately between 10² to 10⁴ keV) the time lags are soft because the monoenergetic photon light curves are dominated by the eletron cooling; see Fig. A.3. Since the energetic electrons cool faster, one expects that harder photons will appear before softer ones, thus the appearance of a soft lag. Note that the range of photon energies corresponds to γ_cool²E₀⁷ and γ_pk²E₀ where γ_cool and γ_pk are the two limiting electron Lorentz factors as shown in Figs A.3 and A.4.

• Below a characteristic photon energy E_cool = γ_cool²E₀ the spectrum will be dominated mainly by the tail of the ICS distribution produced by electrons with Lorentz factors γ_cool. These will have a power law of index +2 in xdℓ/dx units (see, e.g. (Blumenthal & Gould 1970)) and they will oscillate in phase since their radiation is dominated by electrons of Lorentz factors γ ≃ γ_cool. Therefore, their relative time lags will be zero and the plot will have a horizontal branch. In our example this occurs roughly between 10 and 100 keV. Note also that in the vast majority of the cases examined in the present paper, the peaks of this horizontal branch come after the peak of the reprocessed component. This can be seen in Fig. A.6, where the curves at low energies (red, orange, yellow, green) come after the peak of the bolometric luminosity (black line). This has its significance as we show below.

• At even lower energies (in the current example for E < 10 keV), contributions from the reprocessed distribution start to become important. Therefore, the spectrum will consist of two components: the tail of the ICS component mentioned above and the Wien part of the gray-body spectrum. As photon energies decrease, the contribution of the latter becomes gradually more important than that of the former. Thus, since ΔT_i > 0 for the light curves of the horizontal part, as mentioned earlier, the time lag ΔT_i will show a decreasing trend with decreasing photon energy as it has to reach, by definition, the value ΔT_i = 0 at E = E₀. Therefore this part of the spectrum will show a hard lag.

• Finally, it is worth mentioning that at very high energies, there is another hard lag as energization overcomes losses for electrons; however, both electron and photon distributions have very low values there as they are deep in their high energy cutoffs so one should not expect some observational effect, at least for the parameters used in the present paper.

Fig. B.1. Photon spectrum x.dℓ/dx for various photon energies as a function of normalized time around the first peak of the light curve shown in Fig. A.1. The photon energies correspond approximately to the characteristic radiation of the electrons through the ICS relation in the Thomson limit xi = γi2x0. For color codes see Fig. A.5.

The above analysis covers the basic trends to be expected in a plot of ΔT_i vs. E. However, there are some variations which depend on the initial parameters. The most common is that for parameters that favor a strong cooling, the photon energy E_cool which corresponds to γ_cool, moves to the left and, as a result, the horizontal branch will shrink or disappear altogether. In such cases, the soft lag branch gives its place immediately to a hard lag as the photon energy decreases.

The calculation of ΔT_i as outlined above allows us to estimate at once the lags between any energy band. For instance, in Fig. B.1 we have plotted with vertical black lines the energy bands 2-6 keV (soft) and 6-15 keV (hard) that many authors use to derive lags in QPO observations. By integrating over the corresponding energy bins we can find the total energy entering the two bands at each instance and use this to compute time lags (defined as ΔT_sh in the main text) and the corresponding phase lags. Thus, for example, the case examined above shows a hard lag of ΔT_sh = 0.4t_cr. The average values⁸ of the two lags are depicted with black horizontal lines in Fig. B.1. If, for different parameters, the cooling was increased, then the whole ΔT curve would have moved to the left and different parts of the curve would have entered the soft and hard bands bringing a change to ΔT_sh. This can be seen in Fig. B.2 which depicts two such cases. Due to the aforementioned shifting, the orange one shows a positive (hard) lag and the violet one a negative (soft) one.

Fig. B.2. Time lags ΔTi of monochromatic photon energy peaks with respect to the time where the black body peak occurs (red line). The average multiwavelength spectrum for the same cycle has also been plotted with green line. The vertical black lines indicate the "soft" (2-6 keV) and "hard" (6-15 keV) energy bands sometimes used in the literature to measure phase lags. The horizontal black lines indicate the average lags in these two bands. In this case the relative soft-hard lag is positive.

Fig. B.3. Time lags ΔT of monochromatic photon energy peaks with respect to the time where the peak of the reflected component occured. The violet line depicts a negative lag, while the orange one a positive. The vertical black lines indicate the "soft" (2-6 keV) and "hard" (6-15 keV) energy bands. The horizontal colored lines indicate the average lags in these two bands.

The fact that hard lags appear in the present model due to the contributions of the Wien part of the gray-body spectrum of the reprocessed component can be seen in Figure B.3 that shows time lags and multi-wavelength spectra as in Fig. B.1, but this time for a delta function reprocessing component. As it can be readily verified, in this case there is no rising part in the ΔT vs. E plot; instead, the so-called "horizontal" branch reaches exactly at the energy of the reprocessing component.

Fig. B.4. Monoenergetic time lags (red line) and the multi-wavelength spectrum (green line) in the case where the reprocessed photons are a delta-function.

Finally, we should mention that the temperature of the reflected component T_eff plays an important role in the sign of the lag. Figure B.4 compares the time lags when T_eff takes different values, from T_eff = 1 × 10⁶K (red) up to T_eff = 8 × 10⁶K (violet) which correspond to the runs shown in Tab. 7 of the main text. The corresponding curves are shifted with respect to each other, different parts of each enter the soft and hard energy bins, and, as a result, the sign of the lags change from negative to positive as T_eff increases.

Fig. B.5. Time lags ΔT of monochromatic photon energy peaks with respect to the time where the peak of the reflected component occurred. The temperatures of this component were assumed to be 106K (red), 2 × 106K (orange), 4 × 106K (green) and 8 × 106K (violet) and correspond to the cases shown in Tab. 7 of the main text. The vertical black lines indicate the "soft" (2-6 keV) and "hard" (6-15 keV) energy bands.

The above result mean that not only the temperature of the reprocessed component, but its shape as well, might play a role in the sign of the lags. So, for example, one could use a truncated power law instead of a gray-body spectrum and repeat the analysis above. However, we will not deal with such cases here, but we note that Bellavita et al. (2022) have reached the same conclusion.

Appendix C: Effects of synchrotron radiation

As shown in Sect. 6, synchrotron radiation can have a damping effect on the limit cycles, provided that the magnetic field exceeds some critical value given by Eq. 26. This value, however, is only indicative, as it was derived from solving the simplified system of equations presented in that section. In the present Appendix we show the effects of synchrotron on the output parameters when the full code is used. Tab. C.1 tabulates the usual parameters having as the only variable the corona B-field value while keeping the rest parameters fixed. To determine indicative B-values we started from $B_{\mathit{crit}}({\stackrel{\sim }{t}}_{\mathit{en}})$ and we multiplied it by 1/4,  1/2 and 1. As can be seen from the values of ${\stackrel{\sim }{T}}_{10}$ and f_rms in the aforementioned table, damping sets in practically for B = B_crit/2 and quenches oscillations for B = B_crit⁹. Furthermore, both L_bol and L_X are reduced with increasing B, as synchrotron takes effectively electron energy and channels it to lower photon frequencies, which, however, are efficiently absorbed due to synchrotron self-absorption. Finally, for B ≥ B_cr synchrotron losses largely impede electron energization, leading to steady-state and low corona luminosity.

Appendix D: Relation to previous work

We proceed now to show the relation of the present work to MPK. While both works rely on the same principles and essentially come to the same conclusions, they are based on different assumptions.

MPK used for powering the electrons a term which was the luminosity related to the mass accretion rate (see their Eq. 2). Here, on the contrary, we have defined an injected mass accretion rate Ṁ_inj that is not related directly to the luminosity, but was used instead to define the injection rate of electrons in the energization process; see Eq. 9. However, since one of our output parameters is the bolometric luminosity L_bol, we can calculate the mass accretion rate of MPK from the relation Ṁ_MPK = L_bol/L_Edd and obtain Ṁ_MPK a-posteriori. Therefore, this offers a way to compare quantitatively the results of the two papers. Note, however, that for the present setup, the black hole mass is not required; it can be used only indirectly to calculate the expected L_Edd and to verify that R_c > > R_s, where R_s is the Schwarzschild radius.

Our numerical results show that both ν_QPO and L_bol can be satisfactorily represented as power laws of the energization timescale ${\stackrel{\sim }{t}}_{\mathit{en}}$(see Tab. 3 and Fig. 7) and corona radius R_c (see Tab. 5). Thus we can write

$$\begin{array}{c}\hfill {\nu }_{\mathit{QPO}}=C_1{\stackrel{\sim }{t}}_{\mathit{en}}^{{\alpha }_1}R{{}_c}^{{\alpha }_2}=C_1^{\prime }{t}_{\mathit{en}}^{{\alpha }_1}R{{}_c}^{{\alpha }_2-{\alpha }_1}\end{array}$$(D.1)

and, similarly,

$$\begin{array}{c}\hfill L_{\mathit{bol}}=C_2{\stackrel{\sim }{t}}_{\mathit{en}}^{{\beta }_1}R{{}_c}^{{\beta }_2}=C_2^{\prime }{t}_{\mathit{en}}^{{\beta }_1}R{{}_c}^{{\beta }_2-{\beta }_1}\end{array}$$(D.2)

where C₁,  C₁′, C₂,  C′₂ are proportionality constants; we have also made use of the relation ${\stackrel{\sim }{t}}_{\mathit{en}}=t_{\mathit{en}}/t_{\mathit{cr}}$. If we eliminate t_en between the above relations, we find

$$\begin{array}{c}\hfill {\nu }_{\mathit{QPO}}=C_3{L}_{\mathit{bol}}^{{\alpha }_1/{\beta }_1}{R}_c^{{\alpha }_2-{\beta }_2/{\beta }_1}\end{array}$$(D.3)

where C₃ = C′₁C′₂^{−α₁/β₁}. Inserting the values α₁ = −0.49, α₂ = −1.03, β₁ = −0.85 and β₂ = 0.83 obtained from a χ²− fit to the corresponding values tabulated in Tables 3 and 4, we get ${\nu }_{\mathit{QPO}}=C_3{L}_{\mathit{bol}}^{0.57}{R}_c^{-1.52}$ where C₃ = 5.10⁻⁸. Therefore, Eq. D can be written

$$\begin{array}{c}\hfill {\nu }_{\mathit{QPO}}=5.{10}^{-8}{L}_{\mathit{bol}}^{0.58}{R}_c^{-1.52}Hz=1.1L_{38}^{0.58}{R}_9^{-1.52}Hz,\end{array}$$(D.4)

where L₃₈ = L_bol/10³⁸ erg/sec and R₉ = R_c/10⁹ cm. This should be compared with Eq. (17) of MPK where they find ${\nu }_{\mathit{QPO}}^{\mathit{MPK}}\propto {\dot{m}}_{\mathit{MPK}}^{0.5}{R}_c^{-1.5}$. Making the correspondence between L_bol and Ṁ_MPK, we find that the numerical factors in the two expressions differ by a factor of ≃3. We note, however, that for this estimation we have used only one curve of Fig. 7. Using other curves would change slightly this result.

A difference in the results of the two papers lies on the type of oscillations found. In MPK the oscillations were always damped, while in the present paper we found that in some cases the oscillations have negligible damping. Both papers solve essentially Lotka-Volterra type of equations, albeit with a subtle difference which can be attributed to the implicit assumptions made about electron energy gains. In MPK the energy gains depend on ṁ_MPK, and therefore the electron energy density grows linearly with time (see their Eq. 2); in the present approach the energization term used causes the electron energy density to increase exponentially (see Eqs 2 and 5 of the present paper as well as Appendix A). Since the Lotka-Volterra equations initially produce an exponential growth of the preys, the present paper more closely resembles these. However, here we do not solve simply for two equations as in MPK, but we use a system of many equations for prey (electrons) and for predators (photons) which complicates the overall picture. These complications can be seen in Appendix A which shows, for example, that damping starts from photon energies and move towards the lower ones with time.

Appendix E: Energy dependent acceleration and escape

In the present paper we have adopted an energization and escape timescale that are both independent of energy. However, in most acceleration scenarios, such as Fermi-I type, one would expect these timescales to increase with electron energy. In this Appendix, therefore, we assume that both energization and escape timescales depend on energy and repeat the procedure outlined in the main body of the paper. The formal solution of the electron equation can be found in Kirk et al. (1998) where it was shown that the electron distribution function n_e ∝ γ^{−1 − t_en, 0/t_esc, 0} where t_en(γ) = t_en, 0γ and t_esc(γ) = t_esc, 0γ, therefore the electron distribution function is flatter than the energy-independent case studied previously. Table E.1 shows results using the same parameters of Tab. 3. All output parameters are comparable, i.e. both ν_QPO and luminosities differ by less than a factor of 2-3, therefore a more realistic choice of the energization timescale does not change the overall results significantly. We note that the luminosities increase in the present case because the electron distribution function becomes flatter, therefore there is even more energy in the upper cutoff.

Appendix F: Implications for Active Galactic Nuclei

The present work has been concentrated around QPOs in BHXRBs. However, its scale-free results can be readily applied to Active Galactic Nuclei (AGN) as well. Following the arguments given in Sect. 5.3 and 6, one can scale all results obtained in the present work by the following rules

• Scale R_c to a new value R′_c and let κ = R′_c/R_c.

• Use a new value of $\dot{M}{\prime }_{\mathit{inj}}=\kappa {\dot{M}}_{\mathit{inj}}$.

• Keep all other input parameters (${\stackrel{\sim }{t}}_{\mathit{en}}$, ${\stackrel{\sim }{t}}_{\mathit{esc}}$, α, T_eff) unchanged.

Then the output parameters will scale as follows

• The QPO frequency will scale as $\nu {\prime }_{\mathit{QPO}}={\kappa }^{-1}{\nu }_{\mathit{QPO}}$.

• Both L_bol and L_X will scale as L′=κL.

• All other output parameters (${\stackrel{\sim }{T}}_{10}$, f_rms, ΔT_sh/t_cr, E_pk, Γ and η will remain unchanged.

Note that the above show that the results scale with the size of the corona R_c and not with the black hole mass. Based on the above and on the results shown in Tab. 3 we can construct Table F.1 for three R_c values that can characterize AGN corona.

All Tables

All Figures

	Fig. 1. Schematic representation of the feedback loop between electrons and photons in the coronal environment. For artistic reasons the corona is depicted as a cloud, while in the paper its shape is assumed to be spherical. (The figure is AI generated).
In the text

Fig. 2. Light curves of ICS compactness in the 2–200 keV band for various values of the parameter ℓds. All other parameters have been kept constant and have the values $R_c={10}^9\mathrm{cm},{\dot{M}}_{\mathrm{inj}}=3.2\times {10}^{-7}{M}_{\odot }/\mathrm{yr},{\stackrel{\sim }{t}}_{\mathrm{en}}=2.5,t_{\mathrm{esc}}/t_{\mathrm{en}}={10}^3,T_{\mathrm{eff}}=4\times {10}^6\mathrm{K},\alpha =0.1$. The colors of the lines correspond to ℓds = 10−6,  10−5,  10−4,  10−3,  10−2, and 10−1 from red to violet. For the sake of clarity each curve is displaced from the previous one by unity.

In the text

	Fig. 3. Average photon spectra for the same cases shown in the previous figure. The color codes have been kept the same.
In the text

Fig. 4. Set of light curves all sharing the same value of Qe, inj (Minj = 16 × 10−8 M⊙/yr) but having different energization timescales. Here ${\stackrel{\sim }{t}}_{\mathrm{en}}=1.25$ (violet), 2.5 (dark blue), 5 (light blue), 10 (green), 20 (yellow), 40 (orange), and 80 (red), all given in units of tcr. The other parameters are Rc = 109 cm,  tesc/ten = 103,  α = 0.1,  Teff = 4 × 106 K,  ℓds = 10−6. For the sake of clarity each curve is displaced from the previous one by 1.

In the text

	Fig. 5. Average multiwavelength spectrum for the case shown in the previous figure, keeping the color codes the same.
In the text

	Fig. 6. Heat map of the parameter ${\stackrel{\sim }{T}}_{10}$ as a function of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁinj for Rc = 109 cm. Other parameters are α = 0.1, ℓds = 10−6 and Teff = 4 × 106 K. Dark regions correspond to light curves with strong decay.
In the text

Fig. 7. QPO frequency νQPO versus the energization timescale ${\stackrel{\sim }{t}}_{\mathrm{en}}$ for various values of Ṁinj. Values of Ṁinj, in units of 10−8 M⊙/yr, are 64 (violet), 32 (dark blue), 16 (light blue), 8 (green), 4 (yellow), 2 (orange), and 1 (red). Other parameters are Rc = 109 cm, α = 0.1, and Teff = 4 × 106 K. The black line indicates the relation ${\nu }_{\mathrm{QPO}}\propto t_{\mathrm{en}}^{-1/2}$ given by Eq. (21).

In the text

	Fig. 8. QPO frequency νQPO versus the 2 − 200 keV luminosity LX for various values of Ṁinj. Values of Ṁinj, in units of 10−8 M⊙/yr, are 32 (violet), 16 (dark blue), 8 (light blue), 4 (green), 2 (yellow), 1 (orange), and 0.5. Other parameters are Rc = 109 cm, α = 0.1, and Teff = 4 × 106 K.
In the text

	Fig. 9. Heat map of the average peak energy of the photon distribution Epk as a function of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁinj for Rc = 109 cm. Other parameters are α = 0.1 and Teff = 4 × 106 K. Dark blue regions indicate regimes with low values of Epk.
In the text

	Fig. 10. Heat map of the time lag ΔTsh (in units of tcr) as a function of ${\stackrel{\sim }{t}}_{\mathrm{en}}$ and Ṁinj for Rc = 109 cm. Other parameters are α = 0.1 and Teff = 4 × 106 K. Violet color indicates negative lags, i.e., hard photons come before the soft ones in the prescribed energy bands.
In the text

Fig. 11. Phase lag Δϕ versus the QPO frequency νQPO. The violet line corresponds to Ṁinj = 8 ×10−8 M⊙/yr for Rc between 109.3 cm and Rc = 108.55 cm, the dark blue line corresponds to Ṁinj = 16 ×10−8 M⊙/yr for Rc between 109.6 cm and 108.55 cm, while the light blue line corresponds to Ṁinj = 32 ×10−8 M⊙/yr for Rc between 109.9 cm to 108.7 cm. In all cases the radius Rc is decreasing from left to right. The other parameters are ${\stackrel{\sim }{t}}_{\mathrm{en}}=2.5$, α = 0.4, and Teff = 4 × 106 K. The red lines show the results of varying ${\stackrel{\sim }{t}}_{\mathrm{en}}$ by a factor of 2 (i.e., between .25 and 5) when Ṁinj = 16 ×10−8 M⊙/yr and the corona radius is Rc = 109.3 cm,  109 cm, and 108.7 cm (l to r). As ${\stackrel{\sim }{t}}_{\mathrm{en}}$ decreases, the curves move up and to the right. Orange lines repeat the same procedure with Teff varying by a factor of 2 (i.e., between 2 and 8 × 106 K). As Teff increases, the curves move upwards. The data points are taken from Zhang et al. (2020).

In the text

	Fig. 12. Fractional RMS versus the QPO frequency νQPO for the cases depicted in Fig. 11. The fractional RMS increases with decreasing ${\stackrel{\sim }{t}}_{\mathrm{en}}$ (red lines) and increasing Teff (orange lines). The data points are taken from Zhang et al. (2020).
In the text

	Fig. A.1. Generic light-curve showing the first two peaks of the bolometric ICS compactness. The red and green lines indicate the start of the phases where electron cooling and re-energization prevails respectively; for details see text.
In the text

Fig. A.2. Snapshots of a generic case showing the electron distribution ne(γ) multiplied by γ2 (left panel) and the emitted photon spectra in x.dℓ/dx units (right panel) taken every 2tcr during the initial particle energization phase. Here x = ϵ/mec2. The photon light curve of this case is given in Fig. A.1 from t = 0 up to the first red vertical line. Red lines in the present figure indicate the peak of the electron instantaneous distribution (left) and of the ICS photons (right). As time passes both peaks move from left to right.

In the text

Fig. A.3. Same quantities as in Fig. A.2 for the electron cooling phase, i.e. for the time interval between the red and green lines of Fig. A.1. Snapshots are taken at every tcr. The violet lines indicate the distribution functions at the peak of the photon light curve. Light blue lines depict the distribution functions for times before the peak, while the dark blue lines for times after the peak. As time passes, the peaks of both distributions move from right to left.

In the text

Fig. A.4. Same quantities as in Fig. A.2 for the electron re-energization phase, i.e. for the time interval defined between the green and the second red line of Fig. A.1. Snapshots are taken at every 2tcr. Light blue lines depict the photon distribution function for times before the photon minimum, while the dark blue lines for times after the minimum. As time passes, the peaks of both distributions move from left to right.

In the text

Fig. A.5. Electron density function ne multiplied by γ2 for various electron energies as a function of normalized time for the run with the parameters that produced the light curve shown in Fig. A.1. Red line corresponds to γ = 100.3, orange to 100.6, yellow to 100.9, green to 101.2, light blue to 101.5, dark blue to 101.8 and violet to 102.1. The bolometric photon light curve of Fig. A.1 is shown with black line. The right panel shows a zoom of the first peak for clarity.

In the text

	Fig. B.1. Photon spectrum x.dℓ/dx for various photon energies as a function of normalized time around the first peak of the light curve shown in Fig. A.1. The photon energies correspond approximately to the characteristic radiation of the electrons through the ICS relation in the Thomson limit xi = γi2x0. For color codes see Fig. A.5.
In the text

Fig. B.2. Time lags ΔTi of monochromatic photon energy peaks with respect to the time where the black body peak occurs (red line). The average multiwavelength spectrum for the same cycle has also been plotted with green line. The vertical black lines indicate the "soft" (2-6 keV) and "hard" (6-15 keV) energy bands sometimes used in the literature to measure phase lags. The horizontal black lines indicate the average lags in these two bands. In this case the relative soft-hard lag is positive.

In the text

	Fig. B.3. Time lags ΔT of monochromatic photon energy peaks with respect to the time where the peak of the reflected component occured. The violet line depicts a negative lag, while the orange one a positive. The vertical black lines indicate the "soft" (2-6 keV) and "hard" (6-15 keV) energy bands. The horizontal colored lines indicate the average lags in these two bands.
In the text

	Fig. B.4. Monoenergetic time lags (red line) and the multi-wavelength spectrum (green line) in the case where the reprocessed photons are a delta-function.
In the text

Fig. B.5. Time lags ΔT of monochromatic photon energy peaks with respect to the time where the peak of the reflected component occurred. The temperatures of this component were assumed to be 106K (red), 2 × 106K (orange), 4 × 106K (green) and 8 × 106K (violet) and correspond to the cases shown in Tab. 7 of the main text. The vertical black lines indicate the "soft" (2-6 keV) and "hard" (6-15 keV) energy bands.

In the text