# Peeling Back the Layers of Extinction of Dusty Galaxies in the Era of JWST:
Modelling Joint NIRSpec + MIRI Spectra at rest-frame 1.5 - 28 $\mu$m

F. R. Donnan, 1 I. García-Bernete,1 D. Rigopoulou,1,2 M. Pereira-Santaella,3
P. F. Roche,1 A. Alonso-Herrero4

1Department of Physics, University of Oxford, Keble Road, Oxford, OX1 3RH, UK
2School of Sciences, European University Cyprus, Diogenes street, Engomi, 1516
Nicosia, Cyprus
3 Instituto de Física Fundamental, CSIC, Calle Serrano 123, E-28006, Madrid,
Spain
4 Centro de Astrobiología (CAB), CSIC-INTA, Camino Bajo del Castillo s/n,
E-28692 Villanueva de la Cañada, Madrid, Spain E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present an analysis of the combined NIRSpec and MIRI spectra of dusty
galaxies between 1.5 - 28 $\mu$m restframe by implementing a differential
extinction model, where the strength of extinction varies across the spectrum
as different layers of the obscuring dust are probed. Our model is able to
recover a 2D distribution of dust temperature and extinction which allows
inference of the physical nature of the dust in these environments. We show
that differential extinction is necessary to reproduce the spectra of 4 highly
obscured Luminous Infrared Galaxies observed with NIRSpec IFU and MIRI MRS,
where simple screen or uniformly mixed dust distributions fail to fit the
data. We additionally compare the extinction of HII regions in these galaxies
via hydrogen recombination lines, the extinction of molecular gas via the H2
lines, Polycyclic Aromatic Hydrocarbons via the 12.7/11.3 PAH ratio and the
stellar continuum. We find that the molecular gas is deeply buried with the
HII regions in star-forming regions, with a similar extinction to the hottest
dust components. However we find the cooler dust to be less obscured, at a
similar extinction to the stellar continuum and PAHs. The nuclei show a
complex dust distribution with VV114 NE, NGC 3256 S, IIZw96 SW showing a
deeply buried continuum source relative to the molecular gas/HII regions.
Additionally, NGC 3256 S, NGC 7469 and VV114 SW show an isolated hot dust
component, indicative of AGN heating, where NGC 3256 S and NGC 7469 are
previously known AGN.

###### keywords:

galaxies: nuclei – galaxies: evolution – techniques: spectroscopic

††pubyear: 2015††pagerange: Peeling Back the Layers of Extinction of Dusty
Galaxies in the Era of JWST: Modelling Joint NIRSpec + MIRI Spectra at rest-
frame 1.5 - 28 $\mu$m–D

##  1 Introduction

Luminous Infrared Galaxies (LIRGs, LIR $>10^{11}\rm{L}_{\odot}$) and
Ultraluminous Infrared Galaxies (ULIRGs, LIR $>10^{12}\rm{L}_{\odot}$) are
dust and gas rich galaxies, with dust enshrouded star formation and AGN
activity (e.g. Sanders & Mirabel, 1996; Rigopoulou et al., 1999; Lonsdale et
al., 2006). These galaxies play a key role in the evolution of galaxies,
making up the bulk of the star formation rate density at cosmic noon (e.g. Le
Floc’h et al., 2005; Magnelli et al., 2011). In addition $\sim 40\%$ of local
ULIRGs contain Compact Obscured Nuclei (CONs; Falstad et al., 2021; García-
Bernete et al., 2022b; Donnan et al., 2023b), which hide rapid, active
supermassive black hole growth and/or extremely compact starbursts (e.g.
Veilleux et al., 2009; Aalto et al., 2015, 2019; Pereira-Santaella et al.,
2021).

As a consequence of the vast quantities of dust, emission in the optical is
absorbed and reprocessed in the infrared. Therefore observations in the
infrared are required to accurately measure properties of these galaxies such
as the star-formation rate, properties of the ISM etc. In particular, the mid-
infrared is a feature rich region of the spectrum, with dust emission and
absorption signatures, Polycyclic Aromatic Hydrocarbons (PAHs), and numerous
emission lines from ionised and molecular gas.

With all these emission/absorption features, the mid-IR spectra of galaxies
can be extremely complex, therefore extracting the emission features and
determining accurately the level of the continuum emission is a challenge. As
the PAH features are broad and blended they create a pseudo-continuum, under
which the true continuum lies (e.g. Brandl et al., 2006; Smith et al., 2007;
Gallimore et al., 2010; Li, 2020; Rigopoulou et al., 2021). This means a
simple, local continuum around emission features is often not sufficient to
obtain accurate PAH fluxes, particularly around the 7.7 $\mu$m band. This
problem is compounded in more obscured sources, where strong absorption
features from water ice ($\sim 6\mu$m), CH ($\sim 7\mu$m) and silicates ($\sim
9.8\mu$m) create a complex continuum (e.g. Spoon et al., 2007, 2022; Veilleux
et al., 2009; García-Bernete et al., 2024a).

In highly obscured environments it becomes difficult to reproduce the
continuum, particularly over a large wavelength range. This may be a
consequence of differential extinction, where multiple regions of different
dust densities all contribute within a single aperture. The resulting spectrum
is therefore a combination of different intrinsic spectra with varying levels
of obscuration. Such environments have been mapped in the Milky Way (e.g.
Lallement et al., 2022). Differential extinction has an impact on the shape of
the observed continuum, such as changing the ratio of the 18 $\mu$m to
9.8$\mu$m silicate features. This has been observed in AGN tori with clumpy
torus models showing a higher 18 $\mu$m to 9.8 $\mu$m ratio than smooth torus
models (e.g. Hatziminaoglou et al., 2015; Martínez-Paredes et al., 2020;
García-Bernete et al., 2024a).

Different approaches have been used to model the mid-IR spectra of galaxies,
however the majority of these assume a single level of extinction i.e an
extinction law is scaled by a single extinction value over the whole spectrum.
Popular choices include PAHFIT (Smith et al., 2007), which models the full
spectrum with a series of blackbodies subject to extinction and a series of
PAH emission features, with flexible profiles. This works well for relatively
unobscured star-forming galaxies which are dominated by PAH emission but can
struggle to fit AGN and/or obscured galaxies (e.g. Gallimore et al., 2010).
Another code, CAFE (Marshall et al., 2007), similarly uses modified
blackbodies but is restricted to only three where each has a different
extinction.

If an AGN is present, one approach is to model the dust from AGN using a
template, typically produced from radiative transfer simulations of the torus
(e.g. Ramos Almeida et al., 2009; Alonso-Herrero et al., 2011, 2012; Martínez-
Paredes et al., 2015, 2020, 2021; González-Martín et al., 2019a, b; García-
Bernete et al., 2019, 2022d). This approach is valuable in providing
constrains on physical parameters of the torus such as opening angle,
inclination etc. However as this is often highly model dependent, one often
needs to fit multiple libraries to obtain a complete picture of the torus for
a given object. A wide variety of libraries exist in the literature, with
different geometries, dust distributions etc. These include smooth torus
models (e.g. Fritz et al., 2006; Efstathiou et al., 2013), clumpy models (e.g.
Nenkova et al., 2008a, b; Hönig et al., 2010), two-phase models (e.g.
Stalevski et al., 2016) and models incorporating dust outflows (e.g. Hönig &
Kishimoto, 2017). These templates can also struggle to fit the data in certain
cases such as for the most obscured nuclei, where the deep silicate feature is
difficult to reproduce (García-Bernete et al., 2022b), with only the most
obscured smooth tori able to generate the deep silicates. As a result, semi-
empirical templates of highly obscured galaxies have often been used to
reproduce their spectra such as that employed in QUESTFIT (Veilleux et al.,
2009), however these often fail to fit the data due to a lack of model
flexibility.

While the techniques outlined above are extremely valuable, there is a need
for a more general, data driven approach that can fit a wide range of objects
such as star-forming regions, type 1/2 AGN, extremely obscured nuclei, etc,
but still provide inference of interesting/useful physical parameters of the
dusty structure in these galaxies such as the degree of obscuration and dust
temperature. This is especially true in the era of JWST where the much higher
data quality of NIRSpec/MIRI (e.g. Wells et al., 2015; Wright et al., 2015;
Rieke et al., 2015; Jakobsen et al., 2022; Argyriou et al., 2023; Böker et
al., 2023) over previous instruments such as Spitzer or Akari, means the
current models can struggle to fit the data due to an increase in spectral
resolution and sensitivity. Additionally, the significant improvement in
spatial resolution enables better isolation of the nuclei of galaxies, and
therefore the spectra are less diluted by star-formation. As a consequence,
the spectra are more dominated by a complex dust continuum, requiring new
models to reproduce.

In this work we present a new model, where we describe the continuum with a
generalised differential extinction model to fit combined NIRSpec and MIRI
spectra from 1.5 $\mu$m to 28 $\mu$m, where dust of different temperatures are
subject to different levels of extinction. The structure of this paper is as
follows. We first describe the model in detail in Section 2. We then test with
simulated data in Section 3 to verify the model before applying it to JWST
data in Section 4. Subsequently we compare the dust extinction model with
other tracers of extinction, namely the stellar continuum, Hydrogen
recombination lines and H2 lines for the 4 LIRGs observed with NIRSpec and
MIRI and discuss the inferred dusty structure in Section 5.

##  2 Model

The simplest way to predict the observed spectrum of a galaxy, $f_{\nu}$, is
to model the intrinsic spectrum, $f_{\nu}^{\textrm{Int}}$, with a sum of
modified blackbodies in the infrared, which is then subject to a single screen
of extinction of the form

$\frac{f_{\nu}}{f_{\nu}^{\textrm{Int}}}=e^{-\tau_{9.8}\tau(\lambda)}$ (1)

where $\tau(\lambda)$ is an assumed extinction law that is normalised to 1 at
9.8 $\mu$m. The level of extinction is then controlled by the optical depth at
9.8 $\mu$m, $\tau_{9.8}$, which scales the extinction law to increase the
extinction. This is equivalent to $A_{v}$, which is often used for rest-frame
UV/optical data. For reference $\tau_{9.8}\sim 0.1A_{v}$ (e.g. Draine, 1989).
While mathematically simple, the screen geometry is not well motivated, as the
obscuring dust and source dust are assumed to be separate, where the obscuring
dust does not to emit (see Fig. 11).

Alternatively, one can assume a mixed geometry where the emitting and
obscuring dust is uniformly mixed (e.g. Smith et al., 2007), which has the
form

$\frac{f_{\nu}}{f_{\nu}^{\textrm{Int}}}=\frac{1-e^{-\tau_{9.8}\tau(\lambda)}}{\tau_{9.8}\tau(\lambda)}.$
(2)

This is more physically motivated for star-forming regions compared to a
simple screen, as demonstrated in Fig. 11. These simple approaches tend to
fail however when modelling highly obscured galaxies, particularly over a long
wavelength range, where the emergent spectrum at $<5$ $\mu$m may be from a
much less obscured environment (with a different intrinsic spectrum) compared
to longer wavelength emission. A simple solution is to model two components, a
relatively unobscured star-forming component and an obscured nucleus, such as
that implemented in PAHDecomp (Donnan et al., 2023b). This is also the
approach of Hernán-Caballero et al. (2016) and references therein. This
however does not work for sources where there may be a gradient of extinction
as we probe different layers as a function of wavelength. An example where
such a geometry may be present is shown in 11. Therefore in order to be able
to fit the spectrum of any (U)LIRG, we require a more general model, to
capture this “differential extinction” behaviour. This is described in detail
in Section 2.1.1.

In addition to the differential extinction model, we model the PAH emission
features following a similar procedure to that of PAHFIT (Smith et al., 2007)
and include additional molecular absorption features (e.g. Imanishi et al.,
2010; Lai et al., 2020, 2023). Our full model, $f_{\nu}$, is given by

$f_{\nu}=\sum_{i=1}^{N_{\rm PAH}}I_{\nu,\rm
PAH}^{(i)}(\lambda)+C_{\nu}e^{-\tau_{\rm
ices(Dust)}(\lambda)}+S_{\nu}e^{-\tau_{\rm ices(Stellar)}(\lambda)}$ (3)

where the first term is a sum over $N_{\rm PAH}$ PAH profiles $I_{\nu,\rm
PAH}(\lambda)$ and the second is the total continuum comprising of the dust
model incorporating differential extinction, $C_{\nu}$, and a stellar
continuum, $S_{\nu}$, both of which are subject to a screen of the following
absorption features. For the dust continuum we include absorption from water
ice at $\sim$ 3 $\mu$m and 6$\mu$m, CO2 at 4.0 $\mu$m and CH at $\sim$ 7$\mu$m
while for the stellar continuum we only apply the near-IR absorption features
of water ice at $\sim$ 3 $\mu$m and CO2 at 4.0 $\mu$m. These are described in
detail in the following sections.

We show five examples of fits with this model to a variety of spectra in Fig.
1 to demonstrate the ability to fit a range of different objects with varying
obscuration. The components comprising the model in equation 3 are shown.

To fit this model to data we use a Bayesian approach, where we maximise the
following log probability

$\ln\textrm{Prob}=-\chi^{2}-\sum\ln\sigma^{2}-\Gamma P+\textrm{const.},$ (4)

where $\chi^{2}$ is simply the sum of the squared normalised residuals between
the model and the data, $\sigma$ is the data error bars, and $P$ is a penalty
for complex models of the dust distribution, scaled by a regularisation scale,
$\Gamma$. This factor controls the level of regularisation, where a larger
value increases the penalty for non-smooth dust distributions. This term is
described in detail in Section 2.1.2.

