# Obscuration beyond the nucleus: infrared quasars can be buried in extreme
compact starbursts

Carolina Andonie ,1 David M. Alexander,1 Claire Greenwell,1 Annagrazia
Puglisi,2,1 Brivael Laloux,3,1 Alba V. Alonso-Tetilla,2 Gabriela Calistro
Rivera, 4 Chris Harrison,5 Ryan C. Hickox,6 Melanie Kaasinen,4 Andrea Lapi,7
Iván E. López,8,9 Grayson Petter,6 Cristina Ramos Almeida,10,11 David J.
Rosario,5 Francesco Shankar,2 and Carolin Villforth12
1 Centre for Extragalactic Astronomy, Department of Physics, Durham
University, Durham, DH1 3LE, UK
2 School of Physics and Astronomy, University of Southampton, Highfield SO17
1BJ, UK
3Institute for Astronomy & Astrophysics, National Observatory of Athens, V.
Paulou & I. Metaxa 11532, Athens, Greece
4European Southern Observatory (ESO), Karl-Schwarzschild-Straße 2, 85748
Garching bei München, Germany
5School of Mathematics, Statistics and Physics, Newcastle University,
Newcastle upon Tyne, NE1 7RU, UK
6 Department of Physics and Astronomy, Dartmouth College, 6127 Wilder
Laboratory, Hanover, NH 03755, USA
7 SISSA, Via Bonomea 265, I-34136 Trieste, Italy
8Dipartimento di Fisica e Astronomia "Augusto Righi", Università di Bologna,
Via Gobetti 93/2, 40129 Bologna, Italy
9 INAF – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via
Gobetti 93/3, 40129 Bologna, Italy
10 Instituto de Astrofísica de Canarias, Calle Vía Láctea, s/n, 38205 La
Laguna, Tenerife, Spain
11 Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna,
Tenerife, Spain
12 Department of Physics, University of Bath, Claverton Down, Bath BA2 7AY, UK
E-mail<EMAIL_ADDRESS>Fellow

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In the standard quasar model, the accretion disk obscuration is due to the
canonical dusty torus. Here, we argue that a substantial part of the quasar
obscuration can come from the interstellar medium (ISM) when the quasars are
embedded in compact starbursts. We use an obscuration-unbiased sample of 578
infrared (IR) quasars at $z\approx 1-3$ and archival ALMA submillimeter host
galaxy sizes to investigate the ISM contribution to the quasar obscuration. We
calculate SFR and ISM column densities for the IR quasars and a control sample
of submillimeter galaxies (SMGs) not hosting quasar activity and show that:
(1) the quasar obscured fraction is constant up to $\rm SFR\approx
300\>M_{\odot}\>yr^{-1}$, and then increases towards higher SFR, suggesting
that the ISM obscuration plays a significant role in starburst host galaxies,
and (2) at $\rm SFR\gtrsim 300\>M_{\odot}\>yr^{-1}$, the SMGs and IR quasars
have similarly compact submillimeter sizes ($R_{\rm e}\approx 0.5-3\rm\>kpc$)
and, consequently, the ISM can heavily obscure the quasar, even reaching
Compton-thick ($N_{\rm H}>10^{24}\rm\>cm^{-2}$) levels in extreme cases. Based
on our results, we infer that $\approx 10-30\%$ of the IR quasars with $\rm
SFR\gtrsim 300\>M_{\odot}\>yr^{-1}$ are obscured solely by the ISM.

###### keywords:

quasars: general – galaxies: starburst – submillimetre: galaxies

††pubyear: 2023††pagerange: Obscuration beyond the nucleus: infrared quasars
can be buried in extreme compact starbursts–Obscuration beyond the nucleus:
infrared quasars can be buried in extreme compact starbursts

##  1 Introduction

Quasars are the most luminous subset of the Active Galactic Nucleus (AGN)
population and represent the most rapid phase of supermassive black hole
(SMBH) growth. Strong observational evidence supports a connection between
SMBH accretion and the evolution of their host galaxies (Alexander & Hickox,
2012; Kormendy & Ho, 2013; Madau & Dickinson, 2014), which is expected to
occur in an early obscured phase (Sanders et al., 1988; Hopkins et al., 2008;
Lapi et al., 2014; Hickox & Alexander, 2018), with large reservoirs of gas and
dust feeding the SMBH and the galaxy. In the last few years, there has been
growing evidence favouring this synchronized evolution due to the discovery of
fundamental differences in the host galaxy and nuclear properties between
obscured ($N_{\rm H}>10^{22}\rm\>cm^{-2}$) and unobscured quasars ($N_{\rm
H}<10^{22}\rm\>cm^{-2}$), which cannot be easily explained by the standard AGN
orientation model (Antonucci, 1993; Urry & Padovani, 1995); i.e., the observed
differences cannot be only due solely to the orientation of a geometrically
and optically thick dusty torus (Ramos Almeida & Ricci, 2017). A number of
recent studies have found systematic differences between obscured and
unobscured quasars in the clustering properties (Hickox et al., 2011; DiPompeo
et al., 2014; Donoso et al., 2014; Villarroel & Korn, 2014; Powell et al.,
2018; Petter et al., 2023), star-formation (SF) properties (Chen et al., 2015;
Andonie et al., 2022), and galaxy companion environments (Satyapal et al.,
2014; Weston et al., 2017; Dougherty et al., 2023). Combined, these results
suggest that obscured quasars represent a distinct early growth phase closely
connected to a period of enhanced SF.

Recently, Andonie et al. (2022) (hereafter A22) found that obscured quasars
have star-formation rates (SFR) $\approx 3$ times higher than unobscured
quasars, a result inconsistent with the standard AGN orientation model. Most
striking is the finding that an excess of obscured quasars reside in starburst
galaxies with $\rm SFR>300\>M_{\odot}\>yr^{-1}$, providing evidence that the
host galaxy may contribute to the total obscuration of the quasar. A number of
previous studies have investigated the host galaxy contribution to the total
AGN obscuration (e.g., Gilli et al., 2007; Buchner et al., 2017; Circosta et
al., 2019; Gilli et al., 2022; Silverman et al., 2023), but our current study
is the first to directly constrain the host galaxy obscuration as a function
of the SFR. To achieve this we use observations from the Atacama Large
Millimeter/submillimeter Array (ALMA) of a sub sample of the A22 IR quasars.
We use these ALMA data to improve the SFR constraints and explore the host
galaxy submillimeter (submm) sizes to assess the contribution of the
interstellar medium (ISM) to the total quasar obscuration. We also analyse a
control sample of submm galaxies (SMGs), which do not host an AGN.

In this work, we adopt a concordance cosmology (Hinshaw et al., 2013), and a
Chabrier initial mass function (Chabrier, 2003).

##  2 Observations and analysis

###  2.1 SED fitting, photometric catalogues, and the sample

Our parent sample is taken from the quasar catalogue published by A22 from the
Cosmic Evolution Survey (COSMOS, Scoville et al., 2007). A22 identified a
complete sample of 578 infrared (IR) quasars ($L_{\rm AGN,8-1000\>\mu
m}>10^{45}\rm\>erg\>s^{-1}$) at $z<3$, with minimal obscuration bias, using
detailed UV-to-far IR spectral energy distribution (SED) fitting. We find that
24 IR quasars (hereafter the submm-IR quasars) have archival ALMA host galaxy
submm size measurements (see Section 2.2 for details). In this study we re-fit
the rest-frame $0.3-1000\rm\>\mu m$ SEDs of the submm-IR quasars but now
including the ALMA Band 6/7 photometry from the latest version of the A3COSMOS
photometric catalogue (version 20200310, Liu et al., 2019). Our SED fitting is
performed using the multicomponent Bayesian SED fitting code fortesfit
111https://github.com/vikalibrate/FortesFit (Rosario, 2019), and the same
models and SED fitting approach of A22. The SED modelling includes the
unabsorbed stellar emission, the UV emission from the accretion disk, IR
emission from the AGN, and the dust obscured star-formation (for more details,
see Section 2 and 3 of A22). As a comparison, we also fit a control sample of
all the 85 SMGs with reliable host galaxy sizes in COSMOS reported by Fujimoto
et al. (2017), using the same SED fitting approach as for the IR quasars. We
exclude from the SMG sample all the sources either identified as AGN by our
SED fitting, detected in the Chandra COSMOS Legacy survey (Civano et al.,
2016, see Section 2.3), or part of the AGN optical spectroscopic catalogue of
Rosario et al. (2013), ending up with 78 sources. For details on the AGN
identification process and SFR upper-limits estimation, we refer the readers
to Section 3.3 of A22. The SED best-fitting parameters are reported in Table 1
of the online Supplementary Material.

Figure 1 depicts the host galaxy properties of the IR quasar sample of A22,
the submm-IR quasars, and our control sample of SMGs, and also reports the
Kolmogorov Smirnov (KS) test p-value between the submm-IR quasars and SMGs. We
plot the 50th percentile from the posterior distribution for the SED best-
fitting parameters. Panels (a) and (d) show that the three samples have
consistent redshifts and stellar masses. Similarly, panel (c) shows that the
SFRs of the IR quasars span a wide range of values, and that the submm-IR
quasars are a good representation of the IR quasars with the largest SFRs; the
distributions of the redshifts, SFRs, $\rm M_{\star}$, and $L_{\rm AGN,IR}$ of
the submm-IR quasars and the IR quasars with $\rm SFR>100\>M_{\odot}\>yr^{-1}$
always have p-values of $>0.3$. However, we note that the majority of the IR
quasars have low SFRs ($<10\rm\>M_{\odot}\>yr^{-1}$), and we can only provide
an SFR upper limit in these cases; hence, those sources are excluded from the
KS test analysis. We find that the distributions of the SFRs and stellar
masses of the submm-IR quasars (mean and standard deviation values of $\rm
SFR_{submm-IR\>quasars}=610\pm 400\>M_{\odot}\>yr^{-1}$ and
$M_{\rm\star,submm-IR\>quasars}=(8\pm 4)\cdot 10^{10}\rm\>M_{\odot}$) and SMGs
(mean and standard deviation values of $\rm SFR_{SMGs}=700\pm
550\>M_{\odot}\>yr^{-1}$ and $M_{\rm\star,SMGs}=(8.6\pm 4.2)\cdot
10^{10}\rm\>M_{\odot}$) are also similar (in both cases the KS p-value $>0.2$)
and place both populations on the starburst region of the main-sequence of
star-forming galaxies, following Whitaker et al. (2012) and Aird et al.
(2019).

Figure 1: Normalized distributions for the full IR quasar sample (578 sources)
of A22 (green dotted), the 24 IR quasars with ALMA sizes (red solid), and the
control sample of 78 SMGs (black solid). In 3 panels, the K-S $\rm p-value$
between the IR quasars with sizes and the SMGs is quoted. (a): redshifts; (b):
$8-1000\>\rm\mu m$ AGN luminosity (c): SFRs; (d) stellar masses.

