An Analytical Fourier-Transformation Model for the Production of Hard and Soft X-Ray Time Lags in AGNs: Application to 1H 0707-495

The variability of the X-ray emission from active galactic nuclei is often characterized using time lags observed between soft and hard energy bands in the detector. The time lags are usually computed using the complex cross spectrum, which is based on the Fourier transforms of the hard and soft time series data. It has been noted that some active galactic nuclei display soft X-ray time lags, in addition to the more ubiquitous hard lags. Hard time lags are thought to be produced via propagating fluctuations, spatial reverberation, or via the thermal Comptonization of soft seed photons injected into a hot electron cloud. The physical origin of the soft lags has been a subject of debate over the last decade. Currently, the reverberation interpretation is recognized as a leading theory. In this paper, we explore the alternative possibility that the soft X-ray time lags result partially from the thermal and bulk Comptonization of monochromatic seed photons, which in the case of the narrow-line Seyfert 1 galaxy 1H 0707-495, may correlate with fluorescence of iron L-line emission. In our model, the seed photons are injected into a hot, quasi-spherical corona in the inner region of the accretion flow. We develop an exact, time-dependent analytical model for the thermal and bulk Comptonization of the seed photons based on a Fourier-transformed radiation transport equation, and we demonstrate that the model successfully reproduces both the hard and soft time lags observed from 1H 0707-495.


Introduction
It has often been noted that all black hole sources demonstrate similar patterns of Xray variability, with the nature of the variability dependent on whether the source is in the low/hard state or the high/soft state (see van der Klis 2006 for a review). The same general variability patterns are observed for black holes of all mass scales, ranging from Galactic black holes (GBHs) to supermassive active galactic nuclei (AGNs). For example, by investigating the power spectral densities (PSDs) for a large sample of GBHs and AGNs, McHardy et al. (2006) showed that the timescale, t break , associated with the location of the characteristic PSD break frequency scales as t break ∝ (M )(Ṁ /Ṁ E ) −1 , where M is the black hole mass,Ṁ is the accretion rate, andṀ E denotes the Eddington accretion rate. A similar mass scaling law was discovered by De Marco et al. (2013), who observed that the Fourier frequency associated with the onset of negative (soft) time lags, ν neg , scales as ν neg ∝ (Ṁ /Ṁ E )/(M ). These correlations suggest the universality of the underlying accretion process over many decades of black hole mass.
In this paper we investigate the structure of black hole accretion flows using studies of variability based on observations of Fourier X-ray time lags in AGN. Our results successfully replicate an important, and relatively new observation, while simultaneously demonstrating a similar mass scaling relationship to that described by De Marco et al. (2013). In our model, the observed mass scaling laws occur as a natural consequence of bulk and thermal Comptonization occurring near the black hole.

Fourier Time Lags
The idea of computing the time lags between two photon energy bands using Fourier analysis was first suggested by van der Klis et al. (1987), with subsequent development and application by Nowak et al. (1999). The calculation of the time lags is based on the development of the complex cross spectrum, C, defined by where denotes the circular Fourier frequency, ν F is the Fourier frequency in Hz, and S(ω) and H(ω) represent the Fourier transforms of the soft and hard photon energy light curves, corresponding to photon energies s and h , respectively. The asterisk in Equation (1) represents the complex conjugate.
The Fourier transforms S(ω) and H(ω) appearing in Equation (1) where the soft and hard photon energy light curves are denoted by s(t) and h(t), respectively. The two light curves are related to the time-dependent X-ray photon count-rate spectrum observed at the detector, F ( , t), via The corresponding Fourier time lag, δt, is then evaluated using the relationship where the phase lag, φ, is given by According to the sign convention adopted here, a positive time lag δt > 0 corresponds to a hard lag, in which the hard channel signal lags the soft signal for a given Fourier frequency (or inverse variability period). It can be shown that if an exact time delay is introduced between the hard and soft channel light curves h(t) and s(t) so that h(t) = s(t − ∆t), then the resulting time lag is δ(ω) = ∆t as expected.
It can also be easily deduced that the existence of time lags necessarily implies time variability in the source, since steady-state light curves would generate no phase lag, and therefore no time lag. Hence the detection of time lags from a given source necessarily implies the existence of variability in the X-ray signal from the source. This implies that the emission from a given source can be viewed as a superposition of continual steady-state emission, combined with episodic bursts or flashes, which produce time-dependent phenomenon such as the time lags. That is the approach we take here in our analysis of the emission from the narrow-line Seyfert 1 (NLS1) galaxy 1H 0707-495.

Time Lag Observations
Hard Fourier time lags have been detected in the X-ray emission from many AGNs over the past three decades, as discussed by Hasinger et al. (1986), van der Klis et al. (1987), Lyubarskii (1997), Crary et al. (1998), Böttcher & Liang (1999), Kotov et al. (2001), Arévalo & Uttley (2006), and Kroon & Becker (2016). However, a relatively new field of investigation of AGN X-ray variability involves the robust detection of soft lags from AGNs. First reported for 1H 0707-495 (Fabian et al. 2009), broader studies have revealed that soft lags are not uncommon, and are found in all AGN mass scales, ranging from low-mass (Mallick et al. 2021) to high-mass (De Marco et al. 2013). For example, De Marco et al. (2013) detected soft time lags in 15 of the 32 radio-quiet AGNs in their survey, and also determined to high accuracy ( 4σ) that the magnitude of the lags scaled in proportion to the black hole mass over 2.5 orders of magnitude. The universality of the scaling law suggests that the same underlying physical mechanism may be responsible for producing the soft lags in all AGNs, regardless of the black hole mass.
The time scales of the observed high-mass AGN soft lags are on the order of ∼ 10 − 100 seconds (e.g., Fabian et al. 2009;Emmanoulopoulos et al. 2011;Cackett et al. 2013;De Marco et al. 2013), which are on the order of multiples of gravitational radii of AGN, making them extremely useful probes for understanding the structure of the inner regions of the accretion flows, including temperature, density, and accretion configuration. It should be noted that soft Fourier time lags have also been discovered in stellar-mass sources (e.g., De Marco et al. 2015). However, there is a great deal of difficulty in using these sources to study the underlying physical processes due to the extremely small time scales.

