© The Authors 2026

1. Introduction

Compact obscured nuclei (CONs, e.g., Falstad et al. 2021), are bright < 100 pc galaxy nuclei with extremely high column densities (N_H2 ≳ 10²⁵ cm⁻²). They have been mainly identified via the detection of HCN vibrational emission in the millimeter (e.g., Aalto et al. 2015; Falstad et al. 2019), by molecular lines in absorption in the far-infrared (IR, e.g., González-Alfonso et al. 2004, 2015), and through the PAH equivalent width method (García-Bernete et al. 2022a, 2025a). Their power source is, however, hard to establish due to the high nuclear extinction. Even the luminosity arising from these nuclei is difficult to establish from the millimeter continuum and lines observed with high resolution because the trapping of dust-emitted photons increases the dust temperature (T_d) within the nucleus (González-Alfonso & Sakamoto 2019) (hereafter gon19), and the luminosity depends on T_d at the nuclear surface (the “photosphere”), which is better probed with molecular absorption lines. Nevertheless, there is growing evidence that CONs typically harbor an active galactic nucleus (AGN) based on the extremely compact millimeter continuum emission compared to the more extended gas mass distribution, as observed with ALMA (Pereira-Santaella et al. 2021). Still, accurate IR to millimeter continuum models, together with detailed analysis of the hottest dust mid-IR emission from CONs, are scarce.

Here we combine ALMA, James Webb Space Telescope (JWST), and Herschel observations of the eastern (E) nucleus of IRAS 07251−0248 to model its continuum emission from mid-IR to submillimeter (submm) wavelengths. IRAS 07251−0248 is a merging ultra-luminous infrared galaxy (ULIRG) with an IR (6–1500 μm) luminosity of 10^12.45 L_⊙ and nuclear separation of ∼1.8 kpc (Lamperti et al. 2022); its E-nucleus dominates the total luminosity, contributing 91% (Pereira-Santaella et al. 2021). There are no reported observations of HCN vibrational emission from the galaxy, but there are other significant indications that it harbors a CON from the strong absorption in the OH 65 μm doublet (González-Alfonso et al. 2015), the 6.2/3.3 μm PAH equivalent-width ratio (García-Bernete et al. 2025a), and the absorption below the continuum in the blueshifted wing of the CO 2-1 line (Lamperti et al. 2022). IRAS 07251−0248 is very faint in the Chandra soft X-ray band (0.5 − 2 keV) and undetected in the hard bands of Chandra (2 − 7 keV, Iwasawa et al. 2011) and NuSTAR (10 − 24 keV, Ricci et al. 2021). On the other hand, JWST observations have recently shown that the source has a nuclear cosmic-ray dominated chemistry as derived from the high abundances of H₃⁺ (Pereira-Santaella et al. 2024), hydrocarbons (García-Bernete et al. 2025b), and molecular cations (Speranza et al. 2025). We adopt luminosity and angular distances of D_L = 400 Mpc and D_A = 339 Mpc (scale of 1.64 pc/mas), respectively.

2. Observations and results

2.1. ALMA observations of the 667 μm continuum and model

ALMA observations at ν_obs = 413.11 GHz (Band 8, λ_rest = 667 μm) of IRAS 07251−0248 E (program 2023.1.00942.S, PI: M. Pereira-Santaella) were carried out on November 28, 2023 using the 12 m array with high angular resolution (75 × 67 mas² = 123 × 110 pc²) and sensitivity (0.2 mJy/beam). Details of the reduction process are given in Appendix A. The map of the 667 μm continuum (Fig. 1) shows a very bright (peak of ≈25 mJy/beam) and compact emission. To derive its size and flux, we used the method described in Pereira-Santaella et al. (2021) and attempted to fit the observed emission with a point-source, a Gaussian, and a point-source and a Gaussian (pG). The best fit, shown in Fig. 1b, was found for the pG model with parameters listed in Table 1. Appendix A shows that the point-source and the Gaussian fits can be disregarded.

Fig. 1. ALMA 667 μm continuum emission of IRAS 07251−0248 E and model. Left: Observed emission, with colored scale in units of mJy/beam. Middle: Model, composed of an unresolved source (full-width at half maximum, FWHM < 21 mas = 34 pc) and a Gaussian source (FWHM = 89 × 86 mas2 = 146 × 141 pc2). Right: Residuals in μJy/beam.

Due to the high signal-to-noise ratio of the 667 μm emission (> 120), the unresolved source of the pG model was constrained to have an upper limit for the deconvolved size of FWHM ≤ 21 mas = 34 pc, assuming a Gaussian light distribution. Conservatively using the equivalence between a Gaussian source and a disk of uniform brightness (radius ≈ 0.8 × FWHM, Sakamoto et al. 2008), an upper limit for the nuclear source radius is R ≲ 27 pc, yielding a brightness of T_B ≳ 190 K. The 667 μm emission is entirely attributed to dust, with the optical depth given by

$$\begin{array}{c}\hfill {\tau }_{667\mathrm{\mu }\mathrm{m}}=ln{(1-\frac{{10}^{-22}{\lambda }_0^2{D}_L^2{F}_{667\mathrm{\mu }\mathrm{m}}}{2\pi k(1+z)T_d{R}^2})}^{-1}=ln{(1-\frac{7546F_{667\mathrm{\mu }\mathrm{m}}}{T_d{R}^2})}^{-1},\end{array}$$(1)

where k is the Boltzman constant (in cgs units), λ₀ = 667 μm, D_L = 400 Mpc, F_667 μm is in mJy, T_d is the dust temperature in K, and the source radius R is in pc. Using T_d = 500 K, Eq. (1) gives τ_667 μm = 0.5 for R = 27 pc, and has no solution for R ≤ 16 pc. Therefore, the ALMA observations reveal an extremely compact source with a very high column density of hot dust, and will be denoted as the “hot component” (H_C).

For the Gaussian (resolved) continuum source, a size of FWHM = 89 × 86 mas² is found, corresponding to a source (most likely a disk or a torus) seen nearly face-on with a radius of R ≈ 75 pc. Equation (1) gives τ_667 μm ≈ 0.26 for an adopted T_d = 123 K (Section 3.1). This component will be denoted as the “warm component” (W_C).

The spectral energy distribution (SED) of IRAS 07251−0248 E from near-IR to millimeter wavelengths, including JWST (Section 2.2), Spitzer, Herschel/PACS and SPIRE, and ALMA data, is shown in Fig. 2. At 667 μm, we plot separately the flux densities for the H_C (blue) and W_C (light-blue); the black filled circle is the sum of both and lies below the extrapolation of the Herschel/SPIRE data, indicating that a fraction of the total 667 μm emission from the merger is spatially filtered out by ALMA. The continuum measurements at 1.2 mm are included (Table 3 in Pereira-Santaella et al. 2021), where the black open circle indicates the unresolved emission (R < 82 pc, which is large enough to include the emission from the H_C and W_C), and the cross indicates the resolved emission (∼500 × 300 pc²).

Fig. 2. Spectral energy distribution (SED) of IRAS 07251−0248 E from λrest = 2.6 to 1200 μm and model. In black, the full JWST NIRSpec & MIRI/MRS spectrum is shown, together with Spitzer, Herschel/PACS (from 60 to 160 μm, both photometric and spectroscopic) and SPIRE (from 230 to 470 μm), and ALMA (667 and 1200 μm) data. The colored circles at 667 μm indicate the contributions by the unresolved (blue) and resolved (light-blue) components inferred from the model in Fig. 1. At 1.2 mm, the open black circle and the cross indicate the flux densities of the unresolved (< 82 pc) and resolved (540 × 330 pc2) components, respectively (Pereira-Santaella et al. 2021). The model includes the HC in blue (dashed: unattenuated emission; solid: attenuated), the WC in light-blue (extincted), a cold component (CC) in orange, the stellar component in gray, and the “hidden” component (whose attenuated emission is added to that of the WC, see Section 3.2) in dotted-dashed green. Red is total. The insert panels show details of the model results. Upper left: The Tdust profile of the HC. Upper right: The optical depths τ of the HC and the WC as function of wavelength. Lower middle: The Tdust profile of the HC as function of the optical depth from surface at 7 and 14 μm.

2.2. JWST observations of the mid-IR continuum and molecular bands

The JWST observations of IRAS 07251−0248 E, recently presented by García-Bernete et al. (2025b) and Speranza et al. (2025), were carried out as part of the JWST GO Cycle 2 Large Program #ID:3368 (P.I. L. Armus and A. Evans) including NIRSpec and MIRI/MRS integral field spectroscopy of the nucleus. Details of the reduction process are given in Appendix B.

