© ESO 2021

1. Introduction

The X-ray variability in both active galactic nuclei (AGN) and X-ray binaries can be assessed through the use of the power spectral density (PSD). The PSD describes how the X-ray power varies on different timescales (or different temporal frequencies), which depends on the mechanisms responsible for the X-ray production near black holes (e.g. González-Martín & Vaughan 2012). Another timing technique to probe the innermost regions is by measuring X-ray reverberation lags, which are the time delays between the changes in direct continuum and reprocessing echoes from the disc (see Uttley et al. 2014, for a review). Due to the longer distance travelled by the reflected photons, changes in energy bands dominated by reflection, or reprocessing, lag behind changes in the continuum-dominated band. The first hint of such reverberation delays was reported by McHardy et al. (2007) in the AGN Ark 564, followed by the first robust detection by Fabian et al. (2009) in the AGN 1H0707–495.

X-ray reverberation in black hole X-ray binaries was first robustly detected in GX 339–4 by Uttley et al. (2011) when the source was in its hard state. Previous studies of GX 339–4 pointed to the approximate central mass being ≥6 M_⊙ (e.g. Hynes et al. 2003) and a small disc inclination angle (De Marco et al. 2015). Miller et al. (2008) fitted the Suzaku spectra and found that the central black hole has a very high spin, a ∼ 0.998. The X-ray spectroscopic analysis of the hard state spectra from the RXTE archive carried out by García et al. (2015) suggested the black hole spin to be a ∼ 0.95. Spectral fitting of GX 339–4 during its very high flux state using NuSTAR and Swift also suggested a high spin of a ∼ 0.95 (Parker et al. 2016). According to the time-lag analysis, Uttley et al. (2011) found that the disc thermal emission (∼0.3–0.7 keV, soft band) leads the power-law variations (∼0.7–1.5 keV, hard band) on long timescales (> 1s). Mahmoud et al. (2019) assumed that the soft component that leads the power-law emission is a soft Comptontized component. Rapisarda et al. (2016) and Rapisarda et al. (2017) instead modelled it as a variable inner region of the thin disc. However, the disc black-body variations lag behind the power-law variations by a few milliseconds on short timescales (< 1s). This switch from low-frequency hard to high-frequency soft lags is thought to be produced by two distinct mechanisms. While the hard lags are likely due to inward propagating fluctuations (e.g. Kotov et al. 2001; Arévalo & Uttley 2006), the soft lags can be explained by thermal reverberation associated with the longer light-travel time the hard photons take from the central power-law X-ray source to the disc where they are reprocessed into relatively soft black-body emission. The thermal reverberation lags then provide clues to the geometry of the X-ray source and the inner accretion flow close to the event horizon of the central black hole.

The timing analysis of GX 339–4 including additional XMM-Newton observations revealed that the changes in time lags in the high flux states can be characterized as a function of luminosity (De Marco et al. 2015). The luminosity may change with the truncation radius. When the luminosity increases, the truncation radius becomes smaller so that the observed flux increases due to additional disc dissipation and reflection from the inner disc. Meanwhile, the amplitude of reverberation lags become smaller as the photons reflecting off the inner radii of the disc have a smaller average light crossing time. On the other hand, if the truncation radius increases, the luminosity decreases and the reverberation lag increases. Reig et al. (2018) demonstrated that the correlation between the average time lag and the photon index (i.e. the lags increase as the X-ray continuum becomes softer) in the black-hole X-ray binaries could be a result of inverse Comptonization in the base of the jet. Kylafis & Reig (2018) also reported this tight correlation in GX 339–4. Sridhar et al. (2020) found that the temperatures of both inner disc and corona are sensitive to the luminosity of state transition.

The physical properties as well as the geometry of the Comptonizing region in GX 339–4, however, are still unclear. Different fitting techniques and models also lead to different implied values of truncation radius (see Dziełak et al. 2019, and discussion therein). Wang-Ji et al. (2018) investigated the mean spectrum of this source during a failed outburst in 2013 observed by NuSTAR and Swift. They found that a maximum truncation radius could reach ∼37r_g (1r_g = 1GM/c² is the gravitational radius, where G is the gravitational constant, M is the black hole mass, and c is the speed of light). They also reported a smaller truncation radius of ∼3–15r_g during the 2015 outburst. Evolution of the reverberation lag at the end of the 2015 outburst and during the return to quiescence was investigated by De Marco et al. (2017). By converting the amplitude of the lags to light-crossing distance, De Marco et al. (2017) estimated the truncation radius to be ∼40–200r_g. Kylafis & Reig (2018) showed that during the hard and hard-intermediate states the inner disc could extend inwards, shrinking the hot inner flow and the jet base. García et al. (2019) carried out a spectroscopic analysis of this source when it went through the failed outburst in 2017. They found that a dual-lamppost model provides a better fit than the standard single lamppost source and inferred that the truncation radius reaches a few gravitational radii as the luminosity increases in the hard state.

While some of the spectroscopic studies of the refection component in the bright hard state suggested that the inner disc could extend very close to the innermost stable circular orbit (e.g. García et al. 2015; Steiner et al. 2017; Wang-Ji et al. 2018), there is also previous literature that reported a significantly larger truncation radius of for example > 35r_g (Tomsick et al. 2009; Basak & Zdziarski 2016) during the hard state. The truncation of the optically thick, geometrically thin disc in the hard state can be the result of thermal conduction of heat from the corona that causes the inner disc to evaporate (Meyer et al. 2000). Mahmoud et al. (2019) modelled the hard state spectral-timing data of GX 339–4 taking into account the time lags due to reverberation and found evidence of the truncation radius in the order of ∼20r_g.

