Simulating X-Ray Reverberation in the Ultraviolet-emitting Regions of Active Galactic Nuclei Accretion Disks with Three-dimensional Multifrequency Radiation Magnetohydrodynamic Simulations

Although active galactic nuclei (AGN) are among the most luminous objects in the Universe, due to their great distance, it is generally not possible to resolve the microarcsecond scales of their accretion disks, broad-line regions (BLRs), or event horizons. However, because AGN light curves vary over a wide range of frequencies, temporal resolution can be used in place of spatial resolution. A primary technique used to study AGN, called reverberation mapping, makes use of this temporal resolution (Peterson 2014).

Reverberation mapping was first proposed by Blandford & McKee (1982) as a way to measure distances to the BLR by finding a time lag between variability in continuum light curves and variability in broad emission lines (see also Kaspi et al. 2000; Peterson et al. 2004; Bentz & Katz 2015; Grier et al. 2017). Disk continuum reverberation mapping can similarly be used to measure the radial extent of the accretion disk itself because X-rays emitted by the corona move outwards along the disk and are absorbed and reemitted by different temperature regions of the disk at different wavelengths based on the local temperature (e.g., Sergeev et al. 2005; Cackett et al. 2007; De Rosa et al. 2015; Edelson et al. 2015, 2017, 2019; Fausnaugh et al. 2016; Jiang et al. 2017; Starkey et al. 2017; Cackett et al. 2018, 2021; Homayouni et al. 2022). Therefore, there will be a lag on the light-crossing timescale between variability in light curves from short wavelength bands and light curves that have been reprocessed and reemitted in longer wavelength bands. This lag is used to determine the distance between different temperature regions of the accretion disk.

Traditional models for reverberation mapping make the assumption that the variability in UV and optical AGN light curves is driven by the variability of the X-ray irradiation. However, there is significant evidence for additional variability in UV-optical light curves that does not appear to come from X-ray light curves (e.g., Uttley et al. 2003; Arévalo et al. 2009; Edelson et al. 2019; Yu et al. 2022; Cackett et al. 2023; Hagen & Done 2023; Kara et al. 2023). The lack of correlation between X-ray and UV-optical light curves could be due to absorption by intervening material (e.g., Kara et al. 2021) or changes in the height and temperature of the corona on timescales shorter than the observing timescale (Panagiotou et al. 2022a). Another alternative is that this variability could come from fluctuations in the UV-optical regions of the accretion disk due to magnetorotational instability (MRI)–driven turbulence (Balbus & Hawley 1991) or convection resulting from the enhanced opacity in the UV-optical region of the disk (Jiang & Blaes 2020).

While there has been considerable observational work done using reverberation mapping, theoretical work to model the underlying physics of the reprocessing of hard irradiation into the softer radiation emitted by the disk has lagged behind. Most of the previous work modeling this reprocessing is (semi)analytic and makes assumptions about the vertical profiles of various disk parameters, which often remain fixed over time or are not calculated self-consistently to include the interaction of radiation and gas (e.g., Sun et al. 2020; Kammoun et al. 2021; Panagiotou et al. 2022b; Salvesen 2022). The recent development of accurate, high-speed, numerical methods to calculate radiative transfer coupled to magnetohydrodynamics (MHD) makes it possible for the first time to examine how both radiation and magnetic fields have an impact on AGN disk structure and turbulence (Jiang et al. 2013; Jiang & Blaes 2020; Jiang 2021). Adding realistic Rosseland and Planck mean opacities that include line opacities instead of only electron scattering and free–free absorption also has an impact on driving convection in AGN disk models (Jiang & Blaes 2020). However, all of these MHD simulations model radiation using a single frequency group, with the radiation integrated over all frequencies, and do not examine the impact of X-ray irradiation on the disk and the light curves emitted by the disk.

In this paper we provide an early application of a multifrequency radiation MHD code, which instead of integrating radiation over all frequencies, integrates radiation over multiple frequency groups (Jiang 2022). This multifrequency capability allows us to simulate from first principles the reprocessing of high-frequency light incident on an accretion disk into light emitted in the UV based on the local disk temperature to better understand the drivers of variability in UV-optical AGN light curves. By first principles, we mean that the radiation energy and momentum are coupled to the gas evolution and that the temperature, opacity, and energy of the simulation depend self-consistently on the local gas density and pressure. In Section 2 we describe the methods we use for our simulations. In Section 3 we present, compare, and discuss the simulated disks and light curves from our simulations with and without X-ray irradiation. Finally, in Section 4 we summarize our results and present ideas for future work.