Figure 1: Five examples of fits to combined NIRSpec and MIRI spectra for
different types of nuclei. From top to bottom is NGC 7469 (type 1 AGN - see
Section 5.4.1), VV114 SW (Potential Type 2 AGN (Rich et al., 2023) \- see
Section 5.4.2), IIZw96 SW (Obscured AGN (García-Bernete et al., 2024b) \- see
Section 5.4.4), VV114 NE (Obscured AGN (Donnan et al., 2023a) \- see Section
5.4.2) and NGC 3256 S (Type 2 (Compton Thick) AGN (Ohyama et al., 2015) \- see
Section 5.4.3). Note the CO band at $\sim$ 4.6 $\mu$m is masked in the
fitting.

###  2.1 Continuum

####  2.1.1 Differential Extinction

A general model of dust emission and extinction can be written as a 2D
weighted average of modified blackbodies, $\epsilon_{\nu}B_{\nu}$, and screens
of extinction, $e^{-\tau_{9.8}\tau(\lambda)}$, where the weights are
parameterised by a 2D distribution, $\Psi(T,\tau_{9.8})$, of dust temperature,
$T$, and extinction, $\tau_{9.8}$. The modified blackbodies are given by the
product of the Planck function, $B(T)$, and a wavelength dependent emissivity
term, $\epsilon_{\nu}$ (see Section 2.1.3). The dust continuum is therefore
given by

$C_{\nu}(\lambda)=A\int_{0}^{\infty}\int_{35\textrm{K}}^{1500\textrm{K}}\Psi\left(T,\tau_{9.8}\right)e^{-\tau_{9.8}\tau(\lambda)}\frac{\epsilon_{\nu}B_{\nu}(T,\lambda)}{B_{0}(T)}dTd\tau_{9.8},$
(5)

where the sum of the weights, $\Psi(T,\tau_{9.8})$, is unity and $B_{0}(T)$ is
a normalisation factor for each modified blackbody. We choose $B_{0}(T)$ to be
the integrated flux between 1.5$\mu$m and 28$\mu$m with some emissivity,
$\epsilon_{\nu}$. We do not consider temperatures below 35K as dust this cold
will not contribute significantly to the flux below 30 $\mu$m. We use a
maximum temperature of 1500K as this is approximately the average dust
sublimation temperature for a mix of silicate and graphite grains (e.g.
Kishimoto et al., 2007).

In this work we assume a fixed extinction law, $\tau(\lambda)$, using that
presented in Donnan et al. (2023a), which consists of a power law and
absorption from silicates at $\sim 9.8$ and $\sim 18$ microns. We choose this
law as the silicate features include absorption from crystalline features,
necessary for highly obscured spectra (e.g. Spoon et al., 2022; Donnan et al.,
2023a). As the weighted average describes only the shape of the model, the
entire continuum is scaled by a factor $A$ to match the data.

In the simplest case where there is a single screen of extinction, the dust
distribution has a form

$\Psi(T,\tau^{{}^{\prime}}_{9.8})=\rho(T)\delta\left(\tau^{{}^{\prime}}_{9.8}-\tau_{9.8}\right),$
(6)

where $\rho(T)$ describes the distribution of blackbodies of different
temperatures (see Fig. 11). The single screen is represented by a Dirac delta
function which simplifies the integral into

$C_{\nu}(\lambda)=e^{-\tau_{9.8}\tau(\lambda)}\int_{35\textrm{K}}^{1500\textrm{K}}\rho(T)\frac{\epsilon_{\nu}B_{\nu}(T,\lambda)}{B_{0}(T)}dT$
(7)

where the observed continuum, $C_{\nu}(\lambda)$, is the intrinsic multiplied
by a factor of $e^{-\tau_{9.8}\tau(\lambda)}$, recovering equation (1).

For the mixed case there is a uniform dust distribution between 0 and some
$\tau_{9.8}$ (see Fig. 11):

$\Psi(T,\tau^{{}^{\prime}}_{9.8})=\rho(T)\begin{cases}\frac{1}{\tau_{9.8}}&0\leq\tau^{{}^{\prime}}_{9.8}\leq\tau_{9.8}\\\
0&\textrm{otherwise}\end{cases}$ (8)

which simplifies the integral to

$C_{\nu}(\lambda)=\frac{1-e^{-\tau_{9.8}\tau(\lambda)}}{\tau_{9.8}\tau(\lambda)}\int_{35\textrm{K}}^{1500\textrm{K}}\rho(T)\frac{\epsilon_{\nu}B_{\nu}(T,\lambda)}{B_{0}(T)}dT,$
(9)

where the observed continuum is the intrinsic continuum multiplied by a factor
$\frac{1-e^{-\tau_{9.8}\tau(\lambda)}}{\tau_{9.8}\tau(\lambda)}$, recovering
equation (2). This is the original PAHFIT model (Smith et al., 2007), where
the amplitude of different blackbodies is allowed to vary, comprising
$\rho(T)$, and is then subject to the extinction factor.

To implement our differential extinction model, instead of fixing the dust
distribution, $\Psi\left(T,\tau_{9.8}\right)$, we constrain this from the data
by fitting the continuum with a flexible model for
$\Psi\left(T,\tau_{9.8}\right)$. We do this under some prior restrictions, to
prevent nonphysical solutions for $\Psi\left(T,\tau_{9.8}\right)$, by ensuring
the distribution is smooth and only applies differential extinction where
necessary to fit the data.

To implement this model one could choose a parametric form of
$\Psi\left(T,\tau_{9.8}\right)$ such as a 2D Gaussian or a 2D Lorentzian,
however there is no obvious/physically motivated choice for the form of this
distribution. Assuming a given parametric form may therefore bias the results.
We therefore opt for a non-parametric approach where we bin
$\Psi\left(T,\tau_{9.8}\right)$ into a 20x20 grid and allow each grid element
to vary individually. This approach is similar to the stellar population
modelling of pPXF (Cappellari, 2023), where a non-parametric 2D distribution
of stellar metalicity and age is used to generate the model continuum.

We choose the grid carefully to ensure that each grid element gives a similar
weight to the continuum. We bin $\tau_{9.8}$ into 20 evenly spaced values in
natural log space between, $\tau_{9.8}=0.05$ and $\tau_{9.8}=15$. By binning
in $\ln$ space, we are evenly sampling different screens of extinction of the
form $e^{-\tau_{9.8}\tau(\lambda)}$. For the temperature axis we bin in
$\log_{10}$ space between 35K and 1500K. To ensure each grid element has equal
weighting when fitting, we normalise the spectrum of each grid element by the
integrated flux over the wavelength range of the data. After fitting, we re-
scale each grid element such that the normalisation matches that in equation 5
to obtain the $\Psi\left(T,\tau_{9.8}\right)$ in that equation.

The wavelength dependent extinction factor is given by the ratio of the
observed to the intrinsic continuum such as those shown in equation (1) and
(2) for the screen and mixed cases respectively. From our model, the intrinsic
continuum is defined as

$f_{\nu}^{\textrm{Int}}=S_{\nu}^{\textrm{Int}}+A\int_{0}^{\infty}\int_{35\textrm{K}}^{1500\textrm{K}}\Psi\left(T,\tau_{9.8}\right)\frac{\epsilon_{\nu}B_{\nu}(T,\lambda)}{B_{0}(T)}dTd\tau_{9.8},$
(10)

where $S_{\nu}^{\textrm{Int}}$ is the stellar population model defined in
equation (14). From this we calculate the wavelength dependent extinction
factor

$\frac{f_{\nu}}{f_{\nu}^{\textrm{Int}}}=\frac{S_{\nu}+C_{\nu}}{f_{\nu}^{\textrm{Int}}}$
(11)

where $S_{\nu}$ and $C_{\nu}$ are defined in equations (5) and (14)
respectively while $f_{\nu}^{\textrm{Int}}$ is given above. In general
$\frac{f_{\nu}}{f_{\nu}^{\textrm{Int}}}\neq e^{-\tau_{9.8}\tau(\lambda)}$, as
this only holds for a single screen geometry. Rather, the extinction factor
depends on the inferred $\Psi\left(T,\tau_{9.8}\right)$ and therefore can vary
between objects. For example, if there is a buried hot dust component compared
to the cool dust, the extinction factor will drop at $\sim 5$ $\mu$m compared
to long wavelengths. This is demonstrated with simulated data in Section 3.

####  2.1.2 Regularisation

Allowing each grid element on a 20x20 grid to vary, effectively adds 400
parameters to the model (typical NIRSpec+MIRI spectra contain $\sim$ 20000
data points). This provides the flexibility to fit complex, obscured galaxy
spectra, however the best fit solution for $\Psi\left(T,\tau_{9.8}\right)$ may
not be physically meaningful with so many parameters. While we find MCMC
sampling returns reasonably smooth solutions compared to a simple maximisation
of the $\ln\textrm{Prob}$ (equation (4)), we include some regularisation to
favour solutions where the distribution is simple/smooth.

We follow the approach of Cappellari (2017, 2023) where we add a penalty
factor to the log-likelihood to disfavour complex models of
$\Psi\left(T,\tau_{9.8}\right)$. In particular we use the sum of the Laplacian
for each row and column of $\Psi$:

$P=\int_{0}^{\infty}\int_{35\textrm{K}}^{1500\textrm{K}}\Psi^{\prime\prime}\left(T,\tau_{9.8}\right)dTd\tau_{9.8}=\\\
\sum_{i}\left(\Psi_{i-1,j}-2\Psi_{i,j}+\Psi_{i+1,j}\right)^{2}+\sum_{j}\left(\Psi_{i,j-1}-2\Psi_{i,j}+\Psi_{i,j+1}\right)^{2}$
(12)

where the integral of the 2nd derivative of $\Psi$ is calculated numerically
on the 20x20 grid in each axis as shown by the second line of the equation.
This factor must scale with the number of data points to be comparable to the
other terms in equation (4). We therefore choose $\Gamma=1000N$, where there
are $N$ data points in a given spectrum.

In addition to the preference for smooth solutions, we include additional
regularisation on the “effective” $\tau_{9.8}$ as function of wavelength.
During the fit, the effective $\tau_{9.8}$ is calculated for 20 wavelength
bins across the spectrum, evenly spaced in $\log_{10}\lambda$. This is done by
computing the extinction factor
($f_{\nu}^{\textrm{Obs}}/f_{\nu}^{\textrm{Intr}}$) at each wavelength bin and
converting to an effective $\tau_{9.8}$ where
$\tau_{9.8}^{\textrm{Eff}}=-\ln\left(f_{\nu}^{\textrm{Obs}}/f_{\nu}^{\textrm{Intr}}\right)$
at 9.8 $\mu$m. $\tau_{9.8}^{\textrm{Eff}}$ is then converted into units to
match the grid elements of $\Psi\left(T,\tau_{9.8}\right)$, such that we are
smoothing in steps of $e^{-\tau_{9.8}}$. We then place a Gaussian prior on the
differences between successive $\tau_{9.8}^{\textrm{Eff}}$ with a mean of zero
and a standard deviation of 6 grid elements. We chose this prior to be wide
such as to not limit the inference on $\Psi\left(T,\tau_{9.8}\right)$, but
restrictive enough to prevent the model fitting for differential extinction
where it is not needed. Therefore this prior ensures that the fit defaults to
no differential extinction and only deviates from this when necessary to fit
the data.

####  2.1.3 Emissivity

The shape of the blackbody profile for a given temperature, $T$, is modified
by an emissivity factor, $\epsilon_{\nu}$. Typically this factor takes the
form of a power law

$\epsilon_{\nu}\propto\lambda^{-\beta},$ (13)

with values of $\beta=1.5-2$. Codes such as PAHFIT (Smith et al., 2007), adopt
$\beta=2$, however this approach does not account for silicates, which
contribute significantly in the mid-infrared (e.g. Li & Draine, 2001).
Adopting a more physically motivated emissivity that includes silicates (as
well as graphite grains), results in a more realistic spectrum including
intrinsic silicate emission at 9.8 $\mu$m and 18$\mu$m (e.g. Li & Draine,
2001; Draine & Li, 2007; Marshall et al., 2007). However, silicate emission is
mainly observed in type 1 AGN (e.g. Gallimore et al., 2010; García-Bernete et
al., 2017; Martínez-Paredes et al., 2020; García-Bernete et al., 2022a), where
lines of sight of hot silicates, heated by the black hole accretion disk, are
visible without a significant hot background source, requiring clumpy or low-
inclination tori. In star-forming regions, the continuum is typically flat or
shows moderate silicate absorption. To achieve such a continuum with intrinsic
silicate emission requires either a cool dust temperature, where the silicates
contribute less and/or extra extinction to compensate. The problem is
therefore degenerate from an observational perspective, as the observed
continua for star-forming regions can be reproduced by either the power law
emissivity or a physically motivated emissivity.

In this work we explore both possibilities, starting with the simple power law
with $\beta=2$ and then with a emissivity following the curve presented in
Table 6 of Li & Draine (2001). To adjust to the spectral resolution of JWST,
we fit these points with a polynomial and the empirical silicate templates in
Donnan et al. (2023a). This is shown in Fig. 2 in comparison with the simple
power law.

Figure 2: Emissivity used to modify the intrinsic blackbody dust emission. The
blue line shows a simple power law with index $\beta=2$. The black points are
taken from Table 6 of Li & Draine (2001) and include silicate $\&$ carbon
grains. The red line shows our physical emissivity profile based on these
points by fitting a polynomial and a silicate template.

###  2.2 Stellar Continuum

