# The Relation Between AGN and Host Galaxy Properties:
I. Obscured AGN reside in disturbed hosts at $0<z<4$

Nina Bonaventura Steward Observatory, University of Arizona, 933 North Cherry
Avenue, Tucson, AZ 85721, USA  Jianwei Lyu (吕建伟) Steward Observatory,
University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA  George
H. Rieke Steward Observatory, University of Arizona, 933 North Cherry Avenue,
Tucson, AZ 85721, USA  Stacey Alberts Steward Observatory, University of
Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA  Christopher N. A.
Willmer Steward Observatory, University of Arizona, 933 North Cherry Avenue,
Tucson, AZ 85721, USA  Pablo G. Pérez-González Centro de Astrobiología (CAB),
CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain
Andrew J. Bunker Department of Physics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford OX1 3RH, UK  Meredith Stone Steward Observatory,
University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA
Francesco D’Eugenio Kavli Institute for Cosmology, University of Cambridge,
Madingley Road, Cambridge, CB3 0HA, UK  Cavendish Laboratory, University of
Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK  Christina C. Williams
NSF’s National Optical-Infrared Astronomy Research Laboratory, 950 North
Cherry Avenue, Tucson, AZ 85719, USA  Steward Observatory, University of
Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA  Michael V. Maseda
Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St.,
Madison, WI 53706 USA  Chris J. Willott NRC Herzberg, 5071 West Saanich Rd,
Victoria, BC V9E 2E7, Canada  Zhiyuan Ji Steward Observatory, University of
Arizona, 933 N. Cherry Ave, Tucson, AZ 85721, USA  William M. Baker Kavli
Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK  Cavendish Laboratory, University of Cambridge, 19 JJ Thomson
Avenue, Cambridge CB3 0HE, UK  Stefano Carniani Scuola Normale Superiore,
Piazza dei Cavalieri 7, I-56126 Pisa, Italy  Stephane Charlot Sorbonne
Université, UPMC-CNRS, UMR7095, Institut d’Astrophysique de Paris, F-75014
Paris, France  Jacopo Chevallard Department of Physics, University of Oxford,
Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK  Emma Curtis-Lake
Centre for Astrophysics Research, Department of Physics, Astronomy and
Mathematics, University of Hertfordshire, Hatfield AL10 9AB, UK  Daniel J.
Eisenstein Center for Astrophysics — Harvard and Smithsonian, 60 Garden
Street, Cambridge, MA 02138, USA  Kevin Hainline Steward Observatory,
University of Arizona, 933 N. Cherry Ave, Tucson, AZ 85721, USA  Ryan Hausen
Department of Physics and Astronomy, The Johns Hopkins University, 3400 N.
Charles St., Baltimore, MD 21218, USA  Erica J. Nelson Department for
Astrophysical and Planetary Science, University of Colorado, Boulder, CO
80309, USA  Marcia J. Rieke Steward Observatory, University of Arizona, 933 N.
Cherry Ave, Tucson, AZ 85721, USA  Brant Robertson Department of Astronomy and
Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa
Cruz, CA 95064, USA  Irene Shivaei Centro de Astrobiología (CAB), CSIC-INTA,
Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain

###### Abstract

The morphology of a galaxy is a manifestation of the complex interplay of
physical processes occurring within and around it, and therefore offers
indirect clues to its formation and evolution. We use both visual
classification and computer vision to verify the suspected connection between
galaxy merging activity — as evidenced by a close/merging galaxy pair, or
tidal features surrounding an apparently singular system — and AGN activity.
This study makes use of JADES JWST/NIRCam imagery, along with an
unprecedentedly complete sample of AGN built using JWST/MIRI photometry in the
same field. This 0.9-25 $\mu$m dataset enables constraints on the host galaxy
morphologies of the broadest possible range of AGN beyond z$\sim$1, including
heavily obscured examples missing from previous studies. We consider two AGN
samples, one consisting of lightly to highly obscured X-ray-selected AGN (Lyu
et al., 2022), and the other, presumed Compton-thick mid-infrared-
bright/X-ray-faint AGN recently revealed by MIRI (Lyu et al., 2023). Both
samples contain a significant fraction of host galaxies with disturbed
morphologies at all redshifts sampled, and increasingly so towards higher
redshift and AGN bolometric luminosity. The most obscured systems show the
highest fraction of strongly disturbed host galaxies at $95\pm 4\%$, followed
by the moderately and unobscured/lightly obscured subsets at $78\pm 6\%$ and
$63\pm 6.5\%$, respectively. From this pattern of disturbances, we conclude
that mergers are common amongst obscured AGN, and that the obscured AGN phase
may mark a period of significant SMBH growth. This finding presents tension
with the leading model on AGN fueling mechanisms (Hopkins et al., 2014) that
needs reconciling.

##  1 Introduction

In the hierarchical framework for large-scale structure formation in the
Universe, a ‘cosmic cycle’ for galaxy formation and evolution is expected as a
direct result of regularly occurring galaxy mergers (Hopkins & Hernquist,
2006). In this picture, two gas-rich disk galaxies of similar masses merge to
violently drive a nuclear inflow of gas towards the center of the system that
triggers both starbursts and quasars ($L_{bol}\geq 10^{45}$ ergs/s), and
therefore the simultaneous growth of black holes and stellar bulges in the
centers of galaxies. Observational studies support this modeled co-evolution
of galaxies with their central supermassive black holes (SMBHs), most famously
the empirical $M-\sigma_{disp}$ relation showing the tight correlation between
black hole mass and host galaxy velocity dispersion, demonstrating an
undeniable link between the origin of galaxies and SMBHs (Ferrarese & Merritt,
2000; Gebhardt et al., 2000; Tremaine et al., 2002). There is also the
persistent coincidence of intensive star formation and embedded active nuclei
in ultra-luminous infrared galaxies (ULIRGs) commonly associated with mergers
(Sanders & Mirabel, 1996; Jogee, 2006).

It is predicted that for most of the lifetime of a merger, quasar activity is
obscured at optical and (soft) X-ray wavelengths by a screen of gas and dust,
manifesting in the appearance of an obscured, Type II quasar. The obscuring
material could be confined to the nuclear region of the galaxy, and/or re-
distributed to galaxy-wide scales by the triggering merger event (Hickox &
Alexander, 2018). This early obscured phase would then transition into the
relatively brief, luminous, unobscured Type I quasar state once the radiative
power of the black hole increases to the point of ejecting the obscuring
material out and away, halting the ability of the active nucleus to sustain
itself. At this point, a merger remnant would emerge in the form of a normal,
passively evolving galaxy with a central stellar bulge and SMBH, with or
without a surviving (or reformed) disk structure (Hopkins et al., 2009). As
this process propagates across the galaxy population throughout cosmic time
(and sometimes multiple times within the lifetime of a given galaxy/merger
remnant), there is the expectation of a shift towards merger remnants with
lower masses towards lower redshifts given the diminishing gas supply in the
Universe; and, ultimately, to gas-poor elliptical galaxies (Hopkins &
Hernquist, 2006; Hickox & Alexander, 2018; Hopkins et al., 2008; Blecha et
al., 2018)).

The implication of this theoretical evolutionary scenario is that major galaxy
mergers are largely responsible for the bulk of SMBH growth in the Universe
(Hickox & Alexander, 2018; Hopkins et al., 2008; Gilli et al., 2007); the
alternative is that they are predominantly fueled by stochastic processes,
slowly over time, instead of through a series of short and dramatic episodes
(Kormendy & Kennicutt, 2004; Hopkins & Hernquist, 2006; Hopkins et al., 2014).
However, the observational evidence for the former conclusion is ambiguous
(Treister et al., 2012; Zhiyuan et. al., 2022) and the issue remains
unresolved (see, e.g., Shah et al., 2020; Pierce et al., 2023; Hernández-
Toledo et al., 2023; Villforth, 2023). A further observational complexity is
that Type I quasar host galaxies tend to lie on the star-forming fundamental
plane (e.g., Xu et al., 2015; Zhang et al., 2016), rather than emerging from
ULIRGs or lying in post-starburst, passively evolving galaxies, as might be
expected from theoretical predictions. Finally, it is not clear if this
generally accepted model for quasar evolution extends into the realm of
Seyfert galaxies, defined by significantly lower luminosities and black hole
masses than quasars (Hopkins et al., 2014).

Clues to this evolutionary scenario should lie in the host galaxy
morphologies. There have been many programs to image the host galaxies of
nearby and moderate redshift AGN and quasars; results are available for well
over 200 systems (Kim et al., 2017). Although there are disagreements in the
assignment of disturbances in individual cases, in general there is broad
consistency among the studies that roughly 50% of the hosts appear disturbed,
as illustrated in the comprehensive source-by-source summary by Kim et al.
(2017). This agreement holds when different techniques are used, providing
credibility to the assignments.

The great majority of the available studies have focused largely on Type I
unobscured or lightly obscured AGN. However, the discovery of substantial
numbers of Type II AGN in the SDSS (Zakamska et al., 2004) has enabled similar
observations of their host galaxies. For example, Urbano-Mayorgas et al.
(2019) report that, in a sample of 41 Type II quasars, $34^{+6}_{-9}$% of the
host galaxies are highly disturbed; they summarize the existing work as
showing that $44-62$% of Type II hosts are disturbed. Zakamska et al. (2019)
present HST images of 10 extremely red quasars (ERQs)111Zakamska et al. (2019)
suggest that the relatively low number of disturbed hosts detected amongst the
ERQs of their study may be due to the undetectability of faint merger
signatures in the HST imagery utilized, e.g. if these systems are late-stage
merger remnants and/or lie in very dusty host galaxies. and 6 Type II AGN,
finding, respectively, that ($2\pm 1$)/10 and $2/6$ of the hosts are
disturbed. That is, it appears that these obscured AGN also show a significant
fraction of disturbed hosts, as classified by the available methods. However,
even the Type II samples discussed above do not include the most heavily
obscured, Compton-thick ($N_{H}\geq 10^{24}$ $cm^{-2}$) AGN. A comprehensive
study of the relation of host morphology to AGN obscuration needs to draw on
the full range of AGN identification methods (Lyu et al., 2022), including
those centered in the infrared that can detect AGN that are so obscured that
their traces are hidden in other wavelength ranges. Earlier studies with HST
and the _Spitzer_ Infrared Camera (IRAC) attempted to compare the
morphological properties of AGN selected by different methods (Zhiyuan et.
al., 2022), although those studies were only able to focus on AGN much
brighter than those found by Lyu et al. (2023).

