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

Active galactic nuclei (AGN) light curves observed with different wave bands show that the variability in longer wavelength bands lags the variability in shorter wavelength bands. Measuring these lags, or reverberation mapping, is used to measure the radial temperature profile and extent of AGN disks, typically with a reprocessing model that assumes X-rays are the main driver of the variability in other wavelength bands. To demonstrate how this reprocessing works with realistic accretion disk structures, we use 3D local shearing box multifrequency radiation magnetohydrodynamic simulations to model the UV-emitting region of an AGN disk, which is unstable to the magnetorotational instability and convection. At the same time, we inject hard X-rays (>1 keV) into the simulation box to study the effects of X-ray irradiation on the local properties of the turbulence and the resulting variability of the emitted UV light curve. We find that disk turbulence is sufficient to drive intrinsic variability in emitted UV light curves and that a damped random walk model is a good fit to this UV light curve for timescales >5 days. Meanwhile, X-ray irradiation has negligible impact on the power spectrum of the emitted UV light curve. Furthermore, the injected X-ray and emitted UV light curves are only correlated if there is X-ray variability on timescales >1 day, in which case we find a correlation coefficient r = 0.34. These results suggest that if the opacity for hard X-rays is scattering dominated as in the standard disk model, hard X-rays are not the main driver of reverberation signals.


Introduction
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. 2015Edelson et al. , 2017Edelson et al. , 2019;;Fausnaugh et al. 2016;Jiang et al. 2017;Starkey et al. 2017;Cackett et al. 2018Cackett et al. , 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.

Methods
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 7 M e under the local shearing box approximation with the box centered at r 0 = 45R s , where R s is the Schwarzschild radius.The corresponding Keplerian orbital frequency Ω = 4.73 × 10 −6 s −1 and shear parameter = - W = q d d r ln ln 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 12 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 12 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 highfrequency 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 E 10 r,1 4 r,0 mid , where E r,0 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 timesampling 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 4 to increase the Maxwell stress in order to achieve a temporally averaged » ´-  M M 2 10 edd 2 once the simulation reaches steady state.

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 twofrequency Run A, not shown here, are qualitatively the same as in Run B.
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 5 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.

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 −2 for 0.03 < f < 0.2 day −1 , 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.
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 highcadence observations of AGN light curves from Kepler, on timescales around 10 days, there is evidence that the PSD slope is steeper than the f −2 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 MRIand convection-driven fluctuations.Previous global MHD simulations suggest that the fluctuations producing this r 0 = 45 R s = 6.68 × 10 14 cm τ 0 6.26 × 10 4 P g,0 2.77 × 10 6 dyn cm −2 P rat = P r,0 /P g,0 4.37 c rat = c/c s 5.69 × 10 3 τ therm 80 days Notes.H r is the initial total (gas and radiation combined) pressure scale height.r 0 is the fixed location of the center of the shearing box, and τ 0 and P g,0 are the initial midplane optical depth and gas pressure, respectively.P rat and c rat are the initial ratio of the radiation pressure to the gas pressure and the ratio of the speed of light to the sound speed, respectively, at the midplane.We also give the approximate thermal timescale (τ therm ) calculated from the average α viscosity in our simulations once they reach steady state.We show the UV radiation energy above the midplane and the X-ray radiation energy below the midplane.In all panels, the solid white line is z(τ UV = 1) ≈ z(τ XR = 1), and the dashed white line is z(τ a,UV = 1).τ a,XR < 1 throughout our simulation.Left panels: horizontally and time-averaged vertical profiles of gas temperature, density, and turbulent velocity in the top, middle, and bottom panels, respectively, for our two-frequency simulation.All figures are for Run A.