Veledina (2016, 2018) investigated the interference of X-rays under the propagating-fluctuations framework and found that the presence of multiple sources (e.g. disc Comptonization and synchrotron Comptonization sources) could produce a complex PSD shape seen in X-ray binaries. Reverberation signatures can also be imprinted in the PSD profiles (Papadakis et al. 2016; Chainakun 2019). The aim of our work is to develop a PSD model, taking into account both propagating-fluctuations and reverberation effects, to predict the changes in the accretion disc and the inner hot flows of GX 339–4. It is clear that there is a large scatter in the reported values of the truncation radius in GX 339–4 (see e.g. Wang-Ji et al. 2018, and references therein). Here we focus on the very last stages of the soft-to-hard transition as the source returns to quiescence, similar to the observations analysed by De Marco et al. (2017).

The GX 339–4 observations used here and the data reduction are explained in Sect. 2. The theoretical PSD models are presented in Sect. 3. The fitting results are shown in Sect. 4, followed by the discussion in Sect. 5. The conclusions are drawn in Sect. 6.

2. Observations and data reduction

The data analysed in this paper were selected from those that were in the final phase of the GX 339–4 outburst observed during August–September 2015 using the EPIC-pn detector on board XMM-Newton, and were obtained from XMM-Newton Science Archive¹. The selected observations are listed in Table 1. All observational data were cleaned following the standard data reduction method as described on the XMM-Newton data analysis webpage² using Science Analysis Software (SAS) version 18.0.0 with the latest calibration files. We created EPIC-pn light curves using the SAS task EVSELECT with a selection criteria of PATTERN ⩽4 and a time bin size of 6 ms. The source extraction region for the data in timing mode was defined as the data that have 28 ⩽ RAWX ⩽ 48, while the region for those observed in small window mode was defined as the circular area centred on the source position with a radius of 40″.

Since the source is very bright, the light curves should be dominated by the source counts. Here, we followed the data reduction outlined in De Marco et al. (2017). To maximize the data signal to noise ratio, we did not remove background flaring events from the observational data. In fact, the fraction of the useful exposure time affected by the flaring events for each observation is ≲7%, except for the observation O1 in which ∼60% of the exposure time is affected. Nevertheless, we checked that timing products obtained from the observation with and without removing the flaring events, specifically the power spectra, are well consistent. This verifies that the flares did not affect the GX 339–4 timing properties in our timing analysis. We also note here that all observations were not affected by pile-up. In fact, the source count rates during observations O1 – O4 are substantially below the maximum count rate allowed for an observation in the timing mode of the pn detector (< 240 count s⁻¹). For observations O5 – O6, which were observed in small window mode, we ran the SAS task EPATPLOT to check and confirm that the observations were not significantly affected by the pile-up. Therefore, the useful exposure time for each observation after performing the data reduction is shown in Col. 4 of Table 1. The light curves in two different energy bands – 0.3–0.7 keV and 0.7–1.5 keV – were extracted from all pn observations. These bands are dominated by thermal reverberation and power-law continuum, respectively, as previously reported by De Marco et al. (2017).

Finally, based on the light curves obtained, we created the power spectra using the FTOOLS task POWSPEC³. In brief, each light curve was divided into a number of segments with an interval of 41 s and a time bin size of 6 ms; then these segments were converted into power spectra and averaged to create the single power spectrum for each light curve. The power spectra were rebinned by a factor of 1.04 dex for the data below 2 Hz and 1.2 dex for the data above 2 Hz. We ignored data points consistent with zero (or negative) power. The power spectra obtained from this method were then used as the basis for further analysis.

3. Theoretical model

3.1. Geometry setup

To avoid having too many free parameters, we fixed the black hole mass of GX 339–4 to be 10 M_⊙, the inclination angle i = 30°, and the black hole spin a = 0.998 (e.g. Hynes et al. 2003; De Marco et al. 2015; García et al. 2015; Parker et al. 2016; Zdziarski et al. 2019). Therefore, the innermost stable circular orbit (ISCO) is at ∼1.23r_g. The disc is geometrically thin and optically thick (Shakura & Sunyaev 1973) ranging between 400r_g and r_trc, where r_trc is the truncation radius. Inside r_trc the disc is replaced by two hot-flow zones responsible for the emission of soft spectrum (Γ_sz) and hard spectrum (Γ_hz) ranging between r_trc–r_sh and r_sh–1.23r_g, respectively. The Γ_sz and Γ_hz are the photon indices of the X-ray continuum emitted from the inner soft and hard zones, respectively, and r_sh is the transition radius between the two zones. The accretion disc varies on long timescales while the turbulent inner hot flows vary intrinsically on relatively fast timescales. The parameter t_p characterizes the time taken for fluctuations to propagate from r_trc to r_sh. Our geometric setup is presented in Fig. 1.

Fig. 1. Geometry of our model. The accretion disc varies on long timescales and extends from 400rg–rtrc. The inner regions from rtrc–rsh and rsh–1.23rg are turbulent inner flows varying intrinsically on fast timescales that are responsible for the soft and hard spectral emission, respectively. The accretion disc provides seed photons to the inner hot flows, and reprocesses and reflects the Comptonized emission photons coming back from the inner regions. The parameter tp is the time taken for fluctuations to propagate from rtrc to rsh.

We consider two distinct mechanisms that take place in our system: fluctuations in the mass accretion rate and the thermal reverberation. The accretion disc provides seed photons to the inner hot-flow zones, and reflects the Comptonised photons travelling outwards from the hot flows. The continuum fluxes in the energy band of interest E_j emitted from the spectrally soft and the hard hot-flow zones are

