# Liquid-organic time projection chamber for detecting low energy
antineutrinos

T. Radermacher<EMAIL_ADDRESS>J. Bosse  S. Friedrich  M.
Göttsche  S. Roth  G. Schwefer

###### Abstract

The MeV region of antineutrino energy is of special interest for physics
research and for monitoring nuclear nonproliferation. Whereas liquid
scintillation detectors are typically used to detect the Inverse Beta Decay
(IBD), it has recently been proposed to detect it with a liquid-organic Time
Projection Chamber, which could allow a full reconstruction of the particle
tracks of the IBD final state. We present the first simulation-based
statistical analysis of the expected signatures. Their unequivocal signature
could enable a background-minimized detection of electron antineutrinos using
information on energy, location and direction of all final state particles. We
show that the positron track reflects the antineutrino’s vertex. It can also
be used to determine the initial neutrino energy. In addition, we investigate
the possibility to reconstruct the antineutrino direction on an event-by-event
basis by the energy deposition of the neutron-induced proton recoils. Our
simulations indicate that this could be a promising approach which should be
further studied through experiments with a detector prototype.

###### keywords:

Liquid detectors; Time projection Chambers (TPC); Neutrino detectors; Search
for radioactive and fissile materials

††journal: Nuclear Inst. and Methods in Physics Research, A

[1]organization=III. Physikalisches Institut B, RWTH Aachen University,
addressline=Sommerfeldstraße 16, city=Aachen, postcode=52074, country=Germany
[2]organization=Aachen Institute for Advanced Study in Computational
Engineering Science (AICES), RWTH Aachen University,
addressline=Schinkelstraße 2, city=Aachen, postcode=52056, country=Germany

##  1 Introduction

Thanks to important advances in detector technology and analysis techniques,
experimental neutrino research has become a field of precision physics that
continues to push boundaries of sensitivity and accessible energy. Over the
last years, the $\mathrm{MeV}$ region of neutrino energy has been of
particular interest, e.g. for the investigation of reactor, solar and
geoneutrinos.

In addition, the energy domain below $10\,\mathrm{MeV}$ can be useful for
nuclear monitoring purposes where antineutrinos are produced by
$\upbeta$-decaying fission products. With regard to reactors, a number of
studies, for instance [1, 2], have examined antineutrino monitoring to observe
whether reactors are operational, or even to deduce parameters such as fuel
burnup or power. In such applications, detectors would typically be placed in
close vicinity to the reactor building. Another project [3] looks at the
potential for antineutrino monitoring to detect undeclared reactors from some
distance, a significantly more challenging task.

For nuclear waste, both in interim storage facilities and in future geological
repositories, verification activities are required to control the non-
diversion of spent fuel which contains plutonium. For the case that other
verification techniques fail, antineutrino monitoring could act as a
redundancy without common failure modes [4, 5]. The main challenge will be a
low neutrino flux intensity.

Antineutrinos in the $\mathrm{MeV}$-energy range can be detected via their
interaction with a proton producing a neutron and a positron, which had also
been used for the discovery of the neutrino. This reaction, called Inverse
Beta Decay (IBD), has the advantage of a clear final state signature
consisting of a prompt signal of the two $511\,\mathrm{keV}$ photons from the
positron annihilation and a delayed signal from the absorption of the neutron
in a nucleus. Selecting this event signature reduces background sources to a
large extent.

Typically, organic scintillation detectors are used to measure the energy
depositions of the prompt and delayed signals. For high fluxes, e.g. adjacent
to a nuclear reactor, no extensive shielding is necessary [6, 7]. It has been
shown that using IBD is also effective enough for measuring moderate neutrino
fluxes in some distance to the reactors in neutrino oscillation experiments
[8, 9, 10], though requiring much larger detectors and a more significant
overburden that shields against cosmic radiation. Lastly, the same method even
allows the measurement of the very small flux of geoneutrinos that has been
successfully performed by Borexino [11]. Here, not only very large amounts of
overburden are necessary, but also further measures against background sources
such as veto counters and shielding against intrinsic radiation of the
surrounding rocks are required. In addition, the scintillator needs to be of
high purity to avoid radioactive decays inside of the detector.

These elaborate precautions against background are needed due to the fact that
in liquid scintillation detectors the final state of the IBD reaction is not
reconstructed in detail. Instead, the measurements are based on the three
following observables of the event: The prompt energy from the positron
annihilation, the delayed energy caused by the neutron absorption and the time
between both energy depositions. For instance, IBD events cannot be
distinguished from events involving electrons and neutrons (e.g. from
cosmogenically produced ${}^{9}\mathrm{Li}$ or ${}^{8}\mathrm{He}$) as
annihilation photons cannot be identified separately, because the prompt
energy of the event, the kinetic energy of positron and the energy from the
positron annihilation, is usually measured as a whole.

In this paper, we propose a different detection method that could have the
benefit of strongly reduced background, namely detecting IBD events using the
tracking capabilities of a Time Projection Chamber (TPC) to fully reconstruct
the ionization of each particle of the IBD final state. This is to be achieved
using an organic liquid. It contains sufficient free protons to serve as a
target for IBD reactions, as has already been proposed by previous studies
[12, 13]. Specifically, in each IBD event one could seek to identify the
initial positron track, the energy depositions of the two annihilation photons
and of the delayed photon emitted due to the neutron absorption by a nucleus.
Furthermore, ionization clusters induced by the neutron scattering off the
nuclei of the detector material might be detectable.

Reconstructing the tracks and showers from the original individual hits could
provide information about the energy, the location and, with limitations, also
the direction of the final state particles. The unequivocal signature of a
completely resolved IBD final state could then enable background-minimized
detection of electron antineutrinos. This would then allow to use the IBD
antineutrino detection also for sources with smaller neutrino fluxes like
nuclear waste or geoneutrinos and also in locations with less or even without
any overburden.

An additional benefit could be the access to information on the initial
direction of the detected antineutrino. Liquid scintillation detectors have
also looked for final state variables that contain information on the
antineutrino direction [14]. However, their position reconstruction usually
lies in the order of few centimeters due to the time resolution achievable
from the photodetectors surrounding the volume. Additional segmentation of the
detector volume can further improve the position reconstruction and by that
the directional sensitivity as discussed for SANDD [15]. Nonetheless, the
results are based on the statistical analysis of a large number of
antineutrino events. By using the final state as measured with a TPC, the
variables sensitive to the antineutrino direction could be reconstructed with
higher precision. New variables like the position of the neutron absorption or
even neutron scattering with respect to the positron track become available
and could be included in the reconstruction of the antineutrino direction.
Hence, the antineutrino direction could possibly be measured on an event-by-
event basis. If the directionality information can be improved sufficiently by
this technique, one could be able to discriminate antineutrino background from
the signal, e.g. reactor antineutrino background when measuring nuclear waste.
The same holds for the determination of the distribution of radioactive
nuclides within the rock formations when studying geoneutrinos.

In this paper, we present results of the first simulation-based statistical
analysis of the expected IBD signatures in a liquid-organic TPC. The goal is
to identify opportunities and challenges for their reconstruction. We use the
example of antineutrinos emitted from nuclear waste as they present a
significant challenge due to their energies reaching barely above the IBD
detection threshold. Particular emphasis is placed on examining whether the
antineutrino direction might indeed be reconstructable, as this would clearly
be the most challenging reconstruction task. Compared to other experiments it
is clear that there will be practical challenges to realize such detector.
This study looks at potential capabilities of the general detection approach
we present, in order to guide future experimental work to further assess the
concept’s practical feasibility.

##  2 Antineutrinos from Nuclear Waste

In the first few hundreds of years, the dominant radioactivity in nuclear
waste comes from the $\upbeta$-decaying elements like ${}^{90}\mathrm{Sr}$ and
${}^{137}\mathrm{Cs}$. For the evaluation of antineutrino emissions, we
obtained the isotopic concentrations of radioactive isotopes in spent fuel
discharged from a pressurized water reactor, PWR, using a simulation with the
reactor simulation code SERPENT2 [16]. A burn-up of $55\,\mathrm{MWd/kg}$ and
a uranium enrichment of 4% was implemented. We then used the beta spectra of
all beta-decaying nuclides from the software tool BetaShape [17] and converted
them to the antineutrino spectra shown in figure 1 for different ages of the
nuclear waste.

Figure 1: Energy spectra of antineutrinos emitted from one fuel assembly of a
simulated PWR reactor with a burnup of $55\,\mathrm{MWd/kg}$ and enrichment of
4%.

The antineutrinos are emitted at low energies with the emission significantly
decreasing at $2.3\,\mathrm{MeV}$. The maximum energy is reached at
$4.3\,\mathrm{MeV}$ with an emission rate more than ten orders of magnitude
smaller than for the lower energy range up to $2.3\,\mathrm{MeV}$. Considering
the detection via IBD with its threshold at $1.8\,\mathrm{MeV}$, the energy
region that can be detected is fairly small. In this region, the antineutrino
emission is dominated by the decay of the ${}^{90}\mathrm{Y}$ isotope, which
is a decay product of the highly abundant fission product ${}^{90}\mathrm{Sr}$
and stays in equilibrium with it.