The possibility with the JWST Near-Infrared Camera (NIRCam) to study AGN host
galaxy morphology at cosmic noon with unprecedented resolution and sensitivity
at rest-optical wavelengths, provides the motivation for a new approach to
these issues. Furthermore, the thoroughness of the AGN identification in the
GOODS-South (GOODS-S) field, including the identification of a large sample of
previously unknown obscured AGN out to high redshift (Lyu et al., 2022, 2023)
allows for the study of the connection between large-scale morphological
disturbances in an active galaxy and evidence of obscuration by dust and/or
gas of its accretion disk. This paper initiates such a study. Here (Paper I),
we present the details of our analysis of the JWST observations. In Paper II,
we will supplement Paper I with additional multiwavelength data to support an
evolutionary scenario for massive Seyfert galaxies implied by the results of
the initial study.

Figure 1: JADES V1 NIRCam F150W mosaic image containing the AGN sample
targeted in GOODS-S/HUDF. Blue points mark the X-ray-selected AGN sample from
Lyu et al. (2022), and the magenta points the mid-IR-selected MIRI AGN from
Lyu et al. (2023)

Our discussion will develop as follows. Section 2 provides an overview of the
new data we have used, specifically JWST/NIRCam images (Rieke et al., 2023) of
both a pre-JWST X-ray-selected AGN sample (Lyu et al., 2022), and a JWST Mid-
Infrared Instrument (MIRI) mid-IR-selected AGN sample. We also summarize the
AGN identification methods relevant to our study. Section 3 describes the
experimental design of our analysis, consisting of two distinct trials. Trial
1 takes a traditional approach to the analysis of the X-ray-selected subset of
the AGN sample, utilizing the hydrogen absorbing column density, $N_{H}$, as
the measure of AGN obscuration. We evaluate the host galaxy morphology for
unobscured/light ($N_{H}<10^{22}$ $cm^{-2}$), moderate
($10^{22}<N_{H}<10^{23}$ $cm^{-2}$), and heavy X-ray obscuration
($N_{H}>10^{23}$ $cm^{-2}$), and place the results in the context of existing
studies. Given that the use of X-rays as a selection criterion biases the
sample against the more heavily obscured objects (Del Moro et al., 2016; Lyu
et al., 2022), Trial 2 analyzes the behavior of the host galaxies from the
mid-IR-selected sample that should include a better and less-biased
representation of the behavior of the most-obscured (presumed Compton-thick)
AGN. Finally, we present the results of our analysis trials in Section 4, and
a discussion of those results in Section 5.

The galaxy parameter values quoted in this paper are drawn directly from the
SED analyses of Lyu et al. (2022) and Lyu et al. (2023), which adopt a slighty
different cosmology: Lyu et al. (2022) assumes $\Omega_{m}=0.27$ and
$H_{0}=71$ km s-1Mpc-1, while (Lyu et al., 2023) assumes $\Omega_{m}=0.287$
and $H_{0}=69.3$ km s-1Mpc-1.

Figure 2: Histograms showing the distribution of the bolometric luminosities
($L_{bol}$), redshifts (_z_), and stellar masses ($M_{*}$) of the X-ray and
MIRI AGN samples (Lyu et al., 2022, 2023); and the hydrogen absorbing column
densities ($N_{H}$) used to measure obscuration in the X-ray AGN sample (Lyu
et al., 2022). The $L_{bol}$ values for the X-ray AGN sample represent the
K-corrected intrinsic X-ray luminosity ($L_{x,int}$) (Lyu et al., 2022; Duras
et al., 2020); and those for the MIRI AGN sample are based on the intrinsic,
unobscured $4400$ Å luminosity derived from the SED fitting procedure
described in Lyu et al. (2023).

##  2 Sample and Data

This study is enabled by previous work to identify AGN in the GOODS-S region
using ultra-deep _Chandra_ and _Spitzer_ imaging, now complemented by both the
mid-IR identification of obscured AGN in the same region with JWST/MIRI, plus
the ability to examine the host galaxies uniformly using JWST/NIRCam.

###  2.1 AGN sample

We gather sources from two data sets: the comprehensive census of GOODS-S/HUDF
AGN detected at X-ray, optical, infrared, and/or radio wavelengths analyzed in
Lyu et al. (2022) over an area of $\sim$ 170 arcmin2; and a new sample of
JWST/MIRI-detected AGN revealed by the SMILES MIRI GTO program (PI: G. Rieke)
covering $\sim$ 34 arcmin2 Lyu et al. (2023). Both Type I and II AGN are
included in these samples, detected out to $z\sim 6$ via spectroscopic and/or
photometric measurements, and are presented with a full set of SED-derived
physical parameters.

####  2.1.1 X-ray AGN

The X-ray AGN sample was defined as all sources in the catalog of Lyu et al.
(2022) for which the ratio of the 0.5-7 keV X-ray luminosity to the 3 GHz
radio luminosity exceeds the physical threshold for stellar processes in a
galaxy, and is therefore confidently attributed to AGN emission (Alberts et
al., 2020). In fact, Lyu et al. (2022) discovered that the selection of AGN
using the X-ray-to-radio luminosity ratio was the most effective of the
multiwavelength selection methods they considered, until the recent multi-band
JWST/MIRI measurements became available. We also restrict our sample to those
X-ray-bright AGN with an available SED-derived $N_{H}$ measurement from Luo et
al. (2017). The vast majority of the resulting sample (hereafter the ‘pre-
JWST’ AGN sample of Lyu et al., 2022) are optically faint, which conveniently
minimizes AGN contamination in the rest-optical NIRCam imagery we utilize for
morphology analysis. The majority likewise lack a mid-IR counterpart in the
SED analysis of Lyu et al. (2022), which appears to confirm the mostly faint
intrinsic X-ray luminosities inferred from their SEDs, as photometric
detection in the mid-IR requires a stronger intrinsic AGN luminosity than in
the X-ray selection due to a higher competition from the host galaxy light in
this wavelength regime of the SED (Hickox & Alexander, 2018).

The properties of the X-ray AGN sample are shown in Figure 2. Of the total of
169 X-ray AGN meeting the criteria above, 151 were covered by the JADES/NIRCam
mosaic utilized in our analysis.

####  2.1.2 Mid-Infrared AGN

The mid-IR AGN sample constitutes the X-ray-faint subset of a batch of newly
identified, JWST/MIRI-selected AGN from SMILES, a JWST Cycle 1 GTO program
that has targeted the central region of GOODS-S with eight MIRI filters from
5.6 to 26 $\mu$m (Lyu et al., 2023). These sources were classified as AGN
through an SED analysis consisting of JWST MIRI and NIRCam photometry, in
addition to HST photometry (Lyu et al., 2023). As we require the selected
subset of sources to lack a confident identification as X-ray AGN using the
selection methods discussed in Lyu et al. (2022, 2023), these mid-IR sources
are expected to represent the most heavily dust- and gas-obscured AGN in the
sample.

The properties of the mid-IR AGN sample are shown in Figure 2 alongside those
of the X-ray sample. Of the total of 182 mid-IR AGN meeting the criteria
above, 138 were covered by the JADES/NIRCam mosaic utilized in our analysis.

###  2.2 JWST NIRCam Data

The NIRCam imaging utilized in this study was obtained as part of the JADES
survey in JWST program 1180 (PI: Eisenstein), described in full detail in the
JADES Data Release v1 paper (Rieke et al., 2023), as well as the JWST
Extragalactic Medium-band Survey (JEMS) (JWST PID 1963) (Williams et al.,
2023). JADES v1 NIRCam mosaic images of the GOODS-S/HUDF field are available
in the $0.9-4.4$ $\mu$m range in each of the following wide and medium
filters: F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, and F444W.
Each mosaic image achieves $0.03^{\prime\prime}$/pixel resolution and an AB
magnitude depth of 29-30 (Rieke et al., 2023). We chose to work exclusively
with the F150W NIRCam mosaic image for this study, as it proved to be the
deepest and most resolved of the mosaics that also maximizes coverage of the
sky positions of the AGN in our sample (see the Appendix for further details).
It is displayed in Figure 1, showing the overlapping coverage of the AGN
sample.

The NIRCam F150W filter has significant overlap with the HST WFC3 F160W band
that has been widely utilized in various types of galaxy studies, from probing
the peak of galaxy stellar emission at z$\sim$0 to assisting in the capture of
‘Lyman break’ galaxies in the SEDs of $z>9$ ‘ _J_ -band dropout’ galaxies
(see, e.g., Henry et al., 2008). The observed effective wavelength of 1.5
$\mu$m of this filter translates to the rest-optical to rest-near-IR
wavelength range of 0.33-1.5 $\mu$m for our $0<\textit{z}<4$ AGN sample, which
primarily probes stellar emission in the AGN host galaxies, at $\sim$5x
greater sensitivity than the HST F160W imaging band. While dust attenuation of
starlight has been shown to bias measures of galaxy sizes and morphologies
(see, e.g., Suess et al., 2022), we do not expect this effect to significantly
impact the results and conclusions of the present study. The reason for this
lies in our consideration of the non-parametric _shape_ asymmetry of the AGN
host galaxies as measured in the binary detection mask (also known as a
‘segmentation map’) associated with each galaxy image, composed of equally-
weighted pixels (see Section 3.1.2 for details). This contrasts with the more
commonly used, classic asymmetry parameter that also measures variations in
the brightness distribution of the pixel values directly drawn from the galaxy
image (as well as parametric fits to galaxy surface brightness profiles to
determine galaxy sizes, morphologies, and sub-structures). More specifically,
given that the shape asymmetry parameter effectively traces the outline of a
galaxy’s shape in a binary image, without regard to differences in pixel
intensity, it would only be influenced by the effects of dust attenuation if
the outer edges of an AGN host galaxy were significantly dusty, and the
starlight significantly absorbed, fully or partially hiding the full extent
and shape of the host galaxy.