Using athena++ we solve the ideal MHD equations coupled with time-dependent, implicit radiative transfer equations for intensities over discrete angles in the manner described in Jiang (2022), which builds on Jiang (2021) by allowing for multiple frequency groups. We use two frequency groups, one group covering the frequency range [0, 1 keV] to model the radiation field emitted by the disk locally (at UV wavelengths) and the other group covering the frequency range (1 keV, ∞ ) for the X-ray irradiation. Notice that we turn off the frequency shift due to Doppler effect as its effect is small and can cause numerical diffusion in the frequency space for the two groups we adopt.

We model a local patch of the accretion disk around a supermassive black hole (SMBH) with mass M_BH = 5 × 10⁷ M_⊙ under the local shearing box approximation with the box centered at r₀ = 45R_s, where R_s is the Schwarzschild radius. The corresponding Keplerian orbital frequency Ω = 4.73 × 10⁻⁶ s⁻¹ and shear parameter $q=-d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}r=3/2$. In the shearing box approximation, the x- and y-coordinates correspond to the radial and azimuthal dimensions, while the z-coordinate corresponds to the height. We choose our simulation domain to be (x, y, z) ∈ (−1.5H_g , 1.5H_g ) × (−6H_g , 6H_g ) × (−48H_g , 48H_g ) with resolution N_x = 48, N_y = 192 and N_z = 1536. Here H_g = 1.11 × 10¹² cm is the gas pressure scale height of the disk, while the total pressure scale height (including radiation pressure) will be H_r = 4.06 H_g = 4.51 × 10¹² cm. We use the periodic and shear periodic athena++ boundary conditions for the x- and y-boundaries, respectively. For the z-boundaries, we use outflow boundary conditions for the gas. Intensities for the low-frequency group are allowed to leave the simulation domain, but they are not allowed to enter. For the high-frequency group, we set the incoming intensities based on a prescribed X-ray light curve. Outgoing intensities are copied from the last active zone to the ghost zone.

We first perform a one-frequency run to model only the locally emitted UV photons without X-ray irradiation. Then we restart the simulation with two frequency groups to model the locally generated UV photons and X-ray irradiation simultaneously. We use N = 48 discrete angles in each cell for each frequency group to resolve the angular distribution of the radiation field. Inside the simulation box, we initialize the UV radiation isotropically based on the local disk temperature and the X-ray radiation isotropically based on a constant radiation energy, ${E}_{{\rm{r}},1}={10}^{-4}{E}_{{\rm{r}},0}^{\mathrm{mid}}$, where ${E}_{{\rm{r}},0}^{\mathrm{mid}}$ is the midplane radiation energy for the UV frequency group.

We inject two different X-ray light curves for two otherwise initially identical two-frequency runs (Run A and Run B). The time dependence of the intensity of this radiation is determined by a mock X-ray AGN light curve that we model as a damped random walk (DRW) as in Kelly et al. (2009) because the DRW model has been shown to be a decent representation of AGN light curves for timescales of days to years (e.g., Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2012; Zu et al. 2013). For Run A, we use a damping timescale τ_damp = 30 days and a time sampling of dt = 1 day; for Run B, we use τ_damp = 0.1 day and a time sampling of dt = 0.001 day. For Run A, we choose a shorter damping timescale than the thermal timescale of our simulation box and a low time sampling to avoid introducing very high frequency variability. For Run B, we choose a much shorter damping timescale that is closer to the dynamical timescale of the innermost disk where the X-rays originate and a high time-sampling rate to capture this short timescale variability. In summary, Run A will have more variability at low frequencies than Run B, while at high frequencies, Run B will have more variability. The goal is to test and compare two very different injected light curves.

We found that if we normalized the intensity of the injected X-rays to the time-averaged total UV flux, leaving the top and bottom of the simulation box in the one-frequency run, the X-rays do not have a large impact on our simulation, at least over a few thermal timescales. Since we are interested in X-ray reprocessing, we instead normalize the X-ray irradiation to 8 times the average UV emission, which is the value we found sufficient to produce changes in the UV light curve. In the future, a parameter study is needed to determine how this normalization impacts the reprocessing. Finally, we average the injected flux over all angles in our simulation and inject this averaged intensity along all angles pointing into the box. Note that because X-rays are injected directly into the simulation, there will be no lag between the X-ray and UV light curves from the light-travel time.

We determine the timestep-dependent opacity for the UV component in each cell with a bilinear interpolation of the OPAL opacity tables from Iglesias & Rogers (1996) based on each cell's current gas temperature and density. For our X-ray component, we take the proper average of the free–free absorption opacity and electron scattering opacity, which ends up being dominated by the scattering value.