###  2.2 ALMA data

We use submm (i.e., dust continuum) measurements from the Band 6 or 7 ($\sim
1\rm\>mm$, whichever is available) ALMA data of Fujimoto et al. (2017), Chang
et al. (2020), and Lamperti et al. (2021), in addition to measurements from a
combination of CO(5-4), CO(4-3), CO(2-1), and $1.1-1.3\rm\>mm$ continuum
observations from Puglisi et al. (2019, 2021)222There is no bias in the sizes
from Puglisi et al. (2019) as they showed that the CO(5-4), CO(4-3) are
consistent with $\sim 1\rm\>mm$ sizes for their sample.. Overall, there are 24
IR quasars with reliable host galaxy sizes from ALMA observations with $\rm
S/N>5$, and all of them used a circularized Gaussian profile, or a consistent
geometry333Lamperti et al. (2021) fitted different models to measure sizes,
and for the three overlapping sources from our sample, the best-fitting model
sizes are consistent within $1\sigma$ with those derived using a circularised
Gaussian profile., to fit the data. Here, we define the submm effective radius
as $R_{e}=\rm FWHM/2$. The size measurements and references for each source in
our sample are reported in Table 1 of the Supplementary Material.

The angular resolution of the observations used to calculate the sizes covers
$\theta_{\rm ALMA}\approx$ 0.′′18-0.′′9 ($\approx 1.5-7\rm\>kpc$ at $z=1-3$),
including 6 out of 24 objects with $\theta_{\rm ALMA}<0.$′′3\. We note that we
do not expect any of these observations to over-resolve the galaxy, given that
most of the dust continuum emission is typically concentrated into a central
compact component of $\sim$0.′′$15$ (Gullberg et al., 2019).

The AGN-heated dust could affect the host galaxy size measurements if its
contribution to the ALMA Band 6/7 flux is significant. However, on the basis
of our SED decomposition (see Section 2.1 for details), the AGN contribution
to the total rest-frame $850\rm\>\mu m$ luminosity is $>10\%$ for only three
sources and $<5\%$ for 19 sources. For the three sources with the slightly
largest quasar contribution, the host galaxy sizes are $>1.3\rm\>kpc$, which
is larger than the average size of the submm-IR quasars (see Section 3.2 for
details). Hence, we conclude that the quasar emission has a minimal impact on
the host galaxy submm size measurements.

###  2.3 X-ray spectral constraints

In the A22 sample, 63% (365/578) of the IR quasars are X-ray detected in the
Chandra COSMOS Legacy survey (Civano et al., 2016), while the other 37% remain
undetected. In this work, we adopted the X-ray spectral fitting constraints
from Laloux et al. (2023), unlike A22, who used the measurements from Marchesi
et al. (2016) and Lanzuisi et al. (2018). For the overlapping sources, both
catalogues are in broad agreement with the absorbing column densities for
sources with $N_{\rm H}<10^{23}\rm\>cm^{-2}$, and the classification between
obscured/unobscured sources does not change. However, Laloux et al. (2023)
used a Bayesian approach that robustly constrains the X-ray spectral physical
parameters for sources with $<30$ counts, which lie below the lower count
threshold adopted in Marchesi et al. (2016). Analyzing such faint sources is
fundamental in this study, justifying our choice of X-ray catalogue.

Briefly, in Laloux et al. (2023), each source is fitted using the BXA package
(Buchner et al., 2014) with the UXCLUMPY X-ray model (Buchner et al., 2019).
The output of the fit is a parameter posterior distribution that fully
encapsulates the measurement uncertainties. We calculate the probability that
an individual source is obscured by taking the fraction of the line of sight
(l.o.s) column density posterior distribution above $N_{\rm
H,X\>(los)}=10^{22}\rm\>cm^{-2}$, $f_{N_{\rm H}>22}$. In the case of the X-ray
undetected sources, we set $f_{N_{\rm H}>22}=1$ following the results of A22,
who used a combination of X-ray stacking and $L_{\rm 2-10\rm keV}$ – $L_{\rm
6\>\mu m}$ analyses to demonstrate that the majority ($>$ 85%) of the X-ray
undetected sources are heavily obscured by $N_{\rm
H,X\>(los)}>10^{23}\rm\>cm^{-2}$, and many are likely Compton-thick ($N_{\rm
H,X\>(los)}>10^{24}\rm\>cm^{-2}$, CT); see Section 4.3, and Figures 6 and 7 of
A22 for details. We note that we do not expect the stacked X-ray signal to be
contaminated by SF processes444Adopting Ranalli et al. (2003), the X-ray
luminosity due to $\rm SFR>300\>M_{\odot}\>yr^{-1}$ ranges $L_{\rm
2-10\>keV,SF}=(2-9)\cdot 10^{42}\>\rm erg\>s^{-1}$, which is $<5\%$ of the AGN
X-ray luminosity calculated from the $L_{\rm 6\>\mu m}$ following Stern
(2015)., but to further confirm that the X-ray undetected IR quasars with
large SFRs are obscured, we take the same approach of A22 and stack all the
X-ray undetected quasars with $\rm SFR>300\>M_{\odot}\>yr^{-1}$ using the
python code StackFast555The StackFast code is available here
https://www.tonima-ananna.com/research (see section 4.3 of that paper). As in
A22, we detect a strong stacked X-ray signal. The count rates in the soft
($0.5-2\rm\>keV$) and hard ($2-7\rm\>keV$) bands are $\rm S=(9.95\pm
1.84)\times 10^{-6}\>counts\>s^{-1}$ and $\rm H=(7.73\pm 1.82)\times
10^{-6}\>counts\>s^{-1}$, which is equivalent to a hardness ratio of $\rm
HR=\frac{H-S}{H+S}=0.13\pm 0.15$. At the mean redshift of the sample
($z=2.2\pm 0.5$), the HR indicates that, on average, these IR quasars are
obscured by $N_{\rm H}>10^{23}\rm\>cm^{-2}$ (see Figure 6 in A22) and,
consequently, that the majority of the X-ray undetected IR quasars with $\rm
SFR>300\>M_{\odot}\>yr^{-1}$ are X-ray undetected because the X-ray emission
was absorbed by heavily-obscured column densities666The X-ray undetected IR
quasars with $\rm SFR<300\>M_{\odot}\>yr^{-1}$ have a $\rm HR=-0.1\pm 0.2$ and
$z\approx 2$, which also implies $N_{\rm H}>10^{23}\rm\>cm^{-2}$..

For the nine submm-IR quasars that are undetected in the X-ray band, we
calculated $N_{\rm H,X\>(los)}$ lower limits. For this, we adopt the Stern
(2015) relationship to convert the rest-frame 6$\mu m$ AGN luminosity
calculated from the SED fits into the intrinsic $2-10\>\rm keV$ luminosity. We
then utilize the Chandra sensitivity maps generated in Laloux et al. (2023)
and use UXCLUMPY to compute the minimum $N_{\rm H,X\>(los)}$ that each source
needs to be undetected in the X-ray band at its redshift, which corresponds to
the value for what the detection probability drops below 90%. The $N_{\rm
H,X\>(los)}$ ranges $N_{\rm H,X\>(los)}=(0.3-5)\cdot 10^{23}\rm\>cm^{-2}$ and
are reported in Table 1 of the Supplementary Material.

##  3 Results

###  3.1 Connection between the obscured fraction and the star-formation rate

The SFR enhancement in obscured quasars found in A22 suggests that dust-
obscured SF from the host galaxy may cause at least part of the obscuration
towards the quasars. In this section, we demonstrate the connection between
obscuration and SFR.

Panel (a) of Figure 2 shows the obscured quasar fraction as a function of the
SFR for the entire IR quasar sample of A22. The obscured fraction values and
uncertainties correspond to the mean and standard error of the $f_{N_{\rm
H}>22}$ distribution at each SFR bin (see Section 2.3 for details). Similarly,
the x-axis shows the mean and standard deviation of the SFR distribution at
each bin, calculated using
PosteriorStacker777https://github.com/JohannesBuchner/PosteriorStacker. This
code takes into account the upper limits by modelling samples of the posterior
distributions using a Gaussian model.. The total obscured fraction of the
sample is $68\pm 2\%$. However, we find that the obscured fraction has a clear
dependence with SFR, being broadly constant ($\approx 65-70\%$) for $\rm
SFR<300\>M_{\odot}\>yr^{-1}$, increasing to $\approx 80\%$ for $\rm
SFR=300-900\>M_{\odot}\>yr^{-1}$, and reaching $\approx 90\%$ for $\rm
SFR>900\>M_{\odot}\>yr^{-1}$. Panel (b) clearly illustrates this increase by
showing that the obscured fraction at $\rm SFR>300\>M_{\odot}\>yr^{-1}$ is
between $1.2\pm 0.1$ and $1.4\pm 0.1$ times larger than the obscured fraction
at the lowest SFR bin. We note that the number of sources contributing to each
SFR bin is not constant: the majority (311/578) reside in the lowest SFR bin,
while the two largest SFR bins are comprised of $<30$ sources. However, we
still observe the same trend when taking three equally-populated bins
($\approx$ 70 sources) for the quasars with $\rm SFR>50\>M_{\odot}\>yr^{-1}$
(see the grey area on the graph). The observed increase in the obscured
fraction at $\rm SFR>300\>M_{\odot}\>yr^{-1}$ suggests that at those SF
levels, the host galaxy starts to significantly contribute to the total
obscuration of the system.

Panel (c) of Figure 2 plots the mean and standard deviation of the redshift
distributions at each SFR bin, with redshifts ranging $\sim 1.7-2.4$. A
general trend of a modest redshift increase with SFR is seen. However, we find
that the change in the obscured fraction with SFR is unlikely to be explained
by evolution in the intrinsic obscured AGN fraction since this is expected to
increase by only up to $\sim 1.1\times$ between $\rm z=[1.5,2.5]$ (e.g., Gilli
et al., 2022; Laloux et al., 2023). We also note that the stellar mass is
broadly constant (median values ranging $M_{\rm\star}\approx(5-7)\cdot
10^{10}\rm\>M_{\odot}$) across the different SFR bins. Hence, the increment of
the obscured fraction with SFR is unlikely to be due to a redshift or stellar
mass evolution.

