# The Obscured Nucleus and Shocked Environment of VV 114E Revealed by
JWST/MIRI Spectroscopy

F. R. Donnan, 1 I. García-Bernete,1 D. Rigopoulou,1 M. Pereira-Santaella,2 A.
Alonso-Herrero,3 P. F. Roche,1 A. Hernán-Caballero, 4 H. W. W. Spoon5
1Department of Physics, University of Oxford, Keble Road, Oxford, OX1 3RH, UK
2 Observatorio Astronómico Nacional (OAN-IGN)-Observatorio de Madrid, Alfonso
XII, 3, 28014, Madrid, Spain
3 Centro de Astrobiología (CAB), CSIC-INTA, Camino Bajo del Castillo s/n,
E-28692 Villanueva de la Cañada, Madrid, Spain
4 Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan,
1, E-44001 Teruel, Spain
5 Cornell Center for Astrophysics and Planetary Science, Ithaca, NY 14853, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Compact Obscured Nuclei (CONs) potentially hide extreme supermassive black
hole (SMBH) growth behind large column densities of gas/dust. We present a
spectroscopic analysis of the heavily obscured nucleus and the surrounding
environment of the eastern region of the nearby ($z=0.02007$) interacting
galaxy VV 114 with the JWST Mid-InfraRed Instrument (MIRI). We model the
spectrum from 4.9 - 28 $\mu$m to extract Polycyclic Aromatic Hydrocarbon (PAH)
emission and the underlying obscured continuum. We find that the NE nucleus
(A) is highly obscured where the low PAH equivalent width (EW) ratio,
EW(12.7)/EW(11.3), reveals a dust enshrouded continuum source. This is
confirmed by decomposing the continuum into nuclear and star-forming where the
nuclear component is found to be typical of CONs. The 11.3/6.2 PAH flux ratio
is consistent with originating in star-forming regions rather than typical
AGN. The second nucleus (B) is much less obscured, with PAH flux ratios also
typical of star-forming regions. We do not detect any high ionisation lines
such as [Ne V] or [Ne VI] which suggests that if an AGN is present it must be
highly obscured. Additionally, we detect a shock front south of the secondary
nucleus (B) in the [Fe II] (5.34 $\mu$m) line and in warm molecular hydrogen.
The 6.2 PAH emission does not spatially coincide with the low-J transitions of
H2 but rather appears strong at the shock front which may suggest destruction
of the ionised PAHs in the post-shock gas behind the shock front.

###### keywords:

galaxies: nuclei – galaxies: evolution – shock waves – techniques:
spectroscopic

††pubyear: 2015††pagerange: The Obscured Nucleus and Shocked Environment of VV
114E Revealed by JWST/MIRI Spectroscopy–C

##  1 Introduction

Luminous Infrared Galaxies (LIRGs, LIR $>10^{11}\rm{L}_{\odot}$) and
Ultraluminous Infrared Galaxies (ULIRGs, LIR $>10^{12}\rm{L}_{\odot}$) contain
some of the most extreme dusty environments and disturbed morphology, with a
high prevalence of major/minor galaxy mergers (e.g. Rigopoulou et al., 1999).
These obscured environments hide some of the most intense and rapid galaxy
evolution with high star-formation rates and supermassive black hole (SMBH)
growth. Studying such environments is a challenge at optical wavelengths due
to the high levels of dust extinction and therefore requires study at longer
wavelengths. The JWST Mid-InfraRed Instrument (MIRI) enables the study of
these environments in the mid-infrared at a significantly improved
spatial/spectral resolution and sensitivity compared to previous facilities
(Rieke et al., 2015; Wright et al., 2015; Wells et al., 2015).

The most obscured galactic nuclei, so-called Compact Obscured Nuclei (CONs)
show extreme column densities ($N_{\rm{H}}>10^{25}\rm{cm}^{-2},A_{v}>1000$
mag) (e.g. Aalto et al., 2019; Falstad et al., 2021, and references therein.)
where either a compact starburst or an accreting supermassive black hole
(SMBH) is completely obscured. These objects are prevalent among LIRGs and
ULIRGs in the local universe (Donnan et al., 2022; Falstad et al., 2021;
García-Bernete et al., 2022c) and therefore represent a crucial phase of
galaxy evolution where rapid SMBH growth may occur, especially considering the
increased prevalence of dusty galaxies towards cosmic noon (e.g. Magnelli et
al., 2011).

Detecting such objects is a challenge due to their obscured nature with most
studies being confined to the sub-mm regime (e.g. Aalto et al., 2015, 2019;
Falstad et al., 2021). However recent work by García-Bernete et al. (2022c) in
the mid-infrared has found that these objects can be identified by measuring
the impact of obscuration on the mid-IR continuum compared to Polycyclic
Aromatic Hydrocarbon (PAH) emission which traces relatively unobscured star-
formation in the circumnuclear region. The development of this technique is
particularly timely as the recent launch of the James Webb Space Telescope
(JWST) provides data in the mid-infrared with a significant increase in the
spatial and spectral resolution over previous facilities.

To identify CONs using PAHs involves measuring the impact of the 9.8 $\mu$m
silicate absorption feature on the continuum compared to the PAH emission,
through the use of the ratio of the Equivalent Widths (EWs) of different PAH
features. As the PAH emission originates from the relatively unobscured
circumnuclear region, the presence of a highly obscured nucleus will cause
deep silicate absorption at 9.8 $\mu$m and lower the continuum around the 11.3
$\mu$m PAH (García-Bernete et al., 2022c). This will increase its equivalent
width relative to the other PAH features where the continuum is less affected
by the extinction from the nucleus. Therefore the presence of a highly
obscured nucleus can be revealed through lower values of the EW(12.7)/EW(11.3)
and the EW(6.2)/EW(11.3) PAH ratios (García-Bernete et al., 2022c). Typical
star-forming galaxies that do not host a highly obscured nucleus show a
typically higher, approximately constant PAH EW ratio as the intrinsic flux
ratio is approximately constant and the extinction of the PAHs and continuum
cancel when taking the ratio (Hernán-Caballero et al., 2020; García-Bernete et
al., 2022c).

In this work we apply the technique to JWST data for the first time with MIRI
observations of the highly obscured interacting system VV 114 (IC 1623, Arp
236). This is a local ($z=0.02007$, $D=$ 84 Mpc, Strauss et al., 1992) galaxy
with a highly disturbed morphology. This object is a LIRG with $L_{\rm
IR}=10^{11.71}$ L⊙ (Armus et al., 2009) and is in the mid-stages of a major
merger (Stierwalt et al., 2013). The western region of this galaxy is much
more prominent in the UV/optical (e.g. Frayer et al., 1999; Evans et al.,
2022) whereas the eastern region is extremely dust rich, appearing very red in
the optical. Moving towards the infrared, the luminosity of the eastern region
dramatically increases (e.g. Charmandaris et al., 2004) and dominates over the
western region. Therefore the eastern region is ideal for testing the
diagnostic power of JWST/MIRI for extremely obscured environments.

The eastern region contains complex morphology with clumps of star formation,
shocks (Rich et al., 2011) and a potential buried AGN (Saito et al., 2015).
This nucleus is a very good CON candidate considering the heavy obscuration
and potential AGN power source. From Spitzer spectral mapping observations
Donnan et al. (2022) found that the PAH EW ratios approached those that might
indicate a CON as the aperture size was reduced towards the nucleus, however
the low angular resolution meant that the spectrum was heavily diluted by the
circumnuclear star-formation. The low angular resolution of Spitzer (PSF FWHM
of $\sim$ 3.8” at 14 $\mu$m which corresponds to an angular size of $\sim$ 1.6
kpc) meant that the eastern region of VV 114 was completely unresolved in the
spectral mapping mode with no structure visible. The almost factor of $\sim$10
increase in the spatial resolution of JWST/MIRI MRS (PSF FWHM of $\sim$ 0.48”
at 14 $\mu$m which corresponds to an angular size of $\sim$ 200 pc) allows
spatially resolved spectroscopic analysis of the eastern region of VV 114,
where the previously unresolved region now shows multiple “nuclei” and a very
complex morphology.

The eastern region of VV 114 has been observed with JWST as part of Early
Release Science Program 1328 (PI: Lee Armus) with MIRI imaging (Evans et al.,
2022) and MIRI MRS. We use this data to demonstrate the diagnostic power of
mid-infrared spectroscopy to reveal highly obscured nuclei.