In addition to the dust continuum, stellar emission will dominate below 3
$\mu$m. Fitting the stellar component is critical to constraining any possible
contribution from hot dust, $T\sim 1000K$, at around $\sim$ 5 $\mu$m where the
inclusion of NIRSpec to provide data at $<$ 5 $\mu$m provides strong
constraints on the stellar contribution. Particular features of stellar
continua can be seen such as CO absorption from stellar atmospheres at $\sim
2.2$ $\mu$m.

To model the stellar continuum we use two templates with ages of 10 Gyr and
100 Myr, at solar metallicity, provided by FSPS (Conroy et al., 2009; Conroy &
Gunn, 2010), with a Salpeter (1955) Initial Mass Function (IMF). At
wavelengths > 1 $\mu$m, the stellar continuum is dominated by more evolved
stellar populations in LIRGs (e.g. Pereira-Santaella et al., 2015), despite
these galaxies being highly star-forming. We chose these two templates based
on Pereira-Santaella et al. (2015), as combination of these templates provides
a reasonable stellar continuum without doing a full stellar population model,
which would be beyond the scope of this work. We allow each template to be
scaled individually to allow different combinations of young/older populations
which relates to the star-formation history. The scaled sum of the templates
is subject to its own screen of extinction, constrained by the shape of the
spectrum at $\sim$ 1.5 - 2 $\mu$m. Therefore the stellar component is simply

$S_{\nu}=e^{-\tau_{9.8}\tau(\lambda)}S_{\nu}^{\textrm{Int}},$ (14)

where $S_{\nu}^{\textrm{Int}}$ is the sum of each template scaled by its
corresponding scale factor.

###  2.3 Ices

In addition to the extinction present in the dust model, there are absorption
features from H2O ice, CO, CO2 and aliphatic CH present in the spectrum as
seen in equation (3). To account for this, we generate templates of the
optical depth from highly obscured sources and apply them as a screen to the
continuum. The templates for the mid-IR water ice and CH features at $\sim$ 6
$\mu$m are taken from Donnan et al. (2023b), which were generated from Spitzer
IRS spectra of the deeply embedded source, NGC 4418 (see also García-Bernete
et al., 2024a). We find the low resolution Spitzer data is sufficient to
generate a template for fitting the MIRI data as these features are broad.
However, the near-IR water ice at $\sim$ 3 $\mu$m and CO2 at $\sim$ 4.2 $\mu$m
were generated using the NIRSpec spectrum of the southern nucleus of NGC 3256.
This was necessary in order to create a reliable template as the depth of this
feature, with the high signal to noise of NIRSpec, much better captures the
shape compared to Akari data and the higher spatial resolution better isolates
the nucleus.

We do this by first creating a local continuum using a cubic spline with
anchor points at [2.5, 2.6, 2.65, 3.8, 4.18, 4.48, 4.82, 4.95] $\mu$m. This is
shown in red in the left panel of Fig. 3. We calculate an optical depth curve,
$\tau_{\textrm{ices}}(\lambda)$, using this local continuum,
$C_{\textrm{local}}$ where

$\tau_{\textrm{ices}}(\lambda)=\ln\left(f_{\textrm{data}}/C_{\textrm{local}}\right)$
(15)

with $f_{\textrm{data}}$ as the spectra. We mask the main emission features
before computing equation (15), most notably the 3.3 $\mu$m and 3.4 $\mu$m PAH
features. This results in the black points in the right panel of Fig. 3. To
convert this into a usable template we simply smooth the data for the CO2
feature. For the broad water ice feature we fit two 5th order polynomials to
each side of the feature.

Figure 3: Generation of water ice and CO2 templates from the southern nucleus
of NGC 3256. Left: NIRSpec spectrum between 2.5 and 5 $\mu$m shown in black.
The red line shows a local spline continuum, masking the main absorption
features of water ice, CO2 and CO. Right: Optical depth curve of the
absorption features using the local continuum in the left panel and masking
the main emission features in the data. The blue line shows the generated
template for water ice, created using two 5th order polynomials for each side
of the feature. The CO2 template in red was generated by smoothing the data.
We avoid generating a template for the CO feature and simply mask this region
of the spectrum in the fitting.

We do not create a template for the CO features at $\sim$ 4.6 $\mu$m as this
band consists of a variety of absorption/emission features which will change
shape from object to object and thus requires more complex modeling (e.g.
Pereira-Santaella et al., 2023; González-Alfonso et al., 2023; García-Bernete
et al., 2024b; Buiten et al., 2023). We therefore mask this region during the
fitting process. However we find the extrapolated best fit continuum over this
region provides a good fit to an underlying continuum (e.g. see Fig. 1).

We apply the near-infrared ices, namely the $\sim$ 3 $\mu$m water ice and CO2
at $\sim$ 4.2 $\mu$m, to both the dust continuum and the stellar continuum.
For both continua, the depth of the $\sim$ 3 $\mu$m water ice and the CO2
feature are allowed to vary independently. Moreover, the depths of each
feature are not necessarily the same for each continuum as a highly obscured
AGN, such as NGC 3256 S, shows much deeper ices than typical stellar continua
from star-forming regions. Therefore to obtain an adequate fit for NGC 3256 S,
we indeed require different optical depths for the ices applied to the stellar
vs dust continuum. This can be seen in the lower panel of Fig. 1. We therefore
allow the optical depth of the near-IR ices to vary separately as seen in
equation (3).

For most cases, there is not enough information in the data to constrain each
component separately (unlike NGC 3256 S), and therefore we place a Gaussian
prior that the depth of the ices is the same for the dust and stellar
continuum, such that these are only separate components when necessary to fit
the data. The prior is on the ratio of the two components with a mean of 1 and
a standard deviation of 0.2, such that $\pm$5$\sigma$ varies between zero and
two.

###  2.4 PAH Profiles

We model the PAH features with a series of Drude profiles, following the
original approach of Smith et al. (2007). As discussed in Donnan et al.
(2023a), we add additional profiles to be able to fit the higher spectral
resolution data from JWST. To fit the NIRSpec portion of the spectra, we
follow Lai et al. (2020), where the original PAHFIT was extended to fit AKARI
data.

To model asymmetric features, we introduce an asymmetry parameter to a single
Drude profile to better capture the shape of the feature with fewer parameters
than summing multiple components. This is applied to the 3.3, 3.4, 5.2, 6.2,
11.0, 11.3, 13.55 and 14.2 $\mu$m PAH profiles. We follow the approach of
Gordon et al. (2021) and Stancik & Brauns (2008) where we introduce parameter,
$a$, into the following equation for the

$I_{\nu,\rm
PAH}(\lambda)=A\frac{\left(\gamma/\lambda_{0}\right)^{2}}{\left(\lambda/\lambda_{0}-\lambda_{0}/\lambda\right)^{2}+\left(\gamma/\lambda_{0}\right)^{2}}$
(16)

where the original FWHM $\gamma_{0}$ is modified

$\gamma=\frac{2\gamma_{0}}{1+e^{a(\lambda-\lambda_{0})}}$ (17)

where $A$ is the amplitude and $\lambda_{0}$ is the central wavelength.

Similar to previous iterations of PAHFIT, we restrict the width, centre and
now the asymmetry parameters to a small range of values to prevent additional
degeneracies.

For the majority of the features the asymmetry is simply turned off. A full
list of these parameters can be found in Table. 2. The individual Drude
profiles can be seen in Fig. 1 for some example fits with the model.

We have additionally added some restrictions to the relative strengths of
certain PAH features. In particular, the less prominent PAH features, such as
those between the 6.2 $\mu$m and 7.7 $\mu$m cannot be too large relative to
the 6.2 $\mu$m and 7.7 $\mu$m PAHs. This is to prevent these PAHs fitting the
continuum when the PAH emission is weak. We therefore place the following
restrictions on the integrated PAH flux between certain limits: PAH[6.6,
7.3]/PAH[6.0, 6.5] < 1.0, PAH[6.6, 7.3]/PAH[7.4, 8.2] < 0.25, PAH[12.3,
13.0]/PAH[11.7, 12.2] < 2.8, PAH[11.7, 12.2]/PAH[11.2, 11.5] < 1.0, PAH[10.0,
11.1]/PAH[11.2, 11.5] < 0.4, where the values in the square brackets denote
the integration limits. These were chosen based on the star-forming spectra
where the PAH features are bright and well resolved.

###  2.5 Emission Lines

In the Spitzer era it was possible to model emission lines with simple
Gaussian profiles (e.g. Smith et al., 2007). This was necessary, not only to
obtain fluxes for each line but to properly decompose lines that may be
blended with PAH profiles. This was particularly important for the 12.7 $\mu$m
PAH feature which was heavily blended with the strong [Ne II] line at 12.81
$\mu$m.

An early modified version of PAHFIT for MIRI/MRS data was presented in Donnan
et al. (2023a) (also used in García-Bernete et al., 2022c) which modelled each
line with a 2 component profile to capture the more complex line profile
observed with the significantly higher spectral resolution of MIRI/MRS (R
$\sim$ 1500 - 3500 ) over Spitzer/IRS data (R $\sim$ 500) (e.g. Labiano et
al., 2021). Considering the abundance of lines now detected in typical MIRI
spectra (e.g. Pereira-Santaella et al., 2022; García-Bernete et al., 2022a; U
et al., 2022; Lai et al., 2022; Armus et al., 2023; Donnan et al., 2023a; Rich
et al., 2023; Young et al., 2023), and the number of parameters required per
line to capture the profile, this approach is extremely inefficient.
Fortunately, the high spectral resolution means that blending of emission
lines and PAH features is no longer an issue as the narrow widths of the
emission lines makes them easily distinguishable from the much broader PAH
features. We therefore simply mask the emission lines where they are detected
in the spectra, and integrate using the model subtracted data to obtain the
corresponding line fluxes.

Figure 4: Model subtracted residuals for the fit to NGC 3256 SF1 (grey; see
Section 4.2), to be used to integrate to obtain emission line fluxes. The red
dashed line marks a local continuum correction to the H2 S(6) line, which when
subtracted from the residuals, results in the red data points. This process is
applied to every emission line before calculating the integrated flux. The
vertical dashed lines show the range over which the line is integrated.

For some emission lines, particular those with low signal to noise, the model
subtracted data may contain some residuals which when integrated, will cause
inaccurate flux measurements. This is particularly relevant if the line is
faint. An example is shown in Fig. 4, where we show the model subtracted
residuals for two lines from a fit to a star-forming region in NGC 3256 (see
Section 4.2). The continuum around the H2 S(6) line (6.109 $\mu$m) appears
overestimated resulting in negative residuals redwards of the line, whereas
the continuum around first line (H2O $\nu$=1-0 212-101, 6.049 $\mu$m,
González-Alfonso et al., 2023; García-Bernete et al., 2024b) is predicted
accurately.

Before integrating to obtain a line flux, we correct the residuals by
measuring a local continuum around the line and subtracting it from the
residuals. We do this by calculating the median flux at either side of the
line and connecting the points linearly. This creates the red dashed line in
Fig. 4 which is subsequently subtracted to obtained the corrected residuals.
This is done for every emission line before integrating to obtain a flux. We
found that this correction led to underestimated flux errors as it effectively
corrects for scatter between each MCMC sample. To overcome this, we re-sample
the model subtracted residuals, for each MCMC sample, using the flux error
bars before integrating to obtain a flux.

##  3 Simulated Data

To assess the feasibility of this method in recovering the dust distribution
from the data, we test the model with simulated data. In particular, we want
to test that the model recovers accurately the input dust distribution. We
also want to test that the model recovers the simple screen and mixed cases
and thus does not fit for differential extinction where it does not exist.

###  3.1 Generating Mock Data

We generate three sets of mock data, a screen, mixed and a differential
distribution where there is a gradient of obscuration, with the hot dust
appearing more obscured than the cool dust. Specifically, we generated mock
spectra by inputting a dust distribution $\Psi\left(T,\tau_{9.8}\right)$ for
each of the three cases. For the differential case we use three 2D Gaussian’s
with parameters
$\left\\{x_{0},y_{0},\sigma_{x},\sigma_{y}\right\\}=\left\\{4,4,2,4\right\\}$,
$\left\\{8,8,2,4\right\\}$ and $\left\\{12,12,2,4\right\\}$ in units of the
grid elements with a position angle, of -40 ∘. For these tests we adopt the
simple power law emissivity profile (see Section 2.1.3). This is to allow a
more direct comparison to the traditional PAHFIT methods.

The dust distribution for the three sets of simulated data are shown in Fig. 5
with the corresponding continuum and extinction factor in the right panels. We
add a stellar continuum to simulate the near-IR portion of the spectra. We use
the 100 Myr template (we chose the 100 Myr arbitrarily here and test whether
the model recovers the correct stellar continuum in Section 3.2), and apply an
extinction equal to that of the dust continuum for the screen and mixed cases.
For the differential case, we apply a screen of extinction to the stellar
continuum, equivalent to the screen case.

We additionally add a PAH and emission line template to the continuum before
generating the mock data. This template was generated as the average of all
the continuum subtracted spectra of the star-forming regions in NGC 3256 and
NGC 7469 (see Section 4.2) after fitting with the original PAHFIT-like
continuum (Donnan et al., 2023a). We then generate the mock data with
wavelengths equal to that of the real data and noise of 1% of the flux values.
This gives a signal-to-noise comparable to the JWST data.

Each of the three simulated spectra are fitted with a fixed screen, fixed
mixed and the new differential extinction model to investigate how well the
true dust distribution and extinction factor can be recovered. The best fit
extinction factor is shown in the right panels of Fig. 5, in comparison with
the true, input extinction factor. The contours in the left most panel show
the best fit dust distribution from the differential model.