Spent nuclear fuel is usually stored in a sealed containment that shields most
of the ionizing radiation emitted from the waste. This shielding, however, has
no impact on the antineutrino flux due to their low interaction cross
sections. This makes the antineutrinos very interesting to be used to monitor
the nuclear waste. Furthermore, the emitted antineutrino flux is directly
depending on the waste mass, providing the opportunity for re-verification if
the continuity of knowledge is lost.

In A an exemplary scenario of a monitoring application for a geological
repository is shown. It presents the expected event rates in an active medium
of Tetra-Methyl-Silane as studied in this work.

##  3 Liquid-organic time projection chamber

The TPC is a widely used tracking detector due to its good track
reconstruction ability and energy loss resolution. It usually consists of a
large detection volume in which charged particles deposit energy by ionizing
the media and liberating electrons. Due to an externally applied electric
field these ionization electrons drift towards an anode plane. A three-
dimensional reconstruction is provided by a segmented, two-dimensional readout
at the anode and the measurement of the drift time of the electrons.

Originally, TPCs use a gas mixture as sensitive material to measure the
passage of charged particles. The neutral neutrino cannot be detected directly
by an ionization detector, but charged final state particles from neutrino
interactions with the detector material can be measured. Due to the extremely
small cross section of neutrino interactions, a denser material than a gas is
necessary. For this reason, shortly after the invention of the gaseous TPC,
the usage of liquid argon as active medium was suggested because of its higher
density [18]. Nowadays, using liquefied noble gases has become an established
technology for the detection of weakly interacting particles. In neutrino
experiments like MicroBooNE [19] and ArgonNeuT [20] liquid argon TPCs have
been exploited, whereas liquid xenon is used in dark matter searches, e.g. the
experiments XENON100 [21], LUX [22] and PandaX [23]. Currently, the DUNE
neutrino experiment is building huge liquid argon TPCs for its far detector
with a combined 70 kiloton of liquid argon as the target mass [24]. First
results with this technology were gained at the ProtoDUNE setup of the CERN
neutrino platform [25].

One drawback of liquid noble gases is that they need to be cooled down by a
cryogenic system and therefore, it is very challenging to construct and
operate a detector system with this technology. Hence, studying alternative
high-density media that allow drifting electrons at room temperatures has been
of special interest for a long time. One promising material is Tetra-Methyl-
Silane, TMS, which is a non-polar organic liquid with spherically shaped
molecules. Its usage in ionization detectors has been previously investigated
[26] with the focus on its application in calorimetry. For example, KASCADE-
Grande [27] used liquid ionization chambers filled with TMS in their central
calorimeter. In recent years, the use of these organic liquids as sensitive
material of a TPC has been discussed [12, 13]. By operating small test
chambers, the charge transport in organic liquids has been studied by
measuring the charge yield [28], the drift velocity [29] and the mean free
drift length [30]. Still, these studies have to be intensified as open
questions remain. For instance, further parameters relevant for the electron
drift in organic liquids, such as the diffusion coefficients, still have to be
investigated.

Another advantage of an organic liquid with respect to liquid noble gases is
the presence of free protons due to its hydrogen content. Through this, a low
threshold for the IBD process can be reached, making it possible to detect
antineutrinos down to the energy of $1.8\,\mathrm{MeV}$. Event rates should
also be comparable with liquid scintillation detectors. For example, the
number of target nuclei for IBD per volume in TMS is as large as 81 % of the
one in the scintillator material of KamLAND [31].

As explained in the section above, the antineutrinos emitted from spent
nuclear fuel have energies in the range between $1.8\,\mathrm{MeV}$ and
$2.3\,\mathrm{MeV}$. Hence, liquid noble gases can not be used as target
material to detect them via IBD due to their higher energy threshold of this
reaction.

##  4 Simulation setup and event generation

In order to investigate the reconstruction of IBD events in a liquid-organic
(LOr) TPC, we simulated antineutrino interactions within Tetra-Methyl-Silane
as target material using the event generator GENIE [32] and the energy
spectrum of figure 1 at 100 years. Since the flux effectively vanishes above
$2.3\,\mathrm{MeV}$, IBD can only take place at the hydrogen nuclei. The
outgoing positron and neutron from the IBD reaction are then fed to a GEANT4
[33] simulation to investigate their interaction with TMS. From this
simulation we obtain the energy deposition of each particle. These are then
translated into the number of electrons with an estimated electron yield
$G_{if}$ that are used to simulate the spatial charge distribution at an anode
with a readout structure of squared pads with a size of $5\,\mathrm{mm}\times
5\,\mathrm{mm}$ considering the drift velocity $v_{d}$ and the diffusion
coefficient for longitudinal $d_{L}$ and transverse diffusion $d_{T}$.

We examine two extreme, but still reasonable detector configurations with
different maximum drift lengths of $l_{1}=1\,\mathrm{m}$ and
$l_{2}=10\,\mathrm{m}$, respectively, where we keep a sensitive volume of
$80\,\mathrm{m^{3}}$ congruous to our study presented in the A. For both
configurations we assume a cathode high voltage of $-500\,\mathrm{kV}$ which
generates a drift field of $E_{1}=5\,\mathrm{kV/cm}$ and
$E_{2}=0.5\,\mathrm{kV/cm}$, respectively. For the first detector
configuration, detector configuration 1, $G_{if,1}=7\,\mathrm{e^{-}/keV}$ [34]
and $v_{d,1}=5.5\,\mathrm{\upmu m/ns}$ [30] are used. For the second
configuration, detector configuration 2, we assume
$G_{if,2}=5\,\mathrm{e^{-}/keV}$ and $v_{d,2}=0.5\,\mathrm{\upmu m/ns}$ using
a mobility of $\mu=105\,\mathrm{cm^{2}/V/s}$ [35].

Due to the lack of experimental data, we use the Einstein Relation
${D_{L,T}=\mu\cdot k_{B}T/\mathrm{e}}$ to estimate the thermal limit of the
diffusion coefficient $D_{L,T}$ with given mobility $\mu$. To express the
diffusion as a function of the drift distance $l_{d}$ instead of drift time
$t_{d}$ via the formula $\sigma=\sqrt{2\cdot D_{L,T}\cdot
t_{d}}=d_{L,T}\sqrt{l_{d}}$ we transform the diffusion coefficient using
equation 1:

$d_{L,T}=\sqrt{\frac{2D_{L,T}}{v_{d}}}=\sqrt{\frac{2k_{B}T}{\mathrm{e}E}}.$
(1)

Measurements of electron diffusion in liquid argon show that the actual
diffusion coefficient $D_{L,T}$ is a factor 3–4 larger than the thermal limit
[36] which can be explained by an energy increase of the electrons due to the
external electric field [37]. Such an increase should also appear in TMS
although measurements have to be performed to find out the magnitude. However,
for simplicity we assume the same factor 4 for TMS and get a diffusion
coefficient of $d_{L,T,1}\approx 60\,\mathrm{\upmu m/\sqrt{cm}}$ for detector
configuration 1 and $d_{L,T,2}\approx 190\,\mathrm{\upmu m/\sqrt{cm}}$ for
detector configuration 2. We use the maximum drift length to study the maximum
impact of diffusion which in the end limits to which extent tracks can be
separated from each other. A summary of the simulation parameters can be found
in table 1.

Table 1: Simulation parameters for the two detector configurations. In both
cases a high voltage of $500\,\mathrm{kV}$ and a sensitive volume of
$80\,\mathrm{m^{3}}$ is assumed.
|  Detector configuration 1 |  Detector configuration 2
---|---|---
maximum drift length $l$ | $1\,\mathrm{m}$ | $10\,\mathrm{m}$
electric field $E$ | $5.0\,\mathrm{kV/cm}$ | $0.5\,\mathrm{kV/cm}$
drift velocity $v_{d}$ | $5.5\,\mathrm{\upmu m/ns}$ | $0.5\,\mathrm{\upmu m/ns}$
diffusion coefficient $d_{L,T}$ | $60\,\mathrm{\upmu m/\sqrt{cm}}$ | $190\,\mathrm{\upmu m/\sqrt{cm}}$
electron yield $G_{if}$ | $7\,\mathrm{e^{-}/keV}$ | $5\,\mathrm{e^{-}/keV}$

With these detector scenarios we want to evaluate two extreme, but still
reasonable combinations of electric field and drift lengths. Detector
configuration 1 has a higher electric field which comes with a higher drift
velocity and smaller diffusion coefficient. Detector configuration 2 on the
other hand has a larger drift length which results in less readout channels
when keeping the same sensitive volume. Both scenarios are limited by the
maximum high voltage. Even though we know that achieving this is challeging,
we consider $500\,\mathrm{kV}$ as still manageable. In comparison, for the
planned dual-phase detector of the DUNE project a cathode voltage of
$600\,\mathrm{kV}$ was envisaged [38]. However, due to the design evolution of
the ProtoDUNE dual-phase to Vertical Drift only $300\,\mathrm{kV}$ are needed
[39]. Their stability could be demonstrated in a HV campaign [40].

