Changing Accretion Geometry of Seyfert 1 Mrk 335 with {\it NuSTAR}: A Comparative Study

We present a detailed spectral study of the narrow-line Seyfert~1 galaxy, Markarian~335, using eight epoch observations made between 2013 and 2020 with the Nuclear Spectroscopic Telescope Array. The source was variable during this period both in spectral flux and flow geometry. We estimated the height of the Compton cloud from the model fitted parameters for the whole observation period. This allowed us to investigate the underlying physical processes that drive the variability in X-rays. Our model fitted mass varies in a narrow range, between $(2.44\pm0.45-3.04\pm0.56)\times10^{7}M_\odot$, however, given the large error bars, it is consistent with being constant and is in agreement with that known from optical reverberation mapping observations. The disk mass accretion rate reached a maximum of 10\% of the Eddington rate during June 2013. Our study sheds light on mass outflows from the system and also compares different aspects of accretion with X-ray binaries.


Introduction
Active galactic nuclei (AGN) host an accretion disk with an energetic X-ray producing Compton cloud so-called corona around a supermassive black hole at the center. The emitted X-ray radiation also acts as one of the most direct probes of accretion onto AGN. They show rapid flux and polarization variations which indicate that the observed high energy radiation mainly originates from the inner few to tens of Schwarzschild radii (r s = 2GM BH /c 2 , M BH , G, and c are the mass of the black hole, gravitational constant and the speed of light respectively) from the event horizon of the black hole [1]. The power-law (PL) component is widely accepted as the effect of inverse Compton (IC) scattering of thermally produced soft photons from the accretion disk by a corona of hot electrons at the inner edge of the disk [2,3]. However, the geometry and size of the corona as well as the physical mechanisms governing the energy transfer between the disk and the corona are not well understood, which necessitates a physical model to better understand the observed signatures in them.
There are proposed models in the literature [4] which consider the corona to be the region between the truncation radius and the innermost stable circular orbit. However, it is not clear why the disk truncates at a certain radius from the central black hole and how the truncation radius is connected with the corona. Furthermore, most of the models in the literature [5,6] use optical depth or the coronal temperature or the spectral index as a parameter to compute the spectrum from the corona. None of them are the basic physical quantities of accretion. According to two-component advective flow model, TCAF [7,8], at a certain distance from the BH where gravitational force balances with the centrifugal force, shock forms by the low angular momentum, hot, sub-Keplerian halo, and satisfying Rankine-Hugoniot conditions, which is the region of the truncation of the disk. Apart from that as the two forces balances each other accreting matter virtually stops and its velocity goes down, which also increases the optical depth of the corona region. As the gravitational force dominates closer to the BH, velocity again increases and the flow passes through the inner sonic point and falls onto arXiv:2103.16793v1 [astro-ph.HE] 31 Mar 2021 the BH, therefore the flow is transonic [7]. Beyond this shocked region, both the disk and halo matter pile up to decide the optical depth and temperature of the corona region. The same region also upscatters the incoming soft radiation from the disk via IC. A typical X-ray spectrum of AGN in the 2-10 keV shows primarily the signature of PL. Due to IC, the corona cools down and its area decreases, leading to low brightness level and a softer spectrum [9]. These physical properties and the dynamics of the flow makes the TCAF more favourable as a physical model than other existing models in the literature. In addition to the above, TCAF uses mass accretion rates and the mass of the BH as model parameters, which are the basic physical quantities of a flow. Furthermore, the same shock location which can explain the spectral features, can also be able to explain the temporal features from its oscillation [10,11]. Therefore, among the various models available in the literature to understand the observed X-ray emission in AGN, we preferred to consider TCAF model for our study.
TCAF model has mainly five parameters: (i) mass of the black hole (M BH ), (ii) disk accretion rate (ṁ d ), (iii) halo accretion rate (ṁ h ), (iv) location of the shock (X s ), and (v) shock compression ratio (R). Increasingṁ d keeping the other parameters fixed, increases the number of soft photons from the Keplerian disk [12] making the spectrum softer. Higherṁ h increases the corona temperature and hardens the spectrum. Similarly increase of either X s or R lead to hardening of the spectrum. However, in reality all parameters change in multidimensional space. Therefore such variations may not sustain due to non-linear change in optical depth and coronal temperature in different epochs. A detailed parameter study has been made in Ref. [13]. In recent years, TCAF model is implemented in XSPEC and has been used widely to study both X-ray binaries (XRBs) [14][15][16] and AGN [17][18][19]. However, there are physical processes which are not yet incorporated in the current TCAF model e.g. jet, bulk motion Comptonization and the spin effect of the BH.
Apart from the above components, many AGN exhibit extra features in their X-ray spectra that likely originate from the accretion disk [20]. The so-called X-ray reflection features produced by the accretion disk might be because of its illumination from the central compact object, such as a magnetically dominated corona residing above the surface of the disk [e.g. 2,21,22]. A prominent feature in the X-ray spectrum is the Fe Kα emission line at around 6.4 keV which is due to fluorescence in a dense and relatively cold medium [23][24][25]. The shape of the Fe Kα line is a powerful probe of the general relativistic effects close to the BH [26,27], leading to estimates of the spin of the BH [28][29][30]. Indeed, observations over the last decades by X-ray satellites have uncovered relativistically broadened Fe Kα line in the spectra of several AGN [31][32][33] and black hole XRBs [e.g. 34]. Models of X-ray reflection spectra [5,35] show additional features that collectively can constrain the ionization state, heating, and cooling of the illuminated surface [36].
Mrk 335 is a bright narrow-line Seyfert 1 (NLS1) galaxy at a redshift z = 0.026 [37], with a black hole mass, M BH = 2.6 × 10 7 M [38] and showing ionized absorption in the X-ray and ultra-violet (UV) bands [39]. It has been observed to go into extremely low-flux states [40] where the soft X-ray flux dropped by a factor of up to ∼ 30, while the hard flux dropped only by a factor of ∼ 2 in 7 years. Using observations from XMM-Newton, Ref. [41] found signatures of absorption and high reflection fraction in the spectrum of Mrk 335 in its lowest flux state, whereas at the intermediate flux state the observed X-ray spectra is explained well by blurred reflection model without the requirement of variable absorption [42]. Ref. [43] estimated a high-frequency lag of ∼ 150 s between the continuum dominated energy bands and the iron line and soft excess components for this source. This time lag suggests that the continuum source is located very close to the central black hole and that relativistic effects are supposed to be worth considering. Using combined Swift and NuSTAR observations, Ref. [44] interpreted the X-ray flare during 2014 as arising from the vertical collimation and ejection of the X-ray emitting corona at a mildly relativistic velocity, that causes the continuum emission to be relativistically beamed away from the disk. Ref. [45] discussed the effect of the geometry of the corona on the relativistically blurred X-ray reflection arising from the accretion disk, which can also explain the variability in between low and high-flux states. Ref. [46] performed the structure function analysis using long term Swift optical/UV and X-ray data of Mrk 335 and showed that the X-ray low flux state could be due to the physical changes in the corona or absorption, whereas the variability in the optical/UV band is more consistent with the thermal and dynamic timescales associated with the accretion disk. Ref. [47], studied the low-flux state data using XMM-Newton and observed absorption lines from a highly ionized outflowing wind. Several studies on Mrk 335 clearly show that the system is highly variable leaving a range of possibilities in its dynamical behavior and emission features. Mrk 335 is known to show change in geometry [45], such a change in geometry can be found from a disk-based model fit such as TCAF to data. It is not clear to date if there is any change in the mass accretion of the source, and if so, what effects it can have on the geometry of the corona in the source? To investigate the above, in this work we chose this object and carried out an analysis of the X-ray data on the source acquired by NuSTAR during the period 2013-2020. This paper is organized as follows: in the next section, we describe the observation details and the analysis procedures. In §3, we explain our results obtained from model fits to the data and summarize our conclusions in the final section.

