A variable corona for GRS 1915+105

Most models of the low frequency quasi periodic oscillations (QPOs) in low-mass X-ray binaries (LMXBs) explain the dynamical properties of those QPOs. On the other hand, in recent years reverberation models that assume a lamp-post geometry have been successfull in explaining the energy-dependent time lags of the broad-band noise component in stellar mass black-holes and active galactic nuclei. We have recently shown that Comptonisation can explain the spectral-timing properties of the kilo-hertz (kHz) QPOs observed in neutron star (NS) LMXBs. It is therefore worth exploring whether the same family of models would be as successful in explaining the low-frequency QPOs. In this work, we use a Comptonisation model to study the frequency dependence of the phase lags of the type-C QPO in the BH LMXB GRS 1915+105. The phase lags of the QPO in GRS 1915+105 make a transition from hard to soft at a QPO frequency of around 1.8 Hz. Our model shows that at high QPO frequencies a large corona of ~ 100-150 R_g covers most of the accretion disc and makes it 100% feedback dominated, thus producing soft lags. As the observed QPO frequency decreases, the corona gradually shrinks down to around 3-17 R_g, and at 1.8 Hz feedback onto the disc becomes inefficient leading to hard lags. We discuss how changes in the accretion geometry affect the timing properties of the type-C QPO.