Previous Models
To date, two primary scenarios have emerged regarding the mechanism responsible for producing the soft Fourier time lags observed from GBHs and AGNs. The first scenario centers on the reprocessing of radiation in a zone in which the electrons cool on a timescale comparable to the lag time. The second scenario, commonly referred to as the reverberation interpretation, has become somewhat of a consensus in the literature, and focuses on the reflection of hard emission produced near the black hole off a cool electron population located at a larger radius, causing down-scattering with a time lag created by the transit-time propagation delay. Malzac & Jourdain (2000) investigated the first scenario, in which emission produced in the corona is down-scattered over time by electrons in the corona that cool in response to the scattering of the radiation. They predicted that if the flares arise from this mechanism, then the lags should invert from hard to soft if the energy dissipation timescale is between 10 H/c and 100 H/c, where H is the size of the emission region. These authors presented illustrative numerical simulations, but made no attempt to fit the data for any source.
When soft time lags were first discovered in the X-ray emission from 1H 0707-495, Zoghbi et al. (2010) argued that the model of Malzac & Jourdain (2000), in which the spectrum is formed in a cooling corona, is inconsistent with the observed spectral variability around ∼ 1 keV. They proposed an alternative scenario in which the time lags are the result of spatial reverberation. In their model, the soft lags are formed when hard photons from the corona (generated near the event horizon) are reflected off the surrounding, cooler accretion disk, which extends from an inner radius of ∼ 1.3 R g to an outer radius ∼ 400 R g , where R g = GM/c 2 is the gravitational radius for the central black hole. Zoghbi et al. (2010) concluded that this type of reflection scenario naturally produces soft time lags with a magnitude comparable to the light travel time between the corona and the reflecting region of the disk. A time lag of ∼ 30 s would therefore correlate with a black hole mass of M ∼ 2 × 10 6 M , which is consistent with the uncertain mass of this black hole (e.g., Zhou & Wang 2005). Miller et al. (2010) analyzed relative lags between three energy bands and concluded that their results invalidated the reflection model of Zoghbi et al. (2010). As an alternative, they proposed that the soft lags in 1H 0707-495 are caused by a partially opaque reverberating region with a minimum radius of ∼ 1, 000 light-seconds from the black hole corresponding to ∼ 100 R g for a black hole with mass M ∼ 2 × 10 6 M . In Miller's model, this region has decreasing opacity with increasing energy and has scattering or reflection present in all Xray bands. Zoghbi et al. (2011) argue that Miller's results are invalid because they imposed artificial constraints on the model geometry.
Subsequent authors have further explored the application of the reverberation model to the interpretation of the observed soft lags. Gardner & Done (2014) investigated the NLS1 AGN PG1244+026 and proposed that reflection alone is unable to produce the soft lags without additional emission from the accretion flow as well as additional reprocessed emission. Emmanoulopoulos et al. (2016) used a method proposed by Papadakis et al. (2016) to search for evidence of reverberation lags in the PSD of 1H 0707-495. We note that the non-detection of interference fringes in the PSD was consistent with the model, since the observations had insufficient signal-to-noise, coupled with too much uncertainty regarding the shape of the underlying PSD of the direct emission, which made it impossible to rule out the null hypothesis. There have also been numerous studies examining the iron K features in the lag versus energy spectrum for AGNs and GBHs that provides supporting evidence for the reverberation lag hypothesis (e.g., Kara et al. 2016;Kara et al. 2019;Wang et al. 2021).