Another important point that has to be kept in mind is that in liquids,
organics or cryogenic noble gases, high purity is needed to reach a sufficient
free electron lifetime. For TMS this electron lifetime is constant for
electric fields below $E=20\,\mathrm{kV/cm}$ according to [41]. With the drift
velocity rising proportional to the electric field a balance needs to be found
between electric field and drift length. However, measurements have to be
carried out to show if the required electron lifetimes are reachable. Electron
lifetimes required for detector configuration 1 with the short drift length
have already been shown [42]. For detector configuration 2 with the long drift
distance much larger electron lifetimes have to be achieved, comparable to the
electron lifetimes in liquid Argon TPCs. Even though it is not yet proven if
this can be achieved, the simulation of configuration 2 should show the
influence of the diffusion. For this purpose these studies were performed with
infinite electron lifetime.

##  5 Inverse beta decay signature in a LOr-TPC

Figure 2: Simulation of an IBD event projected onto the anode for detector
configuration 1. The left side shows the distribution of the single electrons
in a color-coded scale. The right side shows which particle is associated to
the ionization clusters.

Figure 2 shows the interaction of an antineutrino entering the TPC
perpendicular to the drift direction for detector configuration 1. It shows
all energy depositions of the secondary particles projected on the anode
integrated over the whole event time. On the left side, for each readout pad
the recorded ionization in number of electrons is shown. On the right side,
the associated particle is mapped to each of the ionization clusters. Contrary
to liquid scintillation detectors, which can only distinguish a prompt from a
delayed signal, here also the ionization of the positron and the scattering of
the neutron is directly accessible.

For the comparison of the detector configurations, figure 3 depicts a zoomed
view into the projected charge distributions simulated for the same event with
the respective drift parameters. It shows the ionization cluster of the
positron in the center, the ionization from the neutron-induced proton recoils
on the right, the neutron capture recoil on top and some ionization of one of
the annihilation photons on the bottom right. For detector configuration 2
that is shown on the right side of figure 3, the charge clusters are spread
over more pads as expected due to the larger diffusion coefficient and the
longer drift distance. As a result, the clusters of the positron and the
neutron ionization are more likely to overlap and therefore, lead to a poorer
separation. In addition, the larger spread of the charge reduces the potential
electric pad signals that are especially low for the small ionization signals
of the proton recoils which is slightly lower for detector configuration 2.

Figure 3: Comparison between the two detector configurations using the same
event. The charge distribution of detector configuration 1 is shown on the
left side. The right side shows the distribution of detector configuration 2.
Both show the same zoomed-in view in the center of the event also shown in
figure 2.

We will explain each signal created by the particles of the IBD final state in
the following. The numbers are given for detector configuration 1 first and
for detector configuration 2 in brackets behind. In the simulated energy
region, the positron receives a kinetic energy smaller than
$0.5\,\mathrm{MeV}$. This leads to a typical track length below
$1\,\mathrm{mm}$ before annihilation. This track starts at the primary vertex
of the IBD interaction and can be clearly identified by the up to
$3200\,(2286)\,\mathrm{electrons}$ that are created along its path. The
ionization released along the positron track is a direct measure of the
positrons’ initial kinetic energy, which in turn measures the energy of the
antineutrino as depicted in figure 4.

Figure 4: Number of ionization electrons created by the positron track in
detector configuration 1 as function of the antineutrino energy.

At rest, the positron annihilates with an electron of the detector material
and two $511\,\mathrm{keV}$ photons are created which can be used as a prompt
signal to identify the IBD event as also done in liquid scintillation
detectors. As shown in the event display, these photons produce characteristic
energy depositions due to the fact that they initially travel a few
centimeters losing their energy by Compton scattering before the final
photoelectric effect takes place.

During its propagation, the neutron scatters off the nuclei of the detector
liquid until it is thermalized and finally captured. The nucleon created by
the neutron capture is produced in an excited state and finally releases the
excitation energy via photon emission. From 100,000 simulated events we found
that in 93 % of the cases the neutron is captured by a hydrogen nucleus which
releases a single photon of $2.2\,\mathrm{MeV}$ energy. Again, this photon
loses its energy by Compton scattering, traversing the liquid for a few
centimeters. If the neutron is captured by another nucleus of the detector
liquid, several photons can be emitted with energies according to the
excitation levels of the produced nucleus. On average it takes $(244.2\pm
1.2)\,\mathrm{\upmu s}$ until the neutron is captured and the photon is
released. Detecting the photon from the neutron capture process in addition to
the positron annihilation photons gives a clear event signature.

Even though the photons of the annihilation and the capture of the neutron are
emitted at different times, the arrival time of the two charge signals at the
anode can be ambiguous. First, the spatial distribution of the photon
absorptions results in different arrival times of the ionization charge at the
anode. Second, the electron cloud of this ionization diffuses while drifting
to the anode. This can lead to an overlap in the arrival time of the signal
making it harder to distinguish between the positron and the neutron signal.
For detector configuration 1 28 % of the events have electrons liberated by
the neutron capture arriving at the readout before all electrons of the
annihilation have arrived. For detector configuration 2 82 % of the events
show a concurrent arrival of the electrons due to the larger drift distance
and the slower drift velocity. In most cases, the simultaneous arriving charge
clusters, however, should be assignable to the different processes because of
the spatial separation in the readout plane. Modern intelligent reconstruction
algorithms could then reconstruct the events and may resolve the ambiguities,
although this has to be studied in more detail.

The neutron produced by the IBD has typical kinetic energies between
$1\,\mathrm{keV}$ and $3\,\mathrm{keV}$ due to its larger mass compared to the
positron. Therefore, the positron obtains practically all of the energy of the
incoming antineutrino. At the IBD energy threshold, the positron receives no
kinetic energy and since the proton is at rest before the interaction, the
neutron takes over the complete momentum of the antineutrino. Hence, it also
has the identical direction as the incoming antineutrino. With increasing
neutrino energy and thus rising positron momentum, the deviation between
neutron and antineutrino direction increases due to momentum conservation
[43]. Still, the neutron direction is largely correlated with the antineutrino
direction at the considered energies.

Figure 5: Number of electrons created by neutron-induced proton recoils in
detector configuration 1. The left side shows the number of electrons created
in each proton recoil. The right side shows the integrated number of electrons
from proton recoils in each IBD event.

Afterwards, the neutron changes its direction due to scattering from the
protons in the detection material. The mean distance between the antineutrino
vertex and the first neutron-induced proton recoil is $(9.20\pm
0.04)\,\mathrm{mm}$. Usually, the largest energy transfer occurs in the first
collision. It decreases with each collision as shown on the left side of
figure 5. The first proton recoil creates on average $8.4\,(6)$ electrons
while it is only $3.9\,(2.8)$ for the second and $1.8\,(1.3)$ for the third.
As the right side of figure 5 shows, in each event an integrated signal of
5–30 (4–21) electrons can be expected from the neutron-induced proton recoils.
The most probable value is around 10 (7).

The various physical processes discussed so far yield very different amounts
of ionization to be detected. Especially the detection of the neutron signal
consists of only a few electrons and is therefore very challenging. However,
different from dark matter experiments, this small signal is not be used as a
trigger signal, but is intended to reconstruct a satellite charge cluster near
the cluster created by the positron. To register the small charge amount of
the neutron signal the concept of a dual-phase TPC can be utilized. Here,
adding a gaseous layer on top of the liquid drift volume allows for electron
amplification in the gas due to an avalanche process induced by a high
electric field applied in this region. This concept of dual-phase TPCs has
been successfully achieved in liquid xenon TPCs for dark matter searches [44].

##  6 Investigations on the reconstruction of the antineutrino direction

An additional benefit of the full reconstruction of IBD final states could be
the access to information on the initial direction of the detected
antineutrino on an event-by-event basis.

We investigated IBD events simulated with GENIE for the nuclear waste
antineutrino spectrum. The left side of figure 6 shows the angle between the
incident antineutrino direction and the initial direction of the neutron as
function of the antineutrino energy. Directly at the IBD threshold, the
initial neutron direction is identical to the incident antineutrino direction.
For higher antineutrino energies $E_{\nu}$ the deviation between the neutron
and antineutrino direction increases. The maximum angle can be described
analytically, see equation 2 [43].

$\cos{\phi_{\texttt{max}}}=\frac{\sqrt{2E_{\nu}\Delta-(\Delta^{2}-m_{e}^{2})}}{E_{\nu}}\qquad\text{with}~{}\Delta=m_{n}-m_{p}.$
(2)

Figure 6: Angular distribution of the neutron as a function of the incident
antineutrino direction. The left side shows the angular deviation as a
function of the antineutrino energy. The right side shows the angular
probability distribution of IBD events based on the nuclear waste spectrum
that we consider.

Since the most probable angle between neutron and antineutrino corresponds to
the maximum angle (as can be seen in figure 6), the neutron is initially
directed in a cone around the antineutrino direction. Later the neutron
undergoes several scattering events. The correlation with the initial neutron
direction and the antineutrino direction is reduced by each scattering
process. Therefore, it would be desirable to reconstruct the location of the
first or of the first few neutron-proton scatterings. The connecting line
between this location and that of the positron track would then yield a very
good prediction of the direction of the antineutrino.

On the right side of figure 6, the angle between the neutron and antineutrino
direction is shown for the IBD events based on the nuclear waste antineutrino
spectrum considered by us. The maximum angle is $21^{\circ}$ and the
distribution peaks at around $15^{\circ}$.