###  3.2 Case 1: Screen Dust Distribution

For an input screen dust distribution, the screen model (see the top panel of
Fig. 5) recovers the true extinction factor within 1 $\%$ over the entire
wavelength range as expected, while the mixed model performs the worst,
underestimating the extinction across the wavelength range by $\sim 5\%$ up to
$\sim 15\%$ at $\sim 4$ $\mu$m. The differential model recovers the true
extinction well within 1 $\%$ across the majority of the wavelength range,
only deviating at $\sim 4$ $\mu$m by a maximum of 4 $\%$. For all the fits we
find the stellar continuum is recovered accurately, with the input stellar
population of 100 Myr being recovered, as the near-IR is dominated by this
component and so there is no difficulty in reproducing the input stellar
continuum.

The recovered dust distribution in the left panel of Fig. 5 peaks in the
correct region but shows a larger spread. This larger spread is not seen for
the mixed input data where the best fit dust distribution matches the spread
accurately. Physically this is a bias towards more mixed dust distributions
and so while this does not change the resultant extinction factor/continuum in
any significant way, one must take care when drawing any inference about the
physical dusty structure from the shape of $\Psi\left(T,\tau_{9.8}\right)$.
This bias towards more smooth dust distributions is a consequence of the
regularisation with the non-parametric method. This is a known bias of non-
parametric approaches (e.g. Leja et al., 2019; Iyer et al., 2019), for example
with non-parametric SFH-histories, where sudden bursts in star-formation can
be physically expected but can difficult to reproduce with this approach.

###  3.3 Case 2: Mixed Dust Distribution

For the mixed input case we find the mixed model (see the middle panel of Fig.
5) performs best as expected, recovering the extinction factor within $1\%$.
The screen and differential models also perform well in this case, both
finding the true extinction factor within 1$\%$ with both only deviating by a
maximum of 4$\%$ at $\sim 1.6\mu$m. These tests highlight that the
differential extinction model is able to recover the extinction factor in the
simple cases where there is no differential extinction.

###  3.4 Case 3: Differential Dust Distribution

The third simulated spectrum was produced using a dust distribution with a
clear gradient where the hotter dust is more obscured than the colder dust
(see the bottom panel of Fig. 5). We find that the fixed screen and mixed
models fail to recover the extinction factor, and indeed fail to produce a
satisfactory fit to the continuum. The differential model however reproduces
the input dust continuum accurately and recovers the extinction factor well,
with a deviation at $\sim 5$ $\mu$m by up to 15$\%$. This deviation implies
the model has difficulty recovering the most obscured hot dust, which would
lower the extinction factor at $\sim 5$ $\mu$m. This is a limitation of the
data as these obscured blackbodies will have a minimal contribution to the
flux/shape of the model and thus makes little difference to the fit. Therefore
due to the regularisation of the model, the fit favours more conservative
solutions by not switching on components unnecessary to fit the data, pushing
the recovered extinction factor closer to a case with less extreme
differential extinction. This is therefore a hard limit of the data, where the
information is simply not present to infer the most obscured hot dust that was
included when generating the mock data.

Figure 5: Tests with simulated data for three input dust distributions: screen
(upper), mixed (middle) and differential (lower). The left panels show the
dust distribution, $\Psi(T,\tau_{9.8})$, used to generate the simulated data
overlayed with the best fit model in the grey contours. The colour scale is
linear, with the grid elements summing to one. The 2nd column shows the
resulting continuum from the input dust distribution (solid red line) and the
input stellar continuum (solid blue line). The simulated data is shown with
the grey points. The dashed lines show the different best fit models from the
screen fit (green), mixed (blue) and the differential model (red). The 3rd
column shows just the dust continuum from the simulated data in comparison
with the best fit models. The right panel shows the wavelength dependent
extinction factor for the input (black) and the best fit results when fitting
with a single screen (green dashed), pure mixed (blue dashed) and the new
differential extinction model (red dashed).

For the third case we see that the traditional fixed or screen models do not
properly fit the data and predict the extinction factor incorrectly, which
would lead to incorrect flux ratios from emission features. This highlights a
case where a more complex dust distribution is required and where our new
model is able to reproduce the input conditions. As discussed in Section 4.4,
dust distributions similar to this are present in LIRGs (and therefore also
ULIRGs), and so accounting for the effect of differential extinction is
necessary to model these spectra.

##  4 JWST Data

###  4.1 Targets

We use NIRSpec IFU and MIRI MRS data from Director’s Discretionary Early
Release Science Program 1328 (PI: Lee Armus & Aaron Evans) of four local
LIRGs, namely NGC 7469, NGC 3256, IIZw96 and VV 114. We show RGB NIRCAM images
in Fig. 6 for the 4 galaxies, with the F120W, F200W filters tracing stellar
continuum and the F335M/F356W filters tracing the 3.3 $\mu$m PAH emission. The
fully reduced files for the NIRCAM imaging were downloaded from the MAST
archive.

NGC 7469 is a barred spiral galaxy hosting a type 1 AGN (e.g. Osterbrock &
Martel, 1993) with a compact nuclear starburst and starburst ring (e.g. Díaz-
Santos et al., 2007; García-Bernete et al., 2022a; Lai et al., 2022, 2023).
The AGN drives a wide, highly ionised outflow towards the SE (e.g. García-
Bernete et al., 2022a; U et al., 2022; Bianchin et al., 2023), coinciding with
a gap in the star-forming ring. This object has been subject to much analysis
of the PAH population with JWST (García-Bernete et al., 2022a; Zhang & Ho,
2023b; Lai et al., 2023).

NGC 3256 is a late-stage merger consisting of a face-on spiral galaxy (NGC
3256 N Sakamoto et al., 2014) with a bright nuclear starburst (Lira et al.,
2008). The southern component of the merger (NGC 3256 S) is an edge on galaxy
hosting a highly obscured, Compton Thick AGN (e.g. Ohyama et al., 2015) which
drives a highly collimated outflow (e.g. Sakamoto et al., 2014; Pereira-
Santaella et al., 2023).

VV114 is a mid-stage major merger (e.g. Stierwalt et al., 2013) where the
eastern component appears extremely dusty outputting a large infrared
luminosity (Armus et al., 2009). The region is complex, with multiple nuclei
and a large scale shock resulting from the interaction (Saito et al., 2017;
Donnan et al., 2023a; Rich et al., 2011, 2023; González-Alfonso et al., 2023;
Buiten et al., 2023).

Finally, IIZw96 is also a major merger, with a dusty eastern region consisting
of numerous compact, dusty star-forming clumps (Inami et al., 2010, 2022),
separate from the two main nuclei of the interacting galaxies. The SW nucleus
is thought to host a buried AGN from modelling IR rovibrational lines of the
molecular bands (García-Bernete et al., 2024b).

###  4.2 Data Reduction

We reduced the data using modifications to the default pipeline. The NIRSpec
data were reduced using pipeline version 1.9.4 with the 1063 calibration
context files. The details of the reduction can be found in Pereira-Santaella
et al. (2023).

The MIRI data were reduced using the default pipeline (v 1.9.4) up to the
stage 2 data, using the 1041 calibration context. Before step 3 we run the
residual fringe correction step offline. We run stage 3 without the background
subtraction, instead performing this afterwards. We do this by calculating the
median in each background frame before subtracting it. More details of the
reduction are found in Appendix A of García-Bernete et al. (2024a).

We extract spectra through a 0.3" radius aperture. To ensure the region is
identical between the sub-channels we correct for small spatial misalignment’s
due to errors in the astrometry. We do this by measuring the position of the
aperture we wish to extract the spectrum from relative to the brightest point
source in the data. By measuring relative positions to a point source we avoid
re-sampling the entire cube which may introduce further errors into the data.
After extraction of the spectra, we run an additional residual fringe
correction on the 1D spectra to remove the remaining fringing. Additionally,
we use a PSF wavelength dependent correction factor for a point-like
extraction for the nuclear apertures (see also Pereira-Santaella et al., 2022;
García-Bernete et al., 2022a; Donnan et al., 2023a). This was obtained using a
standard star for both NIRSpec and MIRI where the PSF correction factor is
calculated as the ratio of the spectrum within 0.3" to the total spectrum as
obtained from CALSPEC (Bohlin et al., 2014; Bohlin & Lockwood, 2022). For
NIRSpec we use TYC 4433-1800-1 while for MIRI we use HD 163466.

To estimate errors for the fluxes at each wavelength we avoid those generated
by the pipeline for the MIRI data as we found them to be underestimated. We
instead use the dedicated background exposures and calculate the rms in each
frame. For the NIRSpec data, we found the pipeline produced reasonable errors.

Fig. 6 shows the positions of the apertures used in this work for each object.
We extract six star-forming regions from each of NGC 7469 and NGC 3256, chosen
based on near-IR bright regions. For NGC 7469, SF2, SF3 and SF4 lie on the
redshifted region of the nuclear outflow (e.g. García-Bernete et al., 2022a; U
et al., 2022; Bianchin et al., 2023) while a nuclear bar extends (see Fig. 10
of Díaz-Santos et al., 2007) from SF1 to the southwest.

Figure 6: Images of the 4 LIRGs used in this work and the apertures used to
extract the spectra. The upper panels show RGB images with NIRCAM. The
F335M/F356W filters (red) show the 3.3 $\mu$m PAH emission while the F120W and
F200W filters mainly trace stellar continuum. Most of the nuclei show some
saturation at the brightest points. The field of view of the NIRSpec data is
shown by the black box. The lower panels show the integrated NIRSpec detector
2 G395H/F290LP channel, integrated between 4.08 - 5.27 $\mu$m, tracing
continuum (stellar and/or dust). The apertures used in this work to extract
the spectra are shown and labelled. Each aperture has a radius of 0.3" and
those labelled SF do not have any PSF correction applied as their emission
appears extended.

For IIZw96 SW, there is a secondary nucleus to the west, with a very small
separation. At long wavelengths, the PSFs of each source overlap, leading to
contamination. To mitigate this effect we extract the spectra of each nuclei,
applying the aperture correction as discussed above. We then extract a total
spectrum using a 0.6" aperture centred on the midpoint between the two nuclei.
For this aperture we also apply a PSF correction factor by calculating the
flux loss within a 0.6" aperture centred between two standard stars. We find
the sum of the individual nuclei spectra to be larger than the 0.6" total
spectrum, particularly at long wavelengths, due to contamination by each of
the nuclei. To correct the SW nuclear spectrum, we subtract a fraction of the
excess flux (where excess = sum of the nuclei - total), where the fraction is
the contribution of the second nuclei to the sum of the nuclei. We find the
correction to be small but significant at long wavelengths. For reference this
is shown in Appendix C.

###  4.3 Other Estimates of Extinction

With the rich MIRI MRS and NIRSpec spectra, comes additional diagnostics of
extinction. These can then be compared to the differential extinction model of
the dust continuum to provide insight into the physical conditions of the
regions. In the following sections we discuss three different measures of
extinction using the emission features in the spectra.

####  4.3.1 Hydrogen Recombination Lines (HI)

Hydrogen recombination lines have been used in the optical to obtain
extinction measurements for decades (e.g. Salim & Narayanan, 2020). The so-
called Balmer decrement measures the extinction of HII regions around young
stars. This works by assuming an intrinsic H$\alpha$/H$\beta$ ratio under case
B conditions (Hummer & Storey, 1987) and inferring the optical depth required
to reproduce the observed ratio. In the near and mid-infrared the same
technique can be applied using the Brackett (Br, 4.05, 2.62, 2.17, … $\mu$m),
Pfund (Pf, 7.46, 4.65, 3.74, … $\mu$m) and Humphreys (Hu, 12.4, 7.50, 5.91 …
$\mu$m) series rather than the Balmer series.

While a large number of lines can be used, we focus on the ratio of Br$\beta$
(2.62 $\mu$m) to Br$\gamma$ (2.17 $\mu$m) as these lines have a high signal-
to-noise and are unblended with any other lines. A more detailed analysis
combining all the HI lines may provide further insight but is beyond the scope
of this work.

We predict the extinction for each spectrum by applying a screen with a given
$\tau_{9.8}$ that reproduces the observed flux ratio. We do this for each MCMC
sample to obtain uncertainties for the extinction.

We do not apply this to the nucleus of NGC 7469 as there is significant
contamination of the HI recombination lines from the broad line region of the
AGN. Therefore the intrinsic flux ratio of these lines is no longer due to
case B in HII regions and so would provide inaccurate extinction values.

####  4.3.2 Molecular Gas (H2)

The mid-infrared is populated with numerous H2 lines tracing the molecular
gas. In particular the S(1) to S(8) transitions between 5 $\mu$m and 28
$\mu$m. In Donnan et al. (2023a) and Hernandez et al. (2023), it was noted
that the S(3) J=5 transition is extremely sensitive to extinction as its
central wavelength at 9.665 $\mu$m places it within the trough of the 9.8
$\mu$m silicate absorption (see also Pereira-Santaella et al., 2014). In the
referenced works, this line was largely ignored as a result, however this fact
also makes it ideal for predicting level of extinction affecting the molecular
gas. We do this by first constructing a rotation diagram where we plot the
population of each energy level, $\log(N/g)$, against the energy of the upper
level, $E_{u}/k$, as shown in the left panel of Fig. 7. This is for a typical
star forming region (NGC 3256 SF4), where the fluxes used have not been
corrected for extinction (black points). As a result, the S(3) transition
appears lower than expected. We can therefore estimate the extinction by
varying the correction factor until this dip no longer exists.