Observation and Data analysis
In the present manuscript, we analyzed archival data 1 of NuSTAR observations of Mrk 335 made during 2013-2020. During this time interval Mrk 335 was observed nine times, out of which in one observation, data quality is not good, so we considered the remaining eight observations. The details of the observations are listed in Table 1. The NuSTAR data in the energy range of 3.0−30 keV (as the data is noisy above 30 keV), were extracted using the standard NUSTARDAS V1.3.1 2 software. We ran NUPIPELINE task to produce cleaned event lists and NUPRODUCTS for generation of the spectra. We used a region of 80 for the source and 100 for the background using ds9. The data were grouped by grppha command, with a minimum of 30 counts per energy bin. The same binning was used for all the observations. For spectral analysis of the data we used XSPEC 3 [49] version 12.8.1. We used the absorption model TBABS [48] with the Galactic hydrogen column density fixed at 3.6× 10 20 atoms cm −2 [50] throughout the analysis.

Results and Discussion
We fit the NuSTAR observations of Mrk 335 taken between 2013 to 2020. During this observation period the source showed significant variability in spectral flux (see Figure 1). We did spectral analysis using different models. First, we performed spectral fitting using a simple powerlaw (PL) model as TBABS(GAUSS+PL), where the PL index (Γ) varies from 0.72-1.84, a change by a factor greater than two, which is clearly an indication of significant change in the flow behavior and corona properties. The model fitted parameters are provided in Table 2. Our model fits show that the line energy and width vary in a broad range, which point to complex geometry changes and gravity effects (more especially line broadening) in the source. In addition to the simple PL model, we also modeled the observed spectra using the complex phenomenological model PEXRAV [51] and the physical model TCAF [8]. The results of PEXRAV model fits are given in Table 3. Our PL model fits indicate that the source passed through hard and soft spectral states similar to XRBs, where the timescale for changing from one state to other is in the order of few days to weeks, which is the viscous timescale [52, and references therein]. We found a similar timescale when the corona changed from its maximum size to minimum, which is around 12 days. However, as the viscous timescale for AGN is in the order of few hundred to few Myr (scaled by the mass of the BH), it may not be possible to observe the similar timescale for other AGNs. In future we aim to explore this aspect of accretion to estimate the viscosity of the flow.
From PEXRAV model fits, we found the reflection fraction to be high with R ref > 0.9 for all the epochs. When the Fe abundance (A Fe ) was made as a free parameter, A Fe , was found to vary from sub-Solar to super Solar values during the epochs analysed here, which is quite unphysical. It has been pointed out by Ref. [53], reflection model fits to observed X-ray spectra are generally found to yield high A Fe values. Also, according to Ref. [54], the high A Fe obtained from model fits to data could be due to some unknown physical effects being overlooked in current models. Therefore, for PEXRAV model fits, we froze A Fe to solar value. [55]. Next we applied the physical model TCAF and having the form TBABS(TCAF+GAUSS) to extract the accretion flow parameters. For all observations goodness of the fit is more or less similar to PEXRAV model. The TCAF model fitted parameters are provided in Table 4. For some observations, better fits were achieved after adding some additive model components which are detailed in the table. For majority of the epochs, the χ 2 /do f obtained from PEXRAV and TCAF model fits are in agreement, which for three epochs, the χ 2 /do f from TCAF fits are marginally larger than those obtained from PEXRAV model fits. This might be due to differences in the spectral components between models e.g. TCAF does not include line features, that may give relatively higher χ 2 /do f . In Figure 2, we show different model fitted spectra for the OBSID 60001041002/3 in the left/right panels of the figure. Figure 3 shows the same spectral fitting but for the OBSID 60001041005 (left panel) and 80001020002 (right panel). The rest of the observational fits are shown in Appendix. We used TCAF model fitted parameters to estimate the geometrical height of the shock (H shk ). The height of the corona is basically the height of the shock, as the shock is the boundary layer of the corona. We estimated (H shk ) using the equation below [see also 14]: where, γ is the adiabatic index, X s is the location of the shock or the boundary of the corona, and R is the shock compression ratio (= ρ + /ρ − ), a ratio of the density of the downstream to the density of the upstream of the flow.
In Figure 4, we show the variation of various model parameters obtained from TCAF model fits as well as the variation of hardness ratio (HR) with time. The halo rate was always higher than the disk rate by an order of magnitude or more which implies that the flow was dominated by the low angular momentum matter, which might be accreted from the wind or host galaxies. Theṁ d was found to vary by a factor of ∼40 from 0.24% to 10% of the Eddington rate and theṁ h varied between 1.2-3.1 Eddington rate. The size of the shock location (X s ) or corona changed significantly and reached up to 6 r s during 25 June 2013 (60001041005) and 8r s during 2014 (80001020002). For the OBSID 60001041005, we required an extra broad (1.55 keV) Gaussian component at low energy 3.44 keV, which was also found to be broadened by strong gravitational effects [27]. For the OBSID 80001020002, we required an additive DISKLINE [26] component at line energy 6.12 keV with emissivity powerlaw index (β) ∼ −3.0 and R in = 24.4r g (GM/c 2 ) as the current TCAF model does not include line emission due to gravitational effect in the fit. It is worth noting that the X s is achieved for the highestṁ d , which might be due to the cooling effect inside the corona. As the accretion rate increased more disk photons got intercepted by the corona. This cooled the corona, causing the shock to move inwards and thus causing a significant change in the size of the corona. Such change in the corona could also be due to the activity of the jet, which extracted a significant amount of thermal energy and contracted the corona [58]. Therefore it is also be possible to establish a jet-disk connection using TCAF model fits to AGN spectra similar to XRBs [59-61, and references therein].
Such movement of the shock can give rise to the observed variability in spectral and temporal properties [62]. During 2014,ṁ h further increased and made the corona hotter and shock the receded slightly. During 2018, shock further moved away, the corona has become bigger in size, therefore intercepts more soft photons, generates high energy photons, which can also ionize the disk atmosphere. For this particular observation, we required partial covering fraction model (ZXIPCF), with a covering factor ( f c =0.45) and low ionization parameter (logξ=1.2). For the OBSID 90602619004, during 2020, we required another narrow line (0.29 keV) component, when the shock again moved inward. The model fitted shock compression ratio also varied in a broad range, which also dictates the strength of the shock       Figure 4, shows the estimated height of the corona, which changes significantly for the whole observation period.
Mass of the BH is one of the parameters in TCAF. Therefore, keeping it as a free parameter during TCAF model fits to the observed spectra will yield M BH value for Mrk 335. We kept M BH as a free parameter during our fitting and we found M BH to vary in range between (2.44 ± 0.45 to 3.04 ± 0.56) ×10 7 M . Considering the large error bars the derived BH mass is consistent to be constant. Our derived M BH is consistent with that estimated by [38] using optical reverberation mapping observations, however, an order of magnitude larger than that obtained by [53]. Such a low mass found by [53] could be due to them not considering an expanded corona in their model fit to the data. Our results on M BH also shows that by fitting accretion disk based models such as TCAF to observed X-ray spectra one can estimate mass of BH in AGN. Also, we found that freezing M BH to some fixed value from the above range and redoing the fit has negligible effect on other derived parameters, such asṁ d ,ṁ h , X s , and R.
The viscous time scale in AGN is much longer compared to that of XRBs. At small time scales (in the order of few days or month) significant change in the accretion rate is not expected, similar is the case with the heating and cooling rates. Thus at such timescales, variability is expected to be less. This is in fact observed in the radio-quiet category of AGN and accretion disk based models can explain the observed flux variations in them. However, this might not be the case in the radio loud category of AGN, as relativistic jets play a dominant role in the flux variations observed in them. Even in those sources as the ejection of jets extract thermal energy from the corona, the change in the shock location can be drastic in them though the cooling is less. Results on this will be reported elsewhere.