Furthermore, we studied the potential to reconstruct the antineutrino
direction via the energy depositions using our GEANT4 Monte Carlo simulation.
In this study, we considered five different definitions for a neutron
direction vector that always starts at the IBD interaction origin and compared
it with the incident antineutrino direction vector. We investigated (1) the
first neutron-induced proton recoil, (2) the energy-weighted mean of all
proton recoil ionization clusters, (3) all recoils including that of the
nucleus capturing the neutron, (4) only the recoil of the nucleus capturing
the neutron, and (5) the ionization of the first Compton-scatter of the photon
released after the neutron capture. As depicted in figure 7, the angular
correlation between the direction of the initial antineutrino and the
different definitions (1-5) diminishes significantly with more ionization
energy taken into account. For (4) and (5) the distributions even become
isotropic, so these events can not be used to reconstruct the direction on an
event-by-event basis.

Figure 7: Angle between reconstructed neutron direction and incident
antineutrino direction. Distributions are shown for different reconstruction
methods, as explained in the text.

Therefore, we conclude that for the best possible reconstruction of the
antineutrino direction only the first neutron-induced proton recoil should be
considered. However, as it might not be possible in all cases to isolate this
energy cluster from the signature arriving at the anode, the second best
option would be to use the energy-weighted position of all proton recoils
without the recoil of the capturing nucleus.

As we have already shown above, the deposited energy on average decreases with
each neutron scattering. Hence, automatically only the first few recoils will
result in an measurable energy deposition and would be taken into account. As
our simulations show, the capture recoil is usually $\Delta t\approx
240\,\mathrm{\upmu s}$ delayed and can therefore be excluded from the other
proton recoils.

Finally, we studied the reconstruction of the neutron direction from the
charge distribution created by the drifting electrons at an anode using the
drift simulation described above. The question we pose is to which extent the
neutron recoil cluster(s) can be separated spatially from the positron track,
identifying them as small satellite clusters. Since the positron track is
rather short, it can also be considered as a single ionization cluster from
which we determine the center of charge and its two-dimensional standard
deviation. We then determine the same for both the ionization clusters of all
neutron-induced recoil protons and for only the first one. For IBD events
where the antineutrino direction is perpendicular to the drift direction of
the LOr-TPC we found that for 87 % of the events, the charge distributions are
more than two combined standard deviations apart in detector configuration 1
scenario as depicted in figure 8 using only the first neutron-induced proton
recoil. For 88 % of these events, the neutron direction can be reconstructed
with a deviation of less than $10^{\circ}$ from the true neutron direction.
When using all proton recoils for the neutron reconstruction, only 72 % have a
deviation of more than two combined standard deviation. In 68 % of these
events, the neutron direction can be reconstructed with less than a
$10^{\circ}$ deviation.

Figure 8: Distance of the charge distributions of neutron and positron on an
anode for detector configuration 1 in combined standard deviations on the left
side. On the right side the resulting angular deviation of the reconstructed
and the true neutron direction for events with more than two combined standard
deviation is shown. Here only the charge cluster of the first neutron-induced
proton recoil was used.

For detector configuration 2, the charge distributions are only separated for
29 % of the events with more than two standard deviations considering only the
first neutron-induced proton recoil. From these events, 54 % can be
reconstructed with less than $10^{\circ}$ deviation. In case the charge
distribution is reconstructed by all proton recoils 33 % of the events are
more than the two standard deviations apart in this scenario, of which 55 %
can be reconstructed with less than $10^{\circ}$ deviation.

Detector configuration 1 shows much better performance than detector
configuration 2. This is mainly due to the larger diffusion in the case of
detector configuration 2. On the one hand the diffusion coefficient is more
than three times smaller than for detector configuration 2. On the other hand
its maximum drift distance that we used for our evaluation is ten times
larger. We additionally evaluated the separability for both detector
configurations using a diffusion coefficient of $d_{L,T}=300\,\mathrm{\mu
m/\sqrt{cm}}$. For detector configuration 1, this reduces the number of events
that can be separated with more than two combined standard deviations to 49 %
which would still be in an acceptable range. Due to the large drift distance
in detector configuration 2, the clusters smear out that much, that they can
not be separated anymore. We conclude that measurements of the diffusion
coefficients are needed to specify a maximum drift length.

For our example, assuming a detector volume of dimensions $10\times 4\times
2\,\mathrm{m^{3}}$, the setup of detector configuration 1 with the
$1~{}\mathrm{m}$ maximum drift length can be realized by a central cathode in
parallel to the $10\times 4\,\mathrm{m^{2}}$ large side. This setup would
require a total readout area of $80\,\mathrm{m^{2}}$. The setup of detector
configuration 2 could be realized with a cathode on one of the $2\times
4\,\mathrm{m^{2}}$ large sides. In summary, detector configuration 1 delivers
better results due to the ten times shorter drift length, but needs a ten
times larger readout area instrumented.

##  7 Conclusion and Outlook

We presented the first simulation-based statistical analysis of the expected
IBD signatures in a liquid-organic TPC. This is a newly proposed detection
method for neutrino scattering events using a liquid-organic TPC which aims
for a full reconstruction of the IBD final state. The identification of this
unequivocal event signature has the potential to enable a background-minimized
detection of antineutrino events. The basis for this is the reconstruction of
the positron track. On the one hand, it improves the neutrino energy
resolution due to the correlation of ionization of the positron track to the
initial antineutrino energy. On the other hand, it makes it possible to
distinguish the IBD events containing a positron signal plus two annihilation
photons from background events involving $\upbeta$-electrons.

Based on our results, we conclude that it could be possible to indeed
associate the energy depositions to the final state particles and achieve a
complete reconstruction of IBD events. Further studies on the identification
with dedicated reconstruction algorithms have to be performed. Recent advances
in these analysis techniques are promising in this respect.

In addition, we showed that the direction of the incoming antineutrino could
be estimated from the locations of the energy depositions caused by the recoil
protons produced by the scattering of the final state neutrons. Even though
their ionization signal is very small, it might be detectable as a satellite
cluster to the positron track. Experimental studies have to be undertaken to
decide whether this is achievable, given the fact that the positron track’s
signal is two orders of magnitude larger than that of the recoil protons.
Additionally, the separability is strongly depending on the diffusion
coefficients and the drift distance. Since no measurement data for the
diffusion coefficient are available, we could only estimate them for our
studies. The concept could lead to a measurement of the antineutrino direction
on an event-by-event basis, even when our simulations show that not all proton
recoils have a sufficient distance from the positron track. It would in any
case, however, be a large improvement compared to current capabilities using
scintillation detectors.

Future simulation studies should also address signatures for antineutrinos
entering the detector in parallel to the electron drift direction. Here, the
neutrons and positron tracks would be much closer together in the anode plane.
How they could be distinguished in a three dimensional reconstruction also
using the drift time measurement should be further assessed.

Finally, before a liquid-organic TPC as it is considered in our studies can be
realised, some technical specifications still have to be demonstrated
experimentally. Most important is to reach electron lifetimes that allow drift
lengths in the order of meters. To achieve this goal, the liquid needs to be
purified to a ppb-level of oxygen. Here, methods deployed at the purification
of liquefied noble gas TPC can be utilized. Furthermore, to realise a dual-
phase readout for the signal registration, the extraction of electrons from
the liquid into a gas needs to be shown. Otherwise, only the larger signals
from the IBD signature can be detected.

We will soon start further studies based on a TPC prototype to investigate the
technical realisation of a detector described in the simulation studies in
this paper. We intend to measure drift velocity, electron lifetime and the
diffusion parameters of TMS. Additionally, we want to examine the extraction
of electrons from the liquid into a gaseous amplification structure.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgments

We thank Madalina Wittel for her engagement in starting this project and we
acknowledge the contributions of David Betram, Helge Haveresch, Anike Ohm and
Hagen Weigel. Additionally, we thank Yan-Jie Schnellbach for fruitful
discussions about this publication. The project was funded by the Federal
Ministry of Education and Research (BMBF) and the Ministry of Culture and
Science of the German State of North-Rhine Westphalia (MKW) under the
Excellence Strategy of the Federal Government and the Länder [funding code
OPSF 573]. Further support came from the Federal Ministry for Economic Affairs
and Energy [funding code 02W6281]. Malte Göttsche and Thomas Radermacher are
generously funded by the VolkswagenStiftung Freigeist Fellowship.

##  Appendix A Initial study for monitoring a geological repository with a
liquid-organic TPC

As an initial study we investigated a model loosely related to the now
discarded Gorleben site [45] in Germany. It comprises 2240 containers with ten
spent fuel assemblies each. The production year of the fuel was simulated
corresponding to the German program, though we assumed the last fuel to be
discharged in 2020.

The layout foresaw the storage caverns to be located at a depth of
$870\,\mathrm{m}$. In our assessment, we envisaged that an antineutrino
detector would be placed at the same depth as the repository, but
$50\,\mathrm{m}$ away from the outermost containers. In view of the
transportation and deployment requirements, we considered a detection volume
of $80\,\mathrm{m^{3}}$ of Tetra-Methyl-Silane for the antineutrino detector.
The results of this preliminary study are summarised in table 2. It can be
seen that the expected detection rate diminishes significantly with time due
to the decay of ${}^{90}\mathrm{Sr}$. Note that larger detectors or multiple
detectors at different locations could be employed, including above the
repository [46]. These results show that antineutrino monitoring with a
liquid-organic TPC could be possible, but indeed only using a detection
technique with high background suppression.