Fourier Transform Model
The ongoing controversy regarding the origin and interpretation of the soft time lags, as well as the extreme complexity of the existing theoretical models, has motivated us to develop a new model for the production of the observed hard and soft time lags observed from AGN. In this paper, we explore the possibility that the observed time lags may be the result of thermal and bulk Comptonization occurring in the freely-falling, quasi-spherical inner region (corona) of accretion flows onto supermassive black holes. In this scenario, the time-dependent part of the signal (that produces the observed time lags) is the result of the episodic emission of monochromatic photons generated near the black hole, which are then Comptonized by hot electrons in the quasi-spherical coronal region. Conversely, the steady-state X-ray spectrum in this model represents emission produced continuously by the optically thick, geometrically thin accretion disk. The episodic emission of the seed photons may be due to inhomogeneous "clumps" in the accretion flow (e.g., Merloni et al. 2006;Gutiérrez et al. 2021).
The theoretical approach adopted here is a novel one, based on the fact that the observed time lags are computed using Fourier transformed X-ray light curves from the source AGN at two different X-ray energies (Kroon & Becker 2016). This has motivated us to develop theoretical predictions for the time lags by directly solving the transport equation to determine the Fourier transform of the radiation distribution throughout the quasi-spherical coronal region near the center of the accretion flow. This method is convenient because it facilitates the direct computation of Fourier-based time lags, which can then be compared with observational data.
The time-dependent transport equation we solve includes terms describing bulk and thermal Comptonization, spatial diffusion, advection, and photon injection and escape. We demonstrate that the results obtained for the time lags as a function of Fourier frequency agree reasonably well with the observations for 1H 0707-495, without requiring reflection off the accretion disk, or any special geometrical constraints.
The remainder of the paper is organized as follows. In Section 2, we introduce the transport equation and discuss the dynamical model of the spherical accretion flow. In Section 3, we obtain the exact solution for the Fourier transform of the radiation distribution inside the spherical flow. In Section 4 we develop the formalism for computing the X-ray time lags based on the Fourier transform of the radiation distribution escaping from the outer boundary of the spherical flow. In Section 5 we use the new model to generate simulated time lags for the narrow-line Seyfert 1 galaxy 1H 0707-495, and we compare the simulated lags with the observational data. Finally, in Section 6, we provide an overview of our main results and their astrophysical significance.

Model Overview
In this paper, we explore an alternative model for the production of the observed hard and soft time lags from accreting black holes. The previously explored reverberation and reflection scenarios provide a natural mechanism for generating time lags via light transittime effects, in which the reverberation is a geometrical effect due to the propagation and scattering of photons through the physical space of radius. On the other hand, the new model developed here explores the complementary effect of thermal and bulk Comptonization in creating time lags due to reverberation in the energy space rather than the physical space, as in the previous models. We argue that Comptonization lags and spatial transit-time lags probably both contribute to the observed time lags, and there is an interesting duality between the two mechanisms, since the Comptonization occurs in the energy space and geometrical reverberation occurs in the physical space.
The model developed here involves the analytical solution of a time-dependent transport equation describing the scattering of seed photons injected into an accreting plasma cloud in the region of flow near the central black hole. The seed photons are assumed to be injected as an instantaneous flash of energy driven by a monochromatic photon source.
Following the work of Colpi (1988), the transport equation considered here includes the effects of both bulk and thermal Comptonization, as well as the spatial transport of radiation via diffusion and advection. However, we extend the model of Colpi in two important ways. First, we investigate time-dependent solutions using a Fourier transform technique. Second, we impose free-streaming boundary conditions at finite values for the inner and outer radii of the region.
In the case of 1H 0707-495, the steady-state X-ray continuum spectrum is likely produced by an accretion disk, which may extend to an outer radius ∼ 10 3 R g (Zoghbi et al. 2010). Close to the black hole, the disk may extend to the radius of marginal stability, or, alternatively, it may transition into a freely-falling, quasi-spherical corona, either due to passage through a shock (e.g., Chakrabarti 1989;Das et al. 2009;Becker et al. 2011;Chattopadhyay & Kumar 2016), or due to the accretion of a separate population of matter that is supplied with low angular momentum (e.g., Proga & Begelman 2003;Moscibrodzka et al. 2007), or possibly due the removal of a substantial fraction of the angular momentum via a mildly relativistic wind (e.g., Dauser et al. 2012;Done & Jin 2016). In this work, we assume that the disk either extends inward to the radius of marginal stability, with an overlying corona, or that the disk truncates near the outer radius of the corona due to one of the mechanisms listed above. Hence we assume that the time-dependent Comptonization is occurring in a quasi-spherical inner region, where the radial velocity, v r , and the electron number density, n e , follow the free-fall profiles v r ∝ r −1/2 and n e (r) ∝ r −3/2 , respectively.
We assume that the electrons in the quasi-spherical inner region comprise a Maxwellian distribution with temperature T e . This region is expected to be roughly isothermal due to the "thermostat" effect of thermal Comptonization (e.g., Sunyaev & Titarchuk 1980). The isothermal assumption is also supported by the simulations of Meyer-Hofmeister et al. (2012) who show that thermal Comptonization leads to a weak radial dependence for the electron temperature in the corona. A simplified schematic representation of our model is provided in Figure 1.

Transport Equation
Our radiation transport model expands upon the work of Colpi (1988) and Kroon & Becker (2016) by incorporating the effects of time-dependent bulk and thermal Comptonization in a spherically-symmetric converging flow. We begin by writing down the general time-dependent radiation transport equation describing the evolution of the photon distribution function, f ( r, , t), in an arbitrary geometry, which is given by (Becker 2003) where denotes the photon energy, r is the spatial location, v is the accretion velocity, σ T denotes the Thomson cross section, κ is the spatial diffusion coefficient, k denotes the Boltzmann constant, T e is the (constant) electron temperature, m e is the electron mass, c is the speed of light, and Q denotes the photon source term. The third and fourth terms on the right-hand side represent the effects of bulk and thermal Comptonization, respectively.
The radiation distribution function, f , is normalized so that the photon number and energy densities, n r and U r , respectively, are computed using and In a spherically-symmetric coronal region, the transport equation can be rewritten as where v r < 0 is the radial accretion velocity and f G (r, r 0 , , 0 , t) denotes the Green's function, which is the response to the instantaneous injection of N 0 photons of energy 0 at time t 0 from a source at radius r 0 .
The fundamental transport equation (Equation (10)) is linear, and therefore it follows that once the solution for the Green's function, f G , has been obtained, the corresponding particular solution, f , associated with the general photon source function, Q, is given by the integral convolution In our model, the bounds for the radial integration in Equation (11) are replaced with the inner and outer radii of the quasi-spherical corona.