In Section 2 we describe the data and reduction steps. In Section 3 we outline
our modified version of PAHFIT to fit JWST spectra. By applying this fitting
tool we measure the obscuration, properties of the various regions through
their PAH features. By generating line maps we measure the fine structure
lines and molecular hydrogen (H2) transitions that are tracing a shock in
Section 4. We then discuss a possible evolutionary scenario in Section 5.

##  2 Observations

To illustrate the region of VV 114E we show imaging from MIRI (Evans et al.,
2022) with a colour composite image in the left panel of Fig. 1. We used the
reduced data from the MAST
archive111https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html to
construct this image. There are however some issues with this reduction where
the dithered observations within a single filter have not been properly
aligned, which can be seen by the diffraction spikes splitting into three. As
we are using the image for illustration only for purposes, this is not a major
issue.

Two main point sources are apparent in this image, labelled A and B. Source A
is significantly redder than B and is the suggested buried AGN from ALMA
observations (Saito et al., 2015) where the continuum remains unresolved in
the millimetre. Source B shows a stronger continuum in the F560W band relative
to F1500W compared to source A. There are numerous other small clumps
appearing within the image however, most of these are beyond the field of view
of the MIRI-MRS observations. The F770W filter, centered on the 7.7 $\mu$m PAH
shows much more extended emission and traces areas of star-formation where
young stars radiatively excite PAH molecules.

###  2.1 MIRI MRS Reduction

MIRI MRS data consists of 4 sub-cubes with channel 1 (4.9–7.65 µm), channel 2
(7.51–11.71 µm), channel 3 (11.55–18.02 µm) and channel 4 (17.71–28.1 µm) with
an increasing field of view and pixel size towards longer wavelengths.

We downloaded the MIRI MRS level 2-b reduced data from the MAST archive and
constructed the 4 sub-cubes using the Spec3 pipeline (version 1.6.2). Before
running the Spec3 pipeline we performed the residual fringe subtraction step
separately to remove fringing still present after the fringe correction step
of the pipeline. Using the jwst$\\_$0913.pmap context of the CRDS we ran the
Spec3 pipeline to produce 4 sub-cubes and their corresponding backgrounds. We
estimated the background for each wavelength channel by calculating the median
in each channel of the background cubes before subtracting, following the same
method as in García-Bernete et al. (2022a).

To ensure the extracted spectra are from the same spatial position, we first
correct for small discrepancies in the WCS and spatially align all the sub-
channels. We use the DAOStarFinder algorithm (Stetson, 1987) from the
photutils python package to find the position of point sources in the
integrated intensity (moment 0) maps which trace the continuum of each of the
4 sub-channels. We then adjust the WCS header to align the measured point
sources. Subsequently, we extract the spectrum using a circular aperture with
a 1" diameter centered on the measured position. We do this for the two point
sources in the image (A and B) and apply the correction factor from Pereira-
Santaella et al. (2022) for this aperture (see also García-Bernete et al.,
2022a). We also extract the spectrum from a region where no source is present
to represent the extended emission (point C). We chose an aperture placed away
from the two point sources so as to not be contaminated by the PSF at the
longest wavelengths but close enough such that it is within the field-of-view
in the shortest wavelength channel. As the emission within aperture C is not a
point source we do not apply the aperture correction factor. The three
apertures are shown in the right panel of Fig. 1.

The three spectra are shown in Fig. 2 where the extended emission (source C)
shows the strongest emission features compared to source A and B which both
show a strong continuum. Source A is clearly the more obscured one showing
very deep silicate absorption at 9.8 $\mu$m and 18 $\mu$m. There is also
evidence for absorption due to water ice and aliphatic CH around 6 $\mu$m and
7.5 $\mu$m respectively. The spectrum of source B shows a strong continuum
around 5 - 6 $\mu$m, consistent with the colours seen in the imaging in Fig. 1
and Evans et al. (2022).

Figure 1:  Left: Color composite image of VV 114E using JWST/MIRI. The F560W
(blue) and F1500W (red) trace mostly continuum while the F770W (green) traces
the 7.7 $\mu$m PAH feature. The two main point sources are labelled A and B.
Right: Integrated intensity (moment 0) map of Channel 1 (4.90 $\mu$m – 7.65
$\mu$m) which is dominated by thermal continuum dust emission and the 6.2
$\mu$m PAH which is more spatially extended. The black circles show the
apertures used to extract the spectra shown in Fig. 2, labelled A, B and C.
Figure 2: MIRI/MRS spectra of three regions in VV 114E. Spectra A and B are
extracted as point sources using a 1" circular aperture with a PSF wavelength
dependant correction factor applied. Source C represents extended emission
without the presence of any point sources. The fluxes are normalised to 7.7
$\mu$m. The most prominent PAH features are labelled with the vertical dashed
lines.

##  3 Methods

###  3.1 Spectral Fitting

By modelling the entire spectrum we can obtain emission line/PAH fluxes and
extinction features using all of the available MIR data to constrain the
various components. A widely used tool from the Spitzer era was PAHFIT (Smith
et al., 2007) which decomposes the spectrum into various emission lines, PAH
features, a dust continuum and corresponding extinction. We use our own
modified version of PAHFIT which has been fine tuned to model JWST spectra
with additional PAH features, emission lines and extinction features. The
following sections describe the various components of the model.

####  3.1.1 Emission Lines

The mid-infrared wavelength regime is host to various emission lines including
molecular hydrogen (H2) lines and ionised fine structure lines in addition to
the PAH emission.

Unlike the original PAHFIT which modelled each emission line with a single
Gaussian, we model each line with two components, a Gaussian and a Lorentzian.
This allows asymmetries to be captured by the model and allows the shape to
vary between being Gaussian dominated (narrow wings) or Lorentzian dominated
(broad wings). Each line is therefore modelled as the sum of a Gaussian,
$G\left(\lambda_{0},\gamma_{G}\right)$, and a Lorentzian
$L\left(\lambda_{0}+\delta\lambda_{0},\gamma_{L}\right)$:

$I_{\nu,\rm Line}(\lambda)=A\left[\beta
G\left(\lambda_{0},\gamma_{G}\right)+(1-\beta)L\left(\lambda_{0}+\delta\lambda_{0},\gamma_{L}\right)\right]$
(1)

where $A$ is the total amplitude of the line and $\beta$ the fractional
contribution of the Gaussian component which has a central wavelength
$\lambda_{0}$ and a full width at half maximum (FWHM), $\gamma_{G}$. The
Lorentzian component has a centre that is allowed to vary from the Gaussian
component by $\delta\lambda_{0}$, restricted to $\pm 0.005\mu$m and has a FWHM
$\gamma_{L}$. The centre of each line, $\lambda_{0}$ is allowed to vary by
$0.01\mu$m and each width, $\gamma_{G}$, $\gamma_{L}$ can vary by 50$\%$ of
its initial value. We use initial values from the FWHM of 0.0015 $\mu$m, 0.003
$\mu$m, 0.008 $\mu$m and 0.010 $\mu$m for Channels 1 to 4 respectively. The
full list of lines fitted are shown in Table 1 where the new lines not present
in the original version of PAHFIT are highlited in bold.

We find that the Lorentzian component is required to fit the base of the
lines. This may compensate for issues with the current wavelength calibration
in the data reduction pipeline (Rigby et al., 2022; Pereira-Santaella et al.,
2022; García-Bernete et al., 2022a). As the majority of the flux of the line
is not contained within the wings, this has a negligible impact on the
measured fluxes.

Table 1: Spectral Lines

Line | $\lambda_{0}$ | FWHM ($\pm 50\%$)
---|---|---
|  $\mu$m |  $\mu$m
(1) | (2) | (3)
H2 S(8) | 5.053 | 0.0015
$[$Fe II$]$ | 5.340 | 0.0015
$[$Fe VIII$]$ | 5.447 | 0.0015
H2 S(7) | 5.551 | 0.0015
$[$Mg V$]$ | 5.609 | 0.0015
H2 S(6) | 6.109 | 0.0015
H2 S(5) | 6.909 | 0.0015
$[$Ar II$]$ | 6.985 | 0.0015
HI 6-5 | 7.460 | 0.003
$[$Ar V$]$ | 7.901 | 0.003
H2 S(4) | 8.026 | 0.003
$[$Ar III$]$ | 8.991 | 0.003
H2 S(3) | 9.665 | 0.003
$[$S IV$]$ | 10.551 | 0.003
H2 S(2) | 12.278 | 0.008
HI 7-6 | 12.370 | 0.008
$[$Ne II$]$ | 12.813 | 0.008
$[$Cl II$]$ | 14.322 | 0.008
$[$Ne III$]$ | 15.555 | 0.008
H2 S(1) | 17.035 | 0.008
$[$P III$]$ | 17.885 | 0.010
$[$Fe II$]$ | 17.936 | 0.010
$[$S III$]$ | 18.713 | 0.010
$[$Fe III$]$ | 22.925 | 0.010
$[$Ne IV$]$ | 24.318 | 0.010

  * Column (1): Line name, Column (2): Vacuum wavelength, Column (3): Line FWHM allowed to vary by $\pm$50$\%$ of the value shown

  * Additional emission lines not present in the original PAHFIT (Smith et al., 2007) are shown in bold.

