# Reverberation Mapping of IC 4329A

Misty C. Bentz Department of Physics and Astronomy, Georgia State University,
Atlanta, GA 30303, USA<EMAIL_ADDRESS>Christopher A. Onken Research
School of Astronomy and Astrophysics, Australian National University,
Canberra, ACT 2611, Australia  Rachel Street LCOGT, 6740 Cortona Drive, Suite
102, Goleta, CA 93117, USA  Monica Valluri Department of Astronomy, University
of Michigan, Ann Arbor, MI, 48109, USA

###### Abstract

We present the results of a new reverberation mapping campaign for the broad-
lined active galactic nucleus (AGN) in the edge-on spiral IC 4329A. Monitoring
of the optical continuum with $V-$band photometry and broad emission-line flux
variability with moderate-resolution spectroscopy allowed emission-line light
curves to be measured for H$\beta$, H$\gamma$, and He II $\lambda 4686$. We
find a time delay of $16.3^{+2.6}_{-2.3}$ days for H$\beta$, a similar time
delay of $16.0^{+4.8}_{-2.6}$ days for H$\gamma$, and an unresolved time delay
of $-0.6^{+3.9}_{-3.9}$ days for He II. The time delay for H$\beta$ is
consistent with the predicted value from the relationship between AGN
luminosity and broad line region radius, after correction for the $\sim 2.4$
mag of intrinsic extinction at 5100 Å. Combining the measured time delay for
H$\beta$ with the broad emission line width and an adopted value of $\langle
f\rangle=4.8$, we find a central supermassive black hole mass of $M_{\rm
BH}=6.8^{+1.2}_{-1.1}\times 10^{7}$$M_{\odot}$. Velocity-resolved time delays
were measured across the broad H$\beta$ emission-line profile and may be
consistent with an “M”-like shape. Modeling of the full reverberation response
of H$\beta$ was able to provide only modest constraints on some parameters,
but does exhibit agreement with the black hole mass and average time delay.
The models also suggest that the AGN structure is misaligned by a large amount
from the edge-on galaxy disk. This is consistent with expectations from the
unified model of AGNs, in which broad emission lines are expected to be
visible only for AGNs that are viewed at relatively face-on inclinations.

Seyfert galaxies (1447) — Supermassive black holes (1663) — Reverberation
mapping(2019)

††journal: ApJ††software: LCO Reduction Pipeline (McCully et al., 2018),
JAVELIN (Zu et al., 2011), Galfit (Peng et al., 2002, 2010), CARAMEL (Pancoast
et al., 2014)

##  1 Introduction

IC 4329A is a striking example of a broad-lined active galactic nucleus (AGN)
in an edge-on galaxy. Within the unified model for AGNs (Antonucci, 1993),
this would suggest that the rotation axis of the central supermassive black
hole and the rotation axis of the galaxy disk are misaligned by a considerable
angle. The misalignment of the galaxy and black hole are further supported by
a $\sim 50^{\circ}$ difference between the position angle of the disk of IC
4329A and the position angle of extended radio emission from its nucleus
(Schmitt et al., 1997; Nagar & Wilson, 1999). AGN fueling events can torque
the spin axis of the black hole (e.g., King & Pringle 2007) while black hole
mergers can flip the spin axis (e.g., Lousto et al. 2016). With its membership
within a loose group of galaxies inside the cluster Abell 3574 (Disney, 1973;
Wilson & Penston, 1979; Abell et al., 1989) and close proximity to the giant
lenticular IC 4329, it is possible that either of these scenarios may have
occurred in IC 4329A. Furthermore, there is observational evidence for ongoing
interactions between IC 4329 and IC 4329A (Read & Pietsch, 1998), and IC 4329
is a shell galaxy (Malin & Carter, 1983), which may be evidence of a past
merger.

IC 4329A hosts a highly reddened AGN that is nevertheless apparently bright
across the electromagnetic spectrum, which led Wilson & Penston (1979) to
suggest that it is the “nearest quasar”. Prior studies have concluded that the
reddening arises from the AGN structure itself as well as the dust lane of the
galaxy disk (Wolstencroft et al., 1995; Mehdipour & Costantini, 2018). IC
4329A thus provides a rare opportunity to carry out detailed studies of the
absorbing material around an accreting supermassive black hole, from the far-
infrared through the X-rays (e.g., Mehdipour & Costantini 2018).

Figure 1: Left: The central $10\arcmin\times 10\arcmin$ of the reference image
of IC 4329A, centered on the AGN. The image is oriented with North up and East
to the left, and the giant lenticular galaxy IC 4329 is visible to the west of
IC 4329A. Right: Hubble Space Telescope image of the central $20\arcsec\times
20\arcsec$ in IC 4329A, highlighting the bright central AGN peeking over the
dust lane in the disk of the galaxy.

Currently lacking among the studies of IC 4329A, however, is an accurate
constraint on the black hole mass. Reverberation mapping (Cackett et al.,
2021), which uses spectrophotometric monitoring of an AGN to track changes in
the continuum flux and the echoes of those changes in the broad emission
lines, is often applied to bright and nearby broad-lined AGNs and has resulted
in a sample of $\sim 70$ black hole mass measurements to date in galaxies with
$z\lesssim 0.4$ (Bentz & Katz, 2015).

Strong flux variability in the nucleus of IC 4329A has been known for decades
and an early reverberation mapping attempt was made by Winge et al. (1996).
However, that early program was carried out before the typical size of the
broad line region (BLR) was well understood, and when estimates were about a
factor of 10 too large. With an average (median) sampling rate of 8 (4) days,
and utilizing the CTIO 2D-Frutti camera in the spectrograph — a photon-
counting detector with significantly poorer performance than a CCD — the light
curves were noisy, undersampled, and provided unreliable time delay
measurements (Peterson et al., 2004).

We therefore undertook a new reverberation mapping campaign in the first half
of 2022 with the aim of constraining the black hole mass of IC 4329A, and we
describe the results of that program in this manuscript.

##  2 Observations

IC 4329A is an edge-on spiral galaxy located in the direction of Centaurus at
RA=13h49m19.3s, Dec=$-30^{\circ}18\arcmin 34\arcsec$ (J2000.0), and
$z=0.01605$. Based on fundamental plane studies (Hudson et al., 2001; Bernardi
et al., 2002), the distance to the giant lenticular IC 4329 is $D=59\pm 9$ Mpc
(Tully et al., 2009), which is consistent with the distance of $D=69\pm 18$
Mpc derived from observations of a Type 1a supernova (Stahl et al., 2021). IC
4329A is expected to have a similar distance as IC 4329 given the signatures
of ongoing interactions between them.

###  2.1 Imaging

Photometric monitoring of IC 4329A was carried out with the Las Cumbres
Observatory (LCO; Brown et al. 2013) 1.0 m telescope network under program
NSF-2022A-012 (PI: Bentz)111LCO telescope time was granted by NOIRLab through
the Mid-Scale Innovations Program (MSIP). MSIP is funded by NSF.. Observations
began on 2022 February 03 (UT) and continued throughout the 2022A semester,
ending on 2022 July 31. A 60 s $V-$band image was scheduled to occur every
$\sim 8$ hours using the Sinistro imaging cameras, providing a $26\farcm
5\times 26\farcm 5$ field of view and an angular scale of 0$\farcs$389 per
pixel.

Over the course of the semester, a total of 374 images were acquired from the
LCO facilities at three sites – 164 from Cerro Tololo Interamerican
Observatory (CTIO), 138 from South African Astronomical Observatory (SAAO),
and 72 from Siding Spring Observatory (SSO) – with a median (average) time
sampling of 0.35 (0.47) days. All observations were reduced by the LCO
pipeline (McCully et al., 2018), which applies typical CCD corrections using
bias, flat, and dark frames, and the fully reduced images were downloaded from
the LCO archive222https://archive.lco.global.

The $V-$band covers a relatively emission-line free region of the spectral
energy distribution for nearby AGNs, and is therefore an excellent probe of
the continuum variability. However, nearby AGNs live in optically bright and
spatially resolved nearby galaxies. We therefore relied on image subtraction
methods (Alard & Lupton, 1998; Alard, 2000) to isolate and measure the
variable flux in the nucleus of IC 4329A. Image subtraction has the advantage
over aperture photometry in that nonvariable flux sources, such as the host
galaxy, are removed in the processing. The removal of the host galaxy results
in a reduction of the noise caused by variable seeing affecting the amount of
light within a fixed aperture, and it also removes the damping effect on the
intrinsic variability that is caused by a large constant-flux offset. Because
of the advantages, image subtraction methods are commonly applied in studies
of transients and other time domain phenomena including microlensing (Udalski
et al., 2008, 2015), supernovae (Riess et al., 2001; Miknaitis et al., 2007;
Melinder et al., 2008), tidal disruption events (Holoien et al., 2016; Brown
et al., 2018), as well as AGN variability (Fausnaugh et al., 2016; Grier et
al., 2017; Bentz et al., 2021a).

We first registered all images to a common reference grid with the algorithm
of Siverd et al. (2012). Using the algorithms of Alard & Lupton (1998) and
Alard (2000), we then built a reference image (Figure 1) from the subset of
frames collected at CTIO with the lowest backgrounds and best seeing. For each
observation collected throughout the semester, we then convolved the reference
frame to match the image, subtracted the convolved reference from the observed
image, and measured the residual counts at the position of the AGN using an
aperture with a radius of 9 pixels.

To convert the residual counts to calibrated fluxes, we carried out a 2D
surface brightness decomposition of the reference frame. Using Galfit (Peng et
al., 2002, 2010), we modeled the separate photometric components of IC 4329A
and the nearby lenticular IC 4329. A model of the point spread function (PSF)
of the image was generated by fitting four Gaussian profiles to a bright and
isolated field star. For each of the Gaussian profiles, we allowed the first
mode of the Fourier perturbation series to vary in an effort to better account
for slight asymmetries in the PSF shape. The resultant PSF model image was
then used to model the AGN as well as several field stars with catalogued
$V-$band magnitudes. An exponential profile with a radial truncation was used
to model the edge-on galaxy disk and its extinguishing dust lane, while Sérsic
profiles modeled the bulge and companion galaxy. Finally, the background sky
was modeled with a tilted plane to allow for gradients in both the x and y
directions. The final zeropoint of the flux calibration was set by minimizing
the difference between the fitted magnitudes of the fields stars and their
reported magnitudes in the AAVSO Photometric All Sky Survey catalog (Henden &
Munari, 2014). Once the models were optimized, the final modeled AGN magnitude
was combined with the residual counts to determine the calibrated fluxes for
each $V-$band measurement throughout the campaign.

As we found in our previous monitoring study of NGC 3783 (Bentz et al.,
2021a), the images acquired from the three Southern LCO 1 m sites have slight
photometric offsets between them. We identified observations that were
collected close in time ($<0.5$ day) from CTIO and SSO, as well as from CTIO
and SAAO. We then determined and applied the slight linear scaling needed to
bring the SAAO and SSO measurements into agreement with the CTIO measurements,
as the CTIO data set was used to build the reference frame for image
subtraction and to determine the AGN mean magnitude.

Figure 2: Top: Final calibrated and binned $V-$band light curve of IC 4329A.
Bottom: Magnitude residuals for a non-varying $V=15.3$ mag field star.

Previous work has found that the flux uncertainties associated with image
subtraction methods are often underestimated (Zebrun et al., 2001; Hartman et
al., 2005). We therefore followed the basic procedure outlined by Hartman et
al. (2004) and examined the residual counts left behind by non-varying field
stars, identifying the factor by which the uncertainties needed to be expanded
to account for the scatter in residual flux from stars of a similar magnitude
as the AGN. The uncertainties on the AGN fluxes were then inflated by this
scale factor.