The JWST spectrum (Fig. 2) is characterized by deep absorption of both solid- and gas-phase bands. Besides the silicate 9.7 μm and 18 μm features, there are deep absorption features due to H₂O ice at 3 (stretching), 6 (bending), and 13 μm (libration), and due to hydrogenated amorphous carbon (a-C:H) grains at 6.85 and 7.25 μm. Superposed on this absorbed continuum, strong gas-phase absorption bands due to CO, CO₂, H₂O, HCN, and C₂H₂ are detected, as well as additional weaker bands mostly associated with hydrocarbons (García-Bernete et al. 2025b) and molecular cations (Speranza et al. 2025). In this absorption-dominated mid-IR spectrum, no high ionization lines that would directly probe an AGN are detected (Spoon et al. 2022). Here we focus on the H₂O ν₂ = 1 − 0 band at ∼5.2 − 7.0 μm displayed in Fig. 3, although the HCN, C₂H₂, CH₄, and CO₂ bands are also modeled in Appendix F. The H₂O ν₂ = 1 − 0 band has previously been reported in extragalactic sources by García-Bernete et al. (2024b, 2025b), González-Alfonso et al. (2024), Buiten et al. (2024, 2025). We used the H₂O spectroscopic parameters tabulated in the HITRAN2020 database (Gordon et al. 2022).

Fig. 3. H2O ν2 = 1 − 0 band in IRAS 07251−0248 observed with JWST MIRI/MRS. Panels a and c show the observed spectrum, with the orange line indicating the adopted continuum level and the light-blue curve showing the contribution from the warm component. Panels b and d show the continuum normalized spectrum, with the blue line showing the model result for the hot component and the dashed line indicating its covering factor (fhot). We label the main transitions that potentially contribute to the H2O band as shown in González-Alfonso et al. (2024), but some very highly-excited lines are not detected in IRAS 07251−0248. Note that all H2O absorption features are blueshifted relative to the labels.

The H₂O band in IRAS 07251−0248 E shows some similarities with, but also important differences from, that observed towards the lower luminosity LIRG VV 114 E SW-s2 (González-Alfonso et al. 2024). In both sources, the bands are in pure absorption (in contrast with that observed in the LIRG II Zw96-D1, where the P branch is in emission; García-Bernete et al. 2024b). However, the continuum-normalized absorption due to low-excitation H₂O lines in IRAS 07251−0248 E is stronger than in VV 114 (see, e.g., the blended feature at ≈6.05 μm absorbing ≈70% of the continuum in IRAS 07251−0248 but only ≈35% in VV 114). The opposite happens with the very high-excitation H₂O lines, which are strong in VV 114 but most of them (above ∼1500 K) undetected in IRAS 07251−0248 (see also Fig. D.2), indicating significantly lower excitation in this ULIRG (∼200 K, García-Bernete et al. 2025b). Interestingly, both sources show the H₂O lines blueshifted by ∼160 km s⁻¹ relative to the emission lines (Fig. D.1), but the blueshift in VV 114 increases with higher line excitation while the opposite happens in IRAS 07251−0248 E. Close inspection of the H₂O line profiles in Figs. D.1 and D.2 indicates the presence of two kinematic components, and the most blueshifted one (with absorption at velocities up to −500 km s⁻¹) is only seen in the low excitation lines.

3. A model for the mid-IR to millimeter continuum

3.1. The continuum model

Since ALMA observations at 667 μm reveal an H_C with an extreme column density of hot dust that is expected to emit in the mid-IR, and JWST indeed indicates the presence of a mid-IR component with a column density high enough to generate strong absorption in the molecular bands, our hypothesis is that both observations trace the same continuum nuclear component (with R_h ≲ 27 pc as derived from ALMA). Additional support for this assumption is found in the CO 2-1 230 GHz rotational line profile displayed in Fig. D.1 (Lamperti et al. 2022), which shows an asymmetric profile with apparent blueshifted absorption between −500 and −350 km s⁻¹ falling below the continuum. Assuming that this absorption is produced by the same gas responsible for the low-excitation H₂O absorption found at similar velocities, the continuum behind both CO and H₂O line absorption, at 1.3 mm and ∼6 μm, should arise from the same physical source. In the following, we explore and test this scenario further.

The high optical depth of the H_C in the submm indicates that trapping of continuum photons will raise T_d within the hot nucleus, and we have thus generated a grid of “greenhouse” models using the method described in gon19. AGN models are used here, together with an absorption coefficient of dust at < 10 μm similar to that obtained by Fritz et al. (2011) towards the galactic center, and an absolute calibration in terms of total hydrogen column density (N_H, see Appendix C) consistent with the results by Planck Collaboration XXV (2011). In our continuum models, we include the absorption due to solid-phase silicates, but not due to ices or a:C-H, so that we aim to fit the envelope of the 2.6 − 8 μm continuum (see Fig. 2). For simplicity, flat density profiles are used.

The grid of greenhouse models for the H_C covers log N_H(cm⁻²) = 24.0 − 25.6 and unattenuated luminosities log L_unatt(L_⊙) = 11.0 − 12.6, both with logarithmic steps of 0.2¹. A set of mid-IR 2.6 − 11 μm continuum points unaffected by ice, a:C-H, or gas-phase bands, together with the flux density at 667 μm of the H_C, is used to find the best-fit model. The unattenuated nuclear emission is extinguished by a layer of foreground cold dust, which is required because the predicted silicate absorption at 9.7 μm generated by any relevant greenhouse model is found to be much shallower than the observed feature. The fitted parameters are, in addition to the properties of the greenhouse model (N_H and L_unatt), its radius R_h and the optical depth of the foreground layer evaluated at 6 μm, ${\tau }_{6\mathrm{\mu }\mathrm{m}}^{\mathrm{fg}}$. Our best-fit model for the H_C, with parameters listed in Table 2, is shown with a solid-blue line in Fig. 2; the dashed-blue line indicates the unattenuated emission (i.e., the emission that would be observed without the foreground layer).

Once the best fit model for the H_C is obtained, the W_C is fitted by using JWST continuum points at 11.5 − 26 μm, together with the flux density at 667 μm of the resolved component. A graybody (single-T_d) extincted by the same foreground layer as for the H_C (i.e., with the same ${\tau }_{6\mathrm{\mu }\mathrm{m}}^{\mathrm{fg}}$) is adopted. We add to this graybody an extincted blackbody at 230 K, primarily contributing at < 10 μm, which has the effect of producing a better match to the H₂O band (see Section 3.2). The fitted parameters for the W_C are T_d, the intrinsic optical depth evaluated at 6 μm (${\tau }_{6\mathrm{\mu }\mathrm{m}}^{\mathrm{int}}$), and the source radius R_w. The resulting “warm” continuum, with T_d = 123 K, is complemented by a “hidden” component (discussed in Section 3.2) and shown in light-blue in Fig. 2. The fitted parameters are listed in Table 2. Finally, a graybody “cold” (T_d = 45 K) component (C_C) is added to the model to account for the remaining far-IR emission that cannot be explained with the H_C and W_C (orange curve in Fig. 2). The C_C also includes the fraction of the 667 μm emission that is filtered out by ALMA.

The fit to the SED of IRAS 07251−0248 E indicates that the mid-IR (< 10 μm) continuum envelope and the 667 μm unresolved emission can indeed be matched with a single foreground-extincted nuclear component with unattenuated luminosity L_unatt ≈ 10¹² L_⊙. The sizes obtained for this H_C (R_h = 13.2 pc) and for the W_C (R_w = 71 pc) are consistent with the ALMA 667 μm measurements for the unresolved and resolved components. Our dust-to-gas mass ratio calibration (Appendix C) yields a gas mass (M_H2 in Table 2) dominated by the C_C. The total M_H2 = 5.3 × 10⁹ M_⊙ is comparable with the value obtained from CO 1 − 0 (4.3 × 10⁹ M_⊙, González-Alfonso et al. 2015) and CO 2 − 1 (5.0 × 10⁹ M_⊙, Lamperti et al. 2022), while the CO 2 − 1 value for the central (r < 250 pc) region is consistently lower (2.5 × 10⁹ M_⊙, Pereira-Santaella et al. 2021). The inserts in Fig. 2 indicate extreme mid-IR optical depths for the H_C and also for the W_C, as well as hot dust (> 500 K) within the bulk of the H_C. Nevertheless, the surface temperature of the H_C is much more moderate (≈200 − 250 K), owing to the nearly complete conversion of the luminosity into blackbody emission from the nuclear photosphere (T_d ≈ (L_unatt/(4πR_h²σ_SB))^1/4 = 240 K).