The choice to analyze the AGN host galaxies using only image cutouts from the
observed F150W JADES mosaic – which probes a range of rest-frame wavelengths
as opposed to a single consistent rest-frame wavelength – was motivated by the
following additional factors: 1) the limited overlap of our AGN sample with
the mosaics exposed at the other NIRCam wavelengths sampled by JADES; 2) the
obvious resolution degradation of the galaxy images at wavelengths longward of
F150W due to the increasing PSF size222The observed degradation in image
resolution for a given AGN host galaxy in our sample, as highlighted in the
Appendix, when attempting to probe the same rest-frame wavelength by using
mosaic images exposed at increasingly longer observed wavelengths, led us to
conclude that no information was lost, or final conclusions altered, by
foregoing their inclusion in our analysis. In other words, it was clear that
by utilizing only the F150W-band mosaic, that we were always choosing the best
imagery at hand to constrain galaxy morphology.; and 3) the desire to avoid
the extra sources of uncertainty that would be introduced to our analysis by
comparing galaxy images extracted from mosaics in different filters.

##  3 Methodology

Methods to constrain galaxy morphologies historically range from the manual,
visual classification of galaxy shapes that led most famously to Edwin
Hubble’s ‘tuning fork’ sequence of galaxy classes still in use today, to the
equally famous and successful Galaxy Zoo project that employs citizen
scientists around the world to visually classify galaxy images (Lintott et
al., 2008). Now, with continued advances in telescope technology that will
yield sufficiently large data samples to formally qualify as ‘big data’, most
notably represented by the recently launched Euclid mission and upcoming LSST
survey, sophisticated machine learning computer algorithms must increasingly
be relied upon for such a task. A key advantage, however, of a classical
visual morphological analysis of a humanly manageable galaxy sample, is that
it is robust against the inevitable sources of uncertainty introduced by such
automated computer algorithms. These algorithms are sufficient for
constraining generalized properties, but often fail to consider the nuanced
features of individual galaxies that must not be ignored if one is to hope to
achieve a complete understanding of the physical phenomena leading to their
structures (Reza et al., 2021; Walmsley et al., 2002). Furthermore, the
requirements for machine learning and citizen science, i.e., very large data
samples and/or reliable training sets, are not yet available for the data we
use in this work.

###  3.1 AGN Host Galaxy Morphology Analysis

In the present study, we exploit the robustness of the visual classification
method to identify merger signatures in images of 289 AGN host galaxies in the
form of a galaxy-wide disturbance or asymmetry (see Bignone et al., 2017 and
references therein). The method chosen closely adheres to the visual
classification procedure defined by Kocevski et al. (2015) (which represents a
refinement to the set of visual classification rules initially defined in
Kocevski et al., 2012). We then follow this procedure with a quantitative
computer visual analysis in the form of an automated non-parametric measure of
the asymmetry in each galaxy’s spatial pixel distribution using the
_statmorph_ software (Rodriguez-Gomez et al., 2019).

Figure 3: Examples of figure output from statmorph showing strongly
asymmetric/disturbed (top, $As>0.2$), mildly asymmetric/disturbed (middle,
$0.1<As<0.2$), and symmetric/undisturbed (bottom, $As<0.1$). The galaxy image
is shown along the left column, and along the right column, the binary
detection mask (segmentation map) used to calculate the shape asymmetry
($A_{s}$) parameter utilized in this study. The measurement is performed
within the circular aperture set by $r_{max}$, defined by Rodriguez-Gomez et
al. (2019) as the distance between the point that minimizes the standard
asymmetry (A) and the most distant pixel that belongs to the detection mask
(which differs slightly from the $r_{max}$ value originally defined in Pawlik
et al., 2016, who introduce $A_{s}$). The elliptical Petrosian aperture,
centered on the point in the galaxy image that minimizes A, is used to
calculate the sky background level in the galaxy image that sets the
brightness threshold for the binary detection mask.

####  3.1.1 Visual Classification

Kocevski et al. (2015) present a straightforward procedure for classifying the
optical appearance of a galaxy into the following categories: undisturbed
disk, disturbed disk, spheroid, point source, irregular, or galaxy merger. To
streamline the results of our analysis, we classify as “disturbed” all
galaxies identified as irregularly shaped; an obvious merger, such as a close
pair with two visible nuclei, or a classic ‘trainwreck’ appearance; as well as
an asymmetric disk. Likewise, all point sources, spheroids, and apparently
symmetric disks were classified as “undisturbed.” We also chose to refine the
“disturbed disk” category into “strongly” and “mildly” disturbed
subcategories, as we observed significant numbers of both types in our trials.
Finally, we classified the small number of apparently fully edge-on disk
galaxies as “undisturbed”, under the presumption that a galaxy-wide
disturbance would manifest itself in three dimensions and therefore be
observable from all disk inclination angles relative to the observer.

