Multi-timescale reverberation mapping of Mrk 335
MMNRAS , 1–14 (2019) Preprint 10 September 2020 Compiled using MNRAS L A TEX style file v3.0
Multi-timescale reverberation mapping of Mrk 335
Guglielmo Mastroserio, , (cid:63) Adam Ingram, & Michiel van der Klis Cahill Center for Astronomy and Astrophysics, California Institute of Technology, 1200 California Boulevard, Pasadena, CA 91125, USA Astronomical Institute Anton Pannekoek, University of Amsterdam, Science Park 904, NL-1098 XH Amsterdam, Netherlands Department of Physics, Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Time lags due to X-ray reverberation have been detected in several Seyfert galaxies.The different travel time between reflected and directly observed rays naturally causesthis type of lag, which depends directly on the light-crossing timescale of the systemand hence scales with the mass of the central black hole. Featureless ‘hard lags’ notassociated with reverberation, and often interpreted as propagating mass accretionrate fluctuations, dominate the longer timescale variability. Here we fit our reltrans model simultaneously to the time-averaged energy spectrum and the lag-energy spectraof the Seyfert galaxy Mrk 335 over two timescales (Fourier frequency ranges). Wemodel the hard lags as fluctuations in the slope and strength of the illuminatingspectrum, and self-consistently account for the effects that these fluctuations haveon the reverberation lags. The resulting mass estimate is . + . − . × M (cid:12) , which issignificantly lower than the mass measured with the optical reverberation mappingtechnique ( - million M (cid:12) ). When we add the correlated variability amplitudes tothe time lags by fitting the full complex cross-spectra, the model is unable to describethe characteristic reverberation Fe K α line and cannot constrain the black hole mass.This may be due to the assumption that the direct radiation is emitted by a point-likesource. Key words:
Black hole physics – Relativistic processes – X-rays: galaxies – Galaxies:nuclei
Active Galactic Nuclei (AGN) are thought to be powered bythe accretion of matter onto supermassive black holes ( − M (cid:12) ). The gas forms an accretion disc and emits a multi-temperature blackbody spectrum which peaks in the UVband (Shields 1978; Malkan 1983). Some of these photonsact as seed photons for inverse Compton up-scattering in ahot electron ‘corona’ (Eardley et al. 1975; Thorne & Price1975). The resulting direct Comptonised emission is oftensimply modelled as a power-law with a high energy cut-off, where the power-law index and the cut-off energy arerelated to the optical depth and the electron temperature,respectively.Part of this radiation illuminates the disc and is re-emitted, producing the reflection component in the spec-trum (e.g., Lightman & White 1988). The shape of the over-all time-averaged energy spectrum depends on the propertiesof both the corona, such as optical depth and temperature,and the accretion disc, including both its physical character-istics (e.g., ionisation and iron abundance) and its geometry (cid:63)
E-mail: [email protected] (e.g., inclination and radius of the inner edge). One of themost prominent features in the reflection spectrum is an iron K α line emitted at . keV and broadened by the effects oforbital motion in the disc and relativistic redshift (Fabianet al. 1989). The black hole drags space-time around it fur-ther modifying the emission line profile (Fabian et al. 2000),making it possible to measure the spin of the black hole (seeReynolds 2014 for a review).Although this approach of time-averaged spectral fittingled to very interesting results, there are several aspects thathave remained unclear, in particular, the exact geometry ofthe system. For example, the extent of the corona is un-der debate. It might be either a compact cloud of gas thatcould reasonably (compared to the accretion disc dimen-sions) be modelled as a point source (Haardt & Maraschi1991), or a more extended structure, with various geome-tries being considered for it (e.g. Eardley et al. 1975; Haardt& Maraschi 1993). One of the reasons for this ambiguity isthat in time averaged spectral fitting, degeneracies amongthe model parameters can often not be avoided. In orderto help break such degeneracies, it is possible to addition-ally study the time variability of the spectrum. One effec-tive way to study the spectral variability is to measure the © a r X i v : . [ a s t r o - ph . H E ] S e p G. Mastroserio et al. time lags observed in both supermassive and stellar massblack holes between the variability in different photon energybands. Two types of lag are observed in AGN, the ‘hard lag’that is intrinsic to the direct emission (intrinsic lags) andthe ‘soft lag’ that is due to the differences in light cross-ing time from corona to observer between the direct andthe reflected emission (the reverberation lags). The intrin-sic lags are thought to be generated by mass accretion ratefluctuations propagating through the accretion disc towardsthe black hole (e.g. Lyubarskii 1997; Ar´evalo & Uttley 2006;Ingram & van der Klis 2013). The timescales of the fluctua-tions depend on the viscosity of the gas, thus long timescalesarise in the outer part of the accretion flow and shorttimescales closer to the black hole. The fluctuations reachthe corona causing variations in the electron temperature.At each radius the variability is due to the product of localfluctuations and the generally slower fluctuations that prop-agated in from larger radii. Conversely, the reverberationlags are due to the different path lengths of the variabilitysignals to reach the observer. Hence, the variability in theenergy band dominated by the reflected emission lags behindthe variability of the direct, inverse Compton emission.The hard lags dominate at long timescales, and hencein many sources the lags as a function of Fourier frequencycomputed between (cid:39) . − . keV (soft energy band) and (cid:39) . − . (hard energy band) are positive (i.e., hard fluxlags soft flux) at Fourier frequencies ν ≤ M (cid:12) / M Hz.The lag spectrum as a function of photon energy calculatedin this frequency range is featureless and lag depends ap-proximately linearly on log E (e.g. Papadakis et al. 2001;McHardy et al. 2004). At higher frequencies the hard in-trinsic lags are smaller, and soft (negative) lags have beendetected between the same two energy bands (e.g. Fabianet al. 2009). The lag energy spectrum in this higher fre-quency range is indicative of reverberation, exhibiting a neg-ative slope for E (cid:46) keV (hence the soft lag between thetwo broad energy bands), with a broad iron line feature inthe − keV range (e.g. Kara et al. 2016) and a Comptonhump in the − keV range (Zoghbi et al. 2014). De-tecting this Fe line reverberation feature in the lag energyspectrum is easier in AGN than in stellar mass black holesbecause they are much bigger systems and so the timescalesinvolved are consequently longer. However, because of theirdistance, AGN are usually fainter than stellar mass blackholes, so they are often characterized by a worse signal tonoise. For that reason, we need to use as wide as possible afrequency range to calculate the lag energy spectrum, andanalyse it together with the time-averaged energy spectrumin order to constrain the system parameters.Looking only at the time lags means considering only alimited part of the information in the data, as the correlatedvariability amplitude as a function of energy and frequencyis ignored in this approach. Lags and amplitudes have beenconsidered together in only a few cases (e.g. Uttley et al.2011; Kara et al. 2013a; Rapisarda et al. 2016), however,joint modelling of both has not yet been performed for AGN.So, progress could be made by jointly modelling the time-averaged and variability (time lag and correlated amplitude)spectra taking into account both the effects of mass accretionrate propagating fluctuations and reverberation lags. Thebest way to consider time lags and correlated amplitudestogether is by explicitly modeling the real and imaginary parts of the complex cross spectra (van der Klis et al. 1987;Rapisarda et al. 2016; Ingram et al. 2016; Mastroserio et al.2018), and that is the approach taken in this paper.We have developed the fully relativistic reverberationmapping model reltrans (Ingram et al. 2019) that, usinga transfer function formalism (e.g. Campana & Stella 1995;Reynolds et al. 1999; Cackett et al. 2014), computes the com-plex cross-spectrum for a range of Fourier frequencies andthe time-averaged energy spectrum in a prescription basedon a ‘lamppost’ coronal geometry (e.g. Matt et al. 1991),assuming isotropic emission and a flat thin accretion disc.In the configuration of reltrans that we use, the disc ra-dial ionisation profile is self-consistently calculated from theradial density profile corresponding to a Shakura & Sunyaev(1973) zone A accretion disc. However, the rest-frame reflec-tion spectrum is computed using the model xillver (Garc´ıaet al. 2013), which employs a fixed electron number densityof cm − . Future versions of reltrans will include the vari-able electron density version of xillver , xillverD (Garc´ıaet al. 2016), which will enable us to include self-consistentlya radial density profile as well as a radial ionization pro-file. The reltrans model computes the reverberation lagsand accounts for the different light crossing times of thephotons using proper relativistic ray-tracing both betweencorona and each point on the disc and from each point on thedisc to the observer. We use here a further improved versionof reltrans which accounts for the hard lags through thepivoting of the continuum spectrum by allowing fluctuationsin the index and normalisation of the power law (Mastrose-rio et al. 2019). This serves as a proxy for the cooling andheating of the corona by the fluctuations in the rate of seedphotons coming from the disc and the mass accretion ratefluctuations propagating into the corona. Each patch of thedisc sees a different hardness of the incident emission comingfrom the corona due to the effects of orbital motion and red-shift which are properly evaluated taking into account thedifferent paths of the photons (Mastroserio et al. 2018). Thiscauses non-linear effects in the fluctuations of the reflectedemission (i.e., not just the strength but also the shape ofthe reflection spectrum fluctuates) which are taken into ac-count in the model. reltrans has been successfully testedon Cygnus X-1 (Mastroserio et al. 2019) and previously usedfor a proof of principle of the method (fitting only the re-verberation lags) in Mrk 335 (Ingram et al. 2019).In this paper we present a full exploration of the modelwith an XMM-Newton observations of Mrk 335. The Seyfert1 galaxy Mrk 335 has been extensively studied both in termsof time-averaged spectrum (e.g. Wilkins & Gallo 2015; Keek& Ballantyne 2016) and lag spectrum, which shows bothhard and soft lags (e.g., Kara et al. 2013b; De Marco et al.2013a; Chainakun et al. 2016). The source has been observedat different flux levels (defining ’flux epochs’) (Wilkins &Fabian 2013) with XMM-Newton from 2006 to 2009. Re-verberation lags are detected only in the highest-flux epoch(Kara et al. 2013b) although the reflection component in thespectrum is weaker then than in the low flux epochs (Keek& Ballantyne 2016). We therefore select the high flux epochand use reltrans to simultaneously fit the time-averagedenergy spectrum and the complex cross-spectra as a functionof energy in multiple Fourier frequency ranges. MNRAS000
E-mail: [email protected] (e.g., inclination and radius of the inner edge). One of themost prominent features in the reflection spectrum is an iron K α line emitted at . keV and broadened by the effects oforbital motion in the disc and relativistic redshift (Fabianet al. 1989). The black hole drags space-time around it fur-ther modifying the emission line profile (Fabian et al. 2000),making it possible to measure the spin of the black hole (seeReynolds 2014 for a review).Although this approach of time-averaged spectral fittingled to very interesting results, there are several aspects thathave remained unclear, in particular, the exact geometry ofthe system. For example, the extent of the corona is un-der debate. It might be either a compact cloud of gas thatcould reasonably (compared to the accretion disc dimen-sions) be modelled as a point source (Haardt & Maraschi1991), or a more extended structure, with various geome-tries being considered for it (e.g. Eardley et al. 1975; Haardt& Maraschi 1993). One of the reasons for this ambiguity isthat in time averaged spectral fitting, degeneracies amongthe model parameters can often not be avoided. In orderto help break such degeneracies, it is possible to addition-ally study the time variability of the spectrum. One effec-tive way to study the spectral variability is to measure the © a r X i v : . [ a s t r o - ph . H E ] S e p G. Mastroserio et al. time lags observed in both supermassive and stellar massblack holes between the variability in different photon energybands. Two types of lag are observed in AGN, the ‘hard lag’that is intrinsic to the direct emission (intrinsic lags) andthe ‘soft lag’ that is due to the differences in light cross-ing time from corona to observer between the direct andthe reflected emission (the reverberation lags). The intrin-sic lags are thought to be generated by mass accretion ratefluctuations propagating through the accretion disc towardsthe black hole (e.g. Lyubarskii 1997; Ar´evalo & Uttley 2006;Ingram & van der Klis 2013). The timescales of the fluctua-tions depend on the viscosity of the gas, thus long timescalesarise in the outer part of the accretion flow and shorttimescales closer to the black hole. The fluctuations reachthe corona causing variations in the electron temperature.At each radius the variability is due to the product of localfluctuations and the generally slower fluctuations that prop-agated in from larger radii. Conversely, the reverberationlags are due to the different path lengths of the variabilitysignals to reach the observer. Hence, the variability in theenergy band dominated by the reflected emission lags behindthe variability of the direct, inverse Compton emission.The hard lags dominate at long timescales, and hencein many sources the lags as a function of Fourier frequencycomputed between (cid:39) . − . keV (soft energy band) and (cid:39) . − . (hard energy band) are positive (i.e., hard fluxlags soft flux) at Fourier frequencies ν ≤ M (cid:12) / M Hz.The lag spectrum as a function of photon energy calculatedin this frequency range is featureless and lag depends ap-proximately linearly on log E (e.g. Papadakis et al. 2001;McHardy et al. 2004). At higher frequencies the hard in-trinsic lags are smaller, and soft (negative) lags have beendetected between the same two energy bands (e.g. Fabianet al. 2009). The lag energy spectrum in this higher fre-quency range is indicative of reverberation, exhibiting a neg-ative slope for E (cid:46) keV (hence the soft lag between thetwo broad energy bands), with a broad iron line feature inthe − keV range (e.g. Kara et al. 2016) and a Comptonhump in the − keV range (Zoghbi et al. 2014). De-tecting this Fe line reverberation feature in the lag energyspectrum is easier in AGN than in stellar mass black holesbecause they are much bigger systems and so the timescalesinvolved are consequently longer. However, because of theirdistance, AGN are usually fainter than stellar mass blackholes, so they are often characterized by a worse signal tonoise. For that reason, we need to use as wide as possible afrequency range to calculate the lag energy spectrum, andanalyse it together with the time-averaged energy spectrumin order to constrain the system parameters.Looking only at the time lags means considering only alimited part of the information in the data, as the correlatedvariability amplitude as a function of energy and frequencyis ignored in this approach. Lags and amplitudes have beenconsidered together in only a few cases (e.g. Uttley et al.2011; Kara et al. 2013a; Rapisarda et al. 2016), however,joint modelling of both has not yet been performed for AGN.So, progress could be made by jointly modelling the time-averaged and variability (time lag and correlated amplitude)spectra taking into account both the effects of mass accretionrate propagating fluctuations and reverberation lags. Thebest way to consider time lags and correlated amplitudestogether is by explicitly modeling the real and imaginary parts of the complex cross spectra (van der Klis et al. 1987;Rapisarda et al. 2016; Ingram et al. 2016; Mastroserio et al.2018), and that is the approach taken in this paper.We have developed the fully relativistic reverberationmapping model reltrans (Ingram et al. 2019) that, usinga transfer function formalism (e.g. Campana & Stella 1995;Reynolds et al. 1999; Cackett et al. 2014), computes the com-plex cross-spectrum for a range of Fourier frequencies andthe time-averaged energy spectrum in a prescription basedon a ‘lamppost’ coronal geometry (e.g. Matt et al. 1991),assuming isotropic emission and a flat thin accretion disc.In the configuration of reltrans that we use, the disc ra-dial ionisation profile is self-consistently calculated from theradial density profile corresponding to a Shakura & Sunyaev(1973) zone A accretion disc. However, the rest-frame reflec-tion spectrum is computed using the model xillver (Garc´ıaet al. 2013), which employs a fixed electron number densityof cm − . Future versions of reltrans will include the vari-able electron density version of xillver , xillverD (Garc´ıaet al. 2016), which will enable us to include self-consistentlya radial density profile as well as a radial ionization pro-file. The reltrans model computes the reverberation lagsand accounts for the different light crossing times of thephotons using proper relativistic ray-tracing both betweencorona and each point on the disc and from each point on thedisc to the observer. We use here a further improved versionof reltrans which accounts for the hard lags through thepivoting of the continuum spectrum by allowing fluctuationsin the index and normalisation of the power law (Mastrose-rio et al. 2019). This serves as a proxy for the cooling andheating of the corona by the fluctuations in the rate of seedphotons coming from the disc and the mass accretion ratefluctuations propagating into the corona. Each patch of thedisc sees a different hardness of the incident emission comingfrom the corona due to the effects of orbital motion and red-shift which are properly evaluated taking into account thedifferent paths of the photons (Mastroserio et al. 2018). Thiscauses non-linear effects in the fluctuations of the reflectedemission (i.e., not just the strength but also the shape ofthe reflection spectrum fluctuates) which are taken into ac-count in the model. reltrans has been successfully testedon Cygnus X-1 (Mastroserio et al. 2019) and previously usedfor a proof of principle of the method (fitting only the re-verberation lags) in Mrk 335 (Ingram et al. 2019).In this paper we present a full exploration of the modelwith an XMM-Newton observations of Mrk 335. The Seyfert1 galaxy Mrk 335 has been extensively studied both in termsof time-averaged spectrum (e.g. Wilkins & Gallo 2015; Keek& Ballantyne 2016) and lag spectrum, which shows bothhard and soft lags (e.g., Kara et al. 2013b; De Marco et al.2013a; Chainakun et al. 2016). The source has been observedat different flux levels (defining ’flux epochs’) (Wilkins &Fabian 2013) with XMM-Newton from 2006 to 2009. Re-verberation lags are detected only in the highest-flux epoch(Kara et al. 2013b) although the reflection component in thespectrum is weaker then than in the low flux epochs (Keek& Ballantyne 2016). We therefore select the high flux epochand use reltrans to simultaneously fit the time-averagedenergy spectrum and the complex cross-spectra as a functionof energy in multiple Fourier frequency ranges. MNRAS000 , 1–14 (2019) everberation lag in Mrk 335 We analyse a ks observation of Mrk 335 performed with
XMM-Newton in 2006 (obs ID 0306870101). During the ob-servation the source was in the high-flux epoch (Wilkins &Gallo 2015) with a . − keV flux of . × − erg cm − s − (Grupe et al. 2007). We follow the data reduction pro-cedure described in Ingram et al. (2019), considering onlythe EPIC-pn data and discarding the EPIC-MOS data (fol-lowing Kara et al. 2013b). We use the Science Analysis Sys-tem (SAS) v11.0.0 to extract the signal from a circular re-gion with arcsec radius centred on the maximum of thesource emission. We apply the filters PATTERN ≤ andFLAG == and discard background flares at the beginningand at the end of the observation (considering only times252709714 to 252829414 seconds telescope time). We extractthe time-averaged source and background energy spectrumand also extract light curves with second binning from different energy bands with the same energy resolutionas used by Kara et al. (2013b). The SAS task epiclccorr performs various corrections including subtracting the back-ground signal, which is extracted from a circular region of arcsec radius at some distance from the source.We compute the cross-spectrum between the light curvefor each energy band and the reference light curve which isthe sum of all the light curves except the subject one (Utt-ley et al. 2014), and average it into broad frequency ranges.We consider the full length of the light curves in order toprobe the same frequency ranges as Kara et al. (2013b) andDe Marco et al. (2013a). Before Fourier transforming thelight curves, we fill in small gaps (less than seconds) byrandomly extracting the missing count rates from a Pois-son distribution centered on the value interpolated betweenthe previous and the next bin of the gap. We divide thecross-spectrum into two frequency ranges: . − . mHz, . − . mHz. These two frequency ranges are roughly inthe regimes dominated by the hard and soft lags, respec-tively, as found by De Marco et al. (2013a). We find thatthe intrinsic coherence (Vaughan & Nowak 1997) is consis-tent with unity in these two frequency ranges but drops offsignificantly at higher frequencies, as is seen in other AGN(Zoghbi et al. 2010; Uttley et al. 2014). Since the model as-sumes unity coherence, we ignore the higher frequencies thatdisplay low coherence. Since the Fourier frequency resolutionof the cross-spectrum is . × − Hz, the first frequencyrange contains Fourier frequencies, and the second .The cross-spectral amplitudes averaged over these two fre-quency ranges are therefore sufficiently close to Gaussian-distributed that we can use the χ statistic for the purposesof fitting models (Nowak et al. 1999).Epitropakis & Papadakis (2016) pointed out that con-sidering the full length of the light curve and averaging thecross-spectrum into broad frequency bins (smoothing) canintroduce a bias to the time-lag estimates. Specifically, thetime lag estimated for a given frequency bin is not equalto the expectation value of the true time lag at the centralfrequency of the considered frequency bin. This occurs be-cause the amplitude and phase of the cross-spectrum candepend quite steeply on frequency, and therefore the cross-spectrum can change quite substantially and non-linearlywithin a broad frequency bin. The improved version of rel-trans that we present in this paper (see Section 3) models this bias by first calculating the predicted frequency depen-dent cross-spectrum and then averaging over the frequencyrange corresponding to the observational data. This can beunderstood as the model including exactly the same biasesas the data do, leading to the best fitting model parametersthemselves being unbiased . The old version of the model(Ingram et al. 2019 ) onlyapproximately followed the de-scribed procedure. These approximations are good for thelow frequencies (in terms of c / R g ) probed by fits to X-raybinary data (e.g. Mastroserio et al. 2019), but break downfor the higher frequencies (again in terms of c / R g ) probedby fits to AGN data. We therefore first address the inaccu-racies in the old model (Section 3) before fitting to the Mrk335 data (Section 4). In this Section, we describe the new version of the reltrans model used in this paper. A future paper will detail furtherimprovements in addition to those described here, and willbe accompanied by the public release of a new model version(version 2). The initial release of the model was described inIngram et al. (2019). Subsequently, Mastroserio et al. (2019)included the non-linear effect of variations in the power-lawindex of the illuminating spectrum, which can reproduce theobserved hard lags. Here, we build on the Mastroserio et al.(2019) model by relaxing some assumptions that caused in-accuracies in the high frequency ranges (in terms of c / R g )probed by AGN. Averaging the transfer functions over thefrequency range (the method used in the old version of themodel) is correct only when the transfer function (modulusand phase) and the variability amplitude of the continuumvariations do not depend strongly on frequency. This is notthe case for the higher frequencies at which the reverberationlags dominate. Thus the improved model calculates the fre-quency dependent cross-spectrum before then averaging overthe specified frequency range. We present below the math-ematical details of our improvements to the model, whichinclude a slight change in the model parameters specifyingthe spectral pivoting that makes them easier to interpretphysically. We start by representing the time-dependent spec-trum directly observed from the corona as D ( E , t ) = A ( t ) (cid:96) g Γ ( t ) soz E − Γ ( t ) e − E / E cut , obs , where (cid:96) and E cut , obs are respec-tively the lensing factor and observed high energy cut-off asdefined in Ingram et al. (2019) and g soz ≡ g so /( + z ) , where g so is the blueshift experienced by a photon traveling fromthe corona to an observer at infinity in an asymptoticallyflat and static spacetime and z is the cosmological redshiftexperienced by the photon. The normalisation A ( t ) varies We note that, although this fit procedure is unbiased, it does ig-nore diagnostic information. As an extreme example, if the phaseturns by ◦ in a frequency bin, the cross-spectrum of data andmodel can be zero for that bin. In such an example, the fit wouldnot be biased, but it would be degenerate. The old version of the code is publicly available at
MNRAS , 1–14 (2019)
G. Mastroserio et al. around a mean value A and the power-law index Γ ( t ) varies(with small variability amplitude) around a mean value Γ .The total spectrum is the sum of the directly observed emis-sion and the reflected emission, which we calculate using thetransfer function formalism extensively described in Ingramet al. (2019) and Mastroserio et al. (2019).The Fourier transform (FT) of the full spectrum is S ( E , ν ) = e − τ ( E ) (cid:26) A ( ν ) [ D ( E ) + W ( E , ν )] + A Γ ( ν ) [ ln ( g soz / E ) D ( E ) + W ( E , ν ) + W ( E , ν )] (cid:27) . (1)This expression is very similar to Equation 8 in Mastroserioet al. (2019), except here we explicitly include interstellarabsorption with the e − τ ( E ) factor (this was left as implicitin Mastroserio et al. 2019), and Mastroserio et al. (2019)defined a new term B ( ν ) = − A Γ ( ν ) , which we now refrainfrom doing in order to make the model parameters easier tointerpret. Here, D ( E ) ≡ (cid:96) g Γ soz E − Γ e − E / E cut , obs , and W ( E , ν ) , W ( E , ν ) and W ( E , ν ) are transfer functions that are writtenout in full in Appendix A. W represents the response ofthe reflection spectrum to fluctuations in the normalisationof the illuminating spectrum and W and W represent theresponse to fluctuations in the power-law index ( W accountsfor changes to the emissivity profile and W accounts forchanges to the restframe reflection spectrum).Defining γ ( ν ) e i ∆ BA ( ν ) ≡ A Γ ( ν )/ A ( ν ) , Equation (1) be-comes S ( E , ν ) = A ( ν ) e − τ ( E ) (cid:26) D ( E ) + W ( E , ν ) + γ ( ν ) e i ∆ BA ( ν ) [ ln ( g soz / E ) D ( E ) + W ( E , ν ) + W ( E , ν )] (cid:27) . (2)Here γ ( ν ) = A | Γ ( ν )|/| A ( ν )| is the amplitude of fluctuations in Γ divided by the fractional variability amplitude of fluctua-tions in A , and ∆ BA ( ν ) is the phase lag between fluctuationsin Γ and those in A (positive means Γ lags A ). Since it iseasier to interpret, we adopt ∆ BA ( ν ) in place of φ B ( ν ) , whichwe used in Mastroserio et al. (2018) and Mastroserio et al.(2019).For the case of the observational data, we take theFourier transform (FT) of the count rate in the subject bandand then cross with the FT of the reference band countrate. The model has to follow the same procedure. In pre-vious versions of the model, we averaged the transfer func-tions over the frequency range specified as model parameters α ( ν ) = | A ( ν )| , φ A ( ν ) = arg [ A ( ν )] and φ B ( ν ) averaged over thesame frequency range. This procedure is inaccurate if anyof these three quantities, or any of the transfer functions,change significantly across the frequency range. We thereforenow first calculate the frequency dependent cross-spectrumand then average over the input frequency range.The cross-spectrum is G ( E , ν ) = S ( E , ν ) F ∗ r ( ν ) , where F r ( ν ) is the FT of the reference band count rate time series.The reference band count rate is the sum of the count ratesfrom each of the energy channels I min to I max that make upthe reference band. To calculate the reference band FT, wemust therefore fold S ( E , ν ) around the instrument responseto get S ( I , ν ) and then sum over the reference band channels. The cross-spectrum is therefore G ( E , ν ) = | A ( ν )| S raw ( E , ν ) F ∗ r , raw ( ν ) = | A ( ν )| G raw ( E , ν ) , (3)where S raw ( E , ν ) ≡ S ( E , ν )/ A ( ν ) and F r , raw ( ν ) ≡ F r ( ν )/ A ( ν ) .Note that G ( E , ν ) can be very simply folded around the in-strument response to get G ( I , ν ) .Once we have calculated the cross-spectrum for thedata, we average over a frequency range ν lo to ν hi centeredat ν c = ( ν hi + ν lo )/ . Again the model has to do the same toget (cid:104) G ( E , ν c )(cid:105) = ∫ ν hi ν lo G ( E , ν ) d νν hi − ν lo = ∫ ν hi ν lo | A ( ν )| G raw ( E , ν ) d νν hi − ν lo . (4)We therefore need to assume a form for | A ( ν )| . We rep-resent this as a power-law function of frequency, | A ( ν )| = α ( ν c )( ν / ν c ) − k , giving (cid:104) G ( E , ν c )(cid:105) = α ( ν c ) ∫ ν hi ν lo ( ν / ν c ) − k G raw ( E , ν ) d νν hi − ν lo . (5)The full band power spectrum is P ( ν ) ∼ | A ( ν )| , and obser-vationally tends to consist roughly of a twice broken power-law: k ≈ at very low frequencies ( ν (cid:46) . Hz in X-raybinaries), k ≈ at intermediate frequencies and k ≈ athigh frequencies ( ν (cid:38) Hz in X-ray binaries). The model isinsensitive to k in the low and intermediate frequency rangesfor which k ≈ and k ≈ . This is because for such values of ν c , / ν c is much longer than the light-crossing timescale ofthe disc, and therefore the transfer functions do not changemuch from ν = ν lo to ν = ν hi . The model is sensitive to k for the highest frequencies, because here the cross-spectrum(modulus and phase) can vary steeply with ν . We thereforeadopt k = .The output model for the cross-spectrum is therefore (cid:104) G ( E , ν c )(cid:105) (real and imaginary parts). The model lag-energyspectrum is given by t lag ( I , ν c ) = arg [(cid:104) G ( I , ν c )(cid:105)]/( πν c ) , whichexactly mirrors our procedure to constrain the lag-energyspectrum from the observed cross-spectrum. The model pa-rameters for the hard lags are γ ( ν c ) and ∆ BA ( ν c ) . Note thatthe lag-energy spectrum model is not sensitive to the param-eter α ( ν c ) , but the model for the real and imaginary partsof the cross-spectrum is because α ( ν c ) affects the modulusof the cross-spectrum. Fitting the − keV time-averaged spectrum with an ab-sorbed power-law model reveals evident residuals in the iron K α line energy range (see Fig. 1), consistent with previousanalyses (Keek & Ballantyne 2016; Wilkins & Gallo 2015).We also see the emission line at . keV in the rest frameof the host galaxy (the cosmological redshift to Mrk 335 is z = . ; Huchra et al. 1999) that was reported by Keek& Ballantyne (2016). Throughout this paper, we model thisline in the time-averaged spectrum with a narrow Gaussianfunction with fixed width ( − ) and free centroid (allowedto vary between . and . ). Including a xillver (Garc´ıaet al. 2013) component in the model to account for distantreflection leads to a fit with χ / d . o . f . of . / (Model MNRAS000
G. Mastroserio et al. around a mean value A and the power-law index Γ ( t ) varies(with small variability amplitude) around a mean value Γ .The total spectrum is the sum of the directly observed emis-sion and the reflected emission, which we calculate using thetransfer function formalism extensively described in Ingramet al. (2019) and Mastroserio et al. (2019).The Fourier transform (FT) of the full spectrum is S ( E , ν ) = e − τ ( E ) (cid:26) A ( ν ) [ D ( E ) + W ( E , ν )] + A Γ ( ν ) [ ln ( g soz / E ) D ( E ) + W ( E , ν ) + W ( E , ν )] (cid:27) . (1)This expression is very similar to Equation 8 in Mastroserioet al. (2019), except here we explicitly include interstellarabsorption with the e − τ ( E ) factor (this was left as implicitin Mastroserio et al. 2019), and Mastroserio et al. (2019)defined a new term B ( ν ) = − A Γ ( ν ) , which we now refrainfrom doing in order to make the model parameters easier tointerpret. Here, D ( E ) ≡ (cid:96) g Γ soz E − Γ e − E / E cut , obs , and W ( E , ν ) , W ( E , ν ) and W ( E , ν ) are transfer functions that are writtenout in full in Appendix A. W represents the response ofthe reflection spectrum to fluctuations in the normalisationof the illuminating spectrum and W and W represent theresponse to fluctuations in the power-law index ( W accountsfor changes to the emissivity profile and W accounts forchanges to the restframe reflection spectrum).Defining γ ( ν ) e i ∆ BA ( ν ) ≡ A Γ ( ν )/ A ( ν ) , Equation (1) be-comes S ( E , ν ) = A ( ν ) e − τ ( E ) (cid:26) D ( E ) + W ( E , ν ) + γ ( ν ) e i ∆ BA ( ν ) [ ln ( g soz / E ) D ( E ) + W ( E , ν ) + W ( E , ν )] (cid:27) . (2)Here γ ( ν ) = A | Γ ( ν )|/| A ( ν )| is the amplitude of fluctuations in Γ divided by the fractional variability amplitude of fluctua-tions in A , and ∆ BA ( ν ) is the phase lag between fluctuationsin Γ and those in A (positive means Γ lags A ). Since it iseasier to interpret, we adopt ∆ BA ( ν ) in place of φ B ( ν ) , whichwe used in Mastroserio et al. (2018) and Mastroserio et al.(2019).For the case of the observational data, we take theFourier transform (FT) of the count rate in the subject bandand then cross with the FT of the reference band countrate. The model has to follow the same procedure. In pre-vious versions of the model, we averaged the transfer func-tions over the frequency range specified as model parameters α ( ν ) = | A ( ν )| , φ A ( ν ) = arg [ A ( ν )] and φ B ( ν ) averaged over thesame frequency range. This procedure is inaccurate if anyof these three quantities, or any of the transfer functions,change significantly across the frequency range. We thereforenow first calculate the frequency dependent cross-spectrumand then average over the input frequency range.The cross-spectrum is G ( E , ν ) = S ( E , ν ) F ∗ r ( ν ) , where F r ( ν ) is the FT of the reference band count rate time series.The reference band count rate is the sum of the count ratesfrom each of the energy channels I min to I max that make upthe reference band. To calculate the reference band FT, wemust therefore fold S ( E , ν ) around the instrument responseto get S ( I , ν ) and then sum over the reference band channels. The cross-spectrum is therefore G ( E , ν ) = | A ( ν )| S raw ( E , ν ) F ∗ r , raw ( ν ) = | A ( ν )| G raw ( E , ν ) , (3)where S raw ( E , ν ) ≡ S ( E , ν )/ A ( ν ) and F r , raw ( ν ) ≡ F r ( ν )/ A ( ν ) .Note that G ( E , ν ) can be very simply folded around the in-strument response to get G ( I , ν ) .Once we have calculated the cross-spectrum for thedata, we average over a frequency range ν lo to ν hi centeredat ν c = ( ν hi + ν lo )/ . Again the model has to do the same toget (cid:104) G ( E , ν c )(cid:105) = ∫ ν hi ν lo G ( E , ν ) d νν hi − ν lo = ∫ ν hi ν lo | A ( ν )| G raw ( E , ν ) d νν hi − ν lo . (4)We therefore need to assume a form for | A ( ν )| . We rep-resent this as a power-law function of frequency, | A ( ν )| = α ( ν c )( ν / ν c ) − k , giving (cid:104) G ( E , ν c )(cid:105) = α ( ν c ) ∫ ν hi ν lo ( ν / ν c ) − k G raw ( E , ν ) d νν hi − ν lo . (5)The full band power spectrum is P ( ν ) ∼ | A ( ν )| , and obser-vationally tends to consist roughly of a twice broken power-law: k ≈ at very low frequencies ( ν (cid:46) . Hz in X-raybinaries), k ≈ at intermediate frequencies and k ≈ athigh frequencies ( ν (cid:38) Hz in X-ray binaries). The model isinsensitive to k in the low and intermediate frequency rangesfor which k ≈ and k ≈ . This is because for such values of ν c , / ν c is much longer than the light-crossing timescale ofthe disc, and therefore the transfer functions do not changemuch from ν = ν lo to ν = ν hi . The model is sensitive to k for the highest frequencies, because here the cross-spectrum(modulus and phase) can vary steeply with ν . We thereforeadopt k = .The output model for the cross-spectrum is therefore (cid:104) G ( E , ν c )(cid:105) (real and imaginary parts). The model lag-energyspectrum is given by t lag ( I , ν c ) = arg [(cid:104) G ( I , ν c )(cid:105)]/( πν c ) , whichexactly mirrors our procedure to constrain the lag-energyspectrum from the observed cross-spectrum. The model pa-rameters for the hard lags are γ ( ν c ) and ∆ BA ( ν c ) . Note thatthe lag-energy spectrum model is not sensitive to the param-eter α ( ν c ) , but the model for the real and imaginary partsof the cross-spectrum is because α ( ν c ) affects the modulusof the cross-spectrum. Fitting the − keV time-averaged spectrum with an ab-sorbed power-law model reveals evident residuals in the iron K α line energy range (see Fig. 1), consistent with previousanalyses (Keek & Ballantyne 2016; Wilkins & Gallo 2015).We also see the emission line at . keV in the rest frameof the host galaxy (the cosmological redshift to Mrk 335 is z = . ; Huchra et al. 1999) that was reported by Keek& Ballantyne (2016). Throughout this paper, we model thisline in the time-averaged spectrum with a narrow Gaussianfunction with fixed width ( − ) and free centroid (allowedto vary between . and . ). Including a xillver (Garc´ıaet al. 2013) component in the model to account for distantreflection leads to a fit with χ / d . o . f . of . / (Model MNRAS000 , 1–14 (2019) everberation lag in Mrk 335 [A] in Table 1). This poor fit implies that a relativistic reflec-tion component is also required, for which we use the model reltrans (Ingram et al. 2019). The full expression for ourmodel (Model [B]) is tbabs × ( xillver + reltrans + zagauss ) , (6)where the direct continuum emission is included in the rel-trans model as an exponentially cut-off power-law. We fixthe hydrogen column density to n H = . × cm − follow-ing Kalberla et al. (2005) and assume the relative elementalabundances of Wilms et al. (2000). We fix the ionizationparameter of the distant reflector to log ξ = . In contrast,we model the radial profile of the disc ionization parameterusing a self-consistent calculation of the irradiating flux ina lamppost geometry and assuming the radial density pro-file that corresponds to ‘zone A’ in the Shakura & Sunyaev(1973) disc model . We use 10 radial zones to calculate theionisation profile: one zone between r in and ( / ) r in andthe rest logarithmically separated. r = ( / ) r in is approxi-mately the radius whereby the ionisation parameter peaks ,and the value of log ξ at this radius is a model parameter thatwe leave free in the fit. We choose to only model the . − keV energy range of the spectrum because the . − . keVenergy range can only be adequately described by includ-ing a complicated 3 layer warm absorber model (Longinottiet al. 2013), which is beyond the scope of this paper.Fig 2 shows the best fitting model together with the un-folded data (upper panel). The residuals (bottom panel) donot present any significant structure and the fit has χ / d . o . f . of / . Here, and throughout this paper, we have fixedthe high energy cut-off, as this quantity cannot be con-strained in the < keV energy range of XMM-Newton ,to E cut = keV (in the observer frame), and we have alsofixed the inclination angle to i = ◦ (following for both ofthese parameters Keek & Ballantyne 2016). We fix the di-mensionless black hole spin parameter to a = . and allowthe disc inner radius to be a free parameter (constrained tobe larger than the ISCO radius).The best fitting model parameters are listed in Table 1.We obtain a larger disc inner radius than the fit of Keek& Ballantyne (2016) to the same data, and a smaller (moreplausible) relative iron abundance. Since they used relxill (Dauser et al. 2013; Garc´ıa et al. 2014), the two differencesbetween our model and theirs are that we assume a lamp-post geometry as opposed to parameterizing the reflectionemissivity profile with a broken power-law, and that we self-consistently account for a radial disc ionization profile in-stead of assuming a single value of the ionisation parameter.Our fit requires a large value of the ‘boost’ parameter, /B . This sets the normalization of the relativistic reflectionspectrum relative to the direct emission. For /B = , the We note that relativistic corrections are ignored in the Shakura& Sunyaev model. In our model the radial derivative is the rele-vant quantity because the assumed density profile only influencesthe radial ionization gradient , not the absolute value, which is setby the log ξ max model parameter. However, the relativistic cor-rections can change the density gradient and we plan to includethese corrections in a future version of the model. It is exact if the illuminating flux scales as r − . We use thisradius instead of numerically calculating where the ionization pa-rameter peaks to save computational expense. k e V c m s k e V ( d a t a - m o d e l ) / e rr o r Figure 1.
Upper panel Time-averaged energy spectrum unfoldedwith an absorbed power-law model. Lower panel: residuals of thismodel. Clear residuals around the iron line energy range can beseen. relative normalization of the reflection component is set en-tirely by the lamppost geometry and general relativistic lightbending, whereas /B > returns a stronger reflection com-ponent than expected in the lamppost geometry and /B < corresponds to a weaker than expected reflection (see Ingramet al. 2019 for more details). A deviation from the expectedisotropic lamppost emission can be due to the true sourcegeometry being something other than point-like and/or anintrinsic velocity of the plasma in the corona that beams theemission. Since the physical assumptions of reltrans areonly valid for a boost parameter reasonably close to unity,we impose a hard upper limit in our fits of /B = to avoidunrealistic values. We see that the best fitting value in ourfit is close to this upper limit. Interestingly alongside withthe large boost value the fit requires a low iron abundance.Therefore even though the relativistic reflection componentis fundamental to fit to the data, the iron line does not seemprominent enough to require solar iron abundance. A similarbehaviour has been found by Shreeram & Ingram (2020) fora black hole binary. We now model the lag-energy jointly with the time-averagedspectrum, again using reltrans . We again ignore energies < keV in the time-averaged spectrum because of the com-plicated structure of the warm absorber. For the lags, weinstead consider the extended energy range . − . keV.Although soft lags can be caused by the response of thewarm absorber to changes in the ionizing flux (Silva et al.2016) which is not accounted for in our model, these lagsshould only be relevant for variability on longer timescales The configuration of the radial zones in the model is the sameas that used for the time-averaged spectral fitting.MNRAS , 1–14 (2019)
G. Mastroserio et al. k e V c m s k e V ( d a t a - m o d e l ) / e rr o r Figure 2.
Upper panel: Time-averaged energy spectrum un-folded with the best fitting model accounting for direct emissionand relativistic component ( reltrans ) and the distant reflector( xillver ). The narrow emission line at . keV in the sourcerestframe is modelled with a Gaussian component. Lower panel:Residuals of the best fit model. than the Fourier frequency ranges considered in our analysis.We therefore choose to model the reverberation lags downto . keV in order to maximize the signal to noise of thedata.The first model, hereafter Model [1], considers the time-averaged spectrum and lags calculated in the two Fourierfrequency ranges where the source shows high coherence.Since the higher frequency range ( . − . mHz) seems notto be affected by the hard lags (De Marco et al. 2013b; Karaet al. 2013b) we use the pivoting prescription only in thelower frequency range ( . − . Hz). Thus, in the higherfrequency range only the normalisation of the illuminatingpower law is varying. We note that when pivoting is includedin this frequency range, the hard lags also fit to the reverber-ation lags, biasing our results. This happens because eventhough the iron line signal in the lag energy spectrum is sig-nificantly higher than the underlying continuum signal, it isnot sufficiently strong to dominate the fit.The best fitting parameters of Model [1] are listed inTable 1, and the best fitting time-averaged energy spectrumand lag spectra are shown in panels A, B, C in Fig 3 along-side the data. The dashed lines in panel A are the differentcomponents of the time-averaged spectrum. The fit has anacceptable χ / d . o . f . of / . The model has been plot-ted with higher resolution than the data for clarity, whichallows the small reverberation feature to be seen in the lowfrequency range of the model (see also Mastroserio et al.2018), which is unfortunately impossible to detect with theavailable data resolution. Table 1.
Best fitting parameters for the different models in thepaper. First column lists the free parameters in the model; a dashin the model column indicates that the parameter is not part ofthat model. All models include Galactic absorption with hydrogencolumn density fixed to n H = . × cm − , a high energy cut-off fixed to keV, and a narrow Gaussian line in the time-averaged spectrum fixed at . keV in the galaxy rest frame, withcosmological redshift z = . . The black hole is consideredmaximally spinning ( a = . ) in all relativistic models. Models[A] and [B] fit to the time-averaged spectrum only. Model [A] doesnot include the relativistic reflection component, Model [B] does.Model [1] fits two lag spectra in the frequency range . − . mHz.Models [2] fits the complex cross-spectra in the same frequencyranges simultaneously to the time-averaged spectrum, whereasthe last column model fits just the complex cross-spectrum. Errorsare all confidence.Time-ave spec Time-ave &Lag spec Time-ave &Cross specParameter pw + xil [ A ] rel + xil [ B ] model [ ] model [ ] h [ R g ] - + − . + . − . . + . − . r in [ R g ] - + − . + . a . + . − . Γ pl . + . − . . + . − . . + . − . . + . − . log ξ ( xil ) . + . − . . + . − . . + . − . A Fe . + . b . + . − . . + . b . + . b /B (Boost) - . c − . . c − . . c − . Mass[ M (cid:12) ] - - . + . − . . + . d χ / d . o . f . /
127 116 /
124 128 /
144 144 / a The lower limit of the inner radius is the ISCO ( . R g ). b The lower limit of the iron abundance is . . c The upper limit of the boost parameter is d The lower limit of the mass is M (cid:12) Finally we attempt to account for the variability ampli-tude of the reflection component in addition to the timelags by simultaneously fitting the complex cross-spectrumas a function of energy jointly with the time-averaged spec-trum. We now only fit in the − keV energy range, sincethe warm absorber will affect the energy dependence of thecross-spectrum even if it does not contribute to the timelags.We fit simultaneously to the time-averaged spectrumand the complex cross-spectrum in the same two frequencyranges. This model, hereafter Model [2], is therefore thesame as Model [1], except the correlated variability ampli-tude is considered in addition to the time lags and the time-averaged spectrum. The best fitting model is shown in theleft panels A, B and C of Fig. 4, where the upper panel showsthe unfolded time-averaged spectrum and the bottom panelsthe unfolded real (solid lines) and imaginary (dashed lines)parts of the cross-spectrum. Again the hard lag producedby the pivoting power-law is enabled only for the lower fre-quency range. We note that the iron line signal, representedby a single data point, is too weak in the model in both the MNRAS000
144 144 / a The lower limit of the inner radius is the ISCO ( . R g ). b The lower limit of the iron abundance is . . c The upper limit of the boost parameter is d The lower limit of the mass is M (cid:12) Finally we attempt to account for the variability ampli-tude of the reflection component in addition to the timelags by simultaneously fitting the complex cross-spectrumas a function of energy jointly with the time-averaged spec-trum. We now only fit in the − keV energy range, sincethe warm absorber will affect the energy dependence of thecross-spectrum even if it does not contribute to the timelags.We fit simultaneously to the time-averaged spectrumand the complex cross-spectrum in the same two frequencyranges. This model, hereafter Model [2], is therefore thesame as Model [1], except the correlated variability ampli-tude is considered in addition to the time lags and the time-averaged spectrum. The best fitting model is shown in theleft panels A, B and C of Fig. 4, where the upper panel showsthe unfolded time-averaged spectrum and the bottom panelsthe unfolded real (solid lines) and imaginary (dashed lines)parts of the cross-spectrum. Again the hard lag producedby the pivoting power-law is enabled only for the lower fre-quency range. We note that the iron line signal, representedby a single data point, is too weak in the model in both the MNRAS000 , 1–14 (2019) everberation lag in Mrk 335 ( d a t a m o d e l ) e rr o r T i m e l a g ( s e c ) [B] [0.02 0.2] mHz Energy (keV)2.50.02.5 ( d a t a m o d e l ) e rr o r T i m e l a g ( s e c ) [C] [0.2 0.7] mHz ( d a t a m o d e l ) e rr o r ( k e V c m s k e V ) [A] Figure 3.
Joint fit of the time-averaged energy spectrum and thelag energy spectra in different Fourier frequency ranges (specifiedin the panels). Panel A: unfolded time-averaged spectrum, thedashed curves are the model components; panels B and C: timelag spectra. Below every panel we show the residuals of the fit.The time-averaged spectrum is fitted in the − keV energyrange, the time lag spectra in . − keV. The model of thefirst frequency range accounts for both illuminating continuumand reverberation lags whereas only the reverberation lags areaccounted in the second frequency range. Model curves are plottedwith higher resolution than the data for clarity. complex cross-spectra. Even though the χ / d . o . f . of / is still formally acceptable, the structured residuals in theiron line range indicate our model is not able to properlyfit the time-averaged spectrum and complex cross-spectrasimultaneously. It is worth noting that the parameters ofthe fit are similar to those of the previous fit (Model [1]),apart from the black hole mass which is poorly constrainedin this case with only an upper limit of ∼ . × M (cid:12) . Thisis because the model does not reproduce the reverberationfeature in the lag energy spectrum that provides the massconstraint. As an experiment to determine what is driving our Model [2]fit, we fit the cross-spectrum in both frequency ranges with-out considering the time-averaged spectrum. This is instruc-tive because the time-averaged spectrum tends to dominateover the cross-spectrum in the fit due to its far better statis-tics. We refer to this experiment as Model [3]. Right panelsB ∗ and C ∗ of Fig. 4 show the resulting best fitting model,which has a reduced χ of . / and exhibits a prominentiron line feature in both frequency ranges. In panel AA ofFig. 4, instead, we plot the time-averaged spectrum cor-responding to best fitting parameters of Model [3] (lightgreen dashed line). This line has not been fitted to thedata but only superimposed to the time-averaged spectrum.In the same panel as a comparison we plot the Model [2]time-averaged spectrum (dark green solid line). We see thatthe two spectra are dramatically different, with the cross-spectrum-only fit predicting a much stronger and narroweriron line (resulting from a much larger source height ∼ R g and iron abundance ∼ ). We therefore see that the modelcan either reproduce the observed cross-spectrum or the ob-served time-averaged spectrum. The far better statistics ofthe time-averaged spectrum drive the fit to favour it, andtherefore miss the iron line features in the cross-spectrum.This discrepancy is very similar to what Zoghbi et al. (2020)found in NGC 5560, even though in that case the problemalready appeared fitting the lag spectrum. They found thatusing a static lamppost geometry model over-predicts theiron line feature in the time-averaged spectrum. Althoughour model is able to jointly fit lag and time-averaged spec-tra, when we consider the correlated variability amplitudeswe encounter the same problem (compare their Fig.7 withour panel AA of Fig. 4). In order to fit the cross-spectrumour model maximises the iron abundance in the accretiondisc which leads to an unrealistic iron line feature in thetime-averaged spectrum. In our case it seems that the lagsare compatible with the parameters of the time-averagedspectrum, however, the data requires more variability in theiron line compared to what is predicted by the model. Zoghbiet al. (2020) argue the assumed static lamppost geometrycould be the cause of this discrepancy. We expand on thisargument in Section 5.1 also considering our results on theblack hole mass. Inner radius and ionisation profile
When we consider only the time-averaged spectrum(Model [B]) the inner radius is ± R g . This value is largerthan that found by Keek & Ballantyne (2016), who fixed r in to the ISCO and inferred a black hole spin of a = . ± . from the time-averaged spectrum of the same observationthat we consider here. They used relxill for their fit, asnoted in Section 4.1. The configuration of reltrans we usehere is different in two ways: it uses an emissivity profilecalculated in the lamppost geometry instead of parameter-izing it with a broken power-law function, and it considers aradial ionization profile instead of a constant ionization pa-rameter. Shreeram & Ingram (2020) showed that assuminga radial ionisation profile in the accretion disc improves thespectral fit compared to a constant ionisation and the value MNRAS , 1–14 (2019)
G. Mastroserio et al. ( d a t a m o d e l ) e rr o r R e & I m [B] [0.02 0.2] mHz RealImaginary 10 Energy (keV) ( d a t a m o d e l ) e rr o r R e & I m [C] [0.2 0.7] mHz RealImaginary202 ( d a t a m o d e l ) e rr o r ( k e V c m s k e V ) [A] 202 ( d a t a m o d e l ) e rr o r * ]RealImaginary 10 Energy (keV) ( d a t a m o d e l ) e rr o r * ]RealImaginary 10 ( k e V c m s k e V ) [AA]Model [2]Model [3] Figure 4.
Left panels: unfolded spectra of the time-averaged energy spectrum and the real and imaginary part of the cross-spectrumwith the best fitting Model [2]. Panel A shows the time-averaged spectrum (solid line) and the different components (dotted lines). PanelsB to C show the complex cross-spectrum in 2 ranges of frequencies (specified in the panels). Below every panel there are the residualsof the best fitting model. Right panels: Panels B ∗ to C ∗ show the unfolded real and imaginary part of the cross-spectrum with the bestfitting Model [3] that does not consider the time-averages spectrum (the frequency ranges are the same of the left-side panels). In allthe cross-spectrum panels solid lines model the real part and dashed lines model the imaginary part. The top right panel AA shows theunfolded time-averaged spectrum with Model [2] (solid dark green line) superimposed to the time-averaged spectrum calculated from thebest fitting parameters of Model [3] (dashed light green line). For all the panels the pivoting prescription has been accounted for only inthe lowest frequency range. MNRAS000
Left panels: unfolded spectra of the time-averaged energy spectrum and the real and imaginary part of the cross-spectrumwith the best fitting Model [2]. Panel A shows the time-averaged spectrum (solid line) and the different components (dotted lines). PanelsB to C show the complex cross-spectrum in 2 ranges of frequencies (specified in the panels). Below every panel there are the residualsof the best fitting model. Right panels: Panels B ∗ to C ∗ show the unfolded real and imaginary part of the cross-spectrum with the bestfitting Model [3] that does not consider the time-averages spectrum (the frequency ranges are the same of the left-side panels). In allthe cross-spectrum panels solid lines model the real part and dashed lines model the imaginary part. The top right panel AA shows theunfolded time-averaged spectrum with Model [2] (solid dark green line) superimposed to the time-averaged spectrum calculated from thebest fitting parameters of Model [3] (dashed light green line). For all the panels the pivoting prescription has been accounted for only inthe lowest frequency range. MNRAS000 , 1–14 (2019) everberation lag in Mrk 335 of the inner disc radius depends on the nature of this profile.The radial ionisation profile quite dramatically changes theshape of the iron line, in particular the red wing (see Fig. 9in Ingram et al. 2019), requiring less broadening from therelativistic effects. In order to calculate the radial ionizationprofile, our model assumes the density profile correspond-ing to zone A of the Shakura & Sunyaev (1973) disc model,in which the dominant sources of pressure and opacity are,respectively, radiation and electron scattering. According toShakura & Sunyaev (1973), the outer limit of zone A in anAGN with a . − keV luminosity of × erg/s (Grupeet al. 2007), and assuming viscosity parameter α = . andblack hole mass or × M (cid:12) , is R ab ∼ R g or ∼ R g (increasing α increases R ab ). Thus we can comfort-ably consider most of the emission to come from a zone Adensity profile accretion disc. We note that the disc scaleheight in the zone A regime can be significantly greater thanzero (see e.g. Taylor & Reynolds 2018a), whereas here weassume an infinitely thin disc. Since the disc scale heightcan impact on the reverberation signal (Taylor & Reynolds2018b), we plan to include a more general disc geometryin our model in future. Including the lag spectra in the fit(Model [1]) does not change the best fitting inner radius.However, the lower limits are now compatible with ISCOwithin 90% confidence, because the broad and strong ironline observed in the lag spectrum favours a smaller innerradius. Iron abundance
It is also notable that the iron abundance yielded by ourmodels, A Fe ∼ . − . , differs significantly from the largevalue of A Fe = . ± . returned by the fits of Keek & Ballan-tyne (2016). It is common for X-ray reflection modelling tomeasure significantly super-solar iron abundances both forAGN and black hole X-ray binaries (e.g. Dauser et al. 2012;Parker et al. 2015; Garc´ıa et al. 2015; see Garc´ıa et al. 2018for a brief review). Again it is possible that including a ra-dial ionization profile goes some way to reducing the inferrediron abundance, as speculated by Ingram et al. (2019). MCMC simulation
When we include the variability amplitude by fitting thecomplex cross-spectrum instead of the lag spectrum (Model[2]), although we achieve a statistically acceptable χ , wesee clear residual structure in the iron line region and there-fore do not consider this model to adequately describe thedata. We therefore focus our further analysis on the resultsof fitting the time-averaged spectrum together with the lagspectrum (Model [1]). In order to probe the parameter spaceof this model and look for parameter correlations we run aMarkov Chain Monte Carlo (MCMC) simulation in xspec with walkers and a total length of million steps. Thewalkers start from the Model [1] best fit and we set a burn-inperiod of steps to assure the convergence of the chain.Fig. 5 shows the integrated probability distribution for someof the model parameters, the purple regions represent the σ contour, orange and yellow and σ . The inner radiusis strongly positively correlated with the boost parameter.The reason for this is that both contribute to the strength of the reflection emission in the opposite sense: when theinner radius is closer to ISCO the reflection is stronger anda lower value of the boost is required. However, the best fitvalue of the boost parameter is significantly above whichmeans that the data require a stronger reflection componentthan the one produced by the lamppost model. This can beinterpreted as the velocity of the plasma in the lamppost be-ing non-zero and pointing towards the black hole, boostingthe flux illuminating the disc. Contrary to what has beenseen in previous works the black hole mass is not degener-ate with the height of the lamppost (e.g. Cackett et al. 2014)nor with any other parameter (Mastroserio et al. 2019 founda correlation with the ionisation parameter). This is likelydue to the wide range of timescales and the simultaneousfit to the time-averaged spectrum that we considered in thecurrent work. In Section 5.1 we further study the relationbetween mass and height of the lamppost. In this work we have jointly modelled the time-averaged X-ray spectrum and the short timescale X-ray variability of theSeyfert galaxy Mrk 335 using our relativistic reflection model reltrans . We considered the high-flux XMM-Newton ob-servation from 2006, in which Kara et al. (2013b) found ev-idence of reverberation lags in the ∼ . − . mHz Fourierfrequency range. We first modelled the time-averaged energyspectrum alone (Model [B]) before additionally modellingthe lag-energy spectrum simultaneously in two frequencyranges ( . − . mHz and . − . mHz; Model [1]), with ourmodel for the lowest frequency range including a prescrip-tion to account for the intrinsic hard lags (Mastroserio et al.2018, 2019). We then extended our analysis to fit simultane-ously to the complex cross-spectrum and the time-averagedspectrum (Model [2]). We achieve statistically acceptable fitsfor both Model [1] and Model [2]. However, the Model [2] fitsuffers from structured residuals around the iron line region,and so here we further discuss the details of our more suc-cessful Model [1] fit as well as the possible reasons for thedeficiencies in our Model [2] fit. The fits that consider timing properties are sensitive to blackhole mass. In previous analyses, fits to the lag-energy spec-trum have suffered from a degeneracy between the sourceheight and the black hole mass (Cackett et al. 2014; Ingramet al. 2019). This occurs because the same reverberation lagcan be reproduced if the source height is many small gravi-tational radii or a few large gravitational radii. We see fromthe bottom left panel in Fig. 5, and from the more detailedcontour plot of χ as a function of black hole mass andsource height presented in Fig. 6, that there is no such de-generacy in our Model [1] fit. This is because we considertwo frequency ranges instead of only one. This point can beillustrated by comparison with the reltrans fit to only the . − . mHz lag-energy spectrum of the same observationpresented in Fig. 17 of Ingram et al. (2019). In that case,the best fitting black hole mass was ∼ × M (cid:12) , but othersolutions with mass values as high as ∼ × M (cid:12) were sta-tistically acceptable within σ confidence. These high mass MNRAS , 1–14 (2019) G. Mastroserio et al.
Figure 5.
Probability distribution from the MCMC simulation of Model [1] (see text for the details). The purple, orange and yellow are , , σ contours respectively. The inner radius is expressed in units of ISCO which is . R g in this case. solutions were in the phase-wrapping regime, and are ruledout from our new fit because they predict very large rever-beration lags in the . − . mHz frequency range that arenot observed.Our best-fitting black hole mass of . + . − . × M (cid:12) (errors are 90% confidence) is roughly one order of magni-tude smaller than the optical reverberation mass measure-ments ( M = . ± . × M (cid:12) ; Peterson et al. (2004) and ± × M (cid:12) ; Grier et al. 2012). These higher mass valuesare incompatible with our modelling with high statisticalconfidence since, as discussed above, parameter combina- tions that can reproduce the . − . mHz lag-energy spec-trum with such high masses predict large lags in . − . mHz that are not observed. This discrepancy is not due todilution of the reverberation lag by the presence of contin-uum emission variability with no time lags between energybands (see e.g Uttley et al. 2014) because dilution is fully ac-counted for in our model – although there could be an extrasource of dilution that is not in our model. The discrepancycould be because we assume a lamppost corona, if in re-ality the corona is extended (as was suggested by Wilkins& Gallo 2015 for this very observation). The schematic in MNRAS000
Probability distribution from the MCMC simulation of Model [1] (see text for the details). The purple, orange and yellow are , , σ contours respectively. The inner radius is expressed in units of ISCO which is . R g in this case. solutions were in the phase-wrapping regime, and are ruledout from our new fit because they predict very large rever-beration lags in the . − . mHz frequency range that arenot observed.Our best-fitting black hole mass of . + . − . × M (cid:12) (errors are 90% confidence) is roughly one order of magni-tude smaller than the optical reverberation mass measure-ments ( M = . ± . × M (cid:12) ; Peterson et al. (2004) and ± × M (cid:12) ; Grier et al. 2012). These higher mass valuesare incompatible with our modelling with high statisticalconfidence since, as discussed above, parameter combina- tions that can reproduce the . − . mHz lag-energy spec-trum with such high masses predict large lags in . − . mHz that are not observed. This discrepancy is not due todilution of the reverberation lag by the presence of contin-uum emission variability with no time lags between energybands (see e.g Uttley et al. 2014) because dilution is fully ac-counted for in our model – although there could be an extrasource of dilution that is not in our model. The discrepancycould be because we assume a lamppost corona, if in re-ality the corona is extended (as was suggested by Wilkins& Gallo 2015 for this very observation). The schematic in MNRAS000 , 1–14 (2019) everberation lag in Mrk 335 h [ R g ]10 M a ss ( M ) Figure 6. χ contour plot of black hole mass and sourceheight resulting from Model [1]. Fig. 7 illustrates how an extended corona model could inprinciple accommodate a higher mass than the lamppostmodel. The top panel illustrates our best fit, whereby thereverberation lag roughly results from the ∼ R g distancefrom the corona to the disc. The second panel shows thatthe same physical distance from corona to disc cannot bereproduced if the mass is increased by a factor of 14. This isbecause the physical size of the ∼ . R g disc inner radiusis now a factor of 14 larger (and the inner radius in units of R g is constrained by the iron line profile). The third panelillustrates that a distance of ∼ × R g is possible for an M ∼ × M (cid:12) black hole with an extended corona at aheight of ∼ . R g .Previous X-ray reverberation analyses of Mrk 335 havealso inferred a higher black hole mass than what we findhere. Chainakun et al. (2016) fit to the lag-energy spec-trum of the same observation that we consider here. Theirmodel assumes a lamppost geometry and is largely simi-lar to ours, with some small differences such as the form ofthe radial ionization profile; their best fitting mass valuewas M ∼ . × M (cid:12) . The major differences betweentheir analysis and ours are that they only considered onefrequency range, whereas we consider two (their one fre-quency range corresponds to our higher frequency range),and that we conduct a more thorough exploration of param-eter space. Given that the single frequency fit presented inIngram et al. (2019) included statistically acceptable phase-wrapped solutions with high black hole mass, it could be thatthe Chainakun et al. (2016) fit is also in the phase-wrappingregime. If this is the case, then the high black hole massreturned by their fit would also be ruled out upon consider-ation of a second, lower frequency range. Emmanoulopouloset al. (2014) used a different procedure of fitting the lag-frequency spectrum between the . − keV and . − keVenergy bands and obtained a mass of M ∼ × M (cid:12) . How-ever, they did not simultaneously fit to the time-averaged ∼ 11.4 × 10 R g ⊙ ∼ 17.6 × 10 R g ⊙ ∼ 21 × 10 R g ⊙ M ∼ 10 M ⊙ ∼ 17.6 × 14 × 10 R g ⊙ > 246 × 10 R g ⊙ M ∼ 14 × 10 M ⊙ ∼ 1.5 × 14 × 10 R g ⊙ ∼ 17.6 × 14 × 10 R g ⊙ ∼ 21 × 10 R g ⊙ M ∼ 14 × 10 M ⊙ Figure 7.
Schematic comparing the lamppost geometry with aradially extended corona geometry ( R g (cid:12) = GM (cid:12) / c ). The toppanel represents our Model [1] fit. For the middle panel, all pa-rameters are the same, except the black hole mass is now 14 timeslarger. The physical distance from the corona to the disc now mustbe larger than it was for the top panel, no matter how low thesource is placed. The bottom panel illustrates a radially extendedcorona. We see that the same physical distance from the coronato the disc can now be achieved even for a × M (cid:12) black hole. spectrum, and they simply assumed that respectively and of the flux in the . − keV and . − keV bandsis due to reflection, whereas we include a full model of therest frame reflection spectrum.It is also worth investigating whether the assumptionsgoing into the optical reverberation mass measurements canbe reasonably changed in order to obtain values compati-ble with our X-ray reverberation measurement. Optical re-verberation mapping consists of measuring the time delay τ between variations in the optical continuum, assumed tocome from the accretion disc, and corresponding variationsin the broad optical emission lines, assumed to come fromclouds orbiting the black hole at a distance far larger thanthe accretion disc in the so-called broad line region (BLR).Assuming that τ is a light crossing delay, and assuming thatthe measured (full width at half maximum) velocity widthof the lines, V FWHM , results from virialized rotation aroundthe black hole, the black hole mass is given by M = f c τ ( V FWHM / ) G . (7)Here, f is a factor that depends on the structure, kinemat-ics and orientation of the BLR ( f = for a BLR consisting MNRAS , 1–14 (2019) G. Mastroserio et al. of a thin spherical shell with an isotropic velocity distribu-tion; Onken et al. 2004). Peterson et al. (2004) set f = . for their mass measurement of M = . ± . × M (cid:12) forMrk 335, which is the value that results from the reason-able assumption that AGN and quiescent galaxies follow onaverage the same relationship between black hole mass andbulge stellar velocity dispersion (Onken et al. 2004). It is inprinciple possible to adjust the assumptions of optical rever-beration mapping to obtain M ∼ M (cid:12) for Mrk 335, butthe required adjustment, setting f = . , therefore seemsunlikely. Moreover, such a low mass is in tension with themeasured ˜AˇE flux, from which Peterson et al. (2004)infer a bolometric luminosity of ∼ . × erg s − . This lu-minosity is equal to the Eddington limit for a black hole ofmass M ∼ × M (cid:12) , which roughly coincides with our confidence upper limit. Thus our mass value implies Edding-ton limited or mildly super-Eddington accretion in Mrk 335,which is not impossible but has not previously been thoughtto be the case for Mrk 335. Although we have achieved a good fit for the time lags, weare skeptical of our fit for the full cross-spectrum because itsuffers from structured residuals around the iron line region(although formally the fit is still statistically acceptable). Ittherefore seems that the model has problems reproducingthe energy-dependent correlated variability amplitude. Thiscould again potentially be because we assume a lamppost ge-ometry, if in reality the corona is extended. Wilkins & Gallo(2015) inferred a radial extent of ∼ R g during the high-flux epoch that we consider, followed by a contraction whenthe flux decreased. In such an extended corona, intrinsichard lags may arise due to inward propagation of accretionrate fluctuations (Wilkins et al. 2016), which we instead ac-count for with a fluctuating power-law index. Our prescrip-tion provides a good mathematical description of fluctua-tions propagating through a disc and into a compact corona(Uttley & Malzac in prep.), but likely differs significantlyfrom propagation in an extended corona. The reverberationsignal itself would also differ significantly from the lamppostcase (see e.g. Fig. 7 but also Wilkins et al. 2016).We note that our analysis is the first ever attemptto jointly fit the full energy-dependent cross-spectrum andtime-averaged spectrum of an AGN. However, we have pre-viously applied the same analysis to the black hole X-raybinary Cygnus X-1 (Mastroserio et al. 2019), and in thatcase the model worked very well. We cannot currently ruleout this being due to the two systems having different coro-nal geometries, with Cygnus X-1 having a more compactcorona than Mrk 335 during the observations we analysed.Alternatively, the two sources could have roughly similar ex-tended coronal geometries, with our lamppost model betterable to approximate the true signal for the lower frequen-cies (in terms of R g / c ) typically probed by X-ray binariesthan the higher frequencies (in terms of R g / c ) probed byAGN. High signal to noise NICER observations of X-ray bi-naries, for which timing properties can be constrained up tofrequencies as high as ∼ Hz, may enable this hypothesisto be tested. The signal for AGN may also be more signif-icantly modified by absorption and winds, given the verydifferent environment associated with the two classes of ob- ject and the different opacity in the disc (e.g. Miller et al.2010; Mizumoto et al. 2018)
The results of our X-ray spectral-timing analysis of Mrk 335using the model reltrans can be summarised as follows:(i) Our fits return a best fitting disc inner radius that istruncated outside of the ISCO for any value of black holespin ( R in ≈ R g ), and a sub-solar iron abundance ( A Fe < . ). Although our fits that include timing properties are con-sistent with R in at the ISCO within confidence, ourresults differ from previous spectral fits for the same obser-vation (Keek & Ballantyne 2016), which returned a smalldisc inner radius and a super-solar iron abundance. We at-tribute this to the more realistic treatment of the ionisationstructure of the accretion disc in our model, allowing fora radial profile rather than assuming a constant ionisationparameter. This suggests that models with a radial ioniza-tion profile may return larger disc inner radii and smalleriron abundances for other AGN than the constant ionisa-tion models that have mainly been used in the past.(ii) Our simultaneous fit to the time-averaged spectrumand the lag-energy spectrum (Model [1]) returns a black holemass of M = . + . − . × M (cid:12) , which is not compatible withUV/optical reverberation mapping results ( − × M (cid:12) ).We argue that our use of the lamppost geometry might beresponsible for this low mass value, whereas an extendedcorona model may return a larger mass. Proper modellingof an extended corona is necessary to confirm this claim.(iii) Our model is not satisfactorily able to simultaneouslydescribe the energy dependent correlated variability ampli-tude of Mrk 335 and the time-averaged spectrum. Fittingto the cross-spectrum alone returns parameters that predicta very strong iron line in the time-averaged spectrum thatis not observed. Again this seems to challenge the lamppostgeometry, as it has been already suggested by Wilkins &Gallo (2015). An extended corona could provide more cor-related variability amplitude in the reflection response dueto the different illumination of the accretion disc, withoutsignificantly changing the lags because the distance betweenthe corona and the disc is similar to the lamppost case. ACKNOWLEDGEMENTS
The authors would like to acknowledge the anonymous ref-eree for the very helpful comments and suggestions. G.M. ac-knowledges the support from NASA grant 80NSSC19K1020and together with M.K. the support from NWO Spinoza.A.I. acknowledges the Royal Society.
DATA AVAILABILITY
The data underlying this article are available from theXMM-Newton science archive ( http://nxsa.esac.esa.int/ ). MNRAS000
The data underlying this article are available from theXMM-Newton science archive ( http://nxsa.esac.esa.int/ ). MNRAS000 , 1–14 (2019) everberation lag in Mrk 335 REFERENCES
Ar´evalo P., Uttley P., 2006, MNRAS, 367, 801Cackett E. M., Zoghbi A., Reynolds C., Fabian A. C., Kara E.,Uttley P., Wilkins D. R., 2014, MNRAS, 438, 2980Campana S., Stella L., 1995, MNRAS, 272, 585Chainakun P., Young A. J., Kara E., 2016, MNRAS, 460, 3076Dauser T., et al., 2012, MNRAS, 422, 1914Dauser T., Garcia J., Wilms J., B¨ock M., Brenneman L. W.,Falanga M., Fukumura K., Reynolds C. S., 2013, MNRAS,430, 1694De Marco B., Ponti G., Cappi M., Dadina M., Uttley P., CackettE. M., Fabian A. C., Miniutti G., 2013a, MNRAS, 431, 2441De Marco B., Ponti G., Miniutti G., Belloni T., Cappi M., DadinaM., Mu˜noz-Darias T., 2013b, MNRAS, 436, 3782Eardley D. M., Lightman A. P., Shapiro S. L., 1975, ApJ, 199,L153Emmanoulopoulos D., Papadakis I. E., Dovˇciak M., McHardyI. M., 2014, MNRAS, 439, 3931Epitropakis A., Papadakis I. E., 2016, A&A, 591, A113Fabian A. C., Rees M. J., Stella L., White N. E., 1989, MNRAS,238, 729Fabian A. C., Iwasawa K., Reynolds C. S., Young A. J., 2000,PASP, 112, 1145Fabian A. C., et al., 2009, Nature, 459, 540Garc´ıa J., Dauser T., Reynolds C. S., Kallman T. R., McClintockJ. E., Wilms J., Eikmann W., 2013, ApJ, 768, 146Garc´ıa J., et al., 2014, ApJ, 782, 76Garc´ıa J. A., Steiner J. F., McClintock J. E., Remillard R. A.,Grinberg V., Dauser T., 2015, ApJ, 813, 84Garc´ıa J. A., Fabian A. C., Kallman T. R., Dauser T., ParkerM. L., McClintock J. E., Steiner J. F., Wilms J., 2016, MN-RAS, 462, 751Garc´ıa J. A., Kallman T. R., Bautista M., Mendoza C., DeprinceJ., Palmeri P., Quinet P., 2018, in Workshop on AstrophysicalOpacities. p. 282 ( arXiv:1805.00581 )Grier C. J., et al., 2012, ApJ, 744, L4Grupe D., Komossa S., Gallo L. C., 2007, ApJ, 668, L111Haardt F., Maraschi L., 1991, ApJ, 380, L51Haardt F., Maraschi L., 1993, ApJ, 413, 507Huchra J. P., Vogeley M. S., Geller M. J., 1999, ApJS, 121, 287Ingram A., van der Klis M., 2013, MNRAS, 434, 1476Ingram A., van der Klis M., Middleton M., Done C., AltamiranoD., Heil L., Uttley P., Axelsson M., 2016, MNRAS, 461, 1967Ingram A., Mastroserio G., Dauser T., Hovenkamp P., van derKlis M., Garc´ıa J. A., 2019, MNRAS, p. 1677Kalberla P. M. W., Burton W. B., Hartmann D., Arnal E. M.,Bajaja E., Morras R., P¨oppel W. G. L., 2005, A&A, 440, 775Kara E., Fabian A. C., Cackett E. M., Steiner J. F., Uttley P.,Wilkins D. R., Zoghbi A., 2013a, MNRAS, 428, 2795Kara E., Fabian A. C., Cackett E. M., Miniutti G., Uttley P.,2013b, MNRAS, 430, 1408Kara E., Alston W. N., Fabian A. C., Cackett E. M., Uttley P.,Reynolds C. S., Zoghbi A., 2016, MNRAS, 462, 511Keek L., Ballantyne D. R., 2016, MNRAS, 456, 2722Lightman A. P., White T. R., 1988, ApJ, 335, 57Longinotti A. L., et al., 2013, ApJ, 766, 104Lyubarskii Y. E., 1997, MNRAS, 292, 679Malkan M. A., 1983, ApJ, 268, 582Mastroserio G., Ingram A., van der Klis M., 2018, MNRAS, 475,4027Mastroserio G., Ingram A., van der Klis M., 2019, MNRAS,p. 1682Matt G., Perola G. C., Piro L., 1991, A&A, 247, 25McHardy I. M., Papadakis I. E., Uttley P., Page M. J., MasonK. O., 2004, MNRAS, 348, 783Miller L., Turner T. J., Reeves J. N., Braito V., 2010, MNRAS,408, 1928 Mizumoto M., Done C., Hagino K., Ebisawa K., Tsujimoto M.,Odaka H., 2018, MNRAS, 478, 971Nowak M. A., Vaughan B. A., Wilms J., Dove J. B., BegelmanM. C., 1999, ApJ, 510, 874Onken C. A., Ferrarese L., Merritt D., Peterson B. M., PoggeR. W., Vestergaard M., Wandel A., 2004, ApJ, 615, 645Papadakis I. E., Nandra K., Kazanas D., 2001, ApJ, 554, L133Parker M. L., et al., 2015, ApJ, 808, 9Peterson B. M., et al., 2004, ApJ, 613, 682Rapisarda S., Ingram A., Kalamkar M., van der Klis M., 2016,MNRAS, 462, 4078Reynolds C. S., 2014, Space Sci. Rev., 183, 277Reynolds C. S., Young A. J., Begelman M. C., Fabian A. C., 1999,ApJ, 514, 164Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337Shields G. A., 1978, Nature, 272, 706Shreeram S., Ingram A., 2020, MNRAS, 492, 405Silva C. V., Uttley P., Costantini E., 2016, A&A, 596, A79Taylor C., Reynolds C. S., 2018a, ApJ, 855, 120Taylor C., Reynolds C. S., 2018b, ApJ, 868, 109Thorne K. S., Price R. H., 1975, ApJ, 195, L101Uttley P., Wilkinson T., Cassatella P., Wilms J., Pottschmidt K.,Hanke M., B¨ock M., 2011, MNRAS, 414, L60Uttley P., Cackett E. M., Fabian A. C., Kara E., Wilkins D. R.,2014, A&ARv, 22, 72Vaughan B. A., Nowak M. A., 1997, ApJ, 474, L43Wilkins D. R., Fabian A. C., 2013, MNRAS, 430, 247Wilkins D. R., Gallo L. C., 2015, MNRAS, 449, 129Wilkins D. R., Cackett E. M., Fabian A. C., Reynolds C. S., 2016,MNRAS, 458, 200Wilms J., Allen A., McCray R., 2000, ApJ, 542, 914Zoghbi A., Fabian A. C., Uttley P., Miniutti G., Gallo L. C.,Reynolds C. S., Miller J. M., Ponti G., 2010, MNRAS, 401,2419Zoghbi A., et al., 2014, ApJ, 789, 56Zoghbi A., Kalli S., Miller J. M., Mizumoto M., 2020, ApJ, 893,97van der Klis M., Hasinger G., Stella L., Langmeier A., vanParadijs J., Lewin W. H. G., 1987, ApJ, 319, L13
APPENDIX A: TRANSFER FUNCTIONS
The three transfer functions are W ( E , ν ) = ∫ α,β g (cid:15) e i π ( + z ) τν R( E d ) d α d β (A1) W ( E , ν ) = ∫ α,β g (cid:15) e i π ( + z ) τν ln g sd R( E d ) d α d β (A2) W ( E , ν ) = ∫ α,β g (cid:15) e i π ( + z ) τν ∂ R ∂ Γ ( E d ) d α d β. (A3)Here, R( E ) is the restframe specific intensity of reflectedemission emergent from the disc and g doz = g do /( + z ) ,where g do and g sd are respectively the blueshift experiencedby a photon travelling from the disc to infinity and from thepoint source to the disc in an asymptotically flat and sta-tionary spacetime. α and β are the horizontal and verticalimpact parameters at infinity, such that a given patch of thedisc appears on the image plane to be centered on coordi-nates ( α, β ) and subtends a solid angle d α d β / D where D isthe distance between the observer and the black hole. τ isthe time delay between the arrival at the observer plane ofphotons direct from the corona and those that reflected off adisc patch specified by the coordinates ( α, β ) . (cid:15) is the emis-sivity, which is a function of radius only (see Equation 20 of MNRAS , 1–14 (2019) G. Mastroserio et al.
Ingram et al. 2019). Note that these expressions explicitlyinclude cosmological time dilation due to the factor of ( + z ) in the index of the exponential. This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000