Finally, over the course of $\sim 2$ months in the middle of the campaign, the
fluxes measured for the AGN varied strongly and extremely rapidly from
observation to observation. Inspection of the images acquired during these
dates found that some of this rapid variability coincided with bright and
uneven backgrounds in the images that could not be approximated with a simple
gradient. Many of these dates also coincide with the Full Moon, and the uneven
backgrounds appear to be associated with the reflection of moonlight within
the domes. Rather than discard all of these frames and accept gaps of several
days in the middle of the light curve, we instead discarded only the most
egregious frames and binned all the remaining measurements with a bin size of
1.5 days. This bin size was determined through testing to provide a reasonable
compromise between retaining temporal sampling and decreasing wild flux
variations that cause difficulties in the time delay analysis.

In Figure 2 we display the final calibrated and binned $V-$band light curve of
IC 4329A, along with a plot of the residuals from a field star of similar
brightness for comparison.

###  2.2 Spectroscopy

Figure 3: Mean (top) and rms (bottom) of the IC 4329A spectra collected
throughout the monitoring program. Integration limits for the broad emission
lines are displayed as the horizontal bars, and the adopted continuum windows
are displayed as the blue horizontal bars. The gray (black) lines show the
spectra before (after) narrow-line subtraction.

Spectroscopic monitoring with the FLOYDS spectrographs on the 2.0 m Faulkes
Telescopes North and South was carried out under programs NSF-2022A-012 and
ANU-2022A-002. The FLOYDS spectrographs are cross-dispersed instruments
providing wavelength coverage of $540-1000$ nm with a dispersion of 3.51 Å/pix
in the first order, and wavelength coverage of $320-570$ nm with a dispersion
of 1.74 Å/pix in the second order. Observations began on 2022 March 07 and
were scheduled every $\sim 24$ hours through 2022 July 31. Over the course of
the semester, 45 spectra (23 from FTN and 22 from FTS) were acquired with a
median (average) time sampling of 2.0 (3.2) days.

Each observation consisted of a 1200 s spectrum of IC 4329A through the 6.0″
slit, with an arc lamp and flat field taken at the same position. On several
visits we also acquired a 60 s spectrum of the nearby A0 spectrophotometric
standard star CD-32 9927 (Hamuy et al., 1992, 1994) through the 6.0″ slit.

The LCO pipeline splits the two spectral orders into separate images, after
which it rectifies each order and applies typical CCD reductions including
biases, flats, and darks as well as an initial wavelength and flux calibration
based on cataloged wavelength and sensitivity functions. We began our custom
processing with the 2D flux- and wavelength-calibrated frames produced by the
pipeline, focusing on the 2nd order (blue) spectra with their higher spectral
resolution and coverage of the H$\beta$ and [O III] emission. Working in
IRAF333IRAF is distributed by the National Optical Astronomy Observatories,
which are operated by the Association of Universities for Research in
Astronomy, Inc., under cooperative agreement with the National Science
Foundation. we detected and cleaned the science frames of cosmic rays and then
extracted 1D spectra of IC 4329A, the standard star, and the arc lamps with a
10 pixel extraction width (corresponding to an angular width of 3$\farcs$37).
We then used the arc lamp spectra to improve the wavelength calibration for
each observation, and finally adopted and applied a common dispersion of 1.6
Å/pix.

For the visits that included an observation of the standard star, we fit a low
order polynomial to the standard star spectrum and compared this to the
tabulated spectrum of Hamuy et al. (1992, 1994). We then determined the
scaling needed to match the tabulated standard spectrum, and applied this
scaling to the spectrum of IC 4329A that was collected in the same visit.
Given the lack of standard star observations on several nights and the
nonphotometric conditions under which many spectra were acquired, we were
unable to improve the flux calibration of all the observations of IC 4329A in
this way. We instead focused on determining the absolute calibration of the [O
III] $\lambda 5007$ Å emission line from photometric nights that included
observations of both IC 4329A and CD-32 9927, finding $F(5007)=(2.34\pm
0.13)\times 10^{-13}$ erg s-1 cm-2. This value agrees well with the flux
reported by Morris & Ward (1988) of $F(5007)=2.4\times 10^{-13}$ erg s-1 cm-2
and is slightly higher, but within 2$\sigma$, of the value reported by Winge
et al. (1996) of $F(5007)=(1.87\pm 0.22)\times 10^{-13}$ erg s-1 cm-2.

Finally, we applied the spectral scaling method of van Groningen & Wanders
(1992) to all of the individual spectra of IC 4329A, focusing on the [O III]
doublet. The [O III] emission lines are not variable on the timescales probed
within a single observing semester (Peterson et al., 2013) and may therefore
be used as an internal calibration “lamp”. The algorithm applies small shifts
and smoothing to minimize the differences between an observed spectrum and a
reference spectrum (created by averaging together several of the best-quality
spectra), thus correcting for residual wavelength shifts, offsets in the flux
calibration relative to the value of $F(5007)$ that we determined above, and
differences in the spectral resolution arising from variable seeing throughout
the campaign.

In Figure 3 we display the mean and root-mean-square (rms) of the blue spectra
of IC 4329A collected throughout the monitoring campaign. Strong variability
is evident in the rms of the broad Balmer lines H$\beta$ and H$\gamma$ as well
as the He II $\lambda 4686$ Å line.

Table 1: Continuum and Emission-Line Light Curves

Continuum |  | Emission Lines
---|---|---
| $V$ |  |  | H$\beta$ | H$\gamma$ | He II
HJD | $F_{\lambda}$ | $\sigma_{F}$ |  | HJD | $F$ | $\sigma_{F}$ | $F$ | $\sigma_{F}$ | $F$ | $\sigma_{F}$
9615.1116 | 2.666 | 0.046 |  | 9646.0925 | 3.124 | 0.031 | 0.700 | 0.025 | 0.746 | 0.031
9616.9787 | 2.651 | 0.032 |  | 9647.0745 | 3.075 | 0.029 | 0.680 | 0.023 | 0.747 | 0.029
9618.7941 | 2.636 | 0.031 |  | 9649.1018 | 3.020 | 0.033 | 0.639 | 0.027 | 0.584 | 0.034
9620.5635 | 2.650 | 0.034 |  | 9651.0854 | 3.155 | 0.028 | 0.711 | 0.022 | 0.715 | 0.028
9622.5523 | 2.660 | 0.038 |  | 9657.0617 | 3.022 | 0.032 | 0.685 | 0.030 | 0.521 | 0.036

Note. — Heliocentric Julian dates are provided as HJD$-2450000$ (days).
$V-$band flux densities have units of $10^{-15}$ erg s-1 cm-2 Å-1 while
emission-line fluxes have units of $10^{-13}$ erg s-1 cm-2. Table 1 is
published in its entirety in the machine-readable format. A portion is shown
here for guidance regarding its form and content.

Table 2: Light Curve Statistics

Time Series | $N$ |  $\langle\Delta T\rangle$ (days) |  $\Delta T_{\rm med}$ (days) | $\langle F\rangle$ | $\langle\sigma_{F}/F\rangle$ | $F_{\rm var}$ | $R_{\rm max}$
---|---|---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8)
V | 91 | $2.0\pm 0.4$ | 1.9 | $2.49\pm 0.18$ | 0.014 | 0.07 | $1.332\pm 0.027$
H$\beta$ | 45 | $3.2\pm 3.9$ | 2.0 | $2.57\pm 0.25$ | 0.011 | 0.10 | $1.436\pm 0.028$
H$\gamma$ | 45 | $3.2\pm 3.9$ | 2.0 | $0.52\pm 0.08$ | 0.049 | 0.15 | $1.859\pm 0.126$
He II | 45 | $3.2\pm 3.9$ | 2.0 | $0.44\pm 0.13$ | 0.079 | 0.28 | $6.553\pm 2.661$

Note. — $V-$band flux densities have units of $10^{-15}$ erg s-1 cm-2 Å-1
while emission-line fluxes have units of $10^{-13}$ erg s-1 cm-2.

##  3 Analysis

###  3.1 Emission-Line Light Curves and Time Delays

With our final scaled spectra, we were able to measure the integrated fluxes
of the broad emission lines and create their light curves. We began by setting
a local linear continuum based on line-free regions on either side of a broad
emission line of interest, and then integrating all of the flux above the
continuum within a specified window. Thus the emission-line fluxes include the
narrow line components as constant flux offsets. While simplistic, this method
avoids inserting additional uncertainties into the measurements and has been
shown to work well when applied to high signal-to-noise spectra with
relatively unblended emission lines (e.g., Peterson et al. 2004; Bentz et al.
2009b).

Figure 4: Left: $V-$band and emission-line light curves for IC 4329A. Flux
densities are in units of $10^{-15}$ erg s-1 cm-2 Å-1, and integrated fluxes
are in units of $10^{-13}$ erg s-1 cm-2. Right: Cross-correlation functions
relative to the $V-$band light curve (for the $V-$band light curve, this is
the autocorrelation function).

We provide the light curves in Table 1 and plot them in Figure 4, while we
tabulate several useful statistics for each light curve in Table 2. For each
spectral feature listed in column (1), we give the number of measurements in
column (2), and the average and median temporal sampling, respectively, in
columns (3) and (4). Column (5) lists the mean flux and standard deviation,
while column (6) lists the mean fractional error. Column (7) gives the noise-
corrected fractional variation, which is computed as

$F_{\rm var}=\frac{\sqrt{\sigma^{2}-\delta^{2}}}{\langle F\rangle}$ (1)

where $\sigma^{2}$ is the variance of the fluxes, $\delta^{2}$ is their mean-
square uncertainty, and $\langle F\rangle$ is the mean flux. Finally, column
(8) lists the ratio of the maximum to the minimum flux in the light curve.

To measure the time delays between the emission-line light curves and the
$V-$band light curve, we began by employing the Interpolated Cross Correlation
Function (ICCF) method of Gaskell & Sparke (1986); Gaskell & Peterson (1987)
with the modifications of White & Peterson (1994). The ICCF method determines
the cross-correlation function twice, first when the continuum light curve is
interpolated and then when the emission-line light curve is interpolated, and
then averages the two together. These averaged CCFs are displayed in the
panels on the right side of Figure 4.

From the CCF, the time delay may be reported as $\tau_{\rm peak}$, the time
lag at which the peak ($r_{\rm max}$) of the CCF occurs, or $\tau_{\rm cent}$,
the centroid of the CCF above some value (usually $0.8r_{\rm max}$). We employ
the flux randomization/random subset sampling (FR/RSS) method (Peterson et
al., 1998, 2004) to quantify the uncertainties on the time delays. In short,
the flux randomization accounts for the effects of the measurement
uncertainties on the derived time delay, and the random subset sampling
accounts for the inclusion or exclusion of any specific data point in the
analysis. The FR/RSS method is run in a Monte Carlo fashion, carrying out a
large number of realizations to build up a distribution of CCF measurements.
We ran 1000 realizations of the FR/RSS method and adopted the median and 68%
confidence interval of the distributions to define the time delays, $\tau_{\rm
cent}$ and $\tau_{\rm peak}$, and their uncertainties, which are tabulated in
Table 3.

Table 3: Emission-Line Time Lags and Widths

|  |  |  | mean | rms
---|---|---|---|---|---
Line | $\tau_{\rm cent}$ | $\tau_{\rm peak}$ | $\tau_{\rm jav}$ | FWHM | $\sigma_{\rm line}$ | FWHM | $\sigma_{\rm line}$
| (days) | (days) | (days) | (km s-1) | (km s-1) | (km s-1) | (km s-1)
H$\beta$ | $16.33^{+2.59}_{-2.28}$ | $16.00^{+0.50}_{-2.75}$ | $16.09^{+0.54}_{-1.26}$ | $6944\pm 51$ | $2247\pm 8$ | $4789\pm 869$ | $2112\pm 93$
H$\gamma$ | $16.04^{+4.83}_{-2.60}$ | $13.25^{+4.25}_{-2.75}$ | $15.48^{+0.96}_{-2.23}$ | $4176\pm 162$ | $1703\pm 61$ | $4607\pm 948$ | $1908\pm 149$
He II | $-0.59^{+3.87}_{-3.85}$ | $0.25^{+4.75}_{-4.75}$ | $-1.37^{+2.04}_{-2.37}$ |  $\cdots$ |  $\cdots$ | $8099\pm 1574$ | $3334\pm 134$