Figure 4: The results of the visual classification of the morphologies of the
289 AGN host galaxies in our sample, showing the percentage that are
undisturbed (‘point source’, ‘spheroid’ , and ‘disk’ ) and disturbed
(‘merger’, ‘irregular’, and ‘disturbed disk’). Upper left: Least obscured
X-ray AGN ($N_{H}<10^{22}cm^{-2}$). Upper-right: Moderately obscured X-ray AGN
($10^{22}<N_{H}<10^{23}cm^{-2}$). Lower-left: Heavily obscured X-ray AGN
($N_{H}>10^{23}cm^{-2}$. Note that the orange slice in this plot that
represents ‘unresolved asymmetric disk’ ($17.5\%$) marks cases where the
galaxy image resolution was sufficient to visually detect an asymmetry, but
not fine enough (or the galaxy size large enough) to distinguish between a
mild and strong disturbance. Lower-right: presumed Compton-thick mid-IR-
bright/X-ray-faint AGN.

####  3.1.2 Computer Vision

While visual classification tasks performed by a human have proven to be very
accurate for many types of analyses (Brinchmann et al., 1998; Bundy et al.,
2005; Kampczyk et al., 2007; Darg et al., 2010; Treister et al., 2012;
Kartaltepe et al., 2015; Walmsley et al., 2002), by definition they are
subject to ‘human error’ on account of oftentimes uncontrollable factors such
as imperfect vision, and bias and differences in reasoning. In the context of
our study, such error could have occurred in instances of modest image
resolution or signal-to-noise, hidden evidence of weak galaxy disturbances
such as faint tidal features; or if the strong inclination angle of a disk
galaxy limited the view of mild spatial disturbances. Another important
consideration in this context is the nonlinear response of the human eye to
light intensity, that could have led to a mis-classification of truly
disturbed galaxies as symmetric/undisturbed, e.g. in cases of high brightness
contrast that would hide faint asymmetries/disturbances. Therefore, to
systematically quantify the results of our visual analysis and catch such
cases of subtle host disturbances that may have gone unnoticed or mis-
categorized by the human eye, we utilized the _statmorph_ software package to
perform an independent, computer visual check of the results. This software
package was specifically chosen for its proven ability to detect faint merger
signatures around the edges of galaxies (Rodriguez-Gomez et al., 2019).

Figure 5: A comparison of the results of the _statmorph_ $A_{s}\geq 0.2$ (dark
colors) and visual (light colors, width artificially narrowed and offset for
clarity) classification methods, showing the total percentage of AGN host
galaxies presenting as strongly asymmetric/disturbed by the two independent
analyses. It can be seen that the computer visual analysis identifies all of
the strongly disturbed cases identified by the visual analysis, plus
additional ones that evaded detection by the human eye. (Given that all/most
of the _statmorph_ $0.1<A_{s}<0.2$ mild asymmetry identifications were
classified as undisturbed to the human eye, we show only the strongly
disturbed _statmorph_ results for the most conservative comparison to the
visual classification results.) Bin placement and width (dark colors)
correspond to the mean and standard deviation of the $N_{H}$ bins for the
X-ray test and control sub-samples; while the bin for the mid-IR-bright/X-ray-
faint MIRI AGN does not reflect a physical measurement, but is artificially
placed to represent the Compton-thick regime. Vertical error bars reflect the
$1\sigma$ binomial standard deviation.

The _statmorph_ computer code reads a galaxy image and its associated weight
map (the $1\sigma$ error image, also known as the ‘sigma image’ in _Galfit_
and similar image analyses) to examine both the brightness and spatial
distribution of its pixel values. It then calculates various measures of
galaxy morphology, including the non-parametric measures of galaxy structure:
asymmetry (A), clumpiness (S), concentration index (C), Gini index (G), and
moment of light ($M_{20}$) (Lotz et al. 2004, Wu 1999, Bershady et al. 2000,
Conselice et al. 2000, Schade et al. 1995, Abraham et al. 1994, Abraham et al.
1996, Isserstedt & Schindler 1986, Takamiya et al. 1999, Conselice et al.
2003). For our analysis, we focused on the _shape_ asymmetry parameter,
$A_{s}$, a relatively recent addition to the suite of galaxy morphological
indicators introduced in Pawlik et al. (2016), explicitly to automate the
recognition of galaxies with faint asymmetric tidal features suggestive of an
ongoing or past merger. $A_{s}$ is calculated in the same way as the standard
asymmetry parameter, except the measurement is performed on the binary
detection mask (segmentation map) of a galaxy rather than the galaxy image
itself (see Pawlik et al., 2016 and Rodriguez-Gomez et al., 2019 for details).
Therefore, while _A_ considers asymmetries in both the spatial and brightness
distributions of the pixel values in a galaxy image, $A_{s}$ considers only
those due to a spatial asymmetry, thereby isolating disturbances that can be
attributed to a physically disruptive event such as a major merger, as opposed
to an internal galaxy sub-structure unrelated to merging activity that would
qualify as a brightness asymmetry.

The binary detection mask utilized in the _statmorph_ calculation of $A_{s}$
is defined as the contiguous group of pixels brighter than 1$\sigma$ above the
mode of the pixel values in the background-subtracted, smoothed galaxy image,
that includes the brightest pixel. Examples of the statmorph figure output
showing the details of the $A_{s}$ calculation are displayed in Figure 3.
Pawlik et al. (2016) proved unequivocally that $A_{s}$ far outperforms the
classic A parameter in detecting low surface brightness features around the
edges of a galaxy. In fact, until the addition of $A_{s}$ to the suite of the
classic CAS (concentration-asymmetry-clumpiness) and Gini coefficient model-
independent measures of galaxy morphology, visual inspection was still
considered a more accurate method to detect merger signatures than these
parameters (see, e.g., Lotz et al., 2008, Darg et al., 2010, Treister et al.,
2012).

For each of the 289 galaxies analyzed in our sample, we inspected and verified
the automatically generated statmorph results to ensure proper source
detection and aperture size and location. We identified and remedied instances
where obvious foreground objects or image features/artefacts obviously not
physically associated with the source being analyzed, were being included in
the aperture used to calculate $A_{s}$. This was done by iteratively
redefining the source aperture until such features were no longer erroneously
being included in the $A_{s}$ measurement. Therefore, the ultimate methodology
we employed in our analysis of AGN host morphology was a synergistic union of
human reasoning and computer vision.

###  3.2 Trial Design

Our analysis of AGN host galaxy morphology was conducted in two distinct
trials (hereafter Trial 1 and Trial 2): the first considering only the X-ray-
selected subset of AGN from Lyu et al. (2022) described in Section 2.1.1; and
the second, those identified as mid-IR AGN without X-ray counterparts in Lyu
et al. (2023), described in Section 2.1.2. To test the hypothesis of our
study, we organized the trials in the following way: the X-ray AGN sample in
Trial 1 was split into a test set of highly obscured sources ($N_{H}>10^{23}$
$cm^{-2}$) and two control sets of unobscured-to-lightly obscured
($N_{H}<10^{22}cm^{-2}$) and moderately ($10^{22}<N_{H}<10^{23}$ $cm^{-2}$)
obscured sources, along the lines of the analysis in Kocevski et al. (2015),
using the SED-modeled hydrogen absorbing column density parameter, $N_{H}$, as
the measure of AGN obscuration. Trial 2 was designed to represent an
overlapping extension of Trial 1 in terms of an assumed increasing level of
obscuration into the Compton-thick regime, as we require the mid-IR AGN sample
to be undetected at X-ray wavelengths.

Figure 6: Top: X-ray intrinsic luminosities for the X-ray AGN sample as a
function of redshift. Middle and bottom: X-ray and MIRI AGN sample $M_{*}$ and
$L_{bol}$ as a function of redshift. Vertical gray lines mark the bounds of
the four redshift bins in which the host galaxy morphologies of the AGN sample
were separately analyzed and compared.

While the X-ray AGN sample size is not large enough to simultaneously match
the test set of heavily obscured AGN to the moderately and lightly obscured
control sets in narrow ranges of redshift (_z_), stellar mass ($M_{*}$),
bolometric luminosity ($L_{bol}$), and intrinsic X-ray luminosity
($L_{x,int}$, X-ray $0.5-7$ keV luminosity corrected for absorption), we found
that dividing the sources into redshift bins of width unity between $z=0$ and
$z=4$ had the effect of naturally binning on most of these other parameters.
Where the primary redshift binning unavoidably led to a relatively coarse
secondary bin on some of the other galaxy parameters tested, with ranges
spanning $1-2$ orders of magnitude in a given redshift bin, we repeated our
analysis on each such parameter separately, to ensure no hidden biases in our
final results. The chosen redshift bins, displayed in Figure 6, represent the
maximum number of evenly spaced bins that allow for a statistically
significant sample of galaxies in each, and are expected to be numerous enough
to mitigate the effect of a redshift bias within a given bin. In other words,
the probability of detecting morphological features in galaxies lying at
similar redshifts should be roughly equal, given that any observed differences
in the their brightness and size would be due to intrinsic differences between
the galaxies, and not to vastly different distances from the observer (see,
e.g., Bamford et al., 2009). We also restricted our analysis to AGN in massive
host galaxies with SED-derived (3D-HST photometry) stellar masses of
$M_{*}>10^{9}M_{\odot}$ (Lyu et al., 2022), given the possibility that the
dominant physical processes shaping the evolution of low-mass active galaxies
may differ from those molding more massive ones; and that, in any case, the
sample is likely incomplete and potentially biased for lower-mass hosts (Lyu
et al., 2023).

What is evident in Figure 6 is the selection bias that plagues the X-ray
detection of AGN, that low-luminosity, highly obscured AGN fall below the
detection limits of the deepest X-ray surveys, leaving only the brightest
sources to be selected at each redshift. The effect of this observational bias
on the Trial 1 analysis is that the most obscured sources exhibit a median
intrinsic X-ray luminosity significantly higher than the least obscured
control set between z=0 and 1, and therefore these two subsets may not be
meaningfully compared in this redshift bin. However, the moderately obscured
control sample is well-matched in luminosity to both the highly obscured test
set and least obscured control set in each redshift bin, and therefore
‘bridges the gap’, so to speak, in the Trial 1 analysis.

For Trial 2, we split the mid-IR-bright/X-ray-faint AGN into the same redshift
bins as for Trial 1, and repeated the same procedural analysis.

Figure 7: JADES v1 NIRCam F150W mosaic image cutouts ($6"\times 6"$) showing
rest-optical AGN host galaxy emission of a representative sampling of the
X-ray- and mid-IR-selected sources lying at $z=1-2$. In these images, used in
both the human and computer visual analyses of the present study, it can be
seen that the galaxy morphologies range from visually undisturbed, symmetric
shapes to highly asymmetric/disturbed. For the visual analysis task, the
images were examined using a variety of SAOImage DS9 scales, stretches, and
colors, to ensure the most accurate morphological classification, as
demonstrated.

##  4 Results

The results of the visual classification procedure for all trials are
displayed in Figure 4, showing the percentage of AGN host galaxies with
undisturbed (‘point source’, ‘spheroid’ , and ‘disk’ ) and disturbed
(‘merger’, ‘irregular’, and ‘disturbed disk’) morphologies. It can be seen
that the morphology predominantly characterizing the disturbed X-ray AGN hosts
is the ‘disturbed disk’ morphology, followed by irregular; the undisturbed
fraction of the X-ray AGN hosts was likewise dominated by the (symmetric) disk
morphology. The MIRI-detected sample, on the other hand, appears
overwhelmingly irregular, with a higher fraction of obvious mergers.

In Figure 5, we compare the measured fractions of disturbed AGN host galaxies
recovered by both the human and computer vision analyses. Unsurprisingly, the
results of both largely corroborate one another, and it is clear that both
methods reveal a correlation between disturbed host galaxy morphology and the
level of AGN obscuration. The main difference found was that, in the visual
classification, often what was qualitatively labeled as undisturbed, or a
‘mild’ galaxy asymmetry/disturbance to the human eye (a subcategory ‘learned’
after observing instances of significantly stronger disturbances), was
classified by the statmorph $A_{s}$ parameter as a ‘mild’ and ‘strong’
disturbance, respectively. Such undisturbed and mildly disturbed galaxies
identified in the visual classification almost exclusively showed an $A_{s}$
value between $0.1-0.2$, and $0.2-0.3$, respectively, ranges which define mild
and strong asymmetries in statmorph. As a result, the computer vision analysis
yields a higher fraction of disturbed AGN host galaxies than the visual
classification procedure. We attribute this ultimately inconsequential
difference between the two analysis methods to the numerical precision and
objectivity of the shape asymmetry algorithm over human vision in classifying
objects in an inherently unbiased way, and to its enhanced sensitivity to
morphological asymmetries such as faint tidal features that might go unnoticed
by the human eye. It proved particularly useful in the classification of the
highly obscured systems lying at the highest redshifts, as these were
sufficiently resolved to detect an asymmetry, but not enough for a fine visual
distinction between a ‘mild’ and ‘strong’ (hidden) asymmetry (which may have
led to a bias in the visual classification results where strong disturbances
were more easily detected at low redshifts). We therefore adopt the statmorph
analysis results as ‘final’, and all subsequent discussion in this paper is
based on it.

Figure 8: Histograms depicting the statmorph shape asymmetry parameter values
for the full X-ray-plus-mid-IR AGN sample. From the top to the bottom of the
column of figures, the values for the following subsets are separately
overlaid in order to highlight the increasing fraction of disturbed hosts with
$N_{H}$ level: the X-ray unobscured/lightly obscured ($N_{H}<10^{22}$
$cm^{-2}$), X-ray moderately obscured ($10^{22}<N_{H}<10^{23}$ $cm^{-2}$),
X-ray highly obscured ($N_{H}>10^{23}$ $cm^{-2}$), and Compton-thick mid-IR
AGN, respectively.

We find that both the Trial 1 X-ray and Trial 2 mid-IR AGN samples
significantly reside in host galaxies with disturbed morphologies at all
redshifts sampled, and increasingly so towards higher redshift (in agreement
with the general trend for galaxies to show a higher incidence of disturbances
with increasing redshift (e.g., Conselice et al., 2009). The incidence of
disturbance also increases with bolometric luminosity. The subset exhibiting
the least amount of disturbance overall is the least obscured and least
luminous X-ray AGN ($64\pm 6.5\%$ strongly disturbed); while the subset with
the highest fraction of strongly disturbed hosts is represented by the
presumed most heavily obscured mid-IR-bright/X-ray faint AGN, whose hosts
appear $95\pm 4\%$ disturbed. The most obscured and luminous X-ray AGN appear
just as disturbed at $94\pm 6\%$, while the moderately obscured X-ray AGN show
measurably less disturbance at $78\pm 6\%$.

Figure 9: _statmorph_ shape asymmetry ($A_{s}$) parameter values for the full
AGN sample displayed as a function of redshift. $A_{s}$ values between 0.1-0.2
represent mild asymmetry/disturbance, while values greater than 0.2 represent
strong asymmetry/disturbance.

##  5 Discussion

In both the visual classification and computer vision analyses of the JADES v1
rest-optical NIRCam images of X-ray- and mid-IR-identified AGN in GOODS-S/HUDF
(Lyu et al., 2022, 2023), we find a strong correlation between the SED-
inferred (or assumed, in the case of the mid-IR AGN sample) value of $N_{H}$,
and the statistical fraction of host galaxies with a measurable galaxy-scale
asymmetry/disturbance. This result is directly in line with the observed trend
for $z<2$ X-ray-luminous AGN analyzed in the same way in Kocevski et al.
(2015), and directly extends their analysis out to $z~{}\sim~{}4$. Like the
authors of that study, we compare a test set of heavily obscured X-ray AGN to
moderately obscured and unobscured-to-lightly obscured control counterparts,
allowing us to isolate the effect of obscuration on AGN host galaxy
morphology. Furthermore, by complementing the X-ray AGN sample with a newly
identified sample of MIRI-detected, obscured mid-IR-bright/X-ray-faint AGN, we
were able to carry our analysis into the presumed Compton-thick regime and
begin to test the evolutionary scenario for AGN put forth in Kocevski et al.
(2015) to explain their findings. These authors suggest that the most obscured
X-ray- and mid-IR-detected Seyfert systems occur relatively early in the
lifetime of a major merger, and chronologically precede the less-obscured and
fainter X-ray sources. One way to frame the outcomes of the distinct trials in
which we separately examined the host galaxy morphologies of X-ray- and mid-
IR-identified AGN, would indeed be to propose such a scenario. This appears
supported by the merger simulations of Lotz et al. (2010), that disturbed host
galaxy morphologies should be most noticeable both at the very start of a
major merger, and also in the final stages. They also show that late-stage
major merger remnants ($\sim 1$ Gyr after coalescence) tend to appear as
early-type spiral galaxies, which would suggest that our disturbed X-ray AGN
sample, that is dominated by the disk morphology, may be major-merger remnants
being ‘caught’ at a later stage than the largely irregular/train-wreck mid-IR
AGN (Bignone et al., 2017).

The strength of our analysis lies in the experimental design of our trials,
and the superior NIRCam rest-optical imagery utilized. These points are
expanded upon, below.

1) By partitioning a single and relatively complete sample of X-ray AGN that
were detected using the same method and extracted from the same data set – and
matched on numerous physical traits but exhibiting distinct levels of $N_{H}$
obscuration – we were able to eliminate a maximum number of uncontrolled
variables and observational biases that might otherwise skew the results of
our analysis. Most importantly, this approach allowed us to test a model-
predicted evolutionary scenario in which the obscured and less-obscured phases
of AGN occur in chronological succession of one another. Had we compared a
test set of AGN hosts to a control sample of non-active counterparts without
distinction between the obscured and unobscured subsamples (see, e.g.,
Kocevski et al., 2012), we could easily have missed this connection. In other
words, the experimental design of our trials was aimed at directly testing the
narrowly focused scientific question of whether or not the obscured AGN phase
is strongly associated with mergers, as opposed to measuring the fraction of
active galaxies versus non-active galaxies associated with mergers (which is a
worthwhile and necessary pursuit that fits into our broader scientific study,
but is beyond the scope of the current paper).

2) The capabilities of JWST have allowed the measurement of detailed
morphologies of AGN host galaxies in rest-optical imagery with high signal-to-
noise and out to high redshift. In addition, the follow-up to our visual
classification of AGN host morphologies with a computational analysis using
the relatively new statmorph shape asymmetry parameter, was key in allowing us
to circumvent the obstacles faced by previous authors: we observed firsthand
how a non-negligible number of mildly asymmetric/disturbed disk galaxies in
our sample could go unnoticed by the human eye and initially be classified as
‘undisturbed’ – only to discover in the results of the computer vision
analysis that they are, in fact, significantly asymmetric. This can be
attributed to the inherently unbiased nature of the _statmorph_ $A_{s}$
algorithm in its assignment of equal weights to all pixels in an input image
belonging to a galaxy, using the binary segmentation image of the galaxy, i.e.
without consideration of the brightness variations amongst the pixels, as
traditionally done with classic non-parametric morphology indicators. This
translates to detecting asymmetries resulting purely from
spatial/morphological disturbances, as opposed to an asymmetry in pixel
brightness values, e.g. due to a galaxy sub-structure not caused by or related
to physically disruptive merging activity. At the time of writing, to the
extent of the authors’ knowledge, the present study represents the first to
quantitatively tie ongoing or past merging activity to the obscured phases of
AGN through the measure of AGN host galaxy spatial morphology.