We note from Figure 2 that the X-ray detection fraction does not significantly
vary with the SFR, except at the highest SFR bin, dropping from $\approx
0.55-0.70$ to $\approx 0.25$; therefore, the obscured fraction at the largest
SFR bin is dominated by X-ray undetected sources. In Section 2.3, we provide
strong evidence that the majority of the X-ray undetected sources are obscured
by $N_{\rm H}>10^{23}\rm\>cm^{-2}$; hence, we set $f_{N_{\rm H}>22,\rm
X\>und}=1$ for the X-ray undetected sources. However, we still observe the
same trend between the obscured fraction and the SFR when we conservatively
consider that only 75% of the X-ray undetected quasars are obscured by $N_{\rm
H}>10^{22}\rm\>cm^{-2}$ (i.e., $f_{N_{\rm H}>22,\rm X\>und}=0.75$). We finally
note that the same trend is seen when only including the X-ray detected
quasars for which we have directly measured the obscured fraction from X-ray
spectral fitting (magenta circles in Figure 2); however, the obscured fraction
is very uncertain at the largest SFR bin due to it only containing six
sources.

Figure 2: Panel (a): Obscured quasar fraction versus SFR for the full IR
quasar sample of A22. Each obscured fraction value corresponds to the mean and
standard error of $f_{N_{\rm H}>22}$ at each SFR bin, defined as $\rm
SFR=[<100,100-200,200-300,300-600,600-900,>900]\rm\>M_{\odot}\>yr^{-1}$. Black
points correspond to the IR quasars and magenta points are the X-ray detected
quasars only. The grey area represents 1$\sigma$ limits for the obscured
fraction of three equal-size SFR bins for sources with $\rm
SFR>50\>M_{\odot}\>yr^{-1}$. The black and blue text on top of each point
corresponds to the sample size and X-ray detection fraction for each SFR bin,
respectively. The x-axes are the mean and standard deviation of the SFRs in
each bin. Panel (b): Fractional increase of the obscured fraction with respect
to the obscured fraction at the lowest SFR bin. Panel (c): mean and standard
deviation of the redshift distribution at each SFR bin.

###  3.2 Host galaxy obscuration of highly star-forming quasars

The results of Section 3.1 suggest that the host galaxy can significantly
contribute to the total quasar obscuration in highly star-forming systems.
This section aims to quantify the host-galaxy ISM contribution to the total
quasar obscuration using ALMA submm size constraints for the 24 submm-IR
quasars.

We calculate the column density associated with the gas in the host galaxy,
taking an approach already used in previous works (e.g., Circosta et al.,
2019). To do this, we first calculate the gas mass of the host galaxy, which
is estimated as $M_{\rm gas}=M_{\rm H_{2}}+M_{\rm HI}$. The molecular gas mass
is calculated from the $\rm CO(1-0)$ luminosity as $M_{\rm H2}=\alpha_{\rm
CO}L^{\prime}_{\rm CO}$, where $\alpha_{\rm CO}\approx 0.8\rm\>M_{\odot}/(\rm
K\>km\>s^{-1}\>pc^{2})$ for starburst galaxies (e.g., Magdis et al., 2012;
Carilli & Walter, 2013; CalistroRivera et al., 2018; Birkin et al., 2021). We
calculate $L^{\prime}_{\rm CO}$ using the well-measured empirical relationship
to the rest-frame $850\rm\>\mu m$ luminosity: $L^{\prime}_{\rm CO}[\rm
K\>km\>s^{-1}\>pc^{2}]=3.02\times 10^{-21}\>L_{\rm 850\mu m}[\rm
erg\>s^{-1}\>Hz^{-1}]$ (Scoville et al., 2016). We estimate the rest-frame
$L_{\rm 850\mu m}$ from the best-fitting SF component in our SED fitting,
removing any AGN contribution; when the SFR is unconstrained, we can only
provide an $L^{\prime}_{\rm CO}$ upper-limit, which we calculate using the
99th-percentile of $L_{\rm 850\mu m}$ distribution. We note that this relation
is valid for galaxies with stellar masses $M_{\star}>\rm 2\cdot 10^{10}\>\rm
M_{\odot}$, which is well within the range of our submm-IR quasars and SMGs
(see Figure 1). We then calculate the atomic gas mass following theoretical
models, which predicts $M_{\rm HI}\approx[2-5]\cdot M_{\rm H_{2}}$ (Lagos et
al., 2011; Fu et al., 2012) for the redshifts and stellar masses of our
quasars; hence, we adopt a mean value of $M_{\rm HI}=2.5\cdot M_{\rm H_{2}}$.
The gas masses of the IR quasars and SMGs range $M_{\rm gas,submm-
IR\>quasars}=(0.04-13)\rm\>\times 10^{10}\>M_{\odot}$ and $M_{\rm
gas,SMGs}=(0.3-20)\rm\>\times 10^{10}\>M_{\odot}$, respectively.

To calculate the average ISM column density, we adopt a uniform spherical
geometry. Assuming that half of the gas resides inside $R_{e}$, the numeric
density of hydrogen atoms is calculated as $n_{\rm H}=1/2M_{\rm gas}/(m_{\rm
H}\>4/3\pi R_{\rm e}^{3})$, where $m_{\rm H}$ is the Hydrogen mass. Then,
$\langle N_{\rm H,ISM}\rangle=\int_{0}^{R_{e}}n_{\rm H}\>\rm dr=1/2M_{\rm
gas}/(m_{\rm H}\>4/3\pi R_{\rm e}^{2})$. All the quantities calculated in this
section are reported in Table 1 of the Supplementary Material.

Figure 3: Submm host galaxy size ($R_{\rm e}$) against the SFR, coloured by
the mean ISM column density ($\langle N_{\rm H,ISM}\rangle$). Stars correspond
to SMGs, and circles correspond to submm-IR quasars. SFR upper limits are
plotted at 99% confidence. Empty circles are quasars with an upper limit on
$R_{\rm e}$ and SFR; hence, no ISM column density can be calculated. The
dotted lines represent linear regressions between $R_{\rm e}$ and SFR for
sources with $\langle N_{\rm H,ISM}\rangle<3\cdot 10^{23}\>\rm cm^{-2}$
(green), $\langle N_{\rm H,ISM}\rangle=[3,10]\cdot 10^{23}\>\rm cm^{-2}$
(orange), and $\langle N_{\rm H,ISM}\rangle>10^{24}\>\rm cm^{-2}$ (red). The
text next to each dotted line is the mean $\langle N_{\rm H,ISM}\rangle$ in
each bin. The histograms on the right side of the plot show the size
distributions of the submm-IR quasars (black) and SMGs (red).

Figure 3 shows the effective radius versus SFR of our sources colour-coded by
$\langle N_{\rm H,ISM}\rangle$. We find that the submm sizes of the submm-IR
quasars and SMGs are in broad agreement, with values ranging $R_{e}\approx
0.5-2.7\rm\>kpc$ (mean and standard error of $R_{\rm e,submm-
IR\>quasars}=1.1\pm 0.6\rm\>kpc$ and $R_{\rm e,SMG}=0.95\pm 0.43\rm\>kpc$).
Additionally, the ISM obscuration values of the submm-IR quasars range
$\langle N_{\rm H,ISM}\rangle_{\rm submm-IRquasars}=(0.01-4)\times
10^{24}\>\rm cm^{-2}$, with a mean value consistent with that of the SMGs
(mean and standard error of $\langle N_{\rm H,ISM}\rangle_{\rm submm-
IRquasars}=(0.8\pm 0.1)\times 10^{24}\>\rm cm^{-2}$ and $\langle N_{\rm
H,ISM}\rangle_{\rm SMGs}=(1.5\pm 3.5)\times 10^{24}\>\rm cm^{-2}$), confirming
that the only significant difference between the two populations is the
presence or absence of quasar activity (see also Figure 1).

Figure 3 also shows the results of the linear regressions between $R_{\rm e}$
and SFR for three different ISM obscuration bins, combining the submm-IR
quasars and SMGs. We find that CT ISM obscurations only occur in very compact
host galaxies with $R_{\rm e}\lesssim 0.8\rm\>kpc$. For larger host galaxy
sizes, even for extremely high SFRs ($>1000\>\rm M_{\odot}\>yr^{-1}$), the ISM
obscuration does not exceed $\langle N_{\rm H,ISM}\rangle\approx 5\cdot
10^{23}\rm\>cm^{-2}$. Similarly, galaxies with $R_{\rm e}\sim 2\rm\>kpc$
typically have $\langle N_{\rm H,ISM}\rangle\lesssim 10^{23}\rm\>cm^{-2}$,
irrespective of the SFR value. These results indicate that the ISM can have an
important contribution to the total obscuration of the system in extremely
compact star-forming galaxies and that its contribution increases towards
smaller host galaxy sizes.

Some studies suggest that the dust continuum sizes from submm observations can
underestimate the gas sizes by a factor of $\approx 1.3-2$ (e.g.,
CalistroRivera et al., 2018; Popping et al., 2022), which would imply an
overestimation of our ISM column densities. We find that if we increase the
host galaxy sizes by a factor of two, the ISM column densities decrease by a
factor of four but not changing the main results of the paper, since the
majority of the submm-IR quasars have inferred $\langle N_{\rm
H,ISM}\rangle>10^{23}\rm\>cm^{-2}$. Furthermore, we note that modelling the
ISM with a spherical geometry is a simplification; for example, there are
eight submm-IR quasars in the spectroscopic catalogue of broad-lines (BL)
quasars from Rosario et al. (2013), indicating that we can have an unobscured
view of the accretion disk due to the patchiness and clumpiness of the
obscuring material. If we instead adopt a disk with a covering factor of 50%,
our ISM column densities would increase by a factor of 2.

The X-ray spectra provide direct constraints of the l.o.s column densities
($N_{\rm H,X(los)}$) towards the central SMBH. Although we cannot directly
compare these values with the estimates of the ISM column densities (which are
instead average values), it is insightful to compare both distributions since
they will give us some indication of the ISM contribution to the overall
quasar obscuration. Figure 4 shows the distributions of $\langle N_{\rm
H,ISM}\rangle$ and $N_{\rm H,X(los)}$. We exclude the BL submm-IR quasars
since we know that, by definition, they have a preferential and
unrepresentative unobscured l.o.s view. We also calculate the mean $\langle
N_{\rm H,ISM}\rangle$ and $N_{\rm H,X(los)}$ of our sample by adopting uniform
distributions when we only have limits: for the ISM column density upper
limits, we set the lowest possible value to $\log\langle N_{\rm H,ISM}/\rm
cm^{2}\rangle=21.7$, following the results of this section; while for X-ray
undetected sources, we set the highest possible value of $\log(N_{\rm
H,X(los)}/\rm cm^{2})=24.7$, following Alexander et al. (2005). We find that
the mean ISM and X-ray column density values of our sample are $\log\langle
N_{\rm H,ISM}/\rm cm^{2}\rangle=23.4\pm 0.1$ and $\log(N_{\rm H,X(los)}/\rm
cm^{2})=23.6\pm 0.1$, respectively, implying that the ISM is comparable to
(and potentially dominates) the overall obscuration towards the quasar in many
systems. As expected, the BL submm-IR quasars have a substantially lower X-ray
column density ($\log(N_{\rm H,X(los)}/\rm cm^{2})=23.1\pm 0.1$), despite
their $\langle N_{\rm H,ISM}\rangle$ remaining consistent to the overall
sample, demonstrating that their l.o.s column density is not representative of
the typical column density towards these sources. We caution that the ALMA
observations of AGN in COSMOS are biased toward X-ray sources; hence, the
fraction of BL submm-IR quasars (8/24) might be overestimated.