Dynamics
In the quasi-spherical coronal region under consideration here, the radial velocity, v r , and the electron number density, n e , are assumed to follow the free-fall profiles v r ∝ r −1/2 and n e (r) ∝ r −3/2 , respectively. The radial velocity v r (r) profile can therefore be written as where v ff (r) denotes the local free-fall velocity, and the constant lies in the range 0 ≤ ≤ 1, with = 1 corresponding to exact free-fall (Colpi 1988).
Assuming that the accreting gas is composed of pure, fully-ionized hydrogen, the spatial diffusion coefficient, κ, appearing in Equation (10) is related to the electron number density, The electron number density, n e , is related to the mass accretion rate,Ṁ , viȧ where m p denotes the proton mass. Combining Equations (12), (13), and (14), we note that the radial velocity, v r , the electron number density, n e , and the spatial diffusion coefficient, κ, can be written as v r (r) = −vcr −1/2 , n e (r) =r where R g = GM/c 2 , and we have introduced the dimensionless constantsv andκ and the dimensionless radiusr We note that the constantv introduced in Equation (15) is related to Colpi's constant viav = √ 2 . Hence settingv = √ 2 yields the exact free-fall case with v r (r) = v ff (r) (see Equation (12)). We also note that Equations (14) and (15) can be combined to relatev and κ to the accretion rateṀ viaṀ Transforming the spatial coordinate from r tor in the transport equation and substi-tuting for v r , n e , and κ using Equation (15) yields where we have also introduced the dimensionless photon energy, χ, and the dimensionless temperature, Θ, using

Trapping Radius
In spherically-symmetric accretion flows, the spatial transport of the radiation is dominated by inward-bound advection for r < r t , and by outward-bound diffusion for r > r t , where r t is the "trapping radius," defined by (Payne & Blandford 1981;Colpi 1988) We can also use Equation (17) to express the trapping radius in the form where the dimensionless accretion rate,ṁ, is defined bẏ and the Eddington accretion rate,Ṁ E , is defined bẏ It is also of interest to compute the electron scattering optical depth, τ (r), between radius r and the outer radius of the spherical coronal region, located at r = r out . The result obtained is wherer = r/R g ,r out = r out /R g , and we have utilized Equation (15) to substitute for n e (r).
As noted by Payne & Blandford (1981) and Colpi (1988), the physical significance of the trapping radius motivates a transformation of the spatial variabler in the transport equation to the new spatial variable y, defined by so that the trapping radius corresponds to y = 1, and the radiation is trapped in the region with y > 1 (r < r t ). Making the change of variables in the transport equation fromr to y yields where the dimensionless time, q, is related to the time t via q ≡κ and q 0 , y 0 , and χ 0 denote the dimensionless injection time, injection radius, and injection energy, respectively, given by

Fourier Transform Method
In order to compute theoretical time lags using Equation (5), we must use our transport equation formalism to evaluate the Fourier transforms of the light curves observed at the soft and hard channel energies, s and h , respectively. The required solution for the Fourier transform can be obtained by performing a Fourier transformation of the partial differential transport equation, Equation (26). This yields an ordinary differential equation satisfied by the Green's function Fourier transform, F G , which is related to the Green's function, f G , via the fundamental definition (cf. Equation (3)) We can also express the Fourier transform in terms of the dimensionless time, q, by writing where the dimensionless Fourier frequency,ω, is related to ω viã ω ≡v Performing a Fourier transformation of Equation (26) by applying the operator +∞ −∞ e iωq dq yields an ordinary differential equation for the Fourier transform Green's function, F G . The result obtained is

Separation Functions
In order to solve Equation (32) to determine the solution for the Fourier transform Green's function, F G , we note that for χ = χ 0 , the source term vanishes, and consequently Equation (32) can be separated using the function where λ denotes the separation constant, and the functions G(y) and H(x) represent the spatial and energy separation functions, respectively. Substituting Equation (33) into Equa-tion (32) and rearranging terms yields, we find that for χ = χ 0 , The separation constant λ is independent of the coordinates y and χ, and therefore we can obtain two distinct ordinary differential equations satisfied by the functions G(y) and H(χ). The equations obtained are and 1 χ 2 where

Spatial Boundary Conditions
The spatial boundary conditions utilized in this model are based on the transition to free-streaming that must occur at the inner and outer surfaces of the spherical coronal region. We note that the free-streaming or "absorbing" boundary condition was also discussed by Titarchuk et al. (1997). The specific (spatial) radiation flux, F, also referred to as the "streaming function," is given by (Becker 1992) The inner and outer boundaries of the spherical coronal region are located at r = r in and at r = r out , respectively, and the corresponding values of the dimensionless location parameter y are (see Equation (25)) -15 - The free-streaming boundary condition operative at the inner and outer boundaries can be written as Transforming the free-streaming boundary conditions from r to y yields Fourier transformation of Equations (41) demonstrates that the spatial boundary conditions satisfied by the Green's function Fourier transform, F G , are given by By combining Equations (33) and (42), we conclude that the spatial separation function G(y) satisfies the free-streaming boundary conditions where primes denote differentiation with respect to y. The eigenfunctions, G n (y), represent the discrete set of solutions to Equation (35) that simultaneously satisfy both the inner and outer boundary conditions given by Equations (43) and (44), respectively. The corresponding eigenvalues for the separation constant λ are denoted by λ n . The solution procedure is further discussed below.