Even if we were to take a conservative approach and entirely remove the
statmorph computer vision analysis to use only visual classification, in
keeping with all previous studies on the topic, our main result would still
hold. In this case, we would effectively re-classify a significant number of
the _statmorph_ -detected-but-invisible “mildly disturbed” cases ‘back’ to the
visually undisturbed category, still leaving the “strongly disturbed” cases
that were observed using both the visual and computational methods, and that
represent the bulk of the disturbed fraction.

###  5.1 Comparison to previous observational studies

We confirm and expand upon the main conclusion of Kocevski et al. (2015), that
a majority of the heavily obscured X-ray AGN at sub-quasar luminosities are
found in hosts with a measurable galaxy-scale spatial asymmetry/disturbance.
We additionally observe that the overall fraction of AGN hosts exhibiting a
disturbed morphology in rest-optical imagery grows in step with redshift in
both the X-ray and mid-IR samples of our analysis, along with $N_{H}$ and
intrinsic X-ray luminosity in the X-ray sample, and bolometric luminosity in
the mid-IR-bright/X-ray-faint subset. The hosts of the most obscured X-ray AGN
($N_{H}>10^{23}cm^{-2}$), and the Compton-thick mid-IR-bright/X-ray-faint MIRI
AGN, represent the most disturbed subsets of the entire sample showing $\sim
94\%$ and $\sim 95\%$ strong disturbance, respectively, and rank consistently
among the most disturbed in all data bins considered in our trials with fixed
redshift, stellar mass, and intrinsic X-ray luminosity, between 85% and 100%
percent disturbed. The moderately obscured ($10^{22}<N_{H}<10^{23}$ $cm^{-2}$)
X-ray AGN show a slightly smaller strongly disturbed host fraction of 78%,
while the least/unobscured ($N_{H}<10^{22}$ $cm^{-2}$) show the lowest of the
entire sample at $\sim 64\%$.

In Kocevski et al. (2015), a picture of AGN evolution is suggested that would
unify the results of their analysis through a chronological ordering of the
obscured and unobscured X-ray AGN phases considered in their study. It
proposes that the X-ray AGN exhibiting the highest (sub-quasar) $L_{x,int}$
and $N_{H}$ values, and level of galaxy-scale morphological disturbance,
appear relatively early in the lifetime of a major merger, before a fully
enshrouding dusty torus begins to block all X-ray (and optical) emission from
the accretion disk. This heavily obscured phase, during which the AGN may only
be detected via significant warm and hot dust re-emission in the mid-IR, and
marking the presumed period of significant black hole growth, would then
transition to the briefly visible quasar phase as the increase in radiative
power from the nucleus reactively “blows out” the obscuring dust and gas. This
would ultimately lead to the presentation of a naked, unobscured X-ray AGN
with a steadily diminishing X-ray luminosity and increasingly more ordered and
relaxed morphology. This observed strong correlation of host galaxy
morphological disturbance with $N_{H}$ level in our own analysis of obscured
and unobscured/lightly obscured X-ray AGN, taken together with the trends in
the luminosities of each of the X-ray and mid-IR AGN subsamples discussed in
Section 4, appears to fit into the same evolutionary picture.

###  5.2 Comparison to models and merger simulations

The AGN-fueling model put forth in Hopkins et al. (2014) predicts that quasars
are predominantly merger-fueled, as supported by numerous empirical studies
(e.g. Hopkins & Hernquist, 2006, Glikman et al., 2015, Assef et al., 2016, Fan
et al., 2016, Treister et al., 2012), while stochastic processes dominate SMBH
growth in low-to-moderate luminosity AGN with $M_{BH}<10^{7}$ $M_{\odot}$ (see
Treister et al., 2012). However, given the expectation that the bulk of the
faint/obscured AGN population remains largely hidden below the detection
limits of even the deepest X-ray surveys, there is the possibility that this
model has been significantly under-tested in the mass and luminosity regime of
Seyfert galaxies. The newly uncovered set of obscured MIRI AGN added to our
comprehensive X-ray sample from Lyu et al. (2022), places us at a unique
vantage point to test this very theory. Furthermore, it is also possible, as
Hopkins et al. (2014) directly address, that any significant evidence of
mergers that may exist in the low- and intermediate-luminosity AGN population
has simply gone undetected due to the inherent faintness and short lifetime of
merger signatures such as tidal tails. While the results of our study agree
with theirs superficially in the prediction of a rise in merger-driven SMBH
accretion with redshift and luminosity, reflecting the expectation of a
cosmological merger rate rapidly increasing with redshift, we observe this
trend in the realm of AGN luminosities and SED-inferred black hole masses far
lower than predicted by Hopkins et al. (2014). Therefore, we must consider
that, while this model is physically sound and has been proven to accurately
predict the dominant fueling mechanism of the brightest and most massive (i.e.
most easily detected) AGN (see, e.g. Treister et al., 2012), that its
predictions in the Seyfert regime need to be improved.

####  5.2.1 On the predominance of disk morphologies