3.2. The 6 μm photosphere traced by the H₂O ν₂ = 1 − 0 band

The nuclear H_C, with extreme column densities and accounting for the bulk of the mid-IR < 8 μm continuum, naturally generates strong mid-IR absorption in the gas-phase molecular bands. If the H₂O excitation is due to radiative pumping, the H₂O ν₂ = 1 − 0 band is key to check whether the inferred T_d profile can reproduce the observed line ratios and also provides valuable information on the gas-phase H₂O abundance and the kinematics of the nuclear region.

Using the T_d profile obtained from our best-fit continuum model for the H_C (Fig. 2), we have generated a grid of non-LTE, non-local radiative transfer models for H₂O (e.g., González-Alfonso et al. 2024) by varying only the H₂O abundance relative to H₂ (X_H2O, assumed uniform) and the gas velocity field. Dust and gas are uniformly mixed within the nucleus, and due to the extreme nuclear mid-IR extinction only the 6 μm photosphere can be traced with the H₂O ν₂ = 1 − 0 band. With the redshift obtained from emission lines (z = 0.08778, Speranza et al. 2025, see also Fig. D.1), we interpret the blueshift of the H₂O lines as an expansion of the nuclear surface and apply across the nuclear photosphere an outflow velocity increasing linearly from 75 to 170 km s⁻¹. Absorption at higher velocities (up to ∼ − 500 km s⁻¹) is also observed in low-excitation lines (Section 2.2), and thus an additional H₂O shell flowing with velocities of 150 − 400 km s⁻¹ was placed in front of the nucleus.

Due to the strong absorption associated with H₂O ice around 6 μm and a:C-H grains at 6.85 μm, the model predictions are better compared with the continuum-normalized spectrum in Figs. 3b,d. However, the models that give a good match to most lines systematically overestimate the absorption lines located close to the absorption troughs, particularly at 6.85 μm. This specific model departure indicates the presence of an additional mid-IR component, unaffected by ice-mantle and H₂O absorption, that dilutes the H₂O line absorption at wavelengths where its strength is comparable to the observed continuum. Specifically, the flux density at frequency ν can be written as

$$\begin{array}{c}\hfill F_{\nu }=F_{\nu }^{H_C}{f}_{\nu ,\mathrm{norm}}^{{\mathrm{H}}_2\mathrm{O}}+{\displaystyle \sum _{i\ne H_C}}{F}_{\nu }^i,\end{array}$$(2)

where F_ν^H_C and $F_{\nu }^i$ are the (attenuated) flux densities from the H_C and all other components i (including the additional mid-IR component), and f_ν, norm^H₂O is the continuum-normalized H₂O absorption predicted by our model for the H_C. Since the observed continuum $F_{\nu }^{\mathrm{cont}}$ has the same expression as Eq. (2) but excluding the f_ν, norm^H₂O factor, the continuum-normalized spectrum is given by

$$\begin{array}{c}\hfill \frac{F_{\nu }}{F_{\nu }^{\mathrm{cont}}}=1+f_{\mathrm{hot}}(f_{\nu ,\mathrm{norm}}^{{\mathrm{H}}_2\mathrm{O}}-1),\end{array}$$(3)

where we have defined a λ−dependent covering factor for the H_C, which is the fractional contribution by the H_C to the observed continuum:

$$\begin{array}{c}\hfill f_{\mathrm{hot}}=1-\frac{{\sum }_{i\ne H_C}{F}_{\nu }^i}{F_{\nu }^{\mathrm{cont}}}·\end{array}$$(4)

We use Eqs. (3) and (4) to compare the H₂O models with the observed continuum-normalized spectrum in Figs. 3b,d. At wavelengths where f_hot becomes significantly lower than unity, because the H_C continuum is severely attenuated by ice or a:C-H and thus the expected $F_{\nu }^i$ from any other component i becomes comparable to $F_{\nu }^{\mathrm{cont}}$, the continuum normalized absorption due to any H₂O line in front of the H_C will be reduced according to Eq. (3). The 123 K component for the W_C in Table 2 yields negligible emission below 10 μm, and we find that a blackbody of 230 K, attenuated in the same way as the H_C and W_C (green dotted-dashed curve in Fig. 2), is required to successfully correct the mentioned discrepancies. This 230 K component does not dominate the continuum emission at any wavelength, so we denote it as the “hidden” component. Although its origin is uncertain, we tentatively associate it with the W_C (and add it to the W_C; light-blue curve in Fig. 2). In Figs. 3b,d, f_hot is shown as black dashed lines.

Our best-fit model for the H₂O band (Figs. 3b,d) is obtained for X_H2O = 8 × 10⁻⁵ within the nucleus and N_H2O ≈ 1.6 × 10¹⁸ cm⁻² for the surrounding H₂O shell. A more detailed comparison between the observed line profiles and model predictions is shown in Fig. D.2 and discussed in Appendix D, but here we emphasize the following points: first, that ground-vibrational levels of H₂O are excited through the pumping of the ν₂ = 1 state by the strong mid-IR radiation that bathes the nuclear region and its surroundings. Second, the high X_H2O ∼ 10⁻⁴ is in line with chemical models of hot regions where H₂O ice has been vaporized and nearly all gas-phase oxygen not locked into CO is converted into H₂O (e.g., Nomura & Millar 2004), as well as with models of C-type shocks (Kaufman & Neufeld 1996). Third, despite very high T_d within the bulk of the H_C, the model captures the observed moderate H₂O excitation due to the moderate T_d values of the nuclear photosphere (lower insert in Fig. 2), which is the only region the absorption band can probe due to extinction. We show in Appendix F that the HCN, C₂H₂, CH₄, and CO₂ bands are also reproduced with our model.

4. Discussion and conclusions

The schematic illustration of the model source shown in Fig. 4 reflects the structure inferred from the observations of IRAS 07251−0248 E and the model: the expanding nuclear region (H_C), is surrounded by outflowing gas (the high velocity foreground shell seen in H₂O), and a nearly face-on torus (W_C). The H_C and the W_C are both embedded in an envelope representing the C_C, which dominates the far-IR emission between 50 and 300 μm. In our illustration, the C_C lies in front of the H_C and W_C, and thus the C_C also represents the foreground layer that extinguishes the continuum emission from the H_C and the W_C and is responsible for the 9.7 and 18 μm silicate absorption. The alternative, that the W_C is a spherical shell surrounding the H_C, as Fig. 1 might suggest, and is thus responsible for the foreground extinction of the H_C, is slightly disfavored because the high optical depth of the W_C (right insert in Fig. 2) would totally extinguish the mid-IR continuum emission from the H_C. Nevertheless, it is still possible that the high-velocity outflow has cleared the path to the H_C along the line of sight and the extinction is caused by the residual material of the W_C or by the outflowing material. The C_C should include the far-IR emission from the western galaxy of the merger, with 9% contribution to the overall luminosity (∼2.5 × 10¹¹ L_⊙, Pereira-Santaella et al. 2021). Regardless of the origin of the foreground extinction layer, it will probably have spatially varying chemical properties, since it is fed by outflowing material from the direction of the nucleus, generating the observed H₂O ice and a:C-H absorption as the gas-phase H₂O freezes out and the hydrocarbons attach to carbonaceous grains (García-Bernete et al. 2025b).

Fig. 4. Schematic representation of the model source. The extremely compact and expanding core (HC) is shown in blue (the 6 μm photosphere is its external surface), with gas flowing outward (the surrounding shell outflowing at ∼ − 400 km s−1) in green. The starburst torus shown in magenta represents the WC. The torus also appears to be expanding, because the far-IR lines are blueshifted (Appendix E, Fig. E.1). Both the HC and WC are embedded in an extended cold component (CC), in orange. The model is axisymmetric about the line of sight.