Table 2: Expected detection rates of the antineutrino emissions without
background evaluated with a detector volume of $80\,\mathrm{m^{3}}$ TMS in a
TPC.
Time elapsed since 2020 in years |  Number of IBD events detected in 3 months
---|---
10 | 761
25 | 532
50 | 293
100 | 89
125 | 49
150 | 27

## References

  * [1] R. Carr, et al., Neutrino-Based Tools for Nuclear Verification and Diplomacy in North Korea, Science & Global Security 27 (1) (2019) 15–28.  doi:10.1080/08929882.2019.1603007.
URL https://doi.org/10.1080/08929882.2019.1603007

  * [2] E. Christensen, et al., Antineutrino Monitoring for Heavy Water reactors, Phys. Rev. Lett. 113 (2014) 042503\.  doi:10.1103/PhysRevLett.113.042503.
URL https://link.aps.org/doi/10.1103/PhysRevLett.113.042503

  * [3] M. Askins, et al., The Physics and Nuclear Nonproliferation Goals of WATCHMAN: A WAter CHerenkov Monitor for ANtineutrinos (2015).  arXiv:1502.01132, doi:10.48550/arXiv.1502.01132.
URL https://arxiv.org/abs/1502.01132

  * [4] V. Brdar, P. Huber, J. Kopp, Antineutrino Monitoring of Spent Nuclear Fuel, Phys. Rev. Applied 8 (2017) 054050\.  doi:10.1103/PhysRevApplied.8.054050.
URL https://link.aps.org/doi/10.1103/PhysRevApplied.8.054050

  * [5] M. Wittel, M. Göttsche, Antineutrino Detection Techniques for Monitoring Long-Term Geological Repositories, ESARDA Bulletin 60 (2020) 20–27.  doi:10.2760/840092.
URL https://publications.rwth-aachen.de/record/793192

  * [6] J. Ashenfelter, et al., First Search for Short-Baseline Neutrino Oscillations at HFIR with PROSPECT, Phys. Rev. Lett. 121 (2018) 251802.  doi:10.1103/PhysRevLett.121.251802.
URL https://link.aps.org/doi/10.1103/PhysRevLett.121.251802

  * [7] H. Almazán, et al., Accurate Measurement of the Electron Antineutrino Yield of U235 Fissions from the STEREO Experiment with 119 Days of Reactor-On Data, Phys. Rev. Lett. 125 (20) (2020).  doi:10.1103/physrevlett.125.201801.
URL http://dx.doi.org/10.1103/PhysRevLett.125.201801

  * [8] H. de Kerret, et al., The Double Chooz antineutrino detectors (Jan 2022).  arXiv:2201.13285.
URL https://arxiv.org/abs/2201.13285v1

  * [9] F. An, et al., The detector system of the Daya Bay reactor neutrino experiment, Nucl. Instrum. Meth. A 811 (2016) 133–161.  doi:10.1016/j.nima.2015.11.144.
URL https://www.sciencedirect.com/science/article/pii/S0168900215015636

  * [10] J. Ahn, et al., RENO: An Experiment for Neutrino Oscillation Parameter $\theta_{13}$ Using Reactor Neutrinos at Yonggwang (2010).  arXiv:1003.1391.
URL https://arxiv.org/abs/1003.1391

  * [11] M. Agostini, et al., Comprehensive geoneutrino analysis with Borexino, Phys. Rev. D 101 (2020) 012009.  doi:10.1103/PhysRevD.101.012009.
URL https://link.aps.org/doi/10.1103/PhysRevD.101.012009

  * [12] N. McConkey, Y. Ramachers, Liquid Scintillator Time Projection Chamber Concept, Nucl. Instrum. Meth. A 718 (2013) 459–461.  arXiv:1205.5695, doi:10.1016/j.nima.2012.11.047.
URL https://arxiv.org/abs/1205.5695v1

  * [13] J. Dawson, D. Kryn, Organic Liquid TPCs for Neutrino Physics, JINST 9 (2014) P07002.  arXiv:1405.1308, doi:10.1088/1748-0221/9/07/P07002.
URL https://iopscience.iop.org/journal/1748-0221

  * [14] E. Caden, Studying Neutrino Directionality with Double Chooz (2012).  arXiv:1208.3628.
URL https://arxiv.org/abs/1208.3628v1

  * [15] F. Sutanto, et al., SANDD: A directional antineutrino detector with segmented 6Li-doped pulse-shape-sensitive plastic scintillator, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 1006 (2021) 165409.  doi:https://doi.org/10.1016/j.nima.2021.165409.
URL https://www.sciencedirect.com/science/article/pii/S0168900221003934

  * [16] J. Leppänen, et al., The Serpent Monte Carlo Code: Status, Development and Applications in 2013, Annals of Nuclear Energy 82 (2014) 06021.  doi:10.1051/snamc/201406021.
URL https://doi.org/10.1051/snamc/201406021

  * [17] X. Mougeot, Betashape: A new code for improved analytical calculations of beta spectra, EPJ Web Conf. 146 (2017) 12015.  doi:10.1051/epjconf/201714612015.
URL https://doi.org/10.1051/epjconf/201714612015

  * [18] C. Rubbia, The Liquid Argon Time Projection Chamber: A New Concept for Neutrino Detectors, CERN-EP-INT-77-08, CERN-EP-77-08 (1977).
URL https://cds.cern.ch/record/117852

  * [19] R. Acciarri, et al., Design and construction of the MicroBooNE detector, JINST 12 (02) (2017) P02017–P02017.  doi:10.1088/1748-0221/12/02/p02017.
URL http://dx.doi.org/10.1088/1748-0221/12/02/P02017

  * [20] C. Anderson, et al., The ArgoNeuT detector in the NuMI low-energy beam line at fermilab, JINST 7 (10) (2012) P10019–P10019.  doi:10.1088/1748-0221/7/10/p10019.
URL https://doi.org/10.1088/1748-0221/7/10/p10019

  * [21] E. Aprile, et al., The XENON100 dark matter experiment, Astroparticle Physics 35 (9) (2012) 573–590.  doi:10.1016/j.astropartphys.2012.01.003.
URL http://dx.doi.org/10.1016/j.astropartphys.2012.01.003

  * [22] D. Akerib, et al., The Large Underground Xenon (LUX) experiment, Nucl. Instrum. Meth. A 704 (2013) 111–126.  doi:10.1016/j.nima.2012.11.135.
URL https://www.sciencedirect.com/science/article/pii/S0168900212014829

  * [23] X. Cao, et al., PandaX: A Liquid Xenon Dark Matter Experiment at CJPL, Sci. China Phys. Mech. Astron. 57 (2014) 1476–1494.  arXiv:1405.2882, doi:10.1007/s11433-014-5521-2.
URL https://doi.org/10.1007/s11433-014-5521-2

  * [24] R. Acciarri, et al., Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE): Conceptual Design Report, Volume 1: The LBNF and DUNE Projects, FERMILAB-DESIGN-2016-01 (2016).  arXiv:1601.05471.
URL https://arxiv.org/abs/1601.05471

  * [25] B. Abi, et al., First results on ProtoDUNE-SP liquid argon time projection chamber performance from a beam test at the CERN Neutrino Platform, JINST 15 (12) (2020) P12004–P12004.  doi:10.1088/1748-0221/15/12/p12004.
URL http://dx.doi.org/10.1088/1748-0221/15/12/P12004

  * [26] A. Gonidec, et al., Ionization chambers with room temperature liquids for calorimetry, CERN-EP-88-36 (1988).
URL https://cds.cern.ch/record/186317

  * [27] J. Engler, et al., A warm-liquid calorimeter for cosmic-ray hadrons, Nucl. Instrum. Meth. A 427 (3) (1999) 528–542.  doi:10.1016/S0168-9002(99)00051-0.
URL https://www.sciencedirect.com/science/article/pii/S0168900299000510

  * [28] M. Farradèche, et al., Free ion yield of trimethyl bismuth used as sensitive medium for high-energy photon detection, JINST 13 (11) (2018) P11004–P11004.  doi:10.1088/1748-0221/13/11/p11004.
URL https://doi.org/10.1088/1748-0221/13/11/p11004

  * [29] H. Faidas, L. Christophorou, D. McCorkle, Drift velocities of excess electrons in 2,2,4,4-tetramethylpentane and tetramethylsilane: A fast drift technique, Chem. Phys. Lett. 163 (6) (1989) 495–498.  doi:https://doi.org/10.1016/0009-2614(89)85175-9.
URL https://www.sciencedirect.com/science/article/pii/0009261489851759

  * [30] S. Wu, et al., A Tetramethylsilane TPC with Cherenkov light readout and 3D reconstruction, Nucl. Instrum. Meth. A 972 (2020) 163904.  doi:10.1016/j.nima.2020.163904.