$$\begin{array}{cc}\hfill F_{\mathrm{sz}}(E_j)& \propto {\displaystyle {\int }_{E_{j,\mathrm{l}ow}}^{E_{j,\mathrm{h}igh}}}{E_j^{\prime }}^{-{\mathrm{\Gamma }}_{\mathrm{sz}}}\mathrm{d}{E}_j^{\prime },\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill F_{\mathrm{hz}}(E_j)& \propto {\displaystyle {\int }_{E_{j,\mathrm{l}ow}}^{E_{j,\mathrm{h}igh}}}{E_j^{\prime }}^{-{\mathrm{\Gamma }}_{\mathrm{hz}}}\mathrm{d}{E}_j^{\prime },\hfill \end{array}$$(2)

where E_j,low and E_j,high are the lowest and highest energy of the specific band of interest, respectively.

Recently, Ingram et al. (2019) showed that there could be some misclassified photons in the energy band < 0.7 keV that should belong to higher energy bands due to the fact that the response matrix of XMM-Newton is non-diagonal. Taking into account this effect as well as the absorption column found in GX339–4 could change the flux normalization in a particular soft energy band. In this work, these effects are included indirectly by employing the reflected response fraction obtained from the spectral fits by De Marco et al. (2017) to produce the PSD, which will be described later in Sect. 4. Therefore, the PSD results here should tie to the mean spectrum fits of the data that have already folded through the instrument responses and that should incorporate the flux contamination induced by the effects of the response matrix as well.

3.2. Mass accretion rate fluctuations and PSD

When fluctuations in mass accretion rate, ṁ, propagate through the accretion disc, they induce variations in the seed photon emission. The disc seed photons are then Compton up-scattered inside the hot corona (or hot flow) producing the Comptonized photons whose variability is directly proportional to that of the mass accretion rate fluctuations. The observed continuum flux in the energy band E_j then consists of two components associated with the emission from soft and hard zones, which can be described by

$$\begin{array}{c}\hfill x_0(t,E_j)=x_{\mathrm{hz}}(t,E_j)+x_{\mathrm{sz}}(t,E_j),\end{array}$$(3)

where

$$\begin{array}{cc}\hfill x_{\mathrm{hz}}(t,E_j)& \propto F_{\mathrm{hz}}(E_j)\dot{m}(t),\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill x_{\mathrm{sz}}(t,E_j)& \propto F_{\mathrm{sz}}(E_j)\dot{m}(t+t_{\mathrm{p}})\ast h(t).\hfill \end{array}$$(5)

The variations of mass accretion rate fluctuations is ṁ(t) and h(t) is the filter function that screens out the high-frequency variability of signals in the soft hot-flow zone so that the harder, inner zone produces the hard continuum spectra varying more at higher frequencies, which is the traditional property observed in PSD of AGNs (González-Martín & Vaughan 2012) and X-ray binaries (Cui et al. 1997; Remillard & McClintock 2006). Also, h(t) produces the different X-ray variability in the two hot-flow components, so that the oscillatory interference features in the PSD can be produced. There is no one unique framework to explain distinct emission regions and the peaks seen in the PSD. This interference framework was previously proposed by Veledina (2016, 2018), and it was suggested by Mahmoud & Done (2018) that it can encompass other current models (e.g. Ingram & Done 2012) within the parameter space (see Mahmoud et al. 2019, for further discussion). The variability produced from the soft hot-flow zone is modelled by a convolution (* sign) of the filter function with the mass accretion rate. The parameter t_p is a characteristic propagation time from the r_trc to r_sh where the hard Comptonization operates.

The PSD is estimated by the modulus squared of the discrete Fourier transform of the light curve (e.g. Nowak et al. 1999; Emmanoulopoulos et al. 2013), which can be written as

$$\begin{array}{c}\hfill P_0(f,E_j)\propto {|X_0(f,E_j)|}^2/{\overline{x_0}}^2,\end{array}$$(6)

where upper case letters represent the quantities in the frequency domain corresponding to those in the time domain written in lower case letters. According to Eqs. (3)–(6), the PSD under this framework can be calculated via (e.g. Veledina 2016)

$$\begin{array}{c}\hfill P_0(f,E_j)=\frac{1+{ϵ}^2{H}^2(f)+2ϵH(f)cos(2\pi f{t}_{\mathrm{p}})}{{(1+ϵ)}^2}|\dot{M}(f){|}^2,\end{array}$$(7)

where it is normalized in terms of the squared fractional rms (Miyamoto et al. 1991) and h(t) is a zero-lag function, so that H(f) is real as this is instrumental in the derivation of Eq. (7). The relative flux ratio of the X-ray continuum emitted from two different hot-flow zones is ϵ = N_F(F_sz/F_hz). While F_sz/F_hz depends on Γ_sz, Γ_hz, and the energy band being considered (see Eqs. (1)–(2)); we assume the normalization factor N_F is the same for all energy bands in the same observation. The N_F in the same energy band, however, could be different between different observations. The Fourier form of the filter function is H(f), which can be expressed in the form (Veledina 2016)

$$\begin{array}{c}\hfill H(f)=\frac{1}{{(f/f_{\mathrm{filt}})}^4+1}\end{array}$$(8)

so that the transmitted signals above the frequency f_filt are damped.

The variability power of the mass accretion rate significantly decays at the viscous frequency (Ingram & van der Klis 2013)

$$\begin{array}{c}\hfill f_{\mathrm{visc}}(r_n)=r_n^{-3/2}{(H/R)}^2\alpha /2\pi ,\end{array}$$(9)

where H/R is the disc scale-height ratio and α is the viscosity parameter. The radius r_n is in units of r_g so the viscous frequency f_visc(r_n) is in units of 1/t_g, where t_g = GM/c³. One can convert gravitational units to physical units when the black hole mass is known. For a black hole mass of 10 M_⊙, 1r_g = 1.48 × 10⁴ m and 1t_g = 4.9 × 10⁻⁵ s.