The W_C is a very optically thick structure (right insert in Fig. 2) surrounding the compact nucleus. Its far-IR emission is much stronger than the H_C due to its larger size, which is similar to the size of the warm and optically thick regions that in other ULIRGs generate the far-IR absorption in rotational lines of OH, H₂O, and other molecular species (e.g., Arp 220 and Mrk 231; González-Alfonso et al. 2012, 2014). The far-IR continuum emission from these regions is usually diluted within the emission from colder regions that also contribute to the far-IR, which is the reason why (together with re-emission in the line) the far-IR optically thick lines absorb only a fraction (≲20%) of the total continuum (González-Alfonso et al. 2015). A similar situation applies to IRAS 07251−0248, as the C_C dominates the far-IR and the continuum from the W_C should be diluted (Fig 2). The Herschel/PACS spectra around the wavelengths of the OH doublets at 119, 79, and 65 μm, shown in Fig. E.1, exhibit optically thick absorption in the OH doublets and also in lines of H₂O, CH⁺, and CH, so that their absorption strengths are sensitive to the fractional contribution of the W_C to the far-IR. We show in Appendix E a model for these lines that reasonably fits the absorption troughs, indicating that the modeled 65 − 120 μm continuum level generated by the W_C is consistent with the available far-IR spectroscopic data. All far-IR lines are blueshifted by ∼150 km s⁻¹ (Fig. E.1), indicating that the W_C is also expanding.

4.1. Uncertainties in the inferred parameters

The extreme N_H2 and M_H2 inferred in the H_C (Table 2) could still be lower limits, because the adopted inner radius of the source is 0.15 × R_out, enabling high T_d > 1000 K for r < 8 pc (left insert in Fig. 2). However, if the dust is concentrated in a thinner shell close to the photosphere, as suggested by the outflow observed in the mid-IR molecular bands, or if the greenhouse effect in our models overestimates the actual T_d values (see discussion in gon19), both N_H2 and M_H2 would increase because of the overall lower T_d.

As shown in Table 2, the apparent luminosities (L_att) of the three components yield a total output power of 2.3 × 10¹² L_⊙, but the combined unattenuated luminosities (L_unatt) obviously give a higher value. This indicates that the values of L_unatt are, for both the W_C and the C_C, in part re-emission from more inner components, i.e., they are not only due to the power sources located within their physical regions (e.g., González-Alfonso et al. 2004). The W_C surrounds the H_C and part of its emission will be an effect of the H_C illumination. Similarly, if the C_C surrounds both the W_C and the H_C, a fraction of its emission will be an effect of heating by the innermost components. Defining the “intrinsic luminosity” L_int of a component as the unattenuated luminosity generated by only the power sources within the physical region, and under the assumption that isotropic emission is applicable to the optically thick H_C and W_C components, the values of L_int would be 10¹², ∼10¹², and ∼3 × 10¹¹ L_⊙ for the H_C, W_C, and C_C, respectively. These estimates also assume that nearly the entire luminosity from the H_C is absorbed and re-emitted by the W_C, which would then intrinsically generate only half of its L_unatt value; the C_C would then produce the remaining emission that accounts for the total output power, including the emission from the western galaxy.

Nevertheless, isotropic emission must not necessarily apply to the optically thick H_C and/or W_C components, in which case their L_int will depart from the isotropic values quoted above (see Efstathiou et al. 2014, 2022). If the morphology of the W_C was a flat disk seen nearly face-on, L_unatt(W_C) = 2πR_W²σ_SBT_d⁴ ∼ 10¹² L_⊙ for R_W = 71 pc and T_d = 123 K (Table 2). The actual scenario is most likely intermediate between the spherical and flat morphologies (Fig. 4). Similarly, if the H_C has a structure elongated along the line of sight, its L_int would exceed 10¹² L_⊙, and the W_C emission would be predominantly attributable to external heating. This possibility is supported by the fact that, rather than representing a separate physical structure, the W_C appears to be connected to the H_C, since both exhibit high column densities and similar kinematics (Appendix E) although quite different brightnesses. In Table 2 we provide plausible ranges for L_int, taking these uncertainties into consideration.

4.2. A CON hiding a quasar, surrounded by a starburst

From a comparative perspective, the high brightness of the H_C favors a dominant AGN contribution to the luminosity. The luminosity surface density, Σ_bol ≈ 4.6 × 10⁸ L_⊙ pc⁻², is ∼1 dex higher than the values found in other CONs with comparable sizes and gas mass surface densities (Σ_H2 ≈ 5 × 10⁵ M_⊙ pc⁻²), such as NGC 4418 ((3 − 7)×10⁷ L_⊙ pc⁻², Sakamoto et al. 2013; González-Alfonso et al. 2012), IC 860 ((1 − 8)×10⁷ L_⊙ pc⁻², Aalto et al. 2019; gon19) or Arp 220 W (∼3 × 10⁷ L_⊙ pc⁻², gon19)². Our model for the H_C yields an unattenuated SED that peaks in the mid-IR (Fig. 2), while the continuum models in gon19 consider lower values of Σ_bol and the corresponding SEDs peak in the far-IR.

Nevertheless, the enormous accumulation of gas in the H_C requires consideration of the role of a potential transitory burst of star formation. We first note that the H₂O absorption line shapes do not show any distinct feature at central velocities (Figs. D.1 and D.2), meaning that the whole mid-IR photosphere of the H_C is outflowing and suggesting the entire expulsion of gas from the nuclear region. If responsible for the outflowing gas, a nuclear starburst could not have a duration much longer than the outflow dynamical timescale of t_dyn ∼ R_h/v = 8 × 10⁴ yr, and, indeed, the free-fall time, t_ff = (3π/32Gρ)^1/2 = 4.1 × 10⁴ yr (calculated from the averaged n_H2 = 5.8 × 10⁵ cm⁻³), accomplishes this constraint. However, with these short timescales a pure starburst scenario would imply that most stars are actually in the protostar phase with the luminosity arising from gas accretion, while expelling the inflowing material that feeds them.

Thompson et al. (2005) argued that for an optically thick starburst disk radiating at its Eddington limit for dust, a characteristic maximum Σ_bol ∼ 10⁷ L_⊙ pc⁻² is attained that is independent of Σ_gas (≈1.36 × Σ_H2). This limit for Σ_bol falls too short to explain the inferred value in the H_C, and an AGN would then be responsible for the bulk of the luminosity. Andrews & Thompson (2011) postulated the possibility of hot starbursts with extreme values of Σ_bol ∼ 5 × 10⁸ L_⊙ pc⁻² for Σ_gas ∼ 5 × 10⁵ M_⊙ pc⁻², but attributed them to the parsec-scale region fueling a bright AGN.

More recently, Grudić et al. (2018) proposed that for sufficiently high Σ_H2 > Σ_crit ∼ 3 × 10³ M_⊙ pc⁻², feedback fails due to the strong gravitational force and an extraordinary burst of star formation converts most of the gas into stars within a few t_ff. The key feature of these models is that the momentum deposition per stellar mass is assumed to be independent of environment; that is, the “τ_IR boost” characteristic of radiation pressure in optically thick components is neglected owing to presumable leakage of photons along optically thin paths. We note that the maximum Σ_gas ∼ 1.3 × 10⁴ M_⊙ pc⁻² considered in their simulations is ∼40× lower than the value of the H_C, but here we assume that their results can be extrapolated to our more extreme conditions. Adopting a zero-age main sequence (ZAMS) population with 3 × 10³ erg s⁻¹ g⁻¹, a stellar mass of M_* = 6.3 × 10⁸ M_⊙ is required to account for 10¹² L_⊙. The current star-formation efficiency would then be SFE = M_*/(M_* + M_gas)∼0.63 and, with a formation timescale of ∼2 t_ff, the star-formation rate would be SFR ∼ 7.7 × 10³ M_⊙ yr⁻¹ with a per-freefall SFE of ϵ_ff ∼ 0.6. Within the uncertainties, the values of SFE and ϵ_ff agree with the results of Grudić et al. (2018) simulations.

Nevertheless, the above pure-starburst scenario faces serious drawbacks in the case of the H_C in IRAS 07251−0248. First, the implied values of Σ_* ∼ 10⁶ M_⊙ pc⁻² are found in a few very dense nuclear star clusters (NSCs), but on smaller spatial scales (≲5 pc) and with lower M_* (Grudić et al. 2019; Neumayer et al. 2020). Second, the concentration of the observed gas mass in the H_C would be hard to explain with such a burst of star formation, because accretion of ∼10⁹ M_⊙ onto a volume with radius R_h ∼ 13 pc would also have to be achieved within ∼1 t_ff (see also Grudić et al. 2019), requiring unrealistic accretion rates. Third, protostars do not efficiently generate the cosmic rays (CRs) required to explain the observed H₃⁺, CR-dominated chemistry, and carbonaceous grain erosion in the H_C (Pereira-Santaella et al. 2024; Speranza et al. 2025; García-Bernete et al. 2025b). Fourth, the H₂O line profiles strongly suggest a collective motion of the whole nuclear gas that is better explained by a single feedback event rather than by the effect of many individual protostars. Finally, the “greenhouse effect” required to explain the high submm brightness of the H_C translates into the “τ_IR boost”, which cannot be neglected for the extreme Σ_gas in the H_C, and the observed feedback would be hard to explain without this effect.