Figure 4: Column-density distributions of the submm-IR quasars. Top panel:
distribution of the mean ISM column density ($\langle N_{\rm H,ISM}\rangle$).
The green histogram shows sources with a $\langle N_{\rm H,ISM}\rangle$ upper
limit, due to an unconstrained star formation component in the SED. Bottom
panel: distribution of the l.o.s X-ray column density ($N_{\rm H,X\>(los)}$).
The green histogram represents $N_{\rm H,X\>(los)}$ lower limits for the X-ray
undetected sources. The red dotted lines show the mean $N_{\rm H}$ of each
distribution; see section 3.2.

##  4 Discussion and Conclusions

Overall, in this letter we have two key results: (1) a relationship between
obscuration and SFR for IR quasars, showing that the host galaxy is likely to
make a significant contribution to the increase in the obscured quasar
fraction at $\rm SFR\gtrsim 300\>M_{\odot}\>yr^{-1}$ (see Figure 2); and (2)
estimates of the ISM column densities which shows that quasars can be heavily
obscured by their host galaxies (see Figure 3). The first result is in
agreement with Chen et al. (2015), who found that the obscured quasar fraction
increases from 0.3 to 0.7 for $\rm SFR=200-630\rm\>M_{\odot}\>yr^{-1}$. While
our second result is in agreement with previous works which found that AGN
host galaxies can obscure the central SMBH (e.g., Buchner et al., 2017;
Silverman et al., 2023) and can even be CT in some extreme cases (e.g.,
Circosta et al., 2019; Gilli et al., 2022). However, our study combines these
two results for the first time by exploring the connection between the host
galaxy obscuration and the SFR in IR quasars. In this final section, we
discuss our results and their implications for quasar models.

###  4.1 ISM contribution to the obscured quasar fraction

Section 3.2 shows that submm-IR quasars and SMGs can be heavily obscured by
the ISM and even CT for the most compact sources. Such high levels of host
galaxy obscuration are not unexpected. Indeed, SMGs are known to have very
high gas surface densities, with typical values of
$\rm\Sigma_{gas}\approx(10^{2}-10^{4})\>M_{\odot}\>pc^{-2}$ (Daddi et al.,
2010; Tacconi et al., 2013; Kennicutt & De Los Reyes, 2021, e.g.,, equivalent
to $\rm N_{\rm H}\approx 10^{22}-10^{24}\>cm^{-2}$), and high ISM column
densities, with typical values of $N_{\rm H}=9.8_{-0.7}^{+1.4}\times
10^{23}\rm\>cm^{-2}$ (e.g., Simpson et al., 2017). Extreme host galaxy
obscuration has also been found in the local ultraluminous IR galaxies NGC
4418 and Arp 220, which are estimated to be obscured by $N_{\rm H}\sim
10^{25}-10^{26}\rm\>cm^{-2}$ (e.g., Barcos-Muñoz et al., 2015; Dwek & Arendt,
2020; Sakamoto et al., 2021).

Large ISM column densities are also predicted in simulations. For example,
Trebitsch et al. (2019) found that during periods of rapid growth, the gas at
galactic scales has a column density at least comparable to the gas at nuclear
scales. Similarly, Bournaud et al. (2011) finds that disk instabilities can
trigger mass inflow towards the central engine and that during this phase, the
total column density is dominated by the gas over galactic scales, which can
achieve $\rm\langle N_{\rm H,ISM}\rangle>10^{23}\>cm^{-2}$, and even reach CT
levels. However, in these simulations, unobscured l.o.s. are still expected
towards the central engine 10% of the time, in agreement with the broad-line
(unobscured) quasars residing in starburst galaxies found in our work.

We note that the submm sizes of the submm-IR quasars hosts and SMGs are much
more compact than those of massive main-sequence galaxies at similar
redshifts, ranging between $\approx 3-6\rm\>kpc$ (e.g., Kaasinen et al.,
2020), which implies typical ISM column densities of $\langle N_{\rm
H,ISM}\rangle\approx(1-9)\cdot 10^{21}\>\rm cm^{-2}$ for comparably massive
main sequence galaxies with $\rm SFR\approx 50-90\rm\>M_{\odot}\>yr^{-1}$.
Assuming that the quasars living in normal star-forming galaxies have submm
sizes comparable to typical main-sequence galaxies, we would therefore expect
their ISM column densities to be typically unobscured and, therefore, to have
a negligible contribution to the quasar obscuration, particularly for $\rm
SFR\lesssim 100\>M_{\odot}\>yr^{-1}$, consistent with our results.

Overall, our results indicate that a large component of the quasar obscuration
can come from the ISM, in addition to the obscuration from the dusty torus.
Figure 5 illustrates three scenarios of quasar obscuration along different
lines of sight. At $\rm SFR\lesssim 300\>M_{\odot}\>yr^{-1}$, the obscured
quasar fraction remains constant around $\approx 0.65$; hence, it is
reasonable to assume that the torus dominates the obscuration for these
systems (panel a). Assuming no increase in the torus covering factor with SFR,
the increase in the obscured quasar fraction to $\approx 0.8-0.9$ would be due
to the ISM, implying that the ISM solely obscures $\approx 10-30\%$ of the
quasars with $\rm SFR\gtrsim 300\>M_{\odot}\>yr^{-1}$ (see panels b and c).
Taking a basic probabilistic approach and assuming that the orientation of the
torus and the ISM are independent (see panel b), the ISM dust-covering
fraction would be $\approx 0.4-0.7$. However, if there is some co-alignment
(see panel c), then the ISM dust-covering fraction could be as high as $0.9$
in the highest SFR systems.

Figure 5: Illustration of the dominant sources of obscuration. Panels (a),
(b), and (c) show different relative orientations between the torus (orange
flared disk) and host galaxy (blue flared disk with stars). The dashed and
continuous arrows represent obscured and unobscured l.o.s, respectively. The
cartoons are illustrative, and the sizes and opening angles are not drawn to
scale.

###  4.2 Implications of our results for quasar models

The similarity in the host galaxy properties of the SMGs and submm-IR quasars,
together with their large estimated ISM column densities, suggest that a
simple unification model postulating that the obscuration is only due to the
dusty torus does not hold for these extreme systems. It has been postulated
that SMGs (and consequently, submm-IR quasars) may be experiencing a
compaction phase where the ISM material is falling towards the centre of the
galaxy, triggering a compact starburst, and subsequently activating the SMBH,
which is supported by the fact that both populations have similarly compact
submm sizes. This phenomenon, also known as ’wet compaction’, is triggered by
an extreme loss of angular momentum caused by either violent disc
instabilities, counter-rotating streams or wet mergers (e.g., Dekel & Burkert,
2014; Zolotov et al., 2015; Lapiner et al., 2023). A number of studies have
proposed ’wet compaction’ as the explanation for the high fraction of compact
star-forming galaxies hosting an AGN, as compared to extended co-eval star-
forming galaxies with similar stellar masses and redshifts (e.g., Kocevski et
al., 2017; Chang et al., 2020; Silverman et al., 2019; Lamperti et al., 2021),
an explanation also supported by cosmological simulations (e.g., Habouzit et
al., 2019; Lapiner et al., 2023).

Recent results have shown that AGN in mergers are typically heavily obscured
(e.g., Ricci et al., 2017), that obscured AGN/quasars reside in preferential
close galaxy pairs and massive large scale environments (e.g., Dougherty et
al., 2023; Petter et al., 2023), and have also excess radio emission over that
found for unobscured quasars (A22), consistent with that expected from quasar-
driven outflows and jets interacting with the host-galaxy ISM (i.e., “AGN
feedback”, Klindt et al., 2019; Fawcett et al., 2020; Calistro Rivera et al.,
2021, e.g.,). Our study extends these differences to the star-forming
environments, showing that obscured quasars are preferentially found in
compact starbursts. Combined, these suite of results can be used to argue that
different quasar sub populations represent different phases in the overall
evolution of quasars (e.g., Lapi et al., 2014; Mancuso et al., 2016; Lapi et
al., 2018). A key uncertainty is the submm sizes of the lower SFR quasars
($\rm SFR\lesssim 100\>M_{\odot}\>yr^{-1}$), which can be used to determine
their relationship with the high SFR quasars studied here (e.g., comparably
extended to co-eval main-sequence galaxies versus compact like the high SFR
quasars). A combined ALMA and radio analysis of IR quasars across the full
range of SFR would provide key insights into the relationship between
different quasars sub populations and test quasar evolutionary models.

## Acknowledgements

This work has been supported by the EU H2020-MSCA-ITN-2019 Project 860744
“BiD4BESt: Big Data applications for black hole Evolution STudies”. DMA and AP
thank the Science Technology Facilities Council (STFC) for support from the
Durham consolidated grant (ST/T000244/1). AP acknowledges partial support by
STFC through grant ST/P000541/1. CH acknowledges funding from an United
Kingdom Research and Innovation grant (code: MR/V022830/1). AL is partly
supported by the PRIN MIUR 2017 prot. 20173ML3WW 002 ‘Opening the ALMA window
on the cosmic evolution of gas, stars, and massive black holes’. DJR
acknowledges the support of STFC grant NU-012097

## Data Availability

The datasets generated and/or analysed in this study are available from the
corresponding author on reasonable request.

