Revealing the accretion disc corona in Mrk 335 with multi-epoch X-ray spectroscopy

Active galactic nuclei host an accretion disc with an X-ray producing corona around a supermassive black hole. In bright sources, such as the Seyfert 1 galaxy Mrk 335, reflection of the coronal emission off the accretion disc has been observed. Reflection produces spectral features such as an iron emission line, which allow for properties of the inner accretion disc and the corona to be constrained. We perform a multi-epoch spectral analysis of all XMM-Newton, Suzaku, and NuSTAR observations of Mrk 335, and we optimize our fitting procedure to unveil correlations between the Eddington ratio and the spectral parameters. We find that the disc's ionization parameter correlates strongly with the Eddington ratio: the inner disc is more strongly ionized at higher flux. The slope of the correlation is less steep than previously predicted. Furthermore, the cut-off of the power-law spectrum increases in energy with the Eddington ratio, whereas the reflection fraction exhibits a decrease. We interpret this behaviour as geometrical changes of the corona as a function of the accretion rate. Below ~10% of the Eddington limit, the compact and optically thick corona is located close to the inner disc, whereas at higher accretion rates the corona is likely optically thin and extends vertically further away from the disc surface. Furthermore, we find a soft excess that consists of two components. In addition to a contribution from reflection in low ionization states, a second component is present that traces the overall flux.