We set the initial vertical profile of our simulations using the equations of radiation hydrostatic equilibrium, as in Jiang et al. (2016). Initial and normalization parameters are given in Table 1. For our initial conditions, the ratio of the radiation to the gas pressure is P_rat = 4.4. We initialize the magnetic field in our simulation as in Jiang et al. (2013) as two oppositely twisted flux tubes with zero net flux and set the ratio of the gas pressure to the magnetic pressure to β = 200 in the midplane. We also add a purely vertical magnetic field with β = 2 × 10⁴ to increase the Maxwell stress in order to achieve a temporally averaged $\dot{M}/{\dot{M}}_{\mathrm{edd}}\approx 2\times {10}^{-2}$ once the simulation reaches steady state.

3.1. Disk Structure

We show the horizontally averaged vertical profiles as a function of time for the gas temperature, gas density, and radiation energy for the one-frequency run without X-ray irradiation (to the left of the gray line) and two-frequency Run B with X-ray irradiation (to the right of the gray line) in the right panels of Figure 1. In the top panel, the red arrows represent the time-varying X-ray irradiation we inject into the simulation from the upper and lower boundaries. The solid white lines show the location of the UV photosphere, which is very similar to the location of the X-ray photosphere, because both optical depths are dominated by the scattering opacity near the surface of the disk (z(τ_UV = 1) ≈ z(τ_XR = 1)). The dashed white line shows the absorption photosphere for the UV component z(τ_a,UV = 1). τ_a,XR < 1 throughout the whole disk, but the effective absorption optical depth is τ_a,XR ≈ 1.3. The left panels show the horizontally and time-averaged vertical profiles for the gas temperature, density, and turbulent velocity for the two-frequency Run B. The vertical profiles for two-frequency Run A, not shown here, are qualitatively the same as in Run B.

Zoom In Zoom Out Reset image size

The vertical profiles are very similar between the one- and two-frequency runs, except for the profile for the gas temperature (top-right panel), which is higher throughout the simulation when we inject X-rays. Unsurprisingly, the temperature increases in the photosphere the most and also within regions that are an additional 10 H_g, where the temperature can be as high as ∼ 10⁵ K. In contrast to isothermal simulations, we do no see any sign of disk winds in our simulations (e.g., Bai & Stone 2013).

In the bottom panel of Figure 1, we show the radiation energy for the low- and high-frequency group above and below the midplane, respectively, to the right of the gray line. The midplane radiation energy is dominated by the energy from the low-frequency group. However, at around ±30 H_g, the radiation energy from the high-frequency group becomes the dominant radiation energy. On average, the energy of the X-ray radiation manages to penetrate down to τ_XR ≈ 240 or ±6.5 H_g before it drops by an e-folding.

Overall, the α-viscosity parameter in our simulations ranges from around 0.025 to 0.04. In all runs, the Maxwell stress is about a factor of 6 larger than the Reynolds stress. The average thermal timescale is τ_therm = 1/(αΩ) ≈ 80 days in all runs.

3.2. Intrinsic Variability in UV Light Curves

We show the intrinsic UV light curve from the one-frequency run without X-ray irradiation and the power spectral density (PSD) of this light curve in blue in the top panels of Figures 2 and 3, respectively. The intrinsic UV light curve has variability over timescales of monthly to subdaily that is driven by MRI turbulence and convection. There is evidence that the UV PSD is damped at the lowest frequencies (f ≲ 1/(30 days)) and then goes as f⁻² for 0.03 < f < 0.2 day⁻¹, suggesting that a DRW model might be a good fit for the intrinsic UV light curve at these frequencies. Recent examinations of long-baseline optical AGN light curves have suggested that the DRW damping timescale is related to the thermal timescale of the disk (Kelly et al. 2009; Burke et al. 2021; Stone et al. 2022). The damping timescale of around 30 days that we see evidence of in Figure 3 is very similar to the thermal timescale near the surface of the disk in our simulations (∼30 days), but longer light curves are needed to robustly model this damping timescale.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the DRW model may be a good fit at lower frequencies, the slope of the UV light curve's PSD is steeper, ∼f^−2.5, than the DRW model at f ≳ 1/(5 days). In high-cadence observations of AGN light curves from Kepler, on timescales around 10 days, there is evidence that the PSD slope is steeper than the f⁻² power law of the DRW (e.g., Smith et al. 2018). In our simulations, the break may be at a higher frequency because we are simulating the UV region as opposed to the optical region being probed by Kepler. Tachibana et al. (2020) proposed that this break comes from the fact that variability will be averaged over different regions of the disk. Our simulations are highly localized and still have a steeper slope at high frequencies solely because disk fluctuations are unable to produce much variability on timescales this short. However, the lack of averaging in our local simulations may be another reason why we do not see a break at frequencies as low as in Kepler light curves.