To describe the broken power-law shape usually seen in the PSD, the form of |Ṁ|²(f) is assumed to be

$$\begin{array}{c}\hfill |\dot{M}{(f)|}^2=\frac{1}{1+{(f/f_{\mathrm{break}})}^{\beta }},\end{array}$$(10)

where f_break is the break frequency and β is an arbitrary index. We fix β = 1.5 when investigating other model parameters, but allow it to be free when we fit the data. The Ṁ(f) represented in Eq. (10) is then the frequency domain form of the driving signal that is used to activate the X-ray continuum emission and variability from the hot flows.

There are two important parameters in our model giving the radial positions of the truncation radius, r_trc, and the transition radius where the hot flows change from the spectrally soft to spectrally hard zone, r_sh. We then relate the break and the filter frequency to the viscous timescales at the truncation and transition radius, respectively,

$$\begin{array}{cc}& f_{\mathrm{break}}=f_{\mathrm{visc}}(r_{\mathrm{trc}}),\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill & f_{\mathrm{filt}}=f_{\mathrm{visc}}(r_{\mathrm{sh}}).\hfill \end{array}$$(12)

The characteristic time for fluctuations that originate at a radius r to propagate across a radial distance Δr is Δt = Δr/[rf_visc(r)] (Ingram & van der Klis 2013). Assuming the characteristic time for fluctuations to propagate through the soft zone depends on the viscous timescale at the truncation radius, the value of t_p in Eqs. (5) and (7) then can be estimated as

$$\begin{array}{c}\hfill t_{\mathrm{p}}=(r_{\mathrm{trc}}-r_{\mathrm{sh}})/(r_{\mathrm{trc}}{f}_{\mathrm{break}}).\end{array}$$(13)

Figure 2 shows examples of the model PSD for different energy bands when the photon indices of the continuum emitted from soft and hard zones are Γ_sz = 2.5 and Γ_hz = 2.0, respectively. We assume the disc is truncated at r_trc = 10r_g and the transition radius is r_sh = 5r_g. The disc parameter (H/R)²α is set to be 0.01. The model can produce larger high-frequency power in higher energy bands, which are the traditional PSD properties seen in both AGN and X-ray binaries. Furthermore the dips and humps are naturally produced by the model. We fix N_F = 0.5 so the parameter ϵ = N_F(F_sz/F_hz) that regulates the importance of the filter level depends only on F_sz/F_hz, which is different in each energy band. Generally, the dips and humps could also depend on r_trc and r_sh that determine f_break, f_filt, and the characteristic propagation time t_p.

Fig. 2. Power spectral density in the 0.3–0.7 keV energy band (black line), 0.7–1.5 keV (blue line), and 3–5 keV (red line) with rtrc = 10rg and rsh = 5rg. The soft and hard zones of the inner flow have the spectral photon indices Γs = 2.5 and Γh = 2.0, respectively. We fix NF = 0.5 so that all energy bands are different due to the values of the intrinsic flux ratio Fs/Fh.

How the geometry of the inner flows affects the PSD in the 0.3–0.7 keV band when ϵ = 0.05 and (H/R)²α = 0.01 is shown in Fig. 3. We note that when ϵ is fixed, there can be different pairs of Γ_sz and Γ_hz that produce the same PSD profile for each value of N_F. To illustrate this, Fig. 3 (top panel) shows the Γ_sz versus Γ_hz plots that provide ϵ = 0.05 for different N_F. In principle, we expect Γ_sz >  Γ_hz, so the black solid line marks the upper limit of N_F of which above this value the model returns Γ_sz <  Γ_hz in this energy band, which is unrealistic and can be neglected. Figure 3 (middle panel) shows that stronger high-frequency power is obtained with decreasing r_trc. Since we fix r_sh = 2r_g, larger r_trc results in a longer propagation-time delay, t_p, through the inner soft to the inner hard hot-flow zones. Increasing time t_p results in more humps being imprinted on the PSD due to the effect of the cosine term in Eq. (7). This is also true when we fix r_trc = 20r_g, but vary r_sh and hence also t_p (Fig. 3, bottom panel). Moreover, the smaller r_sh is, the higher the frequency at which the oscillations are significantly filtered. This is because the f_filt is larger for smaller r_sh. Decreasing r_sh means the inner soft zone could produce X-rays on shorter characteristic timescales to interfere with X-rays from the hard hot flows, leading to the oscillation, interference structures on the PSD profiles towards higher frequencies.

Fig. 3. Top panel: Γsz vs. Γhz plots that give ϵ = 0.05 for different NF values for a 0.3–0.7 keV band of interest. Middle panel: corresponding PSD in 0.3–0.7 keV band varying with rtrc. Bottom panel: corresponding PSD in the same band varying with rsh.

3.3. Simulated PSD including reverberation

Since the radial size of the disc is relatively large compared to that of the inner hot flows, the source illuminating the disc can be approximated as a central illuminating source (e.g. Gardner & Done 2014; Mahmoud et al. 2019). Therefore, the time delay of photons travelling from the central source to be reprocessed by the disc at radius r is given by (Welsh & Horne 1991)

$$\begin{array}{c}\hfill \tau =\frac{r}{c}(1-sinicos\varphi ),\end{array}$$(14)

where i is the inclination angle measured from the observer’s line of sight to the disc axis and ϕ is the azimuthal angle between each specific point on the disc and the projection of the line of sight onto the disc. We note that Eq. (14) does not include general relativistic effects but is still acceptable since we do not specify the exact geometry of the flows. The impulse response function due to reverberation, ψ(τ), can be produced by summing the photon counts over a given range of radii r_trc ≤ r ≤ 400r_g and all azimuthal angles 0 ≤ ϕ ≤ 2π as a function of the time delay using Eq. (14). At radii smaller than r_trc, the disc is replaced by the hot flow and X-ray reverberation does not take place.