Figure 7: Inferring the extinction of H2, in this case for NGC 3256 SF4, by
varying the extinction correction and measuring the smoothness of the
resultant rotation diagram. Left: Rotation diagram for H2 for the observed
fluxes (black) and the extinction corrected fluxes (red). Note the change for
the S(3) transition. Right: Penalty factor for the smoothness of the resulting
rotation diagram given an extinction correction with optical depth
$\tau_{9.8}$. The best fit value is given by the minimum of this curve which
is marked with the solid vertical red line. The vertical dashed red liens show
the 1$\sigma$ error range.

We do this by varying the optical depth of a single screen of dust and
calculating the extinction corrected rotation diagram for all the transitions.
For each optical depth, we calculate the smoothness of the resulting rotation
diagram using the sum of the Laplacian, which is the 1D version of equation
(12),

$P_{H_{2}}=\sum_{i}\left[\log(N_{i-1}/g_{i-1})-2\log(N_{i}/g_{i})+\log(N_{i+1}/g_{i+1})\right]^{2}$
(18)

where, $\log(N_{i}/g_{i})$ is the population of transition $i$. The inferred
extinction is that where the smoothness penalty, $P_{H_{2}}$ is minimised.
This is demonstrated in the right panel of Fig. 7 with the corresponding
extinction corrected rotation diagram in the left panel in red. As before we
can perform this for each MCMC sample to get an uncertainty on the inferred H2
extinction.

The gradient of the rotation diagram is related to the temperature of the gas,
where steeper curves result from cooler gas as the higher energy transitions
are relatively less populated. Typically one will fit multiple ($\sim$ 2 or 3)
linear components, each with a different temperature to fit the curve (e.g.
Rigopoulou et al., 2002; Pereira-Santaella et al., 2022) or assume some power
law distribution of temperature (Togi & Smith, 2016; Pereira-Santaella et al.,
2014). To obtain an estimate without assuming a model, we calculate the median
gradient of pairs of transitions to provide an average molecular gas
temperature as traced by the mid-infrared.

####  4.3.3 PAHs

If an intrinsic ratio between two PAH bands is assumed, the extinction can be
inferred by the same method as the Hydrogen recombination lines. This method
was proposed by Hernán-Caballero et al. (2020) where they showed that the
12.7/11.3 PAH ratio is intrinsically constant for star-forming galaxies with
the observed ratio depending only on the extinction.

We calculate the intrinsic 12.7/11.3 PAH ratio using the star-forming regions
of NGC 3256 and NGC 7469 as a calibration sample. This step is required as the
increased spectral resolution of MIRI MRS over Spitzer IRS, leads to a more
resolved 12.7 PAH feature (e.g. Donnan et al., 2023a). We measure an intrinsic
12.7/11.3 PAH ratio of $0.72\pm 0.08$, after correcting for extinction using
the average continuum extinction as given in Table 1. As we rely on an
empirical calibration, the resulting measured PAH extinctions will largely
follow that of the continuum for the star-forming regions by construction.
Therefore the measured PAH extinction is not an completely independent measure
but is valuable as a comparison between objects.

Taking the samples of the intrinsic ratio, we vary the extinction until the
observed 12.7/11.3 PAH ratio is reproduced for each MCMC sample of the
measured ratio, to obtain a posterior for the PAH extinction.

###  4.4 Results

We fit the model to each spectrum extracted as discussed in Section 4.2. We
fit the entire sample with each emissivity, first with the simple power law
with $\beta=2$ and then with the physical emissivity profile which includes a
mixture of silicates and graphite based on Li & Draine (2001). We found the
physical emissivity profile provides a better fit in all cases. Moreover it is
necessary to fit the silicates seen in emission for the nucleus of NGC 7469
(e.g Alonso-Herrero et al., 2020; García-Bernete et al., 2022a). We compare
the choice of emissivity in detail in Section 5.2. The results hereafter are
with the physical emissivity, as it provides a better fit and is better
motivated from a theoretical perspective.

We show the best fit dust distributions for all the nuclei and a
representative star-forming region in Fig. 8. The remaining star-forming
regions are shown in Appendix 13. We use the physical emissivity as we found
this to provide a better fit to the data and is more physically motivated.

We find the shape of the dust distribution, $\Psi\left(T,\tau_{9.8}\right)$,
to be largely consistent between the various star-forming regions. In
particular it shows a relatively low dust temperature ($\sim$ 90 - 100 K) that
is lightly obscured ($\tau_{9.8}\sim 0.2$), resulting in a steep dust
continuum with little to no silicate absorption. There additionally appears to
be a faint, obscured hot dust component with $\tau_{9.8}>1$. A few of the
star-forming regions, NGC 7469 SF1, SF3 and SF4, show a stronger gradient,
with a stronger hot dust component resulting in moderate silicate emission.
Comparing with the other estimates of extinction, the molecular gas appears
the most obscured, even more obscured than the dust continuum. We discuss
possible reasons for this in Section 5.1.

As expected, the nuclei generally show more complex dust distribution, where
multiple dust components contribute to the observed spectrum. We discuss each
individually in detail in Section 5.4.

We display the inferred properties of each fitted spectrum in Table 1. The
temperature, $T$, and extinction, $\tau_{9.8}$, are given for the total
continuum and individual components of the dust distribution for those cases
as described in Section 5.3. To obtain this, we sum over one of the axes of
the dust distribution to collapse the 2D distribution into a 1D distribution
either of temperature, $T$, and extinction, $\tau_{9.8}$ (by summing over the
other axis). We then determine the median and 16$\%$/84$\%$ upper and lower
bounds of the 1D distributions. These values reflect the shape/spread of
$\Psi(T,\tau_{9.8})$ rather than the statistical uncertainties, which are
negligible by comparison.

We additionally show the measured extinction from the other tracers as
outlined in Section 4.3. For the molecular gas we also show the measured gas
temperature.

For NGC 7469 Nuc and IIZw96 SW, the 12.7/11.3 PAH ratio is not well
constrained, and therefore cannot measure the extinction accurately.

We note that the uncertainties on the inferred properties are dominated by the
priors on the intrinsic PAH ratio as the statistical uncertainties are very
small due to the high signal-to-noise of the data. This is reflected in the
errors on the stellar extinction. These are likely underestimated as we have
not folded in errors based on the assumption of the chosen stellar population
templates, which for example, assumed solar metallicity.

Table 1: Inferred properties from each spectrum

Spectrum | Total Cont. | SF Cont. | AGN Warm | AGN Hot | HI | H2 | PAHs | Stellar
---|---|---|---|---|---|---|---|---
| $\tau_{9.8}$ | $T$ | $\tau_{9.8}$ | $T$ | $\tau_{9.8}$ | $T$ | $\tau_{9.8}$ | $T$ | $\tau_{9.8}$ | $\tau_{9.8}$ | $\bar{T}$ | $\tau_{9.8}$ | $\tau_{9.8}$
|  | K |  | K |  | K |  | K |  |  | K |  |
NGC 7469 SF1 | $0.49^{+1.2}_{-0.38}$ | $111.0^{+120.0}_{-50.0}$ | - | - | - | - | - | - | 0.56${}^{+0.17}_{-0.18}$ | 0.86${}^{+0.05}_{-0.074}$ | 716.0${}^{+100.0}_{-100.0}$ | 0.88${}^{+0.27}_{-0.23}$ | 0.3767${}^{+0.0033}_{-0.0031}$
NGC 7469 SF2 | $0.43^{+1.6}_{-0.34}$ | $100.0^{+79.0}_{-41.0}$ | - | - | - | - | - | - | 1.09${}^{+0.13}_{-0.14}$ | 0.72${}^{+0.084}_{-0.072}$ | 664.0${}^{+69.0}_{-74.0}$ | 0.56${}^{+0.28}_{-0.24}$ | 0.3513${}^{+0.0039}_{-0.0039}$
NGC 7469 SF3 | $0.40^{+1.2}_{-0.3}$ | $104.0^{+66.0}_{-40.0}$ | - | - | - | - | - | - | 1.56${}^{+0.29}_{-0.28}$ | 1.11${}^{+0.096}_{-0.11}$ | 654.0${}^{+54.0}_{-56.0}$ | 0.54${}^{+0.25}_{-0.24}$ | 0.2436${}^{+0.0036}_{-0.0037}$
NGC 7469 SF4 | $0.41^{+1.4}_{-0.31}$ | $100.0^{+65.0}_{-38.0}$ | - | - | - | - | - | - | 1.20${}^{+0.20}_{-0.28}$ | 0.87${}^{+0.11}_{-0.13}$ | 660.0${}^{+76.0}_{-70.0}$ | 0.55${}^{+0.27}_{-0.21}$ | 0.2401${}^{+0.0027}_{-0.0027}$
NGC 7469 SF5 | $0.50^{+1.9}_{-0.4}$ | $93.0^{+75.0}_{-38.0}$ | - | - | - | - | - | - | 0.55${}^{+0.15}_{-0.21}$ | 1.1${}^{+0.11}_{-0.084}$ | 629.0${}^{+45.0}_{-46.0}$ | 0.70${}^{+0.24}_{-0.24}$ | 0.2578${}^{+0.0022}_{-0.003}$
NGC 7469 SF6 | $0.53^{+1.7}_{-0.41}$ | $105.0^{+87.0}_{-43.0}$ | - | - | - | - | - | - | 0.55${}^{+0.15}_{-0.16}$ | 1.21${}^{+0.084}_{-0.084}$ | 654.0${}^{+46.0}_{-42.0}$ | 0.49${}^{+0.25}_{-0.24}$ | 0.3233${}^{+0.0059}_{-0.0056}$
NGC 7469 Nuc | $1.02^{+2.8}_{-0.89}$ | $154.0^{+62.0}_{-40.0}$ | - | - | $1.04^{+3.0}_{-0.92}$ | $149.0^{+47.0}_{-37.0}$ | $0.97^{+0.75}_{-0.48}$ | $634.0^{+390.0}_{-160.0}$ | - | 0.34${}^{+0.012}_{-0.024}$ | 712.0${}^{+5.9}_{-5.9}$ | - | -
NGC 3256 SF1 | $0.33^{+1.5}_{-0.25}$ | $87.0^{+68.0}_{-33.0}$ | - | - | - | - | - | - | 0.85${}^{+0.28}_{-0.29}$ | 1.26${}^{+0.048}_{-0.06}$ | 538.0${}^{+92.0}_{-21.0}$ | 0.15${}^{+0.29}_{-0.23}$ | 0.1862${}^{+0.0035}_{-0.0033}$
NGC 3256 SF2 | $0.55^{+2.2}_{-0.44}$ | $90.0^{+68.0}_{-33.0}$ | - | - | - | - | - | - | 1.10${}^{+0.24}_{-0.22}$ | 1.45${}^{+0.048}_{-0.048}$ | 684.0${}^{+29.0}_{-34.0}$ | 0.55${}^{+0.24}_{-0.34}$ | 0.4949${}^{+0.0058}_{-0.0059}$
NGC 3256 SF3 | $0.38^{+1.6}_{-0.29}$ | $91.0^{+62.0}_{-33.0}$ | - | - | - | - | - | - | 0.72${}^{+0.16}_{-0.15}$ | 1.43${}^{+0.036}_{-0.024}$ | 722.0${}^{+20.0}_{-24.0}$ | 0.64${}^{+0.26}_{-0.24}$ | 0.3288${}^{+0.0034}_{-0.0035}$
NGC 3256 SF4 | $0.41^{+1.7}_{-0.33}$ | $89.0^{+69.0}_{-36.0}$ | - | - | - | - | - | - | 0.83${}^{+0.19}_{-0.19}$ | 1.39${}^{+0.072}_{-0.084}$ | 634.0${}^{+58.0}_{-77.0}$ | 0.39${}^{+0.26}_{-0.24}$ | 0.2575${}^{+0.0052}_{-0.0052}$
NGC 3256 SF5 | $1.07^{+1.9}_{-0.81}$ | $116.0^{+84.0}_{-46.0}$ | - | - | - | - | - | - | 1.01${}^{+0.060}_{-0.064}$ | 1.63${}^{+0.024}_{-0.024}$ | 765.0${}^{+14.0}_{-8.9}$ | 0.40${}^{+0.28}_{-0.22}$ | 0.4173${}^{+0.004}_{-0.0041}$
NGC 3256 SF6 | $0.47^{+1.7}_{-0.38}$ | $93.0^{+61.0}_{-33.0}$ | - | - | - | - | - | - | 0.81${}^{+0.17}_{-0.17}$ | 1.33${}^{+0.048}_{-0.048}$ | 658.0${}^{+29.0}_{-52.0}$ | 0.68${}^{+0.25}_{-0.25}$ | 0.3821${}^{+0.0045}_{-0.0043}$
NGC 3256 N | $1.07^{+1.3}_{-0.83}$ | $109.0^{+62.0}_{-29.0}$ | - | - | - | - | - | - | 1.35${}^{+0.036}_{-0.036}$ | 2.38${}^{+0.012}_{-0.012}$ | 436.0${}^{+19.0}_{-3.2}$ | 1.18${}^{+0.26}_{-0.23}$ | 0.2937${}^{+0.0017}_{-0.0012}$
NGC 3256 S | $4.66^{+6.4}_{-3.6}$ | $83.0^{+43.0}_{-27.0}$ | $0.66^{+0.61}_{-0.54}$ | $77.0^{+30.0}_{-22.0}$ | $6.1^{+6.4}_{-3.4}$ | $83.0^{+40.0}_{-27.0}$ | $6.6^{+2.7}_{-1.8}$ | $1274.0^{+230.0}_{-290.0}$ | 2.49${}^{+0.040}_{-0.036}$ | 2.24${}^{+0.012}_{-0.024}$ | 540.0${}^{+38.0}_{-39.0}$ | 2.07${}^{+0.27}_{-0.22}$ | 1.3971${}^{+0.0018}_{-0.0018}$
VV114 NE | $5.27^{+4.4}_{-2.6}$ | $93.0^{+68.0}_{-32.0}$ | $1.02^{+0.46}_{-0.88}$ | $64.0^{+23.0}_{-14.0}$ | $5.48^{+4.4}_{-2.5}$ | $95.0^{+68.0}_{-33.0}$ | - | - | 1.76${}^{+0.058}_{-0.059}$ | 1.98${}^{+0.026}_{-0.012}$ | 591.0${}^{+11.0}_{-7.3}$ | 2.03${}^{+0.27}_{-0.23}$ | 0.7016${}^{+0.00056}_{-0.00056}$
VV114 SW | $0.95^{+2.1}_{-0.81}$ | $84.0^{+100.0}_{-29.0}$ | $0.25^{+0.36}_{-0.17}$ | $71.0^{+27.0}_{-20.0}$ | $2.1^{+2.2}_{-1.0}$ | $93.0^{+66.0}_{-33.0}$ | $2.29^{+1.4}_{-0.83}$ | $732.0^{+260.0}_{-200.0}$ | 1.46${}^{+0.038}_{-0.039}$ | 1.51${}^{+0.024}_{-0.024}$ | 656.0${}^{+92.0}_{-56.0}$ | 1.47${}^{+0.28}_{-0.22}$ | 0.7309${}^{+0.0046}_{-0.0043}$
IIZw96 NE | $0.69^{+1.8}_{-0.55}$ | $72.0^{+70.0}_{-30.0}$ | - | - | - | - | - | - | 1.11${}^{+0.14}_{-0.14}$ | 0.69${}^{+0.24}_{-0.31}$ | 461.0${}^{+120.0}_{-75.0}$ | 1.49${}^{+0.27}_{-0.24}$ | 0.5045${}^{+0.0037}_{-0.0038}$
IIZw96 SW | $2.25^{+3.3}_{-1.7}$ | $91.0^{+58.0}_{-30.0}$ | $0.69^{+0.56}_{-0.54}$ | $76.0^{+35.0}_{-21.0}$ | $3.59^{+3.2}_{-1.5}$ | $104.0^{+60.0}_{-37.0}$ | $13.11^{+1.9}_{-9.5}$ | $774.0^{+340.0}_{-240.0}$ | 1.30${}^{+0.020}_{-0.019}$ | 2.08${}^{+0.036}_{-0.024}$ | 902.0${}^{+9.2}_{-7.8}$ | - | 0.194${}^{+0.0096}_{-0.0090}$