Conclusion
In the present paper, we analyzed NuSTAR data of Mrk 335, observed between 2013-2020 to understand the dynamics of the flow and the accretion behavior of the source. We found that the source was significantly variable both in spectral flux, accretion geometry, and reflection fraction. We also found the geometry of the corona to change between epochs. Such a change can affect the reflection fraction and the spectral energy distribution of the illuminating radiation, and consequently the ionization rate. From the TCAF model fits we obtained the dynamics of the flow along with the geometry change of the corona and accretion rates. In our TCAF model fits, the line profiles due to relativistic effects were taken into account using an additive DISKLINE model where ever needed rather than using throughout for all observations. Among other models, even though the XILLVER model has relativistic and reflection effects incorporated, the accretion dynamics and the origin of corona are still lacking in the model. Therefore we did not take into account such models in our current study. Below we provide the key findings on the variable nature of the source:

1.
During 2013, the corona at the inner region has changed significantly from its elongated stage (∼ r s ) to destruction stage ∼ 6r s , consistent with Ref. [44]. As the compression ratio has changed and the corona has contracted significantly, it can be possible that jet/mass outflow was launched around June 13, 2013, therefore a significant amount of thermal energy has been extracted by the jet/outflow and the corona contracted.

2.
The observation during 2014, required a blurred reflection component along with TCAF to fit the spectra. This is quite natural as the corona contracted, the inner edge of the disk moved significantly inward, therefore the gravitational effect became dominant [26,27] and blurred the Fe Kα line. We also required a broad Gaussian line component at ∼ 2keV. 3.
The HR was roughly constant during 2013, and is also similar to that obtained during 2014.