We introduce the parameter γ to describe the emissivity of the accretion disc so that the number of emitted photons from each annulus of area 2πrdr is proportional to r^−γ. Higher γ means the X-ray reflection is more concentrated at the inner part of the disc. Meanwhile, γ = 0 means that the same number of photons is produced in each area of 2πrdr and that emissivity, in terms of flux per unit area on the disc, is proportional to r⁻¹. Examples of the disc-response functions are presented in Fig. 4. Here we fix M_BH = 10 M_⊙ and i = 30° while allowing the r_trc to be a free parameter. The solid and dotted lines represent the case of γ = 0 and 1.5, respectively.

Fig. 4. Simulated impulse disc responses in 0.3–0.7 keV band varying with the disc truncation radius. The solid and dotted lines represent the cases when the emissivity index γ = 0 and 1.5, respectively. The inclination angle is fixed at i = 30° and MBH = 10 M⊙.

The impulse reverberation response also acts as a filter on the driving signal that affects the shape of the observed PSD (Papadakis et al. 2016; Chainakun 2019). According to the convolution theorem, the Fourier transform of an observed light curve including reverberation effects can be written as (Uttley et al. 2014; Papadakis et al. 2016)

$$\begin{array}{c}\hfill X(f,E_j)=X_0(f,E_j)\mathrm{\Psi }(f,E_j),\end{array}$$(15)

where X₀(f) and Ψ(f) are the Fourier transform of the driving signal and the impulse response, respectively. The observed PSD in the energy band E_j is then given by

$$\begin{array}{c}\hfill P_{\mathrm{obs}}(f,E_j)=\frac{|X(f,E_j){|}^2}{(1+R(E_j){)}^2}=\frac{|X_0(f,E_j){|}^2{|\mathrm{\Psi }(f,E_j)|}^2}{(1+R(E_j){)}^2},\end{array}$$(16)

where R(E_j) is the reflected response fraction defined as the (reflection flux)/(continuum flux) measured in the specific energy band E_j, which can vary between observations.

The 0.3–0.7 keV PSD predicted by the model that includes both effects of disc fluctuations and reverberation are presented in Fig. 5 (solid lines). We fix ϵ = 0.1 and (H/R)²α = 0.01, but vary γ, r_trc, and r_sh. For comparison, the corresponding PSD with reverberation excluded is also shown as the dotted lines. It is clear that the reverberation acts as a filter that reduces the high-frequency power in the observed PSD, which is independent of the input signals (see also the discussion in Papadakis et al. 2016; Chainakun 2019). Interestingly, if the emissivity γ is too high (i.e. the source illumination is too centrally concentrated), most photons would reflect from the inner disc so the reverberation delays would be quite short and, consequently, relatively small effects of reverberation are seen in the PSD profiles towards low temporal frequencies (i.e. long timescales). We note that the reverberation features are imprinted on the PSD at high temporal frequencies associated with the short timescales of the disc reflection. Given M_BH = 10 M_⊙, the significant drop of power due to reverberation at frequencies around 10–30 Hz should be the beginning of the main dip of the oscillatory structures in the PSD produced by reverberation such as in case of AGN (Papadakis et al. 2016; Chainakun 2019). Furthermore, the results in Fig. 5 (bottom panel) suggest that the overall features of the PSD are quite different for different combinations of key model parameters such as r_trc, r_sh, ϵ, and γ.

Fig. 5. Top panel: total PSD in 0.3–0.7 keV band varying with γ when i = 30°, rtrc = 10rg, rsh = 5rg, ϵ = 0.1 and (H/R)2α = 0.01. Bottom panel: total PSD in 0.3–0.7 keV band varying the geometry of the system with γ = 1.5. All dotted lines represents the corresponding cases with reverberation excluded. The reduction of power at high frequencies is produced by reverberation.

4. Fitting results

Our model consists of nine parameters in total: truncation radius, r_trc, soft-to-hard transition radius, r_sh, disc parameter, (H/R)²α, normalization flux, N_F, photon index of the hard zone, Γ_hz, photon index of the soft zone, Γ_sz, low-frequency break index, β, reflected response fraction, R(E_j), and disc emissivity index, γ. To minimize the free parameters, the R(E_j) is fixed to the values obtained from the spectral fits by De Marco et al. (2017), which are approximately 4.0 and 1.0 in the 0.3–0.7 keV and 0.7–1.5 keV bands, respectively. The reflected response fractions employed here describe the effects of contamination flux between cross-components in the reverberation-dominated and continuum-dominated energy bands. In principle, they should also be able to describe the effect of misclassifying photons due to the non-diagonal response matrix of the instrument that plays a role in changing contamination flux as well. The effects of changing the photon index on the PSD are similar to changing N_F, which is a free parameter. Since the photon index of the X-ray continuum in these GX 339–4 observations was previously reported to vary between 1.5–2 (e.g. De Marco et al. 2017), we fix Γ_hz = 1.5 and Γ_sz = 2.0. This helps avoid the degeneracies in the model and also avoid predicting a spectrum that is wildly different from the observed spectrum. The parameter N_F is fixed between two energy bands, but allowed to vary across the observations.