Conversely, an AGN burst responsible for the observed luminosity and feedback would have enabled smoother gas accretion onto the nuclear region, owing to the different spatial scales and solid angles of AGN feedback and nuclear gas accretion, and the stochastic nature of the former. While the previous formation of an NSC in IRAS 07251−0248 E is a natural outcome of the highly concentrated ISM, the extreme Σ_bol strongly suggests a limited star-formation efficiency within a compact obscured nucleus (CON) that is powered by an AGN currently emitting at quasar level, in agreement with Pereira-Santaella et al. (2021).

With uncertainties related to geometry in optically thick sources, it appears that the W_C is forming stars at high rates equivalent to at least ∼5 × 10¹¹ L_⊙, supporting the view of the composite nature of ULIRGs (e.g., Veilleux et al. 2009) with the specificity that the W_C is a nuclear starburst characterized by high column densities. The high gas mass accumulation at the center of the galaxy both drives a starburst in situ, and thus the formation of an NSC, and efficiently feeds the central supermassive black hole (Neumayer et al. 2020), exemplifying the black hole and star formation co-evolution in its most extreme form. The transition from star formation to black hole accretion in the H_C may be smooth under these conditions, and tidal disruption and capture of (proto)stars close to the black hole may be an efficient black hole growth channel (Strubbe & Quataert 2009).

4.3. Outflow energetics: a hot bubble driving the H_C expansion?

With the values of M_H2 and R from Table 2 and a velocity of 160 km s⁻¹, the mass outflow rate associated with the H_C is estimated as Ṁ = μ M_H2 v/R = 2.9 × 10³ M_⊙ yr⁻¹, where μ = 1.36 corrects for species other than hydrogen. Here we have assumed that the entire nuclear source is expanding, although we can only probe the kinematics of its mid-IR photosphere. The momentum rate is then Ṗ = Ṁv = 2.7 × 10³⁶ dyn = 21.5L_AGN/c and the energy flux is Ė = 0.5 Ṁv² = 2.05 × 10⁴³ erg s⁻¹ = 5.4 × 10⁻³ L_AGN, where we have adopted L_AGN = 10¹² L_⊙ (Table 2).

Radiation pressure on dust grains could, in principle, drive the observed outflow if it can overcome the opposite effect of (self-)gravitation on the gas (e.g., Ishibashi & Fabian 2015). Taking into account the backpressure effects in optically thick environments gon19 and with a gas mass of μ M_H2 = 3 × 10⁸ M_⊙ within the H_C (Table 2), the gravitational force would, however, balance the force due to radiation pressure at the H_C surface as long as the combined black hole and stellar mass attains ∼2.5 × 10⁸ M_⊙. This is plausible, and thus the effect of radiation pressure in driving the outflow is uncertain.

Another possibility consists of an inner hot bubble inflating the nucleus such that molecular outflow is generated as the bubble sweeps up the concentrated nuclear gas. We have explored this point using the analytic self-similar solutions developed by Faucher-Giguère & Quataert (2012) (their Appendix A), assuming that the hot shocked wind conserves energy. We first constrain the swept-out gas mass to μ M_H2 = 3 × 10⁸ M_⊙ at 13.2 pc with a flat density profile, and assumed that 1/2 of the mechanical energy injected by the AGN (v_wind L_AGN/2c, where v_wind is the inner wind velocity and L_AGN = 10¹² L_⊙) goes into kinetic motion of the swept-up gas. Then, for v_wind = 2 × 10⁴ km s⁻¹, a shock velocity driven into the ambient gas of 170 km s⁻¹ at 13.2 pc is obtained, in agreement with the blueshift of the peak absorption observed in the molecular lines. However, the flow time of ∼5 × 10⁴ yr is about half the proton cooling time calculated using the two-temperature effects derived from thermal uncoupling of e⁻ and protons (Faucher-Giguère & Quataert 2012), and thus a partially energy-conserving phase is not ruled out.

This hot bubble scenario is an ultracompact version of the (partially) energy-conserving phase of outflows driven by AGN (Tombesi et al. 2015; Feruglio et al. 2015), with much lower shock velocities that are a consequence of the enormous gas masses that are being swept-up at small radii. The bulk of the gas will not escape from the potential well of the galaxy, and thus a series of CON phases are expected to occur over longer timescales, while holes in the structure will be increasingly opened, generating in later times more extended and faster OH and CO outflows. The interesting feature of the proposed hot bubble scenario is that the inner wind shock presumably accelerates cosmic-ray particles (e.g., Zubovas & King 2012), a key ingredient of the strong H₃⁺ absorption (Pereira-Santaella et al. 2024) and unique carbon-rich (García-Bernete et al. 2025b) and molecular cation (Speranza et al. 2025) chemistry observed in IRAS 07251−0248 E. The strength of linking kinematics and chemistry may favor this scenario over radiation pressure as the origin of the observed outflow.

Acknowledgments

We thank the referee for the constructive comments that helped improve the clarity of the manuscript, and acknowledge the DD-ERS teams for developing their observing program with a zero–exclusive–access period. EG-A thanks the Spanish MICINN for support under projects PID2022-137779OB-C41 and PID2023-146667NB-I00. MPS acknowledges support under grants RYC2021-033094-I, CNS2023-145506 and PID2023-146667NB-I00 funded by MCIN/AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR. IGB is supported by the Programa Atracción de Talento Investigador “César Nombela” via grant 2023-T1/TEC-29030 funded by the Community of Madrid. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST; and from the European JWST archive (eJWST) operated by the ESAC Science Data Centre (ESDC) of the European Space Agency. These observations are associated with program #3368. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2023.1.00942.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

References

Appendix A: ALMA observations and models for the 667 μm continuum

Fig. A.1. ALMA 667 μm continuum emission of IRAS 07251−0248 E and alternative models: a point source (upper) and a Gaussian (lower). Left: Observed emission, with colored scale in units of mJy/beam. Middle: Model. Right: Residuals in μJy/beam. These models give residuals much more prominent than our best point-source + Gaussian model fit in Fig. 1.

We obtained high angular resolution ALMA observations of IRAS 07251−0248 E in Band 8 at ∼410 GHz (rest-frame ∼667μm; program 2023.1.00942.S, PI: M. Pereira-Santaella). Four 1.875 GHz bandwidth spectral windows with 10 km s⁻¹ channels were defined for these observations. In this paper, we focus on the continuum detected after combining all the line-free channels in the spectral windows. We used the ALMA reduction software CASA (v6.5.4; McMullin et al. 2007) to calibrate and clean the data using the standard pipeline. For the cleaning, we applied the Briggs weighting with a robustness parameter of 1.0. To improve the quality of the continuum image, we applied a round of phase self-calibration using the emission from the bright H₂O 448 GHz line (Pereira-Santaella et al. 2017). The final continuum image was corrected for the primary beam. The synthesized beam full-width half-maximum (FWHM) for the combined continuum image is 75×67 mas² and the sensitivity is 0.20 mJy beam⁻¹. For Band 8, the absolute flux accuracy is ∼20% (ALMA Technical Handbook).

Similarly to Pereira-Santaella et al. (2021), we modeled the continuum of IRAS 07251−0248 E to determine the flux, size, and position of the detected emission. Briefly, we used three simple models consisting of a point-source, a Gaussian, and a point-source + Gaussian (pG model). These models were convolved with the beam and their parameters varied to minimize the χ² obtained from the comparison with the observed image. Figure A.1 shows the results of the fit for the point-source and Gaussian model, which both show significant residual structures indicating that they fail to reproduce the nuclear emission of this object. Instead, the fit is much better for the pG model, which shows no evident residuals (Fig. 1).

Appendix B: JWST observations and data reduction