4.
There is a significant change in Γ PL between 2013 and 2014 spectra, which can be due to sudden change in accretion rates. During this period disk accretion had increased by a factor of a few and also the size of the corona contracted significantly.

5.
The steepening of the emissivity profile of Mrk 335 indicates that the corona is compact for this source [44]. In our study, we found that the size of the corona is indeed small and compact during 25 th June 2013 and 20 th September 2014 (see Table 4). It is also noticeable that the height of the corona reduced significantly during these two epochs. 6.
During 2014, spectral flux between 3-30 keV changed/increased by a factor of ∼ 3 compared to 2013 and 2018. This could be due to an increase in accretion rates as well as the change in corona. As the accretion rate increased (in 25 th June 2013 and 10 th July 2018), the number of soft photons increased, thereby increasing the cooling rate, i.e reduction of more energy from the corona by the seed photons from the Keplerian disk [see theoretical aspects in 9]. It should be noted that though the shock location changed significantly, the other parameters triggered that change, mainly the increase in disk accretion rate, therefore the cooling rate. This also infers that not only the shock location but other parameters are equally important to explain the observed variability. 7.
During 2018 and onward, the corona and H shk again elongated. During this period the HR also increased. This also implies that there is a correlation between HR and the geometry of the corona. 8.
During 2018, to take into account the disk ionization effects along with TCAF, we required partial covering, ZXIPCF model to better fit the data. The model fit showed an absorption column density with N H = 6.8 × 10 23 cm −2 is present. The fit also required a low ionization with, logξ = 1.2 erg cm s −1 with partial covering fraction of 0.45. This added component also indicates the presence of mass outflow from the system, which is evident in the monitoring observations of the source in optical/UV and X-rays [63,64, and references therein]. 9.
The reflection fraction, R ref , is measured as the ratio of the photon fluxes from the blurred reflection and powerlaw continuum model components. From PEXRAV model fitting we found this value > 0.9. 10. The mass of the black hole which is kept as a free parameter is found to vary in a very narrow range (2.44-3.04)×10 7 M , and considering the error bars is consistent with a constant. This is in agreement with that of Ref. [38]. 11. Earlier studies on observed X-ray variability inferred the origin of flux variations to changes in the primary powerlaw continuum possibly exhibited through intrinsic variations in the corona, or possible changes in its size or location [42,44,65]. This is in agreement with our present findings.

Informed Consent Statement: Not Applicable
Data Availability Statement: We used archival data for our analysis in this manuscript. Appropriate links are given in the manuscript. For the details of the data fitting, one can directly contact to the first author.