Note. — Reported time lags are in the observer’s frame, while line widths are
in the rest frame.

We also investigated the time delays with JAVELIN (Zu et al., 2011), which
fits a damped random walk model to the continuum light curve and then
determines the best parameters for a top-hat function with which to delay and
smooth the continuum light curve to match the emission-line light curve.
JAVELIN is capable of fitting multiple emission lines at the same time, and it
assesses the uncertainties on the reported time delay through a Bayesian
Markov Chain Monte Carlo method. We fit all of the emission lines
simultaneously, as well as individually, and found that both methods recovered
time delays that agree fairly well. We list the time delays, denoted as
$\tau_{\rm jav}$, from fitting all three emission lines simultaneously in
Table 3.

###  3.2 Line Width Measurements

The widths and shapes of the broad emission lines provide information on the
kinematics of the line-emitting gas, and are generally reported as the
velocity dispersion of the line profile, $\sigma_{\rm line}$, and the full-
width at half maximum, FWHM. Previous reverberation studies have often
reported different broad-line widths when measured in the mean vs. the rms
spectrum (for a compilation, see Bentz & Katz 2015), with the rms spectrum
crucially providing a measure of the kinematics of the reverberating gas.
Furthermore, the rms spectrum suppresses the non-variable narrow-line
contribution.

While the rms spectrum suppresses the narrow line contribution, the narrow
lines must be fitted and removed from the mean spectrum before measuring the
broad line widths. We used the [O III] $\lambda 5007$ line as a template that
was shifted and scaled to match the other narrow lines in the relevant
portions of the spectrum — [O III] $\lambda 4959$ and $\lambda 4363$, and the
narrow components of H$\beta$, He II, and H$\gamma$ — and then subtracted. In
Figure 3 we display the mean and rms spectra before and after narrow-line
subtraction in gray and black, respectively. As expected, the rms profiles of
the broad emission lines are relatively unaffected by the narrow-line
subtraction.

To determine the uncertainties on the line widths, we employed a bootstrap
sampling method that selects a random subset of spectra, creates a mean and
rms spectrum, and then measures the line width directly from the data by
fitting a local linear continuum under each emission line and then
constraining the width of the flux in excess of the continuum. We ran 1000
realizations of the method and built up a distribution of $\sigma_{\rm line}$
and FWHM measurements, from which we report the median and 68% confidence
interval as the measurement and its uncertainty, respectively. As expected,
the line widths in the mean spectrum increased after narrow-line subtraction
had been applied (especially FWHM), while the line widths remained the same in
the rms spectrum before and after narrow-line subtraction. We report both FWHM
and $\sigma_{\rm line}$ for the mean and rms spectra, after narrow-line
subtraction, in Table 3.

All reported line widths were corrected for the resolution of the spectrograph
by assuming

$\Delta\lambda_{\rm obs}^{2}\approx\Delta\lambda_{\rm
true}^{2}+\Delta\lambda_{\rm disp}^{2}$ (2)

where $\Delta\lambda_{\rm obs}$ is the measured width of an emission line and
$\Delta\lambda_{\rm disp}$ is the broadening caused by the instrument. We
estimated $\Delta\lambda_{\rm true}$, the intrinsic line width, by adopting
the width of [O III] $\lambda 5007$ Å measured by Whittle (1992) through a
small spectrograph slit and at high resolution, FWHM$=550$ km s-1. Combining
$\Delta\lambda_{\rm true}$ with $\Delta\lambda_{\rm obs}$, we found
$\Delta\lambda_{\rm disp}=13.9$ Å.

##  4 Discussion

###  4.1 Radius vs. Luminosity

The mean and rms flux at $5100\times(1+z)$ Å throughout the monitoring
campaign was measured to be $F_{\lambda}(5100)=(3.49\pm 0.30)\times 10^{-15}$
erg s-1 cm-2 Å-1. However, with the large spectroscopic aperture employed, a
significant fraction of this flux is expected to come from stars in the host
galaxy rather than the AGN itself. IC 4329A was observed with the Hubble Space
Telescope Advanced Camera for Surveys High Resolution Channel through the
F550M filter in 2006 as part of program GO-10516 (see Figure 1). Bentz et al.
(2009a) report a 2D surface brightness decomposition of the image, from which
we determined the host galaxy contribution through the spectroscopic aperture
utilized in this work, $F_{\lambda,\rm gal}(5100)=(2.15\pm 0.22)\times
10^{-15}$ erg s-1 cm-2 Å-1.

Removing the starlight contribution to the flux measured at $5100\times(1+z)$
Å, correcting for Galactic absorption ($A_{B}=0.21$ mag, Schlafly & Finkbeiner
2011), and assuming a distance of $D=59\pm 9$ Mpc (as measured for IC 4329,
Tully et al. 2009), we find an average nuclear luminosity of $\log\lambda
L_{\lambda}(5100)=42.53\pm 0.15$ L⊙ for IC 4329A during this monitoring
campaign. Based on the relationship between BLR radius and AGN luminosity for
local Seyferts (Bentz et al., 2013), we would thus predict an H$\beta$ time
delay of $\sim 6$ days.

However, the observed nuclear luminosity is likely a severe underestimate
because of the intrinsic extinction along our line of sight to the center of
IC 4329A. Mehdipour & Costantini (2018) analyzed the spectral energy
distribution of IC 4329A from the X-rays through the far-infrared, finding
that the internal extinction was best described by a flat or gray extinction
curve, similar to that described by Czerny et al. (2004), with $E(B-V)=1.0\pm
0.1$ mag. These results suggest that the intrinsic AGN flux at
$5100\times(1+z)$ Å is $\sim 9.5\times$ larger than the observed flux.
Adjusting the nuclear luminosity by this factor then predicts an H$\beta$ time
delay of $\sim 18$ days, in good agreement with the measurements we present
here given the uncertainties involved in correcting internal extinction.

###  4.2 Black Hole Mass

The time delay for a broad emission line ($\tau$) gives a measure of the
light-crossing time and thus the physical size of the emission-line region,
while the width of the line ($V$) provides a constraint on the line-of-sight
velocity of the gas in the region. The two may be combined to derive the black
hole mass as:

$M_{\rm BH}=f\frac{c\tau V^{2}}{G}$ (3)

where $c$ is the speed of light, and $G$ is the gravitational constant. The
scale factor $f$ is an order unity term that includes the detailed geometry
and kinematics of the BLR as well as the inclination angle at which the AGN
system is viewed.

Historically, it has been difficult to constrain the BLR structure and
kinematics for most AGNs that have been studied with reverberation mapping.
Thus, a population-average scale factor, $\langle f\rangle$, is generally
adopted to bring the sample of reverberation-based masses into general
agreement with samples of black hole masses from dynamical modeling through
comparison of their $M_{\rm BH}-\sigma_{\star}$ relationships. Previous
studies have found values of $\langle f\rangle=2.8-5.5$, depending on the
details of the adopted samples and the fitting method (Onken et al., 2004;
Graham et al., 2011; Park et al., 2012; Grier et al., 2013a; Batiste et al.,
2017). These values of $\langle f\rangle$ are appropriate for time delays
given by $\tau_{\rm cent}$ and line widths given by $\sigma_{\rm line}$(rms),
the combination of measurements that has been shown to provide the smallest
scatter among black hole masses derived from multiple reverberation
experiments targeting the same AGN (Peterson et al., 2004). We note, however,
that while adopting $\langle f\rangle$ minimizes any bias in a sample of
reverberation masses, it also means that the reverberation mass for any single
AGN may be uncertain by a factor of a few. In particular, the inclinations of
individual AGNs to our line of sight are expected to be random, and the
analysis of Williams et al. (2018) suggests that inclination plays a
significant role in the variance of individual $f$ values determined for
different AGNs (see also, e.g., Collin et al. 2006).

Here we adopt $\langle f\rangle=4.8$ from Batiste et al. (2017) because of
their careful assessment of galaxy morphology on the determination of
$\sigma_{\star}$ among the AGN sample. Combined with $\tau_{\rm cent}$ and
$\sigma_{\rm line}$(rms) for H$\beta$, we find $M_{\rm
BH}=6.8^{+1.2}_{-1.1}\times 10^{7}$$M_{\odot}$. We find a similar mass with
slightly larger uncertainties, $M_{\rm BH}=5.5^{+1.9}_{-1.2}\times
10^{7}$$M_{\odot}$, using the measurements of H$\gamma$ and adopting the same
value of $\langle f\rangle$.

The $M_{\rm BH}-\sigma_{\star}$ relationship is commonly used to estimate
black hole masses in galaxies, although estimates for the black hole mass in
IC 4329A based on the $M_{\rm BH}-\sigma_{\star}$ relationship are likely to
be biased somewhat high because of the rotational contribution of the edge-on
galaxy disk to $\sigma_{\star}$. Using the $M_{\rm BH}-\sigma_{\star}$
relationship of Tremaine et al. (2002), Markowitz (2009) estimated $M_{\rm
BH}=2^{+2}_{-1}\times 10^{8}$$M_{\odot}$ for IC 4329A, which is indeed
somewhat higher than the reverberation mass we present here. Markowitz (2009)
also estimated the black hole mass using the McHardy et al. (2006)
relationship between bolometric luminosity, black hole mass, and X-ray power
spectral density break, finding $M_{\rm BH}=1.3^{+1.0}_{-0.3}\times
10^{8}$$M_{\odot}$. Markowitz (2009) assumed a luminosity distance of 78.6 Mpc
when calculating the bolometric luminosity, so adopting the measured distance
to IC 4329 of $D=59$ Mpc decreases the predicted mass by $\sim 23$%, bringing
it into better agreement with the reverberation mass. Finally, Ponti et al.
(2012) explored the relationship between X-ray excess variance and $M_{\rm
BH}$, using the reverberation sample to define a scaling relationship between
the two. With this scaling relationship, they predict $M_{\rm
BH}=2.1^{+1.3}_{-0.8}\times 10^{8}$$M_{\odot}$ for IC 4329A, although we note
that they adopted a mass estimate for IC 4329A based on $M_{\rm
BH}-\sigma_{\star}$ and included it as one of the few black holes with $M_{\rm
BH}>10^{8}$$M_{\odot}$ in the reverberation sample that defined the scaling
relationship.

Figure 5: Velocity-resolved time delays (top) across the profile of the broad
H$\beta$ emission line (bottom). The narrow-line contributions to the mean
spectrum are separately plotted in gray.

###  4.3 Velocity-Resolved Time Delays

While a population-average scale factor $\langle f\rangle$ is generally
adopted for reverberation masses, it was understood from the initial
development of reverberation mapping that the response as a function of
velocity across an emission line would provide information on the overall
geometry and kinematics of the BLR gas (Blandford & McKee, 1982). Mapping the
BLR via reverberation would then allow the black hole mass to be directly
constrained, without reliance on a scale factor. Initial attempts to constrain
velocity-resolved time delays were generally hampered by observational
uncertainties and systematics, but more recently there have been a growing
number of successes (e.g., Bentz et al. 2008; Denney et al. 2009; Grier et al.
2013b; Fausnaugh et al. 2017; Brotherton et al. 2020; Oknyansky et al. 2021;
Bao et al. 2022).

Therefore, in addition to measuring the time delay for the integrated H$\beta$
emission, we also investigated the possibility of constraining velocity-
resolved time delays across the broad emission line profile. Following the
methods described in Section 3, we divided the H$\beta$ profile into five
velocity bins. We created a light curve for each bin and measured the time
delay with respect to the continuum light curve.