## References

  * Aird et al. (2019) Aird J., Coil A. L., Georgakakis A., 2019, MNRAS, 484, 4360
  * Alexander & Hickox (2012) Alexander D. M., Hickox R. C., 2012, New Astron. Rev., 56, 93
  * Alexander et al. (2005) Alexander D. M., Bauer F. E., Chapman S. C., Smail I., Blain A. W., Brandt W. N., Ivison R. J., 2005, ApJ, 632, 736
  * Andonie et al. (2022) Andonie C., et al., 2022, MNRAS, 517, 2577
  * Antonucci (1993) Antonucci R., 1993, ARA&A, 31, 473
  * Barcos-Muñoz et al. (2015) Barcos-Muñoz L., et al., 2015, ApJ, 799, 10
  * Birkin et al. (2021) Birkin J. E., et al., 2021, MNRAS, 501, 3926
  * Bournaud et al. (2011) Bournaud F., Dekel A., Teyssier R., Cacciato M., Daddi E., Juneau S., Shankar F., 2011, ApJ, 741, L33
  * Buchner et al. (2014) Buchner J., et al., 2014, A&A, 564, A125
  * Buchner et al. (2017) Buchner J., Schulze S., Bauer F. E., 2017, MNRAS, 464, 4545
  * Buchner et al. (2019) Buchner J., Brightman M., Nandra K., Nikutta R., Bauer F. E., 2019, Astronomy & Astrophysics, 629, A16
  * Calistro Rivera et al. (2021) Calistro Rivera G., et al., 2021, A&A, 649, A102
  * CalistroRivera et al. (2018) CalistroRivera G., et al., 2018, ApJ, 863, 56
  * Carilli & Walter (2013) Carilli C. L., Walter F., 2013, ARA&A, 51, 105
  * Chabrier (2003) Chabrier G., 2003, ApJ, 586, L133
  * Chang et al. (2020) Chang Y.-Y., et al., 2020, ApJ, 888, 44
  * Chen et al. (2015) Chen C.-T. J., et al., 2015, ApJ, 802, 50
  * Circosta et al. (2019) Circosta C., et al., 2019, A&A, 623, A172
  * Civano et al. (2016) Civano F., et al., 2016, ApJ, 819, 62
  * Daddi et al. (2010) Daddi E., et al., 2010, ApJ, 714, L118
  * Dekel & Burkert (2014) Dekel A., Burkert A., 2014, MNRAS, 438, 1870
  * DiPompeo et al. (2014) DiPompeo M. A., Myers A. D., Hickox R. C., Geach J. E., Hainline K. N., 2014, MNRAS, 442, 3443
  * Donoso et al. (2014) Donoso E., Yan L., Stern D., Assef R. J., 2014, ApJ, 789, 44
  * Dougherty et al. (2023) Dougherty S. L., Harrison C. M., Kocevski D. D., Rosario D. J., 2023, MNRAS,
  * Dwek & Arendt (2020) Dwek E., Arendt R. G., 2020, ApJ, 901, 36
  * Fawcett et al. (2020) Fawcett V. A., Alexander D. M., Rosario D. J., Klindt L., Fotopoulou S., Lusso E., Morabito L. K., Calistro Rivera G., 2020, MNRAS, 494, 4802
  * Fu et al. (2012) Fu J., Kauffmann G., Li C., Guo Q., 2012, MNRAS, 424, 2701
  * Fujimoto et al. (2017) Fujimoto S., Ouchi M., Shibuya T., Nagai H., 2017, ApJ, 850, 83
  * Gilli et al. (2007) Gilli R., Comastri A., Hasinger G., 2007, A&A, 463, 79
  * Gilli et al. (2022) Gilli R., et al., 2022, arXiv e-prints, p. arXiv:2206.03508
  * Gullberg et al. (2019) Gullberg B., et al., 2019, MNRAS, 490, 4956
  * Habouzit et al. (2019) Habouzit M., et al., 2019, MNRAS, 484, 4413
  * Hickox & Alexander (2018) Hickox R. C., Alexander D. M., 2018, ARA&A, 56, 625
  * Hickox et al. (2011) Hickox R. C., et al., 2011, ApJ, 731, 117
  * Hinshaw et al. (2013) Hinshaw G., et al., 2013, ApJS, 208, 19
  * Hopkins et al. (2008) Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008, ApJS, 175, 356
  * Kaasinen et al. (2020) Kaasinen M., et al., 2020, ApJ, 899, 37
  * Kennicutt & De Los Reyes (2021) Kennicutt Robert C. J., De Los Reyes M. A. C., 2021, ApJ, 908, 61
  * Klindt et al. (2019) Klindt L., Alexander D. M., Rosario D. J., Lusso E., Fotopoulou S., 2019, MNRAS, 488, 3109
  * Kocevski et al. (2017) Kocevski D. D., et al., 2017, ApJ, 846, 112
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Lagos et al. (2011) Lagos C. D. P., Baugh C. M., Lacey C. G., Benson A. J., Kim H.-S., Power C., 2011, MNRAS, 418, 1649
  * Laloux et al. (2023) Laloux B., et al., 2023, MNRAS, 518, 2546
  * Lamperti et al. (2021) Lamperti I., et al., 2021, A&A, 654, A90
  * Lanzuisi et al. (2018) Lanzuisi G., et al., 2018, MNRAS, 480, 2578
  * Lapi et al. (2014) Lapi A., Raimundo S., Aversa R., Cai Z. Y., Negrello M., Celotti A., De Zotti G., Danese L., 2014, ApJ, 782, 69
  * Lapi et al. (2018) Lapi A., et al., 2018, ApJ, 857, 22
  * Lapiner et al. (2023) Lapiner S., et al., 2023, MNRAS, 522, 4515
  * Liu et al. (2019) Liu D., et al., 2019, ApJS, 244, 40
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Magdis et al. (2012) Magdis G. E., et al., 2012, ApJ, 760, 6
  * Mancuso et al. (2016) Mancuso C., Lapi A., Shi J., Cai Z. Y., Gonzalez-Nuevo J., Béthermin M., Danese L., 2016, ApJ, 833, 152
  * Marchesi et al. (2016) Marchesi S., et al., 2016, ApJ, 830, 100
  * Petter et al. (2023) Petter G. C., Hickox R. C., Alexander D. M., Myers A. D., Geach J. E., Whalen K. E., Andonie C. P., 2023, ApJ, 946, 27
  * Popping et al. (2022) Popping G., et al., 2022, MNRAS, 510, 3321
  * Powell et al. (2018) Powell M. C., et al., 2018, ApJ, 858, 110
  * Puglisi et al. (2019) Puglisi A., et al., 2019, ApJ, 877, L23
  * Puglisi et al. (2021) Puglisi A., et al., 2021, MNRAS, 508, 5217
  * Ramos Almeida & Ricci (2017) Ramos Almeida C., Ricci C., 2017, Nature Astronomy, 1, 679
  * Ranalli et al. (2003) Ranalli P., Comastri A., Setti G., 2003, A&A, 399, 39
  * Ricci et al. (2017) Ricci C., et al., 2017, MNRAS, 468, 1273
  * Rosario (2019) Rosario D. J., 2019, FortesFit: Flexible spectral energy distribution modelling with a Bayesian backbone (ascl:1904.011)
  * Rosario et al. (2013) Rosario D. J., et al., 2013, A&A, 560, A72
  * Sakamoto et al. (2021) Sakamoto K., González-Alfonso E., Martín S., Wilner D. J., Aalto S., Evans A. S., Harada N., 2021, ApJ, 923, 206
  * Sanders et al. (1988) Sanders D. B., Soifer B. T., Elias J. H., Madore B. F., Matthews K., Neugebauer G., Scoville N. Z., 1988, ApJ, 325, 74
  * Satyapal et al. (2014) Satyapal S., Ellison S. L., McAlpine W., Hickox R. C., Patton D. R., Mendel J. T., 2014, MNRAS, 441, 1297
  * Scoville et al. (2007) Scoville N., et al., 2007, ApJS, 172, 1
  * Scoville et al. (2016) Scoville N., et al., 2016, ApJ, 820, 83
  * Silverman et al. (2019) Silverman J. D., et al., 2019, ApJ, 887, L5
  * Silverman et al. (2023) Silverman J. D., et al., 2023, arXiv e-prints, p. arXiv:2306.03205
  * Simpson et al. (2017) Simpson J. M., et al., 2017, ApJ, 839, 58
  * Stern (2015) Stern D., 2015, ApJ, 807, 129
  * Tacconi et al. (2013) Tacconi L. J., et al., 2013, ApJ, 768, 74
  * Trebitsch et al. (2019) Trebitsch M., Volonteri M., Dubois Y., 2019, MNRAS, 487, 819
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, PASP, 107, 803
  * Villarroel & Korn (2014) Villarroel B., Korn A. J., 2014, Nature Physics, 10, 417
  * Weston et al. (2017) Weston M. E., McIntosh D. H., Brodwin M., Mann J., Cooper A., McConnell A., Nielsen J. L., 2017, MNRAS, 464, 3882
  * Whitaker et al. (2012) Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012, ApJ, 754, L29
  * Zolotov et al. (2015) Zolotov A., et al., 2015, MNRAS, 450, 2327


# Learning graph-Fourier spectra of textured surface images for defect
localization

ff

###### Abstract

In the realm of industrial manufacturing, product inspection remains a
significant bottleneck, with only a small fraction of manufactured items
undergoing inspection for surface defects. Advances in imaging systems and AI
can allow automated full inspection of manufactured surfaces. However, even
the most contemporary imaging and machine learning methods perform poorly for
detecting defects in images with highly textured backgrounds, that stem from
diverse manufacturing processes. This paper introduces an approach based on
graph Fourier analysis to automatically identify defective images, as well as
crucial graph Fourier coefficients that inform the defects in images amidst
highly textured backgrounds. The approach capitalizes on the ability of graph
representations to capture the complex dynamics inherent in high-dimensional
data, preserving crucial locality properties in a lower-dimensional space. A
convolutional neural network model (1D-CNN) was trained with the coefficients
of the graph Fourier transform of the images as the input to identify, with
classification accuracy of 99.4%, if the image contains a defect. An
explainable AI method using SHAP (SHapley Additive exPlanations) was used to
further analyze the trained 1D-CNN model to discern important spectral
coefficients for each image. This approach sheds light on the crucial
contribution of low-frequency graph eigen waveforms to precisely localize
surface defects in images, thereby advancing the realization of zero-defect
manufacturing.

Keywords: Quality assurance; Surface defect localization; Graph fourier
transforms; Convolutional neural networks; Explainbale AI

##  1 Introduction

In recent years, the field of quality control in manufacturing has witnessed a
significant transformation owing to advancements in imaging systems. A study
by GlobeNewswire, projects a remarkable surge for AI-based visual systems in
manufacturing, anticipating an estimated global market value of $21.3 billion
by 2028, with a compelling compound annual growth rate of 43.4% [3]. The
augmented availability and heightened sophistication of AI-based vision
systems have transcended their role beyond traditional end-of-line product
inspection to encompass on-line, intermittent and in-situ quality monitoring
[18, 22]. Imaging systems have evolved into indispensable tools for
identifying a myriad of morphological and geometric defects, including
scratches, holes, warping, and inclusions, throughout various stages of the
manufacturing process [21].