Hard X-Rays Have Little Impact on UV Emission
We compare the injected X-ray and reprocessed UV light curves from the two-frequency runs that include X-ray irradiation in the bottom panels of Figure 2 and the PSDs of these light curves in the top panels of Figure 3.Note that because we are injecting X-rays directly into the simulation, there will be no lag between the X-ray and UV light curves from the light-travel time.For Run A, the Pearson r coefficient between the injected X-ray light curve with damping timescale τ damp = 30 days and the reprocessed UV light curve is r = 0.34, with a p-value, p < 10 −5 .For Run B, where the damping timescale of the injected X-ray light curve is τ damp = 0.1 day, we do not get a significant correlation between the X-ray and reprocessed UV light curves, p = 0.64.This lack of correlation is despite the fact that the injected X-ray flux is a factor of 8 larger than the emitted UV flux, which is larger than typically observed.The lack of correlation between the X-ray and reprocessed UV light curves in Run B may be due to the short damping timescale of the X-ray light curve, which makes the X-ray light curve essentially white noise for all timescales of interest and therefore hard to cross-correlate.
For both runs, the injected X-ray light curves appear to have no impact on the PSDs of the reprocessed UV light curves, which are nearly identical to the intrinsic UV light curves from the one-frequency run without X-ray irradiation.This lack of impact is particularly surprising for Run B, where at high frequencies, the X-ray irradiation has around 6 orders of magnitude more power than the intrinsic UV light curve.The lack of additional high-frequency variability in the reprocessed UV light curve in Run B is apparent in the light curve itself as well.
The reason the UV PSD from Run B is not enhanced at high frequencies may be related to the fact that the emitted X-ray light curve in Run B has less power at f  1/(1 day) than the injected X-ray light curve (see bottom right panel of Figure 3).The smoothing of high-frequency variability in the injected X-ray light curve occurs because the X-ray irradiation gets scattered instead of absorbed by the disk.As a result, the injected and emitted light curves in Run B are less strongly correlated, r = 0.89, than in Run A, where r = 0.998 and there is less high-frequency variability to be smoothed.
Due to the very low absorption opacity for the X-ray irradiation, there is very little absorption and direct heating, as is typically assumed for reverberation mapping.Instead, our results suggest that the scattering of the incident X-ray photons off of electrons leads to a momentum transfer that eventually will thermalize, producing the small amount of disk heating we see in our simulations (see Figure 1) and the correlation we see between the X-ray and UV light curves in Run A. This momentum transfer and disk heating do have some impact on the reemitted UV light curve in Run B as well, which, in the top-right panel of Figure 2, is clearly not identical to the intrinsic UV light curve.For both two-frequency runs, the reprocessed UV light curves are not strongly correlated with the intrinsic UV light curve from the one-frequency run without X-ray irradiation, with r = −0.28(p < 10 −5 ) and r = 0.17 (p = 0.0032) for Run A and Run B, respectively.The anticorrelation between the UV light curves in Run A results from modeling the injected X-ray light curve in this run as a DRW with a damping timescale on the order of the thermal timescale of the simulation, which introduces some incidental and nonphysical correlation and anticorrelation between portions of the X-ray and intrinsic UV light curves.The Figure 2. Top panels compare the intrinsic and reprocessed UV light curves that come from the one-and two-frequency runs, respectively.Middle and bottom panels compare the reprocessed UV light curves and injected X-rays from the two-frequency runs.τ damp = 30 days for the injected X-ray light curve in the left panels (Run A) and τ damp = 0.1 day for the injected X-ray light curve in the right panels (Run B).In the bottom right panel, the X-ray light curve has been window averaged over 2 days for easier comparison.To make each UV light curve, we sum the total outgoing flux from the low-frequency group over every cell at the top and bottom of the simulation box.For easier comparison, we scale the flux, F, by the standard deviation of the flux, σ(F), and then subtract off the minimum flux for each light curve.Time zero corresponds to when we start injecting the X-rays.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.
overall correlation between the injected X-ray and intrinsic UV light curve in Run A is r = −0.032,p = 0.0032.
The scattering atmospheres in our simulations are similar to those modeled by Gardner & Done (2017), who found that UV emission from NGC 5548 was in good agreement with analytic models assuming a puffed-up photosphere with low X-ray absorption opacities.In addition, AGN disks with scattering atmospheres may help account for the longer-than-anticipated lag timescales of recent AGN reverberation mapping campaigns (Hall et al. 2018).However, other models, such as those in Panagiotou et al. (2022a), suggest that the canonical reprocessing of X-rays through absorption can also produce longer lag timescales and poorly correlated X-ray and UV light curves, like those observed for NGC 5548.In addition, observations using X-ray spectroscopy suggest that AGN disks may not be fully ionized or scattering dominated (e.g., García et al. 2019).
The time-dependent opacities in our simulation are calculated self-consistently from the temperature and density determined by evolving the MHD equations using an initial vertical profile in radiation hydrostatic equilibrium.If our simulation box was located in a different temperature region of the disk or we used different AGN parameters, such as SMBH mass or accretion rate, it is possible that we would have different initial conditions that would lead to an increased absorption or decreased scattering optical depth near the surface of the disk.We have also assumed that the opacity for the X-ray radiation is the proper average of the free-free absorption opacity and electron scattering opacity.This assumption is most accurate if the gas is fully ionized.The disk may not be fully ionized if, for example, the temperature near the disk surface drops below T g ≈ 10 4 K, which it might for different initial conditions.In general, a lower ionization fraction would decrease the scattering opacity and likely lead to increased absorption (García et al. 2013).Using a more sophisticated opacity model for the high-energy photons could also lead to higher absorption opacities.Higher absorption opacities would result in a larger energy transfer from the X-ray irradiation to the gas and likely increase the correlation between the X-ray and UV light curves.Salvesen (2022) found there was a nonzero thermalization time for the reprocessing of X-ray irradiation.We tested whether there is a lag between the X-ray and UV light curves and found that the strongest correlation between the X-ray and UV light curves in Run A is still at <0.08 day.On the other hand, we found a very weak correlation for Run B, r = 0.037 and p = 0.002, occurs for a lag of 32 days, which is roughly the thermal timescale near the surface of the disk.However, a lag of around 30 days for reprocessing has not been observed in disk continuum reverberation mapping campaigns.
Overall, our simulations show that the low absorption opacity for the hard X-ray irradiation makes it unlikely that hard X-rays are the main driver of variability in emitted UV light curves.Several AGN reverberation mapping campaigns have found no correlation (Buisson et al. 2018;Morales et al. 2019, e.g.,) or only weak correlations (e.g., Schimoia et al. 2015;Edelson et al. 2019; We see some evidence that the UV PSDs are damped at f  0.03 day −1 .At 0.03 < f < 0.2 day −1 , the slope of the UV PSDs appear to be well fit by f −2 , shown as the dashed line, which is the DRW slope above the damping frequency.At f > 0.2 day −1 , the UV PSDs are well fit by f −2.5 , shown as the dotted line.Cackett et al. 2023) between X-ray and UV-optical light curves, which agrees with our results here.However, it is possible that this observed lack of correlation in these papers is simply due to absorption of intervening material (e.g., Kara et al. 2021).