INTRODUCTION
Low frequency (LF) Quasi-Periodic Oscillations (QPOs) in stellarmass black hole binaries (BHB) have been known for many years (see reviews by Motta 2016 andMotta 2020). These QPOs are distinct peaks in the Power Density Spectra (PDS) of X-ray light curves of these sources. Based on the strength of the underlying broad-band variability, centroid frequency, , and quality factor, = /Δ , where Δ is the full width at half maximum around the centroid frequency, LF QPOs are divided into three types, type-A, -B and -C (Casella et al. 2004). While all of these types share similar centroid frequencies, type-C are the most frequent QPOs, with high -factors. Although type-C QPOs usually appear in the range of a few mHz to 10 Hz, they have also been detected at frequencies as high as 30 Hz (Revnivtsev et al. 2000). Models that explain the dynamical origin of QPOs in BHBs were proposed about two decades ago when Stella & Vietri (1998) introduced the idea of relativistic Lense-Thirring precession (LTP) as the underlying mechanism. Follow-up theoretical work by Stella et al. (1999), Psaltis & Norman (2000), Fragile et al. (2001), ★ E-mail: karpouzas@astro.rug.nl Schnittman et al. (2006) and Ingram et al. (2009) connected the LF QPOs in BHBs to the LTP frequency. Since then, there have been a plethora of studies that used timing and spectral features to add support to LTP as the origin of LF QPOs (Motta et al. 2014a; Motta et al. 2014b;Ingram et al. 2016 andIngram et al. 2017). Other models, implementing different physical mechanisms were being developed almost at the same time; some of those are the Two-Component Advection Flow model (TCAF, Molteni et al. 1996;Chakrabarti et al. 2008) and the Accretion-ejection instability model (Tagger & Pellat 1999). For recent reviews of observations and theory of LF QPOs, we refer the reader to Done et al. (2007), Motta (2016) and Ingram & Motta (2020).
The picture of the geometry of accretion onto a BH presented in Ingram et al. (2009) consists of two main components. The first component is the classical optically thick geometrically thin accretion disc (Shakura & Sunyaev 1973), which is thought to be truncated at a radius from the BH, typically larger than the inner-most stable circular orbit (ISCO), , depending on the spectral state (Ichimaru 1977;Esin et al. 1997;Poutanen et al. 1997;Gierliński et al. 2008;). The second component is a geometrically thick hot flow that is assumed to reside between the truncated disc and the ISCO, and to be misaligned with respect to the ac-cretion disc plane. General relativistic hydrodynamical simulations suggest that, under certain conditions, such thick inner flows can occur (Fragile & Anninos 2005;Fragile et al. 2007). The core idea of Ingram et al. (2009) is that the type-C QPO is produced by such torus-like inner flow which precesses at the LTP frequency. The outer radius of the torus, , is smaller or equal to the truncation radius of the accretion disc, while the inner radius is set by the limit at which the surface density of the torus becomes sufficiently low, due to the high tilt angle of its inner parts (see Ingram et al. 2009 for a detailed description). The idea of a precessing torus, rather than a precessing disc, is partly motivated by the fact that the energy spectrum that is modulated at the frequency of the type-C QPO (covariance or rms spectrum), is well described by a Comptonised spectrum (Sobolewska & Życki 2006).
The general spectral behaviour of BHBs is usually described by a transition from a hard, less luminous, to a soft luminous state, imprinted as a q-shaped trajectory on the Hardness-Intensity Diagram (HID; Belloni et al. 2005 and, with regimes of intermediate hardness and luminosity. The frequency of the type-C QPO increases as the source moves from harder to softer states. Depending on the model, the change in hardness in the HID can be interpreted as a change of the outer radius of a hot inner flow (Ingram et al. 2009) or, for models with an extended Comptonising medium (Kazanas et al. 1997), as a change in properties of this Comptonising medium such as size, optical depth and temperature. Hereafter, we will refer to this extended Comptonising medium as the corona (Thorne & Price 1975;Sunyaev & Truemper 1979).
The timing properties of the type-C QPO, namely the fractional rms amplitude in a broad energy band and phase lag between a hard and a soft band, have also been shown to depend upon QPO frequency and source inclination. In particular, Motta et al. (2015) suggested that the rms amplitude of the type-C QPO as a function of QPO frequency is systematically higher for high than for low-inclination sources. Similarly, van den Eĳnden et al. (2017) provided evidence that above a certain QPO frequency in low-inclination sources the phase lags of the type-C QPOs are positive (hard lags), meaning that high-energy photons arrive at the observer after the low-energy photons, and these hard lags increase with QPO frequency, whereas for high-inclination sources the phase lags are negative (soft lags) and decrease with QPO frequency. The aforementioned results favour a geometric origin for the type-C QPOs, which was further supported by more quantitative results, such as the modulation at half the QPO frequency of the centroid energy of the iron emission line in H1743 − 322 (Ingram et al. 2016). Many of the observational findings mentioned above still remain unexplained.
Over the years, considerable work has been done to study the dependence of the phase lags between two broad energy bands not only upon QPO frequency, as we mentioned in the previous paragraph, but also upon the frequency of the broad-band noise. At low Fourier frequencies of the broad-band noise, the lags are hard (Miyamoto et al. 1988;Kotov et al. 2001), and thought to be produced by fluctuations of the mass accretion rate that propagate from the outer to the inner part of the accretion disc (Lyubarskii 1997;Arévalo & Uttley 2006;Ingram & van der Klis 2013), and are imprinted on the Comptonised and un-Comptonised emission that is used to calculate these lags. At higher Fourier frequencies of the broad-band noise the magnitude of the lags decreases, and the lags usually become soft, providing strong evidence for reverberation (Uttley et al. 2011;De Marco et al. 2015). Strong reverberation signatures were also recently reported in the broad-band noise component of MAXI J1820 + 070 by Kara et al. (2019), who studied the soft phase lags of the broad-band noise both in the energy and frequency domain. These authors suggest that a compact corona with a height of about 5 , where = / 2 is the gravitational radius, and , and are the gravitational constant, BH mass and speed of light, respectively, illuminates an accretion disc that is truncated at ∼ 2 , thus producing the measured soft lags. As the corona contracts, the soft lags become shorter due to the shorter light travel-time. Further attention to reverberation was given, since an accurate reverberation mapping model was developed by Ingram et al. (2019) that is able to put constraints on the mass of the BH, as was shown by Mastroserio et al. (2019).
Undoubtedly, reverberation is an essential physical mechanism for modeling spectra of BHBs at certain states. However, the treatment of the corona as a point-like source at a fixed height above the accretion disc in reverberation models limits the prospects of a self-consistent explanation of the observed phenomenology. In fact, several previous studies (Miller et al. 2010;Legg et al. 2012;Mizumoto et al. 2018;Mizumoto et al. 2019 have presented the timing analyses of Active Galactic Nuclei (AGNs) where, in general, the phenomenology is similar and requires the presence of an extended corona that partially covers the accretion disc. Moreover, the soft lags of the LF QPOs in stellar mass BHBs, on which we focus from now on, are usually much larger than the lags of the broad-band noise component (Wĳnands et al. 1999), and therefore a compact corona or hot flow with a small scale height is unable to produce such large soft lags. Specifically for stellar mass BHBs, Mizumoto et al. (2016) showed that in the case of GRS 1915+105 the lack of a variation of the iron line at different timescales further supports the idea of an extended corona. In addition, De Marco & Ponti (2016) showed that in order to explain the large soft lags of H 1743−322 during the hard states, when the type-C QPO is present, with reverberation, one has to assume an extended corona as an illuminating source, with a height of a few hundred . We note also that the broad-band noise component in NS and BHC sources consists of a combination of QPO-like components (Nowak 2000;Belloni et al. 2002) and therefore it is natural to study the rms and lag spectra of these QPO-like components to understand the properties of the boradband noise.
Although the soft lags of LF QPOs have been extensively studied in the literature, little work has been done to explain quantitatively their energy or frequency dependence in the context of a radiative mechanism. Nobili et al. (2000) suggested a Comptonisation model that consists of a two-component corona, which produces both hard and soft lags through Compton up-and down-scattering, respectively, in these two different components. Their model successfully explains both the magnitude and the change of sign of the lags as a function of QPO frequency in GRS 1915+105.
Recently, Karpouzas et al. (2020) presented a model that explains the observed soft lags of the lower kHz QPO in neutron star (NS) LMXB, as a delayed heating of the seed photon source by photons previously up-scattered in the corona. This effect is referred to as feedback. In the case of NS LMXBs, it has been shown through Monte Carlo simulations (Kumar & Misra 2016) that feedback should play an important role even among different corona geometries, and that feedback depends on the size of the corona. Although the geometry of the corona and seed photon source are very different between NS and BH LMXBs, one can in principle test how well feedback of up-scattered photons in an extended corona can explain the soft lags observed in BH LMXBs. In fact, based on the model of Karpouzas et al. (2020), García et al. (2021) showed that a two-component extended corona explains the energy dependence of the time lags and rms amplitude of the type-B QPO in MAXI J1348−630, and so in this work we extend our analysis for type-C QPOs. A good candidate to perform such a test is GRS 1915+105 (Castro-Tirado et al. 1992Castro-Tirado et al. 1994). GRS 1915+105 is one of the most observed, and thus best studied, sources in the Rossi X-ray timing explorer (Bradt et al. 1993, RXTE) archive. The compact object in GRS 1915+105, is a BH with a mass = 12.4 +2 −1.8 ⊙ , (Reid et al. 2014) and spin estimated, with spectral techniques, to be between * = 0.68 and * = 0.99 (Šrámková et al. 2015). GRS 1915 + 105 is a peculiar source, in the sense that it never showed a full Q-shaped cycle on the HID. For a dedicated review on GRS 1915+105, we refer the reader to Fender et al. (2004). Its spectral behaviour was separated in 14 classes based on the work of Belloni et al. (2000) (see also Klein-Wolt et al. 2002 andHuppenkothen et al. 2017). Timing studies of GRS 1915+105, which are mainly what we are interested in, date back to Chen et al. (1997) and Morgan et al. (1997) and since then, several others focused on the spectral-timing properties of LF QPOs in this source (Markwardt et al. 1999;Vignarca et al. 2003;Rodriguez et al. 2004). One of the most notable timing properties related to the type-C QPO phase lag between two broad bands of GRS 1915+105 is possibly the switch from hard lags, at QPO frequencies below ∼2 Hz (Reig et al. 2000;Pahari et al. 2013;Zhang et al. 2020), to zero at around 2 Hz, and then to soft lags at QPO frequencies above that limit. The slope of the lag-energy spectrum follows the same behaviour, changing from positive below 1.8 Hz to zero at around 1.8 Hz, and then negative above 1.8 Hz. Alongside with the phase lags, the rms amplitude, both energy-and frequencydependent, of the type-C QPO also shows a particular dependence upon QPO frequency (Zhang et al. 2020).
Recent work from Dutta & Chakrabarti (2016) and Chatterjee et al. (2017), implementing the TCAF model (Chakrabarti et al. 2008), linked the QPO frequency to the size of the Centrifugal pressure supported Boundary Layer (CENBOL), which serves as the corona in the TCAF model. These studies suggested that changes in the size of the CENBOL can explain the change in sign of the time lag. More recently, Dutta et al. (2018) applied the TCAF model to data of GRS 1915+105, and showed how a gradual change in the size of the CENBOL can explain the time evolution of the type-C QPO frequency, above and below the 1.8 Hz limit. Dutta et al. (2018) concluded that, as the QPO frequency decreases from over 5 Hz to below 1.8 Hz, the CENBOL increases in size from around 135 to over 500 , and then it decreases again to about 100 , when the QPO frequency increases beyond 5 Hz.
Here, we use the Comptonisation model of Karpouzas et al. (2020) to explain the dependence of both the rms amplitude and phase lags upon QPO frequency and energy. In sections 2 and 3 we, respectively, introduce the model and explain its application to the data of Zhang et al. (2020). In section 4 we show our results. Finally, in section 5 we discuss our results and provide a new view to how the corona evolves over time, and how the latter affects the timing behaviour of GRS 1915+105.

COMPTONISATION MODEL
In this work we use the model presented by Karpouzas et al. (2020), which is an adaptation of the Comptonisation model proposed by Lee & Miller (1998); Lee et al. (2001); and Kumar & Misra (2014). The main idea behind these models is that any QPO can be described as an oscillation of the time averaged spectrum, 0 , that is cou-pled to oscillations of other properties of the system, such as the corona temperature, , seed-photon source temperature, , and the external heating provided to a finite-sized corona, .
No assumption is made about the origin of the QPO frequency in the model. As far as the model is concerned, the QPO is an oscillation of the X-ray flux available for Compton up-scattering in the corona. As a consequence, the properties of the corona and seed-photon source can react to the QPO and oscillate coherently. Another possibility is that the QPO itself is produced by some instability in and or some oscillatory mode in the corona (see Ingram et al. 2009;Fragile et al. 2016;Fragile 2020).
The model of Karpouzas et al. (2020), assumes that the seed photon source is a blackbody and that the Comptonising corona is a spherically symmetric homogeneous shell with optical depth and thickness around that blackbody. The model also takes into account feedback of up-scattered photons onto the blackbody. When applying the model in practice, first we calculate 0 by solving the Kompaneets equation (Kompaneets 1957) in steady-state. Next, we linearise the Kompaneets equation assuming that , and undergo small oscillations at exactly the QPO frequency, and solve for the complex amplitude, , of the averaged spectrum, 0 . We note that the oscillation of the seed-photon source temperature, , is attributed to feedback photons. Feedback, within the model, is regulated by the feedback fraction, , that is defined as the fraction of the luminosity of the seed source that is solely due to feedback. The complex amplitude, , holds information about the energy-dependent amplitude and phase lags at the QPO frequency, , that was used to perform the linearisation. We note here that only 0 and are functions of photon energy, while , , , and are considered constant for a specific QPO frequency. In Karpouzas et al. (2020) we showed that the model can fit the energydependent time lags and fractional rms amplitude of the kHz QPOs of the NS LMXB, 4U 1636−53. In this paper we make a first attempt to quantitatively explain the timing properties of the type-C QPO in the BH LMXB GRS 1915+105 using this same model. Our approach introduces two caveats. Firstly, the source of seed photons that we use is a blackbody instead of a disc-blackbody. Secondly, the photonloss term in the Kompaneets equation assumes an equal optical depth seen by all seed-photons, which is a crude approximation in the case of a disc whose spatial extent is comparable to the size of the corona. We will expand on these issues in the Discussion section. Zhang et al. (2020) presented an extensive timing analysis of 620 observations of GRS 1915+105 showing type-C QPOs, after analysing all of the RXTE archival data (1996−2012). The QPO frequencies ranged from 0.4 Hz to 6.3 Hz. Zhang et al. (2020) measured the phase lags of all type-C QPOs between two broad bands, 2−5.7 keV and 5.7−15 keV. They also measured the lags between multiple narrow energy bands, which we hereafter refer to as lag-energy spectra. In the paper of Zhang et al. (2020), lag-energy spectra at nine selected QPO frequencies were shown for reference. Similarly, Zhang et al. (2020) measured the fractional rms amplitude both in the full PCA band (2−60 keV), and in separate bands. The latter will be referred to as the rms-energy spectrum. The primary focus of Zhang et al. (2020) was to study the frequency dependence of the phase lags between the two broad bands mentioned above, which they called the lag-frequency spectrum. The authors fitted a broken line to the lag-frequency spectrum and placed a suffi-ciently precise limit of 1.8 ± 0.1 Hz on the QPO frequency at which the measured phase lags change from hard to soft. Their best-fitting broken-line to the lag-frequency spectrum also exhibited a statistically significant difference in slope, with a value of −0.21 ± 0.02 below 1.8 Hz and −0.1 ± 0.01 above 1.8 Hz, indicating a potential change in the physical mechanism that produces the lags. Zhang et al. (2020) showed that the frequency dependence of the rms amplitude in the full PCA band, which we will refer to as the rms-frequency spectrum, also changes behaviour around the same QPO frequency, switching from increasing with QPO frequency for QPO frequencies below 1.8 Hz, to decreasing for QPO frequencies above 1.8 Hz. Although the model of Karpouzas et al. (2020) can predict both energy and frequency dependence of the fractional rms amplitude and phase lags, in this work we first study the frequency-dependent fractional rms amplitude in the full band, and phase lags between the two broad bands used in Zhang et al. (2020). The latter approach allows us to analyze all 620 QPO frequencies, while avoiding the computationally expensive process of fitting the rms-and lag-energy spectra at all QPO frequencies. By studying the frequency dependence, we obtain an initial estimate of how the model parameters depend upon QPO frequency. Then we proceed with fitting the rms-and lag-energy spectra presented by Zhang et al. (2020), at representative QPO frequencies, and compare with our initial estimate for completeness. The challenge of the frequency-dependent analysis, with which we start, is that we only have two measurements at our disposal, namely the fractional rms amplitude and phase lag for each set of the 7 parameters of the model, , , , , , that are needed to produce a pair of simulated fractional rms amplitude and phase lag at a certain QPO frequency, . Thus, one can not perform a fitting in the classical sense. In the following, we describe how we can overcome the aforementioned limitation by matching the predicted, from the model, frequency dependence of the fractional rms amplitude and phase lags to the measured ones. Hereafter, we will refer to the term fractional rms amplitude simply as rms amplitude.

DATA ANALYSIS AND MODELING
As shown by Karpouzas et al. (2020), for a given set of the 7 parameters mentioned above, the model simultaneously predicts the rms amplitude in any band and phase lag between any two bands. However, the amplitude of the external heating rate, , only affects the rms amplitude, acting as a normalisation, while the phase lags are not dependent on this parameter. Therefore, a good approach is to only simulate the lag-frequency spectrum at first, and then re-scale the simulated rms-frequncy spectrum, that was anyway produced alongside the lag-frequency spectrum, by adjusting in order to match the measured rms amplitude values. In the following sections we discuss the modeling of the phase lags first, and then we study the resulting rms amplitude.

Modeling of the lag-frequency spectrum
To explain the change of sign of the phase lag around 1.8 Hz and its dependence upon QPO frequency, we designed a simple computational experiment using the model of Karpouzas et al. (2020). In particular, we generated 10 6 random combinations of the model parameters uniformly for in the range 3−100 keV, in the range 0.1−3 keV, in the range 2−500 , in the range 0−1 and in the range 0.4−6.3 Hz. We then calculated using and a random uniformly distributed value of Γ between 1.5 and 3, to cover a wide enough range of the power-law index. To calculate , for every selected value of and Γ, we used the formula: where is the electron mass. We set the external heating rate, , to a constant arbitrary value of 1%, since it only affects the measured rms amplitude and, as has been shown by Karpouzas et al. (2020), it has a dependence upon QPO frequency, which we do not know a priori. Later, when we compare with the measured rms amplitude, we can re-scale to match the data, and reveal its QPO frequency dependence and actual value.
For each randomly generated QPO frequency, there is a corresponding randomly generated combination of the five model parameters , , , and , that when given as input to the model generates a model estimate of the phase lag between the bands 2−5.7 keV and 5.7−15 keV, at that QPO frequency. The energy bands were chosen such that we maintain consistency with Zhang et al. (2020). If one plots all the values of the model estimated phase lags as a function of the randomly generated QPO frequency, as we have done with gray circles in Figure 1, no obvious correlation between the two is observed. However, the data suggests that there exists a significant correlation between phase lag and QPO frequency. To accommodate for the observed correlation, we created two subsamples from our total simulated sample of phase lags, one with pairs of simulated phase lags and QPO frequencies for which the phase lags where positive, at QPO frequencies below 1.8 Hz and clustered around the best-fitting broken line of Zhang et al. (2020) below 1.8 Hz, and another with pairs of simulated phase lags and QPO frequencies for which the phase lags where negative, at QPO frequencies above 1.8 Hz and clustered around the best-fitting broken line above 1.8 Hz. As a condition for the clustering around the best-fitting broken line, we demanded that the predicted phase lag from the model can deviate from the best-fitting broken line as much as one standard deviation of the data. The aforementioned computational experiment is illustrated in Figure 1, alongside with the original data and best-fitting broken line of Zhang et al. (2020).
For each pair of simulated phase lag and QPO frequency, that we restricted alongside the best-fitting broken line in Figure 1, there is also a pair of the five remaining model parameters , , , and . Our goal was to check whether the aforementioned model parameters, that were used to generate the clustered simulated points around the broken line, have any correlation with QPO frequency that separates them from the total simulated population which is uniformly sampled. We will discuss the dependence of the model parameters upon QPO frequency in the Results section of this paper, where they are presented, but for now we will discuss the other by-product of our model, the rms amplitude, which at this point is already calculated but not compared with the measured rms amplitude.

The rms-frequency spectrum
As we mentioned above, our model can predict the rms amplitude at any QPO frequency, simultaneously with the phase lag, for any combination of the model parameters. However, we remind the reader that to compare with measured rms amplitude we need to take into account the amplitude of the external heating rate , . Up to this point, we have clustered the simulated pairs of phase lags and frequencies alongside the best-fitting broken line ( Figure  Figure 1. The measured lag-frequency spectrum of GRS 1915+105 (black points with error bars) alongside with all of our simulations. Grey circles represent all of the simulated pairs of phase lags and QPO frequencies.
Red represents the simulated pairs that cluster within 1-of the best-fitting broken line (blue dashed line) to the measured phase lags. The break, where the phase lags switch from hard to soft, is denoted by the cyan dotted line at ∼1.8 Hz. The white circles that track the best-fitting broken line show individual observations, at selected QPO frequencies, for which the energy dependence of the phase lags was measured (see section 3.3). 1). In Figure 2, we plot the simulated rms-frequency spectrum, i.e. the pairs of simulated rms amplitude and frequency, only for the points for which the corresponding simulated phase lags are clustered around the broken-line of Figure 1. On top of that, we plot the rms amplitude measurements of Zhang et al. (2020) alongside with their uncertainties. In Figure 2, we re-scaled the simulated rmsfrequency spectrum by multiplying all of the values by a factor of 11, which is equivalent to assuming an external heating rate amplitude of = 11% instead of 1% that we used before, equal at all QPO frequencies. This re-scaling is totally arbitrary and only serves as a way to visually compare the simulated rms-frequency spectrum to the measured values of Zhang et al. (2020).
One can immediately distinguish two clusters of simulated points in Figure 2. These clusters reflect the degeneracy caused by the lack of prior knowledge of the spectral parameters of the source in these observations (Karpouzas et al. 2020), , and . To detect and separate the clusters, we used the DBSCAN algorithm (Ester et al. 1996). One of the clusters (blue rectangle in Figure 2) corresponds to very low values of the rms amplitude (< 1%) and does not populate the higher QPO frequencies (>3 Hz) very well, even when arbitrarily re-scaled. The other cluster (red rectangle in Figure 2) appears at higher rms amplitude but does not fully populate the lower QPO frequencies (<2 Hz) sufficiently well. Interestingly, by forcing the simulated lag-frequency spectrum to follow the broken line in Figure 1, the corresponding simulated rmsfrequency spectrum for the upper cluster (red rectangle in Figure  2), has a significant negative correlation, exactly above the 1.8 Hz level, as observed in the measured rms amplitude.
Since is QPO frequency-dependent, in general, the sim-  Figure 1 (contours). The contours were produced by smoothing the two-dimensional space of the simulated points with a Gaussian kernel. The rms amplitude is re-scaled by a factor of 11% to better compare with the data. The red and blue rectangles showcase the separation of the two prominent clusters in the simulations.
ulated rms-frequency spectrum can not perfectly match the data yet, for any arbitrary re-scaling factor that is constant at each frequency. However, if one matches the measured rms-frequency spectrum with the simulated one, the expected frequency dependence of can be extracted. To perform the aforementioned matching, first we divided the simulated frequency range in 15 bins. In each bin we calculated the mode of the distribution of the simulated rms amplitude values and used that as the expected value of the simulated rms amplitude in that bin. We also calculated the 1 − confidence interval around the mode, and used that value as the uncertainty of the simulated rms amplitude in that bin. After doing that, we have an rms-frequency relation which we can match to the measured rms-frequency relation. To do the latter, we simply found a normalisation, different at each QPO frequency that, when multiplied with the simulated rms amplitude, re-produces the measured rms amplitude at each QPO frequency. The value of the normalisation factor at each QPO frequency represents the required external heating rate fractional amplitude, | |, that best fits the data. We plot the result of this matching in the top panel of Figure 3, while the lower panel shows the values of the normalisation factor required at every QPO frequency.
We must note that we discarded the lower cluster (blue rectangle in Figure 2) of the simulated rms-frequency spectrum because it would require an extreme variability of the external heating rate ( ∼ 100%), to agree with the measured rms amplitude. In Figure 2, the contour levels correspond to a smoothing of the simulated rms-frequency pairs in the two-dimensional space using a Gaussian kernel. Seemingly, the contours in Figure 2 show no simulated rms-frequency pairs below 1 Hz. This happens only because the model, under the priors assumed, predicts too few points below 1 Hz, compared to the ones above that frequency, that they appear insignificant after the smoothing process, and thus not visible. This effect will be discussed below.
To study the physical parameters of the model and their dependence upon QPO frequency, which is discussed in the next section, we separated the simulated QPO frequency domain in two regimes, one below 1.8 Hz and the other above 1.8 Hz. For the regime above 1.8 Hz, we used all the simulated models that belonged to the upper cluster of the rms-frequency spectrum (red rectangle of Figure 2). For the regime below 1.8 Hz, we used all the available simulated models due to the fact that the simulations in this regime were more scarce. The issue of scarce model sampling below 1.8 Hz will be addressed in Section 5.3. Zhang et al. (2020) presented rms-and lag-energy spectra of GRS 1915+105 for selected type-C QPO frequencies in the complete range 0.4−6.3 Hz. Those selected frequencies are plotted as white circles in Figure 1, and can be viewed as snapshots of the timing properties of GRS 1915+105 at a given QPO frequency. They showed that, on average, the initially negative slope of the lag-energy spectra systematically increases following a linear trend, while the QPO frequency decreases down to ∼ 2 Hz. Below 2 Hz, the slope of the lag-energy spectrum, on average, crosses the zero level (switching to hard lags) and continues to increase following a significantly different linear trend than the one formed above the 2 Hz limit. Up to this point, our analysis has omitted any information about the energy dependence of the rms amplitude and phase lags provided by the model. Nevertheless, we arrived at an initial estimate of how the physical parameters depend upon QPO frequency. The dependence of the physical parameters upon QPO frequency are shown as black circles in Figure 4.

MCMC fitting to the rms-and lag-energy spectra
Using this initial estimate as a starting point, we performed an individual Markov-chain Monte Carlo (MCMC) fitting to each rmsand lag-energy spectrum selected by Zhang et al. (2020). For the MCMC fitting and analysis we used the affine-invariant ensemble sampler (Foreman-Mackey et al. 2013) and followed the same approach as in Karpouzas et al. (2020). We must note that in the case of the rms-energy spectrum, for reasons that we discuss in Section 5.3, we neglected the rms amplitude measurements above 25 keV. The best-fitting model, alongside with the data of Zhang et al. (2020), at each QPO frequency is plotted separately for the frequencies above 1.8 Hz in Figure 5 and for the frequencies below 1.8 Hz in Figure  6. As expected, the MCMC fitting gives a separate estimate of the dependence of the physical parameters upon QPO frequency. In the next sections we will discuss these results and how they compare to the analysis that does not take the energy dependence of the rms amplitude and phase lags into account.

RESULTS
In section 3.1, we discussed the clustering of the simulated phase lags alongside the best-fitting broken-line of Zhang et al. (2020). This clustering provides an initial estimate of the frequency dependence of the model parameters, thus revealing what, in the context of the model, is the physical mechanism that drives the lag-frequency relation. In Figure 4 we plot five of the model parameters, namely the corona size, , feedback fraction, , electron temperature of the corona, , optical depth, , and seed photon source temperature, , as a function of the QPO frequency plus the photon power-law index, Γ, which is not an independent parameter but was used to generate the values of for the simulation (see section 3.1).
To estimate the frequency dependence of the parameters, we divided the simulated frequency domain of the cases that clustered alongside the broken line in 15 bins. In each frequency bin we studied the distribution of each one of the five parameters and used the mode of the distribution as a best estimate, and the 1 − confidence interval around that mode as the uncertainty of that estimate. These parameters, and their corresponding uncertainties at each QPO frequency are summarized in Table 1. In Figure 4 the black circles represent the values of the parameters as a function of QPO frequency derived from the clustering alongside the best-fitting broken-line ( Figure 1). Hereafter, for simplicity, we will refer to the results derived from the clustering around the best-fitting broken-line process as the frequency-only-approach. The red circles in the same Figure represent the best-fitting parameters derived from the MCMC fitting to the individual rms-and lag-energy spectra of Zhang et al. (2020) at selected QPO frequencies (see Karpouzas et al. 2020 for details on the parameter estimation using the MCMC method). We will refer to the process of MCMC fitting to the rms-and lag-energy spectra as the energy-dependent-approach. The physical parameters derived from the energy-dependent-approach are summarized separately in Table 2.
Based on the frequency-only-approach we find that the corona size decreases systematically from 115 at around 6 Hz down to about 7.2 at 2 Hz, and then increases again to almost 160 at 0.9 Hz. Alongside with the decreasing corona size, in Figure  4 we plot the evolution of the disc inner radius assuming that the type-C QPO is the manifestation of the LTP frequency around a BH with mass 12.4 ⊙ and spin, * , between 0.68 and 0.99. At the same time, the feedback fraction, , stays at a maximum value of ∼ 1, which means that the disc is mainly heated up by up-scattered photons, for QPO frequencies above 2 Hz and then drops abruptly to zero at frequencies below 2 Hz, where the measured lags also change sign. For the rest of the parameters, , and , we did not find significant correlations with QPO frequency, except in the case of the corona temperature which exhibits, on average, lower values at frequencies below 2 Hz than above 2 Hz. More specifically, the corona temperature, , remains above 70 keV, on average, at frequencies above 2 Hz, while at 1.7 Hz it becomes ∼ 60 keV and decreases down to about 11 keV at 1.2 Hz. The amplitude of the external heating rate, | |, which is derived through comparison with the measured rms-frequency spectrum, decreases from around 10% at 5.9 Hz, to 20% at 2 Hz. Below 2 Hz, | | becomes extremely high (>100%), a fact that will be discussed later in section 5.3.
In the energy-dependent-approach, we fitted the rms-and lagenergy spectra of the QPO at the individual frequencies selected by Zhang et al. (2020). This allows us to better constrain the physical parameters of the model. This approach is equivalent to taking a snapshot of what the parameters would be at the time when the QPO frequency had a particular value. Based on the results of the energydependent-approach, we find that the corona size stays above 100 when the QPO frequency is above 2.5 Hz, while at around 2 Hz it becomes 8.7 with a lower 1− uncertainty of 6 , indicating that around that frequency the corona size, within the uncertainties, can be as low as the lowest value in our simulations, which was 2 . However, we note here that, based on the MCMC, the best-fitting corona size between 2 Hz and 3.5 Hz does not seem to agree with the smooth decrease in size that we derive from the frequency-only approach. This happens because the observations between 2 Hz and 3.5 Hz, for which we had energy-dependent measurements (white circles in Figure 1), happened to be selected such that the phase lags were very low, compared to other observations in the same frequency range, which naturally leads to larger fitted corona sizes, , since that parameter is mostly affected by the magnitude of the phase lag. The feedback fraction, , decreases with QPO frequency smoothly from almost 100%, at 5 Hz to almost zero at 0.5 Hz. The corona temperature, although badly constrained, is high (>40 keV) at QPO frequencies above 2 Hz, and it becomes significantly lower (<13 keV) at QPO frequencies below 2 Hz. The optical depth of the corona shows a systematic increase from 0.7 at 5 Hz to 6.2 at 0.5 Hz. The seed-photon source temperature also shows a decrease, on average, from values above 1 keV at QPO frequencies above 2 Hz, to values below 0.3 keV at QPO frequencies below 2 Hz. Finally, the amplitude of the external heating rate increases from 17% at 5 Hz to ∼ 40% at 2 Hz, and then it becomes higher than 100% at QPO frequencies below 2 Hz.
The last column of Table 2 summarizes the 2 of the fits to the rms-and lag-energy spectra at each QPO frequency. For the calculation of 2 we combined both the rms amplitude and phase lags, using the same approach as in Karpouzas et al. (2020). The high 2 values are, in most cases, because of the bad fit of the model to the rms-energy spectrum at energies higher than 20 keV. Although, as we mentioned before, the measured rms in energy bands above 25 keV was not used in the MCMC fitting, the rms amplitude, predicted by the model, fails to flatten out at high energies. The non-flattening of the model rms amplitude at high energies was already noticed in the case of NS LMXBs (Karpouzas et al. 2020), and attributed to a lack of a non-oscillating flux component in the model, such as reflection, that would contribute to a reduction of the rms amplitude at higher energies. This explanation is also valid in the case of GRS 1915+105, where reflection signatures appear in the energy spectrum when the type-C QPOs are observed (Misra et al. 2020).

DISCUSSION
We propose a novel explanation for the change of the time lags of the type-C QPO in GRS 1915+105 from soft, at high QPO frequencies, to hard, at low QPO frequencies. We applied a Comptonisation model that incorporates feedback of the hard corona photons onto the soft photon source to a large data set of timing measurements of this source with the RXTE satellite. We used the rms amplitude and phase lag measurements of Zhang et al. (2020) combined with the model of Karpouzas et al. (2020) and conclude that a corona with time-varying size, accompanied by a switch in feedback efficiency, is able to explain the data. In particular, we find that, as the QPO frequency decreases from 5 Hz to 2 Hz, the size of the corona decreases from around 120 to a few . During this contraction, the accretion disc is covered by the corona, the feedback of Compton up-scattered photons onto the disc is very efficient (∼ 100%), and most of the disc flux is due to feedback photons. As the QPO frequency decreases, so does the feedback efficiency because the corona covers less and less area of the disc, as its size decreases. At the same time, assuming that the truncation radius of the disc is anti-correlated to the QPO frequency, the inner radius of the disc increases and at some critical frequency (∼ 1.8 Hz) the disc is no longer covered by the corona. From that point on, the efficiency of feedback decreases as the disc inner radius moves further out and the corona size continues to decrease. In Figure 7, we present a schematic representation for the evolution of the size of the corona in four phases as the type-C QPO frequency decreases from 5 Hz down to 1 Hz. Our fits are systematically worse when the QPO frequency is below 1.8 Hz than when it is above that value. Our results indicate that when the QPO frequency is below 1.8 Hz the size of the corona increases to values much larger than 130 , the feedback efficiency decreases to values <50% and down to zero, and also the power of the source of external heating for the corona oscillates with an amplitude larger than 100%. The aforementioned paradigm of a decreasing corona and feedback efficiency fits well with all available timing measurements of the type-C QPO in the frequency range of 1.8−5 Hz, and explains both the decrease in phase lags between two energy bands, and the smooth decrease in the slope of the lag-energy spectrum. The change in the behaviour of the lags at 1.8 Hz and below indicates either a possible switch in the physical mechanism that produces the lags or a change in the geometry of the Comptonising region that the model does not take into account. Nobili et al. (2000) proposed a model that was able to explain the lag-frequency spectrum of GRS 1915+105. The difference between their approach and ours is in the mechanism that produces the soft lags. In the model of Nobili et al. (2000), the soft lags are caused by Compton down-scattering of soft disc photons that have previously suffered saturated Comptonisation in an inner optically thick component of the corona. In the model of Nobili et al. (2000) the Compton down-scattering happens in an outer cooler component of the corona. Table 1. Physical parameters of the model at different frequencies of the type-C QPO in GRS 1915 + 105. The frequencies in the first column represent the center of the selected bins in which the distribution of each parameter was analysed. The amplitude of the external heating rate (last column) is deduced by comparison of the simulated rms-frequency spectra to the measured one (see Figure 3). The uncertainties of the model parameters were calculated as described in Figure 3 and are omitted when only upper limits for the parameters could be constrained.

Comparison to other models
[   Table 2. Physical parameters of the model as a function of the type-C QPO frequency in GRS 1915+105. The best-fitting model parameters are derived via MCMC fitting to the energy-dependent rms amplitude and phase lags at each frequency, while the corresponding uncertainties where taken as the 1− confidence interval around the mode of the MCMC posterior distributions for each model parameter. The values of Γ in the seventh column are derived using the best-fitting temperature, , and optical depth, , while the uncertainties of Γ are computed through propagation, using the uncertainties of the same parameters.
[  In each panel the black circles show the parameter estimates based on the total QPO frequency sample (frequency-only-approach) while the red circles show the parameter estimates from the MCMC fitting applied to individual selected frequencies (energy-dependent-approach). The red shaded area, in the upper left panel shows the range of the predicted evolution of the disc inner radius assuming that the type-C QPO frequency is produced by LTP and assuming a range for the masses and spins of the BH in GRS 1915+105. The dotted cyan vertical lines, in all panels, denote the 1.8 Hz limit.
In the model of Karpouzas et al. (2020), the soft lags are interpreted as a delayed response of the seed photon source to feedback photons that impinge back on the disc after being Compton upscattered in the corona. The differences in the physical assumptions between the two models are reflected on the correlation between the inferred corona size and disc inner radius. In the model of Nobili et al. (2000), the disc inner radius is positively correlated to the size of the corona, whereas in our model the two quantities are anticorrelated down to the limit of ∼2 Hz and then positively correlated for lower QPO frequencies. Our model, however, explains not only the lag-frequency spectrum, but also both the rms-and lag-energy spectra at different QPO frequencies. Dutta et al. (2018) based on the TCAF model, provide an explanation opposite to ours when it comes to the evolution of the Comptonising cloud at different type-C QPO frequencies. More specifically, although the size of our corona is of the same order of magnitude as the size of the CENBOL in Dutta et al. (2018), in their interpretation there is no systematic or rapid decrease in the CENBOL size at around 1.8 Hz. Since a transition of a type-C QPO frequency around 1.8 Hz typically happens in matter of a few days (Dutta et al. 2018), the predicted decrease in the size of the corona should happen on the same time scale. Regardless of whether the decrease in corona size is systematic or very rapid, our fits to both the lag-energy and lag-frequency spectra agree that the size should be significantly lower around 1.8 Hz. If our explanation is correct, a reduction in the size of the corona, while the temperature remains the same or decreases, would lead to a reduction in the hard X-ray intensity at low QPO frequencies. Stiele & Kong (2018) and Bhargava et al. (2019)  Best-fitting energy-dependent models to the rms-and lag-energy spectra at selected type-C QPO frequencies above 1.8 Hz in GRS 1915+105. In the left and right panels the rms-and lag-energy spectra are plotted, respectively, alongside with the best fitting models. For the different QPO frequencies, the data-points and model-lines are plotted with different styles colors. In the bottom sub-panels of each side we plot the residuals between the data and model normalised by the uncertainties in the data. In panel a, at high QPO frequency, a spherically symmetric corona (blue layered circle) covers the inner accretion disc (thick horizontal line) and feedback is very efficient. In panel b the QPO frequency is lower and, assuming it reflects the Lense-Thirring frequency at the inner edge of the disc, the inner-disc radius increases; the size of the corona has decreased enough so that it is only marginally covering the inner disc, while feedback efficiency has decreased. In panel c, where the QPO frequency is at around 2 Hz, the inner-disc radius is higher, the corona has reached its minimum size and is now inside the disc. Finally, in panel d the corona size increases again as the QPO frequency decreases, and the inner edge of the disc increases further. matter from the corona is not advected onto the black hole then, given appropriate conditions in the magnetic field, matter could be directed towards the jet, in which case enhanced radio emission should be present at low QPO frequencies. In support to our expectations, Yan et al. (2013), found that increased radio activity is present when the frequency of the type-C QPO is low.
Although at first glance our explanation appears to be in conflict with the model of Ingram et al. (2009), this is not the case. Ingram et al. (2009) assumed that the inner flow around an accreting black hole consists of a geometrically thick torus located inside the truncation radius of a non-precessing geometrically thin disc.
As the torus precesses, it illuminates different parts of the disc causing the modulation of the X-ray flux that produces the LF QPO. Ingram et al. (2016) noted, however, that the same pattern would result if the torus was fixed and it was the disc the one that precessed at the LT frequency (Schnittman et al. 2006;Tsang & Butsky 2013). As Ingram et al. (2016) explained, their choice of one geometry over the other was based on the fact that the rms spectrum of the QPO is hard and hence the emission at the QPO frequency could not come from the disc.
As we showed here, this argument is not applicable, and a scenario with a precessing disc and a fixed corona is possible. Indeed, in our model the rms spectrum of the QPO is a consequence of inverse Compton scattering of soft disc photons in the corona (the torus in the scenario of Ingram et al. 2009), such that the high rms amplitude values of the QPO at high energies simply reflect the variability of the soft disc emission at the LT frequency that is scattered in the corona. This, plus the feedback from the corona to the disc, naturally explain both the variability of the iron line discussed in Ingram et al. (2016) and the rms spectrum of the QPO.
Our interpretation of a corona with a different size at different QPO frequencies, agrees with recent evidence presented by Kara et al. (2019) 2019) studied the broad-band noise component through a reverberation model that assumes a point-like corona above the accretion disc, and presented evidence that the height of the corona changes with time. Here, we used the LF QPOs and applied a Comptonisation model that assumes an extended corona, and provided evidence that the size of this corona changes with QPO frequency. Since the QPO frequency changes with time, the corona size also changes with time. The results of Kara et al. (2019) and ours could be due to a similar mechanism that acts throughout the different spectral states and depending on the exact state, and its spectral-timing properties can be detected by using different analysis and modelling techniques.

On the possible connection between the corona and the jet
Persistent and transient BH LMXBs exhibit significant radio emission (Fender 2001). The radio activity is related to the existence of a jet (Band & Grindlay 1986;Levinson & Blandford 1996;Georganopoulos et al. 2002;Reig et al. 2003), that can produce hard X-rays through Compton up-scattering. Jet models, such as the ones presented in Reig et al. (2003) and Kylafis & Reig (2018), successfully explain the hard lags as a function of broad-band frequency in the context of Comptonisation of disc seed-photons in the jet. However, to the best of our knowledge, the soft lags of the type-C QPOs, and particularly the transition from hard to soft lags, has yet to be explained by jet models.
The radio emission in GRS 1915+105 is occasionally anticorrelated to the hard X-ray flux (Trudolyubov 2001). The fact that the radio emission of GRS 1915+105 is stronger at lower type-C QPO frequencies (Yan et al. 2013), combined with the fact that the lags are hard at low QPO frequencies and, as we showed here, Comptonisation in a uniform corona provides systematically worse fits to the data, naturally leads to the following scenario: At high QPO frequencies an extended corona with efficient feedback onto the disc can explain both the power-law component in the energy spectrum and the lags and rms-amplitude of the type-C QPO. As the size of the corona decreases, feedback onto the disc and Comptonisation in the corona are no longer the dominant components that drive the timing-spectral properties, and so the jet takes over as the medium that produces the power-law and the hard lags. This scenario, if validated, can explain all the spectral and timing properties at once, as well as the correlation between X-ray and radio luminosity.

Model Caveats
We modelled the rms-frequency, lag-frequency, rms-energy and lag-energy spectra of the type-C QPO in GRS 1915+105, using the Comptonisation model of Karpouzas et al. (2020). As with every model, there are certain caveats that we address here. Firstly, the model tends to produce less simulated phase-lag and rms-amplitude values at QPO frequencies lower than 1.8 Hz (Figures 1 and 2). This happens because in the model the energy-dependent rms-amplitude and lags do not change as smoothly for low values of the feedback fraction, , as for high values ( >0.3). This behaviour is inherent to the model and depends on the assumptions upon which it is built. To produce hard lags, the model requires less feedback, since in the model increasing the feedback means increasing the time-delay of the soft component with respect to the hard one, and so regimes of hard lags will tend to be poorly sampled. A detailed study of how the model depends upon each physical parameter will be presented in an upcoming publication (García et al. in prep.). Secondly, the model of Karpouzas et al. (2020) uses a simple blackbody as the seed-photon source for Comptonisation. Despite this simplification, it is remarkable how good the model works when applied to the data. While a multi-colour disc-blackbody would be a more realistic assumption, this would not only add two additional free parameters (source inclination and disc inner radius), but it would complicate the definition of both the optical depth, which should then have an extra spatial dependence, and of the feedback, since the feedback efficiency should then depend on the inner disc radius. Addressing the aforementioned problems will be part of a future work, especially since source inclination is suspected to play a role in the switch of the time lags behaviour as a function of QPO frequency (van den Eĳnden et al. 2017). Finally, the model yields a seemingly worse fit to the data of the rms-energy spectra above 25 keV, and also in general at QPO frequencies below 1.8 Hz. The first issue was also addressed in Karpouzas et al. (2020) and is likely due to the lack of a reflection component in our model. The second issue, which becomes even more apparent from the very large values of the required external heating rate amplitude (| |>100%), can be due to either a break-down of the assumption of the model at QPO frequencies below 1.8 Hz, or to an actual switch in the physical mechanism, as mentioned in the previous sub-section. To answer as many of these questions as possible and retain applicability to large datasets, a more detailed modeling of the disc-corona-jet geometry and interaction in a semi-analytical framework is needed.

CONCLUSIONS
We applied our newly developed Comptonisation model, initially designed to explain the lower kHz QPO in NS LMXBs, to the type-C QPO of the BH LMXB GRS 1915+105. We used a large data-set of timing measurements of this source obtained using archival data from the RXTE satellite. We showed that a spherically symmetric and uniform corona of a finite size is able to explain these timing properties, under the condition that the size of the corona decreases as a function of QPO frequency, down to the critical frequency of 1.8 Hz, where the lags turn from soft to hard. Furthermore, the switch of the lags from soft to hard is moderated by the feedback efficiency of up-scattered photons onto the accretion disc that lies inside the corona. For the first time, we were able to fit the energy dependent rms amplitude and time lags of type-C QPOs, which supports the idea that Comptonisation in the corona should be an essential component of every spectral and timing model.