One key aspect that contributes to the complexity of defect identification is
the diverse range of textures generated by different manufacturing processes.
In metallic manufacturing, casting processes such as sand casting or
investment casting yield textures characterized by irregularities stemming
from the cooling and solidification of molten metal. In additive
manufacturing, textures emerge as a result of layer-wise deposition of
material, frequently accompanied by porosity and shrinkage defects [26, 32,
25]. Machining processes, including milling and turning, impart distinct
surface finishes with tool marks and cutting patterns, contributing to the
overall texture. Turning or cutting operations performed on a lathe often lead
to defects like burrs, grooves, adhered material, and tearings [23, 24].
Plastic and polymer manufacturing, specifically injection molding, introduces
textures contingent on mold surface finish and polymer type, with defects like
sink marks and weld lines influencing the visual appearance. In the realm of
composite manufacturing, lay-up processes shape textures through the
arrangement of fibers and resin distribution, while compression molding
imparts textures affected by molding surface characteristics. Semiconductor
manufacturing involves intricate processes yielding specific surface textures
influenced by factors such as etching and deposition techniques, and defects
in these processes can profoundly impact the chip’s functionality [27, 28]. In
contrast to open-source datasets such as MVTEC-AD [2] that primarily focus on
defect localization, machining processes generate intricate background
textures, as depicted in Figure 1. The complexity of these textures
underscores the significance of understanding their diversity, representing a
crucial aspect in advancing methodologies for defect identification within
manufacturing contexts.

Figure 1: Images illustrating distinct complex surface textures resulting from
a turning process (left) and ripple formation from laser direct energy
deposition (L-DED) process (right) [32]

From an imaging perspective, defects emerge as irregularities embedded within
the diverse textures, presenting a myriad of forms that pose challenges for
precise classification and localization. The intricacies reside on the subtle
interactions of pixel intensities within textures, and their accurate
identification holds key to unveiling the underlying defect formation
mechanisms across diverse manufacturing processes. Nonetheless, this precision
remains elusive with current image segmentation methods. Traditional
supervised learning methods [26, 29, 21, 28], while effective, often demand
substantial amounts of labeled data for each type of defect. Despite their
prowess in defect classification, these methods encounter limitations in
localizing defects accurately and inferring the underlying defect formation
mechanisms. Additionally, the need for a diverse dataset that spans the range
of possible defects across manufacturing processes poses a practical
challenge. Ganatma et al. [14] developed an AI Toolbox to relate distinct
class of defects with specific image processing methodologies such as Wavelet
Analysis, Morphological Component Analysis (MCA), and Basic Line Detector
(BLD). They proposed an empirical recommendation formula based on three image
metrics$-$entropy, Universal Quality Index (UQI), and Rosenberger’s$-$to
evaluate the performance of a method across a specified class of images.

As an alternative paradigm, graph signal processing (GSP) methods aim to
leverage the relationship between signals and the topology of the graph where
they are supported, to make meaningful inferences [8, 9, 10]. The intuition
behind using graph representations and analysis can be associated to the
approach of Belkin and Niyogi [4], where they highlighted that the graph
Laplacian can capture the complex dynamics of high-dimensional data and
locality preserving properties in a low-dimension space, namely, the graph
$G(V,E)$. Algebraic graph theoretic approaches have been adopted for process
monitoring and anomaly detection in ultraprecision machining using
multivariate time series sensor data [20, 5, 6]. Bukkapatnam et al. [7]
employed a planar random graph representation to monitor the surface
morphology evolution of electron beam printed Ti-6Al-4V samples during
polishing and established an endpoint criterion based on Fiedler number
$\lambda_{2}$, that is sensitive to the neighbourhood asperity structure.

Spectral graph theoretic approaches have been previously used to describe
topological relationships in various physical domains involving image
processing [19, 16, 15]. Hunag et al. [30] utilized GSP framework to analyze
brain activity aligned with the structural brain graph, revealing signal
variability during attention-switching tasks.They emphasized the framework’s
versatility for conducting analyses on both functional and structural
connectivity in functional imaging datasets. However, in most real world
scenarios, the underlying graph topology is not readily available, and needs
to be inferred from the signals itself. The construction of a meaningful graph
topology from the data is crucial for effective representation and inferring
functional connectivity within the graph [11, 12]. Models for developing graph
topology $G$ from data are either defined based on the underlying physical
phenomenon on the graph like heat diffusion, or statistically modelled as a
function that draws realization from a probability distribution over the
variables that are representative of the graph structure [13].

In this study, we harness the principles of Graph Fourier Transform (GFT) to
effectively localize the presence of surface defects in images. Utilizing GFT
introduces a powerful means of handling complex spatial relationships within
images with textured backgrounds. By conceptualizing the image as a signal
distributed across a graph, where pixels are interconnected based on their
spatial relationships, we aim to uncover specific graph spectral coefficients
crucial for distinguishing and localizing defects. The graph spectral
coefficients pave way to a more robust identification of nuanced morphological
characteristics, contributing to superior accuracy in localizing defects
amidst intricate image structures. A convolutional neural network (1D-CNN)
classifier is employed to categorize images with diverse surface textures,
discerning between those devoid of defects and those exhibiting defects with
varying shapes, sizes, and locations, encompassing a spectrum of morphological
characteristics. Inferences from subsequent explainable AI method using SHAP
(SHapley Additive exPlanations), helps identify important spectral
coefficients for defect localization within each image, and forms the basis
for establishing correlations with the underlying physical mechanisms
responsible for defect formation. This graph-based approach enables us to
capture both global patterns and localized variations, making it a powerful
tool for precise defect identification—an essential element in achieving zero-
defect manufacturing.

The rest of this paper is organized as follows: Section 2 provides a
comprehensive overview of graph signal processing and its application to
surface defect localization. In Section 3, we delineate the architecture of
the 1D-CNN model and evaluate its performance. In Section 4, we delve into the
interpretability of our model’s predictions through SHAP analysis to infer how
spectral coefficients contribute to defect localization, followed by
conclusions in Section 5.

##  2 Spectral graph theoretic approach for defect localization

We introduce the formulation of graph and graph signals, i.e., signals whose
samples are indexed by the nodes of an arbitrary graph. Let us define some of
the notations associated with graph signals, graph fourier transform and the
interpretation of graph variation operators.

Consider a weighted graph $G=(V,E,W)$ where $V=\\{1,2,\ldots,N\\}$ is a finite
set of $N$ nodes or vertices, $E\in V\times V$ is a set of edges defined as
ordered pairs $(i,j)$ representing pairwise relations between vertices, and
$W:E\rightarrow\mathbb{R}$ is a mapping from the set of edges to scalar values
$w_{ij}$ representing the edge weight between vertex $i$ to vertex $j$. A
graph signal $X:V\rightarrow\mathbb{C}$ is a mapping from the vertices of the
graph to complex numbers $X\in\mathbb{C}^{N}$. This graph signal can be
represented as a vector $X=\\{x_{1},x_{2},x_{3},\ldots,x_{N}\\}$ where $x_{k}$
indicates the value of the signal at the $k^{th}$ vertex in $V$.

Figure 2: Schematic representation of Graph Fourier Transforms for spectral
coefficients extraction and defect localization

Here, we define a grid-based undirected graph structure on an image (see
Figure 2) by representing each pixel as a node and connecting adjacent pixels
with edges $(E)$ defined as ordered pairs $(i,j)$ where $(i)$ and $(j)$ are
pixel indices. In contrast to traditional image segmentation methods, where
edge weights are determined based on pixel intensities utilizing exponential
or gaussian functions, our approach formulates a graph structure solely based
on spatial relationships, specifically considering immediate neighbors in the
$x$ and $y$ directions. The associated adjacency matrix
$A\in\mathbb{R}^{N\times N}$ for the graph is defined as

$A_{ij}=\begin{cases}w_{ij}=1,&\text{if }(i,j)\in E\\\
0,&\text{otherwise}\end{cases}$ (1)

The degree matrix $D\in\mathbb{R}^{N\times N}$ for the graph is a diagonal
matrix defined as follows:

$D_{ii}=\text{deg}(i)=\sum_{j\in n(i)}w_{ij}$ (2)

where $n(i)$ stands for the neighborhood of vertex $i$. The laplacian matrix
$L\in\mathbb{R}^{N\times N}$ for the above defined graph structure $G$ with
adjacency matrix $A$ and degree matrix $D$ is defined as

$L=D-A$ (3)

For undirected graphs with real and non-negative edge weights, the graph
Laplacian matrix $L$ inherits properties of being real, symmetric, and
positive semi-definite [17].

###  2.1 Eigen decomposition of the graph laplacian matrix

The real and symmetric nature of the graph Laplacian matrix $L$ allows for
eigen decomposition into a complete set of orthonormal eigenvectors
$u_{0},u_{1},\ldots,u_{N-1}$ and corresponding eigenvalues, in ascending order
$0=\lambda_{0}<\lambda_{1}\leq\ldots\leq\lambda_{N-1}$:

$L=U\Lambda U^{T}$ (4)

where $U$ is the eigenvector matrix containing the eigenvectors
$u_{0},u_{1},\ldots,u_{N-1}$ as columns, satisfying $U^{-1}=U^{T}$, and
$\Lambda$ is a diagonal matrix containing the eigenvalues, i.e.,
$\Lambda=\text{diag}(\lambda_{0},\lambda_{1},\ldots,\lambda_{N-1})$. Each
eigenvector $u_{k}$ of the graph laplacian matrix $L$ satisfies
$Lu_{k}=\lambda_{k}u_{k}$. The eigenvalue $\lambda_{k}$ is referred to as the
graph frequency/spectrum, with a smaller eigenvalue corresponding to a lower
graph frequency [17].

The eigen decomposition of the graph Laplacian matrix $L$ into eigenvectors
$U$ and corresponding eigenvalues $\Lambda$ embodies a profound depiction of
the graph’s underlying structure and dynamics. In the context of a grid
structure representing an image, where each edge directly connects neighboring
pixels and bears a unit weight, the derived eigenvectors $U$ offer nuanced
insights into the image’s representation. The eigenvector linked to
$\lambda_{0}$ acts as a fundamental mode, resembling a constant function and
offering a comprehensive perspective on the overall structure of the image
grid. Similarly, the eigenvector associated with $\lambda_{1}$ (second
smallest eigen value), widely recognized as the Fiedler vector, encapsulates
vital information regarding the graph’s connectivity, acting as a
representation of the first non-trivial mode.

Figure 3: Eigen waveforms illustrating select low and high frequency modes of
the defined graph structure

Moving beyond these foundational eigenvectors, the subsequent high-frequency
modes unveil localized variations within the image grid. These high-frequency
eigenvectors delve into fine-grained features, highlighting intricate patterns
and details that might be imperceptible in lower-frequency modes. As
eigenvalues ascend, each associated eigenvector offers a more refined
perspective on the spatial relationships and structural intricacies within the
image. Figure 3 provides a visual representation of few select eigen vectors
having low and high frequency waveforms.

###  2.2 Graph fourier transform