Figure 8: Fits to eight of the spectra in this work. We show a star-forming
region (NGC 3256 SF1) and the 7 nuclei. Each panel shows the best fit dust
distribution, $\Psi\left(T,\tau_{9.8}\right)$, on the left and the
corresponding fit to the data in the right panel. The horizontal coloured
lines show the extinction measurements for the molecular gas, PAHs, HI lines
and the stellar continuum, as described in Section 4.3. The dashed grey lines
show a decomposition of the best fit distribution into star-forming and AGN
torus + polar dust for those case where applied. The resulting continuum for
each component is plotted in the right panels.

##  5 Discussion

###  5.1 Layers of Extinction

From Fig. 8 and Fig. 13, we find that the molecular gas and HII regions are
the most obscured while the PAHs and stellar continuum are less obscured. This
is consistent with a picture where the molecular gas is the most buried within
dust clouds which will collapse to form stars. The HII regions are similarly
obscured, resulting from newly formed stars that still remain deeply buried
within the star-forming region.

The dust continuum appears highly mixed ranging from unobscured
($\tau_{9.8}\sim 0.2$) to obscured ($\tau_{9.8}>1$), with the most buried dust
consistent with the molecular gas and HII regions. This suggests the dust
exists throughout star-forming regions where a significant fraction of the
flux is from the outer layers, where the dust emission is unobscured towards
the line of sight of the observer. A significant portion of the heating of
this unobscured dust may therefore originate from old stars that exist outwith
the obscured molecular clouds where stars are actively forming. This is
confirmed by the extinction from the stellar continuum, which at these
wavelengths is dominated by evolved stars, and appears less buried than the
more obscured molecular gas and HII regions. This is also true for the PAHs,
where previous studies have shown that some excitation can be attributed to
more evolved stellar populations (Zhang & Ho, 2023a), which is an important
consideration when determining star-formation rates.

PAHs experience similar obscuration as the cooler dust continuum. This is
however, largely by construction, as the intrinsic 12.7/11.3 ratio is
empirically calibrated using the average extinction of the continuum (see
Section 4.3.3). However, one would expect the PAHs to be somewhat less
obscured than the HII regions as they are known to exist in
Photodisassociation Regions (PDRs), where the PAHs exist on the boundary
between ionised and molecular gas (e.g. See Fig. 1 of Chown et al., 2023) and
therefore the measured values (in Table 1) are not unexpected.

We find that the molecular gas temperature is generally high, likely due to
these lines being sensitive to the warm molecular gas. If the S(0) transition
was also included, the measured average temperature would likely be lower. We
measure generally similar gas temperatures of $\sim 500-700$ K for all of the
regions with the exception of IIZw96 SW, which shows a higher temperature of
$\sim 900$ K. Interestingly, the known AGN, NGC 7469 Nuc and NGC 3256 S both
show temperatures consistent with the star-forming regions. This is likely
because the mid-IR transitions of H2 originate from gas at temperatures around
those shown in Table 1 and so we are sensitive to molecular gas heating by
stars rather than AGN. This has been found in previous studies (e.g.
Rigopoulou et al., 2002) where AGN and star-forming dominated galaxies show
similar molecular gas temperatures as traced by the mid-IR.

###  5.2 Choice of Emissivity

We found a more physically motivated emissivity, which includes silicates,
provides a better fit to the data over a simple power law. In Fig. 9 we show a
comparison of the inferred dust distribution for a star-forming region
depending on the choice of emissivity.

Figure 9: Comparison of the inferred dust distribution for a typical star-
forming region (NGC 3256 SF3) depending on the chosen emissivity function as
shown in Fig. 2. The left is for $\epsilon_{\nu}\propto\lambda^{-2}$ while the
right is a physically motivated emissivity function, including silicates,
based on Li & Draine (2001).

For the power law emissivity, the dust is almost completely unobscured and
shows no sign of any differential extinction, as to create a smooth continuum
with little silicate absorption requires close to zero extinction. By contrast
the dust distribution for the physically motivated emissivity shows a greater
average extinction compared to the power law emissivity, with evidence of
differential extinction, with a more obscured hot dust component relative to
the cooler dust.

The latter case is physically more reasonable, as the dust extinction is more
consistent with the other extinction tracers, where the buried hot dust is
consistent with the obscured molecular gas where stars are actively forming.
We therefore conclude that the physically motivated emissivity is closer to
reality.

###  5.3 Decomposing the Dust Distribution

For objects that show a complex dust distribution such as IIZw96 SW, we can
decompose this into multiple components. For this example, there is a clear
component of relatively unobscured cold dust that is consistent with the HI,
H2 and stellar continuum extinction measures. Additionally there is an
obscured component with a gradient, where the hotter dust appears more
obscured. This is potentially an AGN torus (García-Bernete et al., 2024b) and
so if we split the dust distribution into two, we can extract individual
continua for each component.

As we find the majority of the dust distribution for star-forming regions lie
at $\tau_{9.8}<1\sim 1.5$, we define a cutoff around these values for NGC 7469
Nuc, VV114 NE, VV114 SW, IIZw96 SW and NGC 3256 S, as these show complex dust
distributions with clear components. Firstly, we select dust from star-
formation by using grid elements below $\tau_{9.8}<1.5$ for NGC 3256 S, VV114
NE and IIZw96 SW and $\tau_{9.8}<1.0$ for VV114 SW. We then note that NGC 7469
Nuc, NGC 3256 S and VV114 SW show a clear separate component in the
$T,\tau_{9.8}$ space at high temperatures. We isolate this with a temperature
cutoff of $T>400$ K. We also apply this to IIZw96 SW, which tentatively shows
a separate component here.

The obscured, high temperature component is likely heated very close to the
AGN, near the dust sublimation zone. This may be the inner rim of a torus or
polar dust, where dusty grains are driven as a wind by radiation pressure from
the AGN accretion disk (e.g. Hönig & Kishimoto, 2017; García-Bernete et al.,
2022d).

We have presented a simple interpretation of the best fit dust distribution,
but the reality may be more complex. For example if there is dust at a similar
temperature and obscuration in both an AGN torus and star-forming regions, the
two cannot be disentangled by a simple cutoff.

###  5.4 Nature of the Individual Sources

From the best fit dust distributions and the extinction predicted from the
different emission features, we can infer the nature of these sources in terms
of layers of extinction. We briefly summarise our findings for each object
here in the context of the literature. Beyond the scope of this work, the rich
data of these targets allows much more analysis into the nature of these
objects (e.g. Inami et al., 2022; García-Bernete et al., 2022c; Donnan et al.,
2023a; Rich et al., 2023; Armus et al., 2023; Lai et al., 2023; Zhang & Ho,
2023b; Pereira-Santaella et al., 2023; González-Alfonso et al., 2023; García-
Bernete et al., 2024b; Buiten et al., 2023).

####  5.4.1 NGC 7469

We find the nuclear dust continuum to be dominated by the AGN, with relatively
weak PAH emission and no significant contribution from stellar continuum,
consistent with previous analysis of this target (e.g. Alonso-Herrero et al.,
2020; García-Bernete et al., 2022a; Armus et al., 2023). The dust distribution
shows a large spread of extinction with an average extinction of
$\tau_{9.8}\sim 1$ with a significant hot dust component with $T>400$ K.
Considering its classification as a type 1 AGN, it is expected to have clear
lines of sight to very hot dust close to the sublimation radius at T $>1000$
K. Interestingly, this dust is moderately obscured which may suggest that
there is some intervening dust that obscures the very hot inner regions of the
AGN torus.

For three of the star-forming regions (see Section 4.2 and Fig. 6) we find the
HII regions to be more obscured than the molecular gas, contrary to the other
SF regions. Namely SF2, SF3 and SF4 show this behaviour and lie on the north-
west section of the star-forming ring. Additionally, these regions show
silicates in emission (along with SF1) due to a combination of low obscuration
and higher dust temperatures. From Díaz-Santos et al. (2007), these regions
show slightly more evolved stellar populations and lack molecular gas (via CO
(2-1)) compared to SF1, SF5 and SF6. This offers a potential explanation as to
why the H2 lines appear less obscured than the HII regions, where the relative
lack of CO but presence of the mid-IR H2 lines would suggest a high molecular
gas temperature.

For NGC 7469 SF1, the comparative lack of cool obscuring dust, resulting in
the silicate emission, may be related to the nuclear bar, as this region is
placed on the boundary between the bar and the star-forming ring.

####  5.4.2 VV 114

VV 114 NE shows a significant highly obscured dust continuum appearing more
obscured than the PAHs, molecular gas or HII regions. This is consistent with
initial analysis of this target in Donnan et al. (2023a) where there is likely
an extremely buried AGN considering its obscuration and compact nature (<100
pc). Additionally, the ro-vibration molecular bands of CO and H2O, show a high
excitation temperature (Buiten et al., 2023). From the dust distribution, high
temperatures are reached up to $\sim$ 800 K but unlike the SW nucleus, there
is no clear isolated hot dust component.

VV 114 SW also shows a strong continuum but it is significantly less obscured
than the NE nucleus. There is a significant isolated hot dust component that
dominates the continuum between 2-10 $\mu$m. This is consistent with the claim
that this is an AGN by Rich et al. (2023), due to its low 3.3 $\mu$m and 6.2
$\mu$m PAH equivalent widths, where the hot dust continuum dominates over the
emission features. From the decomposition, the warm AGN component shows a
gradient where the obscuration increases with temperature. This may indicate a
high covering factor smooth torus. This picture can explain the lack of high
ionisation potential lines (Donnan et al., 2023a; Rich et al., 2023) that
would be expected for an AGN with the observed level of obscuration.

An alternative explanation is that the isolated hot dust component is a result
of the shock front that has passed through the region (e.g. Saito et al.,
2017; Donnan et al., 2023a) and that the nucleus is simply an obscured star-
forming region. This may explain the lack of high-ionisation lines, and the
low excitation temperature of the ro-vibration molecular bands (Buiten et al.,
2023). However, similar to the NE core, the source is very compact (< 100 pc
in the near-IR) and so reaching such a high infrared surface brightness is
difficult without an AGN.

####  5.4.3 NGC 3256

NGC 3256S hosts a low luminosity, Compton Thick AGN (e.g. Ohyama et al., 2015)
which has a collimated molecular outflow (Sakamoto et al., 2014; Emonts et
al., 2014; Pereira-Santaella et al., 2023) detected through CO, HCN, HCO+ and
H2 transitions. From the dust distribution we see a clear signature of a
buried continuum source which shows a gradient of extinction, with the
extinction increasing with temperature in the warm component. This suggests a
smooth torus rather than a clumpy one, as the hot dust is buried and not
leaking through low density lines of sight. In addition, the dust distribution
shows an isolated hot dust component that is also highly obscured. This hot
dust, as shown in Fig. 8, dominates the continuum between 2 - 10 $\mu$m.