We also find, like Kocevski et al. (2015), that the majority of the X-ray AGN
hosts in our study are disks, as opposed to the point sources or spheroids
that might be expected to emerge from a merger. In fact, their finding of a
high disk fraction amongst obscured AGN hosts to a luminosity $10\times$
brighter than predicted by Hopkins & Hernquist (2006), prompted those authors
to revise the merger-induced fueling component of their model in Hopkins et
al. (2014) to better reflect such observations. One of the major claims of
Hopkins et al. (2014), based on both their model predictions and observations
they cite in support of them, is that an AGN host with a disk morphology
should be taken as evidence of being ‘stochastically fueled’, that one would
expect a ‘merger-fueled’ AGN to appear as bulge-dominated in a post-merger
phase. The implication of this assumption on their models for AGN fueling is
that at sub-quasar luminosities and at $M_{BH}<10^{7}M_{\odot}$, SMBH growth
should be dominated by secular, passive processes, at all redshifts, and by
merger-induced growth at quasar luminosities and $M_{BH}>10^{7}M_{\odot}$. Our
findings are in direct disagreement with this model, which would be alarming
had Hopkins et al. (2014) not explicitly highlighted the tremendous difficulty
of detecting inherently faint tidal signatures around a galaxy of a recent, or
even ongoing merger, even in local systems, compounded by surface brightness
dimming for more distant galaxies (see, e.g., Lotz et al., 2008). As a result
of this, there is a record in the literature of major mis-classifications of
such galaxies as “relaxed” using HST imagery, which was primarily and
necessarily relied upon before the advent of JWST and its superior sensitivity
and resolution (see, e.g. Hopkins et al., 2009b, Lotz et al., 2004).

It is also important to consider that Hopkins et al. (2009) thoroughly
demonstrate in simulations of mergers of varying mass ratios, gas content, and
orbital parameters, that it is quite common for disks to survive or reform
following a merger, based on simple laws of gravitational physics and
independent of the details of ISM gas physics or feedback processes. Darg et
al. (2010) take their analysis further and provide empirical evidence that the
timescale of detectability for mergers involving spiral galaxies is longer
than those involving ellipticals, as theoretically predicted by Lotz et al.
(2008). Furthermore, the merger simulation study of Bignone et al. (2017) that
uses Illustris TNG mock galaxy data to constrain the non-parametric
morphologies of mergers, finds that rotationally supported gas-rich structures
can form from residual angular momentum after the final coalescence of the
merging galaxies, despite the Illustris galaxy formation model including
feedback from quasar- and radio-mode AGN, and stellar supernovae; this finding
is also supported by Robertson et al. (2006) and Snyder et al. (2015).
Finally, the merger simulations of Lotz et al. (2008) show that the remnants
of equal-mass, gas-rich isolated mergers appear disk-like when observed at
$t>1$ Gyr after the final coalescence. Therefore, the fact that the members of
the X-ray AGN sample studied here lie mostly in disk morphologies, as opposed
to the spheroidal, bulge-dominated systems historically expected to emerge
from a merger, should not be as surprising as Hopkins et al. (2014) would
indicate; it cannot be used as direct evidence against their formation in a
past merger (Hopkins & Hernquist, 2006). It could still be the case, however,
that AGN residing in fully relaxed or very mildly disturbed disks represent a
phase sufficiently long after a potential past merger that their primary
fueling mechanism has transitioned to more passive processes.

####  5.2.2 Redshift-dependent disturbance

Several features of our analysis appear to support the model-predicted
increase in merger-driven accretion with redshift, if we can safely assume
that merging activity is correlated with obscuration and AGN bolometric
luminosity. In other words, we observe an increase in the level of AGN host
galaxy disturbance and level of $N_{H}$ obscuration – which we assume to be
merger-driven – with redshift, which also traces an increase in $L_{bol}$ in
our sample. We do not observe many unobscured/lightly obscured X-ray AGN at
$z>2$, most in our sample reside at $z<1$; conversely, very few of the highly
obscured and most disturbed AGN lie at the lowest redshifts. This trend
reflects the X-ray survey detection limit where we miss the faintest
high-$N_{H}$ AGN at high redshift, but it is also not clear if we can even
expect a significant number of obscured systems below the detection limit at
these high redshifts (see, e.g. Hickox & Alexander, 2018). In other words, we
may be witnessing a true redshift-dependence of the AGN level of obscuration
and luminosity tied to an increasing level of host galaxy disturbance, and not
merely an artificial result from selection bias.333However, the explanation
for why we do not observe luminous unobscured AGN in our sample is less clear,
but may very well reflect the rarity of catching a snapshot of the relatively
short-lived optical quasar phase ($\sim 10^{7}$ years), in comparison to the
longer pre- and post-quasar phases (see Hopkins & Hernquist, 2006 and
references therein). There is also the expectation of a decreasing number of
massive gas-rich mergers with decreasing redshift and associated predominance
of dry mergers at $z<\sim 1.5$ (Treister et al., 2012), reflecting an
antihierarchical black hole mass growth in which quasar peak luminosities and
(final black hole masses) markedly decline at $z<2$ (Hopkins & Hernquist,
2006).

Finally, of note is that in the redshift bin of our analysis sampling AGN at
$z=1-2$, all subsets appear with a high number of disturbed hosts, including
the low-$N_{H}$, unobscured-to-lightly-obscured X-ray AGN. While we observe
that these systems are found in nearly equal numbers of disturbed and
undisturbed hosts at $z<1$, in agreement with numerous other studies and the
expectations of the AGN fueling model of Hopkins & Hernquist (2006), the
fraction of disturbed hosts among this subset jumps from 64% at $z<1$ to 91%
between $z=1-2$, in the same realm as their more obscured counterparts at
these redshifts – though we caution that this data bin may be skewed by small-
number statistics. Even so, we observe the fraction of mildly disturbed
unobscured hosts increases from 40% at $z<1$ to 72% at $z>1$, which may be
significant. Likewise, the disturbed host fraction of the moderately obscured
counterparts sharply rises from 65% to a maximum of 95% between $z=1-2$, with
a constant ratio of mild-to-strong disturbance. This observed uptick in the
disturbed host fraction in our sample at z$\sim$1 matches the findings of
Treister et al. (2012), who concluded that the dominant source of SMBH growth
appears to transition from secular processes to mergers at this redshift,
based on an analysis of ten archival surveys of moderate-luminosity AGN at
$z<3$.

##  6 Outstanding Questions

###  6.1 Major, minor, or accumulated mini mergers?

While it is accepted in the literature that galaxy disturbances measured using
conventional non-parametric morphology indicators are to be interpreted as
merger signatures, the asymmetry measurements alone cannot generally be used
to determine whether the disturbance is caused by a major or minor merger, or
a steady accumulation of mini mergers. The distinction is important in the
context of our study, given the assumption that only gas-rich major mergers
can significantly contribute to the growth of SMBHs. Lotz et al. (2008) and
Bignone et al. (2017) correlate key non-parametric morphology measures with
the recent merger history of galaxies to show that the $Gini-M_{20}$ criterion
is effective at identifying both major and minor mergers over a significant
lifetime of the merger. However, it is most complete and effective in the pre-
coalescence phase while a double nucleus is visible, given the nature of the
measurement. These authors also show that the classic asymmetry (A) parameter
can be used to loosely distinguish major from minor mergers where two nuclei
are still visible. Even so, Bignone et al. (2017) emphasize that these metrics
can be quite incomplete in their identification of mergers, especially in the
post-coalescence phase. In addition, Pawlik et al. (2016) showed that $A_{s}$
is much more sensitive to faint tidal signatures indicative of a merger than
A, which was the very motivation for their invention of $A_{s}$. Furthermore,
the simulations in Bignone et al. (2017) showed that A classifies as disturbed
all galaxies with stellar masses less than $10^{10.5}$ $M_{sun}$, independent
of a merger event. Therefore, we also consider this parameter to be unreliable
in identifying merger remnants. As such, we are not able to show with the
statmorph non-parametric measures alone if the AGN hosts in our sample are
likely to be the remnants of major or minor mergers, given that only a small
handful of our entire sample meet the double-nucleus criterion. However, the
presentation of a significant number of strongly disturbed disk morphologies
in the X-ray sample suggests they could likely be the remnants of the gas-rich
mergers that are believed to fuel SMBH in the Universe, as 1) the presence of
a disk indicates the presence of gas; 2) merger simulations (Lotz et al.,
2008; Bignone et al., 2017) consistently show that spiral morphologies can
reform after major gas-rich mergers and appear as disturbed disks to $\sim 1$
Gyr post-coalescence; and 3) at any given snapshot along the cosmic timeline
during which merging galaxies are expected to occur in the Universe, we are
more likely to observe the remnants of major mergers than minor mergers, given
the significantly longer timescale of observability ($\sim 2$ Gyr) of the
former (Lotz et al., 2010; Bignone et al., 2017). However, the simulations of
Lotz et al. (2008, 2010) did not include the effects of AGN feedback, and
therefore it is not clear if a gaseous disk could survive or reform post-
major-merger in more realistic conditions. They also did not consider mergers
involving mixed-morphologies, e.g. a disk merging with a spheroid.

The Illustris merger simulations of Bottrell et al. (2023) find that so-called
mini mergers ($0.01\leq\mu<0.1$) are expected to occur significantly more
frequently than minor and major mergers, and therefore may represent a
significant, overlooked cause of galaxy morphological disturbance. While the
simulations show that individual mini mergers are not likely to induce a
significant observable disturbance, they reveal that the integrated (classic
A, light-weighted) asymmetry due to mini mergers over $0.01\leq z\leq 0.7$
constitutes $70\%$ of the asymmetric structure in the mergers of their study.
They therefore suggest that the aggregation of asymmetry enhancements by a
series of mini merger events for a given galaxy may be driving the build-up of
light asymmetry and lopsidedness in galaxies, as opposed to major or minor
mergers. However, their study examined only TNG50 star-forming galaxies and
adopted the classic asymmetry parameter as the measure of galaxy disturbance,
while we consider active galaxies and utilize the shape asymmetry parameter
that ignores asymmetries in the galaxy surface brightness distribution. As a
result, it is not immediately clear if the results of their study have direct
implications for ours.

##  7 Conclusion