The graph Fourier transform (GFT) is a mathematical tool applied to signals
defined on graphs. For a given graph Laplacian operator $L=U\Lambda U^{T}$,
where $U$ is the matrix of eigenvectors and $\Lambda$ is the diagonal matrix
of corresponding eigenvalues, the GFT of a signal
$X=(x_{1},x_{2},\ldots,x_{N})$ is defined as

$S=U^{T}X,\quad s_{k}=\langle X,u_{k}\rangle=\sum_{i=1}^{N}x_{i}\cdot
u_{k}(i-1)$ (5)

where spectral coefficient $s_{k}=\langle X,u_{k}\rangle$ is the inner product
between the signal $X$ and eigenvector $u_{k}$, representing the projection of
the signal onto the $k$th graph frequency mode.

Anomalous patterns or defects within an image often manifest as distinct
localized variations, stemming from diverse physical mechanisms encountered
throughout various stages of the manufacturing process. The spectral
coefficients $S=\\{s_{0},s_{1},s_{3},\ldots,s_{N-1}\\}$ derived by projecting
the signal onto the eigenvectors not only accentuate the presence of defects
in the graph frequency domain but also provide a distinctive signature for
precise defect localization.

Leveraging the orthogonality property of the eigenvector matrix ($U\cdot
U^{T}=I$), the original image can be precisely reconstructed using $X=US$. In
sections 3 and 4, we employ a 1D-CNN binary classifier using spectral
coefficients from an image dataset characterized by diverse background
textures, including images with defects and non-defect scenarios. Employing
SHAP explanations, we identify important spectral coefficients $S^{*}$ for
each image. The localization of defects becomes visually discernible through
the reconstruction process $X^{*}=US^{*}$, capitalizing on the essential
spectral coefficients that encapsulate distinct signatures of the defects.

##  3 1D Convolutional neural network classifier

We construct a binary classifier using 1D-CNN, taking the spectral
coefficients as inputs from an image dataset characterized by diverse
background textures, including images with defects, and non-defect scenarios.
We utilized a dataset of 256 $\times$ 256 surface images obtained from the
turning process, generously provided by our collaborator Marposs-STIL. This
dataset comprises 228 instances of defects and 800 non-defective images. To
prepare for the subsequent classification task, we extracted 1000 defect and
non-defect patches, each measuring 64 $\times$ 64, from the original dataset,
as illustrated in Figure 4. Employing set rotation and flipping for
augmentation, we expanded the dataset to include 6000 defect and 6000 non-
defect patches. Following this, Graph Fourier Transform (GFT) features were
extracted from each patch, resulting in a final dataset sized 12000 $\times$
4096\. To assess the effectiveness of our methodology, an 80-20 split was
employed, allocating 80% of the dataset for training a 1D Convolutional Neural
Network (CNN) model on the GFT features. This model was purposefully designed
to classify patches into defect and non-defect categories, leveraging GFT
features as the basis for our exploration into surface image defect
localization.

Figure 4: Top row containing images with diverse background textures and no
defects, and bottom row showcasing images with defects.

The proposed model architecture for defect classification, illustrated in
Figure 5, involves a sequence of two 1D Convolution layers with ReLU
activation functions, utilizing 32 and 64 kernels and a stride of 3. Following
the convolutional layers, the output is flattened and subjected to a dropout
layer. Subsequently, two fully connected layers are integrated, featuring
hidden neurons arranged as (128, 2), accompanied by ReLU and Sigmoid
activation functions, respectively. The model is compiled using categorical
cross-entropy as the chosen loss function, the Adam optimizer (learning rate =
0.001), and accuracy as the primary evaluation metric. The training regimen
spans 50 epochs with a batch size of 64, utilizing labeled image vectors for
both training and testing phases. This architectural configuration, along with
the specified training parameters, is strategically designed to enhance the
model’s proficiency in accurately classifying defects and non-defects within
the surface image dataset. Figure 6 shows that the model achieves a training
accuracy of 99.94% and a testing accuracy of 99.91%, underscoring its robust
performance in accurately classifying defects and non-defects within the
surface image dataset.

Figure 5: 1D-CNN model architecture Figure 6: Performance of 1D-CNN binary
classifier

To validate the significance of spectral coefficients in classification, we
employed a Random Forest Classifier with these coefficients as features,
achieving notable training and testing accuracies of 99.89% and 96.87%,
respectively. This underscores the discriminative power of spectral
coefficients in defect identification. Subsequently, our focus shifts to SHAP
explanations from the trained 1D-CNN classifier in the upcoming section. This
analysis aims to pinpoint essential spectral coefficients crucial for defect
localization, providing key insights into the model’s decision-making process.

##  4 Results and discussion

After successfully classifying an image as defect or non-defect using the
trained 1D CNN model $f$, we aim to identify the GFT spectral coefficients
that exert the most influence on the model’s output. These dominant
coefficients provide insights into the defect localization process, enabling
us to move beyond binary classification and pinpoint the specific spectral
coefficients associated with defects.

Shapley value is a game theory-based metric widely used for quantifying the
contribution of each player in a cooperative game. In our context, Shapley
values measure the importance of each spectral coefficient in determining
whether an image contains a defect. Specifically, the contribution of a
coefficient $s_{i}$ to the model output for an input $X$ is determined by the
difference between outputs when including and excluding coefficient $s_{i}$,
expressed as $f({X}_{R\cup\\{s_{i}\\}})-f({X}_{R})$. Here, $R$ is a subset of
$S\setminus\\{s_{i}\\}$ (the entire set of coefficients $S$ excluding the
coefficient $s_{i}$). ${X}_{R}$ represents a modified input that only contains
coefficients in $R$ and takes the same coefficient values as ${X}$ does.
Overall, the marginal contribution (or Shapley value) of coefficient $s_{i}$
regarding the input $X$ (denoted as $\phi_{s_{i}}^{X}$) is computed as a
weighted average of the output differences over all possible coefficient
subsets, that is

$\phi_{s_{i}}^{X}=\sum_{R\subseteq
S\setminus\\{s_{i}\\}}\frac{|R|!(|S|-|R|-1)!}{|S|!}[f({X}_{R\cup\\{s_{i}\\}})-f({X}_{R})]$
(6)

wherein $|\cdot|$ is the size of a set. $\phi_{s_{i}}^{\bf x}$ satisfies

$\phi_{0}+\sum_{i=1}^{|S|}\phi_{s_{i}}^{X}=f({X})$ (7)

where $f({X})$ is the original model output on ${X}$. $\phi_{0}$ denotes the
model output when all spectral coefficients are zero, i.e., $\phi_{0}$ is a
constant based on the model’s historical training data regardless of ${X}$.
Based on Eq. (7), a positive Shapley value for a particular coefficient of a
defect image indicates an increased likelihood of the image being classified
as defective, and vice versa.

However, the computation complexity of Eq.(6) can be extrememly large when
considering all possible coefficient subsets [31]. Therefore, SHAP,
particularly the DeepSHAP algorithm in our case, approximates the Shapley
value using selected representative subsets. DeepSHAP also approximates
$f({X}_{R\cup\\{s_{i}\\}})-f({X}_{R})$ for a given coefficient subset $R$ by
leveraging the gradient information obtained through backpropagation [31]. The
resulting estimate of the Shapley value for coefficient $s_{i}$ is referred to
as the SHAP value, denoted by $\hat{\phi}_{s_{i}}^{X}$. When applied to
individual images, SHAP unveils important spectral coefficients to each image,
contributing to effective defect localization. In addition, SHAP also provides
a global inference perspective, which considers the collective impact of
coefficients across the entire dataset. Specifically, after obtaining the SHAP
value of coefficient $s_{i}$ for every image ${X}^{(j)}$, $j=1,2,...,M$, the
algorithm aggregate the SHAP values over all $M$ images to derive the global
contribution of coefficient $s_{i}$, denoted by $G_{s_{i}}$, as

$G_{s_{i}}=\sum_{j=1}^{M}|\hat{\phi}_{s_{i}}^{{X}^{(j)}}|,\qquad s_{i}\in S$
(8)

Two curated data subsets, each comprising 15 images, encompass instances of
defects and non-defects. SHAP values for each of the 15 images are
subsequently derived using SHAP DeepExplainer, trained on 500 images from the
overall training dataset, which includes both classes. The outcomes
demonstrate a pronounced discriminative capability between images with defects
and those without, particularly emphasizing the importance of low-frequency
spectral coefficients. Cumulative summaries of the top 10 important spectral
coefficients for classifying each data class are visually presented in Figures
7 and 8, offering insights into the key features influencing the
classification process.

Figure 7: Summary plot of SHAP values for a subset of correctly classified 15
defect images Figure 8: Summary plot of SHAP values for a subset of correctly
classified 15 non-defect images

Figure 9: Marginal distribution of low-frequency spectral coefficients for
defct and non-defect images

Figure 10: Root mean square measures of the first 500 low-frequency spectral
coefficients for defect and non-defect images

Subsequent to this analysis, we examine the marginal distribution of the first
16 low-frequency spectral coefficients, denoted as ${s_{0}}$ to ${s_{15}}$,
across the entire dataset for each class. As emphasized by SHAP, a noteworthy
shift in median values is observed, particularly for ${s_{4}}$, ${s_{5}}$, and
${s_{2}}$. Additionally, Fig 9 illustrates the distinct dissimilarity in the
marginal distribution of spectral coefficient ${s_{0}}$ between the two
classes. The defects examined in this study demonstrate a tendency to disturb
the low-frequency eigen waveforms in contrast to non-defect images with
distinctive background textures as observed in Figure 9. The pronounced
discriminative capability between defect and non-defect images, as clarified
by SHAP through low-frequency spectral coefficients, receives additional
validation through the Root Mean Square plot of the spectral coefficients
across the entire test dataset (Fig 10).

Figure 11: Defect localization outputs: Defect Input Image, Ground truth,
localization (Our method), and localization (SimpleNet)

Figure 12: Defect localization performance: Our Model vs SimpleNet

Defect localization on the 256 $\times$ 256 surface image is conducted by
classifying each 64x64 patch with a stride of 8 across the entire image,
leveraging the capabilities of the trained 1D CNN model. The outcome is a
sparse defect localization output of dimensions 49x49. The defect localization
results from our model are visually presented in the third column of Figure
11, providing a representation of the spatial localization of identified
defects across the surface image.

As discussed in Section 1 and illustrated in Figure 1, the surfaces derived
from various manufacturing processes exhibit complex background patterns,
complicating the development and training of algorithms tailored for defect
localization unique to the process. The proposed method revolves around the
core idea of capturing these complexities through the utilization of well-
established GFT features. Our defect localization approach is then assessed
through a comparative analysis with SimpleNet [1], a widely recognized state-
of-the-art methodology which stands out as a simple yet application-friendly
network. This comparative study highlights the challenge of formulating
effective and customized methods for defect localization on intricate
machining textures. On the MVTec AD benchmark, SimpleNet [1] demonstrates
notable performance, achieving an anomaly detection AUROC of 99.6% and
reducing errors by 55.5% compared to the next best-performing model. The
SimpleNet is trained using hyperparameters consistent with those employed by
its original authors for their study on the MVTEC-AD dataset. Figure 11
visually compares the defect localization outputs from both models, juxtaposed
with the ground truth. Further quantitative insights are presented in Figure
12, featuring boxplots comparing the models’ performance using AUROC, AUPR,
and F1 score metrics. Our observations reveal that our model excels in the
sparse localization and detection of defects within surface images.