URL https://www.sciencedirect.com/science/article/pii/S0168900220303971

  * [31] S. Abe, et al., Precision Measurement of Neutrino Oscillation Parameters with KamLAND, Phys. Rev. Lett. 100 (22) (2008).  doi:10.1103/physrevlett.100.221803.
URL http://dx.doi.org/10.1103/PhysRevLett.100.221803

  * [32] C. Andreopoulos, et al., The GENIE Neutrino Monte Carlo Generator, Nucl. Instrum. Meth. A 614 (2010) 87–104.  arXiv:0905.2517, doi:10.1016/j.nima.2009.12.009.
URL https://doi.org/10.1016/j.nima.2009.12.009

  * [33] S. Agostinelli, et al., Geant4—a simulation toolkit, Nucl. Instrum. Meth. A 506 (3) (2003) 250–303.  doi:10.1016/S0168-9002(03)01368-8.
URL https://www.sciencedirect.com/science/article/pii/S0168900203013688

  * [34] R. Holroyd, Electrons in nonpolar liquids, in: Charged Particle and Photon Interactions with Matter, Taylor&Francis Group, 2002, pp. 175–206.  doi:10.1201/9780203913284.
URL https://doi.org/10.1201/9780203913284

  * [35] J. Engler, Liquid ionization chambers at room temperatures, J. Phys. G: Nucl. Part. Phys. 22 (1) (1996) 1–23, cited by: 14.  doi:10.1088/0954-3899/22/1/002.
URL {https://iopscience.iop.org/article/10.1088/0954-3899/22/1/002}

  * [36] Y. Li, et al., Measurement of longitudinal electron diffusion in liquid argon, Nucl. Instrum. Meth. A 816 (2016) 160–170.  doi:10.1016/j.nima.2016.01.094.
URL http://dx.doi.org/10.1016/j.nima.2016.01.094

  * [37] E. Shibamura, et al., Ratio of diffusion coefficient to mobility for electrons in liquid argon, Phys. Rev. A 20 (1979) 2547–2554.  doi:10.1103/PhysRevA.20.2547.
URL https://link.aps.org/doi/10.1103/PhysRevA.20.2547

  * [38] B. Abi, et al., Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume II: DUNE Physics, FERMILAB-PUB-20-025-ND (2020).  arXiv:2002.03005, doi:10.48550/arXiv.2002.03005.
URL https://arxiv.org/abs/2002.03005

  * [39] D. Autiero, et al., Yearly progress report on NP02 ProtoDUNE Dual Phase (2021), Tech. rep., CERN, Geneva (2021).
URL https://cds.cern.ch/record/2759609

  * [40] D. Autiero, et al., Yearly progress report on NP02 (2022), Tech. rep., CERN, Geneva (2022).
URL https://cds.cern.ch/record/2805710

  * [41] M. Lopes, K. Masuda, T. Doke, Estimation of the electron lifetime in tetramethylsilane, Nucl. Instrum. Meth. A 271 (3) (1988) 464–470.  doi:https://doi.org/10.1016/0168-9002(88)90306-3.
URL https://www.sciencedirect.com/science/article/pii/0168900288903063

  * [42] S. Ochsenbein, et al., Purity in room temperature liquid ionization chambers, Nucl. Instrum. Meth. A 273 (2) (1988) 654–656.  doi:https://doi.org/10.1016/0168-9002(88)90073-3.
URL https://www.sciencedirect.com/science/article/pii/0168900288900733

  * [43] P. Vogel, J. Beacom, Angular distribution of neutron inverse beta decay $\bar{\nu}_{e}+p\to e^{+}+n$, Phys. Rev. D 60 (1999) 053003.  doi:10.1103/PhysRevD.60.053003.
URL https://link.aps.org/doi/10.1103/PhysRevD.60.053003

  * [44] M. Schumann, Dual-Phase Liquid Xenon Detectors for Dark Matter Searches, JINST 9 (2014) C08004.  arXiv:1405.7600, doi:10.1088/1748-0221/9/08/C08004.
URL https://doi.org/10.1088/1748-0221/9/08/C08004

  * [45] W. Bollingfehr, et al., Endlagerkonzepte - Bericht zum Arbeitspaket 5, GRS 272 (2011).
  * [46] Y. Schnellbach, et al., Antineutrino detector cconcepts for safeguarding spent nuclear fuel repositories, in: INMM/ESARDA Joint Annual Meeting, 2023.


# Simulating X-ray Reverberation in the UV-Emitting Regions of Active Galactic
Nuclei Accretion Disks with 3D Multi-Frequency Magnetohydrodynamic Simulations

Amy Secunda Department of Astrophysical Sciences, Princeton University, Peyton
Hall, Princeton, NJ 08544, USA  Yan-Fei Jiang (姜燕飞) Center for Computational
Astrophysics, Flatiron Institute, New York, NY 10010, USA  Jenny E. Greene
Department of Astrophysical Sciences, Princeton University, Peyton Hall,
Princeton, NJ 08544, USA

###### Abstract

Active galactic nuclei (AGN) light curves observed with different wavebands
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 multi-frequency radiation magnetohydrodynamic (MHD) simulations
to model the UV-emitting region of an AGN disk, which is unstable to the
magnetorotational instability (MRI) 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 (DRW) model is a good fit to this UV light curve for
timescales $>5$ days. Meanwhile, the injected X-rays have almost no impact on
the power spectrum of the emitted UV light curve. In addition, 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.52$. These results suggest that hard X-rays with scattering
dominated opacity are likely not the main driver of the reverberation signals.

##  1 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 re-emitted 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, 2018,
2021; De Rosa et al., 2015; Edelson et al., 2015, 2017, 2019a; Jiang et al.,
2017; Fausnaugh et al., 2016; Starkey et al., 2017; 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 re-emitted 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.,
2019b; Yu et al., 2022; Hagen & Done, 2023; Cackett et al., 2023; Kara et al.,
2023). Instead 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; Salvesen, 2022; Panagiotou et al., 2022).
The recent development of accurate, high-speed, numerical methods to calculate
radiative transfer coupled to magnetohydrodynamics (MHD) make 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 multi-frequency radiation
MHD code, which instead of integrating radiation over all frequencies,
integrates radiation over multiple frequency groups (Jiang, 2022). This multi-
frequency 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 injected X-rays. Finally, in Section 4 we summarize our results and
present ideas for future work.

##  2 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\text{kev}]$ to model the radiation field emitted by
the disk locally (at UV wavelengths) and the other group covering the
frequency range $(1\text{keV},\infty)$ 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
with mass $M_{\rm BH}=5\times 10^{7}~{}M_{\odot}$ under the local shearing box
approximation with the box centered at $r_{0}=45R_{\rm s}$, where $R_{\rm s}$
is the Schwarzschild radius. The corresponding Keplerian orbital frequency,
$\Omega=4.73\times 10^{-6}$ s-1 and shear parameter, $q=-d\ln\Omega/d\ln
r=3/2$. We choose our simulation domain to be
$(x,y,z)\in(-1.5H_{g},1.5H_{g})\times(-6H_{g},6H_{g})\times(-48H_{g},48H_{g})$
with resolution $N_{x}=48$, $N_{y}=192$ and $N_{z}=1536$. Here $H_{\rm
g}=$1.11\text{\times}{10}^{12}$$ cm is the gas pressure scale height of the
disk, while the total pressure scale height (including radiation pressure)
will be $H_{\rm r}=6.89H_{\rm g}=$7.65\text{\times}{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 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}^{\rm
mid}$, where $E_{\rm r,0}^{\rm 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. (2009a), 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., 2009a; Kozłowski et al., 2010;
MacLeod et al., 2012; Zu et al., 2013). For Run A, we use a damping timescale
$\tau_{\rm damp}=30$ days and a time-sampling of $dt=1$ day, for Run B we use
$\tau_{\rm damp}=0.1$ days and a time-sampling of $dt=0.001$ days. Run A will
have more variability at low frequencies than Run B, while at high frequencies
Run B will have more variability.

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 eight times
the average UV emission. 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.