This study directly builds upon and reaffirms the results of Kocevski et al.
(2015), who compare the morphologies of a significant number of obscured
X-ray-selected AGN to their unobscured countertparts, in an attempt to test
the suspected connection between obscuration and merger activity. Like
Kocevski et al. (2015), we split an X-ray sample of AGN into ‘test’ highly
obscured and ‘control’ moderately and lightly obscured subsets, and match them
on stellar mass, bolometric luminosity, and redshift. However, we were able to
extend the analysis into the Compton-thick regime using a newly obtained
sample of mid-IR-bright/X-ray-faint JWST/MIRI AGN; and also into the high-
redshift regime to z$\sim$4, using highly resolved, high-signal-to-noise
JWST/NIRCam rest-optical imagery from the JADES v1 data release. We find a
significant correlation between host galaxy disturbance and level of AGN
obscuration amongst a complete, multi-wavelength sample of obscured Seyfert
AGN at $z=0-4$, using both a qualitative visual classification and
quantitative measure of the shape asymmetry of each galaxy in question. From
the observed pattern of disturbances, we conclude that mergers are common
amongst the obscured AGN population, and therefore that the obscured AGN phase
may mark a period of significant supermassive black hole (SMBH) growth – in
direct tension with the leading model on AGN fueling mechanisms (Hopkins et
al., 2014). Specifically, there is the possibility that the merger-fueled AGN
phase predicted by Hopkins et al. (2014) to exclusively reside in the quasar
regime, extends to the significantly lower Seyfert luminosities and black hole
masses sampled here. This may be the result of their model missing a
significant fraction of faint obscured AGN that are now ‘coming to light’, and
that may represent relics of quasar activity in the distant past (so-called
‘dead quasars’). If such systems were in fact born from mergers and not other,
more passive processes – as might be indicated by the observation that the
Seyfert luminosity function smoothly joins onto the quasar luminosity function
(Hao et al., 2005) – they could be expected to be sporadically revived by such
secular processes, leading to their presentation as ‘non-mergers’.

This theory appears to be supported by McAlpine et al. (2020), a study that
utilizes the cosmological hydrodynamical EAGLE simulation to investigate the
degree to which black hole (BH) activity is enhanced during the period of a
major merger. This study finds that the galaxies embodying the most optimal
conditions for triggering an AGN via a merger have the relatively low stellar
masses, black hole masses, and luminosities typical of the AGN in our sample.
The second key discovery of their study is that 50-75% of major-merger-
triggered BH activity occurs in post-merger systems, after the two galaxies
have coalesced. This implies that the majority of observational merger
studies, which focus on the earliest stages of a merger, may have missed a
significant fraction of the peak of BH growth in their samples, and mistakenly
classified the post-merger remnants as isolated, non-merger AGN.

The strength of our conclusion not only lies in the corroboration of the
results of Kocevski et al. (2015) and expansion to higher redshifts and levels
of AGN obscuration, but also in the agreement with the various merger
simulations discussed. In Paper II, we will expand upon the analysis of this
paper with supplementary data from the literature to further investigate if
the Seyfert population sampled may evolve according to a lower-luminosity,
lower-mass extension of the accepted quasar evolution models put forth in
Hopkins & Hernquist (2006), Hopkins et al. (2009), and Hopkins et al. (2014)
as the data would suggest.

Here we show the JADES v1 mosaic NIRCam cutout images of a randomly chosen and
representative selection of the AGN in our sample at the effective wavelengths
available (listed in the image captions), to demonstrate the progressive
resolution degradation inevitably caused by an increasing PSF size towards
longer wavelengths. This instrumental effect motivated our decision to include
only the F150W mosaic imagery in our analysis, as it represents the mosaic
at/near peak NIRCam resolution and sensitivity that also maximizes the
coverage of our AGN sample in the GOODS-S field. We observed, after repeatedly
visually examining a statistically significant number images of the AGN host
galaxies in a number of trials, that no information would be gained or lost,
or conclusions altered, by measuring the AGN host galaxy morphology at a
constant rest-frame wavelength, requiring our analysis procedure to be
repeated using mosaic image cutouts at a range of observed wavelengths other
than (longward of) F150W. It can be witnessed here in all images shown, that
the morphological identification of each galaxy as measured in the F150W image
as symmetric/undisturbed or asymmetric/disturbed, remains unchanged in the
images exposed at other wavelengths; all galaxy images exposed at wavelengths
longward of F150W appear identical to the galaxy image at F150W, except at
poorer resolution. We also note a similar observance of the images taken at
wavelengths shortward of F150W: while they appear equally or slightly more
resolved than the F150W image (but noisier), the galaxy morphology appears
identical and therefore no new information would be gleaned from inclusion of
these images in our analysis.

Figure 10: JADES v1 mosaic image cutouts of lightly obscured
($N_{H}<10^{22}cm^{-2}$) X-ray-selected AGN at $0<z<1$ (Lyu et al., 2022) in
the following NIRCam bands (left to right): F090W, F115W, F150W, F200W, F277W,
F356W, F444W.

Figure 11: JADES v1 mosaic image cutouts of lightly obscured
($N_{H}<10^{22}cm^{-2}$) X-ray-selected AGN at $1<z<2$ (Lyu et al., 2022) in
the following NIRCam bands: (Top) F090W, F115W, F150W, F200W, F277W, F356W,
F444W; (Middle, Bottom) F070W, F090W, F115W, F150W, F200W, F277W, F356W, F444W

Figure 12: JADES v1 mosaic image cutouts of moderately obscured
($10^{22}<N_{H}<10^{23}$ $cm^{-2}$) X-ray-selected AGN at $0<z<1$ (Lyu et al.,
2022) in the following NIRCam bands: (Top, Middle) F090W, F115W, F150W, F200W,
F277W, F356W, F444W; (Bottom) F070W, F090W, F115W, F150W, F200W, F277W, F356W,
F444W

Figure 13: JADES v1 mosaic image cutouts of moderately obscured
($10^{22}<N_{H}<10^{23}$ $cm^{-2}$) X-ray-selected AGN at $1<z<2$ (Lyu et al.,
2022) in the following NIRCam bands: F090W, F115W, F150W, F200W, F277W, F356W,
F444W

Figure 14: JADES v1 mosaic image cutouts of moderately obscured
($10^{22}<N_{H}<10^{23}$ $cm^{-2}$) X-ray-selected AGN at $2<z<3$ (Lyu et al.,
2022) in the following NIRCam bands: F090W, F115W, F150W, F200W, F277W, F356W,
F444W

Figure 15: JADES v1 mosaic image cutouts of highly obscured
($N_{H}>10^{23}cm^{-2}$) X-ray-selected AGN at $0<z<1$ (Lyu et al., 2022) in
the following NIRCam bands: (Top, Middle) F090W, F115W, F150W, F200W, F277W,
F356W, F444W; (Bottom) F070W, F090W, F115W, F150W, F200W, F277W, F356W, F444W

Figure 16: JADES v1 mosaic image cutouts of highly obscured
($N_{H}>10^{23}cm^{-2}$) X-ray-selected AGN at $1<z<2$ (Lyu et al., 2022) in
the following NIRCam bands: F090W, F115W, F150W, F200W, F277W, F356W, F444W

Figure 17: JADES v1 mosaic image cutouts of highly obscured
($N_{H}>10^{23}cm^{-2}$) X-ray-selected AGN at $2<z<3$ (Lyu et al., 2022) in
the following NIRCam bands: F115W, F150W, F200W, F277W, F356W, F444W

Figure 18: JADES v1 mosaic image cutouts of highly obscured
($N_{H}>10^{23}cm^{-2}$) X-ray-selected AGN at $3<z<4$ (Lyu et al., 2022) in
the following NIRCam bands: F115W, F150W, F200W, F277W, F356W, F444W

This work was supported by NASA grants NNX13AD82G, 1255094, and the NASA
contract for NIRCam to the University of Arizona, NAS5-02015. This work is
based on observations made with the NASA/ESA/CSA James Webb Space Telescope.
SA, NB, DJE, KH, ZJ, BR, MR, and CNAW acknowledge support from the NIRCam
Science Team contract to the University of Arizona, NAS5-02015. DJE is further
supported as a Simons Investigator. SA, JL, IS, and GHR acknowledge support
from the JWST Mid-Infrared Instrument (MIRI) Science Team Lead, grant
80NSSC18K0555, from NASA Goddard Space Flight Center to the University of
Arizona. These observations are associated with JWST GTO programs #1180 and
#1207, and GO program #1963. AJB acknowledges funding from the “FirstGalaxies”
Advanced Grant from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (Grant agreement No.
789056). The work of CCW is supported by NOIRLab, which is managed by the
Association of Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation. WB and FDE
acknowledge support by the Science and Technology Facilities Council (STFC)
and by the ERC through Advanced Grant 695671 ”QUENCH”. JWST

