GOALS-JWST: Mid-infrared Molecular Gas Excitation Probes the Local Conditions of Nuclear Star Clusters and the Active Galactic Nucleus in the LIRG VV 114

The enormous increase in mid-IR sensitivity and spatial and spectral resolution provided by the JWST spectrographs enables, for the first time, detailed extragalactic studies of molecular vibrational bands. This opens an entirely new window for the study of the molecular interstellar medium in luminous infrared galaxies (LIRGs). We present a detailed analysis of rovibrational bands of gas-phase CO, H2O, C2H2, and HCN toward the heavily obscured eastern nucleus of the LIRG VV 114, as observed by NIRSpec and the medium resolution spectrograph on the Mid-InfraRed Instrument (MIRI MRS). Spectra extracted from apertures of 130 pc in radius show a clear dichotomy between the obscured active galactic nucleus (AGN) and two intense starburst regions. We detect the 2.3 μm CO bandheads, characteristic of cool stellar atmospheres, in the star-forming regions, but not toward the AGN. Surprisingly, at 4.7 μm, we find highly excited CO (T ex ≈ 700–800 K out to at least rotational level J = 27) toward the star-forming regions, but only cooler gas (T ex ≈ 200 K) toward the AGN. We conclude that only mid-infrared pumping through the rovibrational lines can account for the equilibrium conditions found for CO and H2O in the deeply embedded starbursts. Here, the CO bands probe regions with an intense local radiation field inside dusty young massive star clusters or near the most massive young stars. The lack of high-excitation molecular gas toward the AGN is attributed to geometric dilution of the intense radiation from the bright point source. An overview of the relevant excitation and radiative transfer physics is provided in an appendix.