Spatial Eigenfunctions
The global solutions for the spatial eigenfunctions G n (y) are obtained using a threestep process. The first step is to numerically solve the ordinary differential equation, Equation (35), starting with the inner boundary condition (Equation (43)), imposed at y = y in . This yields the inner solution, denoted by G in (y). The second step is to carry out the numerical integration again, this time starting with the outer boundary condition (Equation (44)), imposed at y = y out , which yields the outer solution, G out (y).
For general, arbitrary values of the separation constant, λ, the inner and outer solutions are linearly independent functions. However, for certain special values of λ, the Wronskian of the inner and outer solutions vanishes, i.e., where and y * is located anywhere in the computational domain, so that y in ≤ y * ≤ y out . The eigenvalues λ n are the roots of Equation (45). We note that the Fourier frequency,ω, appears in Equation (35), and therefore it follows that a unique set of eigenvalues λ n is obtained for each value ofω.

Eigenfunction Orthogonality
The final closed-form solution for the Fourier transform F G can be expressed as an infinite series if we can demonstrate that the spatial eigenfunctions, G n (y), form an orthogonal set. We can establish the orthogonality of the spatial eigenfunctions as follows. First we rewrite Equation (35) in the Sturm-Liouville form d dy S(y) dG n dy + Q(y)G n + λ n 2 Ω(y)G n = 0 , where the weight function, Ω(y), is defined by and the functions S(y) and Q(y) are defined by S(y) = y −3/2 e −y , Next, we assume that λ n and λ m represent distinct eigenvalues associated with the spatial eigenfunctions G n and G m , respectively. Equation (47) is then multiplied by G m and the result is subtracted from the same equation with the indices n and m interchanged. After some algebra, the result obtained is Integration by parts between the inner and outer radii y in and y out yields Utilizing the spatial free-streaming boundary conditions at the inner and outer radii (Equations (43) and (44)), we find that the left-hand side of Equation (51) vanishes. Hence we obtain yout y in This establishes that the spatial eigenfunctions G n and G m are orthogonal with respect to the weight function Ω(y) in the computational domain located between the inner and outer free-streaming boundaries located at y = y in and y = y out , respectively.

Energy Eigenfunctions
The energy separation functions, H(χ), satisfy the second-order ordinary differential equation given by Equation (36). In addition, in order to obtain convergent results for the photon number and energy densities using Equations (8) and (9), we note that H must not increase faster than −3 as → 0. Likewise, H must decrease faster than −4 as → ∞. Additionally, H must be continuous at the injection energy χ = χ 0 in order to avoid an infinite diffusive flux in the energy space. Becker & Wolff (2007) obtained the exact solution to Equation (36) satisfying the required boundary and continuity conditions. The result is given by their Equation (48), which can be written as where M κ,µ , and W κ,µ denote the Whittaker functions, and we have introduced the parameters κ ≡ 1 2 (δ + 4) , µ ≡ 1 2 (3 − δ) 2 + 4 δλ n 1/2 .

Final Closed-Form Solution
Having established that the spatial eigenfunctions form a complete orthogonal set, it is now possible to express the Fourier transform Green's function, F G , using the infinite series where y 0 indicates the injection location, and N max denotes the maximum index included in the sum, which must be large enough to ensure convergence of the numerical results for the time lags. The computation of the time lags is further discussed in Section 5.
The expansion coefficients, C n , can be determined by exploiting the orthogonality of the spatial eigenfunctions and applying the derivative jump condition which is obtained by integrating Equation (32) over a small region in energy space around χ = χ 0 . Equation (56) is now substituted into Equation (57) where the Wronskian, W(χ 0 ), is defined by (Becker & Wolff 2007) and primes denote differentiation with respect to χ. The final result in Equation (59) is obtained using Equation (55) from Becker & Wolff (2007).
The solution for the expansion coefficients, C n , is obtained by multiplying Equation (58) by G m (y)Ω(y) and integrating with respect to y from y in to y out . Utilizing the orthogonality of the spatial eigenfunctions (Equation (52)), after some algebra, the final expression obtained for C n is C n = δN 0κ 3/2 G n (y 0 )Ω(y 0 )e χ 0 /2 e iωq 0 y 0 5/2 Γ(µ − κ + 1 2 ) 2πv 5/2 R 2 g c(m e c 2 ) 3 Θ 3 χ κ 0 I n Γ(1 + 2µ) , where the quadratic normalization integrals, I n , are defined by Combining Equations (56) and (60), we find that the final closed-form solution for the Fourier transform Green's function, F G , can be written as where χ min , χ max , κ, and µ are computed using Equations (54) and (55). Equation (62) gives the final solution for the Fourier transform, F G , of the radiation Green's function, f G , inside the corona, evaluated at dimensionless energy χ, dimensionless location y, and dimensionless Fourier frequencyω, which correspond to the dimensional quantities , r, and ω, respectively. The solution describes the response to the impulsive injection of N 0 photons with energy 0 at radius r 0 and time t 0 , which correspond to the dimensionless location y 0 , dimensionless energy χ 0 , and dimensionless time q 0 .
Based on the linearity of the fundamental transport equation (Equation (10)), it is straightforward to show that the Fourier transform associated with any desired photon source distribution, Q, can be evaluated using the integral convolution (cf. Equation (11)) In practice, the upper and lower bounds for the radial integration in Equation (63) are replaced with the inner and outer radii of the quasi-spherical corona. In Section 4, we use Equation (62) to compute the Fourier transform of the escaping radiation field, which is the basis for the time lag calculations.