In Figure 5 we show the time delays ($\tau_{\rm cent}$) measured as a function
of velocity across the H$\beta$ emission line. For comparison, the gray
horizontal band marks the time delay and its uncertainties for the integrated
emission line. A slight M-shaped pattern may be apparent in the time delays,
reminiscent of that seen in the velocity-resolved time delays of the broad UV
and optical emission lines in NGC 5548 (De Rosa et al., 2015; Pei et al.,
2017). Detailed analysis of the velocity-resolved signatures in NGC 5548 by
Horne et al. (2021) revealed the M-shaped pattern to be the signature of an
inclined Keplerian disk. IC 4329A has long been known (Disney, 1973; Wilson &
Penston, 1979) as one of the few nearby AGNs to show clear double-peaked broad
emission line profiles. Storchi-Bergmann et al. (2017) suggest that double
peaks in optical broad emission lines are more common among local Seyferts
than generally recognized, and they include NGC 5548 among their sample.
Double-peaked broad emission lines are often interpreted as indicating a disk-
like broad line region structure, so the potential similarities between the
velocity-resolved responses of NGC 5548 and IC 4329A are perhaps unsurprising
if they are both double-peaked AGNs. Nevertheless, while velocity-resolved
time delays may provide hints, a more thorough analysis is required to
constrain the details of the BLR geometry in IC 4329A.

###  4.4 CARAMEL Modeling

There are two independent approaches for analyzing velocity-resolved
reverberation responses: as an ill-posed inverse problem or through forward
modeling. The inverse approach attempts to extract the transfer function,
which describes the time delay distribution as a function of velocity,
directly from the data (Horne, 1994; Skielboe et al., 2015; Anderson et al.,
2021). Forward modeling, on the other hand, uses a framework of self-
consistent models to explore the relevant parameter space and determine the
set of models that best agree with the observations (Pancoast et al., 2011,
2014). Forward modeling has the advantage that the results are relatively
simple to interpret, although it is limited by the flexibility and
completeness of the models. The inverse approach, on the other hand, makes few
initial assumptions and is generally more flexible, but produces results that
may be difficult to interpret and generally rely on comparison with models.

We explored the possibility of constraining more details regarding the
geometry and kinematics of the BLR in IC 4329A with the forward modeling code
CARAMEL, following a similar approach to that we previously employed (Bentz et
al., 2021b, 2022). In this case, most of the model parameters were not well
constrained. Parameters with modest constraints included the black hole mass,
$\log($$M_{\rm BH}$/$M_{\odot}$$)=7.64^{+0.53}_{-0.25}$ or $M_{\rm
BH}=4^{+10}_{-2}\times 10^{7}$ $M_{\odot}$, and the average time delay,
$\tau_{\rm mean}=13.6^{+9.2}_{-3.0}$ days, which agree with our results in
Sections 3.1 and 4.2.

Figure 6: Posterior distributions of constraints on the H$\beta$-emitting BLR
geometry.

While the dynamics of the system were not constrained, there were weak
constraints on some of the geometric parameters of the H$\beta$-emitting BLR
(Figure 6), suggesting an average radius of $r_{\rm mean}=14.2^{+7.2}_{-3.7}$
light-days, with a minimum radius of $r_{\rm min}=3.1^{+8.5}_{-2.5}$ light-
days and a radial extent of $\sigma_{\rm r}=12.0^{+22.5}_{-8.0}$ light-days.

Finally, modest constraints were derived for the inclination angle of the BLR
relative to our line of sight, $\theta_{i}=40.1^{+26.6}_{-18.1}$ deg, as well
as the opening angle of the BLR, $\theta_{o}=69.2^{+15.3}_{-29.6}$ deg, where
$\theta_{o}=0\arcdeg$ is a thin disk and $\theta_{o}=90\arcdeg$ is a sphere.
The inclination angle suggests that the AGN system and the galaxy disk are
indeed misaligned by a significant amount ($\sim 25\arcdeg-70\arcdeg$).
Furthermore, we note that the moderate inclination angle constrained by the
models suggests that the value of $\langle f\rangle$ adopted in Section 4.2 is
unlikely to be an underestimate for this particular AGN (cf. Williams et al.
2018).

Visualization of a random sampling of allowed BLR models included wide and
puffy inclined annuli, filled inclined biconical structures, and inclined
“hamburger buns” that resemble the intersection of a bicone with a shell.
Thus, while these models demonstrate consistency with several previous results
including those that we present above, they are limited in their ability to
elucidate a more clear picture of the BLR in IC 4329A. Additional work will be
required to determine whether these limitations are a result of the quality of
the reverberation data, which are a bit noisy but do show a clear velocity-
resolved response, or are instead due to a mismatch between the modeling
assumptions in CARAMEL and the physical conditions of the BLR in IC 4329A.

##  5 Summary

We have presented a new reverberation mapping campaign focused on the AGN in
the center of the edge-on spiral galaxy IC 4329A. With photometric and
spectroscopic monitoring covering 2022 Feb–July, we were able to constrain the
time delay of the broad H$\beta$ and H$\gamma$ emission lines:
$16.3^{+2.6}_{-2.3}$ days for H$\beta$ and $16.0^{+4.8}_{-2.6}$ days for
H$\gamma$. He II $\lambda 4686$ also exhibited variability during these
observations, but the time delay of $-0.6^{+3.9}_{-3.9}$ days is unresolved
and consistent with zero days. The H$\beta$ time delay is consistent with that
expected from the relationship between AGN luminosity and BLR radius for local
Seyferts, after correcting the AGN luminosity for the expected $\sim 2.4$ mag
of intrinsic extinction at 5100 Å.

Combining the time delay and width of H$\beta$, we find a black hole mass of
$M_{\rm BH}=6.8^{+1.2}_{-1.1}\times 10^{7}$$M_{\odot}$ when adopting $\langle
f\rangle=4.8$. Based on the detection of velocity-resolved time delays across
the H$\beta$ broad line profile, we carried out forward modeling of the
reverberation data presented here. Many of the model parameters were
relatively unconstrained, but the black hole mass and average H$\beta$ time
delay agreed with our measurements. Finally, the models suggest that the BLR
is indeed severely misaligned with the disk of the host galaxy, in agreement
with expectations from the unified model for AGNs given that broad emission
lines are clearly visible in the spectrum of the nucleus of IC 4329A, and yet
the galaxy is oriented edge-on to our line of sight.

We thank the referees for suggestions that improved the presentation of this
work. We thank Tommaso Treu and Peter Williams for providing the CARAMEL
modeling code as well as assistance in working with it. MCB gratefully
acknowledges support from the NSF through grant AST-2009230. CAO was supported
by the Australian Research Council (ARC) through Discovery Project
DP190100252. MV was supported by the NSF through AST-2009122. LCOGT