Introduction
The L IR = 4.5 × 10 11 L e luminous infrared galaxy (LIRG) VV 114 (IC 1623, Arp 236) is a merging system undergoing intense star-forming activity.It consists of two gas-rich galaxies, with a projected nuclear separation of 8 kpc at a distance of 80 Mpc.At optical wavelengths, its western component (VV 114W) shows bright spiral arms and many luminous young star clusters, while its eastern component (VV 114E) is heavily obscured by dust lanes.At mid-and far-IR wavelengths, however, the eastern component dominates the emission.With a 7 μm-to-UV flux density ratio of ∼800 for VV 114E and ∼10 for VV 114W (Charmandaris et al. 2004), VV 114 has the most extreme color contrast between galaxies taking part in the merger observed among local LIRGs.
Early radio and infrared images were already able to resolve the dusty eastern nucleus into an NE and SW component (Condon et al. 1991;Knop et al. 1994).The latter has a complex structure, consisting of a bright unresolved point-like source and several secondary emission peaks (Scoville et al. 2000;Soifer et al. 2001).The point source is bright at infrared wavelengths, but does not coincide with any of the peaks in the 8.4 GHz radio image.High-resolution (θ FWHM ≈ 0 2) Very Large Array (VLA) 33 GHz imaging by Song et al. (2022) resolved the SW core into four distinct components.The 33 GHz continuum, which traces the thermal free-free emission from H II regions, is much fainter at the location of the infrared point source than at both the NE core and at an extended knot S of the infrared peak, suggesting relatively weak star-forming activity at the point source (Evans et al. 2022).
Atacama Large Millimeter/submillimeter Array (ALMA) observations reveal abundant cold molecular gas in VV 114.CO is ubiquitous in the system, while the dense gas tracer HCN is found to be concentrated in the nucleus of VV 114E; HCO + is found in the eastern nucleus and in the overlap region between the two galaxies (Saito et al. 2015).An enhanced methanol abundance, indicative of kiloparsec-scale shocks, is found in the overlap region as well (Saito et al. 2017).Highresolution 0 2 ALMA Band 7 observations of HCN(4-3) and HCO + (4-3) were able to resolve the VV 114E nucleus into the NE core and a complex SW structure (Saito et al. 2018), with comparatively little emission coming from the bright infrared point source.
The nature of the nuclear components of VV 114E has been a subject of debate.Mid-infrared Infrared Space Observatory studies (Le Floc'h et al. 2002) and X-ray observations (Grimes et al. 2006) found no evidence for an energetically important active galactic nucleus (AGN), but could not exclude a heavily obscured one either.Iono et al. (2013) found a significantly enhanced HCN(4-3)/HCO + (4-3) ratio toward the NE core, and attributed this to the presence of a deeply embedded AGN.Near-and mid-infrared JWST imaging and spectroscopic observations, however, show a different picture.Evans et al. (2022) found mid-IR colors that are consistent with a starburst toward the NE core, while those toward the SW core are in line with an AGN.Moreover, the extremely low equivalent widths of 3.3 and 6.2 μm polycyclic aromatic hydrocarbon (PAH) features found by Rich et al. (2023) are strong evidence for the presence of an AGN in the SW core.These JWST observations indicate that the IR-bright point source in the SW core harbors an AGN, while the other nuclear continuum peaks are starburst regions.
The integral field spectrographs on board JWST provide a new method of probing local gas conditions.Their high sensitivity and spatial resolution and moderate spectral resolution enable a detailed analysis of rovibrational molecular bands in external galaxies.Toward U/LIRGs, the intense infrared background radiation causes these lines to be seen in absorption, allowing for direct measurements of column density.Since the high-J lines will not be affected by cold foreground gas, the excitation conditions of the gas can be studied in detail.Such studies have previously only been possible in Galactic star-forming regions (e.g., Lahuis & van Dishoeck 2000;González-Alfonso et al. 2002).
The Q-branches of C 2 H 2 , HCN, and CO 2 , blends of closely spaced pure vibrational lines, have been detected by Spitzer toward a handful of local U/LIRGs (Lahuis et al. 2007), but at lower spatial and spectral resolution.Rich et al. (2023) previously reported the detection of the ν 2 14.0 μm band of HCN and the ν 5 13.7 μm band of C 2 H 2 in absorption toward the NE core in VV 114E.The fundamental band of CO at 4.67 μm, tracing transitions between the ground state and the first vibrationally excited state, has been observed in absorption toward several heavily obscured U/LIRGs from the ground (e.g., Spoon et al. 2003;Geballe et al. 2006;Shirahata et al. 2013;Baba et al. 2022;Ohyama et al. 2023).Pereira-Santaella et al. (2024) demonstrated JWST's capabilities for studying the CO band in emission.
In this paper, we present JWST NIRSpec and Mid-InfraRed Instrument (MIRI) absorption spectra of gas-phase CO, H 2 O, C 2 H 2 , and HCN observed toward the NE core and the SW complex in the nucleus of VV 114E.Using rotation diagram analysis (CO and H 2 O) and spectral fitting (C 2 H 2 and HCN), we probe the excitation conditions of the molecular gas, and investigate the responsible excitation mechanism.

Observations and Data Reduction
VV 114E was observed with MIRI (Rieke et al. 2015;Labiano et al. 2021) in medium resolution spectroscopy (MRS) mode (Wells et al. 2015) on 2022 July 5, and with NIRSpec in integral field unit (IFU) mode (Jakobsen et al. 2022;Böker et al. 2022) on 2022 July 19 as part of the GOALS Early Release Science Program 1328 (co-PIs Armus and Evans;doi:10.17909/vhw3-g317).This work is focused on the spectroscopic data, but we include the results of MIRI imaging to indicate the regions of spectral extraction (Figure 1).For a description of the MIRI imaging data, see Evans et al. (2022).

NIRSpec Data
NIRSpec IFU observations cover the wavelength range from 0.97 to 5.27 μm, employing three combinations of gratings and filters: G140H/F100LP, G235H/F170LP, and G395H/ F290LP.This last configuration covers the fundamental band of CO at 4.7 μm.A four-point dither pattern was used, as well as dedicated Leakcal images to correct for failed open microshutter assembly (MSA) shutters.The field of view (FOV) of the data products is ≈3 6 × 3 8 at a position angle (counterclockwise from N) PA ≈ 32°.
We downloaded the uncalibrated data-processed by Science Data Processing (SDP) version 2023 1a-from the MAST portal.We reduced the data using JWST Science Calibration Pipeline version 1.11.1.dev1+g9975346and calibration reference data system (CRDS) reference file jwst_1097.pmap.This is an updated version of the pipeline with respect to the one applied by Rich et al. (2023).
Both the science and Leakcal exposures were first put through the Detector1 pipeline, which applies detector-level calibrations to produce count rate files from nondestructive ramp readouts.The Spec2 pipeline then processes these rate files to produce fully calibrated individual exposures.At this step, additional instrumental corrections are applied, as well as flux and wavelength calibration.The Leakcal count rate files are used here to correct for failed open MSA shutters.The Spec3 pipeline finally combines the Stage 2 calibrated science exposures into a fully calibrated data cube.Stage 3 processing includes an outlier detection step.After some experimentation, we decided on an outlier flag threshold of 99.8% and a kernel size of 7 × 7. The final data cube spaxel size is 0 1, corresponding to ≈40 pc for VV 114.The NIRSpec pointspread function (PSF) has an FWHM of ∼0 17 at the longwavelength end of 5.3 μm (Jakobsen et al. 2022).

MIRI MRS Data
The MIRI MRS observations cover the full 4.9-27.9μm wavelength range of the four IFU channels, using the three grating settings SHORT (A), MEDIUM (B), and LONG (C) to do so.Dedicated off-target background observations were taken, and a four-point dither pattern was employed.Each channel has a different FOV; it is smallest for channel 1 with an FOV of ∼5″ × 4″ at a position angle PA = 255°.The FWHM of the MRS PSF varies between 0 35 in channel 1A and 0 9 in channel 4C (Argyriou et al. 2023).
The uncalibrated data products-processed by JWST SDP version 2022 4a-were downloaded through the MAST portal.These were then reduced using developmental version 1.11.1.dev1+g9975346 of the JWST Science Calibration Pipeline (Bushouse et al. 2022) and CRDS reference file jwst_1094.pmap.
The MIRI MRS data were processed in a similar fashion to the NIRSpec data.The Detector1 pipeline produced count rate files, and we applied the Spec2 step to subsequently turn these into calibrated individual exposures.We also applied the residual fringe correction (rfc) at this stage.For the background subtraction, we applied the master background subtraction step in the Spec3 pipeline.Outlier detection was performed at this stage as well, with a threshold of 99.8% and a kernel size 11 × 1.Finally, the Spec3 pipeline combined the Stage 2 exposures into data cubes.

Astrometric Corrections
To ensure that the regions from which we extract spectra are the same between the NIRSpec and MIRI cubes, we calibrate the astrometry to align with Gaia stars, employing the simultaneous MIRI imaging data.We find an offset of 0 2 between the native MIRI astrometry and the detected Gaia stars in the parallel F560W image.After applying this correction to the MIRI MRS cubes, we collapse both the MRS and the NIRSpec cubes in the (observer-frame) 5.0-5.1 μm spectral region, and shift the NIRSpec astrometry such that the two collapsed images match.The correction required for the NIRSpec cube is found to be 0 3.

Spectral Extraction
The nuclear region of VV 114E contains a large number of mid-IR luminous cores (Evans et al. 2022).One of these is most likely an obscured AGN, as suggested by Rich et al. (2023).Two other prominent mid-IR peaks display strong, most likely thermal radio emission at 33 GHz (Evans et al. 2022;Song et al. 2022), and therefore represent regions of intense star formation.Together, these three cores dominate the 7.7 μm emission from the VV 114E nucleus, and they form the subject of study of the present paper.For easy reference, we will henceforth refer to these as AGN, SF-NE, and SF-SW, and their positions are indicated in Figure 1.These cores correspond to the "SW core," "NE core," and "SW-S knot" respectively, in the notation of Evans et al. (2022, their Figure 3).In the notation of Rich et al. (2023), they correspond closely to apertures c, a, and d respectively; they also align with the HCO + (4-3) emission peaks denoted by E2, E0, and E1 respectively by Iono et al. (2013).
We extract spectra from apertures of 0 322 (130 pc) in radius, centered on the apparent intensity peaks in the 5.0 μm slice of the NIRSpec data cube.The position and sizes of the extraction apertures are indicated in Figure 1.Although the JWST Science Calibration Pipeline contains a 1D rfc for extracted spectra, in addition to the rfc performed in the Spec2 step (see Section 2.2), we do not apply this 1D correction for the majority of this work, as this defringer can partially remove the rovibrational lines of interest due to their periodic nature.We only apply it for the spectral fits for C 2 H 2 and HCN (see Section 4.3).No aperture corrections are applied as our analysis will be based on normalized spectra.
To correct for residual systematic errors in the wavelength solution of the instruments, we manually calibrate the redshift of the extracted spectra.We use H 2 emission lines to set the systemic velocity, as these lines trace the bulk of the molecular gas in the nuclear regions.Since the MIRI MRS wavelength solutions are specific to the grating setting used, we only consider H 2 lines close in wavelength to the molecular features under consideration, and we do so separately for each species studied.We do not, however, differentiate between the spatial regions, as the difference in radial velocity between regions is smaller than the difference in radial velocity between different H 2 lines.

Results
We present the 1D spectra extracted from the three apertures over the full wavelength range of the instruments in Figure 2 (NIRSpec) and Figure 3 (MIRI MRS).Here, we briefly discuss the continuum shape and the most prominent features to contextualize our analysis.
The overall view of the spectra immediately reveals several striking differences between the AGN position on the one hand and the SF-NE and SF-SW knots on the other hand.The AGN spectrum shows a steeply rising continuum from 1 to 4 μm, which remains relatively flat at longer wavelengths.This behavior is characteristic of dust-embedded AGNs, and typically attributed to hot dust emission directly associated with the circumnuclear torus (e.g., Barvainis 1987;Pier & Krolik 1992;Hönig & Kishimoto 2010;González-Martín et al. 2019).The SF-NE and SF-SW spectra, on the other hand, exhibit very similar continuum shapes, which are flatter than the AGN spectrum between 1 and 4 μm but have a stronger rise at wavelengths above about 18 μm.The 9.7 μm silicate absorption feature is very deep in the two starburst regions, indicating deeply dust-embedded star formation; the silicate feature is considerably shallower at the AGN position.Furthermore, the PAH emission features at 3.3, 6.2, 7.7, and 11.3 μm are strong toward the SF-NE and SF-SW regions, but nearly vanish in the AGN spectrum (see also Rich et al. 2023).These general properties, which are in agreement with the ideas of Laurent et al. (2000) and Marshall et al. (2018), all suggest the presence of a moderately obscured AGN at the position "AGN," and more deeply embedded star formation at positions SF-NE and SF-SW.
Several broad ice absorption features are detected as well, toward all three regions under consideration.At 3 μm, the broad H 2 O stretch (e.g., Boogert et al. 2015) is observed.The CO 2 stretch is seen at 4.27 μm; interestingly, it is weakest in the most deeply obscured star-forming region (SF-NE).None of the spectra show considerable CO ice absorption at 4.7 μm,  but all show strong CO gas-phase absorption at this wavelength.

Fundamental CO Band
The fundamental band of CO, arising from absorption from the vibrational ground state (v = 0) to the first excited state (v = 1), appears prominently as a sequence of deep lines in all three regions.These bands are shown in detail in Figure 4.The CO bands are comprised of an R-branch bluewards of 4.67 μm, corresponding to transitions with ΔJ = +1 in absorption, and a P-branch on the red side, where transitions satisfy ΔJ = −1.
Here, ΔJ = J u − J l , where J u and J l are the rotational quantum numbers of the upper and lower level, respectively.Lines are labeled by their branch and the rotational quantum number of the lower level involved in the transition, e.g. the transition between v = 0, J = 0 and v = 1, J = 1 is denoted by R(0).At the SF-NE and SF-SW positions, we observe broad, deep bands, indicating highly excited CO gas against a bright background continuum.The AGN region, however, only exhibits a narrow absorption band, with some weak emission present in the form of P-Cygni profiles in the P-branch.
The R-branch is contaminated by two weak emission lines, partially filling in the R(5) and R(24) lines.Furthermore, the Pf-β emission line at 4.654 μm calls for some care in using the R(0) and R(1) lines.The P-branch suffers from more contaminants: the H 2 0-0 S(9) emission line lies at 4.695 μm, but the fundamental bands of 13 CO and C 18 O run through the 12 CO P-branch as well.

H 2 O
We detect rovibrational lines of water vapor between approximately 5.5 and 7.2 μm toward all three regions under consideration; they are shown in detail in Figure 5.The H 2 O lines are weak and narrow toward the AGN and more clearly visible toward the SF-NE and SF-SW regions, where they are significantly broader and deeper.
This spectral region also encompasses the strong 6.2 μm PAH feature, as well as prominent absorptions at 6.85 and 7.25 μm.We detect the latter only in the deeply embedded SF-NE core.These features are typically attributed to aliphatic or hydrogenated amorphous hydrocarbons (Spoon et al. 2001(Spoon et al. , 2004;;Rich et al. 2023); for further discussion of these assignments, see Boogert et al. (2015).

C 2 H 2 and HCN
At longer wavelengths, our spectral and spatial resolutions decrease, impeding the detection of shallow individual lines.Figure 6 shows the three spectra between 13.6 and 14.2 μm.As reported previously by Rich et al. (2023), we detect the blended Q-branches of both C 2 H 2 and HCN at 13.7 and 14.04 μm at the SF-NE position.They may be present in the SF-SW spectrum as well, but the signal-to-noise ratio (S/N) here is too low to confirm a detection.
The abundances of both HCN and C 2 H 2 can be greatly enhanced by high-temperature chemistry in the inner envelopes around young massive stars (e.g., Doty et al. 2002;Rodgers & Charnley 2003).HCN is of particular interest as its rotational transitions in the ν 2 = 1 state have been detected toward several U/LIRGs (e.g., Sakamoto et al. 2010;Aalto et al. 2015aAalto et al. , 2015b;;Falstad et al. 2021).With MIRI MRS, we observe the absorption due to infrared pumping by the strong 14 μm continuum, causing the excitation to the ν 2 = 1 state.For C 2 H 2 , which has no dipole-allowed rotational transitions, mid-IR observations are a unique detection method.

CO Overtones
In addition to the fundamental CO band from the v = 0 − 1 transitions at 4.7 μm, the first overtones-rovibrational  legend).The lines of the 12 CO fundamental band are labeled in dark blue; some 13 CO lines are indicated in red.Toward the SF-NE and SF-SW regions, we see broad, deep bands, indicating that the molecular gas is highly excited.The AGN position exhibits a narrow 12 CO band; the R(1), R(0), and P(1) lines of 13 CO are also detected here.
transitions for which Δv = ± 2-at 2.3 μm lie in the NIRSpec spectral range as well.This part of the spectrum is shown in Figure 7.We detected these absorption bands toward the SF-NE and SF-SW positions, but not toward the AGN position.The bands display prominent bandheads at the blue side, a shape that is characteristic of the spectra of the atmospheres of cool stars, such as red giants and supergiants.The latter are ubiquitous in starburst regions (e.g., Oliva et al. 1995;Armus et al. 1995).The 2.3 μm CO bands appear a few Myr after the onset of star formation, and remain prominent for more than 10 8 yr (Origlia et al. 1999;Leitherer et al. 1999;Origlia & Oliva 2000).The lack of these CO bands toward the AGN position therefore unambiguously confirms its nature as an AGN.

CO
The CO band spectra presented in Figure 4 display a remarkable dichotomy.The SF-NE and SF-SW positions show broad vibrational bands, revealing high rotational excitation temperatures.In contrast, the narrow vibrational band observed at the AGN position indicates a much lower rotational excitation.CO rotational ladders observed at lower spatial resolution in the (sub)millimeter and far-infrared, however, show exactly the opposite behavior, with the highest rotational excitation observed for AGN, and lower rotational excitation for starbursts (van der Werf et al. 2010;Rosenberg et al. 2015;Lu et al. 2017).As discussed previously, the star formation nature of SF-NE and SF-SW, and the AGN nature of the so labeled position are supported by strong evidence from infrared and radio observations (Evans et al. 2022;Rich et al. 2023).The broad CO bands at the two SF positions and the narrow band at the AGN position are therefore very surprising.
Determining the cause of this behavior is the focus of the present paper.To proceed, we first analyze the various molecular bands using rotational population diagrams, in order to quantify the excitation.We then discuss the physical origin of these results, the underlying physics, and the implications in Section 5.
We begin by constructing rotational population diagrams based on our three fundamental CO band spectra.The principle has been discussed in the literature before, but usually in treatments tailored to rotational line emission (e.g., Goldsmith & Langer 1999).Since here we will be deriving rotational excitation temperatures from vibrational absorption lines, we provide the relevant equations, and details of our methods, in Appendix A.1; a summary is provided here.
We treat the spectra as pure absorption and convert them to optical depth spectra through I obs = I cont e − τ , using a spline fit to spectral regions that appear to be line-free to estimate the local continuum I cont .We then fit a Gaussian profile to each line in optical depth space and take the integral of this model to obtain an integrated optical depth.Using Einstein A values from HITRAN23 (Rothman et al. 2013; with values for CO from Guelachvili et al. 1983;Farrenq et al. 1991;Goorvitch 1994), we calculate the lower-level column density for each detected line using Equation (A7).
The distribution of these vibrational ground state level populations allows us to probe the physical conditions of the CO gas by comparing to the Boltzmann distribution.We employ a rotation diagram analysis to characterize the rotational excitation temperature T ex , defined as follows: Here, N i , g i , and E i denote the column density, statistical weight, and energy of level i; k is the Boltzmann constant.In the rotation diagram, we plot the logarithm of N i /g i as a function of the level energy in temperature units.The slope is then governed by the excitation temperature of the gas (see Appendix A.1).
The CO rotation diagrams for the three regions under consideration are shown in Figures 8 (SF-NE and SF-SW) and 9 (AGN).For the rotation diagram analysis, we only use the R-branch lines, as some P-branch lines appear to be affected by emission and other contaminants (Figure 4).We also remove the R(5) and R(24) lines from the fits, as these lines are contaminated by faint emission from other species.
At the SF-NE position, the CO rotation diagram exhibits a low-excitation component in the first 3-5 J-levels, followed by an extended tail of seemingly thermalized gas that continues out to at least J = 23.To quantify the excitation, we fit a model combining two local thermodynamic equilibrium (LTE) components (Equation (A8)) to the measurements.The highexcitation tail is well fit by a rotational temperature T ex ≈ 700 K.
For the SF-SW region, the CO rotation diagram similarly shows a low-excitation component up to J = 4 and an extended high-excitation tail.However, here, the tail has another kink at E l /k B ∼ 800 K, indicating two distinct high-temperature components.Fitting a three-component LTE model, we find rotational temperatures of T ex ≈ 200 K and T ex ≈ 800 K for the J > 5 levels.
At the AGN position, the picture is quite different: the rotation diagram is well fit by two relatively cool components.We do not detect a high-excitation component.Given the sensitivity of these current data, we exclude the presence of a T ex ≈ 700 K component with column density N(CO)  3 × 10 16 cm −2 along the line of sight-2 orders of magnitude lower than the column densities inferred for the SF-NE and SF-SW regions.
Resulting parameters for all fits are displayed in the figures and summarized in Table 1.All positions show a relatively cool (T ex = 15-32 K) gas component dominating the lowest rotational levels.In all likelihood, this component corresponds to the bulk molecular gas in the system (Saito et al. 2015), probably subthermally excited at the SF-NE and SF-SW positions.This gas is probably "foreground" and not directly associated with the nuclear mid-IR cores.The fitting results confirm the presence of highly excited CO gas in two presumed star-forming regions and absence thereof at the AGN position.We will return to this important issue in Section 5.2.
We extract the velocity dispersion and radial velocity with respect to the velocity of warm H 2 from the optical depth line profile fit for each of the lines.We remove the contaminated R(5) and R(24) lines from this analysis.The measured velocity dispersions are corrected for the instrumental broadening, assuming the nominal resolving power of R = 2700 (Jakobsen et al. 2022).The results are shown in Figure 10.
Figure 10 shows that the CO gas observed toward SF-NE is effectively at the velocity of the H 2 lines throughout the band.Toward the AGN, however, the CO exhibits a clear blueshift of V rad ≈ −60 km s −1 .At the SF-SW position, the radial velocity transitions from V rad ≈ −50 km s −1 at the low-excitation CO lines to V rad ≈ −140 km s −1 at the highest J-lines, suggesting the presence of slightly blueshifted cool gas and more rapidly outflowing highly excited gas.
Toward the star-forming regions, the velocity dispersions reveal a distinct lower-dispersion component that corresponds  to the cold component seen in the rotation diagrams (Figures 8  and 9).The SF-SW region exhibits a peak in velocity dispersion of σ V ≈ 100 km s −1 at the R(4)-R(7) lines-where the rotation diagram transitions from the foreground component to the T ex ≈ 230 K component, and the radial velocities start to transition to V rad ≈ −140 km s −1 .The peak in measured velocity dispersion is therefore explained by the blending of two lines arising from CO gas at two distinct velocities, with velocity dispersions of ∼40 km s −1 for the individual kinematic components.These estimates are roughly consistent with dispersions found for rotational emission lines (Saito et al. 2015(Saito et al. , 2018)).Although the rotation diagram is indicative of three temperature components, we find no evidence of a third velocity component.
We model the three CO spectra using the column densities and excitation temperatures of the various components inferred from the rotation diagram, and the velocity dispersions of the lines dominated by specific components, as collected in Table 1.These models are shown together with the continuum-extracted observed spectra in Figure 11.The model spectra of the individual components are shifted by the measured radial velocities listed in Table 1.
For the star-forming regions, the correspondence between data and model is good in the R-branch-from which the rotation diagram was constructed-but less so in the P-branch.Toward the SF-NE core, the observed P-branch absorption lines are considerably shallower than predicted by the model.This asymmetry has significant implications, which we shall discuss in Section 5.1.We note that imperfections in the SF-SW model are most likely due to the blending of at least two kinematic components that were not explicitly separated in the construction of the rotation diagram.
In the AGN region, the observed absorption lines match the model well in both branches, but they are blueshifted with a radial velocity of −55 km s −1 .Most notably, the P(4), P(5), P(6), P(7), and P(8) lines clearly appear in emission as well, at the systemic velocity.These P-Cygni profiles suggest the presence of a molecular outflow.At the SF-SW position, just SE of the AGN, the high-excitation lines are blueshifted at a radial velocity up to −140 km s −1 , while the low-excitation lines are shifted by only −50 km s −1 .
The above analysis makes the implicit assumption that the absorbing gas fully covers the background continuum source.If this is not the case, our calculation underestimates optical depths and column densities.The derived temperatures are, however, not strongly affected if the peak optical depths of the high-excitation lines remain small compared to 1, since in this case the column density ratios remain unaltered.High optical depths would lead to a situation where the strongest lines saturate and have approximately the same depth.An inspection of Figures 4 and 11 reveals that this is clearly not the case at the AGN position and at position SF-SW.However, at position SF-NE, a number of the deepest lines indeed have similar depths, and the possibility of saturation must be considered.
To assess the effect that incomplete coverage of the background continuum might have on the inferred parameters for the SF-NE region, we assume the minimum covering factor of f c = 0.4-a lower covering factor would not be able to produce the apparent absorption depth of 40% that the deepest line observed here reaches-and correct the observed absorption lines for this covering factor.We then construct a rotation diagram from this corrected spectrum.The implied column densities increase by a factor of approximately 3.5 compared to the case where we assume f c = 1; the inferred temperature of the high-excitation component decreases by 19% to T ex = 559 ± 16 K.A similar exercise for the SF-SW region with covering factor f c = 0.7 leads to a column density increase of a factor 1.5 and a temperature decrease of 6% to T ex = 717 ± 48 K. Thus, even in the most extreme possible  case, the high-excitation temperature measured is not strongly diminished due to saturation of the lower −J lines.
The CO column density corresponds to a 4.5 μm dust optical depth τ 4.5 , which can be calculated from the following: Using a standard molecular cloud CO abundance [CO]/ [H 2 ] = 10 −4 , the Milky Way infrared extinction curve determined by Schlafly et al. (2016), and A V /N(H) = 5 × 10 −22 mag cm 2 /H (Bohlin et al. 1978;Rachford et al. 2009), we find that τ 4.5 = 1 corresponds to N(CO) = 3.3 × 10 18 cm 2 .The dust optical depths implied by the CO column densities in Table 1 are then τ 4.5 = 0.5 for SF-NE, and τ 4.5 = 0.9 for SF-SW.Applying the proposed increase of the column density by a factor 3.5 in SF-NE would result in τ 4.5 ≈ 1.8 for that region.
This number exceeds the value derived from the depth of the 9.7 μm silicate absorption feature (Donnan et al. 2023, their Figures 4 and A1) by less than a factor of 2. With the corrected column density for SF-NE, the relative depths of the 9.7 μm silicate feature for SF-NE and SF-SW are also reproduced.

H 2 O
For water, constructing a rotation diagram is more difficult than for CO as many apparent lines are actually blends of several lines.In particular, for the SF-SW region, where the H 2 O lines are considerably broader, many lines had to be removed from the fit.Modeling these additional lines would require a more elaborate procedure (see González-Alfonso et al. 2024), which is not necessary for the purposes of the present paper.The R-branch H 2 O rotation diagrams for the two starburst regions are presented in Figure 12.We fit a singlecomponent LTE model to the data.We do not attempt a fit for the weak H 2 O lines at the AGN position.
For the lines selected for the rotation diagram analysis, we measure the radial velocities and velocity dispersions as described in Section 4.1, now assuming an instrumental resolving power R = 3400 (Labiano et al. 2021).The results, summarized in Table 1, again show gas observed toward the SF-NE core that is effectively at systemic velocity, while the H 2 O absorption at the SF-SW position is considerably blueshifted.In this region, the radial velocity of the water lines matches that of the highly excited CO lines, with V rad ≈ − 110 km s −1 .Interestingly, in both cores, the velocity dispersion of the H 2 O lines is significantly higher than those of CO.
Using the rotation diagram results obtained from the R-branch lines, and the median velocity dispersions, we model the water spectra for the SF-NE and SF-SW regions with a single-component LTE model; they are presented together with the observed spectra in Figure 13.At both positions, the match is reasonable in the R-branch, but we see considerably shallower absorption than predicted in the P-branch.The same P-R branch asymmetry was noted for CO (Section 4.1).We  discuss the origin of this asymmetry between branches in Section 5.1.

C 2 H 2 and HCN
For the strong Q-branches of C 2 H 2 and HCN at ∼14 μm, we cannot use the rotation diagram method because we do not resolve the individual lines.Instead, we use a direct spectral fitting analysis, employing a grid of single-component LTE models and χ 2 minimization.For both molecules, we use a 30 × 50 grid with column density spanning a range Î and temperature spanning T/ K ä (50, 500).The column densities are logarithmically distributed.We fix the velocity dispersion of the individual lines at σ V = 50 km s −1 -the value we found for the water lines-but through experimentation, we find that this value has little effect on the inferred parameters.We limit the spectral range used in the fits to the fundamental Q-branch, and we only model the strong C 2 H 2 and HCN lines at the SF-NE position.
To account for the noise in the spectrum, we apply a bootstrap procedure: we resample the spectral pixels used in the fit 5000 times, and take the best-fit model from each iteration.The resulting corner plots are shown in Figure 14.To visualize the spread in the best-fit models, we draw 50 random bootstrap samples and plot the corresponding models along with the observed spectrum-these fits are shown in Figure 15.Fit results are again displayed in the figures and in Table 1.

Excitation Mechanism
The rotation diagram for CO in the SF-NE core (Figure 8) shows a well-constrained straight line for J  7, implying that the gas is in equilibrium at a rotational temperature of 700 K. Remarkably, the gas appears to be fully thermalized out to at least J = 23.These extremely high levels of thermalization put very strong requirements on local densities if the levels are collisionally excited.We, therefore, need to carefully consider the excitation mechanism at play.In this discussion, we focus on CO first, because for this molecule we have characterized the level populations in the most detail.

Collisional Excitation
Collisional LTE occurs if collisional deexcitation dominates over radiative downward transitions, i.e., if the density of the main collisional partner species-typically H 2 for neutral molecules-is much greater than the relevant critical density.
To establish whether the CO rotation diagrams in Figure 8 can be produced by collisional excitation, we estimate the critical density required to thermalize the highest levels that follow the straight line in the rotation diagram.For the SF-NE core, the levels appear to be thermalized up to at least J = 23 (E/k ≈ 1500 K)-and likely up to J = 27-at a best-fit rotational temperature of T rot ≈ 700 K.At the SF-SW position, the highest-excitation CO component displays LTE conditions at a temperature of T rot ≈ 800 K up to at least J = 30 (E/k ≈ 2600 K).
For the rovibrational transitions, the critical densities are extremely high, as the Einstein A coefficients are of the order A CO,rovib ∼ 10 s −1 , while the collisional deexcitation rates are at most k CO,rovib ∼ 10 −17 cm −3 s −1 .Thus, if collisions dominate the excitation of CO, it must be through the pure rotational transitions.Therefore, we consider the critical densities of rotational transitions in the vibrational ground state.
For the J = 24 − 23 transition of CO, the Einstein A coefficient is A 24,23 = 1.281 × 10 −3 s −1 , while the sum of the collisional deexcitation rates amounts to ∼5 × 10 −10 cm 3 s −1 .Thus, in the absence of a radiation field, the critical density is 3 .In order for the gas to be thermalized up to this level, as the rotation diagram suggests, we, therefore, need densities n(H 2 )  10 7 cm −3 .We combined non-LTE modeling using DESPOTIC (Krumholz 2014) with Markov Chain Monte Carlo fitting using emcee (Foreman- Mackey et al. 2013) to test whether the level populations that are not dominated by foreground gas (J 3) in Figure 8 could be produced in the absence of a radiation field.In these models, statistical equilibrium level populations were calculated, considering rotational transitions with Einstein A coefficients and collisional rates taken from LAMDA (Schöier et al. 2005;Yang et al. 2010) and using the escape probability formalism assuming a spherical geometry.The resulting posteriors suggest a similar density requirement of n(H 2 )  6 × 10 6 cm −3 .
Although this density is not impossible on small scales such as the central regions of actively star-forming cores, such regions are expected to be subparsec in size.In the present case, however, the thermalized high-temperature components dominate the column densities over R ≈ 130 pc apertures.The absorption signal arises from a likely smaller region set by the projected area of the background continuum source, but radio observations show structures on the scale of several parsecs.Song et al. (2022) estimate deconvolved radii of ≈40 pc for SF-NE and SF-SW from their 33 GHz VLA observations.ALMA observations of the 350 GHz continuum also reveal structures that are considerably larger than the 0 2 (80 pc) beam (Saito et al. 2018).Iono et al. (2013) used HCN(4-3)/(1-0) and HCO + (4-3)/(1-0) line ratios to estimate densities of n (H 2 ) = 10 5 − 10 6 cm −3 , which is still insufficient to thermalize the CO out to J = 23.
Therefore, it is highly unlikely that the CO we observe toward the SF-NE and SF-SW cores is actually in collisional LTE, unless the 4.7 μm continuum totally originates in small ultradense (n(H 2 ) > 10 7 cm 3 ) regions.However, even then, collisional excitation will be irrelevant, since at such high densities the gas temperature will equal the dust temperature, which is set by the local radiation field.Mapping of the [Fe II] 5.34 μm emission line has revealed the presence of a large-scale bow shock covering the SF-SW region (Donnan et al. 2023).It is of interest to consider the possible role of this shock for the excitation of the molecular gas in SF-SW, since C-type shocks in molecular clouds can easily produce postshock temperatures of several times 10 3 K without dissociating the molecules (Draine et al. 1983).In the present case, however, the highly excited molecular gas observed in absorption cannot be identified with this postshock gas, as the H 2 O and high-J CO lines are blueshifted by >100 km s −1 with respect to the [Fe II] line.Furthermore, the detection of prominent CO 2 ice absorption features (Figure 2) shows that no significant grain mantle destruction, which always occurs in such shocks, has taken place.In addition, shocks produce a range of temperatures in the cooling postshock gas, while the observed CO lines are well represented by a single high-temperature component for 19 J 31 (Figure 8).We, therefore, conclude that the large-scale bow shock traced by [Fe II] is not responsible for the high CO excitation in SF-SW, and that the role of more localized shocks is likely minor at most.Donnan et al. (2023) also found strong [Fe II] emission toward the SF-NE region.Here, the molecular lines are all at rest with the [Fe II] line, so an association between the two cannot be ruled out.However, strong absorption by CO 2 ice is clearly present, and in this region, the CO gas appears to be completely dominated by a single excitation temperature for 7 J 27 (Figure 8).For these reasons, shock excitation is unlikely here.

Far-IR Radiative Excitation
We next consider radiative excitation by the far-IR radiation field.This model is suggested by the fact that the relevant CO pure rotational lines have frequencies between 1037 GHz (289 μm) for J = 9-8-where the 700 K component begins to dominate-and 3098 GHz (97 μm) for J = 27 − 26.At these far-IR wavelengths, U/LIRGs have intense local radiation fields, and can become optically thick (resulting in a blackbody local radiation field) down to wavelengths of typically ∼100 μm (Solomon et al. 1997).Since radiative excitation of the CO rotational lines has hardly been discussed in the literature, we summarize the relevant equations in Appendix A.2.
To equilibrate all the rotational levels probed, the radiation field must have the shape of a blackbody of at least the excitation temperature of the gas, across the indicated wavelength range.A blackbody of T = 700 K peaks at 4.1 μm; one of T = 800 K peaks at 3.6 μm.However, the observed spectra of SF-NE and SF-SW peak at much longer wavelengths, as an inspection of Figure 3 immediately confirms.Furthermore, at the column densities (Saito et al. 2015) and local extinctions (Donnan et al. 2023) estimated for SF-NE and SF-SW, the dust will be optically thin at the wavelengths of the relevant CO lines.
We conclude that radiative excitation by the far-IR radiation field must be ruled out as an excitation mechanism.

Mid-IR Pumping
Finally, we consider excitation by mid-IR radiation, through the vibrational transitions at approximately 4.7 μm.This mechanism is attractive since, as noted in Section 5.1.2,the ∼700 K radiation field required to produce the observed level of radiative excitation, if interpreted as a blackbody, peaks at ∼4 μm, so a strong radiation field is available.A close inspection of the CO band spectra presented in Figures 4 and 11 reveals the presence of weak vibrational emission lines, in addition to the much stronger absorption lines, at several positions.As noted in Section 5.1.1,critical densities for the vibrational transitions are extremely high, so these emission lines cannot result from collisional excitation.Their presence therefore directly shows that radiative excitation through the 4.7 μm band plays a role.An additional argument supporting the importance of mid-IR pumping is provided by the pronounced observed asymmetry between the P-and R-branches of both the CO and the H 2 O lines, which was noted in Sections 4.1 and 4.2.In a situation where only absorption occurs, this asymmetry is impossible to explain, since the ratio of P-and R-branch lines originating from the same lower level is determined solely by the relevant Einstein A coefficients.In the presence of radiatively excited emission, this situation changes.Consider a 4.7 μm continuum photon that is absorbed in the n g n g J J J J 2 2 , absorption in the R( ¢ J ) line will be stronger than that in the P( ¢ + J 2) line.This mechanism, which was analyzed by González-Alfonso et al. (2002), is illustrated in Figure 16.The stronger absorption in the R-branch line compared to the P-branch, while the emission lines have similar strengths, leads to the observed asymmetry between P-and R-branches.This process also results in a net population transfer to higher rotational levels, and may therefore affect the rotational temperature in the vibrational ground state (see also Carroll & Goldsmith 1981;Sakamoto et al. 2010).
Depending on the relative strength of emission and absorption, the asymmetry can lead to a situation where both branches are in absorption (as in the cases studied in this paper), or both are in emission, or where the R-branch is in absorption, but the P-branch appears in emission (as illustrated in Figure 16).An example of the latter situation is given by Pereira-Santaella et al. (2024, their Figure 3, middle panel) in the nearby LIRG NGC 3256.We have explored our data cube, and likewise find extended regions where the R-branch is seen in absorption, and the P-branch is seen in emission, confirming the importance of mid-IR radiative excitation.
If mid-infrared pumping indeed dominates the level populations in the vibrational ground state, the implications for the interpretation of the rotational temperature are profound.Since the rovibrational lines, unlike the pure rotational lines, all lie at effectively the same wavelength, a blackbody is not required to produce equilibrium level populations.Instead, the excitation temperature traces the almost monochromatic brightness temperature of the exciting radiation field at the wavelength of the transitions.Thus, the excitation temperature found for CO can be interpreted as a measure of the intensity of the local radiation field at 4.7 μm, and is in fact its (Planck) brightness temperature (see Appendix A.2).
The above discussion assumes that the mid-IR radiation field not only excites the vibrational levels but also determines the high-J rotational levels in the ground vibrational state.If the excitation is purely by radiation, this is guaranteed to be true, since in this case radiative equilibrium at the excitation temperature T rad applies.However, collisions, if sufficiently frequent, can modify the populations of the rotational states.This situation was found in the extended outflow regions in NGC 3256 studied by Pereira-Santaella et al. (2024).In contrast, here, we are studying the most intense luminous cores, where local radiation fields are much more intense.We evaluate the relative importance of radiative and collisional excitation of the rotational levels in Appendix A.3.This calculation shows that, in the intense mid-IR radiation fields encountered in the regions studied here, radiative pumping through the mid-IR vibrational levels also dominates the population of the rotational levels.Only for the lowest rotational levels, where collisional rates are significantly higher, collisional excitation may play a role, and this may account for the low-excitation temperatures derived at levels well below J = 10 (see Figures 8 and 9   Figure 16.Schematic representation of the rovibrational excitation process in the case of a P-and R-branch.The density of each level i is denoted by n i ; A denotes the Einstein A coefficient, which does not vary strongly between adjacent rovibrational lines.As long as there is no population inversion between levels 1 and 2 (i.e., n 1 /g 1 > n 2 /g 2 ), the R-branch absorption is stronger than the P-branch absorption, while the corresponding emission lines -both proportional to An 3 -are equally strong.The net effect is an asymmetry between the appearance of the P-and R-branches.
Turning now to the other molecules, the same branch asymmetry that we noted for CO was found for H 2 O (see Section 4.2), confirming that radiative excitation plays a role.Furthermore, the quality of the LTE fit (see Figure 13) suggests equilibrium conditions, despite the high rotational critical densities of n crit ∼ 10 8 cm −3 for this molecule.Thus, the H 2 O observed toward the SF-NE core is likely radiatively excited as well.
For HCN and C 2 H 2 , we can only analyze the Q-branch, and the S/N of our data is not sufficient to conclude whether a branch asymmetry is present in these cases or not.Due to the limited S/N and spectral resolution, it is difficult to assess whether non-LTE effects play a role.We note, however, that HCN has high rotational critical densities (n crit ∼ 10 7 -10 8 cm −3 ) and is strongly susceptible to both far-IR radiative excitation and mid-IR pumping (Sakamoto et al. 2010).Given the high 14 μm intensity observed toward the SF-NE core (Figure 3), it is likely that HCN is excited through the mid-IR radiation field as well.For C 2 H 2 , no critical densities are available, and we cannot with certainty establish the dominant excitation mechanism.
We note that the wavelengths of the molecular bands all lie close to the peak of a blackbody radiator at the excitation temperature derived for that band.This result is natural in the case of radiative excitation through a molecular band, since this mechanism is most energy-efficient if the exciting radiation field peaks close to the wavelength of that band.In the case of collisional excitation, however, the lowest-excitation temperatures would be expected for H 2 O, which has the highest critical densities.The observed trend in excitation temperatures thus provides further support for radiative excitation through the vibrational bands.

AGN versus Starburst
We have established that the observed molecular bands are dominated by radiative excitation.This conclusion provides an explanation for the conundrum described in Section 5.1: why do the molecular bands display much higher excitation in the star-forming regions than at the AGN position?
The crucial observation is that the radiation field that provides the excitation is the angle-averaged radiation field (see Equation (A9) in Appendix A.2).This leads to a pronounced difference in the excitation between the starforming regions and the AGN environment, which is illustrated in Figure 17.In the star-forming regions, the exciting radiation field is provided by hot dust, heated by the young stars that form throughout the regions (plus possible contribution from starlight).This geometry is found for instance in the eponymous hot core in the Orion-KL region, which is externally heated by the surrounding young star cluster (Blake et al. 1996;Zapata et al. 2011).Since the gas is mixed with the hot dust, there is no geometric dilution of the exciting radiation field.The core of the star-forming region will produce the most intense radiation, allowing gas in front of it-which is the gas probed by our data-to be observed in absorption.Hence, the local mid-IR radiation field excites the molecular gas, and the more intense radiation from deeper layers of the core provides a background continuum.
The situation is different in the AGN case: the 4.7 μm radiation field is now provided by the (parsec-size) circumnuclear torus, which serves both as a background source and provides the exciting radiation field.The solid angle from which the absorbing cloud is excited now depends on the size of the torus and the distance of the cloud from the torus, but is in any case much smaller than 4π sr.As a result, the exciting radiation field is strongly geometrically diluted at the location of the absorbing gas.The result is a strongly reduced angleaveraged exciting radiation field, and therefore, the CO does not reach the high degree of excitation observed toward the intense star-forming regions.
As the CO line detections toward the AGN only extend to J = 8, for which the critical densities are only n crit ∼ 10 5 cm −3 , and the rotational temperature is relatively low, we cannot definitively establish whether the T ex ≈ 200 K component is dominated by collisions or radiation.We note, however, that the detection of rovibrational CO emission lines confirms that a substantial amount of CO near the AGN is infrared-pumped to the v = 1 state, and therefore, the pumping rate must be considerable.
Intense local infrared radiation fields can also be produced very close to single, massive, dust-embedded stars.CO located in the centrally heated dusty envelope will be exposed to a For the AGN region, molecular gas in a fiducial shell (orange) surrounding the central AGN + torus (black) is irradiated from only a small solid angle toward the AGN.Therefore, the molecular gas will not be highly excited, despite the high intensity of the continuum emission from the AGN.In the case of a star-forming region, the gas, dust, and starlight are effectively mixed, and therefore, the sketched patch of molecular gas is irradiated isotropically by infrared continuum emission.
strong locally generated 4.7 μm radiation field from every direction, if the region is optically thick at that wavelength.While this appears to be the case for SF-NE and SF-SW, we note that the AGN environment must have moderate 4.7 μm optical depth, since the AGN is detected already at 2.3 μm, and is not associated with an extinction peak (Donnan et al. 2023, their Figure 4).Any CO in the AGN envelope will then be exposed to a much weaker 4.7 μm radiation field than near the more deeply embedded stars, since the envelope is now optically thin.The physical principle of this model is the same as that above, and only the adopted geometry is different.
We finally note that the column density implied by the CO data toward the AGN is more than an order of magnitude lower than typical column densities of the geometrically thick AGNobscuring circumnuclear torus (N(H 2 ) ∼ 10 23 cm 2 , e.g., Hönig & Kishimoto 2010).Such a torus could only produce the low column densities inferred if it were seen nearly face-on, which would result in a type 1 AGN, contrary to what is observed.It is therefore not possible to associate this gas with the torus itself.

Local Conditions in the Star-forming Regions
For SF-NE, we have derived a maximum possible CO column density N(CO) = 6 × 10 18 cm 2 .The local gas density as derived from CO at this position is about 10 5 cm 3 , but considering also HCN and HCO + yields = n log 5.0 5.9 H 2 (Saito et al. 2015).Combining these estimates gives a typical path length through the high-excitation CO component of at most 0.06 pc.This small typical dimension is not surprising, since the required exciting radiation field can only be generated at the cores of the densest star clusters or close to very bright stars.The emerging picture thus has luminous dense clusters of young stars distributed throughout the star-forming regions SF-NE and SF-SW, with the high-excitation CO located in regions of high radiation intensity within such clusters, as illustrated in Figure 17.
In order to estimate whether the required luminosity densities can be achieved, we adopt parameters of well-studied super star clusters.The densest young massive star cluster in the Milky Way is the Arches cluster with a central density of approximately 2 × 10 5 M e pc −3 (Espinoza et al. 2009).Rapid dynamical evolution will lead to even higher central densities (Harfst et al. 2010).A central density as high as 10 7 M e pc −3 has been proposed for the R136 cluster in the Large Magellanic Cloud, although this number is not undisputed (Selman & Melnick 2013).In any case, it seems not unreasonable to suppose that the massive young star clusters in SF-NE and SF-SW can reach central densities as high as about 10 6 M e pc −3 .Using Starburst99 (Leitherer et al. 1999), we find that an instantaneous starburst with a Salpeter initial mass function (IMF) cutoff at 100M e has a bolometric luminosity-to-mass ratio (L/M) of about 1.7 × 10 3 L e /M e during the first few million years.Using simple energy density arguments, we can estimate the typical radiation intensity within these clusters and find that it falls short of the required value by about an order of magnitude.
However, this estimate is probably overly conservative, for several reasons.First, the star formation rates of SF-NE and SF-SW are much higher than the most extreme values found in Local Group galaxies; therefore, the clusters in SF-NE and SF-SW are likely significantly more extreme in mass, central density, and stellar content than the densest young clusters in the Local Group, in particular if they undergo rapid dynamical evolution.Second, the IMF of the Starburst99 models used above does not consider stars more massive than 100M e , while in reality the IMF in these regions is expected to be fully populated, extending well above 100M e .Finally, it has been proposed that in the most actively star-forming regions of (U) LIRGs the IMF may be top heavy (Romano et al. 2017;Sliwa et al. 2017;Zhang et al. 2018).All of these effects would lead to a significantly more intense local radiation field than in the most extreme star clusters in Local Group galaxies.Indeed, "Super Hot Cores" associated with the youngest super star clusters have been found in the nearby starburst galaxy NGC 253 (Rico-Villas et al. 2020).
Nevertheless, it is clear that the required luminosity densities are extreme.We, therefore, also investigate the alternative model, where the molecular gas is simply located close to luminous young stars.Using the models for dust-embedded stars by Scoville & Kwan (1976), we find that a 120M e star, with a luminosity of 2 × 10 6 L e (Ekström et al. 2012), can produce the required radiation intensity out to a distance of about 0.005 pc.Assuming that the high-excitation CO is located in a layer of this thickness then implies a local density ñ 10 cm H 6 3 2 , in excellent agreement with observationally estimated values.Since the most massive stars form in dustembedded ultracompact and hypercompact H II regions, where most of the Lyman continuum is absorbed by dust (Churchwell 2002;Hoare et al. 2007, and references therein), the molecular gas can be shielded from rapid destruction in these harsh environments.
In summary, although the implied local radiation intensities are very high, they can still be produced by stellar sources in the regions of intense star formation, or close to the most luminous young stars, that are expected to exist in LIRGs, and our spectra provide ideal probes of local conditions in such regions.We note that an alternative model for generating the intense local radiation field, involving an accreting intermediate-mass black hole, has been suggested by González-Alfonso et al. (2024).In this case, a compact nonthermal radio point source is expected to be present.Radio VLBI observations will thus be able to distinguish between these two models.

Kinematics and the Role of Shocks
The CO spectrum observed toward the AGN (Figure 11) shows P-Cygni profiles for the P(4)-P(8) lines: emission from infrared-pumped CO gas is observed at the rest-frame wavelength, and blueshifted lines with a radial velocity of V rad ≈ 10 −50 km s −1 with respect to H 2 are seen in absorption.Absorption in these lines traces the "warm" T ex ≈ 200 K component (see Figures 9 and 11) closer to the AGN.Thus, this outflow is directly associated with the nucleus, and these observational signatures can be produced by an expanding shell around the AGN.
Although the blueshifted absorption lines reveal the presence of outflowing molecular gas, they only hold information about the line-of-sight velocity and column density of the gas inside a narrow line of sight against the AGN background continuum.An estimate of the mass outflow rate is precluded by the absence of information on the dimensions and geometry of the outflow.The P-branch emission lines originate from a larger projected area around the AGN.However, due to the partial blending with the absorption lines, we cannot reliably fit their line profiles, and therefore, we cannot robustly estimate the corresponding column densities or velocity dispersions.Thus, further characterization of the CO outflow is not possible with the present data.
While the excitation conditions of the observed molecules in the SF-NE and SF-SW star-forming regions are similar, we find significant differences in the kinematics.Toward the SF-NE region, the absorption lines of CO, H 2 O, C 2 H 2 , and HCN are all effectively at rest with the bulk of the warm molecular gas, as traced by H 2 .The SF-SW region, however, exhibits strong blueshifts and generally broader lines for both CO and H 2 O.The cold foreground CO observed toward this region has a similar radial velocity as that observed toward the nearby AGN; the more highly excited CO reaches a radial velocity of V rad ≈ −150 km s −1 .The H 2 O is strongly blueshifted as well, with V rad ≈ −100 km s −1 .These findings indicate that, in the SF-SW core, the radiatively excited gas close to the young stars is driven out by stellar winds or radiation pressure.

Relation to Galactic Star-forming Regions, ULIRGs, and CONs
Since the highly excited vibrational bands probe the environments of the most massive recently formed stars, it is not surprising that the spectra show many of the features typically found in Galactic high-mass star-forming regions, such as strong absorption features from gas-phase and solidstate species, and a rising continuum toward long wavelengths (e.g., Evans et al. 1991;Boonman et al. 2003;Boonman & van Dishoeck 2003;Barentine & Lacy 2012;Li et al. 2022;Beuther et al. 2023).However, closer inspection reveals a number of interesting differences.The most luminous dust-embedded young star known in the Milky Way is AFGL 2591VLA 3 with a luminosity of about 2.3 × 10 5 L e and a mass of 40M e .Molecular vibrational bands between 4 and 13 μm have been studied by Barr et al. (2020Barr et al. ( , 2022)), who find a fairly constant excitation temperature of 600-700 K for the warm gas component across a number of different molecular species.This result points to thermalized collisional excitation, possibly in a circumstellar disk, and indeed, the estimated density ñ 10 cm H 9 3 2 is much higher than the density estimates in VV 114.
Another interesting comparison can be made with the externally heated hot in the Orion-KL region.It displays CO vibrational bands in absorption toward the object IRc2 (González-Alfonso et al. 2002), which display striking similarity to the bands observed in SF-NE and SF-SW.These and more recent observations (Beuther et al. 2010;Nickerson et al. 2023) show an excitation temperature of 100-200 K and a density of a few times 10 8 cm 3 It, thus, appears that conditions in SF-NE and SF-SW correspond to the highest excitation temperatures in Galactic high-mass star-forming regions but not require the very high local often invoked for these regions.
The average infrared surface luminosity of SF-NE about 10 13 L e / kpc −2 .This number typical that of ULIRGs (Thompson et al.Therefore, we may expect in ULIRGs local similar to those derived in the present paper.Our thus, confirm earlier suggestions that the star formation surface densities of U/LIRGs are similar to those of the most active star formation in the Milky Way, but scaled up to the entire interstellar medium of these galaxies (Solomon et al. 1997;Downes & Solomon 1998;Papadopoulos al. 2012).display a large range of extinctions (e.g., Spoon et al. 2007Spoon et al. , 2022)), and we may expect to be able to the same excitation diagnostics as in the present paper in the low-extinction ULIRGs, but may not be accessible in the more obscured Compact obscured nuclei (CONs) are defined by their luminous rotational HCN emission lines in the ν 2 = 1 vibrational state (Falstad et al. 2021).They contain highly obscured (N(H 2 ) > 10 24 cm 2 ), compact (10-100 pc), and luminous nuclei and are found in about 50% of all ULIRGs and 20% of LIRGs.No ν 2 = 1 HCN emission lines have been detected from VV 114 (Saito et al. 2018), and the column densities derived for VV 114 are much lower than typical CON values.Furthermore, the AGN identified in VV 114 is only moderately obscured.However, the HCN absorption band discussed in the present paper does provide the population of the ν 2 = 1 level.Studies of this band will thus provide a new diagnostic of local conditions in CONs, as long as foreground extinction is not so high that it prevents detection of the region where the ν 2 = 1 lines originate.

Conclusions
In this work, we have made use of JWST NIRSpec and MIRI MRS spectra in R ≈ 130 pc apertures to study the rovibrational absorption spectra of CO, H 2 O, C 2 H 2 , and HCN in the heavily obscured nucleus of VV 114E.Two of the regions studied (SF-NE and SF-SW) have been previously identified as regions of intense star formation, and one region (labeled "AGN" in this work, or "SW core" in previous studies) has been reported as an AGN.We reach the following conclusions: 1.The detection of CO bandhead absorption at 2.3 μm toward regions SF-NE and SF-SW, but not toward the AGN, is further evidence of the stellar origin of the former and AGN nature of the latter (Section 3.4).2. We detect highly excited CO in the star-forming regions, with excitation temperatures as high as T ex ≈ 700 and 800 K. Toward the AGN, we do not detect such a highexcitation component, and find only weakly excited CO with T ex ≈ 200 K (Section 4.1).3. Strong absorption lines of gas-phase H 2 O are detected toward the SF-NE and SF-SW regions, at excitation temperatures of T ex ≈ 300 and 200 K respectively.Toward the AGN, H 2 O lines are detected, but much weaker (Section 4.2). 4. We confirm the detection of C 2 H 2 and HCN toward the SF-NE region, and find excitation temperatures of T ex ≈ 220 and 150 K respectively.These absorption features are potentially seen toward the SF-SW region as well, but at much lower S/N (Section 4.3). 5. We observe an asymmetry between the P-and R-branches of the vibrational bands of both CO and H 2 O, which is characteristic of radiative excitation of the vibrational levels.Toward the AGN, we observe P-Cygni profiles in the CO P(4)-P(8) lines, indicative of a molecular outflow at radial velocity V rad ≈ 50 km s −1 with respect to H 2 lines (Sections 4.2 and 5.4).6.We conclude that mid-infrared radiative pumping can best account for the equilibrium conditions at high rotational temperature found for CO in the star-forming regions (Section 5.1).
7. We conclude that both the P-/R-branch asymmetry and the highly excited CO are indicative of the geometry of a star-forming region, where molecular gas is mixed with the dust that radiatively pumps it to higher vibrational states.For gas surrounding the AGN, the radiation field is geometrically diluted, leading to lower excitation temperatures (Section 5.2). 8.We find that, in the star-forming regions, the CO vibrational bands probe the highest local radiation fields, associated with young super star clusters or the most massive young stars.
In this work, we have demonstrated the power of JWST NIRSpec/MIRI MRS IFU observations of rovibrational bands to characterize the excitation conditions of molecular gas in highly obscured U/LIRG nuclei, as well as the strong effects mid-infrared pumping may have on the molecular gas in intense starburst regions.Future studies of these molecular bands may explore this new avenue to potentially distinguish between pure starburst effects and AGN activity in similarly obscured systems.
the mid-infrared spectra are typically presented on a wavelength axis.In practice, we fit Gaussian profiles to the wavelength-space lines, and use the analytical solution of the integral to determine the column densities.For an amplitude τ 0 and a line width σ τ,λ , we compute the column density as follows: If stimulated emission is actually important, the derived column densities N l can be (approximately) corrected by dividing them by a factor - l e 1 hc k lT u vib .The derived column densities can now be used to construct a rotation diagram, where we plot the log of weighted column densities ( ) N g ln l l as a function of the corresponding level temperature E l /k.If the gas is in LTE at an excitation temperature T ex the levels will a Boltzmann distribution: Here, N tot is the total column density of the species (as far as it is at this excitation temperature), and Z(T) is the partition function.In the rotation diagram, a medium in LTE at excitation temperature T ex will display a straight line with absolute slope proportional to 1/T ex .If there are several distinct excitation temperature components along the line of sight, or if the medium is not in LTE, the measurements will deviate from this simple linear relation.However, one may still infer one or several excitation temperatures to characterize the system.
The rotation diagrams in this paper are constructed from rovibrational lines from levels in the vibrational ground state.Thus, the excitation temperatures found in these rotation diagrams are rotational temperatures, T rot .The vibrational temperature T vib describes the distribution of vibrational levels, which cannot be directly determined from the fundamental band.For absorption data, the presence of hot bands would provide constraints on the vibrational temperature.

A.2. Combined Collisional and Radiative Excitation of a Twolevel System
In order to best illustrate the physical concepts, we focus our description on a two-level system.The concepts and methods are easily generalized to multilevel systems.Our treatment follows largely the methods of D11.
The system can undergo collisional transitions, spontaneous radiative decay, and induced radiative transitions.We parameterize the radiation field through the angle-averaged photon occupation number n γ and the radiation temperature T rad : Here, u ν is the specific energy density of the radiation field.
, where n Ī is the angle-averaged specific intensity, n γ simply expresses n Ī in a dimensionless form.The radiation temperature T rad , which is defined by Equation (A9), represents the temperature that a Planck function would have, in order to generate a radiation field with energy density u ν .It, thus, expresses the radiation field as a (Planck) brightness temperature.
Following D11 (his Equation (17.5)), the upper level (1) is now populated and depopulated as follows: Here, k 01 and k 10 are the collisional excitation and deexcitation coefficients, A 10 is the Einstein A coefficient of the transition, n i and g i are the density and degeneracy of level i, and n c is the density of the dominant collision partner-H 2 in the case of molecular gas.
In statistical equilibrium, the level populations are now as follows: A 1 1 The critical density in the presence of radiation is expressed as In the absence of a radiation field, this expression reduces to the usual ratio of the Einstein A coefficient and the collisional deexcitation coefficient.If there is a strong radiation field, the critical densities are increased due to the addition of stimulated emission; a higher density is needed to compete with the increased rate of radiative decay.
The collisional excitation and deexcitation coefficients are related by Here, E 10 = E 1 − E 0 is the energy difference between the two levels; T kin is the kinetic temperature of the gas.We can use this relation to simplify Equation (A11) to the following: where we have used the fact that the radiation field involved is that at the transition energy hν = E 10 .Equation (A15) provides a completely general expression for the excitation of a twolevel system excited by both collisions and radiation.Which mechanism dominates the excitation is determined in the first place by the local density n c : if n c ? n crit and no strong radiation field is present, collisions dominate, and the excitation temperature will approach T kin .If on the other hand n c = n crit , and a strong radiation field is present, radiative excitation will dominate, and the excitation temperature will approach T rad .
If the populations are governed by the radiation field, we obtain In other words, if the radiation field dominates the level populations, we retrieve a Boltzmann distribution at the radiation temperature, i.e., radiative LTE.The excitation temperature of the gas is then equal to T rad .If a radiationdominated system has several lines at various widely spaced frequencies, a single excitation temperature will be found only if the exciting radiation field is a blackbody over that entire frequency range.

A.3. Population of the Rotational Levels in the Vibrational
Ground State Here, we estimate the relative importance of radiative and collisional excitation of the rotational levels of the ground vibrational state.Denoting the particle density in level j by n j , the rate of collisional population R j,coll (in units of per cubic centimeters per second) is given by = S ( ) R n k n, A 1 7 where k ij is the collisional excitation coefficient for E i < E j , the collisional deexcitation coefficient for E i > E j , and k ii = 0.The radiative population rate for level j can be calculated as where λ is the central wavelength of the vibrational band.
In order to establish whether collisional or radiative excitation dominates the rotational levels, Equations (A17) and (A21) can be evaluated using population ratios from Figures 8 and 9.We illustrate the process by considering the J = 10 level, which is one of the lowest levels that appears to be definitely in the high-excitation components.CO data from ALMA analyzed by Saito et al. (2015) imply T kin = 25-90 K, and = n log 3.5 5.0 H 2 (valid for region R21 in their Table 12, which corresponds most closely to the SF-NE region studied here).We choose T kin = 90 K, in order to calculate the maximum possible collisional population rate.With these parameters, collision coefficients from Yang et al. (2010), and expressing all level populations relative to n 9 (for reasons that will become clear shortly), we find In order to evaluate Equation (A21), we use the fact that A rovib ∼ 10 s −1 for CO (from the HITRAN database; see Section 4.1 for further references).With T rad = 700 K, and λ = 4.67 μm, we find from Equation (A9) that n γ = 0.012, so that stimulated emission can be ignored.We then obtain In this expression, n 11 can easily be expressed in terms of n 9 using the population ratios displayed in Figures 8 and 9, which then allows the radiative and collisional population rates to be compared directly.This calculation shows that, in the strong radiation fields in the star-forming regions considered here, radiative excitation of the rotational levels dominates over collisional excitation of the high-temperature components by several orders of magnitude, even at the highest densities proposed for these regions.

Figure 1 .
Figure 1.MIRI images of VV 114, where N is up, and E is to the left.The field of view of NIRSpec and MIRI channel 1 are indicated in light blue and dark blue respectively.Top: MIRI F770W image of the full system.In the eastern nucleus, saturated pixels appear as white dots.Bottom: MIRI F770W SUB128 image.The apertures of 0 32 in radius considered in this work are indicated by green circles.The association of these regions with an AGN and star-forming regions is motivated by the work of Evans et al. (2022) and Rich et al. (2023; see text).

Figure 2 .
Figure2.NIRSpec spectra toward the three nuclear regions under consideration.The spectral regions close to 1.4, 2.4, and 4.05 μm rest-frame wavelengths are affected by the NIRSpec detector gaps; these parts of the spectrum are masked out in this figure.Several prominent features are indicated.The AGN spectrum rises much more steeply between 1 and 4 μm than the spectra of the star-forming regions.

Figure 3 .
Figure3.MIRI MRS spectra toward the three regions under consideration.Prominent PAH emission features and the deep 9.7 μm silicate absorption band are labeled.The AGN spectrum is relatively flat and exhibits only a moderately deep 9.7 μm silicate absorption band compared to the spectra of the starburst regions.

Figure 4 .
Figure4.The fundamental CO band as observed toward three nuclear regions in VV 114E.The intensities are offset by a factor a for visibility (see legend).The lines of the 12 CO fundamental band are labeled in dark blue; some 13 CO lines are indicated in red.Toward the SF-NE and SF-SW regions, we see broad, deep bands, indicating that the molecular gas is highly excited.The AGN position exhibits a narrow 12 CO band; the R(1), R(0), and P(1) lines of 13 CO are also detected here.

Figure 5 .
Figure 5. Water vapor absorption lines in the MIRI MRS spectra.They are most prominent at the SF-SW position.The location of the H 2 O R-and P-branches are indicated by black arrows.Several emission lines are marked: dashed lines indicate the H 2 S(5), S(6), and S(7) lines, a dotted line indicates the [Ar II] line, and dashed-dotted lines indicate hydrogen recombination lines.The intensities are offset by a factor a for visualization purposes.

Figure 6 .
Figure 6.MIRI MRS spectra at ∼14 μm.The spectra for the different apertures are offset by a term a (see legend) for visibility.No 1D residual fringe correction was applied here.The Q-branches of HCN and C 2 H 2 are clearly detected toward the SF-NE region, and are tentatively detected toward the SF-SW region as well.

Figure 7 .
Figure 7. CO overtone bands at 2.3 μm.From left to right, the dashed lines indicate the v = 0-2 band, the v = 1-3 band, and the v = 2-4 band; higher bands could not be observed due to the detector gap.The CO absorption features are detected toward the SF-NE and SF-SW positions, but not toward the AGN position.

Figure 8 .
Figure 8. Rotation diagrams from the R-branch lines of CO toward the star-forming cores, including a two-component (SF-NE)/three-component (SF-SW) LTE model fit to the points.Gray dotted lines indicate the LTE models for each individual component.In both regions, the R(5) and R(24) lines are partially filled in by weak emission from other species, leading to an underestimation of the CO level column density; the corresponding points are marked in purple on the rotation diagram and were left out of the fit.

Figure 9 .
Figure 9. Rotation diagram from the R-branch lines of CO toward the AGN position.Gray dotted lines indicate the LTE models for each individual component.The J 3 levels are likely dominated by cool foreground gas.

Figure 10 .
Figure 10.Radial velocities (top) and instrument-corrected velocity dispersions (bottom) of the CO R-branch lines as a function of the energy of the lower level involved in the transition.Gray dashed lines indicate the 0 km s −1 line with respect to H 2 emission lines.Toward the AGN, the CO lines are unresolved; these are indicated by triangles in the bottom panel.

11.
The continuum-extracted observed CO spectra vs. the model spectra corresponding to the best-fit two-component (SF-NE, AGN) or threecomponent (SF-SW) LTE model.The model spectra are shifted to match the measured radial velocities.The model-subtracted spectra are presented as well.The AGN spectrum (bottom panel) is shown on a smaller spectral scale to highlight the blueshift of the absorption lines; an enlarged view of the P-branch P-Cygni profiles is included.The weak, broad feature at ≈4.68 μm is due to CO ice.

Figure 12 .
Figure 12.Rotation diagrams of water vapor R-branch lines toward the SF-NE and SF-SW positions.The single-component LTE fit and its inferred excitation temperature are indicated.

Figure 13 .
Figure13.The continuum-extracted H 2 O spectra observed toward the NE and SW-E cores vs. the LTE model spectra derived from the rotation diagrams in Figure12.The models match the observed spectra quite well, taking into account errors in the continuum fits.

Figure 14 .
Figure 14.Visualization of the distribution of bootstrap samples for C 2 H 2 (left) and HCN (right) toward the SF-NE core.The velocity dispersion of individual lines was fixed at 50 km s −1 , but this exact value does not significantly affect the band shape.Dashed lines indicate 68% confidence intervals on the parameters.The column densities and excitation temperatures are reasonably well constrained, but the distributions are somewhat asymmetrical.

Figure 15 .
Figure15.Model spectra for C 2 H 2 (left) and HCN (right), based on the grid bootstrap fits.The red curves represent 50 random samples from the bootstrap; the shaded gray area indicates the spectral region used in the fit.The sampled models are reasonably good fits, but exhibit significant scatter.

Figure 17 .
Figure 17.Schematic illustration of the geometry and excitation around an AGN (top) and in a star-forming knot (bottom).Solid arrows indicate observed radiation; dashed arrows illustrate undetected radiation.For the AGN region, molecular gas in a fiducial shell (orange) surrounding the central AGN + torus (black) is irradiated from only a small solid angle toward the AGN.Therefore, the molecular gas will not be highly excited, despite the high intensity of the continuum emission from the AGN.In the case of a star-forming region, the gas, dust, and starlight are effectively mixed, and therefore, the sketched patch of molecular gas is irradiated isotropically by infrared continuum emission.
rad .We can therefore express the level populations as follows: account both P-and R-branch transitions; the n γ term corrects for stimulated emission.This expression makes use of the fact that the Einstein coefficients A rovib for rovibrational transitions depend only weakly on the rotational quantum number.Since the population of the vibrational levels is purely radiative (see Section 5.1), we can write and since the vibrationally excited levels in Equation (A18) are close enough in energy that the difference can be ignored, we obtain to sufficient accuracy

Table 1
Summary of the Derived Physical Characteristics of All Considered Species and Regions Summary of Derived Physical Characteristics of the Molecular Gas Note.The excitation temperature T ex , total species column density N tot , mean radial velocity V rad , and mean (instrument-corrected) velocity dispersion σ V are listed separately for each temperature component.For C 2 H 2 and HCN, the radial velocity and dispersion are not constrained by the Q-branch fits; for the CO cold components, some lines are unresolved.All values were derived assuming a complete coverage of the background continuum source (see text).
and Table 1).More complex modeling of the CO and H 2 O bands by González-Alfonso et al. (2024) supports this conclusion.