Time Lag Computation
The time-dependent X-ray photon count-rate spectrum observed at the detector, F ( ), is simply related to the radiation distribution escaping from the outer edge of the spherical region, located at radius r = r out . Since we have employed a free-streaming boundary condition at the outer radius (Equation (40)), the computation of the escaping radiation spectrum is self-consistent, and we can therefore write the time-dependent X-ray photon count-rate spectrum observed at the detector as where D is the distance to the source. In order to utilize the formalism developed here to compute the observed time lags, we need to evaluate the Fourier transform of the observed count-rate spectrum, denoted byF ( , ω), given bỹ or, equivalently (cf. Equation (30)), Applying a Fourier transformation to both sides of Equation (64) where F G is evaluated using Equation (62).
The theoretical soft and hard energy light curves, s(t) and h(t), respectively, required to compute the time lags using Equation (5) are related to the photon count-rate spectrum, F ( , t), via Equations (4). Applying a Fourier transformation to Equations (4) yields the corresponding expressions where S and H denote the soft-and hard-energy light curve Fourier transforms, respectively, and the functionF is evaluated using Equation (67).
In order to compute X-ray time lags using our theoretical model, we must initially specify five parameters, namely: (1) the black hole mass, M ; (2) the dimensionless accretion rate,ṁ =Ṁ /Ṁ E ; (3) the dimensionless velocity parameter,v; (4) the dimensionless inner radius,r in = r in /R g ; and (5) the dimensionless outer radius,r out = r out /R g . We can then determine the value ofκ using Equation (21), which givesκ =v/(3ṁ). This allows us to generate a corresponding vector of values for the dimensionless Fourier frequency,ω, using Equation (31).
After determining the set of eigenvalues and spatial eigenfunctions for each selected value of the dimensionless Fourier frequency,ω, we must then vary the remaining theory parameters in order to obtain an acceptable qualitative fit to the observed time lags. The remaining four theory parameters comprise (1) the energy of the injected seed photons, 0 ; (2) the dimensionless electron temperature Θ = kT e /(m e c 2 ); (3) the hard photon channel energy, h ; and (4) the soft photon channel energy, s .

Application to 1H 0707-495
Fabian et al. (2009) and Zoghbi et al. (2010) have presented and discussed observations of time lags from the narrow-line Seyfert 1 galaxy 1H 0707-495. The time lags comprise an interesting patten, with the source displaying hard lags of magnitude ∼ 100 s for Fourier frequencies ν F 10 −3 Hz, and soft lags on the order of ∼ 10 s for Fourier frequencies ν F 10 −3 Hz. We use our model to replicate the hard and soft lags observed from 1H 0707-495. The simulated time lags, δt, can be computed using Equation (5) once the Fouriertransformed soft and hard channel energy light curves, S(ω) and H(ω), respectively, have been evaluated using Equation (68).

Model Parameters
The mass of the black hole at the center of 1H 0707-495 is estimated to be M = 2 × 10 6 M (Zhou & Wang 2005). Assuming that the steady-state emission from the source is isotropic, the X-ray luminosity, L X , is computed using where F X = 1.44 × 10 −10 erg cm −2 s −1 is the observed X-ray flux in the energy range 0.3-4 keV, calculated using data from Zoghbi et al. (2010) and D = 181 Mpc is the luminosity distance (Sani et al. 2010). The X-ray luminosity obtained using Equation (69) is therefore L X = 5.63 × 10 44 erg s −1 .
We will assume here that the inner edge of the spherical region, at radius r = r in , is located at the radius of the innermost stable prograde circular orbit, r ISCO , so that The value of r ISCO is computed using (Shapiro & Teukolsky 1983) where a = J/(M cR g ) denotes the dimensionless spin parameter for a black hole with mass M and angular momentum J. In the case of a non-rotating black hole, a → 0, and Equation (71) reduces to the Schwarzschild result, r ISCO = 6R g . The radius of the event horizon for prograde orbits, r H , is given by which reduces to r H = 2R g in the Schwarzschild limit.
We set the spin parameter for the black hole using a = 0.98 (Fabian et al. 2009), in which case the inner edge of the spherical region is located at radius r in = r ISCO = 1.61 R g (Equation (71)), and the event horizon is located at r H = 1.20 R g (Equation (72)). The outer edge of the spherical region is set at r out = 120 R g . Recall that in our model, the inner and outer radii are associated with the spherical corona, rather than the accretion disk (see Section 2). The value used for the dimensionless velocity parameter isv = √ 2, which corresponds to exact free-fall.
According to Equation (21), this yields for the dimensionless diffusion parameterκ = 0.429. The values of the inner and outer dimensionless location parameters, y in and y out , can be computed using Equation (39), which yield y in = 2.04 and y out = 0.028. We note that according to Equation (21), the trapping radius is located at r t = 3.3R g forṁ = 1.1. The primary model parameters are listed in Table 1. [Ht] The electron scattering optical depth measured inward from the outer edge of the spherical region, τ (r), is computed using Equation (24), and plotted in Figure 2. We note that the optical depth down to dimensionless radiusr = r/R g = 1.86 is close to unity, which is required in order to ensure that the relativistically blurred iron line emission near the black hole can be resolved spectroscopically, in reasonable agreement with the conclusions of (Wilkins et al. 2014). We find that τ = 1.08 at r = r ISCO , and decreases monotonically out to τ = 0 at r = 120 R g .