According to Veledina (2016), the disc and viscous parameters of GX 339–4, via modelling the PSD in 2–15 keV, were found to be H/R >  0.23 and 0.01 <  α <  0.5. We consider the cases when 0.0005 ≲ (H/R)²α ≲ 0.1, which is approximately between the lower limit and a reasonable upper limit (H/R) = 1 and α = 0.1, but still in the range of what was previously reported (Veledina 2016). Firstly, the global grids of the model have been produced and the PSD data in two energy bands are simultaneously fitted in ISIS, the Interactive Spectral Interpretation System (Houck & Denicola 2000). We find that for the majority of fits the value of (H/R)²α is 0.005 and γ is 0. To improve the fit while limiting computation time, finer local grids of the model are produced independently for each observation, but this time both (H/R)²α and γ are fixed at the values constrained by the global grids, under the assumption that they do not vary much during these observations. The fitting is repeated with the finer, local grids for each observation with fixed (H/R)²α = 0.005 and γ = 0. The fitting results are presented in Fig. 6. The corresponding best-fit parameters are listed in Table 2. We found an increase in the truncation radius from ∼10r_g to ∼55r_g during O1–O6 as the source is decreasing in flux towards the end of the outburst. An increase of the soft–hard transition radius and a variation of N_F are also found during O1–O6. Meanwhile, β seems to decrease towards the end of the outburst.

Fig. 6. Data, model and residuals from fitting the PSD model to 0.3–0.7 keV (blue) and 0.7–1.5 keV (red) PSD data with (H/R)2α fixed at 0.005 and γ fixed at 0 as obtained by the global fits. The PSD displayed here are Poisson noise-subtracted and normalized in terms of squared fractional rms (Miyamoto et al. 1991). Data points consistent with zero (or negative) power have been excluded. The best-fit model parameters for each observation are listed in Table 2.

5. Discussion

For GX 339–4, there is strong evidence that the central black hole has a spin parameter close to the maximal value of a = 0.998 (e.g. Miller et al. 2008; García et al. 2015; Parker et al. 2016). During the high soft state, the inner disc could extend into the ISCO (Plant et al. 2014). Recently, Sridhar et al. (2020) found that the truncation radius remains low (< 9r_g) throughout the bright intermediate state (i.e. the width of the Fe K emission line is nearly constant throughout this transition indicating a quasi-static truncation radius). The evolution of the accretion disc during the intermediate state then should be driven by the variation in accretion rate rather than by the changes in truncation radius. The data analysed here captured the very last stages of the soft-to-hard transition (O1–O2) and the return to quiescence where the source was decreasing its luminosity (O3–O6).

The data analysed here are similar to those used for the reverberation lag analysis of De Marco et al. (2017), where they reported a decrease of the reverberation lag as a function of the source luminosity. Our results suggest that the truncation radius increases from ∼10r_g to ∼55r_g (see Table 2) as the outburst proceeds and the flux decreases from O1–O6. The truncation radius versus the Eddington-scaled luminosity is shown in Fig. 7. The truncation radii from De Marco et al. (2017) are also shown in Fig. 7; these are estimated by converting lag amplitudes to the light-crossing distance from the centre to the inner edge of the disc, assuming a black hole mass of 10 M_⊙. The truncation radii being directly converted from the lag amplitude are significantly larger than those constrained by our PSD model. This is what we expected because the estimates by De Marco et al. (2017) are highly simplistic whereas our model accounts for the dilution. Our results support the commonly agreed framework that the disc truncation radius increases as the luminosity decreases (e.g. Wang-Ji et al. 2018, and references therein).

Fig. 7. Comparison of the truncation radius vs. Eddington-scaled luminosity for GX 339–4 inferred by this work (red diamonds), De Marco et al. (2017) via approximate reverberation measurement (circles), and Mahmoud et al. (2019) via spectral-timing fits (crosses).

Mahmoud & Done (2018) proposed a framework to explain the time-averaged spectrum, the PSD in different energy bands, and the time lags, simultaneously. Their model includes at least three distinct Compton components associated with three specific radii: the truncation radius, the inner radius of the hot flow, and the jet-launch radius. Later on, Mahmoud et al. (2019) used their spectral-timing model, taking into account the effects of reverberation, to fit the O1 data, and found a truncation radius of ∼19.5r_g. Their result is also plotted in Fig. 7. The truncation radius for O1 constrained by our PSD model is ∼10r_g, which is smaller but still comparable to Mahmoud et al. (2019). Spectral fits of the 2015 outburst data observed by NuSTAR and Swift were carried out by Wang-Ji et al. (2018), where the truncation radius was found to be ∼3–15r_g. The accretion disc implied by the spectroscopic analysis (e.g. Wang-Ji et al. 2018; García et al. 2019) is likely to be truncated at smaller radii than all of those suggested by timing analysis.

In principle, the truncation radius, r_trc, and the X-ray source height cannot be well constrained independently via reverberation lag measurements only since both of them directly determine the measured light-travel distance. Furthermore, from Table 2, the frequency break index β decreases as the outburst proceeds, suggesting that there is a variation in accretion rate as well as in r_trc. Decreasing β means the variability power in the frequency domain induced by the mass accretion rate fluctuations varies, from O1–O6, in a way that produces relatively high power at frequencies f >  f_break (see Eq. (10)). Recently, Mushtukov et al. (2018) described the propagation of the fluctuations in the disc using the diffusion equation whose solution was derived by the method of Green functions. The suppression of variability at high frequencies was produced by the Green functions and was stronger for a lower parameter associated with radial kinematic viscosity. Our parameter β that suppresses variability at high frequencies above the viscous frequencies could then describe the dependence of kinematic viscosity and possibly tie back to the Green function treatments such as in Mushtukov et al. (2018). If the changes in luminosity are also driven by the variation in the viscous parameters, the frequency of f_break is changed as well (see Eqs. (9)–(12)). We note that to limit the number of free parameters and avoid model degeneracy we fix (H/R)²α = 0.005 (i.e. assuming this parameter does not significantly vary during O1–O6). This value is obtained from the fits using the global grid and is in the acceptable range reported in the literature (e.g. Veledina 2016).