##  5 Conclusion and Future Work

This work is one of the first investigations of graph Fourier analysis for
detection and localization of defects and such artifacts from images with
highly textured background. The approach harnesses the power of graph
representations to encapsulate intricate dynamics inherent in high-dimensional
data, specifically image pixels, while simultaneously preserving vital
locality properties within a lower-dimensional space, here, the graph. Through
the application of Graph Fourier Transforms, the extraction of essential
spectral coefficients from the image facilitates localization of defects and
artifacts within images characterized by highly textured backgrounds.

A convolutional neural network model (1D-CNN) trained on the graph Fourier
spectral coefficients was able to classify, with classification accuracy
exceeding 99%, whether the image contains a defect. The CNN classifier was
able to exceed the performance of perhaps the most contemporary image-based
defect (artifact) detection method (called SimpleNet) by 70% in terms of the
area under the receiver operating characteristic (AUROC), and F1 score. The
subsequent application of SHAP, an explainable AI (XAI) method enabled the
identification of the key components of the graph Fourier spectrum that inform
the defects, including its location and other morphological characteristics.
The SHAP analysis suggests that the eigen waveforms 4 and 5 (that capture the
lower frequency components) are most sensitive to the presence of the defects.

The current research assumes a constant graph structure to derive spectral
coefficients for each patch. A promising extension involves formulating a
dataset-specific graph structure through the utilization of Graph Neural
Networks (GNNs). By employing GNNs, we can dynamically construct a unique
graph topology tailored to the dataset characteristics. This innovative
approach would enable the derivation of spectral coefficients for patches,
potentially enhancing the efficiency of defect localization in surface images.
This extension holds the potential to further refine the model’s adaptability
to diverse datasets, contributing to more effective and nuanced defect
identification methodologies in manufacturing processes.

## Acknowledgements

This work was supported by Rockwell International Professorship and U.S.
National Science Foundation under Grant NSF-ECCS-1953694. The authors would
also like to express their sincere gratitude to MARPOSS-STIL for generously
providing all the image data.

## References

  * [1] Zhikang Liu, Yiming Zhou, Yuansheng Xu, Zilei Wang, _Simplenet: A simple network for image anomaly detection and localization_ , Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 20402–20411, 2023\.
  * [2] Paul Bergmann, Michael Fauser, David Sattlegger, Carsten Steger, _MVTec AD–A comprehensive real-world dataset for unsupervised anomaly detection_ , Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 9592–9600, 2019\.
  * [3] GlobeNewswire, _Global Artificial Intelligence in Manufacturing Market Size, Share and Industry Trends Analysis Report, Outlook and Forecast, 2022-2028_ , Report Linker, 2022\.
  * [4] Mikhail Belkin, Partha Niyogi, _Laplacian eigenmaps for dimensionality reduction and data representation_ , Neural computation, 15(6), pp. 1373–1396, 2003, MIT Press
  * [5] Mohammad Samie Tootooni, Prahalad K Rao, Chun-An Chou, Zhenyu James Kong, _A spectral graph theoretic approach for monitoring multivariate time series data from complex dynamical processes_ , IEEE Transactions on Automation Science and Engineering, 15(1), pp. 127–144, 2016, IEEE
  * [6] Prahalad K Rao, Satish Bukkapatnam, Zhenyu Kong, Omer Beyca, Kenneth Case, Ranga Komanduri, _Quantification of ultraprecision surface morphology using an algebraic graph theoretic approach_ , Procedia Manufacturing, 1, pp. 12–26, 2015, Elsevier
  * [7] Satish TS Bukkapatnam, Ashif Sikandar Iquebal, Soundar RT Kumara, _Planar random graph representations of spatiotemporal surface morphology: Application to finishing of 3-D printed components_ , CIRP Annals, 67(1), pp. 495–498, 2018, Elsevier
  * [8] Xiaowen Dong, Dorina Thanou, Laura Toni, Michael Bronstein, Pascal Frossard, _Graph signal processing for machine learning: A review and new perspectives_ , IEEE Signal Processing Magazine, 37(6), pp. 117–127, 2020, IEEE
  * [9] Antonio Ortega, Pascal Frossard, Jelena Kovačević, José MF Moura, Pierre Vandergheynst, _Graph signal processing: Overview, challenges, and applications_ , Proceedings of the IEEE, 106(5), pp. 808–828, 2018, IEEE
  * [10] David I Shuman, Sunil K Narang, Pascal Frossard, Antonio Ortega, Pierre Vandergheynst, _The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains_ , IEEE Signal Processing Magazine, 30(3), pp. 83–98, 2013, IEEE
  * [11] Xiaowen Dong, Dorina Thanou, Michael Rabbat, Pascal Frossard, _Learning graphs from data: A signal representation perspective_ , IEEE Signal Processing Magazine, 36(3), pp. 44–63, 2019, IEEE
  * [12] Gonzalo Mateos, Santiago Segarra, Antonio G Marques, Alejandro Ribeiro, _Connecting the dots: Identifying network structure via graph signal processing_ , IEEE Signal Processing Magazine, 36(3), pp. 16–43, 2019, IEEE
  * [13] Daphne Koller, Nir Friedman, _Probabilistic graphical models: principles and techniques_ , MIT Press, 2009
  * [14] Tapan Ganatma Nakkina, Man Vinayaka, Amr Masad, Mohamed El Mansori, Satish Bukkapatnam, _Smart defect identification for manufacturing applications_ , Surface Topography: Metrology and Properties, 10(3), 2022, IOP Publishing
  * [15] Jianbo Shi, Jitendra Malik, _Normalized cuts and image segmentation_ , IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(8), pp. 888–905, 2000, IEEE
  * [16] Ashif Sikandar Iquebal, Satish Bukkapatnam, _Consistent estimation of the max-flow problem: Towards unsupervised image segmentation_ , IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(5), pp. 2346–2357, 2020, IEEE
  * [17] Wei Hu, Jiahao Pang, Xianming Liu, Dong Tian, Chia-Wen Lin, Anthony Vetro, _Graph signal processing for geometric data and beyond: Theory and applications_ , IEEE Transactions on Multimedia, 24, pp. 3961–3977, 2021, IEEE
  * [18] Jorge F Arinez, Qing Chang, Robert X Gao, Chengying Xu, Jianjing Zhang, _Artificial intelligence in advanced manufacturing: Current status and future outlook_ , Journal of Manufacturing Science and Engineering, 142(11), 2020, American Society of Mechanical Engineers
  * [19] Gene Cheung, Enrico Magli, Yuichi Tanaka, Michael K Ng, _Graph spectral image processing_ , Proceedings of the IEEE, 106(5), pp. 907–930, 2018, IEEE
  * [20] Mohammad Montazeri, Prahalada Rao, _Sensor-based build condition monitoring in laser powder bed fusion additive manufacturing process using a spectral graph theoretic approach_ , Journal of Manufacturing Science and Engineering, 140(9), 2018, American Society of Mechanical Engineers
  * [21] Je-Kang Park, Bae-Keun Kwon, Jun-Hyub Park, Dong-Joong Kang, _Machine learning-based imaging system for surface defect inspection_ , International Journal of Precision Engineering and Manufacturing-Green Technology, 3, pp. 303–310, 2016, Springer
  * [22] Fotios K Konstantindis, Antonios Gasteratos, Spyridon G Mouroutsos, _Vision-based product tracking method for cyber-physical production systems in industry 4.0_ , 2018 IEEE international conference on imaging systems and techniques (IST), pp. 1–6, 2018, IEEE
  * [23] A Bordin, S Bruschi, A Ghiotti, _The effect of cutting speed and feed rate on the surface integrity in dry turning of CoCrMo alloy_ , Procedia Cirp, 13, pp. 219–224, 2014, Elsevier
  * [24] Durul Ulutan, Tugrul Ozel, _Machining induced surface integrity in titanium and nickel alloys: A review_ , International Journal of Machine Tools and Manufacture, 51(3), pp. 250–280, 2011, Elsevier
  * [25] Andrew Triantaphyllou, Claudiu L Giusca, Gavin D Macaulay, Felix Roerig, Matthias Hoebel, Richard K Leach, Ben Tomita, Katherine A Milne, _Surface texture measurement for additive manufacturing_ , Surface topography: metrology and properties, 3(2), 2015, IOP Publishing
  * [26] V Akhil, G Raghav, N Arunachalam, DS Srinivas, _Image data-based surface texture characterization and prediction using machine learning approaches for additive manufacturing_ , Journal of Computing and Information Science in Engineering, 20(2), 2020, American Society of Mechanical Engineers
  * [27] Takeshi Nakazawa, Deepak V Kulkarni, _Anomaly detection and segmentation for wafer defect patterns using deep convolutional encoder–decoder neural network architectures in semiconductor manufacturing_ , IEEE Transactions on Semiconductor Manufacturing, 32(2), pp. 250–256, 2019, IEEE
  * [28] Muhammad Saqlain, Qasim Abbas, Jong Yun Lee, _A deep convolutional neural network for wafer defect identification on an imbalanced dataset in semiconductor manufacturing processes_ , IEEE Transactions on Semiconductor Manufacturing, 33(3), pp. 436–444, 2020, IEEE
  * [29] Alessandra Caggiano, Jianjing Zhang, Vittorio Alfieri, Fabrizia Caiazzo, Robert Gao, Roberto Teti, _Machine learning-based image processing for on-line defect recognition in additive manufacturing_ , CIRP annals, 68(1), pp. 451–454, 2019, Elsevier
  * [30] Weiyu Huang, Thomas AW Bolton, John D Medaglia, Danielle S Bassett, Alejandro Ribeiro, Dimitri Van De Ville, _A graph signal processing perspective on functional brain imaging_ , Proceedings of the IEEE, 106(5), pp. 868–885, 2018, IEEE
  * [31] Scott M Lundberg, Su-In Lee, _A unified approach to interpreting model predictions_ , Advances in Neural Information Processing Systems, 30, 2017
  * [32] Himanshu Balhara, Bhaskar Botcha, Sarah J Wolff, Satish TS Bukkapatnam, _Ripple formations determine the heterogeneous microstructure of directed energy deposition (DED)-printed 316L components_ , Materials & Design, 227, 2023, Elsevier