Simulated Time Lags
Since the eigenvalues, λ n , in our calculation are functions of the Fourier frequency, ν F , we must select a sample of discrete Fourier frequencies and evaluate the time lags at those frequencies. The reader should keep in mind that a continuous range of Fourier frequencies is not required in our model, since we do not need to perform the inverse Fourier transformation. Instead, the required time lags are computed directly from knowledge of the Fourier transform itself.
The first 10 complex eigenvalues, λ n , obtained for each discrete value of the Fourier frequencyω adopted here are plotted in Figure 3, with the value ofω increasing from bottom to top. The eigenvalues are indicated by the points, and the colored lines are added to clarify the eigenvalue sequence for each frequency. As noted in Section 4, the computation of the simulated time lags, δt, using our theoretical model also requires specification of the injection radius r 0 , the photon injection energy 0 , the electron temperature Θ = kT e /(m e c 2 ), and the hard and soft channel energies, h and s , respectively. We set the value of the source energy using 0 = 0.89 keV, which notably correlates with the broadened and skewed energy of the iron L-line emission that Fabian et al. (2009) reported for 1H 0707-495. We set the value of the injection radius using r 0 = 16R g , which yields y 0 = 0.206 according to Equation (28). The value adopted for the dimensionless electron temperature is Θ = 0.05, so that kT e = 25.6 keV. This yields for the electron temperature T e = 2.96 × 10 8 K, which is consistent with the expected temperature range for the corona (e.g., Wilkins et al. 2014;Kara et al. 2017).
A remaining issue is how we set the soft and hard channel energies in our model, s and h , respectively. In their time lag computations, Fabian et al. (2009) use s = 0.3 − 1 keV for the soft energy band and h = 1 − 4 keV for the hard energy band. However, our model is not able to accommodate an energy range for the soft and hard bands, and instead we need to specify exact values for s and h . We set our hard and soft channel energies to the values s = 0.33 keV and h = 1.76 keV in order to approximate the variability found in these energy bands.
We have analyzed the convergence of the results obtained for the time lags as a function of the value of the maximum index, N max , appearing in the sum in Equation (62). We find that setting N max = 5 is sufficient to ensure convergence to a relative error of ∼ 1%, and therefore we generally adopt the value N max = 10 as a conservative value when computing the time lags for each Fourier frequency ω.
The resulting time lag distribution, δt, computed using Equation (5), is plotted as a function of the Fourier frequency ν F = ω/(2π) in Figure 4, and compared with the observational data for 1H 0707-495 reported by Fabian et al. (2009) andZoghbi et al. (2010). Note the good qualitative agreement between the theoretical model and the time-lag data from Fabian et al. (2009) andZoghbi et al. (2010), including the change in sign of δt from positive (hard lag) below Fourier frequency ∼ 5 × 10 −4 Hz to negative (soft lag) at higher frequencies.