The stellar continuum, PAHs and HI extinction show a similar pattern to the
star-forming regions and VV114 NE, where these trace different phases of the
circumnuclear star-formation. However, the H2 appears less obscured as it is
likely tracing the molecular outflow dominating the H2 emission over
circumnuclear star-forming molecular clouds.

The northern nucleus shows a dust distribution that is hotter with higher
extinction than the star-forming regions but does not show any significant hot
component that would be consistent with AGN heating. This is consistent with
the literature, that this is a star-forming nucleus (Lira et al., 2008).

Interestingly, the star-forming region to the southwest, SF5, shows a higher
obscuration and dust temperature unlike the other star-forming regions as well
as deeper near-IR ices. This may suggest an additional heating source such as
supernovae and/or X-ray binaries. This is not unexpected considering LIRGs are
believed to have high supernovae rates (many per year) owing to their high
star-formation rates. Multiple supernovae have been detected in NGC 3256 (e.g.
Kankare et al., 2018), however none match this specific source.

With the exception of SF5, all the star-forming regions in NGC 3256 are
remarkably similar, with all showing the same pattern of obscuration where the
stellar continuum, PAHs, HII regions and molecular gas exist from least
obscured to most obscured respectively. The star-forming regions in NGC 3256
appear slightly cooler than those in NGC 7469 and therefore do not show any
silicate emission.

A detailed analysis of the PAH properties, comparing to theoretical spectra,
will provide further insight (Rigopoulou et al. in prep.).

####  5.4.4 IIZw96

The spectrum of the SW nuclei of IIZw96 is highly unusual (Inami et al., 2022;
García-Bernete et al., 2024b), with very weak PAH features and a plethora of
emission features in the near-infrared (< 3 $\mu$m). Indeed, Shipley et al.
(2016) noted this object as an outlier in the PAH vs SFR correlation, due to
its weak PAH emission as probed by Spitzer.

From the dust distribution we find evidence for a cold, buried continuum
source, which shows a gradient, which would again suggest a smooth obscuring
structure. The star-forming component is strong, which dilutes this obscured
source, filling up the silicate feature. Additionally, the star-forming
component is hot, leading to silicate emission at $9.8$ $\mu$m, consistent
with the higher molecular gas temperature of $\sim 900$ K.

The $\sim 4.6\mu$m CO band appears in emission rather than absorption (García-
Bernete et al., 2024b). This is the case for NGC 3256 S, where the nucleus
shows CO absorption while the outflow shows emission (Pereira-Santaella et
al., 2023). We tentatively detect a very hot dust component in the dust
distribution and thus contributing to the continuum at $\sim$ 3 $\mu$m. Unlike
NGC 3256 S and VV114 SW, this is not as clearly isolated from the warm
component and so the temperature/extinction is not well constrained, however
this component dominates the continuum in the near-IR band. The values for the
extinction and temperature for this hot component in Table 1 are consequently
very uncertain. In particular the extinction, as the other dust components
dominate in the silicate bands.

The NE nuclei appears to be a relatively unobscured star-forming region,
however the PAH features appear a-typical. Most notably is the extremely
strong 8.6 $\mu$m feature. Additionally, the molecular gas appears less
obscured than the HII regions and PAHs.

###  5.5 Effective Extinction Curve

For each object, the extinction correction factor, defined as the ratio of the
observed to the intrinsic spectrum (equation (11)), varies depending on the
shape of the best fit dust distribution, $\Psi\left(T,\tau_{9.8}\right)$. In
general this dust distribution is not a screen and therefore cannot be
converted into an optical depth as $\frac{f_{\nu}}{f_{\nu}^{\textrm{Int}}}\neq
e^{-\tau_{9.8}\tau(\lambda)}$. However, to compare the inferred dust
correction factors with known extinction curves from the literature, we
calculate an “effective” extinction curve by inverting this equation. This can
be thought as “locally” assuming a screen at each wavelength and stitching
together an extinction curve. It is worth noting that we are not deriving a
new extinction law, as in general $\frac{f_{\nu}}{f_{\nu}^{\textrm{Int}}}\neq
e^{-\tau_{9.8}\tau(\lambda)}$ and therefore the variation we see between
targets is purely due to the shape of $\Psi\left(T,\tau_{9.8}\right)$. In this
work we have assumed a fixed law of that given in Donnan et al. (2023a), which
governs the intrinsic properties of the dust grains.

Figure 10: Comparison of the inferred extinction factor (eqaution (11)) for
each target, converted to an “effective” optical depth, normalised at 9.8
$\mu$m. We show the targets in Fig. 1 as well as the average for the star-
forming regions ing NGC 3256 and NGC 7469. We compare to the extinction law
used in this work from Donnan et al. (2023a) and others from the literature,
namely Ossenkopf et al. (1992); Kemper et al. (2004); Fritz et al. (2011);
Jones et al. (2013).

Fig. 10 shows the “effective” optical depth normalised at 9.8 $\mu$m for a
variety of our targets.

All of the targets show a shallower gradient in the near-infrared compared to
a single screen from any of the extinction laws. This highlights the effect of
differential extinction, where the mid-infrared traces more obscured activity
than the near-infrared. This is particularly true for VV114 NE, where we
observe a turnover at $\lambda<4\mu$m, where the relatively less obscured
stellar continuum dominates over the highly buried dust continuum. IIZw96 also
shows an unusual shape, where the silicate feature is “filled up” by the hot
dust emission and thus the shorter wavelength portion appears at a relatively
higher optical depth compared tro the silicate band.

Extinction curves with flatter optical depths at $\sim$ 5 $\mu$m, have been
observed in the literature (see Fig. 4 of Draine (2003) or Fig. 2 of Li et al.
(2015)), deviating from the expected power law. Differential extinction offers
a potential explanation, where the mid-infrared traces more obscured dust
while the near-infrared traces less obscured stellar-continuum, flattening the
“effective” extinction curve. When the near-infrared is dominated by buried
hot dust instead, such as for NGC 3256 S, the “effective” extinction curve
becomes steeper, better resembling a power law.

The ratio of the 9.8 $\mu$m and 18 $\mu$m silicate bands also varies between
targets, depending on the relative obscuration and temperature of different
dust components. This can relate to the clumpiness of AGN tori (e.g.
Hatziminaoglou et al., 2015; Martínez-Paredes et al., 2020; García-Bernete et
al., 2024a).

##  6 Conclusions

We have presented a model to reproduce the joint NIRSpec and MIRI spectra
between 1.5 - 28 $\mu$m for a variety of LIRGs. Our main findings are

  * •

By testing with simulated data, we show that the model is able to reproduce
the dust distribution and thus extinction correction factor for input data
with and without differential extinction, where in the latter case, the
traditional screen/mixed approaches fail, highlighting the need for
differential extinction to model the spectra of obscured (U)LIRGs.

  * •

We find evidence for differential extinction from fitting the spectra of LIRGs
where star-forming regions show dust distribution dominated by relatively
unobscured cool dust with some buried hot dust.

  * •

We find a deeply buried continuum component in the dust distribution for VV114
NE, NGC 3256 S and IIZw96 SW, with the latter two showing a hot isolated dust
component at T$\sim 1000$ K, likely due to AGN heating. Similarly, VV114 SW
shows a significant hot component, either due to AGN heating or a shock front
passing through as a result of the galaxy-galaxy merger.

  * •

We compare the dust extinction with that of the stellar continuum, PAHs, HII
regions (HI lines) and molecular gas (H2) lines. For star-forming regions we
find that the molecular gas, and HII regions are the most obscured consistent
with the hottest dust. The stellar continuum and PAHs are significantly less
obscured, consistent with the cooler dust.

  * •

By constructing an “effective” extinction curve, we find the near-infrared to
be much shallower than the mid-infrared which is indicative of differential
extinction, where the near-infrared traces less obscured emission than the
mid-infrared. This can explain flatter extinction curves in the literature.

We have shown that by modelling the distribution of dust we are able to
reproduce the spectrum for a variety of objects from star-forming regions to
highly obscured AGN. Such detailed modelling is imperative considering the
high frequency of dusty galaxies being observed with JWST.

## Acknowledgements

The authors are grateful to the DD-ERS team (Program 1328, PI: Lee Armus &
Aaron Evans) for the observing program with a zero–exclusive–access period.
FRD acknowledges support from STFC through grant ST/W507726/1. DR and IGB
acknowledge support from STFC through grant ST/S000488/1 and ST/W000903/1. AAH
acknowledges support from grant PID2021-124665NB-I00 funded by
MCIN/AEI/10.13039/501100011033 and by ERDF A way of making Europe.

## Data Availability

All the data used in this work is publicly available as part of DD-ERS Program
1328, downloadable from the MAST archive.