Furthermore, Mahmoud et al. (2019) fitted the spectral-timing data of O1 and reported the transition radius where the spectrally soft flow changes to hard inner flow to be ∼6.6r_g. We estimate r_sh for O1 to be ∼7r_g, which is in agreement with their result. Starting from r_sh = 7r_g at O1, our fitting also suggests that the soft-hard transition radius keeps increasing to r_sh = 29r_g towards O6. The size of the spectrally soft inner flow and hard inner flow zones can be estimated as r_trc − r_sh and r_sz − 1.235r_g, respectively. The comparative size of these flows is shown in Fig. 8. As the flux decreases from O1–O6, the size of the hard zone expands from ∼5r_g to 27r_g, which is ∼1.1 − 2.2 times the radial extent of the soft zone that expands from ∼3r_g to 26r_g.

Fig. 8. Comparative size of inner hot-flow zones between different observations (O1–O6) implied from the best-fit parameters. Error bars are small and are excluded in this plot.

Since we fixed Γ_sz = 2.0 and Γ_hz = 1.5 (consistent with the values reported by De Marco et al. 2017), variations in N_F (see Table 2) could actually be due to the variations in spectral indices of the hot flows between different observations. The viscous timescale at the truncation radius, or the outer radius of the hot flow, sets the low-frequency break of the PSD. A low-frequency quasi-periodic oscillation (QPO) due to Lense-Thirring precession of the hot flow (Ingram et al. 2009) is sometimes observed in black hole X-ray binaries (BHXBs) in such a way that the QPO frequency moves with the low-frequency break. However, the QPO was previously found to be weak throughout O1–O6 (De Marco et al. 2017; Mahmoud et al. 2019), so they are not included in our fits. Additional Lorentzian or Gaussian components associating with these QPOs could definitely improve the fits, especially for the O1, but the trend of the PSD model and the implied parameters should not be significantly different.

Our model could allow the emissivity indices to be free but the exact values should depend on the assumed source geometry. Different geometries of the flow produce different illumination patterns on the disc and hence γ could also vary between observations. However, we find that γ ∼ 0 is obtained from the global fits in all six observations. It should be noted that γ = 0 produces the same number of photons in each annulus of area 2πrdr and that the emissivity in terms of flux per unit area on the disc is r⁻¹. This may suggest that the inner flows are significantly vertically extended in such a way that the number of photons illuminating the outer part of the disc is large enough to make the emissivity profile likely flat. It may also require the inner flow to be more vertically extended than the outer flow so that some photons that are supposed to illuminate the thin disc near r_trc are obscured by the outer flow, flattening the emissivity especially at the inner parts of the disc. We caution that this is just a rough approximation and interpretation suggested by the global fits. In principle, γ should not be thought of as an independent parameter but instead should be tied to the exact geometry of the flow. An additional function may be required to tie γ to the characteristic radii of the flow. In that case, the exact geometry of the flow and also the disc (e.g. in the case of a flared disc) must be explicitly specified, which makes the model much more complicated. We then choose to fix γ to the values obtained from the global fits in this work. Nevertheless, the disc emissivity profile is worth investigating in the future, using high density disc models (Jiang et al. 2019).

Dual lamppost sources illuminating the disc were found to be able to explain the spectroscopic data of GX 339–4 in some observations compared to the single lamppost source (García et al. 2019). Therefore, the X-ray source could be extended and the realistic impulse responses need to be produced through the ray-tracing simulations (Wilkins et al. 2016; Chainakun & Young 2017; Chainakun et al. 2019; García et al. 2019). Here, we assume the size scale of the source, or the hot flows, to be relatively small compared to the radial size of the accretion disc (Gardner & Done 2014; Mahmoud et al. 2019). However, Chainakun (2019) reported that the reverberation signatures on the PSD profiles (e.g. the main dip and the oscillation features) would be relatively weak and become more difficult to probe if the source is significantly more extended. Using more realistic, ray-traced impulse responses from different spatially extended source geometries is beyond the scope of this paper, but is planned for the future.

Although we consider two temperature inner flows plus thermal reverberation, the model can naturally produce the wiggles in the PSD profile that are usually seen in the observed data. The bumpy PSD could also be produced by other realistic treatments such as when the propagating fluctuations are solved through the diffusion equation using the Green function (Mushtukov et al. 2018) and when extra variabilities are injected at different characteristic radii (e.g. Rapisarda et al. 2016, 2017). Since the number of free parameters have been minimized to avoid degeneracies, we see that the model provides very good constraints on each key parameter with small statistical errors. The power extracted from observational data towards the high-frequency end at some points turns out to be zero or negative, which values are excluded in our analysis. New data delivered by upcoming X-ray observatories could provide higher quality data (e.g. mean spectrum, lags, and PSD) that could be robustly constrained with the model. High quality PSD extending to higher frequencies may reveal the clear dips and oscillatory structures produced by reverberations that are imprinted on the PSD.

At the moment, we do not directly consider the effects of the response matrix of XMM-Newton that can misclassify photons at low energies (see Ingram et al. 2019, and discussion therein). This, however, could affect the timing data in the same way in which the contamination flux in the thermal reverberation dominated and continuum dominated bands is changed, which has been taken into account in this work by employing the reflected response fraction to determine the contamination ratio between these cross-components. Applying the response matrix directly to the model would be more straightforward and is worth investigating in the future. Last but not least, our model can be further developed to predict the lag-frequency spectrum between these two bands. Fitting the PSD and the lags simultaneously (or the full cross-spectrum) is challenging and optimal statistics must be investigated to account for significantly different numbers of bins between different data sets, otherwise the fitting might be biased towards the best-fit parameters for the PSD data that have significantly larger numbers of data points. The development of a more self-consistent model to simultaneously explain multi-timing data including the effects of the instrument response, however, is planned for the future.

6. Conclusion