## References

  * Abell et al. (1989) Abell, G. O., Corwin, Harold G., J., & Olowin, R. P. 1989, ApJS, 70, 1, doi: 10.1086/191333
  * Alard (2000) Alard, C. 2000, A&AS, 144, 363, doi: 10.1051/aas:2000214
  * Alard & Lupton (1998) Alard, C., & Lupton, R. H. 1998, ApJ, 503, 325, doi: 10.1086/305984
  * Anderson et al. (2021) Anderson, M. D., Baron, F., & Bentz, M. C. 2021, MNRAS, 505, 2903, doi: 10.1093/mnras/stab1394
  * Antonucci (1993) Antonucci, R. 1993, ARA&A, 31, 473, doi: 10.1146/annurev.aa.31.090193.002353
  * Bao et al. (2022) Bao, D.-W., Brotherton, M. S., Du, P., et al. 2022, ApJS, 262, 14, doi: 10.3847/1538-4365/ac7beb
  * Batiste et al. (2017) Batiste, M., Bentz, M. C., Raimundo, S. I., Vestergaard, M., & Onken, C. A. 2017, ApJ, 838, L10, doi: 10.3847/2041-8213/aa6571
  * Bentz et al. (2013) Bentz, M. C., Denney, K. D., Grier, C. J., et al. 2013, ApJ, 767, 149, doi: 10.1088/0004-637X/767/2/149
  * Bentz & Katz (2015) Bentz, M. C., & Katz, S. 2015, PASP, 127, 67, doi: 10.1086/679601
  * Bentz et al. (2009a) Bentz, M. C., Peterson, B. M., Netzer, H., Pogge, R. W., & Vestergaard, M. 2009a, ApJ, 697, 160, doi: 10.1088/0004-637X/697/1/160
  * Bentz et al. (2021a) Bentz, M. C., Street, R., Onken, C. A., & Valluri, M. 2021a, ApJ, 906, 50, doi: 10.3847/1538-4357/abccd4
  * Bentz et al. (2009b) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2009b, ApJ, 705, 199, doi: 10.1088/0004-637X/705/1/199
  * Bentz et al. (2021b) Bentz, M. C., Williams, P. R., Street, R., et al. 2021b, ApJ, 920, 112, doi: 10.3847/1538-4357/ac19af
  * Bentz et al. (2022) Bentz, M. C., Williams, P. R., & Treu, T. 2022, arXiv e-prints, arXiv:2206.03513.  https://arxiv.org/abs/2206.03513
  * Bentz et al. (2008) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2008, ApJ, 689, L21, doi: 10.1086/595719
  * Bernardi et al. (2002) Bernardi, M., Alonso, M. V., da Costa, L. N., et al. 2002, AJ, 123, 2990, doi: 10.1086/340463
  * Blandford & McKee (1982) Blandford, R. D., & McKee, C. F. 1982, ApJ, 255, 419, doi: 10.1086/159843
  * Brotherton et al. (2020) Brotherton, M. S., Du, P., Xiao, M., et al. 2020, ApJ, 905, 77, doi: 10.3847/1538-4357/abc2d2
  * Brown et al. (2018) Brown, J. S., Kochanek, C. S., Holoien, T. W. S., et al. 2018, MNRAS, 473, 1130, doi: 10.1093/mnras/stx2372
  * Brown et al. (2013) Brown, T. M., Baliber, N., Bianco, F. B., et al. 2013, PASP, 125, 1031, doi: 10.1086/673168
  * Cackett et al. (2021) Cackett, E. M., Bentz, M. C., & Kara, E. 2021, iScience, 24, 102557, doi: https://doi.org/10.1016/j.isci.2021.102557
  * Collin et al. (2006) Collin, S., Kawaguchi, T., Peterson, B. M., & Vestergaard, M. 2006, A&A, 456, 75
  * Czerny et al. (2004) Czerny, B., Li, J., Loska, Z., & Szczerba, R. 2004, MNRAS, 348, L54, doi: 10.1111/j.1365-2966.2004.07590.x
  * 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
  * Denney et al. (2009) Denney, K. D., Peterson, B. M., Pogge, R. W., et al. 2009, ApJ, 704, L80, doi: 10.1088/0004-637X/704/2/L80
  * Disney (1973) Disney, M. J. 1973, ApJ, 181, L55, doi: 10.1086/181184
  * 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
  * Fausnaugh et al. (2017) Fausnaugh, M. M., Grier, C. J., Bentz, M. C., et al. 2017, ApJ, 840, 97, doi: 10.3847/1538-4357/aa6d52
  * Gaskell & Peterson (1987) Gaskell, C. M., & Peterson, B. M. 1987, ApJS, 65, 1, doi: 10.1086/191216
  * Gaskell & Sparke (1986) Gaskell, C. M., & Sparke, L. S. 1986, ApJ, 305, 175, doi: 10.1086/164238
  * Graham et al. (2011) Graham, A. W., Onken, C. A., Athanassoula, E., & Combes, F. 2011, MNRAS, 412, 2211, doi: 10.1111/j.1365-2966.2010.18045.x
  * Grier et al. (2013a) Grier, C. J., Martini, P., Watson, L. C., et al. 2013a, ApJ, 773, 90, doi: 10.1088/0004-637X/773/2/90
  * Grier et al. (2017) Grier, C. J., Pancoast, A., Barth, A. J., et al. 2017, ApJ, 849, 146, doi: 10.3847/1538-4357/aa901b
  * Grier et al. (2013b) Grier, C. J., Peterson, B. M., Horne, K., et al. 2013b, ApJ, 764, 47, doi: 10.1088/0004-637X/764/1/47
  * Hamuy et al. (1994) Hamuy, M., Suntzeff, N. B., Heathcote, S. R., et al. 1994, PASP, 106, 566, doi: 10.1086/133417
  * Hamuy et al. (1992) Hamuy, M., Walker, A. R., Suntzeff, N. B., et al. 1992, PASP, 104, 533, doi: 10.1086/133028
  * Hartman et al. (2004) Hartman, J. D., Bakos, G., Stanek, K. Z., & Noyes, R. W. 2004, AJ, 128, 1761, doi: 10.1086/423920
  * Hartman et al. (2005) Hartman, J. D., Stanek, K. Z., Gaudi, B. S., Holman, M. J., & McLeod, B. A. 2005, AJ, 130, 2241, doi: 10.1086/462405
  * Henden & Munari (2014) Henden, A., & Munari, U. 2014, Contributions of the Astronomical Observatory Skalnate Pleso, 43, 518
  * Holoien et al. (2016) Holoien, T. W. S., Kochanek, C. S., Prieto, J. L., et al. 2016, MNRAS, 463, 3813, doi: 10.1093/mnras/stw2272
  * Horne (1994) Horne, K. 1994, in Astronomical Society of the Pacific Conference Series, Vol. 69, Reverberation Mapping of the Broad-Line Region in Active Galactic Nuclei, ed. P. M. Gondhalekar, K. Horne, & B. M. Peterson, 23
  * Horne et al. (2021) Horne, K., De Rosa, G., Peterson, B. M., et al. 2021, ApJ, 907, 76, doi: 10.3847/1538-4357/abce60
  * Hudson et al. (2001) Hudson, M. J., Lucey, J. R., Smith, R. J., Schlegel, D. J., & Davies, R. L. 2001, MNRAS, 327, 265, doi: 10.1046/j.1365-8711.2001.04786.x
  * King & Pringle (2007) King, A. R., & Pringle, J. E. 2007, MNRAS, 377, L25, doi: 10.1111/j.1745-3933.2007.00296.x
  * Lousto et al. (2016) Lousto, C. O., Healy, J., & Nakano, H. 2016, Phys. Rev. D, 93, 044031, doi: 10.1103/PhysRevD.93.044031
  * Malin & Carter (1983) Malin, D. F., & Carter, D. 1983, ApJ, 274, 534, doi: 10.1086/161467
  * Markowitz (2009) Markowitz, A. 2009, ApJ, 698, 1740, doi: 10.1088/0004-637X/698/2/1740
  * McCully et al. (2018) McCully, C., Volgenau, N. H., Harbeck, D.-R., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10707, Software and Cyberinfrastructure for Astronomy V, ed. J. C. Guzman & J. Ibsen, 107070K, doi: 10.1117/12.2314340
  * McHardy et al. (2006) McHardy, I. M., Koerding, E., Knigge, C., Uttley, P., & Fender, R. P. 2006, Nature, 444, 730, doi: 10.1038/nature05389
  * Mehdipour & Costantini (2018) Mehdipour, M., & Costantini, E. 2018, A&A, 619, A20, doi: 10.1051/0004-6361/201833706
  * Melinder et al. (2008) Melinder, J., Mattila, S., Östlin, G., Mencía Trinchant, L., & Fransson, C. 2008, A&A, 490, 419, doi: 10.1051/0004-6361:20079206
  * Miknaitis et al. (2007) Miknaitis, G., Pignata, G., Rest, A., et al. 2007, ApJ, 666, 674, doi: 10.1086/519986
  * Morris & Ward (1988) Morris, S. L., & Ward, M. J. 1988, MNRAS, 230, 639, doi: 10.1093/mnras/230.4.639
  * Nagar & Wilson (1999) Nagar, N. M., & Wilson, A. S. 1999, ApJ, 516, 97, doi: 10.1086/307109
  * Oknyansky et al. (2021) Oknyansky, V. L., Brotherton, M. S., Tsygankov, S. S., et al. 2021, MNRAS, 505, 1029, doi: 10.1093/mnras/stab1138
  * Onken et al. (2004) Onken, C. A., Ferrarese, L., Merritt, D., et al. 2004, ApJ, 615, 645, doi: 10.1086/424655
  * Pancoast et al. (2011) Pancoast, A., Brewer, B. J., & Treu, T. 2011, ApJ, 730, 139, doi: 10.1088/0004-637X/730/2/139
  * Pancoast et al. (2014) —. 2014, MNRAS, 445, 3055, doi: 10.1093/mnras/stu1809
  * Park et al. (2012) Park, D., Kelly, B. C., Woo, J.-H., & Treu, T. 2012, ApJS, 203, 6, doi: 10.1088/0067-0049/203/1/6
  * Pei et al. (2017) Pei, L., Fausnaugh, M. M., Barth, A. J., et al. 2017, ApJ, 837, 131, doi: 10.3847/1538-4357/aa5eb1
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266, doi: 10.1086/340952
  * Peng et al. (2010) —. 2010, AJ, 139, 2097, doi: 10.1088/0004-6256/139/6/2097
  * Peterson et al. (2013) Peterson, B. M., Denney, K. D., De Rosa, G., et al. 2013, ApJ, 779, 109, doi: 10.1088/0004-637X/779/2/109
  * Peterson et al. (1998) Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110, 660
  * Peterson et al. (2004) Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682, doi: 10.1086/423269
  * Ponti et al. (2012) Ponti, G., Papadakis, I., Bianchi, S., et al. 2012, A&A, 542, A83, doi: 10.1051/0004-6361/201118326
  * Read & Pietsch (1998) Read, A. M., & Pietsch, W. 1998, A&A, 336, 855.  https://arxiv.org/abs/astro-ph/9805214
  * Riess et al. (2001) Riess, A. G., Nugent, P. E., Gilliland, R. L., et al. 2001, ApJ, 560, 49, doi: 10.1086/322348
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schmitt et al. (1997) Schmitt, H. R., Kinney, A. L., Storchi-Bergmann, T., Antonucci, & Robert. 1997, ApJ, 477, 623, doi: 10.1086/303744
  * Siverd et al. (2012) Siverd, R. J., Beatty, T. G., Pepper, J., et al. 2012, ApJ, 761, 123, doi: 10.1088/0004-637X/761/2/123
  * Skielboe et al. (2015) Skielboe, A., Pancoast, A., Treu, T., et al. 2015, MNRAS, 454, 144, doi: 10.1093/mnras/stv1917
  * Stahl et al. (2021) Stahl, B. E., de Jaeger, T., Boruah, S. S., et al. 2021, MNRAS, 505, 2349, doi: 10.1093/mnras/stab1446
  * Storchi-Bergmann et al. (2017) Storchi-Bergmann, T., Schimoia, J. S., Peterson, B. M., et al. 2017, ApJ, 835, 236, doi: 10.3847/1538-4357/835/2/236
  * Tremaine et al. (2002) Tremaine, S., et al. 2002, ApJ, 574, 740
  * Tully et al. (2009) Tully, R. B., Rizzi, L., Shaya, E. J., et al. 2009, AJ, 138, 323, doi: 10.1088/0004-6256/138/2/323
  * Udalski et al. (2008) Udalski, A., Szymanski, M. K., Soszynski, I., & Poleski, R. 2008, Acta Astron., 58, 69.  https://arxiv.org/abs/0807.3884
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, Acta Astron., 65, 1\.  https://arxiv.org/abs/1504.05966
  * van Groningen & Wanders (1992) van Groningen, E., & Wanders, I. 1992, PASP, 104, 700, doi: 10.1086/133039
  * White & Peterson (1994) White, R. J., & Peterson, B. M. 1994, PASP, 106, 879, doi: 10.1086/133456
  * Whittle (1992) Whittle, M. 1992, ApJS, 79, 49, doi: 10.1086/191644
  * Williams et al. (2018) Williams, P. R., Pancoast, A., Treu, T., et al. 2018, ApJ, 866, 75, doi: 10.3847/1538-4357/aae086
  * Wilson & Penston (1979) Wilson, A. S., & Penston, M. V. 1979, ApJ, 232, 389, doi: 10.1086/157298
  * Winge et al. (1996) Winge, C., Peterson, B. M., Pastoriza, M. G., & Storchi-Bergmann, T. 1996, ApJ, 469, 648, doi: 10.1086/177812
  * Wolstencroft et al. (1995) Wolstencroft, R. D., Done, C. J., Scarrott, S. M., & Scarrott, R. M. J. 1995, MNRAS, 276, 460, doi: 10.1093/mnras/276.2.460
  * Zebrun et al. (2001) Zebrun, K., Soszynski, I., & Wozniak, P. R. 2001, Acta Astron., 51, 303.  https://arxiv.org/abs/astro-ph/0110612
  * Zu et al. (2011) Zu, Y., Kochanek, C. S., & Peterson, B. M. 2011, ApJ, 735, 80, doi: 10.1088/0004-637X/735/2/80


# LAMMPS Implementation of Rapid Artificial Neural Network Derived Interatomic
Potentials

D.Dickel1,2, M. Nitol2, C. D. Barrett1,2

###### Abstract

While machine learning approaches have been successfully used to represent
interatomic potentials, their speed has typically lagged behind conventional
formalisms. This is often due to the complexity of the structural fingerprints
used to describe the local atomic environment and the large cutoff radii and
neighbor lists used in the calculation of these fingerprints. Even recent
machine learned methods are at least 10 times slower than traditional
formalisms. An implementation of a rapid artificial neural network (RANN)
style potential in the LAMMPS molecular dynamics package is presented here
which utilizes angular screening to reduce computational complexity without
reducing accuracy. For the smallest neural network architectures, this
formalism rivals the modified embedded atom method (MEAM) for speed and
accuracy, while the networks approximately one third as fast as MEAM were
capable of reproducing the training database with chemical accuracy. The
numerical accuracy of the LAMMPS implementation is assessed by verifying
conservation of energy and agreement between calculated forces and pressures
and the observed derivatives of the energy as well as by assessing the
stability of the potential in dynamic simulation. The potential style is
tested using a force field for magnesium and the computational efficiency for
a variety of architectures is compared to a traditional potential models as
well as alternative ANN formalisms. The predictive accuracy is found to rival
that of slower methods.

Keywords. Machine Learning; Artificial Neural Networks; Molecular Dynamics;

##  1 Introduction

Artificial neural networks (ANNs) and other machine learning techniques are
gaining prominence as a method for developing classical interatomic potentials
for use in molecular statics (MS) and molecular dynamics (MD) simulation [1,
2, 3, 4, 5, 6, 7, 8, 9, 10]. Because of their ability to accurately fit large
and complex data sets and because they can be rapidly optimized, ANNs have
shown the potential to predict atomic energies and forces more accurately than
most existing classical formalisms [6, 7]. By using density functional theory
(DFT), which is relatively expensive computationally, to produce large
training databases for the ANN, potentials which reproduce this data to high
precision at much faster computational speeds can be created. While only a
limited number of large scale dynamic simulations have been performed using
these potentials, static simulations for a variety of materials have been able
to reproduce ab initio forces and energies with chemical accuracy [11, 2, 12,
5].
ANNs used to model interatomic potentials behave similarly to traditional
formalisms with an energy being defined for each individual atom based on the
local atomic arrangement and the atomic species of it and its neighbors,
mapping the local environment directly to numerical values for the energy and
forces. This encoding of the local environment is referred to as the
structural fingerprint. Thompson et al. developed the Spectral Neighbor
Analysis Method (SNAP) [13] which uses machine-learning techniques on
components of the local neighbor density. Behler and Parrinello [14, 15]
defined a set of functions based on Gaussians which could be used to generate
the fingerprint. This fingerprint formalism has since been used by a number of
authors as the basis for an ANN. The celebrated Gaussian Approximation
Potentials (GAP)[16] have been used to develop potentials with chemical
accuracy for a variety of systems including tungsten [17], lithium-carbon[18],
iron[19] and water[20]. Artrith and Urban used the Behler-Parrinello basis
functions to create a potential for titanium dioxide using an ANN with two
hidden layers with 10 neurons each which accurately described the various
cystal phases of TiO2[2]. Artrith, Urban, and Ceder have also recently
proposed an efficient scheme for the generation of machine-learned potentials
for composite materials with many atomic species[21]. Takahahsi et al [11]
used a combination of the Gaussian formalism along with insights from the
modified embedded atom method (MEAM) [22] to generate a fingerprint of several
thousand values which was then used to fit a potential for titanium using
linear ridge regression. Recently, a method using the embedding functions from
MEAM has been used to generate potentials for titanium and zinc using a single
hidden layer [23, 24]. Singraber, Behler, and Dellago [25] recently presented
the n2p2 neural network package which uses the Behler symmetry functions in
the construction of the fingerprint. This package has been used to develop a
number of potentials including for Al-Cu [26] and Mg [27].
While the optimal structural fingerprint and network architecture for an ANN
is still unclear and indeed is probably strongly dependent on the material
system being examined, there is a need to be able to perform large scale
dynamic calculations with these potentials. Additionally, while machine
learned potentials greatly outperform ab initio calculations in terms of
efficiency, they still lag behind conventional formalisms such as the embedded
atom method (EAM)[28] and MEAM, which have been the standard for analysis of
metals since their introduction over 25 years ago. Additionally, MEAM takes
advantage of angular screening to limit the size of the neighbor list and
incorporate the known shielding effect to improve performance. To this end, we
have created an implementation of a rapid artificial neural network (RANN)
potential formalism for use in LAMMPS [29]. This implementation accepts
arbitrary Multi-layer Perceptron (MLP) architectures and structural
fingerprints of arbitrary size. Sample potentials have been generated from the
MEAM based structural fingerprints with angular screening, though novel
functions or fingerprint types can be added independently. Section 2 of this
paper reviews the formalism for ANN interatomic potentials and provides the
specific structural fingerprints implemented including a discussion of the use
of angular screening. Section 3 describes the framework which has been
implemented in LAMMPS. Section 4 demonstrates a number of simple static and
dynamic calculations using an ANN potential. In particular, conservation of
energy is shown along with accurate calculation of the forces and pressure,
and a comparison is done between the computational efficiency of this
potential style as compared to other common styles for metal systems. It is
demonstrated that, for the network architectures capable of reproducing a
magnesium training database to within 1meV/atom, performance approximately one
third as MEAM can be achieved. Additionally, the performance of the potential
as a function of neural network architecture and neighbor list size is
considered.