####  3.1.2 PAH Emission Features

The original version of PAHFIT modelled PAH features with a series of Drude
profiles as shown in Table 3 of Smith et al. (2007). With the increased
spectral resolution of JWST/MIRI the PAH emission shows more resolved
structure than with Spitzer data (also noted by García-Bernete et al., 2022a).
We therefore update the PAH model with additional components which we show in
Table. 2. For each Drude feature we allow the central wavelength to vary by
$\pm 0.05\mu$m to account for uncertainties in the wavelength calibration with
the JWST pipeline (Rigby et al., 2022). We also allow the width of each
feature to vary. In this case we allow the width to vary between $+10\%$ and
$-60\%$ of the values shown in Table. 2 where this range allows the PAH
features to be narrower than the previous version of the tool due to the
increased spectral resolution afforded by MIRI/MRS.

Table 2: PAH Features

$\lambda_{0}$ ($\pm 0.05\mu$m) | Average $\lambda_{0}$ | FWHM ($+10\%,-60\%$) | Average FWHM
---|---|---|---
$\mu$m |  $\mu$m |  $\mu$m |  $\mu$m
(1) | (2) | (3) | (4)
5.24 | 5.253 | 0.058 | 0.032
5.27 | 5.266 | 0.179 | 0.138
5.70 | 5.699 | 0.200 | 0.153
5.87 | 5.856 | 0.200 | 0.119
6.00 | 6.019 | 0.200 | 0.085
6.18 | 6.214 | 0.100 | 0.105
6.30 | 6.329 | 0.187 | 0.206
6.69 | 6.705 | 0.468 | 0.515
7.42 | 7.427 | 0.935 | 0.942
7.60 | 7.587 | 0.334 | 0.249
7.85 | 7.841 | 0.416 | 0.432
8.33 | 8.293 | 0.417 | 0.423
8.61 | 8.602 | 0.336 | 0.296
10.68 | 10.630 | 0.214 | 0.229
11.00 | 10.978 | 0.100 | 0.063
11.15 | 11.191 | 0.030 | 0.023
11.20 | 11.211 | 0.030 | 0.028
11.22 | 11.244 | 0.100 | 0.043
11.25 | 11.291 | 0.135 | 0.077
11.33 | 11.380 | 0.363 | 0.362
11.99 | 11.956 | 0.540 | 0.594
12.62 | 12.635 | 0.530 | 0.573
12.69 | 12.738 | 0.120 | 0.111
13.48 | 13.439 | 0.539 | 0.593
14.04 | 13.990 | 0.225 | 0.248
14.19 | 14.207 | 0.355 | 0.390
14.65 | 14.639 | 0.500 | 0.384
15.90 | 15.851 | 0.318 | 0.276
16.45 | 16.459 | 0.230 | 0.156
17.04 | 17.023 | 1.108 | 0.877
17.375 | 17.389 | 0.209 | 0.190
17.87 | 17.887 | 0.286 | 0.258
18.92 | 18.937 | 0.359 | 0.395

  * Column (1): Initial central wavelength for the fit with range $\pm 0.05\mu$m. Column(2): Average best fit central wavelength. Column (3): Initial FWHM allowed to vary within $+10\%$ and $-60\%$). Column (4): Average best fit FWHM.

  * Additional Drude profiles not present in the original PAHFIT (Smith et al., 2007) are shown in bold.

####  3.1.3 Continuum

The continuum is modelled as a sum of modified blackbodies with a range of
fixed temperatures [35, 40, 50, 65, 90, 135, 200, 300, 500 K] and an
emissivity of $\nu^{2}$. An additional blackbody at 5000K is included to
account for the presence of a potential stellar continuum. This continuum is
subject to extinction. The original PAHFIT used a modified version of the
Kemper et al. (2004) Galactic centre extinction curve, while subsequent
versions of the code used various other extinction curves (e.g. Gallimore et
al., 2010). For our new version of PAHFIT we use an extinction curve generated
using templates from extremely obscured sources. We do this as many of the
extinction curves in the literature (e.g. Ossenkopf et al., 1992; Kemper et
al., 2004; Chiar & Tielens, 2006) do not include crystalline features at 11
$\mu$m, 16 $\mu$m, 19 $\mu$m or 23 $\mu$m, which Spoon et al. (2022) found
were a common component of the ISM of galaxies. Moreover, their presence can
be seen by eye in the spectrum for source A (e.g absorption at 23 $\mu$m) and
so are an imperative for the model to achieve a good fit to the data. The
creation of this extinction curve is described in Appendix A. In addition we
use the water ice and CH templates from Donnan et al. (2022) derived from NGC
4418 as these features are very strong for this source. We assume a screen of
dust as the optical depth template is derived with this assumption. Therefore,
the extinction has the form

$e^{-\tau(\lambda)}=e^{-\tau_{\rm Ice}-\tau_{\rm CH}-\tau_{\rm Ext}}$ (2)

where $\tau_{\rm Ice}$ is the optical depth of the ice template, $\tau_{\rm
CH}$ is the optical depth of the CH and $\tau_{\rm Ext}$ is the optical depth
of the extinction curve created in Appendix A. $\tau_{\rm Ext}$ includes the
silicate features at 9.8 $\mu$m, 18 $\mu$m and a power law component.

Combining all these components gives the total model as

$f_{\nu}(\lambda)=\sum_{i=1}^{N_{\rm Lines}}I_{\nu,\rm
Line}^{(i)}(\lambda)+\sum_{i=1}^{N_{\rm PAH}}I_{\nu,\rm
PAH}^{(i)}(\lambda)+\\\ \left(\sum_{i=1}^{N_{\rm
Dust}}A_{i}\frac{B_{\nu}(T_{i},\lambda)}{\lambda^{2}}+A_{*}B_{\nu}(5000\rm{K},\lambda)\right)e^{-\tau_{\lambda}}$
(3)