INTRODUCTION
One of the most direct probes of the physics of the inner accretion disc around an active galactic nucleus (AGN) is through its X-ray emission. As indicated by its rapid variability, the high-energy radiation is produced in the innermost regions of the accretion flow, typically only several to tens of gravitational radii (rg = GM/c 2 , where M is the mass of the black hole and c is the speed of light) from the event horizon of the black hole (e.g., Grandi et al. 1992;McHardy et al. 2006;Uttley 2007;Zoghbi et al. 2013). This size scale is < ∼ 1(M/10 7 M⊙) AU, which shows that it will be impossible to image the X-ray emitting regions of AGNs for the foreseeable future. X-ray spectroscopy is, therefore, the best tool to test theories of accretion physics close to the central black hole.
Fortunately, many AGNs exhibit features in their Xray spectra that likely originate from the accretion disc (Fabian & Ross 2010). These X-ray 'reflection' features are produced by the accretion disc as a result of its illu-⋆ E-mail: l.keek@gatech.edu mination from an external X-ray source, such as a magnetically dominated corona residing above the surface of the disc (e.g., Galeev et al. 1979;Haardt & Maraschi 1991, 1993Haardt et al. 1994;Di Matteo 1998;Merloni & Fabian 2001). The strongest feature in the X-ray reflection spectrum is typically the Fe Kα emission line at 6.4 keV that is due to fluorescence in a dense and relatively cold medium (e.g., Barr et al. 1985;Nandra et al. 1989;Pounds et al. 1990;George & Fabian 1991;Nandra & Pounds 1994). The shape of the Fe Kα line can be a powerful probe of the curvature of space-time Laor 1991), leading to measurements of the spin of the super-massive black hole (e.g., Brenneman & Reynolds 2006;Reynolds 2013). Indeed, observations over the last decade by Chandra, XMM-Newton, Suzaku and NuSTAR have uncovered several examples of relativistically broadened Fe Kα lines in the spectra of AGNs (e.g., Ballantyne et al. 2003;Reynolds & Nowak 2003;Miller 2007;Nandra et al. 2007;Fabian et al. 2009;Emmanoulopoulos et al. 2011;Nardini et al. 2012;Walton et al. 2012Walton et al. , 2013Risaliti et al. 2013;Parker et al. 2014).
Models of X-ray reflection spectra (e.g. Ross García et al. 2013) show many features in addition to the Fe Kα line that collectively can constrain the ionization state of the illuminated surface (Ross et al. 1999), often parametrized by ξ = 4πF/nH, where F is the Xray flux irradiating the reflector and nH is the hydrogen number density of the reflecting surface. Since the Xray emission in AGNs likely originates from the corona, a measurement of ξ may provide information on the density structure of the disc as well as the illuminating conditions produced by the corona. Coronal generation models suggest that the fraction of the accretion power dissipated in the corona, f , can depend on the accretion rate through the disc (Stella & Rosner 1984;Merloni & Fabian 2002;Blackman & Pessah 2009). Thus, there could be a dependence of ξ on the accretion rate, which, if measured, can test ideas on the geometry of the corona and how energy flows into it. Indeed, Ballantyne et al. (2011) considered αdisc accretion theory and predicted a specific relationship between ξ and the AGN Eddington ratio that depended on various disc-corona properties such as f , the viscosity parameter α (Shakura & Sunyaev 1973), and the black hole spin, a. Therefore, if ξ can be estimated for a range of Eddington ratios from a number of AGNs, significant insight into the workings of AGN coronae could be obtained. Ballantyne et al. (2011) examined measurements of ξ culled from the literature and showed that, although there is significant scatter, in general objects with larger Eddington ratios tend to exhibit ionized inner accretion discs. While suggestive, there was too much potential object-to-object variability in the data to draw strong conclusions from the result presented by Ballantyne et al. (2011). A significant improvement can be made by determining ξ for several observations of a single AGN, as this should reduce the scatter and allow for a clearer interpretation of any relationship with ξ. Fortunately, the XMM-Newton, Suzaku, and NuSTAR archives now contain many AGNs that have multiple observations with which to perform such an experiment. This paper presents the first results of this project by analysing 12 X-ray spectra of the Seyfert 1 galaxy Mrk 335 that span nearly a factor of 10 in flux. This AGN has several properties that make it an ideal first target for this investigation: it has been observed several times by XMM-Newton, Suzaku and NuSTAR with 3 observations exceeding 100 ks of exposure time (Gondoin et al. 2002;Crummy et al. 2006;Larsson et al. 2008;Grupe et al. 2008Grupe et al. , 2012Parker et al. 2014;Gallo et al. 2015), is known to have a relativistically broadened Fe Kα line (e.g., , an Eddington ratio that can approach unity (Vasudevan & Fabian 2007, 2009, and a reverberation-mapped black-hole mass estimate (Peterson et al. 2004;Grier et al. 2012b). As is seen below, Mrk 335 provides an excellent test case for the predicted ξ-Eddington ratio correlation, and a window into the physics of AGN coronae.
After describing the employed observations and the extraction of the X-ray spectra from the data (Section 2), we first analyse the spectra with a simple power-law model (Section 3). These fits indicate that spectral properties such as the shape of the iron line correlate with flux. Detailed spectral models including ionized reflection are used to quantify these correlations, including the ξ-Eddington ratio correlation (Section 4). Finally, we interpret the observed behaviour as changes in the coronal optical depth and geometry as a function of the Eddington ratio (Section 5).

Observations and data products
In order to explore ionized reflection over a wide range of source fluxes, we employ all available X-ray observations of Mrk 335 performed with the XMM-Newton (Jansen et al. 2001), Suzaku (Mitsuda et al. 2007), and NuSTAR (Harrison et al. 2013) observatories (Table 1). All three observatories host a focusing telescope with grazingincidence mirrors and CCD imaging detectors that are sensitive around 6.4 keV, where the Fe Kα line appears as a prominent reflection feature. With the exception of the most recent NuSTAR pointing (XII; see Table 1), all observations have been analysed before and shown to exhibit the iron line. The combined exposure time is 0.941 Ms : the most used in any detailed X-ray spectral study of Mrk 335. NuSTAR exposures IX and X are fully covered by Suzaku observation VII.
Creation of the data products for our spectral analysis is similar for all observatories. First we apply the latest calibration to the observation data. Next, an image and a light curve are created for the full detector in the full instrumental band-pass. We select source events from within a circle around the source position, and background events from an off-source region with a similar area that is devoid of other X-ray sources. By comparing the source to the detector light curve, we select times without strong background emission. This is only an issue for the XMM-Newton observations, and reduces their effective exposure time. Excluding the time intervals with background flares, we extract the source and background spectra, and we generate the accompanying response files. Neighbouring bins of the source spectra with fewer than 25 counts are grouped, such that χ 2 statistics are applicable. For most observations the source flux and the instrument's data mode was such that pile-up was no issue, observation I being the exception. Below we describe details of this procedure for the different observatories.

XMM-Newton
We use data from the pn camera (Strüder et al. 2001), which is part of the European Photon Imaging Camera (EPIC). We do not include data from the MOS cameras, because of their smaller effective area, nor from the RGS, which is not sensitive in the energy band that we will consider (> 3 keV). Data products are created from the observation data files using the Science Analysis System (SAS) v14.0.0. At high count rates, multiple photons may be miscounted as one event (e.g., Ballet 1999). We check for so-called pile-up by comparing the distribution of event types as a function of energy to the distribution predicted by SAS task epatplot. Substantial pile-up is only present in observation I, when the source exhibited a high flux and the detector was in "Full Window" mode with a long CCD read-out time. For that observation we exclude the piled-up part of the point-spread function, which is the inner 15 ′′ around the source position. The outer radius of the source extraction region is 37.5 ′′ , and we extract spectra in the 0.3 − 12 keV energy range. The pn spectra typically have the highest quality among our sample, with the ratio of the source and background count rates being 135 to 453 (depending mostly on the source flux).

Suzaku
Data products for the different Suzaku instruments are created using the software package HEASOFT v6.16. For the X-Ray Imaging Spectrometer (XIS) we extract spectra in the energy range 0.3 − 10 keV from a 261 ′′ region around the source location. We combine the spectra of all operational front-illuminated CCDs (three for observation VI and two for VII and VIII), and we also use the back-illuminated CCD. The calibration of the front-and back-illuminated CCDs produces fluxes that differ by at most a few per cent, 1 which is unnoticeable given the data quality. The ratio of source and background count rates for the XIS spectra ranges from 3 to 33.
The Hard X-ray detector (HXD) is a non-imaging instrument with a collimated 4.5 • × 4.5 • field of view. It has two sensors with a combined band-pass of 10 − 600 keV. Above 40 keV, however, the the background dominates completely, and we only consider data from the PIN silicon diodes in the 15 − 40 keV range. The background is strongly dominated (∼ 95%) by the particle background at the satellite. 1 We use the simulated non-X-ray background provided by the instrument team for the respective observations, and we employ the model of Boldt (1987) 1 to account for the small contribution from the cosmic X-ray background. From cross-calibration of Suzaku's instruments, 1 we include an enhancement factor for the PIN spectra with respect to XIS, whose value depends on the pointing of the instruments: 1.18 for observation VI and 1.16 for VII and VIII.

NuSTAR
We employ both Focal-Plane Modules: FPMA and FPMB. Using HEASOFT v6.16 and calibration files with time stamp 7/2/2015, spectra are extracted in the energy range 3.0 − 50 keV from a source region with radius 65 ′′ . Whereas Suzaku's high energy coverage is provided by a non-imaging instrument, NuSTAR's detectors generate images for their entire band-pass. NuSTAR's spectra above 10 keV have, therefore, a substantially lower background, and the background can be much better constrained by using an off-source part of the detector. The ratios of source to background count rates are between 4 and 10.

Light curves and hardness ratio
We investigate the time variability of the count rate and the hardness ratio (Fig. 1). For the latter we take the ratio of the count rate in the 6 − 10 keV and 3 − 6 keV bands. We show data from XMM-Newton's EPIC pn, CCD 0 of Suzaku's XIS, and NuSTAR's FPMA, which all are sensitive in the two mentioned energy bands. Because of the different collecting areas and energy responses of the instruments, the count rates and hardness ratios cannot be directly compared between the different instruments. For easier comparison of the count rates, we scale the count rates such that they are comparable to XMM-Newton observations II-V. The count rate of XMM-Newton observation I is scaled up, since its effective area was reduced due to pile-up. Suzaku count rates are multiplied by a factor 14.3 such that VI matches I, as the two have similar X-ray fluxes (Section 3). We scale NuSTAR by a factor 16.7, such that its count rates match on average with Suzaku during the period when observations IX and X overlap with VII. Because of the different energy responses of the instruments, these scaling corrections are only approximate. For the same reason, no scaling is applied to the hardness ratios (nor are they derived from scaled rates).
A 1 hr time-scale is used for XMM-Newton and 2 hr for Suzaku and NuSTAR, because of their lower effective areas. The light curves exhibit some time variability. The largest variation in the count rate is by a factor 4.9 for observation VIII. The strongest variability is observed at lower count rates, where the uncertainties are correspondingly larger. Most importantly, the significance of the variations in the hardness ratio within each observation is low. This indicates that the spectral shape does not change substantially within each observation. This is consistent with previous studies, that found that the average spectra of each observation can be well fit (e.g., Parker et al. 2014). We, therefore, refrain from splitting observations for a fluxresolved study (e.g., Parker et al. 2014).

Bolometric Eddington ratio
We are interested in how the properties of the accretion disc change as a function of the mass accretion rate, and we use the bolometric luminosity as a tracer of the accretion rate. From the best-fitting model to a spectrum, the 2 − 10 keV luminosity in the rest frame of the host galaxy is calculated using a cosmological redshift of z = 0.025785 (Huchra et al. 1999) and assuming a flat universe with Hubble constant H0 = 70 km s −1 Mpc −1 and cosmological constant Λ0 = 0.73. As most of the disc's thermal emission is outside the X-ray band, we multiply the X-ray flux by a bolometric correction factor. We use the empirical relation derived by  Table 1. They are presented back-to-back, and are not always in time-order. Count rates are scaled for easier comparison (see Section 2.2), but H cannot be compared between instruments because of the different energy responses. Variations in H are small in each observation compared to the uncertainty, so the spectral shape does not evolve significantly during one observation, and we can fit the average spectrum in each observation. (2010), which provides the bolometric correction for the unabsorbed 2 − 10 keV luminosity depending on the shape of the continuum, specified by the photon index Γ of the underlying power law. Under the assumption of isotropic emission, this provides the bolometric luminosity, L bol . We express the luminosity as the Eddington ratio: the ratio of L bol to the Eddington limited luminosity, L Edd . The latter is determined as L Edd = (2.1 ± 0.5) × 10 45 erg s −1 using a black-hole mass of (14 ± 4) × 10 6 M⊙ (Peterson et al. 2004, see also Grier et al. 2012a) as well as a solar hydrogen abundance of X = 0.71.

Zhou & Zhao
The Eddington ratio's uncertainty is dominated by systematics. The relation of Zhou & Zhao (2010) introduces a substantial systematic uncertainty from the spread in the observations that they employed. Including this uncertainty and the uncertainty in the Eddington luminosity, the mean 1σ uncertainty in log L bol /L Edd is 0.77. Without these systematic uncertainties, and only including the uncertainties in Γ and the 2 − 10 keV luminosity, the mean error in log L bol /L Edd is 0.10. As the in-band flux is well measured, the uncertainty is mostly set by the error in Γ. We are interested in the correlation between parameters of a single source, so we omit the systematics from the uncertainty reported in each data point. The systematic error would only shift the observed trends with Eddington ratio.

SIMPLE POWER-LAW FIT
Previous analyses of these observations required a complex spectral model that includes ionized absorption at energies 3 keV as well as photoionized reflection (e.g., Parker et al. 2014;. Before we apply a similar detailed spectral model in Section 4, however, we first fit the data with a simple power law. The observed correlation of the photon index Γ and the X-ray flux (Sarma et al. 2015), suggests that we can gain some insight into the qualitative behaviour of the source with a simple spectral model. Following Parker et al. (2014), we fit a power law to the spectra in the energy intervals 3 − 4, 8.5 − 10, and 40 − 50 keV. This avoids absorption at low energies as well as the Fe Kα line and the Compton hump, which are prominent reflection features. These fits provide a relatively clean look at the underlying power law. We include absorption by the Galactic interstellar medium using the the Tübingen-Boulder model (Wilms et al. 2000) with an equivalent hydrogen column of NH = 3.6 × 10 20 cm −2 (Kalberla et al. 2005), and we use solar abundances as reported by Grevesse & Sauval (1998). A cosmological redshift is applied (Section 2.3). The spectral fits are performed using XSPEC v.12.8.2 (Arnaud 1996), and all presented errors in the fit parameters are at 1σ.
For each spectrum we determine Γ and the 3 − 10 keV power-law flux. We reproduce the logarithmic relation of Γ and the X-ray flux with a break at higher fluxes (Sarma et al. 2015), which has also been observed in other AGN (e.g., Perola et al. 1986;Done et al. 2000;Shih et al. 2002), and we will quantify this relation with the detailed spectral analysis in Section 4. The behaviour of the spectral components that were not explicitly modelled are visualized by the ratio of the spectra to their best-fitting power law (Fig. 2). This provides a consistent picture of how the spectra evolve as a function of flux. The Fe Kα line and the Compton hump are visible as excesses above the power law. Considering, for example, the Suzaku spectra (VI-VIII), the strength of the line is reduced at higher flux, which is expected for a more highly ionized disc. At the low-energy side there is a strong soft excess at E 1 keV as well as absorption in the 1 − 3 keV band. The relative strength of the soft excess appears largest at low fluxes. Also, absorption is most prominent at the lowest fluxes, whereas at the highest fluxes the soft excess extends up to 3 keV. Furthermore, the data points near 50 keV are consistent with the power law, except at the lowest flux values. This suggests that at low flux, the power law has a high-energy cut-off at an energy 10 -1 10 0 10 1 10 2 E (keV) F 3−10 (10 −11 erg s −1 cm −2 ) Γ Figure 2. For all observations identified with a Roman numeral (see Table 1) we indicate the 3 − 10 keV power-law flux F 3−10 , the power-law index Γ, and we show the ratio of the spectra to the best-fitting power law model. The different observations are off set, and horizontal lines indicate unity. We omit instruments with overlapping energies bands: the Suzaku XIS back-lit CCD and NuSTAR's FPMB. Clearly visible are the Fe Kα line around E = 6.4 keV, the Compton hump at E 15 keV, as well as a soft excess and absorption at low energies (E 3 keV).
near ∼ 50 keV, whereas at higher flux there is no cut-off near the considered energy band.

DETAILED SPECTRAL ANALYSIS
Although the picture sketched above is crude, it gives us the confidence that correlations exist between various parameters and the X-ray flux, even when ignoring the absorbed part of the spectrum ( 3 keV; see Sarma et al. 2015 as well as the discussion in Section 5.3). Next we use a more sophisticated spectral analysis to unveil these correlations. First we describe the detailed spectral model, which accounts for photoionized reflection off the disc. Next, we detail our strategy to mitigate model degeneracies, which forces consistency in the results of the different spectra. Finally, the fit results are presented, and we quantify the correlations of log ξ and other model parameters with flux.

Spectral model
The X-ray spectrum of Mrk 335 is rather complex, and consists of several components (e.g., . Thermal emission of the accretion disc dominates the spectrum in the optical and UV. This emission is reprocessed by Compton scattering off hot electrons in a corona to produce a power law that extends well into the (hard) X-ray band. In turn, the power law is reprocessed by scattering off the photoionized accretion disc. The reflection spectrum includes a strong fluorescent Fe Kα emission line close to 6.4 keV, accompanying absorption edges at slightly higher energy, and a "Compton hump" at energies 15 E 40 keV. Furthermore, ionized gas along the line of sight produces strong absorption below 3 keV. Depending on the flux level, as many as three separate absorption components have been detected simultaneously in this source (Longinotti et al. 2013).
We employ the spectral model relxill v0.2h Dauser et al. 2014), which combines relline (Dauser et al. 2010) and xillver (García et al. 2013). Xillver provides the illuminating power law with photon index Γ and a high-energy exponential cut-off at energy E cutoff . Furthermore, it includes the photoionized reflection variant of the power law for an accretion disc with iron abundance AFe relative to the solar value and ionization parameter ξ. In this paper we use log ξ, with ξ in units of erg s −1 cm. The relative strength of the illumination and reflection components is expressed as a reflection fraction, which is calculated in the 20 − 40 keV energy range.
The second part of the relxill model, relline, describes smoothing of the reflection spectrum due to relativistic effects, such as gravitational redshift and relativistic Doppler broadening. It takes into account the black-hole spin, a, the inclination angle of the accretion disc with re-spect to the observer's line-of-sight, and the disc's emissivity profile from an inner radius up to an outer radius.
A second reflection component likely originates at large distance from the black hole, and its main contribution to the spectrum is a narrow Fe Kα line. Instead of including a full reflection model, we describe the line with a narrow unresolved Gaussian profile at 6.4 keV, redshifted to the reference frame of the host galaxy. This simplified model is tested using NuSTAR observation XI, which has a clear detection of the narrow Fe Kα line and broad energy coverage. Replacing the Gaussian line by a xillver component with log ξ = 0, we find that outside of the Fe Kα line it contributes at most 13% of the total flux at any given energy. The best-fitting values of the model parameters are changed by at most 0.3 σ, and a similar goodness of fit is obtained.
Ionized absorption is significant at energies below ∼ 3 keV, and the simple analysis in Section 3 indicates that correlations of the spectral properties and the X-ray flux can be retrieved even when ignoring the absorbed part of the spectrum (see Section 5.3 for further discussion). Furthermore, prominent reflection features, such as the Fe Kα line and edge as well as the Compton hump, are at higher energies, so the reflection spectrum can be constrained without the low-energy part. For simplification we, therefore, ignore the spectrum below 3 keV, and omit spectral components for the soft excess and ionized absorption. Absorption by an ionized gas may introduce spectral lines above 3 keV. Where needed, we include narrow Gaussian profiles to fit these lines.
We apply the same cosmological redshift and absorption by the Galactic interstellar medium as in Section 3. For the absorption model we use solar abundances that match those used in calculating the xillver models (Grevesse & Sauval 1998).

Reducing model degeneracies
When the above-described model is applied to spectra of limited quality, there are multiple solutions with similar goodness of fit. For example, below 3 keV the uncertain nature of the soft excess and the strong absorption pose a challenge for determining the power-law parameters. We, therefore, restrict the fits to E > 3 keV. Furthermore, the ionized and the distant reflectors can be hard to distinguish by a fit: we include a full description of the dominant component with ionized reflection and from distant reflection only the narrow Fe Kα emission line.
Another source of degeneracy in the model parameters is the Fe Kα line (e.g. Dauser et al. 2013). Because of the strong dependence of the line's strength, shape, and energy on the ionization state of the disc, it is crucial for determining log ξ. Its properties are, however, also influenced by the reflection fraction, the disc's iron abundance, AFe, and the relativistic effects described with relline. The latter depend on the the black-hole spin a, the inclination angle, and the disc's emissivity profile. In total a combination of eight parameters determines the properties of the iron line. Consequently, multiple combinations of parameter values produce a similar goodness of fit. Indeed, significantly different parameter values have been obtained for individual spectra of Mrk 335 (e.g., , whereas one would not expect, e.g., AFe or the inclination angle to change between observations. To break these degeneracies, we constrain several parameters. We fix the inclination angle to 30 • , as one typically expects a relatively small angle for Seyfert 1 galaxies (the fit results are insensitive to the precise value). We take the disc emissivity to decrease with the third power of the radius (e.g., Miniutti et al. 2003). The outer radius of the disc is set to a large value of 400 rg, and the inner radius is placed at the inner-most stable circular orbit, which depends on a. We expect a as well as AFe to be the same for all observations. Therefore, we perform a combined fit of all XMM-Newton spectra and tie a and AFe between the different spectra, whereas all other free parameters are different for each spectrum. The XMM-Newton spectra typically have the highest quality, and including more spectra leads to an exceedingly problematic fitting procedure. We find a = 0.89 ± 0.05 and AFe = 3.9 ± 0.7. When fitting the Suzaku and NuSTAR spectra, we fix a and AFe to these values.
This procedure substantially reduces the model degeneracies. Furthermore, as the fixed and tied parameters have the same values across all spectra, correlations among the other parameters are more easily retrieved.

Fit results
Using the procedure described above, we fit all spectra. For back-to-back observations with the same instrument where the spectral shape does not change significantly (Fig. 1), we combine the spectra. This is the case for XMM-Newton pointings IV and V, as well as Suzaku observations VII and VIII. NuSTAR's IX and X were observed simultaneously with Suzaku's VII (Fig. 1), and we fit IX and X each together with the parts of VII that overlap in time. The overlap is only a small part of the total exposure time of VII. For all spectra we obtain good fits with χ 2 ν ≤ 1.17 (Fig. 3). Unresolved red-shifted Gaussian profiles model spectral lines when needed, and we find emission lines at (7.01 ± 0.04) keV (in the rest frame of the host galaxy) for II and (7.60 ± 0.05) keV for VII and VIII as well as an absorption line at (8.2501 ± 0.02) keV for IV and V. The presence of these lines does not significantly influence the goodness of fit.
For all observations we determine the best-fitting values and 1σ uncertainties of all model parameters as well as the flux in the 2 − 10 keV band (Table 2, Fig. 3 and 4). We estimate the bolometric correction depending on Γ (Fig. 3), and calculate the Eddington ratio L bol /L Edd (Fig. 4; Section 2.3).
The values of Γ that we obtain are consistent with those from the simple power-law fit in Section 3 (Fig. 5). Importantly, this means that the simple fits were sufficient for constraining the underlying power-law continuum, and that the qualitative behaviour of deviations from the power law (Fig. 2) are well described by the detailed reflection model. Therefore, the reflection model can quantify this behaviour. For example, the changes in the Fe Kα line apparent from the power-law fits (Fig. 2) now translate into an evolving log ξ measured with the reflection model (Fig. 3, 4).
At high Eddington ratio, log L bol /L Edd −1 (equivalent to log F2−10 −11.3), no high-energy cut-off of the power law is found, and we use a value of E cutoff = 300 keV, which is far outside the energy band of our spectra. 2 The fits of those spectra are insensitive to the precise value of E cutoff : taking E cutoff = 1000 keV changes the best-fitting values of other parameters by less than 1 σ. At lower Eddington ratios, however, we measure E cutoff within the energy band of Suzaku and NuSTAR. Investigation of confidence regions in χ 2 space around the best fits shows that E cutoff is wellconstrained, and does not exhibit strong degeneracies with, e.g., the reflection fraction (see Appendix A). The lowest value, for observations VII and VIII, is at 20 keV. This is the lowest value provided by the xillver model, and the actual value may therefore be somewhat smaller. XMM-Newton's band-pass does not allow us to constrain E cutoff . This may be an issue for the XMM-Newton observation with the lowest flux: III. The low flux, however, produces relatively large uncertainties, such that any deviations in the parameter values due to an incorrect E cutoff are not significant.
The non-ionized reflector is likely located at a large distance from the black hole, e.g., the outer region of the disc, a torus, or the broad-line region (e.g., Nandra 2006). The timescale for changes to propagate to this region is ∼years: similar to the interval covered by the observations considered in our study (Table 1). We, therefore, expect the normalization of the narrow Fe Kα line, KGauss, to be constant. Indeed, for four observations where KGauss is largest, its weighted mean 2 For spectra with a higher signal-to-noise ratio, simulations find that E cutoff can be measured even if it is outside the instrument band (García et al. 2015).  Order of the observations with log L bol /L Edd is the same as with flux in Fig. 3, except that log L bol /L Edd is largest for Suzaku's VI, followed by XMM-Newton's I (with the largest errors) and II.
is (4.2 ± 0.2) × 10 −6 c s −1 cm −2 keV −1 (Fig. 4). For all other observations KGauss is lower, which may be attributed to difficulty in separating this line from the one in the ionized reflection component. There is no evident correlation with Eddington ratio: the correlation coefficient is r = 0.29 with a p-value for no correlation of 0.57, and a linear fit yields a slope that is consistent with 0 within 1σ.

Correlations with Eddington ratio
We quantify the correlations between the logarithm of the Eddington ratio and log ξ, Γ, and the logarithm of the reflection fraction. For log ξ we also investigate the correlation with the logarithm of the 2 − 10 keV flux. The correlation coefficient, r, and the p-value for the null hypothesis (no correlation) are determined. Furthermore, a linear function is fit to the two parameters, taking into account the uncertainties in both, and we give the slope and intercept in Table 3. Often the goodness of fit, χ 2 ν , is substantially larger than unity, due to the small errors in certain data points. In those cases we include a systematic error, which is added in quadrature to the uncertainties of the data points (in the "vertical" direction), such that a "perfect fit" with χ 2 ν = 1 is obtained.
Γ correlates strongly with log L bol /L Edd . We fit a line for log L bol /L Edd < −0.5 . The data points at highest log L bol /L Edd are below this trend, which is consistent with the previously observed break in the trend at high flux (e.g., Sarma et al. 2015). A strong correlation is found between log ξ and the Xray flux, log F2−10, as well as with the bolometric Eddington ratio. The χ 2 ν values of the linear fit are, however, substantially larger than 1. The addition of a systematic error of ∼ 12 % yields χ 2 ν = 1 for both the flux and the Eddington ratio. The obtained slope is consistent within 1 σ with the previous fit. Interestingly, the slope of the linear fit is consistent within 1σ with the measurement of Ballantyne et al. (2011), who combined literature values of 10 AGNs. Also, the intercept is consistent within 2σ (1σ when including the extra systematic error; Table 3). This suggests that the measured relation may be a universal property of the behaviour of these sources. Our measurement is not affected by systematic uncertainties due to differences between sources, and yields a higher correlation coefficient. Therefore, we firmly establish the presence of the correlation of log ξ and flux as well as Eddington ratio for Mrk 335.
Although the reflection fraction clearly appears to decrease for increased values of the Eddington ratio (Fig. 4), the correlation is weaker because of substantial scatter in the data points. To obtain a χ 2 ν = 1 fit, a high systematic error of 94% is needed, and the correlation is no longer significant.
The cut-off energy of the power law transitions from a low value of ∼ 30 keV for log L bol /L Edd < −1, to E cutoff 300 keV at higher Eddington ratio. The transition is, however, not sampled well enough to quantify the correlation.

DISCUSSION
We have performed a multi-epoch analysis of all available XMM-Newton, Suzaku, and NuSTAR observations of Mrk 335 to study the behaviour of the accretion disc corona as a function of Eddington ratio, by inferring correlations of Table 3. Correlations of the ionization parameter log ξ, photon index Γ, and reflection fraction with respect to log L bol /L Edd , and for log ξ also with respect to log F 2−10 (first row). Provided are the correlation coefficient r and the p-value for the null hypothesis (no correlation), as well as the slope, intercept, and χ 2 ν from a linear fit. Furthermore, we include linear fits with an extra systematic error in the parameter measurements such that χ 2 ν ≡ 1. the flux and the spectral parameters. Detailed models of ionized reflection can present degenerate solutions when applied to spectra of limited quality. Other authors have found acceptable fits to the individual spectra of Mrk 335 with alternative models, and found values of Γ, ξ (e.g., , and E cutoff (e.g., Parker et al. 2014) that differ significantly from the trends that we uncover. We optimized our fitting procedure explicitly to infer consistent behaviour across multiple observations at different flux levels. This allows us to quantify strong correlations of the Eddington ratio and Γ (see also Sarma et al. 2015), ξ, as well as correlations with the reflection fraction and E cutoff .

The ionization parameter-Eddington ratio correlation in Mrk 335
The previous section showed that Mrk 335 exhibits a strong and statistically significant correlation between ξ and L bol /L Edd . The Ballantyne et al. (2011) α-disc prediction of this relationship, assuming a radiation-pressure dominated accretion disc and a geometrically thick corona (i.e., H/R = 1, where H is the height of the X-ray source above the disc and R is the radial distance along the disc 3 ) is where η is the radiative efficiency of the disc, f is the fraction of accretion energy dissipated in the corona, α is the viscosity parameter (Shakura & Sunyaev 1973), and (RR, Rz, RT ) encompass the general relativistic effects and are dimensionless functions of a and (R/rg) (e.g., Novikov & Thorne 1973). The spectral analysis of Mrk 335 found good fits when the inner radius of the reflecting region was at the ISCO of a black hole with a = 0.89 (i.e., R = RISCO = 2.39 rg). This removes some of the freedom present in comparing the observations to the prediction of Eq. 1, and potentially eases the physical interpretation of the results. However, Eq. 1 becomes inaccurate for R close to RISCO, as the functions Rz, RT and RR begin to diverge (Krolik 1999). This causes the disc density to increase rapidly, and ξ therefore drops precipitously at radii close to the ISCO (Fig. 6). Such behaviour is unphysical and is not seen in numerical simulations of accretion flows close to the ISCO where the gas density is found to fall quickly as material begins to plunge towards the event horizon (e.g., . Thus, unlike what is seen in Fig. 6, ξ should continue to rise as R moves toward RISCO. To correct this issue, the run of ξ with L bol /L Edd is first computed at R = 4 rg where (RR, Rz, RT ) have not begun to strongly diverge (Fig. 6), and then this result is scaled to R = 2.4 rg using the (R/rg) −7/2 scaling predicted from Eq. 1. This procedure results in the two lines seen in Fig. 7, where the model predictions are plotted over the measurements of Mrk 335. The dashed line assumes a constant f = 0.45 (Vasudevan & Fabian 2007) and the solid line assumes f ∝ (L bol /L Edd ) −0.77 (Stella & Rosner 1984). A constant α = 0.1 and η = 0.1 (Davis & Laor 2011) is also assumed for both these curves. It is remarkable that the straightforward theoretical prediction of Eq. 1 provides such a decent description of the observations of log ξ from Mrk 335. This basic level of agreement suggests that our simple understanding of accretion discs and coronae is on the right track. Yet, it is also clear that the theoretical slope is too steep, with ξ underpredicted at low L bol /L Edd and overpredicted at large L bol /L Edd . This result was also found by Ballantyne et al. (2011) using literature results, and indicates that Eq. 1 is missing one or more important physical effects. Ballantyne et al. (2011) speculated how changing some of the quantities in Eq. 1 (e.g., α) or the underlying disc density may help bring the predicted ξ more in agreement with the measurements, but there were no additional constraints with which to narrow down the possibilities. However, our spectral fits of Mrk 335 provide significant new information that may allow for a clearer picture of the corona in this system.

Piecing it together: a vision for the corona of
Mrk 335 Fig. 4 shows that there is evidence for an anti-correlation between the reflection fraction and L bol /L Edd in Mrk 335. The strongest reflection is found when the source is in the lowest flux state, which is consistent with light-bending enhancing the irradiation of the disc from a compact corona (e.g., Miniutti & Fabian 2004). Our spectral model indicates a disc that is strongly irradiated at the ISCO of a spinning black hole, so, depending on the height of the X-ray source, light bending could be important in illuminating the disc. The enhancement of the irradiating flux due to lightbending is not included in the predictions of Eq. 1, so this is a natural explanation for why log ξ is underpredicted at low Eddington ratios -the flux is too low and therefore the disc is not as ionized as observed. However, this effect would cause ξ to increase at all Eddington ratios, exacerbating the overprediction at large L bol /L Edd , unless there are changes to the coronal geometry that reduce the irradiating the flux at these Eddington ratios. The fact that the reflection fraction seems to drop at higher values of L bol /L Edd indicates that some change in coronal geometry is occurring. The trends observed in Γ and E cutoff (Fig. 4) together support this interpretation. The primary X-ray power-law is thought to be generated by thermal Comptonization in the corona (e.g., Stern et al. 1995;Poutanen & Svensson 1996;Petrucci et al. 2001), and depends on the coronal optical depth and temperature; e.g., Γ ≃ 9 4 + mec 2 kTeτ (1 + τ /3) where τ 1 is the optical depth and kTe the electron temperature (e.g., Shapiro et al. 1976). The upper two panels of Fig. 8 show estimates of kTe and τ (from Eq. 2) as a function of L bol /L Edd for Mrk 335. The coronal temperature is estimated using kTe = E cutoff /3 for L bol /L Edd < 0.1 and kTe = E cutoff /2 for larger Eddington ratios (e.g., Petrucci et al. 2001). These plots clearly suggest that the coronal optical depth diminished by about an order of magnitude as the Eddington ratio increased in Mrk 335. Once the Eddington ratio exceeded ∼ 10 %, the optical depth and the temperature remained roughly constant. One implication of this behaviour is seen in the lower panel of Fig. 8. Here, the Compton y-parameter, defined as (Petrucci et al. 2001) y ≃ 4 4kTe is plotted versus L bol /L Edd . For a given geometry, a Comptonizing corona in equilibrium will adjust so that it has a fixed y (e.g., Stern et al. 1995). The y parameter of the Comptonizing corona of Mrk 335 changes substantially with Eddington ratio and is only relatively constant at larger Eddington ratios (Fig. 8). These results strongly imply that the coronal geometry of Mrk 335 evolves from a high optical depth configuration to a lower one as the Eddington ratio increases. Further evidence for a change in coronal geometry can be seen by considering the compactness of the X-ray source, which is proportional to the ratio of the X-ray luminosity and the size of the corona (e.g., Guilbert et al. 1983;Fabian et al. 2015). Considerations of energy balance and radiation processes in the corona lead to various relation-ships between the compactness and kTe, all of which point to an inverse relationship between the two; that is, independent of the scattering processes ongoing in the corona, the compactness is smaller for larger temperatures. Observations by NuSTAR and other missions are now showing that multiple AGNs follow such a trend . Assuming the radiative processes ongoing within the corona do not change, the fact that the temperature rises with Eddington ratio in Mrk 335 implies that the compactness will fall. Since the luminosity is also increasing with a larger L bol /L Edd , then the compactness can only fall if the size scale of the radiative source increases faster than the luminosity rise.
Putting all these pieces together paints the following portrait of the evolving corona of Mrk 335. At L bol /L Edd < ∼ 0.1 the corona is warm, compact, optically thick and probably situated close to the black hole. This leads to significant light bending and illumination of the inner disc, increasing the reflection fraction and ionization state of the reflecting region. As the Eddington ratio increases, the corona expands and its optical depth drops, allowing the temperature to rise. The reflection fraction also begins to fall indicating that the disc is not being irradiated as strongly as before; this argues that the change in coronal geometry is more in the vertical direction away from the disc, rather than covering more disc area which would normally increase the reflection fraction. If the X-ray source is now further from the black hole, light bending would not be as important, and the flux on the disc surface would be significantly reduced. In that way, the predicted ξ can be brought closer to the observations at L bol /L Edd > ∼ 0.1.

Soft excess and absorption
The ratio of the spectra to simple power-law fits exhibits an evolving soft excess as well as absorption at low energies (E < 3 keV in Fig. 2). Previously, some studies have interpreted the soft excess as part of the reflection signal (Crummy et al. 2006), since reflection models predict a bump in the spectrum due to bremsstrahlung and a multitude of emission lines (e.g., García et al. 2013). Although the peak is at low energies, E 10 −2 keV, the tail extends into the XMM-Newton EPIC pn and Suzaku XIS bands. At large values of log ξ the feature is weak, but at low log ξ its shape is similar to the soft excess.  used the reflionx reflection model (Ross & Fabian 2005) to fit the XMM-Newton and Suzaku spectra of Mrk 335, and included the data at E < 3 keV. As the X-ray count rate peaks in the soft band, it dominates the fit, and the soft excess likely forces small values of log ξ. When they repeated their analysis with similar reflection models as in our study, no good fits were obtained. As these models present an improved description of reflection spectra, it suggests that the bump at low energies is not fully described by the reflection signal. By ignoring the soft band in our fits, our measurement of log ξ depends more strongly on the Fe Kα line and (for Suzaku and NuSTAR) the Compton hump. We find that log ξ correlates with flux: as expected, stronger illumination produces a higher log ξ (Fig. 3). For a cleaner view of the soft excess and absorption, we take the ratio of the spectra to the best-fitting reflection models from Section 4 ( Fig. 9). We show the ratio for XMM-Newton observations down to F 3−10 (10 −11 erg s −1 cm −2 ) Figure 9. Ratios of XMM-Newton and Suzaku spectra and the best-fitting reflection models, similar to Fig. 2 (the flux values are taken from that figure). The spectral fits were performed for E > 3 keV, whereas here we illustrate deviations from the bestfitting models at lower energy due to absorption and a soft excess.

keV and for
Suzaku XIS to 0.6 keV, whereas NuSTAR is not sensitive below 3 keV. The qualitative picture is the same as detailed in Section 3: absorption is stronger at low flux (see also Longinotti et al. 2013), and a soft excess is present. Reflection, therefore, does not fully account for the soft excess. Furthermore, at low flux, observations III, VII, and VIII exhibit an additional feature just below 1 keV, which may signify the presence of an absorption edge. Previous studies found absorption columns of up to NH ∼ 10 23 cm −2 (Grupe et al. 2012;Longinotti et al. 2013), which could be large enough to impact spectral parameters such as Γ and log ξ. These studies, however, typically did not include a reflection component, and obtained poor fits when including E 2 keV data. As a test, we take into account warm absorption when fitting spectrum III, which may have relatively strong absorption as the source is in a low-flux state. Columns of NH ≤ 10 22 cm −2 do not significantly change the fit results, whereas for NH = 10 23 cm −2 the predicted absorption is stronger than observed. Furthermore, Sarma et al. (2015) analysed the XMM-Newton spectra of Mrk 335, and obtained qualitatively and quantitatively similar results when either fitting the E > 3 keV data or when including an ionized absorption component in fits of the full instrument band. This supports our assumption. If one nevertheless includes a strong absorption component in the spectral model, the power law component compensates with a larger photon index. For example, NuSTAR observation XII exhibits a relatively high flux, where we do not expect substantial absorption (Fig. 9).  analyse this pointing in combination with contemporaneous Swift XRT data while including strong absorption (NH ∼ 10 23 cm −2 ), and obtain a photon index that is 29 per cent higher than our best fitting value. Our results can form the basis of an improved study into the nature of the soft excess and ionized absorption, which may also be responsible for some the observed spectral lines (Section 4.3).
To quantify the strength of the soft excess we consider the ratio of the spectra to their best-fitting models at 0.6 keV, which is the lowest energy that we include for Suzaku. For the power-law fits, the ratio is largest at low Eddington ratio, and decreases with luminosity (grey points in Fig. 10). This is consistent with models of reflection spectra, that predict soft emission to be stronger for low values of ξ, which we find at low Eddington ratio (Fig. 4). Indeed, when we take the ratio of the spectra to the best-fitting reflection models, the ratios at low Eddington ratio are substantially reduced (black points in Fig. 10). They are, however, not reduced to zero, but a component remains. This component is present at both low and high Eddington ratio, and it is relatively constant with a mean value of 1.9. A constant ratio implies that it is tracks the overall flux, rather than, e.g., the reflection fraction. It may originate from the source of seed photons for coronal X-ray emission: the accretion disc. For example, the source could be a warm scattering disc atmosphere (e.g., Magdziarz et al. 1998;Done et al. 2012;Petrucci et al. 2013;Różańska et al. 2015).

CONCLUSIONS
Applying multi-epoch spectroscopy to XMM-Newton, Suzaku, and NuSTAR observations of Seyfert 1 galaxy Mrk 335, we have uncovered consistent behaviour of multiple spectral parameters as a function of flux and, therefore, as a function of the Eddington ratio. In particular, we find a strong correlation between the ionization state of the inner accretion disc and the Eddington ratio. The slope is consistent with a previous measurement for literature values of 10 AGNs (Ballantyne et al. 2011), suggesting that the relation describes behaviour that is common for these sources. We firmly establish the correlation with a higher correlation coefficient, as our result is not affected by systematic uncertainties between sources. Furthermore, this correlation is close to the theoretical relation predicted by Ballantyne et al. (2011), but the measured slope is less steep, suggesting changes in the geometry of the corona. We infer the following behaviour of the evolving corona. At low values of the Eddington ratio, the corona is optically thick, compact, and close to the inner accretion disc where light bending is strong. As the accretion rate (and the Eddington ratio) increases, the corona expands, likely in the vertical direction. Above ∼ 10 % of the Eddington ratio, the corona is more extended, optically thin, and light bending is less important.
Our fits to the spectra above 3 keV indicate that at lower energies the soft excess consists of 2 components. One is due to reflection and its strength depends on the ionization state, whereas the other does not. The strength of the latter component is proportional to the overall flux, and possibly originates from the accretion disc.
Compared to previous work that considered individual spectra, we find that multi-epoch spectroscopy is essential for breaking degeneracies in the spectral fits and for obtaining accurate spectral parameters. Furthermore, we show that this method provides a powerful tool to study coronal evolution. The rich archives of the previously mentioned observatories provide the opportunity to extend this investigation to include several other bright AGN, which will reveal whether the behaviour that we found is common or unique to Mrk 335. Finally, future theoretical work needs to investigate the evolution of the coronal geometry with Eddington ratio to explain the observed behaviour of the disc's ionization state.  Figure A1. Contours of the goodness of fit, χ 2 , around the best fitting values of two parameters: the power-law cut-off energy and the reflection fraction. From inside out, the concentric lines represent the 1 σ, 2 σ, and 3 σ confidence levels. Top: Suzaku observation VII (best-fitting result at domain boundary, 1 σ contour not visible); middle: NuSTAR observation IX combined with overlap in Suzaku VII; bottom: NuSTAR observation X with overlap in Suzaku VII.
generacy between the two parameters influences our results in the three cases where we could measure E cutoff (observations VII, IX, and XII; see Fig. 4). Contours of confidence intervals in χ 2 -space around the best fitting values exhibit only minor asymmetries (Fig. A1), and no evidence of strong degeneracies is visible. For Suzaku observation VII, the best fit lies on the domain boundary of E cutoff . Otherwise, our fit results are well localised. Furthermore, when we determined the uncertainties in the fit parameters, we took any asymmetries into account, such that our reported errors include the full extent of the contours in Fig. A1.
This paper has been typeset from a T E X/L A T E X file prepared by the author.