Nonetheless, we find that our intrinsic UV light curve has variability on timescales of monthly to subdaily, due to MRI- and convection-driven fluctuations. Previous global MHD simulations suggest that the fluctuations producing this variability would be accreted inwards (Hogg & Reynolds 2016), leading to inward-propagating lags (i.e., shorter wavelengths lagging longer wavelengths) on the inflow timescale, otherwise known as long negative lags. Evidence for these long negative lags have been observed for the AGN Fairall 9 (Hernández Santisteban et al. 2020; Yao et al. 2023; Secunda et al. 2023).

3.3. Hard X-Rays Have Little Impact on UV Emission

Using the multifrequency radiation MHD code from Jiang (2022) allows us to simulate the reprocessing of X-ray irradiation into UV emission from first principles. This reprocessing is a key component of reverberation mapping, one of the primary tools used to study AGN disks. The traditional model for reverberation mapping assumes that X-rays are the main driver of variability in UV and optical light curves emitted by the AGN disk. However, the results from our simulations do not agree with this assumption. Instead we find that there is significant intrinsic UV variability generated locally by the MRI and convection. This variability can be described by a DRW, with τ_damp ≈ 30 days on timescales greater than 5 days, but shows less power than the DRW at higher frequency, in agreement with AGN light curves observed at high cadence by Kepler (Smith et al. 2018).

The X-rays we inject into our simulations have almost no impact on the PSDs of the reprocessed UV light curves in our two-frequency runs because the short timescale variability of the X-ray irradiation gets smoothed when it is scattered in the outer regions of the disk. As a result, the X-ray and UV light curves from our two-frequency runs are not well correlated if there is not sufficient X-ray variability on timescales greater than a day. Even with longer timescale variability, we only find a correlation of r = 0.34 between the X-ray and UV light curves, again because of the low absorption opacity for X-rays in our simulations. This correlation agrees with observational results that suggest that X-ray and UV light curves are not as well correlated as we might expect if X-rays are the main driver of UV and optical variability as in traditional reprocessing models (e.g., Schimoia et al. 2015; Edelson et al. 2019; Cackett et al. 2023). However, canonical reprocessing models could produce lower X-ray to UV correlations if the X-ray light curve is nonstationary, for example due to changes in the height or temperature of the corona on timescales shorter than the observing campaign (Panagiotou et al. 2022a).

There are many radiation MHD simulations of AGN disks, but most are one-frequency simulations, and none focus on the impact of X-ray irradiation on the UV-optical region of AGN disks and the light curves emitted by the disk. On the other hand, there are numerous sophisticated (semi)analytical calculations for disk reprocessing of X-ray irradiation (e.g., Sun et al. 2020; Kammoun et al. 2021; Panagiotou et al. 2022b; Salvesen 2022). However, they all make assumptions about the vertical profiles of AGN disk parameters and do not evolve the gas turbulence in three dimensions nor the impact that radiation has on the gas and vice versa as a function of time. On the other hand, the multifrequency radiation MHD simulations presented here provide detailed information on the nature of light-curve variability produced by disk turbulence and how X-ray light curves with different PSDs transfer energy and momentum to the AGN disk, affecting the variability in UV light curves emitted by the disk.

Increasing our understanding of the physical processes that influence reverberation mapping is crucial because there is growing evidence that the standard Shakura & Sunyaev (1973) thin disk model does not accurately model an AGN disk (e.g., Antonucci 1988, 2023; Antonucci et al. 1989). For example, radiation MHD simulations show that pressure support from magnetic fields and radiation can lead to puffed-up disks (e.g., Gaburov et al. 2012; Jiang et al. 2016, 2019; Jiang & Blaes 2020; Hopkins et al. 2024), while recent observations provide additional evidence for thicker disks than predicted by the standard thin disk model (Yao et al. 2023; Secunda et al. 2023). In addition, both microlensing and reverberation mapping campaigns suggest that the radial extent of the disk is larger than predicted by the standard disk model (e.g., Morgan et al. 2010; Jiang et al. 2017; Guo et al. 2022).

Improved intuition about the physical processes involved in reverberation mapping is essential to improving our models for AGN disks. Future work should look at how these results will depend on different amounts of absorption opacity for X-ray photons. Future work should also examine a wider range of parameters for AGN disks in order to determine how the physics of reverberation mapping changes as a function of parameters such as SMBH mass and accretion rate. Finally, expanding these shearing box simulations to global simulations is necessary to understand how the fluctuations producing this intrinsic variability are propagated throughout the disk, affecting the reverberation mapping signal. These simulations will allow us to make our reverberation mapping techniques more accurate and mine more information out of AGN light curves about the structure and internal physics of AGN disks.

The authors would like to thank Dr. Erin Kara for insightful comments. A.S. is supported by a fellowship from the NSF Graduate Research Fellowship Program under grant No. DGE-1656466. J.E.G. is supported in part by NSF grants AST1007052, AST1007094, and 22624. The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation. This work is also part of the TCAN collaboration supported by the grant 80NSSC21K0496.