In this study we developed a PSD model that is then fit to the data of GX 339–4 observed during the end of the 2015 outburst. We modelled a truncated accretion disc from r_trc–400r_g. Inside r_trc, spectrally soft and hard hot flows take place down to the ISCO. The model incorporates both disc-fluctuation and reverberation signals. The fluctuations inside the hot flows are propagated inwards on the viscous timescale of the truncation radius. The model can qualitatively reproduce the traditional PSD profiles exhibiting increased high-frequency power for higher energy bands. Furthermore, stronger high-frequency power could be produced with decreasing r_trc, which is also expected. Including the reverberation signals produces a dip at the high frequency end that should be the beginning of oscillatory structures imprinted in the PSD profiles as discussed by Papadakis et al. (2016) and Chainakun (2019).

To produce the model grids, the reflected response fractions of the 0.3–0.7 and 0.7–1.5 keV bands are fixed at the values constrained by the time-averaged spectral analysis (De Marco et al. 2017), as well as the photon index of the X-ray continuum associated with the soft and hard hot-flow zones. In doing this, the PSD models are tied to the realistic parameters obtained from the spectral fitting and the contribution of continuum flux in each energy band as the dilution of reverberation is properly taken into account. The PSD data for both energy bands are fit simultaneously for each observation. We find that the disc parameter (H/R)²α = 0.005 could provide a good fit for the majority of observations using the global grid of the model, hence it is fixed when the finer, local grids are produced independently for each observation and the fitting is repeated.

The fitting results suggest that the truncation radius moves outwards from 10–55r_g as the source luminosity decreases from O1–O6. Although the trend of increasing r_trc with decreasing luminosity is in agreement with previous studies, the values we obtain are smaller than previous reverberation lag analysis (De Marco et al. 2017) and spectral-timing modelling (Mahmoud et al. 2019), but are larger than some of those constrained using spectroscopic data alone (e.g. Wang-Ji et al. 2018; García et al. 2019). We find that the transition radius also increases from 7–29r_g during O1–O6, meaning that the sizes of inner hard hot flows increase from ∼5–27r_g; these spectrally hard hot flows always span a slightly larger radial range than the spectrally soft hot flows by a factor of ∼1.1–2.2. The current PSD model can be straightforwardly adapted for different source and inner flow geometries, such as dual lamppost cases, which may suit the unique data of different X-ray binaries and AGNs.

Acknowledgments

The calculations in this work were carried out using the BlueCrystal supercomputer of the Advanced Computing Research Centre, University of Bristol, UK. All data analyzed are based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. PC thanks the Thailand Research Fund (TRF) for support under grant number MRG6280086, and acknowledges useful discussions with Utane Sawangwit. WL thanks for the financial support from the Faculty of Science, Srinakharinwirot University. We thank the anonymous referee for their comments which led to a clarification of some important points in the paper and has inspired some interesting work for the future.

References

All Tables

All Figures

Fig. 1. Geometry of our model. The accretion disc varies on long timescales and extends from 400rg–rtrc. The inner regions from rtrc–rsh and rsh–1.23rg are turbulent inner flows varying intrinsically on fast timescales that are responsible for the soft and hard spectral emission, respectively. The accretion disc provides seed photons to the inner hot flows, and reprocesses and reflects the Comptonized emission photons coming back from the inner regions. The parameter tp is the time taken for fluctuations to propagate from rtrc to rsh.

In the text

	Fig. 2. Power spectral density in the 0.3–0.7 keV energy band (black line), 0.7–1.5 keV (blue line), and 3–5 keV (red line) with rtrc = 10rg and rsh = 5rg. The soft and hard zones of the inner flow have the spectral photon indices Γs = 2.5 and Γh = 2.0, respectively. We fix NF = 0.5 so that all energy bands are different due to the values of the intrinsic flux ratio Fs/Fh.
In the text

	Fig. 3. Top panel: Γsz vs. Γhz plots that give ϵ = 0.05 for different NF values for a 0.3–0.7 keV band of interest. Middle panel: corresponding PSD in 0.3–0.7 keV band varying with rtrc. Bottom panel: corresponding PSD in the same band varying with rsh.
In the text

	Fig. 4. Simulated impulse disc responses in 0.3–0.7 keV band varying with the disc truncation radius. The solid and dotted lines represent the cases when the emissivity index γ = 0 and 1.5, respectively. The inclination angle is fixed at i = 30° and MBH = 10 M⊙.
In the text

	Fig. 5. Top panel: total PSD in 0.3–0.7 keV band varying with γ when i = 30°, rtrc = 10rg, rsh = 5rg, ϵ = 0.1 and (H/R)2α = 0.01. Bottom panel: total PSD in 0.3–0.7 keV band varying the geometry of the system with γ = 1.5. All dotted lines represents the corresponding cases with reverberation excluded. The reduction of power at high frequencies is produced by reverberation.
In the text

Fig. 6. Data, model and residuals from fitting the PSD model to 0.3–0.7 keV (blue) and 0.7–1.5 keV (red) PSD data with (H/R)2α fixed at 0.005 and γ fixed at 0 as obtained by the global fits. The PSD displayed here are Poisson noise-subtracted and normalized in terms of squared fractional rms (Miyamoto et al. 1991). Data points consistent with zero (or negative) power have been excluded. The best-fit model parameters for each observation are listed in Table 2.

In the text

	Fig. 7. Comparison of the truncation radius vs. Eddington-scaled luminosity for GX 339–4 inferred by this work (red diamonds), De Marco et al. (2017) via approximate reverberation measurement (circles), and Mahmoud et al. (2019) via spectral-timing fits (crosses).
In the text

	Fig. 8. Comparative size of inner hot-flow zones between different observations (O1–O6) implied from the best-fit parameters. Error bars are small and are excluded in this plot.
In the text