## References

  * Aalto et al. (2015) Aalto S., et al., 2015, A&A, 584, A42
  * Aalto et al. (2019) Aalto S., et al., 2019, A&A, 627, A147
  * Alonso-Herrero et al. (2011) Alonso-Herrero A., et al., 2011, ApJ, 736, 82
  * Alonso-Herrero et al. (2012) Alonso-Herrero A., Pereira-Santaella M., Rieke G. H., Rigopoulou D., 2012, ApJ, 744, 2
  * Alonso-Herrero et al. (2020) Alonso-Herrero A., et al., 2020, A&A, 639, A43
  * Argyriou et al. (2023) Argyriou I., et al., 2023, A&A, 675, A111
  * Armus et al. (2009) Armus L., et al., 2009, PASP, 121, 559
  * Armus et al. (2023) Armus L., et al., 2023, ApJ, 942, L37
  * Bianchin et al. (2023) Bianchin M., et al., 2023, arXiv e-prints, p. arXiv:2308.00209
  * Bohlin & Lockwood (2022) Bohlin R. C., Lockwood S., 2022, Update of the STIS CTE Correction Formula for Stellar Spectra, Instrument Science Report STIS 2022-7, 11 pages
  * Bohlin et al. (2014) Bohlin R. C., Gordon K. D., Tremblay P. E., 2014, PASP, 126, 711
  * Böker et al. (2023) Böker T., et al., 2023, PASP, 135, 038001
  * Brandl et al. (2006) Brandl B. R., et al., 2006, ApJ, 653, 1129
  * Buiten et al. (2023) Buiten V. A., et al., 2023, arXiv e-prints, p. arXiv:2312.01945
  * Cappellari (2017) Cappellari M., 2017, MNRAS, 466, 798
  * Cappellari (2023) Cappellari M., 2023, MNRAS, 526, 3273
  * Chown et al. (2023) Chown R., et al., 2023, arXiv e-prints, p. arXiv:2308.16733
  * Conroy & Gunn (2010) Conroy C., Gunn J. E., 2010, FSPS: Flexible Stellar Population Synthesis, Astrophysics Source Code Library, record ascl:1010.043 (ascl:1010.043)
  * Conroy et al. (2009) Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
  * Díaz-Santos et al. (2007) Díaz-Santos T., Alonso-Herrero A., Colina L., Packham C., Radomski J. T., Telesco C. M., 2007, in Guzmán R., ed., Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 29, Revista Mexicana de Astronomia y Astrofisica Conference Series. pp 92–94 (arXiv:astro-ph/0610616), doi:10.48550/arXiv.astro-ph/0610616
  * Donnan et al. (2023a) Donnan F. R., García-Bernete I., Rigopoulou D., Pereira-Santaella M., Alonso-Herrero A., Roche P. F., Hernán-Caballero A., Spoon H. W. W., 2023a, MNRAS, 519, 3691
  * Donnan et al. (2023b) Donnan F. R., et al., 2023b, A&A, 669, A87
  * Draine (1989) Draine B. T., 1989, in Böhm-Vitense E., ed., Infrared Spectroscopy in Astronomy. p. 93
  * Draine (2003) Draine B. T., 2003, ARA&A, 41, 241
  * Draine & Li (2007) Draine B. T., Li A., 2007, ApJ, 657, 810
  * Efstathiou et al. (2013) Efstathiou A., Christopher N., Verma A., Siebenmorgen R., 2013, MNRAS, 436, 1873
  * Emonts et al. (2014) Emonts B. H. C., Piqueras-López J., Colina L., Arribas S., Villar-Martín M., Pereira-Santaella M., Garcia-Burillo S., Alonso-Herrero A., 2014, A&A, 572, A40
  * Falstad et al. (2021) Falstad N., et al., 2021, A&A, 649, A105
  * Fritz et al. (2006) Fritz J., Franceschini A., Hatziminaoglou E., 2006, MNRAS, 366, 767
  * Fritz et al. (2011) Fritz T. K., et al., 2011, ApJ, 737, 73
  * Gallimore et al. (2010) Gallimore J. F., et al., 2010, ApJS, 187, 172
  * García-Bernete et al. (2017) García-Bernete I., Ramos Almeida C., Landt H., Ward M. J., Baloković M., Acosta-Pulido J. A., 2017, MNRAS, 469, 110
  * García-Bernete et al. (2019) García-Bernete I., et al., 2019, MNRAS, 486, 4917
  * García-Bernete et al. (2022a) García-Bernete I., Rigopoulou D., Alonso-Herrero A., Pereira-Santaella M., Roche P. F., Kerkeni B., 2022a, MNRAS, 509, 4256
  * García-Bernete et al. (2022b) García-Bernete I., Rigopoulou D., Aalto S., Spoon H. W. W., Hernán-Caballero A., Efstathiou A., Roche P. F., König S., 2022b, A&A, 663, A46
  * García-Bernete et al. (2022c) García-Bernete I., et al., 2022c, A&A, 666, L5
  * García-Bernete et al. (2022d) García-Bernete I., et al., 2022d, A&A, 667, A140
  * García-Bernete et al. (2024a) García-Bernete I., et al., 2024a, A&A, 681, L7
  * García-Bernete et al. (2024b) García-Bernete I., Pereira-Santaella M., González-Alfonso E., Rigopoulou D., Efstathiou A., Donnan F. R., Thatte N., 2024b, A&A, 682, L5
  * González-Alfonso et al. (2023) González-Alfonso E., García-Bernete I., Neufeld D., Pereira-Santaella M., Fischer J., Donnan F. R., 2023, A&A, submitted
  * González-Martín et al. (2019a) González-Martín O., et al., 2019a, ApJ, 884, 10
  * González-Martín et al. (2019b) González-Martín O., et al., 2019b, ApJ, 884, 11
  * Gordon et al. (2021) Gordon K. D., et al., 2021, ApJ, 916, 33
  * Hatziminaoglou et al. (2015) Hatziminaoglou E., Hernán-Caballero A., Feltre A., Piñol Ferrer N., 2015, ApJ, 803, 110
  * Hernán-Caballero et al. (2016) Hernán-Caballero A., Hatziminaoglou E., Alonso-Herrero A., Mateos S., 2016, MNRAS, 463, 2064
  * Hernán-Caballero et al. (2020) Hernán-Caballero A., et al., 2020, MNRAS, 497, 4614
  * Hernandez et al. (2023) Hernandez S., et al., 2023, ApJ, 948, 124
  * Hönig & Kishimoto (2017) Hönig S. F., Kishimoto M., 2017, ApJ, 838, L20
  * Hönig et al. (2010) Hönig S. F., Kishimoto M., Gandhi P., Smette A., Asmus D., Duschl W., Polletta M., Weigelt G., 2010, A&A, 515, A23
  * Hummer & Storey (1987) Hummer D. G., Storey P. J., 1987, MNRAS, 224, 801
  * Imanishi et al. (2010) Imanishi M., Nakagawa T., Shirahata M., Ohyama Y., Onaka T., 2010, ApJ, 721, 1233
  * Inami et al. (2010) Inami H., et al., 2010, AJ, 140, 63
  * Inami et al. (2022) Inami H., et al., 2022, ApJ, 940, L6
  * Iyer et al. (2019) Iyer K. G., Gawiser E., Faber S. M., Ferguson H. C., Kartaltepe J., Koekemoer A. M., Pacifici C., Somerville R. S., 2019, ApJ, 879, 116
  * Jakobsen et al. (2022) Jakobsen P., et al., 2022, A&A, 661, A80
  * Jones et al. (2013) Jones A. P., Fanciullo L., Köhler M., Verstraete L., Guillet V., Bocchio M., Ysard N., 2013, A&A, 558, A62
  * Kankare et al. (2018) Kankare E., et al., 2018, The Astronomer’s Telegram, 11156, 1
  * Kemper et al. (2004) Kemper F., Vriend W. J., Tielens A. G. G. M., 2004, ApJ, 609, 826
  * Kishimoto et al. (2007) Kishimoto M., Hönig S. F., Beckert T., Weigelt G., 2007, A&A, 476, 713
  * Labiano et al. (2021) Labiano A., et al., 2021, A&A, 656, A57
  * Lai et al. (2020) Lai T. S. Y., Smith J. D. T., Baba S., Spoon H. W. W., Imanishi M., 2020, ApJ, 905, 55
  * Lai et al. (2022) Lai T. S. Y., et al., 2022, ApJ, 941, L36
  * Lai et al. (2023) Lai T. S. Y., et al., 2023, arXiv e-prints, p. arXiv:2307.15169
  * Lallement et al. (2022) Lallement R., Vergely J. L., Babusiaux C., Cox N. L. J., 2022, A&A, 661, A147
  * Le Floc’h et al. (2005) Le Floc’h E., et al., 2005, ApJ, 632, 169
  * Leja et al. (2019) Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019, ApJ, 876, 3
  * Li (2020) Li A., 2020, Nature Astronomy, 4, 339
  * Li & Draine (2001) Li A., Draine B. T., 2001, ApJ, 554, 778
  * Li et al. (2015) Li A., Wang S., Gao J., Jiang B. W., 2015, in , Lessons from the Local Group: A Conference in honor of David Block and Bruce Elmegreen. pp 85–103, doi:10.1007/978-3-319-10614-4_8
  * Lira et al. (2008) Lira P., Gonzalez-Corvalan V., Ward M., Hoyer S., 2008, MNRAS, 384, 316
  * Lonsdale et al. (2006) Lonsdale C. J., Farrah D., Smith H. E., 2006, in Mason J. W., ed., , Astrophysics Update 2. p. 285, doi:10.1007/3-540-30313-8_9
  * Magnelli et al. (2011) Magnelli B., Elbaz D., Chary R. R., Dickinson M., Le Borgne D., Frayer D. T., Willmer C. N. A., 2011, A&A, 528, A35
  * Marshall et al. (2007) Marshall J. A., Herter T. L., Armus L., Charmandaris V., Spoon H. W. W., Bernard-Salas J., Houck J. R., 2007, ApJ, 670, 129
  * Martínez-Paredes et al. (2015) Martínez-Paredes M., et al., 2015, MNRAS, 454, 3577
  * Martínez-Paredes et al. (2020) Martínez-Paredes M., et al., 2020, ApJ, 890, 152
  * Martínez-Paredes et al. (2021) Martínez-Paredes M., et al., 2021, ApJ, 922, 157
  * Nenkova et al. (2008a) Nenkova M., Sirocky M. M., Ivezić Ž., Elitzur M., 2008a, ApJ, 685, 147
  * Nenkova et al. (2008b) Nenkova M., Sirocky M. M., Nikutta R., Ivezić Ž., Elitzur M., 2008b, ApJ, 685, 160
  * Ohyama et al. (2015) Ohyama Y., Terashima Y., Sakamoto K., 2015, ApJ, 805, 162
  * Ossenkopf et al. (1992) Ossenkopf V., Henning T., Mathis J. S., 1992, A&A, 261, 567
  * Osterbrock & Martel (1993) Osterbrock D. E., Martel A., 1993, ApJ, 414, 552
  * Pereira-Santaella et al. (2014) Pereira-Santaella M., Spinoglio L., van der Werf P. P., Piqueras López J., 2014, A&A, 566, A49
  * Pereira-Santaella et al. (2015) Pereira-Santaella M., et al., 2015, A&A, 577, A78
  * Pereira-Santaella et al. (2021) Pereira-Santaella M., et al., 2021, A&A, 651, A42
  * Pereira-Santaella et al. (2022) Pereira-Santaella M., et al., 2022, arXiv e-prints, p. arXiv:2208.04835
  * Pereira-Santaella et al. (2023) Pereira-Santaella M., González-Alfonso E., García-Bernete I., García-Burillo S., Rigopoulou D., 2023, arXiv e-prints, p. arXiv:2309.06486
  * Ramos Almeida et al. (2009) Ramos Almeida C., et al., 2009, ApJ, 702, 1127
  * Rich et al. (2011) Rich J. A., Kewley L. J., Dopita M. A., 2011, ApJ, 734, 87
  * Rich et al. (2023) Rich J., et al., 2023, ApJ, 944, L50
  * Rieke et al. (2015) Rieke G. H., et al., 2015, PASP, 127, 584
  * Rigopoulou et al. (1999) Rigopoulou D., Spoon H. W. W., Genzel R., Lutz D., Moorwood A. F. M., Tran Q. D., 1999, AJ, 118, 2625
  * Rigopoulou et al. (2002) Rigopoulou D., Kunze D., Lutz D., Genzel R., Moorwood A. F. M., 2002, A&A, 389, 374
  * Rigopoulou et al. (2021) Rigopoulou D., et al., 2021, MNRAS, 504, 5287
  * Saito et al. (2017) Saito T., et al., 2017, ApJ, 834, 6
  * Sakamoto et al. (2014) Sakamoto K., Aalto S., Combes F., Evans A., Peck A., 2014, ApJ, 797, 90
  * Salim & Narayanan (2020) Salim S., Narayanan D., 2020, ARA&A, 58, 529
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Sanders & Mirabel (1996) Sanders D. B., Mirabel I. F., 1996, ARA&A, 34, 749
  * Shipley et al. (2016) Shipley H. V., Papovich C., Rieke G. H., Brown M. J. I., Moustakas J., 2016, ApJ, 818, 60
  * Smith et al. (2007) Smith J. D. T., et al., 2007, ApJ, 656, 770
  * Spoon et al. (2007) Spoon H. W. W., Marshall J. A., Houck J. R., Elitzur M., Hao L., Armus L., Brandl B. R., Charmandaris V., 2007, ApJ, 654, L49
  * Spoon et al. (2022) Spoon H. W. W., et al., 2022, arXiv e-prints, p. arXiv:2203.03071
  * Stalevski et al. (2016) Stalevski M., Ricci C., Ueda Y., Lira P., Fritz J., Baes M., 2016, MNRAS, 458, 2288
  * Stancik & Brauns (2008) Stancik A. L., Brauns E. B., 2008, Vibrational Spectroscopy, 47, 66
  * Stierwalt et al. (2013) Stierwalt S., et al., 2013, ApJS, 206, 1
  * Togi & Smith (2016) Togi A., Smith J. D. T., 2016, ApJ, 830, 18
  * U et al. (2022) U V., et al., 2022, ApJ, 940, L5
  * Veilleux et al. (2009) Veilleux S., et al., 2009, ApJS, 182, 628
  * Wells et al. (2015) Wells M., et al., 2015, PASP, 127, 646
  * Wright et al. (2015) Wright G. S., et al., 2015, PASP, 127, 595
  * Young et al. (2023) Young J., et al., 2023, ApJ, 958, L5
  * Zhang & Ho (2023a) Zhang L., Ho L. C., 2023a, ApJ, 943, 60
  * Zhang & Ho (2023b) Zhang L., Ho L. C., 2023b, ApJ, 953, L9

##  Appendix A Example Dust Geometries

In Fig 11, we show three diagrams demonstrating possible dust geometries we
have described in this work.

Figure 11: Diagrams demonstrating three different dust geometries discussed in
this work. The left panels show the dust distribution,
$\Psi\left(T,\tau_{9.8}\right)$ corresponding to the geometry in the right
panels. In each case the observer is on the right. The differential case is
just one example where differential extinction is present.

##  Appendix B PAH Features

We show in Table 2, the parameters governing the Drude profiles used to model
the PAH emission (equations (16) and (17)).

Table 2: PAH Feature Parameters of equations (16,17)

$\lambda_{0}$ | $\gamma_{0}$ | $a$
---|---|---
$\mu$m |  $\mu$m |
(1) | (2) | (3)
3.29 | 0.04 | 0.52
3.40 | 0.03 | -10.0
3.47 | 0.10 | -0.80
5.18 | 0.05 | 0.00
5.24 | 0.10 | -3.00
5.45 | 0.15 | 0.00
5.53 | 0.10 | 0.00
5.64 | 0.10 | 0.00
5.70 | 0.10 | 0.00
5.76 | 0.10 | 0.00
5.87 | 0.15 | 0.00
6.00 | 0.20 | 0.00
6.20 | 0.15 | -6.00
6.69 | 0.40 | 0.00
7.10 | 0.40 | 0.00
7.42 | 0.94 | 0.00
7.55 | 0.30 | 0.00
7.61 | 0.10 | 0.00
7.82 | 0.40 | 0.00
8.33 | 0.20 | 0.00
8.50 | 0.20 | 0.00
8.61 | 0.34 | 0.00
10.60 | 0.10 | 0.00
10.74 | 0.10 | 0.00
11.00 | 0.10 | -1.10
11.20 | 0.10 | 0.00
11.26 | 0.30 | -10.0
11.99 | 0.54 | 0.00
12.60 | 0.50 | 0.00
12.77 | 0.15 | 0.00
13.15 | 0.50 | 0.00
13.55 | 0.20 | -5.00
14.04 | 0.20 | 0.00
14.19 | 0.20 | -5.00
15.90 | 0.32 | 0.00
16.45 | 0.23 | 0.00
17.04 | 1.11 | 0.00
17.38 | 0.21 | 0.00

  * Column (1): Initial central wavelength for the fit with range $\pm 0.05\mu$m. Column (2): Initial FWHM allowed to vary within $+10\%$ and $-60\%$. Column (3): Initial asymmetry parameter $a$, allowed to vary within $+50\%$ and $-50\%$. Those with $a=0.00$ are fixed to a symmetric profile.

  *

##  Appendix C IIZw96 SW Contamination Correction

Fig. 12 shows the correction of the extracted spectrum of IIZw96 SW for
contamination, due to the nearby secondary nucleus. This process is described
in detail in Section 4.2.

Figure 12: Correction of IIZw96 SW spectrum for contamination from nearby
nucleus. The black line shows the total spectrum of the two nuclei from a 0.6"
aperture, while the green shows the sum of the individual nuclei which
displays an excess flux due to contamination. The blue spectrum shows the
individual nucleus of interest which once corrected for contamination, returns
the spectrum shown in red.

##  Appendix D Additional SF Regions

Fig. 13 shows the remaining fits and inferred dust distributions for the
spectra fitted in this work. These are the remaining star-forming regions from
NGC 3256 and NGC 7469.

Figure 13: Same as Fig. 8 but for the remaining spectra, all of which are
star-forming regions.