## References

  * Abraham et al. (1996) Abraham, R. G., Tanvir, N. R., Santiago, B. X., Ellis, R. S., Glazebrook, K., & van den Bergh, S. 1996, MNRAS, 279, L47. doi=10.1093/mnras/279.3.L47
  * Abraham et al. (1994) Abraham, R. G., Valdes, F., Yee, H.K.C., & van den Bergh, S. 1994, ApJ, 432, 75. doi=10.1086/174550
  * Alberts et al. (2020) Alberts, S., Rujopakarn, W., Rieke, G. H., et al. 2020, ApJ, 901, 168. doi=10.3847/1538-4357/abb1a0
  * Assef et al. (2016) Assef, R. J., Walton, D. J., Brightman, M., et al. 2016, ApJ, 819, 111. doi=10.3847/0004-637X/819/2/111
  * Bamford et al. (2009) Bamford, S., et al. 2009, MNRAS, 393, 4. doi:10.1111/j.1365-2966.2008.14252.x
  * Bershady et al. (2000) Bershady, M., Jangren, A., & Conselice, C. 2000, AJ, 119, 2645. doi=10.1086/301386
  * Bignone et al. (2017) Bignone, L. A., et al. 2017, MNRAS, 465, 1. doi=10.1093/mnras/stw2788
  * Blecha et al. (2018) Blecha, L., Snyder, G. F., Satyapal, S., et al. 2018, MNRAS, 478, 3056. doi=10.1093/mnras/sty1274
  * Bottrell et al. (2023) Bottrell, C., et al. 2023, MNRAS, doi=10.1093/mnras/stad2971
  * Brinchmann et al. (1998) Brinchmann, J., et al. 1998, 499, 1. doi=10.1086/305621
  * Bundy et al. (2005) Bundy, K., Ellis, R., & Conselice, C. 2005, ApJ, 625, 2. doi=10.1086/429549
  * Darg et al. (2010) Darg, S. et al. 2010, MNRAS, 401, 2
  * Conselice et al. (2000) Conselice, C., Bershady, M., & Jangren, A. 2000, ApJ, 529, 886. doi=10.1086/308300
  * Conselice et al. (2003) Conselice, C. 2003, ApJS, 147, 1. doi=10.1086/375001
  * Conselice et al. (2009) Christopher J. Conselice, C.J., Yang, C., & Bluck., A.F.L. 2009, MNRAS, 394, 4, doi=10.1111/j.1365-2966.2009.14396.x
  * Del Moro et al. (2016) Del Moro, A., Alexander, D. M., Bauer, F. E., et al. 2016, MNRAS, 456, 2105. doi=10.1093/mnras/stv2748
  * Duras et al. (2020) Duras, F., et al. 2020, å, 636, A73. doi=10.1051/0004-6361/201936817
  * Fan et al. (2016) Fan, L., et al. 2016, ApJ, 822:L32. doi:10.3847/2041-8205/822/2/L32
  * Ferrarese & Merritt (2000) Ferrarese, L. & Merritt, D. 2000, ApJ, 539, L9. doi:10.1086/312838
  * Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13. doi:10.1086/312840
  * Gilli et al. (2007) Gilli, R., Comastri, A., & Hasinger, G. 2007, A&A, 463, 79. doi:10.1051/0004-6361:20066334
  * Glikman et al. (2015) Glikman, E., et al. 2015, ApJ, 806, 218. doi=10.1088/0004-637X/806/2/218
  * Hao et al. (2005) Hao, L., et al. 2005, AJ, 129, 1795, doi=10.1086/428486
  * Henry et al. (2008) Henry, A. L., et al. 2008, ApJ, 680, L97. doi=10.1086/589944
  * Hernández-Toledo et al. (2023) Hernández-Toledo, H. M., Cortes-Suárez, E., Vázquez-Mata, J. A., et al. 2023, MNRAS, 523, 4164. doi:10.1093/mnras/stad1425
  * Hickox & Alexander (2018) Hickox, R. C. & Alexander, D. M. 2018, ARA&A, 56, 625. doi:10.1146/annurev-astro-081817-051803
  * Hopkins & Hernquist (2006) Hopkins, P. F., Somerville, R. S., Hernquist, L., et al. 2006, ApJ, 652, 864. doi:10.1086/508503
  * Hopkins et al. (2008) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2008, ApJS, 175, 356. doi:10.1086/524362
  * Hopkins et al. (2009) Hopkins, P. F., Cox, T. J., Younger, J. D., Hernquist, L., et al. 2009, ApJ, 691, 2, doi = 10.1088/0004-637X/691/2/1168
  * Hopkins et al. (2009b) Hopkins, P. F., & Hernquist, L., 2009, ApJ, 694, 1, doi=10.1088/0004-637X/694/1/599
  * Hopkins et al. (2014) Hopkins, P.F., Kocevski, D. D., Bundy, K., et al. 2014, MNRAS, 445, 1. doi:10.1093/mnras/stu1738
  * Isserstedt & Schindler (1986) Isserstedt, J., & Schindler, R. 1986, å, 167, 11.
  * Zhiyuan et. al. (2022) Ji, Z., et al. 2022, ApJ, 925, 74. doi: 10.3847/1538-4357/ac3837
  * Jogee (2006) Jogee, S. 2006, Physics of Active Galactic Nuclei at all Scales, 143. doi:10.1007/3-540-34621-X_6
  * Kampczyk et al. (2007) Kampczyk, P., et al. 2007, ApJS, 172 329. doi=10.1086/516594
  * Kartaltepe et al. (2015) Kartaltepe, J. S., Mozena, M., Kocevski, D., et al. 2015, ApJS, 221, 11. doi:10.1088/0067-0049/221/1/11
  * Kim et al. (2017) Kim, M., Ho, L. C., Peng, C. Y., et al. 2017, ApJS, 232, 21. doi:10.3847/1538-4365/aa8a75
  * Kocevski et al. (2012) Kocevski, D. D., Faber, S. M., Mozena, M., et al. 2012, ApJ, 744, 148. doi:10.1088/0004-637X/744/2/148
  * Kocevski et al. (2015) Kocevski, D. D., Brightman, M., Nandra, K., et al. 2015, ApJ, 814, 104. doi:10.1088/0004-637X/814/2/104
  * Kormendy & Kennicutt (2004) Kormendy, J., & Kennicutt, R. C., Jr. 2004, ARA&A, 42, 603. doi:10.1146/annurev.astro.42.053102.134024
  * Lintott et al. (2008) Lintott, Chris J., et al. 2008, MNRAS, 389, 3. doi:10.1111/j.1365-2966.2008.13689.x
  * Lotz et al. (2004) Lotz, J. M., Primack, J., Madau, P., 2004, AJ, 128, 163, doi=10.1086/421849
  * Lotz et al. (2010) Lotz, J. M.. et al. 2010, MNRAS, 404. doi=10.1111/j.1365-2966.2010.16268.x
  * Lotz et al. (2008) Lotz, J. M., et al. 2008, ApJ, 672, 1. doi=10.1086/523659
  * Luo et al. (2017) Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2. doi:10.3847/1538-4365/228/1/2
  * Lyu et al. (2022) Lyu, J., Alberts, S., Rieke, G. H., et al. 2022, ApJ, 941, 191. doi:10.3847/1538-4357/ac9e5d
  * Lyu et al. (2023) Lyu, J., Alberts, S., Rieke, G. H., et al. 2023, arXiv:2310.12330. doi:10.48550/arXiv.2310.12330
  * McAlpine et al. (2020) McAlpine, S., et al. 2020, 494, 4. doi=10.1093/mnras/staa1123
  * Pawlik et al. (2016) Pawlik, M. M., Wild, V., Walcher, C. J., et al. 2016, MNRAS, 456, 3032. doi:10.1093/mnras/stv2878
  * Pierce et al. (2023) Pierce, J. C. S., Tadhunter, C., Ramos Almeida, C., et al. 2023, MNRAS, 522, 1736. doi:10.1093/mnras/stad455
  * Reza et al. (2021) Reza, M., et al. 2021, Astronomy and Computing, 37, 100492. doi:10.1016/j.ascom.2021.100492
  * Rieke et al. (2023) Rieke, M. J., Robertson, B. E., Tacchella, S., et al. 2023, arXiv:2306.02466. doi:10.48550/arXiv.2306.02466
  * Robertson et al. (2006) Robertson, B., et al., ApJ, 645, 986. doi=10.1086/504412
  * Rodriguez-Gomez et al. (2019) Rodriguez-Gomez, V., et al. 2019, MNRAS, 483, 4140. doi:10.1093/mnras/sty3345
  * Sanders & Mirabel (1996) Sanders, D. B. & Mirabel, I. F. 1996, ARA&A, 34, 749. doi:10.1146/annurev.astro.34.1.749
  * Schade et al. (1995) Schade, D., Lilly, S. J., Crampton, D., Hammer, F., LeFevre, O., & Tresse, L. 1995, ApJ, 451, L1. doi=10.1086/309677
  * Shah et al. (2020) Shah, E. A., Kartaltepe, J. S., Magagnoli, C. T., et al. 2020, ApJ, 904, 107. doi:10.3847/1538-4357/abbf59
  * Snyder et al. (2015) Snyder, G. F., et al. 2015, MNRAS, 451, 4290. doi=10.1093/mnras/stv1231
  * Suess et al. (2022) Suess, K. A., et al. 2022, ApJ, 937, L33. doi=10.3847/2041-8213/ac8e06
  * Takamiya et al. (1999) Takamiya, M. 1999, ApJS, 122, 1. doi=10.1086/313216
  * Treister et al. (2012) Treister, E., Schawinski, K., Urry, C. M., et al. 2012, ApJ, 758, L39. doi:10.1088/2041-8205/758/2/L39
  * Tremaine et al. (2002) Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, ApJ, 574, 740. doi:10.1086/341002
  * Urbano-Mayorgas et al. (2019) Urbano-Mayorgas, J. J., Villar-Martín, M., & Buitrago, F. 2019, Highlights on Spanish Astrophysics X, 232
  * Villforth (2023) Villforth, C. 2023, The Open Journal of Astrophysics, 6, 34. doi:10.21105/astro.2309.03276
  * Walmsley et al. (2002) Walmsley, M., et al. 2022, MNRAS, 509, 3. doi:10.1093/mnras/stab2093
  * Williams et al. (2023) Williams, C.C., et al. 2023, ApJ, 268, 64. 10.3847/1538-4365/acf130
  * Wu (1999) Wu, K. 1999, Ph.D. thesis, Univ. California, Santa Cruz
  * Xu et al. (2015) Xu, L., Rieke, G. H., Egami, E., et al. 2015, ApJ, 808, 159. doi:10.1088/0004-637X/808/2/159
  * Zakamska et al. (2004) Zakamska, N. L., Strauss, M. A., Heckman, T. M., et al. 2004, AJ, 128, 1002. doi:10.1086/423220
  * Zakamska et al. (2019) Zakamska, N. L., Sun, A.-L., Strauss, M. A., et al. 2019, MNRAS, 489, 497. doi:10.1093/mnras/stz2071
  * Zhang et al. (2016) Zhang, Z., Shi, Y., Rieke, G. H., et al. 2016, ApJ, 819, L27. doi:10.3847/2041-8205/819/2/L27