DISCUSSION AND CONCLUSION
The short time lags of ∼ 10 − 100 seconds detected from many AGNs imply that the physical processes are occurring in the inner region of the accretion flow onto the black hole. Hence, the analysis and interpretation of the time lags provides a unique opportunity to probe the detailed structure of the innermost portion of the accretion flow, including the temperature and density profile, as well as the accretion morphology.
This has motivated us to develop a new model for the production of the observed hard and soft time lags from AGNs based on a combination of thermal and bulk Comptonization occurring in the quasi-spherical inner region of the accretion flow onto the black hole. In our model, the time-dependent emission that is responsible for generating the time lags is produced in the quasi-spherical region, and the steady-state emission, comprising the bulk of the X-ray signal, is produced in the surrounding, geometrically thin but optically thick Time Lag δt (s) accretion disk (Shakura & Sunyaev 1973).
We have applied our new model to the interpretation of the X-ray time lags observed from 1H 0707-495, as reported by Fabian et al. (2009) andZoghbi et al. (2010). The theoretical results for the time lags plotted in Figure 4 indicate good qualitative agreement with the time lag data for 1H 0707-495, suggesting that the observed soft lags may be produced via thermal and bulk Comptonization occurring within the quasi-spherical region of the accretion flow. Hence in our model, the hard and soft time lags are both produced in a single region, via the action of a unified physical mechanism.
Our physical model is based on a rigorous theoretical framework describing the reprocessing of an instantaneous flash of monochromatic photons injected at radius r 0 = 16 R g . The value for the injection energy adopted here, 0 = 0.89 keV, is approximately at the peak value for the energy of the broad iron L-line reported by Fabian et al. (2009), and hence may be due to fluorescence of the iron line, perhaps driven by inhomogeneities or "clumps" in the accretion flow (e.g., Merloni et al. 2006;Gutiérrez et al. 2021).
In the reverberation interpretation, the broadened iron L-line feature observed in the steady-state X-ray spectrum from 1H 0707-495 is a spatial reflection feature (e.g., Zoghbi et al. 2010). This is similar to our model, except that in our model, the seed photons may correlate with broadened iron L-line emission generated inside the corona. Hence the location of the source for the L-line emission in the steady-state spectrum is not necessarily identical with the source location for the impulsive injection of the seed photons used in our calculation of the time lags, although the physical emission mechanism may be the same.
In our application to 1H 0707-495, we have approximated the broadened iron L-line emission using a monochromatic source. This approximation can be generalized in future work using an integral convolution of the Green's function (Equation (63)) with the broadened source distribution. However, we do not expect this modification to substantially alter the results presented here.
The size of the corona assumed here bears further discussion. Morgan et al. (2012) published microlensing results that constrain the X-ray corona of lensed quasars to be smaller than 10 R g , which would seem to be a potential conflict with the size of the corona in our model, which is r out ∼ 120 R g . This problem can be resolved by noting that the overall scattering optical depth of the corona in our model is τ ∼ 1.08, and therefore the distribution of the last scattering over space indicates that many photons will escape promptly from deep layers without much scattering. Hence it follows that the apparent size of the X-ray emitting region in our model is not necessarily equal to the size of the corona implied by the microlensing observations. We assume that the flash of seed photons is produced simultaneously at all points on a spherical shell with radius r 0 . Of course, relativity precludes any coherence in the emission process occurring on timescales shorter than the light-crossing time for the emitting region. In this regard, it's important to note that one would obtain the same result for the Fourier transform solution that we obtain here even if the iron L-line photons where emitted at random locations on the spherical surface. This is because the time window for the Fourier transform "sweeps up" many such emission episodes. Hence, provided they are distributed in a spherically symmetric way, the total Fourier transform would remain in agreement with that computed here, and there would be no causality violation.
The simplicity of our model is in marked contrast to previous models, which generally rely on unusual geometrical configurations, or multiple thermal zones with different properties, in order to reproduce the observed time lags. For instance, Zoghbi et al. (2010) and Miller et al. (2010) invoke different physical mechanisms to explain the production of the hard time lags at small Fourier frequencies and the soft time lags at large Fourier frequencies.
Conversely, in our model, both the hard and soft lags are produced as a natural consequence of thermal and bulk Comptonization occurring in the quasi-spherical inner region of the accretion flow.
It is interesting to further explore the behavior of the time lag δt as a function of the Fourier frequency ν F = ω/(2π). In particular, we note the change in sign of the time lags, at what we term the "critical frequency," ν c ∼ 5 × 10 −4 Hz (see Figure 4). The complex behavior of δt as a function of Fourier frequency, and the sign change between hard and soft lags at the critical frequency, stem from the fact that the energy of the injected iron L-line emission, 0 = 0.89 keV, lies between the soft and hard detector band energy ranges used to compute the time lags. The existence of hard lags for ν ν c reflects the fact that on long time scales, recoil losses are more rapid than stochastic energization. Hence the drift towards high energies is slower than the drift towards lower energies, and this naturally produces hard time lags. Conversely, the appearance of soft lags for ν ν c reflect the rapid energization of photons that occurs in a single collision with an electron, which overwhelms recoil losses on short time scales (Lewis et al. 2016).
De Marco et al. (2013) noted that the Fourier frequency associated with the onset of negative lags (i.e., soft lags), denoted by ν neg , scales as ν neg ∝ (Ṁ /Ṁ E )/(M ). For sources accreting near the Eddington limit, withṀ ∼Ṁ E ∝ M , the negative lag frequency scales as ν neg ∝ M −1 . Within the context of our model for bulk and thermal Comptonization, the negative lag frequency of De Marco et al. (2013) is roughly equivalent to our "critical frequency," ν c , at which the time lags transition from soft to hard. We therefore find that the universality of the mass scaling law discovered by De Marco et al. (2013) can then be understood as a natural consequence of our model, since, according to Equation (31), where ω c = 2πν c denotes the critical Fourier frequency, and the dimensionless critical frequencyω c = 0.131 is a constant in our model. Hence Equation (73) implies that the critical frequency ω c is inversely proportional to the mass of the black hole, in agreement with the results of De Marco et al. (2013). This strongly suggests that thermal and bulk Comptonization play an important role in spectral formation in all AGNs displaying soft time lags, regardless of mass scale.
In this paper we have presented a new model for the production of both hard and soft time lags in AGN via the thermal and bulk Comptonization of seed photons, which may be injected via fluorescence of the iron L-line. We developed a specific detailed calculation applicable to 1H 0707-495 and demonstrated good agreement with the time lag observa-tions for this source. Although our model does not include an explicit reflection component, nonetheless we believe it is complementary to the previous models because it includes reverberation in the energy space via bulk and thermal Comptonization, as well as reverberation in the physical space, due to photon diffusion through the corona. In this sense, the model explored here explains all of the time lags observations for 1H 0707-495 using a single unified physical mechanism.
We note that in other sources there are soft lag features observed at higher energies than those observed from 1H 0707-495, which may be related to iron K-line emission. In future work we plan to extend our model to treat additional photon sources covering a broader energy range in order to more fully explore the role of thermal and bulk Comptonization in producing both hard and soft time lags from AGNs.
Although the effects of general relativity have not been fully incorporated here, we expect that the same fundamental results would be obtained in an equivalent, fully relativistic model, because the fundamental physical process leading to the production of the time lags is electron scattering in the relatively cool, spherical region of the accretion flow, which is rigorously treated in our model.
We are grateful to the anonymous referee for providing numerous suggestions for clarification that led to substantial improvements in the manuscript.