The JWST observations of IRAS 07251−0248 E have been described in detail by García-Bernete et al. (2025b) and Speranza et al. (2025). Data were taken with the integral-field spectrograph NIRSpec, using the grating filter pair G395H (2.9 − 5.3 μm with resolution R ∼ 2700; Jakobsen et al. 2022; Böker et al. 2022), and MIRI MRS (4.9 − 28.1 μm, R ∼ 3700 − 1300; Rieke et al. 2015; Wells et al. 2015; Wright et al. 2015). The data were reduced following the standard MRS pipeline procedure (Labiano et al. 2016), additionally using extra steps to identify and correct hot and cold pixels (Pereira-Santaella et al. 2024; García-Bernete et al. 2024a). The spectrum was extracted according to the method described in Pereira-Santaella et al. (2022) and García-Bernete et al. (2022b, 2024a,b).

Appendix C: The absorption coefficient of dust

Our adopted 6 μm-normalized absorption coefficient of dust from mid-IR to submm wavelengths is shown in Fig. C.1 and compared with the mid-IR values inferred towards the Galactic Center (GC) by Fritz et al. (2011). The latter includes mid-IR features due to ices and aliphatic hydrocarbons that we have ignored, but otherwise our curve fits the overall GC curve for λ < 10 μm. However, we have modified the GC curve for λ = 10 − 20 μm according to Rosenthal et al. (2000), which is in turn based on the work of Draine & Lee (1984). This modification gives a better fit to the IRAS 07251−0248 spectrum at the quoted wavelengths. At λ > 20 μm, we use an emissivity index similar to that of the GC curve (β = 1.85) all the way to millimeter wavelengths.

Planck Collaboration XXV (2011) have given an absolute calibration of τ_250μm/N_H = 2.32 × 10⁻²⁵ cm², which for β = 1.85 yields μ × DGR × k_670μm = 2.2 × 10⁻² cm² g⁻¹, where DGR is the dust-to-gas ratio by mass and μ = 1.36 corrects for all species other than hydrogen. We have adopted this calibration to calculate the H₂ masses in Table 2, which roughly agree with independent methods based on CO emission (Section 3.1). This indicates a suitable calibration for at least the C_C component that dominates the gas mass budget, but the results are more uncertain for the W_C and mostly for the H_C. At 6 μm we obtain τ_6μm/N_H = 1.3 × 10⁻²³ cm², which is used to calculate the molecular abundances in our models where dust and molecules are evenly mixed. The hydrogen column density in terms of the optical depth at 667 μm is N_H = 2.6 × 10²⁵ τ_667μm cm⁻².

Appendix D: Details of the H₂O band

Figure D.1 clearly illustrates the blueshift of the H₂O ν₂ absorption lines relative to several emission lines detected with MIRI/MRS. (Note that the absorption features close to some emission lines are due to R branch lines of the also blueshifted HCN 2ν₂ band.) The CO 2-1 rotational line (see also Fig. C.1 in Lamperti et al. 2022) shows an asymmetric shape with the blueshifted part of the profile weaker than the redshifted part. In view of the absorption of the continuum detected in the velocity range [−570,−350] km s⁻¹, the asymmetry of the line core is also attributable to the absorption of the continuum by the outflowing approaching CO gas, most likely on spatial scales of the W_C. In the lower two panels of Fig. D.1, the profiles correspond to low-excitation H₂O lines that also display absorption up to ∼ − 600 km s⁻¹; other more excited lines do not show this high-velocity absorption (Fig. D.2).

A detailed comparison between a number of observed H₂O line shapes and the best-fit model predictions is presented in Fig. D.2. As shown in the inserted energy level diagram, the displayed lines (mostly within the R branch) span a broad range of level excitation in the ground vibrational state, from the ground ortho level (1_0, 1) to the 10_1, 10 one at ∼1600 K, including both backbone and non-backbone levels.

Fig. D.1. Adopted absorption coefficient of dust normalized at 6 μm (black), compared to the mid-IR values obtained by Fritz et al. (2011) toward the Galactic Center (red).

Fig. D.2. Line profiles of some emission lines as compared with two low-excitation H2O ν2 = 1 − 0 absorption lines. The velocity scale on the abscissa is calculated with a redshift of z = 0.08778. All mid-IR molecular absorption lines have peak absorption blueshifted by 100 − 200 km s−1 relative to the peak emission lines, as exemplified by the H2O lines. Note that the rotational CO 2-1 line shows absorption below the continuum at velocities from −350 to −570 km s−1 as highlighted in the insert, consistent with the H2O wing absorption at similar velocities.

As described in Section 3.2, we include two H₂O components along the line of sight to the mid-IR continuum from the H_C: one is evenly mixed with the hot dust responsible for this continuum, with velocities between 75 and 170 km s⁻¹, and the other is a broad layer extending from the H_C photosphere outwards with velocities between 150 and 400 km s⁻¹. In combination with a turbulent velocity of 90 km s⁻¹, this velocity field can approximately match the observed profiles. Radiative pumping fully dominates the excitation of the first component, so that the H₂ density used in the models (n_H2 = N_H2/(R_h − R_int)∼6 × 10⁵ cm⁻³) has no impact on the results as long as the H₂O abundance relative to the dust remains unchanged. The foreground layer is used to account for the low-excitation absorption at the highest velocities, and is shown to produce significant absorption up to the 5_0, 5 − 4_1, 4 line (E_low ≈ 300 K). In higher excitation lines, the blueshifted absorption is restricted to velocities >  − 300 km s⁻¹ (e.g., 7_1, 6 − 6_2, 5). Because of the strong nuclear mid-IR radiation field, the excitation of this layer component is also dominated by radiative pumping.

Despite the large column densities of H₂O mixed with the dust in the H_C, the predicted absorption troughs are still sensitive to the H₂O abundance relative to the dust, mostly in the highest excitation lines. This is because the dust opacity restricts the H₂O that can generate absorption to the external layers of the source (i.e., to τ_6μm ≲ 1.5 from the surface), and N(H₂O)∼5 × 10¹⁸ cm⁻² in this shell is much lower than the total value. Our best-fit H₂O abundance relative to H₂ of 8 × 10⁻⁵ is similar to the value inferred by García-Bernete et al. (2025b).

The profiles displayed in Fig. D.2 show that our best fit model overpredicts the absorption in the low excitation 1_1, 0 − 1_0, 1 and 1_0, 1 − 1_1, 0 lines. The discrepancy is attributable to potential resonant scattering of the mid-IR continuum over larger scales, as H₂O is also present in the W_C (Appendix E). The model also overpredicts the absorption in some high excitation lines (9_0, 9 − 8_1, 8, 10_1, 10 − 9_0, 9, 9_1, 8 − 8_2, 7, 7_3, 5 − 6_2, 4), but accounts rather well for the rest of the lines.

Appendix E: The far-IR molecular absorption

Due to the dilution of the W_C + H_C far-IR continuum discussed in Section 4, we show in Fig. E.1 the continuum-subtracted far-IR spectra observed with Herschel/PACS around the wavelengths of the OH 119, 79, and 65 μm doublets. If the far-IR continuum level predicted by our model for the W_C is correct, we will be able to match the observed absorption troughs in this absolute scale.

The far-IR absorption features shown in Fig. E.1 are characteristic of very optically thick regions. In panel a, the ground-state OH 119 μm doublet exhibits a flat profile with no dip in the absorption between the two l−doubling components. We also find clear detection of the red component of ¹⁸OH 120.15 μm, as well as strong absorption due to CH⁺ 3-2 at 119.8 μm and to the CH (N = 3, J = 7/2)−(2, 5/2) doublet around 118.5 μm (E_low ≈ 100 K; for a level diagram of CH see Davidson et al. 2001). To our knowledge, this is the only ULIRG where the CH⁺ 3-2 absorption is nearly as strong as that due to OH 119 μm. In panel b, two H₂O lines (4_2, 3 − 3_1, 2 and 6_1, 5 − 5_2, 4 at 78.7-78.9 μm) are detected, together with a broad feature that is in part due to the ground-state OH 79 μm doublet but requires additional absorption at velocities more blueshifted than −350 km s⁻¹ from OH. In panel c, strong absorption in the OH 65 μm doublet (E_low ≈ 300 K) is seen, which has the second highest equivalent width among all galaxies observed with Herschel/PACS (González-Alfonso et al. 2015). ¹⁸OH at 65.6 μm is also detected. As indicated in Fig. E.1 for the OH lines, all features are blueshifted by ∼150 km s⁻¹.