##  2 Artificial Neural Network Architecture

Most machine-learned interatomic potentials make use of MLPs to predict the
energy of individual atoms in a system given their environment. Details of the
training and optimization of these networks in reference to ab initio data can
be found elsewhere [2, 23, 30], but the general structure of a network as an
interatomic potential is as follows. The energy of a particular atom $i$,
determined by its environment, is the last of N layers of the neural network.
The values for any particular layer, $\mathbf{A}^{n}$, after the first is
determined by the previous layer and the weight and bias matrices
$\mathbf{W}^{n}$ and $\mathbf{B}^{n}$:

$\displaystyle
Z^{n}_{l_{n}}=\sum_{l_{n-1}}{W^{n}_{l_{n}l_{n-1}}A^{n-1}_{l_{n-1}}+B^{n}_{l_{n}}}$
(1) $\displaystyle A^{n}_{l_{n}}=g^{n}(Z^{n}_{l_{n}})$ (2)

Where $l_{n}$ is the number of neurons in layer $n$ and $g^{n}(x)$ is a
nonlinear activation function. The input layer, $\mathbf{A}^{0}$ is given by
the structural fingerprint of the local atomic environment. The ouput layer,
$\mathbf{Z}^{N}$, will always contain a single node, representing the energy,
and so $\mathbf{W}^{N}$ will always be a row vector and $\mathbf{B}^{N}$
always a single number.
A number of different functions can be used to describe this structural
fingerprint. Behler and Parrinello [15] consider two types of fingerprint
functions: radial symmetry functions which are a pairwise sum of Gaussians,
and angular terms over all triplets of atoms. For the RANN style, however, we
use the fingerprint style introduced by Dickel, Francis, and Barrett,
motivated by the Modified Embedded Atom Method (MEAM) formalism [22, 31, 32],
with the addition of angular screening. This style seems particularly relevant
for metal potentials for the same physically motivated reasons that the MEAM
formalism is effective. In this style, two different kind of input
fingerprints are considered. First, simple pair interactions are considered
and summed over all the neighbors of a given atom. For a given atom labeled
$i$, we define a set of pair potentials interactions with the form

$F_{n}=\sum_{j\neq
i}{(\frac{r_{ij}}{r_{e}})^{n}e^{-\alpha_{n}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta
r})S_{ij}}$ (3)

Where $j$ labels all the neighbors atoms of $i$ within a cutoff radius
$r_{c}$, $n$ is an integer, different for each member of the pairwise
contributions to the fingerprint, $r_{e}$ is the equilibrium nearest neighbor
distance, $S_{ij}$ is an angular screening term described below, and
$\alpha_{n}$ are metaparameters, which can be tuned to better optimize the
potential. In principle, $\alpha_{n}$ are related to the dimensionless
$\alpha$ used in the MEAM formalism which is related to the bulk modulus of
the material [32]. $f_{c}(x)$ is a cutoff function which smoothly varies the
weight of the interaction from 1 for atoms inside the cutoff radius to 0 once
$r_{ij}>r_{c}$ where $\Delta r$ defines the width of the transition. The
cutoff function used here is the same as employed in MEAM and is given by

$r_{c}(x)=\left\\{\begin{array}[]{lr}1,&x>1\\\ (1-(1-x)^{4})^{2},&0\leq x\leq
1\\\ 0,&x<0\end{array}\right\\}$ (4)

The second kind of fingerprint function considers three body terms, with a
form similar to the partial electron densities used in MEAM. There are two
parameters, $m$ and $k$ which can be varied to generate more inputs for the
fingerprint:

$G_{m,k}=\sum_{j,k}{\cos^{m}{\theta_{jik}}e^{-\beta_{k}\frac{r_{ij}+r_{ik}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta
r})f_{c}(\frac{r_{c}-r_{ik}}{\Delta r})S_{ij}S_{ik}}$ (5)

where $\theta_{jik}$ is the angle between $r_{ij}$ and $r_{ik}$, m is a non-
negative interger and $\beta_{k}$ is a set of metaparameters which determine
the length scale of the different terms. For large numbers of neighbors, it
can be more efficient to calculate this as a sum over a single list of
neighbors as follows.

$\displaystyle
G_{0,k}=\left(\sum_{j}{e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(r_{ij}})S_{ij}\right)^{2}$
(6) $\displaystyle
G_{1,k}=\sum_{\alpha_{1}}\left(\sum_{j}\frac{r^{\alpha_{1}}_{i,j}}{r_{i,j}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta
r})S_{ij}\right)^{2}$ (7) $\displaystyle
G_{2,k}=\sum_{\alpha_{1},\alpha_{2}}\left(\sum_{j}\frac{r^{\alpha_{1}}_{i,j}r^{\alpha_{2}}_{i,j}}{r_{i,j}^{2}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta
r})S_{ij}\right)^{2}$ (8) $\displaystyle
G_{3,k}=\sum_{\alpha_{1},\alpha_{2},\alpha_{3}}\left(\sum_{j}\frac{r^{\alpha_{1}}_{i,j}r^{\alpha_{2}}_{i,j}r^{\alpha_{3}}_{i,j}}{r_{i,j}^{3}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta
r})S_{ij}\right)^{2}$ (9) $\displaystyle
G_{m,k}=\sum_{\alpha_{1}}\sum_{\alpha_{2}}...\sum_{\alpha_{m}}\left(\sum_{j}\frac{r^{\alpha_{1}}_{i,j}r^{\alpha_{2}}_{i,j}...r^{\alpha_{m}}_{i,j}}{r_{i,j}^{m}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta
r}))S_{ij}\right)^{2}$ (10)

where $r^{\alpha_{m}}_{ij}$ is the $x$, $y$, or $z$ component of $r_{ij}$ for
$\alpha_{m}=1,2$ or $3$, respectively. A simplification can be made by noting
that many of the terms in the three body interactions are redundant. In
$G_{2,k}$ for example, there will be a calculation for xy as well as for yx.
By taking advantage of this redundancy, we can reduce the number of terms that
need to be calculated for $G_{m,k}$ from $3^{n}$ to $\frac{(n+1)(n+2)}{2}$.
The simplified expression for $G_{m,k}$ now reads:

$G_{m,k}=\sum_{\alpha_{1}\geq\alpha_{2}\geq...\geq\alpha_{m}}\left(\frac{n!}{n_{x}!n_{y}!n_{z}!}\sum_{j}{\frac{r^{\alpha_{1}}_{i,j}r^{\alpha_{2}}_{i,j}...r^{\alpha_{m}}_{i,j}}{r_{i,j}^{m}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta
r})}\right)^{2}$ (11)