Summary and Future Work
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 1988Antonucci , 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. 2016Jiang et al. , 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.
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).

Figure 1 .
Figure 1.Top-right panel: schematic of our simulation and spacetime diagram of the gas temperature in the one-and two-frequency simulations in the regions left and right of the gray line, respectively.The red arrows represent the time-varying X-ray radiation we inject in the upper and lower boundaries in our two-frequency simualtions.Middle right panel: spacetime diagram of the gas density in the one-and two-frequency simulations in the regions left and right of the gray line, respectively.Bottom right panel: same as middle panel but for the radiation energy.We show the UV radiation energy above the midplane and the X-ray radiation energy below the midplane.In all panels, the solid white line is z(τ UV = 1) ≈ z(τ XR = 1), and the dashed white line is z(τ a,UV = 1).τ a,XR < 1 throughout our simulation.Left panels: horizontally and time-averaged vertical profiles of gas temperature, density, and turbulent velocity in the top, middle, and bottom panels, respectively, for our two-frequency simulation.All figures are for Run A.

Figure 3 .
Figure3.The top panels show the PSDs of the intrinsic UV, reprocessed UV, and injected X-ray light curves in Figure2.The bottom panels compare the PSDs of the injected and emitted X-rays.The left (right) panels show the PSDs from Run A (Run B), where τ damp = 30 days (τ damp = 0.1 day) for the injected X-ray light curve.We see some evidence that the UV PSDs are damped at f  0.03 day −1 .At 0.03 < f < 0.2 day −1 , the slope of the UV PSDs appear to be well fit by f −2 , shown as the dashed line, which is the DRW slope above the damping frequency.At f > 0.2 day −1 , the UV PSDs are well fit by f −2.5 , shown as the dotted line.

Table 1 Ω
0 , ρ 0 , T 0 , and H g are Fiducial Units We Use as Normalization Parameters in Our Simulation for Time, Mass, Temperature, and Distance, Respectively