Fig. E.1. H2O ν2 = 1 − 0 line profiles in IRAS 07251−0248 observed with JWST MIRI/MRS. Filled histograms are the observed profiles and blue lines are model results. Abscissa represents velocity (in km s−1) and ordinate is continuum-normalized. All detected lines have a peak absorption that is blueshifted by 100 − 200 km s−1. The selected transitions are R or Q branch lines that are not extremely blended and cover a representative range of level energies. Note that some displayed lines are not detected (right column). The energy level diagram corresponds to the ground state, with levels involved in the displayed profiles shown in red.

Our model for the far-IR lines includes the quoted species (OH, ¹⁸OH, H₂O, CH, CH⁺) together with H₂O⁺, which has several rotational lines in these spectral regions. Spectroscopic data used for OH, CH, CH⁺, and H₂O⁺ were obtained from the JPL (Pickett et al. 1998) and CDMS (Müller et al. 2001, 2005) catalogs. Since the W_C is very optically thick (τ_100μm ≈ 10, Fig. 2), warm (T_d = 123 K, Table 2), and has a size (R_w = 71 pc) larger than the H_C, it dominates the absorption in the observed lines. The inferred abundances are listed in Table E.1. An OH abundance of ≈5 × 10⁻⁶, similar to the value found in the nuclear region of NGC 4418 (with similar T_d, González-Alfonso et al. 2012) is required to generate absorption troughs in the OH 119 and 65 μm doublets comparable to the observed features (light-blue curves in Fig. E.1). However, the observed very broad OH 79 μm feature is still underpredicted with just the W_C, and we have thus included in the model the expected OH contribution from the H_C, as well as a model for H₂O⁺ in both components. The 79.1 μm absorption is then better reproduced, although its unpredicted blueshifted shoulder at 79.0 μm suggests that an additional unidentified species further contributes to the absorption. With a similar abundance of ≈5 × 10⁻⁶ in the W_C, the H₂O lines are matched as well, thus indicating a decline of its abundance of ∼1 dex from the compact H_C to the more extended W_C. We also note the high abundances inferred for both CH and CH⁺, suggesting that the hydrocarbon-rich chemistry observed in the mid-IR towards the nuclear region (García-Bernete et al. 2025b) applies to some extent to the W_C.

In summary, the reasonable match to the far-IR absorption lines indicates that the continuum level of the W_C in the far-IR predicted by our model is accurate within ∼20%. The similar kinematics found for the far-IR and mid-IR absorption lines and the high abundances found for CH and CH⁺ indicate that the physical and chemical processes in the W_C are also affected by the AGN at the core of the H_C.

Fig. E.2. Herschel/PACS continuum-subtracted spectra of IRAS 07251−0248 around 119, 79, and 65 μm (filled histograms), and model fit. The upper scale in each panel is velocity in km s−1 relative to the blue component of each OH doublet, and the two vertical dotted lines indicate the positions of the two OH l−doubling components (rest wavelengths are calculated for a redshift z = 0.08778). The contribution to the modeled spectrum by different species and components is shown with different colors as labeled in panel (a); red is total.

Appendix F: The HCN, C₂H₂, CH₄, and CO₂ bands

The JWST MIRI/MRS spectrum of IRAS 07251−0248 at 6.8 − 8.0 μm in Fig. F.1a shows a forest of absorption lines that are identified with the gas-phase HCN 2ν₂ 7.1 μm, C₂H₂ν₄ + ν₅ 7.5 μm, and CH₄ν₄ 7.7 μm bands (García-Bernete et al. 2025b). Spectroscopic parameters for these species were obtained from the HITRAN2020 (Gordon et al. 2022) and ExoMol (Chubb et al. 2020; Yurchenko et al. 2024) databases. The strengths of the individual lines (panel b) are comparable to those of the H₂O ν₂ lines, which is surprising given that these bands are intrinsically weak. For the HCN 2ν₂ and C₂H₂ν₄ + ν₅ bands, the Einstein coefficients for absorption of radiation (B_lu) are ∼10⁶ and ∼2.5 × 10⁶ cm² erg⁻¹ s⁻¹, respectively, while the B_lu values for the backbone-backbone transitions of the H₂O ν₂ lines are factors of ∼5 and ∼2 higher.

We have used our model for the H_C, which dominates the continuum at these wavelengths, to fit these bands, with the goals of checking the excitation predicted by the continuum model and deriving the abundances of these species. Our model uses, as for H₂O, a λ-dependent covering factor f_hot as defined in Eq. (4) (Section 3.2). The dilution of the molecular absorption at 7 − 8 μm is also due to the “hidden” 230 K component (Fig. 2). The best-fit model is shown in Fig. F.1b. To explain the band strengths we require extraordinary abundances relative to H₂O, indicating an extreme C-rich chemistry: [HCN]/[H₂O] ≈ 0.56, [C₂H₂]/[H₂O] ≈ 0.54, and [CH₄]/[H₂O] ≈ 0.22. These are within a factor of 2 consistent with the values inferred by García-Bernete et al. (2025b).

The HCN ν₂ 14 μm and C₂H₂ν₅ 13.7 μm bands are also clearly detected (Fig. F.1c). These bands have B_lu values higher than the HCN and C₂H₂ 7.1 − 7.5 μm bands by factors of 5 − 10, and thus the P and R branch lines from the H_C at 14 μm, predicted by our model for the ∼7.5 μm bands, would be expected to be strong. However, the W_C continuum emission at 14 μm, here dominated by the W_C 123 K component (Table 2), is expected to account for about half of the observed continuum (Figs. 2 and F.1c) and thus the 14 μm absorption from the H_C is diluted with f_hot ∼ 0.5 (Fig. F.1d). This dilution has the effect of reducing the predicted 14 μm high-J molecular absorption troughs at levels comparable with the observations. Using a common model for the H_C, our composite continuum model thus yields consistent results for the 7.5 and 14 μm bands, which is required as the high-J lines of the R and P branches in both bands are produced by the same gas.

However, the Q branches of both 14 μm bands are saturated, and the model for the H_C alone cannot reproduce them as they absorb more than 50% of the continuum. We have therefore added a model for the W_C in Fig. F.1d that accounts for the remaining Q branch absorption, using a covering factor for the W_C of f_warm = 1 − f_hot (note that the C_C contributes negligibly to the 14 μm continuum, Fig. 2). The absorbing column densities of both HCN and C₂H₂ in the W_C are 10¹⁷ cm⁻², generating absorption in the Q branches and in the low-J lines of the P and R branches. The low-J lines of the HCN 14 μm band have an additional contribution from the low-excitation shell in front of the H_C, and indeed absorption is seen up to ∼ − 500 km s⁻¹ in the uncontaminated HCN ν₂ R4 line (this has the effect of lowering the apparent T_d). The fit to the 14 μm bands tends to overestimate the absorption in high-excitation lines, most likely reflecting the uncertainties of f_hot.

From the results for the HCN abundance found here and by García-Bernete et al. (2025b), and those for the HCO⁺ abundance by Speranza et al. (2025), a very high [HCN]/[HCO⁺] abundance ratio of ∼100 is obtained in IRAS 07251−0248 E. This result can be related to the variety of HCN/HCO⁺ rotational line ratios observed in galaxies, and widely discussed in the literature (e.g., Krips et al. 2008; Imanishi et al. 2009; Izumi et al. 2016). While a high HCN/HCO⁺ line ratio was first attributed to the chemistry involved in X-ray dominated regions around AGNs, further observations showing high ratios in non-AGN galaxies and low ratios in AGN have dismissed this relationship (e.g., Costagliola et al. 2011; Privon et al. 2020). Recently, Nishimura et al. (2024) have identified the highest HCN/HCO⁺ line ratios in shell-like structures kinematically associated with outflows or inflows, and their chemical models indeed indicate a strong enhancement of [HCN]/[HCO⁺] in shock conditions. The nuclear conditions in IRAS 07251−0248 E are more extreme than those modeled by Nishimura et al. (2024), which are appropriate for the more extended regions probed by the rotational lines, but their shock models illustrate the high [HCN]/[HCO⁺] ratio that can be attained in the outflow associated with the H_C. Nevertheless, the additional effect of cosmic rays in enhancing the ionization rate and processing the carbonaceous grains and PAHs are most likely required to explain the extreme [HCN]/[H₂] ∼ 4 × 10⁻⁵ we find in the nucleus of IRAS 07251−0248 E.