where $N_{\rm Dust}=8$ for fixed $T_{i}=35,40,50,65,90,135,200,300$K with
scale factors $A_{i}$, $B_{\nu}$ is the Planck function and $A_{*}$ is the
scale factor for the stellar blackbody at a fixed temperature of 5000 K. The
emission lines are a sum of each line function, $I_{\nu,\rm
Line}^{(i)}(\lambda)$ (equation (1) for line $(i)$ and the PAH features are a
sum of each Drude function, $I_{\nu,\rm PAH}^{(i)}(\lambda)$ for PAH feature
$(i)$.

For each spectrum we sample the posterior probability using the No-U-Turn
Sampling (NUTS) Markov-Chain Monte Carlo (MCMC) from numpyro (Phan et al.,
2019). This allows posterior probabilities to be obtained from properties
inferred from the spectra.

For each spectrum the integrated flux and equivalent width of each emission
feature is calculated as a numeric integral. The integrated flux is simply

$f^{\rm PAH}=\int f^{\rm PAH}_{\nu}d\nu,$ (4)

for a PAH profile $f^{\rm PAH}_{\nu}$. Before calculating this integral we
divide the feature profile by the extinction factor $e^{-\tau_{\textrm{Ext}}}$
to recover extinction corrected fluxes. We do not include the ice + aliphatic
hydrocarbon absorption in this extinction factor as these are assumed to only
affect the continuum and are not typical of star-forming regions in galaxies
(e.g. Spoon et al., 2022).

To calculate the equivalent width of the PAH features we do not correct by
extinction where we just use the observed flux:

$\rm{EW}=\int\frac{f^{\rm PAH}_{\nu}}{f^{\rm cont}_{\nu}}d\lambda,$ (5)

where $f^{\rm cont}_{\nu}$ is the continuum, corrected by the ice+CH
absorption. This means that the EW of the 6.2 $\mu$m PAH uses a continuum only
subject to the extinction curve containing the power law + silicate
absorption, inline with previous work (e.g. Spoon et al., 2022; Donnan et al.,
2022, and references therein.).

We also measure the strength of the extinction via the strength of the 9.8
$\mu$m silicate feature as defined in (Spoon et al., 2007):

$S_{\rm{sil}}=\ln\left(\frac{f_{9.8,\rm{obs}}^{\rm
cont}}{f_{9.8,\rm{int}}^{\rm cont}}\right),$ (6)

which is simply the natural log of the observed to the underlying continuum at
9.8 $\mu$m. It is worth noting here that this definition differs slightly from
(Spoon et al., 2007) as we use the intrinsic continuum rather than an
interpolated continuum. The intrinsic continuum here contains the power law
component as well as the silicate absorption, which raises the magnitude of
$S_{\rm{sil}}$ by 10$\%$ assuming our extinction curve as discussed in
Appendix A. This value is in-fact equal to the negative of the optical depth
at 9.8 $\mu$m as we model the extinction as a screen.

We fit the spectra of the three regions as described in Section 2. We show
these fits in Appendix B. In columns (2) and (4) of Table 2 we show the
average PAH centres and average FWHMs from these three fits which can be used
to update future PAH models in various spectral fitting codes. We also fit
various star-forming clumps from the the MRS spectra of NGC 7469 (García-
Bernete et al., 2022c) to provide a small but representative reference sample
of pure star-forming regions.

In Fig. 3 we show the 6.2 $\mu$m, 7.7 $\mu$m, 11.3 $\mu$m, 12.7 $\mu$m and 17
$\mu$m PAH bands for region C with the best fit model. We show region C as it
showed the strongest PAH emission features. With the increased spectral
resolution of JWST/MIRI we can reveal new resolved features in the PAH
profiles. In particular the 6.2 $\mu$m PAH is now asymmetric which we now fit
with two Drude profiles rather than a single one. This asymmetry can be
somewhat seen in ISO (Infrared Space Observatory) observations (e.g.
Verstraete et al., 2001), which had a higher spectral resolution than Spitzer
in its low resolution mode. The 11.3 $\mu$m feature also shows more resolved
structure which we fit with four profiles compared to two in the original
PAHFIT (Smith et al., 2007).

Figure 3: Panels showing the main PAH bands at 6.2 $\mu$m, 7.7$\mu$m, 8.6
$\mu$m, 11.3 $\mu$m, 12.7 $\mu$m and 17 $\mu$m for region C. The best fitting
PAHFIT model is shown in red with it constituent components labelled in the
legend. The emission lines are identified in Table 1.

We typically found that the errors of measured properties from the MCMC
routine were very small. The standard deviation of the normalised residuals
was much greater than 1 for all of the model fits which suggests the errors of
the data are underestimated and thus our parameter errors are too small. We
therefore scaled the measured parameter errors by that required to achieve
normalised residuals with a standard deviation of 1 to give a more
conservative estimate for the uncertainties.

###  3.2 Emission Line Maps

To study how the emission line fluxes are spatially distributed we generate
maps using a local continuum and integrating to obtain the line flux. For each
spaxel we extract a spectrum at wavelengths around the emission line. We mask
the data points of the emission line, smooth the data using a Gaussian kernel
and interpolate across the masked region with a cubic spline to produce a
local continuum. We then subtract the local continuum from the data and
integrate to measure a line flux. This is done for each spaxel to produce a
map of the emission line flux.

To estimate error maps for these lines, we first estimate the noise in the
flux data ($f_{\nu}(\lambda)$) by finding the standard deviation of the data
after subtracting the local continuum and masking the emission line. We do
this rather than using the flux error bars as these are likely underestimated
as noted in the previous section. We then propagate the errors in $f_{\nu}$ to
find the error in the integrated flux of the emission line for each spaxel to
obtain an error map.

To ensure accurate line fluxes in such a heavily obscured environment
extinction corrections need to be applied. Obtaining accurate extinction
corrections is a challenge in such a complex environment as different
components of the spectrum may be affected by different levels of extinction.

To do this we generate a map of the optical depth at each spaxel by measuring
the silicate depth at 9.8 $\mu$m. This is done by linearly interpolating
between 6.7 $\mu$m and 13 $\mu$m to approximate an underlying continuum from
which the silicate depth was measured using equation (6). We assume a 10$\%$
contribution at 9.8 $\mu$m from a power law component and increase the
measured silicate depth by a factor 1/0.9 to estimate $\tau_{9.8}$. This
results in the map shown in Fig. 4, which we use to correct any emission line
map by extinction using our extinction curve as discussed in Appendix A.

Figure 4: Extinction map. Optical depth at 9.8 $\mu$m derived from the
silicate depth using an interpolated continuum as discussed in Section 3.2.
This map is used to correct any emission line maps by extinction. The white
crosses show the position of the continuum from source A and B.

##  4 Results

###  4.1 Obscuration

####  4.1.1 PAH EW Ratios

To investigate the nature of the two nuclei and the extended emission, in
particular the presence of a highly obscured nucleus, we measure the PAH EW
ratios of the 12.7 $\mu$m over the 11.3 $\mu$m and the 6.2 $\mu$m over 11.3
$\mu$m PAH features in order to construct a PAH EW diagram (García-Bernete et
al., 2022c). These ratios are shown in Fig. 5. For star-forming regions we
expect an approximately constant PAH EW ratio (Hernán-Caballero et al., 2020;
García-Bernete et al., 2022c). Because of the enhanced spectral resolution
afforded by MIRI/MRS, we had to redefine the CON selection region using JWST
reference spectra of star-forming regions. We do this as the increased
spectral resolution of MIRI MRS compared to Spitzer IRS results in different
values of the PAH EW ratios, and thus recalibration of the CON selection
criteria is required. In particular, the 12.7 $\mu$m PAH is better resolved
and thus has a higher EW than in the Spitzer data. To recalibrate the CON
selection region we use the clumps observed in the starburst ring of NGC 7469
(García-Bernete et al., 2022a). These are shown as the grey points in Fig. 5.
As more targets are observed with MIRI MRS, the sample of star-forming regions
will expand and thus improve the CON selection criteria.

From these star-forming clumps we measure the mean and standard deviation,
$\sigma$, of the EW(12.7)/EW(11.3) and EW(6.2)/EW(11.3) ratios where the
obscured nuclei selection region is defined as values lower than $n\sigma$
from this mean. We show 3$\sigma$, 4$\sigma$ and 5$\sigma$ thresholds in Fig.
5.

Figure 5:  Left: PAH EW diagram (García-Bernete et al., 2022c) for the three
regions in VV 114E, labeled A, B and C. The grey points show the values found
for star-forming clumps in NGC 7469 (García-Bernete et al., 2022a) with means
shown as the grey lines. The selection criteria for highly obscured nuclei are
shown as the purple box defined as either 3$\sigma$, 4$\sigma$ or 5$\sigma$
from the mean of the star-forming clumps. Right: Plot of the 12.7/11.3 PAH
flux ratio as a function of extinction as measured by the strength of the 9.8
$\mu$m absorption feature as given by equation 6. The grey points show the
star-forming clumps for NGC 7469, from which the mean flux ratio subject to
extinction from our custom extinction curve (see App. A) is plotted as the
black line with $\pm 20\%$ intervals shown as the dashed lines. This line
shows how the flux ratio is expected to change as a function of extinction.

We also show the expected 12.7/11.3 PAH flux ratio as a function of the
extinction as measured by the silicate strength (equation 6) in the right
panel of Fig. 5. To construct this, we first measure the intrinsic flux ratio
from the star-forming clumps of NGC 7469 as these show practically no
extinction with a silicate strength of very close to zero (García-Bernete et
al., 2022a). We then apply extinction to this ratio using our extinction curve
and measure the resulting PAH flux ratio. This is shown as the black line with
the dashed lines denoting $\pm$20$\%$. Obscured star-formation is expected to
lie along this line assuming the PAH emission traces star-formation (Hernán-
Caballero et al., 2020; García-Bernete et al., 2022c).

Source A shows lower EW(12.7)/EW(11.3) and EW(6.2)/EW(11.3) with both values
below 3$\sigma$ of the star-forming clumps. This indicates the presence of an
obscured nucleus. This is further shown in the right panel of Fig. 5 where the
12.7/11.3 PAH flux ratio is lower than expected for the level of extinction
(indicated by the silicate strength) measured on the continuum. This suggests
additional extinction on the continuum compared to the PAH emission and thus
the presence of a deeply buried nucleus.

Source B shows a raised EW(12.7)/EW(11.3) but a significantly lower
EW(6.2)/EW(11.3). This is where typical AGN are observed in the PAH EW ratio
diagram in García-Bernete et al. (2022c), which the authors attributed to a
flatter continuum around the 11.3 $\mu$m PAH, due to contribution from
silicates in emission. The low EW(6.2)/EW(11.3) is also expected for AGN as
the presence of hot dust elevates the continuum around 6 $\mu$m resulting in a
low 6.2 PAH EW. This is consistent with the blue MIRI image which Evans et al.
(2022) point out is typical of AGN.

Region C shows an EW(12.7)/EW(11.3) consistent with the star-forming clumps as
expected as this aperture traces the extended PAH emission which is likely
pure star-formation. It does however show a significantly lower
EW(6.2)/EW(11.3) which could be due to a bluer continuum, raising the EW of
the 6.2 $\mu$m PAH.

####  4.1.2 Spectral Decomposition

To further assess whether nucleus A contains a highly obscured source, we
model the spectrum by decomposing it into a star-forming and nuclear
component. We follow the same method as from Donnan et al. (2022) to do this.
This model fits the 5.2 - 14.2 $\mu$m spectrum with a star-forming continuum
using templates from a calibration sample (Hernán-Caballero et al., 2020) and
a nuclear continuum that is heavily obscured. The nuclear continuum is
composed of a smooth spline and a silicate template (derived from IRAS
08572+3915) where the optical depth, $\tau_{N}$, of this template is a free
parameter.

To properly decompose the spectra, we use a Bayesian process where priors are
used to break the degeneracy between the two components. We assume that the
PAH emission is associated with star-formation (the PAH ratios in Section 4.2
are consistent with this) and therefore we place a tight prior on the total
PAH flux, $f_{\rm PAH}$, to the star-forming continuum flux, $f_{\rm SF}$,
between 5.2 - 14.2 $\mu$m. We use a normal prior of $f_{\rm PAH}/f_{\rm
SF}=1.92\pm 0.56$ as inferred from the Hernán-Caballero et al. (2020)
calibration sample in Donnan et al. (2022). We use the silicate profile of
IRAS 08572+3915 for the nuclear component, where this component is constrained
through its contribution to the total flux where the star-forming continuum is
tied to the total PAH flux via the aforementioned prior.

This results in the decomposition in Fig. 6 where we measure a nuclear optical
depth of $\tau_{N}=5.581^{+0.006}_{-0.007}$ and a nuclear fraction of
$\beta=0.772^{+0.001}_{-0.001}$ ($\beta$ is defined as the fractional
contribution to the total continuum flux between 5.2 - 14.2 $\mu$m). It is
worth noting that these uncertainties are likely underestimated by a factor of
$\sim 10$ as mentioned previously. These results show that the spectrum of
source A can be explained by a low contribution from less obscured star-
formation in the foreground in addition to a highly obscured nucleus with a
deep silicate absorption feature, which dominates the spectrum. To place this
optical depth in context we show a modified plot from Fig. B1 of Donnan et al.
(2022) where the measured nuclear optical depths from Spitzer spectra of the
CONquest sample (Falstad et al., 2021) are plotted against the HCN-vib surface
brightness as measured in the millimetre, using the same IRAS 08572+3915
silicate template used here. This is shown in Fig. 7. Placing the
$\tau_{N}=5.6$ onto this plot finds a HCN-vib surface brightness of
$\Sigma_{\rm HCN-vib}\gtrsim 1L_{\odot}$ pc-2, which above the threshold to be
considered a CON (Falstad et al., 2021).

Figure 6: Fit to the spectrum of source A using the decomposition model
described in Section 4.1 and Donnan et al. (2022). This model fits the
spectrum up to 14.2 $\mu$m by decomposing the continuum into star-forming and
nuclear where the relative contributions are constrained by restricting the
total PAH flux to the total star-forming continuum flux. The total continuum
is given by the solid line and is the sum of the two components.  Figure 7:
Plot modified from Donnan et al. (2022) which shows the nuclear optical depth
of the CONquest sample (Falstad et al., 2021) measured using the spectral
decomposition method and plotted against the surface brightness of HCN-vib. We
show the CON definition of $\Sigma_{\rm HCN-vib}>1L_{\odot}$ pc-2 with the
vertical dashed grey line. The horizontal red line shows the value of
$\tau_{N}=5.6$ found for VV 114E A, which predicts a $\Sigma_{\rm HCN-
vib}>1L_{\odot}$, meeting the CON definition.

###  4.2 PAH Properties

The relative fluxes between different PAH emission features are known to
change with the properties of the PAH molecules and the incident radiation
field (see Li, 2020, for a recent review)..

In the immediate vicinity of AGN, PAH molecules are thought to be destroyed
(e.g. Roche et al., 1991; Voit, 1992; Siebenmorgen et al., 2004) however there
is evidence that they can survive close to the AGN, in particular the 11.3
$\mu$m feature (Hönig et al., 2010; Alonso-Herrero et al., 2014, 2016; Jensen
et al., 2017). Most recently García-Bernete et al. (2022a) found that PAHs can
survive in the innermost region of luminous Seyfert galaxies using JWST/MRS
data but that these PAHs are different to typical star-forming galaxies, where
they are observed to be larger and more neutral (García-Bernete et al., 2022b,
a).

For measuring the PAH properties in VV 114E, we avoid using the 7.7 $\mu$m PAH
as the continuum around this feature is difficult to constrain properly due to
both the CH absorption at $\sim 7.3\mu$m and the tail of the silicate
absorption feature. We therefore employ the 11.3/6.2 ratio as a tracer of the
ionisation state of the PAH population (e.g. Rigopoulou et al., 2021). This
was shown in García-Bernete et al. (2022b) to show distinct differences in the
ratios between AGN and pure star-forming galaxies.

In Fig. 8 we show the 11.3/6.2 $\mu$m flux ratio after correcting for
extinction. This ratio is expected to decrease with a greater fraction of
ionised PAH molecules (Rigopoulou et al., 2021; García-Bernete et al., 2022b).
We plot the measured ratios against the [Ne III]/[Ne II] line ratio which
traces the hardness of the radiation field (these lines are investigated
further in the following section) and compare against the average ratios for
star-forming galaxies and AGN from García-Bernete et al. (2022b). We find all
three sources to be consistent with typical star-forming galaxies where the
PAH molecules are less neutral (higher fraction of ionised molecules) than
typical AGN. This means that if any AGN is present then it must be extremely
obscured leaving the PAH emission unaffected.

One concern would be that the extinction is not quite correct for the obscured
source A. As discussed in Section 4.1, source A likely contains a buried
continuum source which means our PAH fluxes are potentially over-corrected for
extinction. This would change the 11.3/6.2 ratio, where over-correcting for
extinction would boost the 11.3 more than the 6.2. This would imply that the
true 11.3/6.2 ratio is slightly lower than plotted which moves it more inline
with the other points, further from the AGN average.

Figure 8: PAH flux ratios of the 11.3 $\mu$m PAH over the 6.2 $\mu$m against
the [Ne II]/[Ne II] line ratio which traces the hardness of the radiation
field. A greater fraction of ionised PAHs will result in a lower value. The
values labelled A, B and C show those obtained from the three regions in VV
114E while the grey points show 4 star-forming clumps from NGC 7469 (García-
Bernete et al., 2022a). The solid grey lines show the average value for star-
forming galaxies ($\sim 1$) and AGN ($\sim 2$) from García-Bernete et al.
(2022b). These ratios have been corrected by extinction.

###  4.3 Fine Structure Lines

To investigate the hardness of the radiation field withing VV 114E we
generated maps of [Ne II] (12.813 $\mu$m, 21.56 eV) and [Ne III] (15.555
$\mu$m, 40.96 eV). The lower ionisation potential lines [Ne II] and [Ne III]
are excellent tracers of star-formation due to their high relative brightness.
These maps are shown in Fig. 9. The errors in these lines are very low as seen
in the error maps in the same figure.

We do not detect any of the higher ionisation potential lines such as [Ne V]
(14.322 $\mu$m, 97.1 eV) or [Ne VI] (7.652 $\mu$m, 126.2 eV), which are
typical tracers of AGN activity (e.g. Sturm et al., 2002; Armus et al., 2006;
Inami et al., 2013). This suggests that there is no opening angle to any
potential AGN present and therefore if an AGN does exist, it must be highly
obscured.

The [Ne III]/[Ne II] line ratio is higher for source B compared to source A
and the immediate surrounding region. This points towards a harder radiation
field within this cluster where the lack of [Ne V] suggests that this is
likely a star-forming clump with a young stellar population rather than an
AGN. This is consistent with the much bluer colours seen in the imaging in
Fig. 1 and Evans et al. (2022). The [Ne III]/[Ne II] ratio for source A may
trace the foreground star-formation as found by the spectral decomposition in
Section 4.1 as it appears similar to the surrounding environment.

We additionally detect the faint [Cl II] (14.37 $\mu$m, 12.97 eV) line which
has a very low ionisation potential. Due to its faintness it was difficult to
detect with Spitzer however the line has been observed in NGC 4945 (Pérez-
Beaupuits et al., 2011) where it appeared to be more diffuse than [Ne II].
This is indeed the case with VV 114E where in Fig. 9 the emission is more
extended. The [Ne III]/[Cl II] ratio is also plotted and clearly shows an
increase for source B, consistent with the [Ne III]/[Ne II] ratio. It is worth
noting that the error in this line ratio is significantly larger than for [Ne
III]/[Ne II] due to the relative faintness of the [Cl II] line.

In Fig. 10 we show a map of [Fe II] (5.34 $\mu$m) which is a strong tracer of
shocks (Allen et al., 2008). There is a clear increase in [Fe II] emission
south of source B which indicates a shock front. The errors in this map are
very low as seen in the lower panel of Fig. 10 and therefore the spatial
variations in the flux of this line are statistically significant. This map
has not been corrected by extinction, however if it were corrected the shock
front would be even more prominent as this region has a higher optical depth
than behind the shock, resulting in a larger correction than behind the shock.

Figure 9: Low ionisation potential line maps and their ratios. The top panels
shows [Ne II] (12.813 $\mu$m, 21.56 eV), [Ne III] (15.555 $\mu$m, 40.96 eV),
and [Cl II] (14.37 $\mu$m, 12.97 eV) from left to right. The bottom panels
show ratios of these lines. For each map the associated error map is shown
directly below. Flux units are in 10-20 Wm-2. The continuum positions of
sources A and B are shown with the white crosses.  Figure 10: Map of the shock
tracer [Fe II] (5.34 $\mu$m) line. We labelled the shock front with the dashed
black line and the potential direction as discussed in Section 4.4. The
continuum positions of sources A and B are shown with the white crosses. This
map has not been corrected for extinction. In the lower panel we show the
error map.

###  4.4 Molecular Hydrogen (H2)

There are numerous transitions of molecular hydrogen from 5 $\mu$m to 17
$\mu$m. We generated maps of these emission lines in Fig. 11 where we show the
S(1) to S(8) transitions. These have been corrected for extinction using the
map in Fig. 4. It is worth noting that this likely leads to an over-correction
of the extinction for source A as the extinction map was inferred from
silicate depth and thus measures the extinction from the continuum. As
discussed in section 4.1 the continuum of source A is more embedded and
therefore using this extinction we are over-correcting the H2 line fluxes, in
particular the S(3) J=5, 9.665 $\mu$m transition, which sits in the 9.8 $\mu$m
silicate absorption band. This is seen in Fig. 11 where for S(3), source A
appears relatively brighter than in the other lines which are less sensitive
to the extinction correction. We therefore also show the uncorrected line maps
in Appendix C.

Fig. 11 illustrates a clear shock with a front south of source B as suggested
by the [Fe II] emission in Fig. 10. Moving to the lower energy transitions the
peak moves behind the shock front tracing the cooler gas.

Figure 11: Maps of the integrated flux of the molecular Hydrogen transitions
from the highest energy H2 S(8) (5.053 $\mu$m) to the lowest energy H2 S(1)
(17.0 $\mu$m). The white crosses mark the continuum positions of sources A and
B. These maps have been corrected by extinction using the optical depth map
shown in Fig. 4. The contours show the 6.2 $\mu$m PAH emission.

We also show a map of the 6.2 $\mu$m PAH emission as contours in Fig. 11 which
traces regions of star-formation. For most of the structure the molecular
hydrogen is consistent with the PAH emission as stars are expected to form
from the molecular gas and young stars excite the PAH molecules. However
behind the shock front there is a clear region of molecular hydrogen which
devoid of any 6.2 $\mu$m PAH emission.

Shocks are expected to destroy PAH molecules (e.g. Micelotta et al., 2010)
where the shock heated gas produces high energy collisions ($\sim$ 1 KeV) with
the PAH molecules that destroys the smaller molecules first. Interestingly we
do not see destruction of the 6.2 $\mu$m PAH band at the shock front but
rather we see strong PAH emission here with the destruction behind the front
where the PAH emission is lacking. A possible explanation is that the
destruction is delayed such that the PAH molecules are destroyed after the
shock has passed through the region in the hot post-shock gas.

####  4.4.1 H2 Excitation

Under the assumption of local thermodynamic equilibrium (LTE), the relative
fluxes of these transitions can be modelled where the molecular hydrogen
column density and temperature determine the relative population of the energy
levels and thus the relative fluxes observed. We use a two temperature gas
model to measure the temperature and the column density of molecular hydrogen
for each component. At each spaxel we fit the hydrogen transitions using the
PhotoDissociation Region Toolbox code222https://github.com/mpound/pdrtpy
(Kaufman et al., 2006; Pound & Wolfire, 2008, 2011) with the ortho-to-para
ratio fixed at 3. As the S(3) transition occurs at 9.665 $\mu$m which sits in
the peak of the silicate absorption, it is highly sensitive to the extinction
correction and so we include this point but with a large (50%) error to ensure
it does not bias the fits. As the S(0) J=2 transition lies outwith the
observing range of MIRI (28.21 $\mu$m), we use the anchor point from Spitzer
observations (Petric et al., 2018) to act as a guide for the model to prevent
nonphysical solutions. To match the flux with the JWST/MIRI data we scale the
S(0) flux by matching the S(1) measurement from Spitzer to the median of the
S(1) map from MIRI. We give the value a 50$\%$ error to account for this large
uncertainty and to make it act as simply an anchor point. For the two point
sources, A and B, we fit the sum of the line fluxes in a 3x3 pixel grid around
the continuum position to avoid PSF effects.

In Fig. 12 we show the measured total column density of molecular hydrogen,
the temperature of the cold gas component and the temperature of the hot gas
component from these fits for each spaxel. Along the shock front there appears
to be a cavity in the column density coinciding with an increased temperature
in both the cold and hot gas. It is worth noting that these two quantities are
somewhat degenerate from the fitting process and so should not be taken as two
independent measures of the gas properties. Behind the shock the molecular
hydrogen column density increases approaching $N_{H_{2}}\xrightarrow{}10^{21}$
cm-2 as the temperature drops. This is expected for a fast molecular shock
where the gas behind the shock reforms (See Fig. 1 of Hollenbach & McKee,
1989). The H2 column densities are lower but similar ($\sim 10^{21}$ cm-2) to
those found by Saito et al. (2015) using CO observations however their value
is sensitive to the $\alpha_{\rm CO}$ conversion factor used.

The lower panels of Fig. 12 show the fits for the two nuclei A (left) and B
(right). As can be seen in these plots the S(0) (J=2) is not predicted by the
model but does act to prevent the solution from being nonphysical. This
suggests that with just two components the model is only able to fit the
hotter gas and an additional cold component may be required to fit this
transition. Adding this additional component to the model would add 2 more
parameters which may be too many as the J=2 transition is highly uncertain and
therefore the constraints on this cooler component would be very weak. For
source A we also see the over-correction of the extinction where the S(3) J=5
transition appears too high.

The temperatures we measure are high ($T_{\rm hot}>1000$ K) as we are
sensitive to the high energy transitions in the fitting. Calculating a
temperature from the S(7) and S(5) fluxes (following Lambrides et al., 2019)
using a 4.5” aperture from Spitzer spectral mapping data (Donnan et al.,
2022), finds a temperature of 1078 K. Considering this is over a larger
spatial scale and is still $>1000$ K, reinforces this result.

Figure 12: Results from modelling the S(1) to S(8) molecular hydrogen
transitions. We show the total molecular hydrogen column density (left), the
temperature of the cold gas component (middle) and the hot gas component
(right). The continuum positions of sources A and B are shown with the white
crosses. In the lower panels we show the fits to the two point sources A
(left) and B (right) where each is fitted to a 3x3 region around the continuum
position. In these plots the best fit is is shown with the solid line composed
of the cold (green) and hot (orange) gas components. These fits assume an
ortho-to-para ratio of 3. The S(0) (J=2) point for these fits is from Spitzer
observations (Petric et al., 2018) and acts as a loose anchor to prevent
nonphysical solutions.

##  5 Discussion

CONs are rare amongst all galaxy types which suggests this is a short-lived
phase in galaxy evolution, however the high obscuration may hide extreme SMBH
growth and/or compact starburst activity. CONs are commonplace among LIRGs and
ULIRGs (Falstad et al., 2021; Donnan et al., 2022) and are prevalent in major
mergers. The evolutionary picture of these objects is relatively unknown, in
particular its origin and subsequent evolution. The obscured source A in VV
114E is an interesting object to study to ask this question.

Our analysis in Section 4.1 finds that the PAH EW ratios for source A are low,
indicative of a buried continuum source. The decomposition model finds a
nucleus that has a sufficient optical depth $\tau_{N}=5.6$ to reach the CON
definition (Donnan et al., 2022) and pass the threshold. This object therefore
may be anywhere between the start or the end of the CON phase. The lack of
high ionisation lines such as [Ne V] suggests if an AGN is present, the
radiation field is completely obscured and unable to affect the surrounding
environment. This is confirmed by looking at the imprint on the PAH molecules
that the central source makes. AGN are known to destroy PAH molecules very
close to the central engine or alter their properties in the circumnuclear
region (García-Bernete et al., 2022b, a) where neutral, large molecules are
observed. This suggests the possibility that the hard radiation field from the
AGN destroys the small molecules first and so the imprint of AGN feedback can
be seen through the PAH properties.

For source A we found the 11.3/6.2 PAH flux ratio to indicate a population of
ionised PAHs, typical of star-forming regions. This suggests that there may
not be any AGN present or that it must be extremely obscured such that the
PAHs have not affected by the AGN. The latter explanation is required to be
consistent with the prediction from Saito et al. (2015) which suggested the
possible presence of a hidden AGN through an elevated HCN/HCO+ ratio.

X-ray emission has been observed from this source however only upper limits
have been placed on the X-ray luminosity of an intrinsic AGN (Ricci et al.,
2021) with $\log L_{2-10}<41.11$ erg s-1 and $\log L_{10-24}<40.95$ erg s-1.
These values are typical of CON dominated sources.

In summary, we detect a completely obscured nucleus for source A which may be
a hidden AGN as suggested by Saito et al. (2015). This object therefore may be
the progenitor of a quasar where the AGN is yet to influence its host galaxy,
as we see no evidence for any feedback on its surroundings. Indeed ULIRGs have
long been thought to be such progenitors (Sanders et al., 1988). To evolve
towards an unobscured quasar requires clearing of the obscuring gas/dust which
will require kinematic analysis to identify possible outflows responsible.
Further study of these objects are therefore necessary to understand the
possible onset of AGN feedback and how this fits in the evolution of galaxies.

From [Fe II] emission and analysis of the transitions of molecular hydrogen,
H2, we detect a shock front south of nucleus B (illustrated in Fig. 10). The
origin of this shock is unclear. Our analysis suggests that source B is a
compact star-forming region and so supernovae (e.g. Reach et al., 2019) or
winds from massive stars (Hollenbach, 1997) may be responsible if the shock
originates from this region. However it appears more extended which suggests
that this shock may simply be a product of the galaxy-galaxy interaction. This
may form part of a larger shock, where in Fig. 1 the 7.7 $\mu$m PAH emission
traces this region and extends to the west. Indeed, Saito et al. (2017) detect
CH3OH emission across $\sim$ 3 kpc which extends between the two merging
components of VV 114. They therefore suggest this is due to a large scale
shock from the galaxy-galaxy interaction. It is also interesting to note that
the shock front appears very dusty as seen in Fig. 4 where the optical depth
increases.

##  6 Conclusions

We have used spatially resolved mid-infrared spectroscopy from JWST/MIRI MRS
observations to investigate the obscured environment of VV 114E. We model the
spectrum from 5 - 28 $\mu$m within three apertures (shown in Fig. 1) from
which we measured the properties of the PAH molecules. We also construct maps
of various molecular and ionised lines. Our main findings are:

  * •

From the low 12.7/11.3 and 6.2/11.3 PAH EW ratios we find that source A is
heavily obscured and hosts a buried nucleus. The PAH EW ratios are elevated
for source B while source C is typical of star-forming regions. By decomposing
the spectrum into nuclear and star-forming we measure a nuclear optical depth
of $\tau_{N}=5.6$, which is large enough to meet the threshold required for
CONs.

  * •

The 11.3/6.2 PAH flux ratios suggests an ionised population of PAHs typical of
star-forming regions for all three sources.

  * •

We do not detect any high ionisation potential lines such as [Ne V] or [Ne VI]
which further rules out the possibility that source B is an AGN however the
high obscuration of source A means a hidden AGN is consistent with this lack
of high ionisation lines.

  * •

We detect clear evidence of shock south of source B with elevated [Fe II]
emission. The molecular hydrogen lines trace this shock where the higher
energy transitions peak at the shock front and move behind towards the lower
energy transitions.

  * •

The 6.2 PAH emission does not trace the low energy H2 emission lines which
suggests the ionised PAH molecules are destroyed in the post-shock gas. This
may be a delayed process as the PAH emission is strong at the shock front and
only subsides behind the shock.

  * •

From modelling the transitions of molecular hydrogen we detect clear evidence
of a shock front south of source B with a ridge of increased gas temperatures
resulting in higher column densities behind the shock front where gas is
expected to reform.

These results suggest that nucleus B is likely a star-cluster while source A
contains an extremely obscured nucleus. This nucleus shows typical
characteristics of CONs and may be powered by an AGN. This initial analysis
shows the power of the JWST for studying these highly obscured environments,
where future kinematics studies are needed to understand the role of outflows
in the evolution of these objects.

## Acknowledgements

FRD acknowledges support from STFC through grant ST/W507726/1. DR and IGB
acknowledge support from STFC through grant ST/S000488/1. DR also acknowledges
support from the University of Oxford John Fell Fund. AAH acknowledges support
from grant PGC2018-094671-B-I00 funded by MCIN/AEI/ 10.13039/501100011033 and
by ERDF A way of making Europe, and grant PID2021-124665NB-I00. The authors
are grateful to the ERS team for the observing program with a
zero–exclusive–access period.

## Data Availability

All the raw data used in this work is available through the MAST
archive333https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html.

## 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
  * Allen et al. (2008) Allen M. G., Groves B. A., Dopita M. A., Sutherland R. S., Kewley L. J., 2008, ApJS, 178, 20
  * Alonso-Herrero et al. (2014) Alonso-Herrero A., et al., 2014, MNRAS, 443, 2766
  * Alonso-Herrero et al. (2016) Alonso-Herrero A., et al., 2016, MNRAS, 455, 563
  * Armus et al. (2006) Armus L., et al., 2006, ApJ, 640, 204
  * Armus et al. (2009) Armus L., et al., 2009, PASP, 121, 559
  * Charmandaris et al. (2004) Charmandaris V., Le Floc’h E., Mirabel I. F., 2004, ApJ, 600, L15
  * Chiar & Tielens (2006) Chiar J. E., Tielens A. G. G. M., 2006, ApJ, 637, 774
  * Donnan et al. (2022) Donnan F. R., et al., 2022, arXiv e-prints, p. arXiv:2211.09628
  * Evans et al. (2022) Evans A. S., et al., 2022, arXiv e-prints, p. arXiv:2208.14507
  * Falstad et al. (2021) Falstad N., et al., 2021, A&A, 649, A105
  * Frayer et al. (1999) Frayer D. T., Ivison R. J., Smail I., Yun M. S., Armus L., 1999, AJ, 118, 139
  * Gallimore et al. (2010) Gallimore J. F., et al., 2010, ApJS, 187, 172
  * García-Bernete et al. (2022a) García-Bernete I., et al., 2022a, arXiv e-prints, p. arXiv:2208.11620
  * García-Bernete et al. (2022b) García-Bernete I., Rigopoulou D., Alonso-Herrero A., Pereira-Santaella M., Roche P. F., Kerkeni B., 2022b, MNRAS, 509, 4256
  * García-Bernete et al. (2022c) García-Bernete I., Rigopoulou D., Aalto S., Spoon H. W. W., Hernán-Caballero A., Efstathiou A., Roche P. F., König S., 2022c, A&A, 663, A46
  * Hernán-Caballero et al. (2020) Hernán-Caballero A., et al., 2020, MNRAS, 497, 4614
  * Hollenbach (1997) Hollenbach D., 1997. pp 181–198
  * Hollenbach & McKee (1989) Hollenbach D., McKee C. F., 1989, ApJ, 342, 306
  * 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
  * Inami et al. (2013) Inami H., et al., 2013, ApJ, 777, 156
  * Jensen et al. (2017) Jensen J. J., et al., 2017, MNRAS, 470, 3071
  * Jones & Merrill (1976) Jones T. W., Merrill K. M., 1976, ApJ, 209, 509
  * Kaufman et al. (2006) Kaufman M. J., Wolfire M. G., Hollenbach D. J., 2006, ApJ, 644, 283
  * Kemper et al. (2004) Kemper F., Vriend W. J., Tielens A. G. G. M., 2004, ApJ, 609, 826
  * Lambrides et al. (2019) Lambrides E. L., Petric A. O., Tchernyshyov K., Zakamska N. L., Watts D. J., 2019, MNRAS, 487, 1823
  * Li (2020) Li A., 2020, Nature Astronomy, 4, 339
  * 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
  * Micelotta et al. (2010) Micelotta E. R., Jones A. P., Tielens A. G. G. M., 2010, A&A, 510, A36
  * Ossenkopf et al. (1992) Ossenkopf V., Henning T., Mathis J. S., 1992, A&A, 261, 567
  * Pereira-Santaella et al. (2022) Pereira-Santaella M., et al., 2022, arXiv e-prints, p. arXiv:2208.04835
  * Pérez-Beaupuits et al. (2011) Pérez-Beaupuits J. P., Spoon H. W. W., Spaans M., Smith J. D., 2011, A&A, 533, A56
  * Petric et al. (2018) Petric A. O., et al., 2018, AJ, 156, 295
  * Phan et al. (2019) Phan D., Pradhan N., Jankowiak M., 2019, arXiv e-prints, p. arXiv:1912.11554
  * Pound & Wolfire (2008) Pound M. W., Wolfire M. G., 2008, in Argyle R. W., Bunclark P. S., Lewis J. R., eds, Astronomical Society of the Pacific Conference Series Vol. 394, Astronomical Data Analysis Software and Systems XVII. p. 654
  * Pound & Wolfire (2011) Pound M. W., Wolfire M. G., 2011, PDRT: Photo Dissociation Region Toolbox (ascl:1102.022)
  * Reach et al. (2019) Reach W. T., Tram L. N., Richter M., Gusdorf A., DeWitt C., 2019, ApJ, 884, 81
  * Ricci et al. (2021) Ricci C., et al., 2021, MNRAS, 506, 5935
  * Rich et al. (2011) Rich J. A., Kewley L. J., Dopita M. A., 2011, ApJ, 734, 87
  * Rieke et al. (2015) Rieke G. H., et al., 2015, PASP, 127, 584
  * Rigby et al. (2022) Rigby J., et al., 2022, arXiv e-prints, p. arXiv:2207.05632
  * 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. (2021) Rigopoulou D., et al., 2021, MNRAS, 504, 5287
  * Roche et al. (1991) Roche P. F., Aitken D. K., Smith C. H., Ward M. J., 1991, MNRAS, 248, 606
  * Saito et al. (2015) Saito T., et al., 2015, ApJ, 803, 60
  * Saito et al. (2017) Saito T., et al., 2017, ApJ, 834, 6
  * Sanders et al. (1988) Sanders D. B., Soifer B. T., Elias J. H., Neugebauer G., Matthews K., 1988, ApJ, 328, L35
  * Siebenmorgen et al. (2004) Siebenmorgen R., Krügel E., Spoon H. W. W., 2004, A&A, 414, 123
  * 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
  * Stetson (1987) Stetson P. B., 1987, PASP, 99, 191
  * Stierwalt et al. (2013) Stierwalt S., et al., 2013, ApJS, 206, 1
  * Strauss et al. (1992) Strauss M. A., Huchra J. P., Davis M., Yahil A., Fisher K. B., Tonry J., 1992, ApJS, 83, 29
  * Sturm et al. (2002) Sturm E., Lutz D., Verma A., Netzer H., Sternberg A., Moorwood A. F. M., Oliva E., Genzel R., 2002, A&A, 393, 821
  * Verstraete et al. (2001) Verstraete L., et al., 2001, A&A, 372, 981
  * Voit (1992) Voit G. M., 1992, MNRAS, 258, 841
  * 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

##  Appendix A Extinction Curve

There are many extinction curves presented in the literature (e.g. Chiar &
Tielens, 2006; Kemper et al., 2004) however we found that none of these
produced good fits for highly obscured objects and in particular the obscured
nucleus (A) from VV 114E. We therefore created an empirically based extinction
curve using a heavily obscured galaxy that has no apparent PAH emission,
namely IRAS 08572+3915.

We use the Spitzer IRS spectra of IRAS 08572+3915 from the IDEOS database
(Spoon et al., 2022), where we find a local continuum using a cubic spline
interpolation with anchor points 5.5, 7.8, 13.0, 14.5 and 26.5 $\mu$m. This
underlying continuum is shown in Fig. 13 with the blue dashed line. We find
the optical depth as a function of wavelength by calculating the log of the
ratio of the underlying continuum to the data. This results in the middle
panel of Fig. 13. We take the silicate absorption features as a smoothed
version of the optical depth data. We then construct an extinction curve by
adding a power law component with exponent 1.7 making up 10% of the optical
depth at 9.8$\mu$m following the original PAHIFT (Smith et al., 2007) i.e

$\tau_{\lambda}=0.1\left(9.8/\lambda\right)^{1.7}+0.9\tau_{\rm Sil}(\lambda),$
(7)

where $\tau_{\rm Sil}$ is the extracted silicate profile from IRAS 08572+3915.
The extinction curve is shown in the right panel of Fig. 13. Comparing the
extinction curve to others from the literature it is clear that ours contains
crystalline absorption features at 11 $\mu$m, 16 $\mu$m, 19 $\mu$m, 23 $\mu$m
and 28 $\mu$m which are not present in the others. These absorption features
are typically weak in less obscured systems however are common in the ISM of
galaxies (Spoon et al., 2022). Our extinction curve also features a
significantly lower 18 $\mu$m silicate absorption feature relative to the 9.8
$\mu$m silicate feature than the other curves. This may be due to the presence
of emission from cold silicates that emit at 18 $\mu$m but not 10$\mu$m,
filling up the 18 $\mu$m silicate absorption band relative to the 9.8 $\mu$m
band (e.g. Jones & Merrill, 1976). The final model continuum is a combination
of numerous blackbodies at different temperatures and the chosen extinction
curve which will compensate for these differences. Therefore this is not an
issue in this work as we are interested in the resulting full continuum rather
than accurately diagnosing the dust constituency based on the shape of the
extinction curve.

Figure 13: Left: Spitzer IRS spectrum of IRAS 08572+3915 in black with the
interpolated cubic spline continuum shown as the dashed blue line. Middle:
Optical depth profile extracted from the log ratio of the interpolated
continuum to the spectrum. The inferred silicate absorption is shown with the
gold line. Right: Extinction curve constructed from a power law + silicate
template. This is compared to three popular choices from the literature Chiar
& Tielens (2006), Kemper et al. (2004), Ossenkopf et al. (1992).

##  Appendix B Model Fits

Figure 14: Fits to the spectra of VV 114E with our modified PAHFIT model from
regions A (top), B (middle) and C (bottom) as shown in Fig. 1. The full model
is shown in red with its constituent parts labelled in the legend. The right
axis of each plot shows the extinction factor, $e^{-\tau_{\lambda}}$, with the
dashed black line.

##  Appendix C Uncorrected H2 Maps

Figure 15: Same as Fig. 11 but without any extinction correction applied.