We determine the timestep-dependent opacity for the UV component in each cell
with a bi-linear 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 to be 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_{\rm 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 $\beta=200$ in the midplane. We also add
a purely vertical magnetic field with $\beta=2\times 10^{4}$ to increase the
Maxwell stress in order to achieve a temporally averaged $\dot{M}/\dot{M}_{\rm
edd}\approx 2\times 10^{-2}$ once the simulation reaches steady state.

Parameter | Unit
---|---
$\Omega_{0}$ |  $4.73\text{\times}{10}^{-6}$ s-1
$\rho_{0}$ |  $10^{-7}$ g cm-3
$T_{0}$ |  $2\text{\times}{10}^{5}$ K
$H_{\rm g}$ |  $1.11\text{\times}{10}^{12}$ cm
$H_{\rm r}$ | 6.89 $H_{\rm g}=$7.65\text{\times}{10}^{12}$$ cm
$F_{\rm max}$ |  $5.80\text{\times}{10}^{12}$ erg s-1 cm-2
$\tau_{0}$ | $6.26\text{\times}{10}^{4}$
$P_{\rm g,0}$ |  $2.77\text{\times}{10}^{6}$ dyn cm-2
$P_{\rm rat}=P_{\rm r,0}/P_{\rm g,0}$ | 4.37
$c_{\rm rat}=c/c_{\rm s}$ | $5.69\text{\times}{10}^{3}$
$r_{0}$ |  $r_{0}=45$ R${}_{\rm s}=$6.68\text{\times}{10}^{14}$$ cm
Table 1: $\Omega_{0}$, $\rho_{0}$, $T_{0}$, and $H_{\rm g}$ are fiducial units
we use as normalization parameters in our simulation, for time, mass,
temperature, and distance, respectively. $H_{\rm r}$ is a fiducial unit for
the total (gas and radiation combined) pressure scale height. $F_{\rm max}$,
$\tau_{0}$, and $P_{\rm g,0}$, are the initial $F_{r,z}(\tau<1)$, midplane
optical depth, and gas pressure, respectively. $P_{\rm rat}$ and $c_{\rm 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.
$r_{0}$ is the fixed location of the center of the shearing box.

Figure 1: Top panel: Schematic of our simulation and space-time diagram of the
gas temperature in the one- and two-frequency simulations in the region left
and right of the grey line, respectively. The red arrows represent the time-
varying X-ray radiation we inject in through the upper and lower boundaries.
Middle panel: Space-time diagram of the gas density in the one- and two-
frequency simulations in the region left and right of the grey line,
respectively. Bottom 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(\tau_{\rm UV}=1)\approx z(\tau_{\rm XR}=1)$ and the dashed white line is
$z(\tau_{a,\rm UV}=1)$. $\tau_{a,\rm XR}<1$ throughout our simulation.

##  3 Results

###  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 (to the left of the grey line) and two-frequency run B (to the right of
the grey line) in Figure 1. In the top panel the red arrows represent the time
varying X-ray irradiation we inject into the top and bottom of our simulation
box in our two-frequency runs. 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(\tau_{\rm UV}=1)\approx z(\tau_{\rm
XR}=1)$). Absorption optical depth for the X-ray component is smaller than 1
throughout the whole disk, but the effective absorption optical depth is $\sim
1.3$. The vertical profiles for two-frequency 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 panel), which is
higher throughout the simulation when we inject X-rays. Unsurprisingly, the
temperature increases the most in the two-frequency run in the photosphere,
where on short timescales the temperature can exceed the midplane temperature,
and inwards an additional $10~{}H_{\rm g}$, where the temperature can nearly
reach $T_{\rm g}\sim 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 grey line. The midplane radiation energy is dominated by the energy
from the low-frequency group. However, at around $\pm 30~{}H_{\rm 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 $\tau_{\rm XR}\approx 240$ or $\pm 6.5$ $H_{\rm g}$ before it drops by
an e-folding.

Overall, the $\alpha$-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 $\tau_{\rm
therm}=1/(\alpha\Omega)\approx 80$ in all runs.

###  3.2 Intrinsic Variability in UV Light Curves

Figure 2: Top panels compare the UV light curves from the one- and two-
frequency runs. Middle and bottom panels compare the UV light curves and
injected X-rays from the two-frequency runs. $\tau_{\rm damp}=30$ days for the
injected X-ray light curve in the left panels and $\tau_{\rm damp}=0.1$ days
for the injected X-ray light curve in the right panels. 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, $\sigma(F)$, and then subtract off the minimum flux for
each light curve. Time zero corresponds to when we start injecting the X-rays.

Figure 3: The top panels show the PSDs of the one-frequency UV, the two-
frequency UV, and injected X-ray light curves in Figure 2. The bottom panels
compare the PSDs of the injected and emitted X-rays for both two-frequency
runs. The left (right) panels show the PSDs from run A (run B), where
$\tau_{\rm damp}=30$ days ($\tau_{\rm damp}=0.1$ days) for the injected X-ray
light curve. We see some evidence that the UV PSDs are damped at $f\lesssim
0.03$ days-1. At $0.03<f<0.2$ days-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$ days-1 the UV PSDs are well fit by
$f^{-2.5}$, shown at the dotted line.

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 clearly has variability over timescales of monthly to sub-daily that is
driven by MRI turbulence and convection. There is evidence that the UV PSD is
damped at the lowest frequencies ($f\lesssim 1/(30~{}\rm{days})$), and then
goes as $f^{-2}$ for $0.03<f<0.2$ days-1, suggesting that a DRW model might be
a good fit for the intrinsic UV light curve at these frequencies. Recent
examinations of long cadence optical AGN light curves have suggested that the
DRW damping timescale is related to the thermal timescale of the disk (Kelly
et al., 2009b; 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 ($\sim 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, $\sim f^{-2.5}$, than the DRW model at
$f\gtrsim 1$/(5 days). In recent high cadence observations of AGN light curves
from Kepler, on timescales around ten 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, it is clear that
our intrinsic UV light curve has variability on timescales of monthly to sub-
daily, due to MRI- and convection-driven turbulence.

###  3.3 Hard X-rays Have Little Impact on UV Emission

We compare the X-ray and UV light curves from the two-frequency runs with
X-ray irradiation in the bottom panels of Figure 2 and the PSDs of these light
curves in the top panels of Figure 3. For run A, the Pearson r-coefficient
between the injected X-ray light curve with damping timescale $\tau_{\rm
damp}=30$ days and the reprocessed UV light curve is $r=0.52$, with a p-value,
$p<10^{-5}$. For run B, where the damping timescale of the injected X-ray
light curve is $\tau_{\rm damp}=0.1$ days, we do not get a significant
correlation between the X-ray and reprocessed UV light curves, $p=0.54$. 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. However, the lack of correlation 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\gtrsim 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 does
have some impact on the re-emitted UV light curves 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, with $r=0.46$ and $r=0.54$, for run A and run B, respectively.
This picture agrees with 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. Furthermore 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).

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 supermassive
black hole (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 gas is fully ionized, which may not be true if the
temperature near the disk surface drops below $T_{\rm g}\approx 10^{4}$ K,
which it might for different initial conditions. A lower ionization fraction
would decrease the scattering opacity and likely lead to more absorption
(García et al., 2013).

Salvesen (2022) found there was a non-zero 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$ days. On the other
hand, we found a very weak correlation for run B, $r=0.048$ and $p=0.017$,
occurs for a lag of 27 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 simulation shows 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 (Morales et al., 2019; Buisson et al.,
2018, e.g.,) or only weak correlations (e.g., Schimoia et al., 2015; Edelson
et al., 2019a; 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).

##  4 Summary and Future Work

Using the multi-frequency radiation MHD code from Jiang (2022) allows us to
simulate the reprocessing of X-ray light curves into UV light curves 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 locally generated by
MRI and convection. This variability can be described by a DRW, with
$\tau_{\rm damp}\approx 30$ days on timescales greater than five 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 will not be 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.52$
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 (e.g., Schimoia et al., 2015; Edelson et al., 2019b;
Cackett et al., 2023).

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; Salvesen, 2022; Panagiotou et al., 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. Doing so in the
simulations presented here, provides detailed information on the nature of
light curve variability produced by disk turbulence, and the ability of X-ray
light curves with different PSDs to transfer energy and momentum to the AGN
disk, affecting the variability of its emission.

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. For example, radiation MHD simulations show that pressure support
from magnetic fields and radiation can lead to puffed up disks (e.g., Jiang et
al., 2016, 2019; Jiang & Blaes, 2020; Hopkins et al., 2023), while recent
observations provide additional evidence for thicker disks than predicted by
the standard thin disk model (Yao et al., 2022; Secunda et al., 2023). In
addition, both micro-lensing 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 whether lower energy soft X-ray or far-UV irradiation, which will have
larger absorption opacities compared to hard X-ray irradiation, are more
likely to be absorbed and reprocessed by the disk producing stronger
correlations between injected and reprocessed light curves. Future work should
also look at 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.

A.S. is supported by a fellowship from the NSF Graduate Research Fellowship
Program under Grant No. DGE-1656466. JEG is supported in part by NSF grants
AST1007052 and AST1007094. The Center for Computational Astrophysics at the
Flatiron Institute is supported by the Simons Foundation.

## References

  * Arévalo et al. (2009) Arévalo, P., Uttley, P., Lira, P., et al. 2009, MNRAS, 397, 2004, doi: 10.1111/j.1365-2966.2009.15110.x
  * Bai & Stone (2013) Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76, doi: 10.1088/0004-637X/769/1/76
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Bentz & Katz (2015) Bentz, M. C., & Katz, S. 2015, Publications of the Astronomical Society of the Pacific, 127, 67, doi: 10.1086/679601
  * Blandford & McKee (1982) Blandford, R. D., & McKee, C. F. 1982, ApJ, 255, 419, doi: 10.1086/159843
  * Buisson et al. (2018) Buisson, D. J. K., Lohfink, A. M., Alston, W. N., et al. 2018, MNRAS, 475, 2306, doi: 10.1093/mnras/sty008
  * Burke et al. (2021) Burke, C. J., Shen, Y., Blaes, O., et al. 2021, Science, 373, 789, doi: 10.1126/science.abg9933
  * Cackett et al. (2021) Cackett, E. M., Bentz, M. C., & Kara, E. 2021, iScience, 24, 102557, doi: 10.1016/j.isci.2021.102557
  * Cackett et al. (2018) Cackett, E. M., Chiang, C.-Y., McHardy, I., et al. 2018, ApJ, 857, 53, doi: 10.3847/1538-4357/aab4f7
  * Cackett et al. (2007) Cackett, E. M., Horne, K., & Winkler, H. 2007, MNRAS, 380, 669, doi: 10.1111/j.1365-2966.2007.12098.x
  * Cackett et al. (2023) Cackett, E. M., Gelbord, J., Barth, A. J., et al. 2023, arXiv e-prints, arXiv:2306.17663, doi: 10.48550/arXiv.2306.17663
  * De Rosa et al. (2015) De Rosa, G., Peterson, B. M., Ely, J., et al. 2015, ApJ, 806, 128, doi: 10.1088/0004-637X/806/1/128
  * Edelson et al. (2015) Edelson, R., Gelbord, J. M., Horne, K., et al. 2015, ApJ, 806, 129, doi: 10.1088/0004-637X/806/1/129
  * Edelson et al. (2017) Edelson, R., Gelbord, J., Cackett, E., et al. 2017, ApJ, 840, 41, doi: 10.3847/1538-4357/aa6890
  * Edelson et al. (2019a) —. 2019a, ApJ, 870, 123, doi: 10.3847/1538-4357/aaf3b4
  * Edelson et al. (2019b) —. 2019b, ApJ, 870, 123, doi: 10.3847/1538-4357/aaf3b4
  * Fausnaugh et al. (2016) Fausnaugh, M. M., Denney, K. D., Barth, A. J., et al. 2016, ApJ, 821, 56, doi: 10.3847/0004-637X/821/1/56
  * García et al. (2013) García, J., Dauser, T., Reynolds, C. S., et al. 2013, ApJ, 768, 146, doi: 10.1088/0004-637X/768/2/146
  * Gardner & Done (2017) Gardner, E., & Done, C. 2017, MNRAS, 470, 3591, doi: 10.1093/mnras/stx946
  * Grier et al. (2017) Grier, C. J., Trump, J. R., Shen, Y., et al. 2017, The Astrophysical Journal, 851, 21, doi: 10.3847/1538-4357/aa98dc
  * Guo et al. (2022) Guo, H., Barth, A. J., & Wang, S. 2022, arXiv e-prints, arXiv:2207.06432.  https://arxiv.org/abs/2207.06432
  * Hagen & Done (2023) Hagen, S., & Done, C. 2023, MNRAS, 521, 251, doi: 10.1093/mnras/stad504
  * Hall et al. (2018) Hall, P. B., Sarrouh, G. T., & Horne, K. 2018, ApJ, 854, 93, doi: 10.3847/1538-4357/aaa768
  * Homayouni et al. (2022) Homayouni, Y., Sturm, M. R., Trump, J. R., et al. 2022, ApJ, 926, 225, doi: 10.3847/1538-4357/ac478b
  * Hopkins et al. (2023) Hopkins, P. F., Squire, J., Su, K.-Y., et al. 2023, arXiv e-prints, arXiv:2310.04506, doi: 10.48550/arXiv.2310.04506
  * Iglesias & Rogers (1996) Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943, doi: 10.1086/177381
  * Jiang (2021) Jiang, Y.-F. 2021, ApJS, 253, 49, doi: 10.3847/1538-4365/abe303
  * Jiang (2022) —. 2022, ApJS, 263, 4, doi: 10.3847/1538-4365/ac9231
  * Jiang & Blaes (2020) Jiang, Y.-F., & Blaes, O. 2020, ApJ, 900, 25, doi: 10.3847/1538-4357/aba4b7
  * Jiang et al. (2019) Jiang, Y.-F., Blaes, O., Stone, J. M., & Davis, S. W. 2019, ApJ, 885, 144, doi: 10.3847/1538-4357/ab4a00
  * Jiang et al. (2016) Jiang, Y.-F., Davis, S. W., & Stone, J. M. 2016, ApJ, 827, 10, doi: 10.3847/0004-637X/827/1/10
  * Jiang et al. (2013) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2013, ApJ, 778, 65, doi: 10.1088/0004-637X/778/1/65
  * Jiang et al. (2017) Jiang, Y.-F., Green, P. J., Greene, J. E., et al. 2017, ApJ, 836, 186, doi: 10.3847/1538-4357/aa5b91
  * Kammoun et al. (2021) Kammoun, E. S., Papadakis, I. E., & Dovčiak, M. 2021, MNRAS, 503, 4163, doi: 10.1093/mnras/stab725
  * Kara et al. (2021) Kara, E., Mehdipour, M., Kriss, G. A., et al. 2021, ApJ, 922, 151, doi: 10.3847/1538-4357/ac2159
  * Kara et al. (2023) Kara, E., Barth, A. J., Cackett, E. M., et al. 2023, ApJ, 947, 62, doi: 10.3847/1538-4357/acbcd3
  * Kaspi et al. (2000) Kaspi, S., Smith, P. S., Netzer, H., et al. 2000, The Astrophysical Journal, 533, 631, doi: 10.1086/308704
  * Kelly et al. (2009a) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009a, ApJ, 698, 895, doi: 10.1088/0004-637X/698/1/895
  * Kelly et al. (2009b) —. 2009b, ApJ, 698, 895, doi: 10.1088/0004-637X/698/1/895
  * Kozłowski et al. (2010) Kozłowski, S., Kochanek, C. S., Udalski, A., et al. 2010, ApJ, 708, 927, doi: 10.1088/0004-637X/708/2/927
  * MacLeod et al. (2012) MacLeod, C. L., Ivezić, Ž., Sesar, B., et al. 2012, ApJ, 753, 106, doi: 10.1088/0004-637X/753/2/106
  * Morales et al. (2019) Morales, A. M., Miller, J. M., Cackett, E. M., Reynolds, M. T., & Zoghbi, A. 2019, ApJ, 870, 54, doi: 10.3847/1538-4357/aaeff9
  * Morgan et al. (2010) Morgan, C. W., Kochanek, C. S., Morgan, N. D., & Falco, E. E. 2010, ApJ, 712, 1129, doi: 10.1088/0004-637X/712/2/1129
  * Panagiotou et al. (2022) Panagiotou, C., Papadakis, I., Kara, E., Kammoun, E., & Dovčiak, M. 2022, ApJ, 935, 93, doi: 10.3847/1538-4357/ac7e4d
  * Peterson (2014) Peterson, B. M. 2014, Space Sci. Rev., 183, 253, doi: 10.1007/s11214-013-9987-4
  * Peterson et al. (2004) Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682, doi: 10.1086/423269
  * Salvesen (2022) Salvesen, G. 2022, ApJ, 940, L22, doi: 10.3847/2041-8213/ac9cdd
  * Schimoia et al. (2015) Schimoia, J. S., Storchi-Bergmann, T., Grupe, D., et al. 2015, ApJ, 800, 63, doi: 10.1088/0004-637X/800/1/63
  * Secunda et al. (2023) Secunda, A., Greene, J. E., Jiang, Y.-F., Yao, P. Z., & Zoghbi, A. 2023, arXiv e-prints, arXiv:2306.05455, doi: 10.48550/arXiv.2306.05455
  * Sergeev et al. (2005) Sergeev, S. G., Doroshenko, V. T., Golubinskiy, Y. V., Merkulova, N. I., & Sergeeva, E. A. 2005, ApJ, 622, 129, doi: 10.1086/427820
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Smith et al. (2018) Smith, K. L., Mushotzky, R. F., Boyd, P. T., et al. 2018, ApJ, 857, 141, doi: 10.3847/1538-4357/aab88d
  * Starkey et al. (2017) Starkey, D., Horne, K., Fausnaugh, M. M., et al. 2017, ApJ, 835, 65, doi: 10.3847/1538-4357/835/1/65
  * Stone et al. (2022) Stone, Z., Shen, Y., Burke, C. J., et al. 2022, MNRAS, 514, 164, doi: 10.1093/mnras/stac1259
  * Sun et al. (2020) Sun, M., Xue, Y., Guo, H., et al. 2020, ApJ, 902, 7, doi: 10.3847/1538-4357/abb1c4
  * Tachibana et al. (2020) Tachibana, Y., Graham, M. J., Kawai, N., et al. 2020, ApJ, 903, 54, doi: 10.3847/1538-4357/abb9a9
  * Uttley et al. (2003) Uttley, P., Edelson, R., McHardy, I. M., Peterson, B. M., & Markowitz, A. 2003, ApJ, 584, L53, doi: 10.1086/373887
  * Yao et al. (2022) Yao, P. Z., Secunda, A., Jiang, Y.-F., Greene, J. E., & Villar, A. 2022, arXiv e-prints, arXiv:2210.13489.  https://arxiv.org/abs/2210.13489
  * Yu et al. (2022) Yu, W., Richards, G. T., Vogeley, M. S., Moreno, J., & Graham, M. J. 2022, ApJ, 936, 132, doi: 10.3847/1538-4357/ac8351
  * Zu et al. (2013) Zu, Y., Kochanek, C. S., Kozłowski, S., & Udalski, A. 2013, ApJ, 765, 106, doi: 10.1088/0004-637X/765/2/106