Fig. F.1. HCN, C2H2, and CH4 bands in IRAS 07251−0248 observed with JWST MIRI/MRS, and model fit. Abscissae are rest wavelengths in μm. Panels a and c show the observed spectra around 7.5 and 14 μm (flux densities in mJy), with the orange line indicating the adopted continuum level (i.e., the baseline) and the light-blue curve showing the contribution to the continuum by the WC. Panels b and d show the continuum normalized spectrum, with the red line showing the model result and the dashed line indicating the covering factor of the hot component (fhot). The contributions to the model by HCN, C2H2, CH4, and CO2 (and their transition labels) are indicated in light-blue, magenta, green, and brown, respectively (and vertically shifted for clarity). Note that all absorption features are blueshifted relative to the labels.

Fig. F.2. CO2ν3 and ν2 bands in IRAS 07251−0248 observed with JWST NIRSpec, and model fit. The upper panel a shows the observed spectrum around 4.25 μm (flux densities in mJy), with the orange line indicating the adopted continuum level (i.e., the baseline) and the light-blue curve showing the contribution to the continuum by the WC. The lower panel shows the continuum normalized spectrum, with the blue line showing the model result for the HC and the dashed line indicating the covering factor of the hot component (fhot). We have also included the required contribution to the absorption by CO2 ice (exp{−τice} in magenta); red is the resulting spectrum. The insert shows the Q branch of the CO2ν2 band at ≈15 μm, with the blue line also indicating the model result for the HC.

Finally, we have also applied the H_C model to the CO₂ν₃ (4.25 μm) and ν₂ (14.97 μm) bands (Fig. F.2). The best fit model to both the CO₂ν₃ R branch and the ν₂ Q branch requires [CO₂]/[H₂O] ≈ 0.03 at the H_C photosphere, with an additional contribution by the surrounding outflowing shell with [CO₂]/[H₂O] ≈ 0.014. The fit, however, grossly underestimates the strength of the CO₂ν₃ P branch up to J ∼ 20, so we added foreground CO₂ ice with H₂O : CO₂ = 10 : 1 (Ehrenfreund et al. 1997; Rocha et al. 2022). The resulting modeled spectrum appears to overestimate the absorption in the high-J lines of the ν₃ P branch, but fits both the ν₃ R branch at 4.23 μm and the ν₂ Q branch at 15 μm rather well. This match validates our fit to the H_C continuum.

All Tables

All Figures

	Fig. 1. ALMA 667 μm continuum emission of IRAS 07251−0248 E and model. Left: Observed emission, with colored scale in units of mJy/beam. Middle: Model, composed of an unresolved source (full-width at half maximum, FWHM < 21 mas = 34 pc) and a Gaussian source (FWHM = 89 × 86 mas2 = 146 × 141 pc2). Right: Residuals in μJy/beam.
In the text

Fig. 2. Spectral energy distribution (SED) of IRAS 07251−0248 E from λrest = 2.6 to 1200 μm and model. In black, the full JWST NIRSpec & MIRI/MRS spectrum is shown, together with Spitzer, Herschel/PACS (from 60 to 160 μm, both photometric and spectroscopic) and SPIRE (from 230 to 470 μm), and ALMA (667 and 1200 μm) data. The colored circles at 667 μm indicate the contributions by the unresolved (blue) and resolved (light-blue) components inferred from the model in Fig. 1. At 1.2 mm, the open black circle and the cross indicate the flux densities of the unresolved (< 82 pc) and resolved (540 × 330 pc2) components, respectively (Pereira-Santaella et al. 2021). The model includes the HC in blue (dashed: unattenuated emission; solid: attenuated), the WC in light-blue (extincted), a cold component (CC) in orange, the stellar component in gray, and the “hidden” component (whose attenuated emission is added to that of the WC, see Section 3.2) in dotted-dashed green. Red is total. The insert panels show details of the model results. Upper left: The Tdust profile of the HC. Upper right: The optical depths τ of the HC and the WC as function of wavelength. Lower middle: The Tdust profile of the HC as function of the optical depth from surface at 7 and 14 μm.

In the text

Fig. 3. H2O ν2 = 1 − 0 band in IRAS 07251−0248 observed with JWST MIRI/MRS. Panels a and c show the observed spectrum, with the orange line indicating the adopted continuum level and the light-blue curve showing the contribution from the warm component. Panels b and d show the continuum normalized spectrum, with the blue line showing the model result for the hot component and the dashed line indicating its covering factor (fhot). We label the main transitions that potentially contribute to the H2O band as shown in González-Alfonso et al. (2024), but some very highly-excited lines are not detected in IRAS 07251−0248. Note that all H2O absorption features are blueshifted relative to the labels.

In the text

Fig. 4. Schematic representation of the model source. The extremely compact and expanding core (HC) is shown in blue (the 6 μm photosphere is its external surface), with gas flowing outward (the surrounding shell outflowing at ∼ − 400 km s−1) in green. The starburst torus shown in magenta represents the WC. The torus also appears to be expanding, because the far-IR lines are blueshifted (Appendix E, Fig. E.1). Both the HC and WC are embedded in an extended cold component (CC), in orange. The model is axisymmetric about the line of sight.

In the text

	Fig. A.1. ALMA 667 μm continuum emission of IRAS 07251−0248 E and alternative models: a point source (upper) and a Gaussian (lower). Left: Observed emission, with colored scale in units of mJy/beam. Middle: Model. Right: Residuals in μJy/beam. These models give residuals much more prominent than our best point-source + Gaussian model fit in Fig. 1.
In the text

	Fig. D.1. Adopted absorption coefficient of dust normalized at 6 μm (black), compared to the mid-IR values obtained by Fritz et al. (2011) toward the Galactic Center (red).
In the text

Fig. D.2. Line profiles of some emission lines as compared with two low-excitation H2O ν2 = 1 − 0 absorption lines. The velocity scale on the abscissa is calculated with a redshift of z = 0.08778. All mid-IR molecular absorption lines have peak absorption blueshifted by 100 − 200 km s−1 relative to the peak emission lines, as exemplified by the H2O lines. Note that the rotational CO 2-1 line shows absorption below the continuum at velocities from −350 to −570 km s−1 as highlighted in the insert, consistent with the H2O wing absorption at similar velocities.

In the text

Fig. E.1. H2O ν2 = 1 − 0 line profiles in IRAS 07251−0248 observed with JWST MIRI/MRS. Filled histograms are the observed profiles and blue lines are model results. Abscissa represents velocity (in km s−1) and ordinate is continuum-normalized. All detected lines have a peak absorption that is blueshifted by 100 − 200 km s−1. The selected transitions are R or Q branch lines that are not extremely blended and cover a representative range of level energies. Note that some displayed lines are not detected (right column). The energy level diagram corresponds to the ground state, with levels involved in the displayed profiles shown in red.

In the text

Fig. E.2. Herschel/PACS continuum-subtracted spectra of IRAS 07251−0248 around 119, 79, and 65 μm (filled histograms), and model fit. The upper scale in each panel is velocity in km s−1 relative to the blue component of each OH doublet, and the two vertical dotted lines indicate the positions of the two OH l−doubling components (rest wavelengths are calculated for a redshift z = 0.08778). The contribution to the modeled spectrum by different species and components is shown with different colors as labeled in panel (a); red is total.

In the text

Fig. F.1. HCN, C2H2, and CH4 bands in IRAS 07251−0248 observed with JWST MIRI/MRS, and model fit. Abscissae are rest wavelengths in μm. Panels a and c show the observed spectra around 7.5 and 14 μm (flux densities in mJy), with the orange line indicating the adopted continuum level (i.e., the baseline) and the light-blue curve showing the contribution to the continuum by the WC. Panels b and d show the continuum normalized spectrum, with the red line showing the model result and the dashed line indicating the covering factor of the hot component (fhot). The contributions to the model by HCN, C2H2, CH4, and CO2 (and their transition labels) are indicated in light-blue, magenta, green, and brown, respectively (and vertically shifted for clarity). Note that all absorption features are blueshifted relative to the labels.

In the text

Fig. F.2. CO2ν3 and ν2 bands in IRAS 07251−0248 observed with JWST NIRSpec, and model fit. The upper panel a shows the observed spectrum around 4.25 μm (flux densities in mJy), with the orange line indicating the adopted continuum level (i.e., the baseline) and the light-blue curve showing the contribution to the continuum by the WC. The lower panel shows the continuum normalized spectrum, with the blue line showing the model result for the HC and the dashed line indicating the covering factor of the hot component (fhot). We have also included the required contribution to the absorption by CO2 ice (exp{−τice} in magenta); red is the resulting spectrum. The insert shows the Q branch of the CO2ν2 band at ≈15 μm, with the blue line also indicating the model result for the HC.

In the text