where $n_{x}$, $n_{y}$, and $n_{z}$ are the total number of x-, y-, and
z-components in the set $(\alpha_{1},\alpha_{2},...,\alpha_{m})$,
respectively. Note that which of these forms (Eq. 5 or Eq.11 will be more
efficiently calculated will depend on the length of the neighbor list and the
magnitude of $m$. Structural fingerprints of arbitrary size can be constructed
by varying the number of different $n$, $m$, and $k$ used giving more or less
information to the ANN. For the potentials tested here, we have used
consecutive values for n and m and totaling $n_{tot}$ and $m_{tot}$
respectively ($n\in(-1,0,1,...,n_{tot}-2)$, $m\in(0,1,...,m_{tot}-1)$). A
particular ANN potential will consist of a fixed structural fingerprint,
number and length of each hidden layer, weight and bias matrices for each
layer and activation functions for each layer. For the potentials considered
here we have used the following activation functions:

$\displaystyle g^{N}(x)=x$ (12) $\displaystyle
g^{n}(x)=\frac{x}{10}+\frac{9}{10}\log(e^{x}+1)\;\textrm{for}\;n<N$ (13)

The size of the weight and bias matrices will depend on the length of the
fingerprint and the length of each hidden layer.
As discussed above, one of the major computational barriers in the
implementation of ANN methods is the large neighbor list included in the
calculation of the structural fingerprint. In order to obtain good convergence
with DFT results, cutoff radii used can be larger than 12 Å[27], which for the
ground state of Mg includes over 300 atoms. By comparison, recent MEAM
potentials for Mg [33] have used a cutoff distance of 5.875 Å, which includes
only 38 atoms in the ground state. Since the computation time scales at least
linearly with the number of neighbors, minimizing the length of the neighbor
list without sacrificing performance should be a key goal for optimal
efficiency. The radial screening function $f_{c}$ can accomplish this, but
leads to nonphysical results at large separations (see, for example [33]).
Angular screening, whereby the effective interaction between atoms is reduced
or eliminated by the presence of an atom located between them can effectively
limit the neighbor list in a similar way without the unphysical results. Such
a method has been utilized by MEAM [34], and the same screening method has
been employed here. Briefly, the screening between two atoms is determined
from the product of all the screening interactions by other atoms in the
neighborhood:

$S_{ij}=\prod_{k\neq i,j}S_{ikj}$ (14)

where $S_{ikj}$ is calculated from a geometric construction considering the
ellipse formed by the 3 atoms with $r_{i,j}$ one of the axes. The screening
parameter, $C_{ikj}$ is then given by:

$C_{ikj}=1+2\frac{r^{2}_{ij}r^{2}_{ik}+r^{2}_{ij}r^{2}_{jk}-r^{4}_{ij}}{r^{4}_{ij}-(r^{2}_{ik}-r^{2}_{jk})^{2}}$
(15)

and then screening value is

$S_{ikj}=f_{c}\left(\frac{C_{ikj}-C_{min}}{C_{max}-C_{min}}\right)$ (16)

where $f_{c}$ is the same cutoff function used for the radial cutoff and
$C_{max}$ and $C_{min}$ are metaparameters which can be tuned to determine
which neighbors can be excluded from calculations.
The effect of including angular screening can be demonstrating by considering
the change of an individual fingerprint as the length scale is changed
continuously. For this example, we consider the value of a single fingerprint
for a Mg atom in a perfect bulk hcp lattice as the lattice parameter is
changed continuously. In the absence of angular screening, the value of the
fingerprint will change rapidly at particular values of the lattice constant
as new neighbors enter the radial screening distance. If angular screening is
included with metaparameters such that, for example, only 3rd nearest
neighbors are ever included regardless of lattice parameter, the change in the
value is considerably smoother. Figure 1 demonstrates the difference between
the two cases. Not only is the unscreened value more expensive to calculate as
it requires more neighbors for most cases, it can also be seen that the
predictive power of the model is hampered as the relation between two states
which are physically similar – ideal 3D lattices with different lattice
constants, given markedly different trends and values as the number of
neighbors changes. This can be overcome through the action of the neural
network and precise cancellation among different fingerprint terms, but the
quality of fit and the extrapolative power of the screened functions appears
more natural. Indeed unphysical deviations at relatively moderate isotropic
compressions, as seen for example in [27] can be avoided entirely through
angular screening as the fingerprint values change continuously in the
screened case even at high compression.

Figure 1: A sample value for a fingerprint in the presence (orange) and
absence (blue)of angular screening as the volume of a cell is changed
isotropically. Because the neighbor list frequently changes without angular
screening, the curve is seen to be less smooth, requiring more data to provide
predictive information to the ANN.

For the full network, forces are calculated by taking a gradient of the energy
over the entire system. This includes the individual ANN contributions from
every atom in the system, so the force on atom $i$ is determined not only by
the ANN output of that atom, but of all of its neighbors as well.

##  3 Numerical Implementation

The ANN potential file written for LAMMPS includes all of the details of the
neural network. It begins with a header which lists the types and numbers of
fingerprints to be used. While we have only included the two varieties, the
bond power and radial power terms, $F_{n}$ and $G_{m,k}$, in this work, the
LAMMPS potential style can be expanded in a straightforward way to include
other descriptors. For each of the types used, the potential file then
specifies the number used and any metaparameters. For $F_{n}$ this includes
the range of values for $n$ as well as the metaparameters $\alpha_{n}$,
$r_{c}$, $C_{max}$, $C_{min}$ and $r_{e}$. For $G_{m,k}$, $m_{tot}$ and
$k_{tot}$ are given along with the metaparameters $\beta_{k}$, $r_{e}$,
$r_{c}$, $C_{min}$, and $C_{max}$. Note that while we use the same cutoff
radius and angular screening parameters for both varieties, they can in
principle be different as every fingerprint type has its own set of
metaparameters. Next, the network architecture is specified by giving the
total number of layers and the length of each layer. Figure 2 shows a sample
header for the potential file. Following the header is the weight and bias
matrices for each layer and finally the activation functions used for each
layer. Using the potential file, the ANN subroutine then calculates the
feature list for every atom and uses these to calculate the energy and force
for every atom in the usual way. The subroutine and potential file also have
the capability of assigning multiple ANNs to different atoms types to be used
for multi-element simulations. In this way, the framework can be used for
arbitrary structural fingerprints as inputs to arbitrary MLP architectures for
systems with an arbitrary number of element types.
The subroutine is included as a normal pair-style in LAMMPS called ‘rann’ and
can be called in the usual way:

pair_style rann
pair_coeff * * potential_file.nn Element1 Element2 …

Figure 2: Sample header for the neural network potential style. The header
specifies the different element types as well as the style and number of
fingerprints used for each type. Metaparameters needed for each style are also
included as well as the overall network architecture.

##  4 Validation and Benchmarking

With the pair style implemented in LAMMPS, several potentials were fit for a
reference database for magnesium (Mg). The database, containing over 300,000
atomic environments was constructed using the open-source density functional
theory (DFT) software, Quantum Espresso [35]. The configurations used to train
the networks are shown in Table 1. These configurations included simple
triaxial strain perturbations of the low energy SC, FCC, BCC, and HCP phases
of Mg along with larger supercells with atoms displaced from equilibrium to
mimic the effect of thermal excitation. Free surface data and vacancy data,
with atoms also perturbed from equilibrium positions were also included to
improve stability. Effective temperatures up to 1000K were included in the
thermally excited database. Training of the network was carried out by
minimizing the RMSE of the system energies compared to the DFT results using
the Levenerg-Marquadt method [36, 37]. 10% of the training data from each
configuration type was left out of training and used for validation. The most
significant variation however, came in the size and construction of the
fingerprint used as the input layer. Along with the number and type of
fingerprints included, the most important factor in determining accuracy and
efficiency was the number of atoms included in the neighbor list. This is a
function of the cutoff radius $r_{c}$ and the angular screening parameter
$C_{min}$. Table 2 summarizes the different fingerprints considered and the
RMSE for both the training and validation datasets for each potential. As can
be seen, even the smallest architectures had RMSE of less than 3 meV/atom and
the best were less than 1 meV/atom. Table 2 also includes computational
performance for a sample calculation which will be described in detail below.

Table 1: DFT simulations used to generate training and validation data for the
Al potential. Training database attempts to consider all relevant low energy
structures which may be encountered during MD simulations.

Sample Description | Atoms per simulations | Number of Simulations | Total atomic environments
---|---|---|---
SC cubic cell w/ strains up to $\pm$15% | 1 | 500 | 500
FCC primitive cell w/ strains up to $\pm$15% | 1 | 500 | 500
BCC primitive cell w/ strains up to $\pm$15% | 1 | 500 | 500
HCP unit cell w/ strains up to $\pm$10% | 2 | 2700 | 5400
FCC 2x2x2 orthogonal supercell | 32 | 500 | 16000
BCC 3x3x3 orthogonal supercell | 54 | 500 | 27000
HCP 3x3x3 primitive supercell | 54 | 2000 | 108000
HCP 4x3x3 primitive supercell with vacancy | 67 | 100 | 6700
HCP 0001 free surface | 54 | 500 | 27000
HCP $10\overline{1}0$ free surface | 54 | 100 | 5400
Totals |  | 8100 | 307000

Table 2: Input Fingerprints considered in this work. $k_{tot}$ and $m_{tot}$ are the total number of values uses for $k$ and $m$ in Eq. 11. For the pair-style input, $n_{tot}$ was always 5. All ANNs contained a single hidden layer of 20 neurons. #SF is the total number of structural fingerprints for a particular potential and $r_{c}$ and $C_{m}in$ are the radial cutoff and angular screening parameter that determine the number of neighbors considered in calculation of the fingerprints. #NN is the number of neighbors included in a calculation of bulk HCP Mg at the equilibrium lattice constant. RMSE values are for the full training and validation databases and are given in meV/atom. Calculation speed is given in timesteps/second. The starred entry was used for validation tests. $k_{tot}$ | $m_{tot}$ | #SF | $r_{c}$ | $C_{min}$ | #NN | Training RMSE | Validation RMSE | Calculation speed
---|---|---|---|---|---|---|---|---
3 | 4 | 13 | 6 | 0.49 | 20 | 0.90 | 2.25 | 21.77
3 | 4 | 13 | 8 | 0.49 | 20 | 0.67 | 1.64 | 15.74
3 | 4 | 13 | 10 | 0.49 | 20 | 0.69 | 1.71 | 10.71
2 | 4 | 13 | 12 | 0.25 | 38 | 0.89 | 1.71 | 5.56
2 | 4 | 13 | 12 | 0.49 | 20 | 1.09 | 1.05 | 8.26
2 | 4 | 13 | 12 | 0.70 | 18 | 0.93 | 1.08 | 9.34
3 | 4 | 17 | 12 | 0.25 | 38 | 0.80 | 1.03 | 5.07
3 | 4 | 17 | 12 | 0.49 | 20 | 0.71 | 1.38 | 8.04
3 | 4 | 17 | 12 | 0.70 | 18 | 0.80 | 1.06 | 8.75
5 | 3 | 20 | 12 | 0.25 | 38 | 1.57 | 2.95 | 4.91
5 | 3 | 20 | 12 | 0.49 | 20 | 2.80 | 3.08 | 7.88
5 | 3 | 20 | 12 | 0.70 | 18 | 2.88 | 3.10 | 8.32
MEAM [33] |  |  |  |  |  |  |  | 61.963
N2P2 [27] |  |  |  |  |  |  |  | 5.12

###  4.1 Validation of the potential

Tests of the validity of the potential and its implementation in LAMMPS were
all performed using the 13 member input fingerprint with a single hidden layer
containing 20 neurons. The metaparameters for this potential are shown in
Table 3. $C_{min}$ was chosen exclude atoms beyond the third nearest neighbors
in the equilibrium hcp Mg structure. First, to confirm agreement with first
principles results beyond the RMSE in energy, the bulk properties of hcp
magnesium were tested for this potential. Table 4 shows a comparison of these
bulk properties among the sample ANN, the first principles database, and
experimental results. Agreement can be seen on the simple bulk properties
targeted.

Table 3: Metaparameters of the benchmark potential used for validation. $n_{tot}$ | 5
---|---
$k_{tot}$ | 3
$m_{tot}$ | 4
$\alpha$ | 5.52
$\beta_{k}$ | 1,2,9
$r_{c}$ | 6 $\AA$
$C_{min}$ | 0.49
$C_{max}$ | 2.90

As an initial test of the fidelity of the implementation, a number of simple
structures were generated and tested to confirm that the forces and energies
were correct and in agreement with one another. As a simple test, we displace
a single atom from its equilibrium location in a minimized HCP Mg lattice.
Figure 3 shows the force on the displaced atom as a function of displacement
as predicted by the ANN and by numerical differentiation of the energy.
Agreement is seen between the two measures demonstrating that forces are
correctly predicted by the algorithm.

Table 4: Properties of bulk magnesium as determined by experiment, the DFT database used for the fitting of the potential, and the 13-20-1 neural network used in this study. Property | Experiment | DFT | ANN
---|---|---|---
$a$ (nm) | 3.209 [38] | 3.194 | 3.196
$c$ (nm) | 5.211 [38] | 5.184 | 5.189
$c/a$ (nm) | 1.623 | 1.623 | 1.623
$E_{c}$ (eV) | 1.51 [38] | 1.45 | 1.46
$C_{11}$ (GPa) | 63.5 [39] | 61.6 | 63.4
$C_{12}$ (GPa) | 25.9 [39] | 23.8 | 24.3
$C_{13}$ (GPa) | 21.7 [39] | 21.2 | 18.5
$C_{33}$ (GPa) | 66.5 [39] | 64.3 | 62.9
$C_{44}$ (GPa) | 18.4 [39] | 17.3 | 18.4
Figure 3: The force as determined by the numerical derivative of the system
energy (solid red) and the calculated force (dashed black) as a single atom is
displaced in the ($\overline{1}\overline{1}20$) direction in a 4x4x4 hcp Mg
lattice. (inset) Difference between derivative and direct force calculation.

As a second test, a periodic fcc cubic unit cell was generated and allowed to
relax following a conjugate gradient minimization scheme. Symmetry is
maintained by the structures throughout minimization. The simulated cell is
then uniaxially compressed 10% in the x direction and slowly pulled in tension
to a total positive strain of 10% and the changes in energy and stress tensor
are examined. The results are shown in Figure 4. We see, again, that the
numerical derivative of the energy as a function of strain agrees exactly with
the stress calculated by the algorithm.
Having demonstrated that the pair style correctly reproduces the energy,
forces, and stresses in static configurations, we now confirm that dynamic
behavior can also be correctly produced. A dynamic simulation was carried out
with a 256 atom periodic hcp system (4x4x4 orthogonal unit cells). The system
was initialized from the equilibrium ground state with atoms given a random
velocity distribution producing an average temperature of 300K. The system was
then allowed to evolve at constant energy and volume. Figure 5 shows the
evolution of the kinetic, potential and total energy of the system. We see
that total energy is conserved.

Figure 4: The xx component of the stress tensor (solid red) and the derivative
of the system energy (dashed black) as uniaxial tension/compression is applied
in the x direction. (inset) Difference between derivative and stress
calculations. Figure 5: The potential, (red) kinetic, (black) and total (blue)
energy relative to the 0K minimum for a 4x4x4 hcp lattice of magnesium atoms
with initial velocity distribution producing a temperature of 300K. The system
was allowed to evolve under NVE conditions. Conservation of energy can be
observed for this system.

###  4.2 Computational Efficiency

While ANN potentials have the capability to reproduce experimental or first
principles results with greater accuracy than traditional formalisms, if the
computational performance suffers too greatly, the benefits become negligible.
The speed at which the ANN pair style will perform depends on the size and
complexity of the structural fingerprint, the number of neighbors within the
cutoff radius and the network architecture.
As a benchmark of the performance, we compare the performance of the Mg ANN
potential to an magnesium Modified Embedded Atom potential[33] and a recent
n2p2 style neural network potential [27]. To test the speed of the MEAM
potential and the various NN potentials, a sample system of hcp magnesium
containing 4000 atoms was minimized and each atom given a random velocity
matching a distribution with a temperature of 300K. The systems was then
allowed to evolve at constant energy and volume for 1000 timesteps. The number
of timesteps calculated per second is shown in Table 2. Calculations were
performed on one core of a Intel Xeon E5-2690 2.0GHz processor. We see that
all of the potentials are slower than MEAM, with a strong dependence on the
number of neighbors considered and on the cutoff radius $r_{c}$, but the
fastest potential, used in the validation study, runs approximately 1/3 as
fast as MEAM and four times faster than competitive n2p2 potentials. As can be
seen the first columns of Table 2, there is a strong dependence on the cutoff
radius even when the same number of atoms are included due to the angular
screening parameter. This is due to the computational time required to
calculate the angular screening needed to remove atoms from the neighbor list.
This is still seen to be more efficient than considering all the atoms within
the larger cutoff radius.
For the size of networks considered here, the computation time is dominated by
the calculation of the structural fingerprint with the propagation through the
network only a limiting factor for the smallest fingerprints and largest
architectures. This suggests that, assuming the structural fingerprint
contains sufficient information to describe the relevant configurations, it is
more efficient to increase the size of the network to increase the accuracy of
the potential, using this formalism, rather than expanding the fingerprint.
This is particularly true for extending the length of the cutoff radius. While
a large cutoff radius is often desirable and is used in existing ANN
potentials as it should improve the smoothness of the potential as atoms move
away from the target atom, the high computational cost of increasing the
neighbor list suggests that ideally only first, second, and possibly third
neighbors should normally be considered for metals. Angular screening combines
both of these aspects, limiting the number of neighbors in bulk settings while
still allowing relatively long-range interactions for studies involving voids,
free surfaces, or fracture.
The RANN formalism presented here has the advantage of requiring similar
computer time as MEAM to run simulations using the developed potentials and,
thanks to the additional of angular screening, faster than other existing ANN
formalisms, but the development time should be considerably smaller, as
fitting to an arbitrary number of targets is automated with many free
parameters available to tune to all available data.

##  5 Conclusions

We have presented an implementation of the Rapid Artificial Neural Network
(RANN) interatomic potential pair style for use in the LAMMPS molecular
dynamics code. The potential is suitable for static and dynamic calculations,
conserving total energy and correctly calculating the pressure and forces of
bulk solids. Such potential styles have been shown to be capable of
reproducing training data from first principles calculations at a level which
exceeds previous semi-empirical formalisms [6, 7, 10]. The developed pair
style is capable of accommodating network architectures of arbitrary dimension
and activation function, and the calculation of the structural fingerprint
requires only singly looped summations over the neighbor list of a target
atom, improving computational efficiency. Angular screening has also been
introduced both to improve efficiency by limiting the number of included
neighbors and to improve predictive power by introducing the physically
motivated phenomenon of shielding and smoothing fingerprints as atoms move
through the radial cutoff. Computational efficiency is of particular
importance for this formalism as the improvements in accuracy over existing
formalisms such as MEAM are greatly diminished if the runtime is closer to
that of first principles calculations. As such, we demonstrate that for a
network with 13 input fingerprints and a cutoff radius of 6.0 $\AA$ for the
neighbor list, and angular screening restricting interactions to third nearest
neighbors in bulk HCP the implementation performs at one third the speed of
MEAM. Larger networks, fingerprints, or neighbor lists will hamper this
performance, with the most significant cost associated with increasing the
neighbor list through the length of the radial cutoff. Ultimately, a flexible,
scalable formalism for ANNs, which can utilize the native parallel processing
available in LAMMPS has been demonstrated with a formalism which produces high
accuracy potentials which operate on speeds comparable to those of MEAM.

## Acknowledgements

## References

## References

  * [1] Hongxiang Zong, Ghanshyam Pilania, Xiangdong Ding, Graeme J Ackland, and Turab Lookman.  Developing an interatomic potential for martensitic phase transformations in zirconium by machine learning.  npj Computational Materials, 4(1):48, 2018.
  * [2] N. Artrith and A. Urban.  An implementation of artificial neural-network potentials for atomistic materials simulations: Performance for tio2.  Computational Materials Science, 114:135–150, 2016.
  * [3] Gabriele C Sosso, Giacomo Miceli, Sebastiano Caravati, Jörg Behler, and Marco Bernasconi.  Neural network interatomic potential for the phase change material gete.  Physical Review B, 85(17):174103, 2012.
  * [4] Jörg Behler.  Representing potential energy surfaces by high-dimensional neural network potentials.  Journal of Physics: Condensed Matter, 26(18):183001, 2014.
  * [5] Daniele Dragoni, Thomas D Daff, Gábor Csányi, and Nicola Marzari.  Achieving dft accuracy with a machine-learning interatomic potential: Thermomechanics and defects in bcc ferromagnetic iron.  Physical Review Materials, 2(1):013808, 2018.
  * [6] Atsuto Seko, Akira Takahashi, and Isao Tanaka.  First-principles interatomic potentials for ten elemental metals via compressed sensing.  Physical Review B, 92(5):054113, 2015.
  * [7] Tran Doan Huan, Rohit Batra, James Chapman, Sridevi Krishnan, Lihua Chen, and Rampi Ramprasad.  A universal strategy for the creation of machine learning-based atomistic force fields.  NPJ Computational Materials, 3(1):37, 2017.
  * [8] Volker L Deringer and Gábor Csányi.  Machine learning based interatomic potential for amorphous carbon.  Physical Review B, 95(9):094203, 2017.
  * [9] Albert P Bartók, James Kermode, Noam Bernstein, and Gábor Csányi.  Machine learning a general-purpose interatomic potential for silicon.  Physical Review X, 8(4):041048, 2018.
  * [10] Ryo Kobayashi, Daniele Giofré, Till Junge, Michele Ceriotti, and William A Curtin.  Neural network potential for al-mg-si alloys.  Physical Review Materials, 1(5):053604, 2017.
  * [11] A. Takahashi, A. Seko, and I. Tanaka.  Conceptual and practical bases for the high accuracy of machine learning interatomic potentials: Application to elemental titanium.  Physical Review Materials, 1:063801, 2017.
  * [12] G. P. P. Pun and Y. Mishin.  Optimized interatomic potential for silicon and its application to thermal stability of silicene.  Physical Review B, 95(22):224103, 2017.
  * [13] Aidan P Thompson, Laura P Swiler, Christian R Trott, Stephen M Foiles, and Garritt J Tucker.  Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials.  Journal of Computational Physics, 285:316–330, 2015.
  * [14] Jörg Behler.  Atom-centered symmetry functions for constructing high-dimensional neural network potentials.  The Journal of Chemical Physics, 134:074106, 2011.
  * [15] Jörg Behler and Michele Parrinello.  Generalized neural-network representation of high-dimensional potential-energy surfaces.  Physical review letters, 98(14):146401, 2007.
  * [16] Albert P Bartók, Mike C Payne, Risi Kondor, and Gábor Csányi.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  Physical review letters, 104(13):136403, 2010.
  * [17] Wojciech J Szlachta, Albert P Bartók, and Gábor Csányi.  Accuracy and transferability of gaussian approximation potential models for tungsten.  Physical Review B, 90(10):104108, 2014.
  * [18] So Fujikake, Volker L Deringer, Tae Hoon Lee, Marcin Krynski, Stephen R Elliott, and Gábor Csányi.  Gaussian approximation potential modeling of lithium intercalation in carbon nanostructures.  The Journal of chemical physics, 148(24):241714, 2018.
  * [19] D. Dragoni, T. D. Daff, G. Csányi, and N. Marzari.  Achieving dft accuracy with a machine-learning interatomic potential: Thermomechanics and defects in bcc ferromagnetic iron.  Physical Review Materials, 2(1):013808, 1991.
  * [20] Thuong T Nguyen, Eszter Székely, Giulio Imbalzano, Jörg Behler, Gábor Csányi, Michele Ceriotti, Andreas W Götz, and Francesco Paesani.  Comparison of permutationally invariant polynomials, neural networks, and gaussian approximation potentials in representing water interactions through many-body expansions.  The Journal of chemical physics, 148(24):241725, 2018.
  * [21] Nongnuch Artrith, Alexander Urban, and Gerbrand Ceder.  Efficient and accurate machine-learning interpolation of atomic energies in compositions with many species.  Physical Review B, 96(1):014112, 2017.
  * [22] M. I. Baskes.  Modified embedded-atom potentials for cubic materials and impurities.  Physical Review B, 46(5):2727, 1992.
  * [23] D Dickel, DK Francis, and CD Barrett.  Neural network aided development of a semi-empirical interatomic potential for titanium.  Computational Materials Science, 171:109157, 2020.
  * [24] Mashroor S Nitol, Doyl E Dickel, and Christopher D Barrett.  Artificial neural network potential for pure zinc.  Computational Materials Science, 188:110207, 2021.
  * [25] Andreas Singraber, Jörg Behler, and Christoph Dellago.  Library-based lammps implementation of high-dimensional neural network potentials.  Journal of chemical theory and computation, 15(3):1827–1840, 2019\.
  * [26] Daniel Marchand, Abhinav Jain, Albert Glensk, and WA Curtin.  Machine learning for metallurgy i. a neural-network potential for al-cu.  Physical Review Materials, 4(10):103601, 2020.
  * [27] Markus Stricker, Binglun Yin, Eleanor Mak, and WA Curtin.  Machine learning for metallurgy ii. a neural-network potential for magnesium.  Physical Review Materials, 4(10):103602, 2020.
  * [28] M. S. Daw and Baskes M. I.  Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals.  Physical Review B, 29:6443, 1984.
  * [29] S. J. Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  Journal of Computational Physics, 117(1):1–19, 1995.
  * [30] Jörg Behler.  Perspective: Machine learning potentials for atomistic simulations.  The Journal of chemical physics, 145(17):170901, 2016.
  * [31] M. I. Baskes.  Determination of modified embedded atom method parameters for nickel.  Mater. Chem. Phys., 50:152, 1997.
  * [32] B.-J. Lee and M. I. Baskes.  Second nearest-neighbor modified embedded atom method potential.  Physical Review B, 62:8564, 2000.
  * [33] Z Wu, MF Francis, and WA Curtin.  Magnesium interatomic potential for simulating plasticity and fracture phenomena.  Modelling and Simulation in Materials Science and Engineering, 23(1):015004, 2015.
  * [34] Michael I Baskes.  Determination of modified embedded atom method parameters for nickel.  Materials Chemistry and Physics, 50(2):152–158, 1997.
  * [35] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, et al.  Quantum espresso: a modular and open-source software project for quantum simulations of materials.  Journal of physics: Condensed matter, 21(39):395502, 2009.
  * [36] Kenneth Levenberg.  A method for the solution of certain non-linear problems in least squares.  Quarterly of Applied Mathematics, 2(2):164–168, jul 1944.
  * [37] Donald W. Marquardt.  An algorithm for least-squares estimation of nonlinear parameters.  Journal of the Society for Industrial and Applied Mathematics, 11(2):431–441, jun 1963.
  * [38] N Saunders, AP Miodownik, and AT Dinsdale.  Metastable lattice stabilities for the elements.  Calphad, 12(4):351–374, 1988.
  * [39] G Simmons and H Wang.  Single crystal elastic constants and calculated aggregate properties: A handbook 2nd ed., 370, 1971.