© The Authors 2023

1. Introduction

Recent observational evidence has revealed the presence of highly compact (< 100 pc), dusty nuclei in a subset of luminous (10¹¹L_IR < L_⊙10¹² L_⊙) and ultra-luminous (L_IR > 10¹² L_⊙) infrared galaxies (LIRGs and ULIRGs, respectively; e.g., Aalto et al. 2015; Falstad et al. 2021). These objects represent a crucial phase in galaxy evolution where the nucleus builds mass rapidly through star formation and/or supermassive black hole (SMBH) growth (e.g., García-Burillo et al. 2012; Gao et al. 2021). As SMBHs are known to co-evolve with their host galaxy (e.g., Kormendy & Ho 2013), these nuclei may play a key role in this co-evolution.

Compact obscured nuclei (CONs) exhibit high column densities of gas (N_H > 10²⁵ cm⁻²) and dust that are highly opaque at short wavelengths (UV/optical). This has prevented the nature of the power source from being uncovered. It is unclear whether the nuclear power source is accretion onto a SMBH and/or a compact nuclear starburst with a top-heavy initial mass function (e.g., Hickox & Alexander 2018; Tacchella et al. 2018; Kocevski et al. 2017).

Rotational-vibrational lines of HCN have been used to identify and probe the radiation field hidden behind the dust (e.g., Aalto et al. 2015, 2019; Falstad et al. 2021). These lines are a result of a heat-trapping effect (e.g., Kaufman et al. 1998; Rolffs et al. 2011; González-Alfonso & Sakamoto 2019), where the extreme density of dust and gas traps radiation from a central heating source, increasing its internal temperature and populating the vibrational states of HCN while the external dust remains cool. Falstad et al. (2021) carried out a search for CONs in the local Universe (z < 0.08) using a representative sample of 46 galaxies, finding CONs more frequently in ULIRGs than LIRGs. However, this technique relies on a relatively weak emission feature that is extremely challenging to detect at higher redshifts, prompting the need for alternate techniques.

García-Bernete et al. (2022a) showed that polycyclic aromatic hydrocarbon (PAH) features in the mid-infrared can be used to detect highly obscured nuclei via their equivalent width (EW) ratios. The presence of a highly obscured nucleus will cause deep absorption at 9.8 μm from silicates, which will suppress the nuclear continuum around the 11.3 μm PAH feature, increasing its EW relative to the 6.2 μm and 12.7 μm PAH features, which are less affected by extinction. Therefore, highly obscured nuclei can be revealed through lower EW ratios of the 12.7 μm feature relative to the 11.3 μm and the 6.2 μm to the 11.3 μm PAH (García-Bernete et al. 2022a). Pure star-forming galaxies show near-constant PAH EW ratios as the intrinsic flux ratio is approximately constant and the PAHs are subject to the same extinction as the continuum (Hernán-Caballero et al. 2020).

Additional techniques of identifying highly obscured nuclei in the mid-infrared include the 14 μm HCN absorption feature (Lahuis et al. 2007) and the various crystalline absorption bands (e.g., Spoon et al. 2022). However, the former was difficult to observe with Spitzer, as this feature is faint and easily diluted, and the latter relies on absorption at 28 μm and 33 μm, which are outside the observing range of the James Webb Space Telescope (JWST), limiting future application.

In this work we investigate the properties of deeply obscured nuclei based on mid-infrared spectroscopy and evaluate the efficacy of different methods in identifying such objects. Utilising Spitzer InfraRed Spectrograph (IRS; Houck et al. 2004) spectra, we introduce a new spectral decomposition method that splits the continuum into a relatively unobscured star-forming component and a nuclear component, which is subject to higher extinction. We apply the new decomposition method to samples of LIRGs and ULIRGs to determine the properties of their nucleus. We also use spatially resolved spectra to test how dilution from the host galaxy affects the measured PAH EW ratios.

The paper is structured as follows: In Sect. 2 we describe the samples of the IRS staring mode data and the construction of cubes for the spectral mapping data. In Sect. 3 we describe our new spectral decomposition method. In Sect. 4 we explore the properties derived from the spectral decomposition and analyse the spectral mapping data. Finally, in Sect. 5 we estimate the number of CONs in ULIRGs and LIRGs and discuss the effect of galaxy inclination and prospects for future work using the JWST.

2. Observations

In this work we used Spitzer spectroscopy for a number of galaxy samples. The majority of the data presented are low-resolution staring mode spectra. The reduced IRS staring mode spectra were obtained from the Infrared Database of Extragalactic Observables from Spitzer (IDEOS; Spoon et al. 2022).

2.1. Staring mode data

For the present study we used two main samples of LIRGs and ULIRGs. For a representative sample of ULIRGs, we used the Herschel Ultra luminous infrared galaxy Survey (HERUS) sample (Farrah et al. 2013), which consists of the 42 local ULIRGs that were observed with Spitzer. For a complete sample of LIRGs we used the Great Observatories All-sky LIRG Survey (GOALS; Armus et al. 2009) sample, which consists of 179 LIRGs and 22 ULIRGs. As not all of these have been observed with Spitzer, this brings the sample down to 143 LIRGs and 15 ULIRGs.

We also used a reference sample of purely star-forming galaxies from Hernán-Caballero et al. (2020) where full details of the sample selection can be found. In addition we used the CONquest sample, which consists of 44 declination and distance limited objects selected from the IRAS revised bright galaxy sample (Sanders et al. 2003). Full details can be found in Falstad et al. (2021). Of these 44 galaxies, 29 have Spitzer spectra with the majority of the objects without spectra being sub-LIRGs.

2.2. Spectral mapping data

In addition to the staring mode spectra we also employed a sub-sample of LIRGs from the GOALS sample that had Spitzer IRS spectral mapping observations with the SL1 and SL2 modules, providing spectra between ∼5 and 14 μm (Alonso-Herrero et al. 2009, 2012; Pereira-Santaella et al. 2010). We chose targets without staring mode spectra with the exception of NGC 6926, which was selected as a CON candidate by García-Bernete et al. (2022a) to provide an additional candidate to investigate. We also excluded a few objects (IRAS 08339+6517, NGC 6052, NGC 6240, and NGC 6621) that had poor signal to noise ratio or issues with the data reduction. The resulting sample presented in this work is shown in Table 1 with some basic properties.

The basic calibrated data were obtained from the Spitzer Heritage Archive¹ and subsequently combined to construct spectral cubes using CUBISM. This process combines the different slit exposures, subtracts the sky background and rejects bad pixels yielding a spectral cube for each module and an associated cube containing the flux errors.

To accurately extract spectra from the two SL modules, we determined the centre location of the object of interest and spatially aligned the two modules. We use the DAOStarFinder algorithm (Stetson 1987) from the PHOTUTILS python package, to obtain positions of sources in the image. We applied this technique to integrated intensity (moment 0) maps of the 5.4 μm continuum (5.47–5.44 μm), to identify the position of the nucleus. In order to ensure the cubes are spatially aligned we generated integrated intensity maps of the overlapping channels (7.53–7.6 μm), and found the position where the flux peaks. The spatial offset between the two modules was generally found to be small but non-negligible in some cases (up to ∼1″). This offset correction was then applied when extracting the spectra from each cube.

As the point spread function for these observations is large (full width at half maximum of ∼3.8″ at 14 μm), we extracted spectra in circular apertures each of which provide different dilutions of the galaxy disc with respect to the nuclear region. Specifically, they were circular apertures of radius 2.5, 4.5, 6.5, and 8.5 arcsec, which are shown in Fig. 1 for Arp 299 A.

Fig. 1. Example of the spectral extraction, in this case for Arp 299 A, at different spatial scales. Left: integrated intensity map of the SL1 cube (7.46–14.29 μm) overlaid with circular apertures of radii 2.5, 4.5, 6.5, and 8.5″, shown as the dashed lines. Right: extracted spectra from the apertures in the left panel with the corresponding colour. The innermost spectrum (purple) has had an aperture correction applied.

The smallest aperture was corrected for slit losses using a standard star to achieve a point source extraction for the nuclear region. The correction factor was obtained by calculating the flux ratio between an aperture containing the total flux of the star (7.5″) and the nuclear aperture (2.5″). A smooth function of the correction factor as a function of wavelength was obtained by fitting a 4th order polynomial to the measured flux ratios. This factor was subsequently applied to the 2.5″ spectral extractions.

The majority of spectra analysed in this work consists only of the SL1 & SL2 modules. This is because the longer wavelength LL1 & LL2 modules probe larger spatial scales and so galaxies containing an active galactic nucleus (AGN) or a CON will have different relative contributions of the host at long wavelengths compared to short wavelengths if the full spectral range was used. Therefore, in order to avoid making physical assumptions when applying a correction factor we simply restrict our analysis to the SL modules.

3. Spectral fitting

There are a number of approaches to modelling the mid-infrared spectrum of galaxies. A popular tool from the Spitzer era is PAHFIT (Smith et al. 2007), which models the continuum with a series of black bodies subject to extinction and PAH emission features with a series of Drude profiles. The tool is well suited to model pure star-forming galaxies; however, it is less well suited to objects that contain highly obscured nuclei and/or AGN. Alternative techniques have been used to address this problem, such as QUESTFIT (Veilleux et al. 2009) and DeblendIRS (Hernán-Caballero et al. 2015); however, both of these assume a set of fixed templates for the PAH emission, which prevents us from measuring the properties of the PAHs independently for each object.

The presence of a distinctly different nuclear continuum in heavily obscured objects makes it difficult to model these sources with a single extinction component. We therefore created a new model that retains the PAH emission and spectral lines of the PAHFIT tool but models the continuum with two components, a nuclear and a star-forming one, each subjected to different levels of extinction. In what follows we describe the parameterisation of each component in detail.

3.1. Star-forming component

To model the star-forming component of the continuum we generated templates by fitting the star-forming sample of (Hernán-Caballero et al. 2020) with PAHFIT. The left panel of Fig. 2 shows the resulting continuum templates ordered by their flux ratio between 8 and 13 μm. Rather than selecting an individual template, we linearly interpolate between the templates ordered in this way with a single parameter between 0 and 1 where 0 corresponds to the steepest template. The resulting template, C_SF(λ), is then subject to extinction (assuming mixed geometry) of the form $\frac{1-{\mathrm{e}}^{-{\tau }_{\lambda }}}{{\tau }_{\lambda }}$, where we restrict the optical depth, τ_λ, using a half-normal prior based on the resulting fits from the star-forming sample. This is shown in the right panel of Fig. 2, where the histogram shows the optical depths at 9.8 μm inferred from the star-forming sample and the red line shows the half-normal prior used for the extinction of the star-forming component for the spectral decomposition.

Fig. 2. Star-forming continuum templates derived from the star-forming sample. Left: unobscured continua, CSF(λ), obtained from the fits to the star-forming sample, used as templates for the spectral decomposition. The colour indicates the steepness (8.0 μm–13 μm): the darker the colour, the steeper the continuum. Right: histogram of the optical depths, assuming mixed geometry, for the star-forming sample. The red line shows the half-normal prior used for the spectral decomposition based on these results.

3.2. Nuclear component

The nuclear component in our model consists of a continuum with some extinction applied. The continuum, C_Nuc(λ), is modelled as a quadratic spline with three evenly spaced anchor points at 5.5 μm, 9.85 μm and 14.2 μm. The y-value of these knots is allowed to vary as a free parameter. As the continuum is normalised by the total flux and scaled by the nuclear fraction, β, the spline has effectively only two free parameters. Extinction is then applied to the continuum as a screen using the ice template and silicate template derived from NGC 4418 and has the form e^{−(τ_Ice(λ)+τ_N(λ))}. The silicate template is normalised to 1 at 9.8 μm, and therefore the optical depth, τ_N, is that measured at 9.8 μm. However, it is worth noting that this is an underestimation as it assumes τ = 0 at the anchor points when in reality a power law component to the extinction curve causes τ > 0 at the anchor points.

The ice template (including the CH absorption) is allowed to vary separately to the silicate template, resulting in a total of four free parameters (two for the shape of the spline, C_Nuc(λ), and two for the extinction, τ_Ice and τ_N) determining the shape of the nuclear component. The creation of these templates is discussed in Appendix A. For the silicate profile we used a template derived from a highly obscured source, namely NGC 4418, because the composition of silicates is likely different (resulting in a different profile) compared to extinction laws derived from the Milky Way (e.g., Tsuchikawa et al. 2021) or the line of sight is more complex than towards the galactic centre.

It is also work noting that the optical depth peak of the silicate profile of NGC 4418 is ∼4.4 from the Spitzer spectroscopy whereas ground based observations with a smaller beam find a higher peak ∼7 (Roche et al. 2015). This means that the silicate profile used is still contaminated by some amount of star formation and so the nuclear optical depth values inferred should be taken as lower limits for the true optical depth. In Appendix B we compared this silicate template to another one – that of IRAS 08572+3915. From this testing we find the results from this work to be unchanged.

For objects where the nuclear optical depth is minimal, the continuum decomposition is unreliable as there is nothing to differentiate the star-forming component from the nuclear one. However, the total continuum should be reliable. Therefore, the decomposition is particularly sensitive to cases where the optical depth is high (as in the case of a CON) and the contribution of the nuclear component is above a certain threshold. We investigate this further in Sect. 3.5.

3.3. Full model

The resulting full model is given by

$$\begin{array}{cc}\hfill f_{\nu }(\lambda )& ={\displaystyle [\sum _{i=1}^{N_{\mathrm{Lines}}}{I}_{\nu ,\mathrm{L}ine}^{(i)}(\lambda )+\sum _{i=1}^{N_{\mathrm{PAH}}}{I}_{\nu ,\mathrm{P}AH}^{(i)}(\lambda )+(1-\beta )C_{\mathrm{SF}}(\lambda )]\frac{1-{\mathrm{e}}^{-{\tau }_{\lambda }}}{{\tau }_{\lambda }}}\hfill \\ \hfill & +\beta C_{\mathrm{Nuc}}(\lambda ){\mathrm{e}}^{-({\tau }_{\mathrm{Ice}}(\lambda )+{\tau }_N(\lambda ))},\hfill \end{array}$$(1)

where C_SF(λ) is the star-forming unobscured continuum, β is the nuclear fraction, and C_Nuc(λ) is the nuclear unobscured continuum. These unobscured continua multiplied by their respective extinction factors are normalised by their area to ensure that the scale factor, β, gives the fractional contribution of the nucleus to the total continuum flux (between 5 and 14 μm). The area is calculated numerically during the fitting process.

To aid in the decomposition, we enforced a prior on the ratio of the total PAH flux, f_PAH, to the integrated continuum flux of the star-forming continuum, f_SF, between 5 and 14 μm. From the PAHFIT results for the star-forming calibration sample we found this ratio to be f_PAH/f_SF ≈ 1.92. We use a very wide normal prior with a standard deviation of 10.0. This wide prior is designed to discourage solutions with zero contribution from the star-forming component as the large apertures of Spitzer will contain at least some extended contribution. The large width of this prior was chosen arbitrarily such that the prior discourages zero star-forming component where $f_{\mathrm{PAH}}/f_{\mathrm{SF}}\stackrel{\infty }{\to }$ but does not heavily bias the results.

3.4. Bayesian inference

To fit these models to the spectra, we used a Bayesian approach to provide accurate posterior probabilities for the galaxy properties of interest. We used Markov chain Monte Carlo (MCMC) sampling from NUMPYRO (Phan et al. 2019) to sample the relative posterior probability of the model given the data, Pr(M|D), which is given by

$$\begin{array}{c}\hfill ln\text{Pr}(M|D)=ln\text{Pr}(M)+ln\text{Pr}(D|M)+\mathrm{const}.,\end{array}$$(2)

where the prior, Pr(M), is uniform between sensible limits as discussed in the previous sections and the log-likelihood, lnPr(D|M), is effectively just the chi-squared as the error bars are fixed by the data:

$$\begin{array}{c}\hfill ln\text{Pr}(D|M)={\displaystyle \sum _{i=1}^N\left(\frac{{(f_i-f_{\nu }(\lambda ))}^2}{{\sigma }_i^2}\right)+\mathrm{const}.,}\end{array}$$(3)

where there are N flux data f_i, with error σ_i. From NUMPYRO, we used the No-U-Turn Sampling (NUTS), a Hamiltonian Monte Carlo method that allows for faster sampling through parallelisation compared to traditional MCMC methods. After a burn-in of 2000 samples we take another 2000 samples to obtain marginalised posteriors for each of the parameters.

From the MCMC samples, posterior probability distributions for the properties derived from the spectra were constructed. We are primarily interested in the EW of the PAH features, which were calculated numerically for each sample of the parameters using

$$\begin{array}{c}\hfill \mathrm{EW}={\displaystyle \int }\frac{f_{\nu }^{\mathrm{PAH}}}{f_{\nu }^{\mathrm{cont}}}\mathrm{d}\lambda ,\end{array}$$(4)

where $f_{\nu }^{\mathrm{PAH}}$ is the PAH profile and $f_{\nu }^{\mathrm{cont}}$ is the continuum. The EW ratios of the various PAH features can then be calculated with associated error bars. The continuum used is the ice-corrected continuum to allow comparisons with the IDEOS sample (Spoon et al. 2022), and to better determine how the silicates may affect the PAH features. The integrated flux of the PAH features is also calculated using

$$\begin{array}{c}\hfill f^{\mathrm{PAH}}={\displaystyle \int }{f}_{\nu }^{\mathrm{PAH}}\mathrm{d}\nu ,\end{array}$$(5)

where posteriors for the PAH flux ratios can be calculated. We do this with and without extinction applied to obtain estimates for the intrinsic ratios as well.

We also estimate the strength of the 9.8 μm silicate feature defined as (Spoon et al. 2007)

$$\begin{array}{c}\hfill S_{\mathrm{sil}}=ln\left(\frac{f_{9.8,\mathrm{obs}}^{\mathrm{cont}}}{f_{9.8,\mathrm{int}}^{\mathrm{cont}}}\right),\end{array}$$(6)

which is the log of the ratio of the observed to the intrinsic continuum at 9.8 μm. It should be noted that this is different from the nuclear optical depth, τ_N, in the spectral decomposition method, which allows the nuclear silicate feature to be ‘filled up’ by the star-forming dust continuum. Therefore, objects with a buried nucleus but with a strong extended star-forming component may have a low observed silicate depth as defined by Eq. (6) but a deep nuclear optical depth, τ_N.

3.5. Testing the method

The primary motivation for modelling the spectra using this new decomposition tool is to provide a realistic continuum to enable us to accurately infer the properties of the emission features. However, it also allows us to better constrain the physical properties of the two components (nuclear and star-forming), particularly the optical depth of the nuclear region, the fraction of the nuclear to the total continuum and, consequently, the shape of the nuclear continuum. It is therefore instructive to test the effectiveness of this decomposition method before inferring the properties of the deeply obscured galaxies studied here.

3.5.1. Accuracy of inferred properties

We first generated simulated data of a typical star-forming galaxy hosting a CON (with a given τ_N) with varying degrees of dilution from the host galaxy (different values of β). This was done by constructing a spectrum using a star-forming galaxy and a CON template, where different fractional contributions, β, of the CON were used. We then ran the decomposition tool on the spectra and tested its ability to recover the true nuclear fraction and nuclear optical depth.

We used NGC 1797 for the star-forming component as this shows a typical spectrum of a star-forming galaxy with some silicate absorption, and the nucleus of NGC 4418 as the CON template. The star-forming continuum was normalised to the integrated continuum flux between 5 and 14 μm. The PAH and emission lines from NGC 1797 were also included as part of the star-forming component and were normalised by the same factor as the continuum. The CON continuum from NGC 4418 was normalised by its integrated flux over the same wavelength range. The star-forming component was then scaled by (1 − β) and the CON component by β for 20 values of β between 0 and 1. We generated two sets of simulated data, the first with a CON component with τ_N = 4.5 and the second with a lower nuclear optical depth of τ_N = 3.0 to test the reliability of the method for different levels of nuclear obscuration.

Figure 3 shows two example fits to the simulated data with the left panel showing a CON highly diluted by the host galaxy (β = 0.42) and the right shows a CON-dominated source (β = 0.89). In both cases the nuclear optical depth, τ_N, and nuclear fraction, β, are successfully recovered by the model. At values of β ≲ 0.4, we find that the method is unable to recover the nuclear fraction and/or the optical depth as reliably, although in the case of the optical depth, the errors do a reasonably good job of accounting for this. This is demonstrated in Fig. 4, where the measured nuclear fraction is plotted against the true value of the nuclear fraction in the left panel. The red points shown the data set with a τ_N = 4.4 whereas the blue show a lower τ_N = 3.0. This plot demonstrates that even when the nuclear component has a lower optical depth, the nuclear component is still recovered by the model. The right panel shows the measured optical depth for each simulated spectrum with the true value of τ_N = 4.5 displayed for the red points and τ_N = 3.0 for the blue points.

Fig. 3. Example of our spectral decomposition fit (described in Sect. 3) to simulated data of a galaxy that hosts a CON. The left panel shows a highly diluted nucleus (β = 0.42) with a spectrum dominated by star formation, and the right panel shows a CON-dominated source (β = 0.89). The full model, shown in red, is composed of the various components shown in the legend. The continuum is shown as the solid gold line and is a sum of the nuclear and star-forming (SF) components, which are shown as dashed lines. The flux is in arbitrary units. The equation of the full model is given in Eq. (1).

Fig. 4. Testing the recovered properties of simulated data using the spectral decomposition fitting method. Left: measured nuclear fraction, β, against the true value. The solid black line shows where the measured equals the true value. Right: measured nuclear optical depth, τN, against the true nuclear fraction. The solid horizontal lines show the true optical depth of τN = 4.5 in red and τN = 3.0 in blue. The vertical dashed black lines show a cutoff of β = 0.4, below which the measured values become unreliable.

In the cases of maximum dilution β ≲ 0.15, the optical depth recovered is significantly below the true value but with a high measured value of β. This is likely due to the nuclear component being used to fit the star-forming component as the two components become degenerate when β is low. We therefore conclude that the optical depth is accurately recovered for values of τ_N > 2.5 (> 5σ of the silicate depths of the star-forming calibration sample) and nuclear fractions of β > 0.4.

3.5.2. Multiple-aperture fitting

From the spectral mapping data, each aperture yields a spectrum (see Fig. 1) that we would like to fit. As each of these apertures has the same centre but different radii they contain the same nucleus but different star-forming fractions. It is therefore important to fit all these spectra simultaneously assuming the same nuclear component for all the spectra but different star-forming components for each.

To test this idea, we used the spectral map of a known CON to measure the nuclear properties from the spectra of each aperture individually and then compare the results to those from fitting all the apertures simultaneously. We used ESO 320-G030, which is known to be a CON but has a strong star-forming component (e.g., Alonso-Herrero et al. 2006; González-Alfonso et al. 2021).

The left panel of Fig. 5 compares the measured nuclear fraction for each aperture when fitting all the spectra simultaneously sharing the same nucleus against fitting each of the spectra independently. The measured values are consistent within 1σ between the two methods. The right panel of Fig. 5 shows the measured optical depth of the nucleus measured from each aperture individually compared to those measured from all apertures when fit simultaneously. Again, we find that the measurements are consistent between the two methods and crucially, even in the presence of significant contribution from the host galaxy. In this case we find that the value of τ_N does not change, instead the error bars expand accordingly. This test also demonstrates the added value in fitting all the apertures simultaneously, with tighter constraints on the nuclear optical depth. Having proved the robustness of our method, we proceed and deploy the tool to fit the spectral mapping data.

Fig. 5. Comparison between fitting spectra from individual apertures vs. fitting all apertures simultaneously with shared nuclear parameters for ESO 320-G030. Left: nuclear fraction, β, measured from the simultaneous fitting on the y-axis vs. the individual spectra on the x-axis. The solid black line shows where they are equal. Right: measured optical depth of the nucleus for each aperture individually and for all apertures simultaneously. Note the much smaller error when using all the spectra simultaneously.

4. Results

Based on the spectral fits using the new decomposition tool, we first refine the CON selection criteria (Sect. 4.1) before investigating the physical properties of the nucleus in Sects. 4.2 and 4.3. We then present results using the spectral mapping data in Sect. 4.4.

4.1. Refining the CON selection criteria

To select CONs using PAH EW ratios, García-Bernete et al. (2022a) defined two ratios: EW(12.7)/EW(11.3) and EW(6.2)/EW(11.3). Both are the ratio of the EW of a PAH feature outside the 9.8 μm silicate absorption (12.7 μm and 6.2 μm) relative to the 11.3 μm PAH that is located within this absorption band. Using the star-forming sample of Hernán-Caballero et al. (2020); García-Bernete et al. (2022a) showed that both of these ratios were constant for the sample although the 6.2/11.3 ratio shows a larger scatter. As the presence of a heavily obscured nucleus results in a lower continuum around the 11.3 μm PAH feature compared to the 6.2 μm or 12.7 μm PAH, the CON selection region was defined by values of the EW(12.7)/EW(11.3) and EW(6.2)/EW(11.3) ratios below those found for the star-forming sample.

As the previous PAH EW thresholds for selecting highly obscured nuclei from García-Bernete et al. (2022a) were based on IDEOS fits using Pearson-IV profiles for the PAH features, we need to refine those criteria in line with our new fitting method. For this purpose we therefore re-fitted the Hernán-Caballero et al. (2020) star-forming sample using our new spectral decomposition method.

The resulting EW ratios are shown in the top left panel of Fig. 6 where the 12.7/11.3 PAH EW ratio has a mean value of ⟨EW_12.7/EW_11.3⟩ = 0.427 with a scatter of σ(EW_12.7/EW_11.3) = 0.077. Subtracting the average uncertainty of the data (0.030) in quadrature, we obtain an intrinsic scatter of ∼17%. The 6.2/11.3 PAH EW ratio shows a much larger scatter with a mean ⟨EW_6.2/EW_11.3⟩ = 0.99 and a scatter σ(EW_6.2/EW_11.3) = 0.49.

Fig. 6. PAH EW diagram (García-Bernete et al. 2022a) of the 12.7/11.3 PAH features against the 6.2/11.3 PAH features obtained from fitting the star-forming sample (top left; Hernán-Caballero et al. 2020), CONquest (top right; Falstad et al. 2021), HERUS (bottom left; Farrah et al. 2013), and GOALS (bottom right; Armus et al. 2009). Points in red have a nuclear optical depth τN > 3.5, and for the CONquest sample the triangle points are those identified as CONs by ΣHCN-vib > 1 L⊙ pc−2. The grey lines show the mean values for the star-forming sample, with the 2σ values for the 12.7/11.3 PAH shown by the dashed lines. The purple shaded area shows the CON selection region.

Comparing our fitting (Drude PAH profiles) to that of IDEOS (Pearson-IV PAH profiles), we find a larger mean value in the EW_12.7/EW_11.3 and a higher scatter. They found a mean of 0.346 and ∼5% intrinsic scatter (Hernán-Caballero et al. 2020).

The larger scatter in the 6.2/11.3 PAH EW ratio may be due to differences in the intrinsic properties of the continuum. Different dust temperatures will more strongly affect the continuum between 6.2 and 11.3 than 12.7 and 11.3. This can be seen in the left panel of Fig. 2 where the presence of hotter dust (> 200 K) will increase the continuum around 6.2 μm.

Using the newly measured PAH EW values we defined the CON selection region as follows: we take the 2σ lower boundary as the threshold for the 12.7/11.3 and the 1σ lower limit for the 6.2/11.3. These criteria result in the purple shaded region shown in Fig. 6, where the thresholds are EW_12.7/EW_11.3 = 0.271 and EW_6.2/EW_11.3 = 0.503 for our spectral decomposition tool.

Two objects in the Hernán-Caballero et al. (2020) star-forming sample fall in the CON region namely VV 283a and UGC 01845. Both of these objects are LIRGs with log(L_IR/L_⊙) = 11.68 and 11.12, respectively (Armus et al. 2009). These two objects are not selected by García-Bernete et al. (2022a). From our optical depth selection discussed in Sect. 4.2 we identify a further 5 potential CON candidates from this sample (all of them are LIRGs). Since this sample was selected only by excluding AGN, identifying CON candidates is therefore not unexpected.

4.2. Physical properties of CONs

Fitting the mid-infrared spectra with our physically motivated model allows us to recover the continuum of the nuclear component, thus providing insight into the physics of CONs, in particular the nuclear optical depth at 9.8 μm, τ_N, and the level of dilution from the host galaxy in form of the nuclear fraction, β.

It is therefore instructive to investigate further properties of the CONs. First, we examined whether there is a possible relation between the model-derived quantities and the properties of the HCN-vib line. Using our method to fit the CONquest sample from Falstad et al. (2021), we measured the nuclear optical depth, τ_N. The measured PAH EW ratios are shown in the top right panel of Fig. 6. In Fig. 7 we show τ_N as a function of the HCN surface brightness in the left panel, and the ratio to the L_IR in the right panel. In both plots there is a clear trend for a higher optical depth correlating with stronger HCN-vib emission. A Pearson correlation test results in coefficients of 0.91 and 0.84 for the HCN-vib surface density plot and HCN-vib to L_IR plot, respectively. This correlation provides good evidence that the heat-trapping effect required to populate the vibrational states of HCN is consistent with the presence of large quantities of dust with a high internal temperature producing the mid-infrared radiation field and a cool dusty exterior to provide the high optical depths of the silicate absorption feature.

Fig. 7. Optical depth of the nuclear component from the spectral decomposition fits to the CONquest sample (Falstad et al. 2021) against the surface brightness of HCN-vib (left) and the strength of HCN-vib to LIR (right). Linear fits are shown in the dashed black lines, with 1σ confidence intervals shown in grey. The upper limits of HCN-vib measurements are shown for objects with τN > 2.5 and β > 0.4. The dashed grey lines show the CON identification threshold from Falstad et al. (2021) on the x-axis and the corresponding τN on the y-axis.

Unlike the nuclear optical depth, the apparent silicate strength as measured using Eq. (6) shows no trend with HCN-vib as this property is highly sensitive to dilution from the host-galaxy. We show this in Appendix C for reference. We also show the 12.7/11.3 PAH EW ratio against the strength of HCN-vib, which does show a trend, albeit one that is weaker than that of the nuclear optical depth.

The Falstad et al. (2021) definition of a CON is based on the strength of the mm vibrational transitions of HCN. They used two definitions, the surface brightness of HCN: Σ_HCN-vib > 1 L_⊙ pc⁻² and the ratio to the infrared luminosity: L_HCN-vib/L_IR > 10⁻⁸. Based on the correlation found in Fig. 7 we can therefore attempt to derive a new CON definition based on the optical depth values determined through our fitting method.

To do this we fitted a straight line in log-log space using SCIPY orthogonal distance regression to account for the errors in both x and y. This allowed us to find an optical depth corresponding to the definitions of Falstad et al. (2021). Taking Σ_HCN-vib > 1 L_⊙ pc⁻² requires τ_N > 3.5 and L_HCN-vib/L_IR > 10⁻⁸ requires τ_N > 3.3. As the surface brightness is a more robust definition and gives a stricter threshold, we adopt τ_N > 3.5 as a criterion to select deeply obscured nuclei. In addition, we impose the extra condition that the nuclear contribution must be greater than 40% to ensure the value of τ_N is reliable from the fitting as discussed in Sect. 3.5.

This now means we have an additional method for selecting CON candidates in the mid-infrared based on the decomposition technique presented in this work.

4.3. CON spectral shape

In Fig. 8 we show the continua of the nuclear components of CONs in LIRGs and ULIRGs from the GOALS and HERUS samples, respectively. These were selected based on the optical depth as described in the following paragraphs. While all show a deep silicate absorption, the slope of the continuum between 5 and 8 μm varies significantly between objects, with some displaying flat continua, such as ESO 374-IG 032, and others showing steeper spectra, such as IRAS 17578−0400. This may reflect differences in the amount of hot dust present since hotter dust will result in a flatter continuum towards shorter wavelengths. A possible explanation for this might be the presence of hotter dust where a more direct line of sight to the central engine allows the dust to reach higher temperatures (e.g., Efstathiou et al. 2022; Lyu & Rieke 2018). This may suggest an inclination dependence where more face-on sources show this hotter dust; however, it remains difficult to reconcile this picture with the large levels of obscuration required to produce the deep silicate absorption.

Fig. 8. Continua of the nuclear component for LIRGs (black) and ULIRGs (red) selected as CONs from the GOALS and HERUS samples, respectively. The spectrum is normalised by the integrated flux between 5 and 15 μm.

4.4. Spectral mapping sample

In this section we present the results of the fits the spectra from from the spectral mapping sample, where we fitted all the apertures simultaneously (and the shape of the nuclear spectra is the same for all of them, as discussed in Sect. 3.5).

Figure 9 shows the PAH EW ratios of this sample and how these vary with aperture size. From this figure it is clear that the PAH EW ratio values from the inner aperture deviate from the others for a number of objects, for example Arp 236 A and Mrk 938. Before analysing further we investigated the effect of the aperture correction on the spectra extracted from this aperture. The slope of the continuum is strongly dependent on this correction function and thus the PAH EWs will depend on this. Therefore, due to the aperture corrections, the CON selection box that is calibrated from staring mode spectra may not be directly applicable to spectral maps.

Fig. 9. EW ratios of the 12.7/11.3 PAH and the 6.2/11.3 PAH for the spectral mapping sample. For each galaxy, values are coloured by the aperture radius used to extract the spectra; a larger aperture will contain more emission from star formation in the galactic disc. The inner purple box shows the CON selection region as defined for the staring mode spectra, and the outer box shows a larger selection region that accounts for aperture correction effects. The grey lines show the mean EW ratios for the star-forming calibration sample. The two bona fide CONs are shown with a red title.

4.4.1. Aperture correction effects

To quantify if the aperture correction has a significant effect we compared the PAH EW ratios measured from the spectral mapping mode compared with those from the staring mode. We used NGC 5990 to do this as it has a significant nuclear component (AGN fraction of 0.49 from Díaz-Santos et al. 2017). Figure 10 shows the PAH EW ratios for each aperture of the spectral maps with the staring mode added as the black point. The aperture of radius 2.5″ (in purple) has the aperture correction applied and shows a deviation from the interpolated curve when including the staring mode point but excluding this point. This suggests that the aperture correction may be responsible for changing the PAH EW ratio values for the innermost aperture.

Fig. 10. PAH EW ratios for NGC 5990. The coloured points are from spectral mapping observations and are the same as in Fig. 9. The additional black point is the staring mode spectra, which have a slit radius of approximately 1.8″. The purple point has an aperture radius of 2.5″ and has an aperture correction applied as described in Sect. 2. The black line shows a cubic spline interpolation that excludes this point. The inner purple box shows the CON selection region as defined for the staring mode spectra, and the outer box shows a larger selection region that accounts for aperture correction effects.

To account for this uncertainty we extended the CON selection criteria to a larger region using the EW ratios for NGC 5990. We used a cubic interpolation of each EW ratio against aperture radius excluding the 2.5″ aperture and measured the deviation of the interpolated vs. measured PAH EW ratios at 2.5″. To do this we used a radius of 1.8″ for the staring mode spectra as this roughly corresponds to the slit width. Using the measured deviation for the 2.5″ aperture we created a larger box shown in Figs. 9 and 10 where the 12.7/11.3 EW threshold is increased by 0.058 and the 6.2/11.3 EW threshold is increased by 0.074. This creates an ‘error’ region where the EW ratios inferred from the innermost aperture are consistent with a CON.

4.4.2. CONs in the spectral mapping sample

Using the adjusted PAH EW CON selection criteria we find two objects that meet both the 12.7/11.3 and the 6.2/11.3 PAH EW criteria: ESO 320-G030 and ZW 049.057, both of which are known CONs with a HCN-vib surface brightness of Σ_HCN-vib > 1 L_⊙ pc⁻² (Falstad et al. 2021).

There are additional objects that may be CONs, selected by only one of the PAH EW criteria. The strongest candidate of these is Arp 299 A as it has a low enough 12.7/11.3 PAH EW ratio, which is more reliable than the 6.2/11.3. This object also almost meets the optical depth criterion of τ_N > 3.5 with τ_N = 3.3. Other objects that could be CONs are Arp 236 A and Mrk 938, which are both within 1σ of the adjusted selection region. Additionally, NGC 6926 meets both PAH EW ratios with the staring mode but this spectrum suffers from some reduction issues.

For the two known CONs (ESO 302-G030 and ZW049.057) we see a clear trend with how the EW ratios change with aperture size. As more of the galaxy disc is included within the aperture, the EW ratio moves away from the CON selection region and towards the mean values of the star-forming calibration sample. This shows that dilution from the host galaxy can cause objects to be missed by the PAH EW method; however, higher spatial resolution will overcome this issue. This last point is important considering the factor of ∼10 increase in the spatial resolution of the Mid Infrared Instrument (MIRI) on the JWST compared to the Spitzer spectral maps.

A few other objects show trends with aperture size such as the other potential CON candidates Mrk 938 and NGC 6926 but also NGC 7552 and UGC 03410. Another interesting trend with aperture size can be seen in Arp 299 A, where the 12.7/11.3 PAH EW ratio increases as the aperture decreases. The high optical depth measured for this object means that the continuum will cause the PAH EW ratio to move in the opposite direction therefore changes in the 12.7/11.3 PAH flux ratio must drive this change and so the PAH properties may be different in the nuclear region.

The flux ratio of PAHs are known to be sensitive to changes in the properties of PAHs such as the size, charge and hardness of the incident radiation field (see Li 2020, for a review). There is evidence that the PAH molecules in the vicinity of AGN are larger and more neutral (García-Bernete et al. 2022b,c) and may even be excited by AGN on very small scales (Jensen et al. 2017). It, therefore, is plausible that Arp1299 A represents an AGN that used to be completely obscured and has expelled some dust exposing PAH molecules to radiation from the central engine and therefore changing the properties of PAH molecules in the nuclear region. Arp 299 B is a known AGN (Alonso-Herrero et al. 2009) that also shows a strong increase in the 12.7/11.3 PAH EW ratio with the smallest aperture, although this is likely strongly affected by the aperture correction, as discussed in Sect. 4.4.1.

5. Discussion

5.1. How many CONs exist in the local Universe?

The first systematic search for CONs was done by Falstad et al. (2021) using the HCN-vib line; they found CONs in ${38}_{-13}^{+18}\%$ of ULIRGs, ${21}_{-6}^{+12}\%$ of LIRGs, and $0_{-0}^{+9}\%$ of sub-LIRGs. This result is limited by the small sample size, hence the large errors. Using the PAH EW technique, García-Bernete et al. (2022a) found CONs in 30% of ULIRGs but only 7% of LIRGs, with the discrepancy likely due to dilution of the nuclear continuum emission by star formation in the disc of the host galaxy. To build on this work we can now use the optical depth selection criteria to recover those CONs that were missed by the PAH EW technique to achieve a more accurate estimate of the number of CONs in ULIRGs and LIRGs.

For the ULIRG sample we used HERUS (Farrah et al. 2013), as outlined in Sect. 2. For the LIRG sample we used GOALS (Armus et al. 2009). To estimate 1σ uncertainties on the fraction of CONs in each of the sample, we followed Falstad et al. (2021) and used the Cameron (2011) method to construct a beta distribution of the CON fraction, which depends on the sample size.

From the HERUS sample of ULIRGs we identify ${29}_{-7}^{+7}\%$ as CONs from the PAH EW criteria as shown in the bottom left panel of Fig. 6, in accordance with García-Bernete et al. (2022a). From the nuclear optical depth definition of τ_N > 3.5 we select slightly more objects with ${36}_{-7}^{+8}\%$ as CONs.

The EW ratios of the full GOALS sample are shown in the bottom right panel of Fig. 6 including the ULIRGs in that sample. From this study, we exclude objects with only spectral maps available as the required aperture corrections may bias the results. Out of the LIRGs in the GOALS sample the EW ratios select $7.7_{-2.0}^{+2.3}\%$ as CONs, consistent with García-Bernete et al. (2022a). The optical depth definition selects ${17}_{-3}^{+3}\%$ as CONs, which is consistent with the CONquest results. In the case of LIRGs, the difference between the optical depth and PAH EW ratios techniques is larger.

To understand this discrepancy we show in Fig. 11 the measured nuclear fraction, β, for the LIRGs selected as hosting CONs with each technique. The plot shows that the nuclear optical depth method selects additional sources with lower nuclear fractions and thus greater dilution from the host galaxy, which the PAH EW method misses. The discrepancy being larger for LIRGs compared to ULIRGs is a consequence of the fact that LIRGs contain more extended star-forming components and thus are more susceptible to dilution from the host galaxy. This is consistent with what was observed with the spectral mapping data in Sect. 4.4.

Fig. 11. Histogram of the nuclear fraction, β, for LIRGs in the GOALS sample selected as CONs via the PAH EW method (red) and the nuclear optical depth (black). This shows that the nuclear optical depth identifies additional sources that are more heavily diluted by the host galaxy.

For both LIRGs and ULIRGs the fraction of CONs is consistent with the results of the CONquest investigation (Falstad et al. 2021) with overlapping 1σ intervals. A summary comparing the various methods is given in Table 2. The larger sample sizes afforded by using mid-infrared observations results in tighter constraints on these numbers. It is worth noting that the peak of the distribution for the number of CONs in LIRGs is slightly lower than the CONquest results. This may just be statistical error but it could also be due to objects where the nucleus is so diluted by the host galaxy that they are excluded by our cut requiring a nuclear contribution of 40%.

A comprehensive table of all the galaxies analysed in this work can be found in Appendix E, where measured PAH EW ratios and nuclear fraction/optical depths are reported. We also report on the detection of the HCN 14 μm absorption line (Lahuis et al. 2007). Finally, we also include information on the detections of the 23 μm, 28 μm or 33 μm crystalline silicate absorption features from Spoon et al. (2022). These features in absorption indicate heavily obscured nuclei via two methods. Method I indicates sources with the 23 μm feature and 33 μm feature in absorption (i.e., s₂₃ < 0.0 and s₃₃ < 0). Method II indicates sources with s₂₃ < −0.09 and s₂₃ < −0.02. Spoon et al. (2022) also includes a method III for detection of the blue wing of the 33 μm feature where the full feature is cut due to the redshift of the source. No objects analysed in this work are classified by method III.

5.2. Effect of galaxy inclination

As mentioned previously, the PAH EW method can select CON candidates independent of galaxy inclination as extinction from the star-forming disc will affect both the PAH flux ratio and the continuum ratio leaving the EW ratio unchanged. However, this is not necessarily the case for the nuclear optical depth selection. In Fig. 6 there is one object, NGC 3628, with a τ_N > 3.5 but is not selected as a CON by either the PAH EW ratios or the HCN-vib surface brightness. This is a highly inclined galaxy, which suggests that the nuclear optical depth obtained from our spectral decomposition fitting may be sensitive to such galaxies with strong dust lanes obscuring the line of sight, which may bias our estimates of the number of CONs.

To quantify if this is an issue we obtained estimates for the inclination of galaxies in the GOALS sample from the NASA/IPAC Extragalactic Database (NED)² using the ratio of the minor axis, b, to the major axis, a, to estimate the inclination, i where cos i = b/a. For the majority of the sample, the axis ratios were from 2MASS imaging in the K_s band (Jarrett et al. 2000) or r band imaging from SDSS (York et al. 2000). We selected a sub-sample with τ_N > 3.5 and compared to the full sample. This is shown in Fig. 12.

Fig. 12. Inclination of galaxies in the GOALS sample obtained from axis ratios. The full sample is shown in black, with objects selected as obscured nuclei from the optical depth shown in red and those not as CONs in blue.

Visually, the two samples look similar; however, there may be some extra galaxies selected at the highest inclinations (i ≳ 65°). To determine if there is any statistically significant difference between the samples we performed a two-sample Kolmogorov–Smirnov (K–S) test (Karson 1968) from the SCIPY.STATS package. This tests the null-hypothesis that the two samples are drawn from the same underlying distribution, which typically requires a p-value < 0.05 to reject. We performed the K–S test on the sample with τ_N < 3.5 and τ_N > 3.5 to test whether there is any significant difference. We find a p-value of 0.43, which suggests there is no statistically significant difference between the two samples. It is worth noting, however, that the small sample size of the CONs may make the K-S test less likely to find a difference even if one exists and so it is worth inspecting visually. We would also only expect a difference at high inclinations where the dust lanes would obscure the line of sight towards the nucleus and so if we only compare the highest inclinations (i ≳ 65°), we do see a bias where sources may be falsely selected as CON candidates. The case of NGC 3628 is one such example; however, this case also shows the strength of the PAH EW method for selecting CON candidates.

While there may be a bias falsely selecting highly inclined sources as CON candidates, this bias may actually counteract the exclusion of the most diluted CONs by their host galaxy (β < 0.4) and so while individual objects may be misclassified, the total number of CONs likely remains accurate.

5.3. Future prospects

To really test the impact of CONs on galaxy evolution requires observations across cosmic time, particularly at cosmic-noon (z ∼ 1 − 3). The advent of the JWST will be at the forefront of finding and understanding CONs beyond the local Universe (García-Bernete et al. 2022a). In particular the MIRI for JWST (Rieke et al. 2015) will enable the imaging (Wright et al. 2015) and spectroscopy (Wells et al. 2015) of galaxies between 5 and 28 μm.

We found in this work that dilution of the nuclear emission by relatively unobscured star formation is a challenge for identifying CONs with Spitzer. As the JWST provides a nearly 10 times increase in angular resolution, this issue will become apparent at z ∼ 0.5 where the physical scale probed will surpass that of local galaxies with Spitzer.

Future surveys in the mid-infrared, such as PRIMA (Glenn et al. 2021), are also important as they will enable analysis of large numbers of these objects at cosmic noon.

6. Conclusions

In this work we have investigated the physical properties of CONs by decomposing the mid-infrared spectra into nuclear and star-forming components and evaluated how to identify such objects using PAH EW ratios. We have made the code publicly available³. Our main findings are:

• From our spectral decomposition, the optical depth of the nuclear component at 9.8 μm, τ_N, is strongly correlated with the surface brightness of HCN-vib emission at millimetre wavelengths, suggesting the same physics is responsible for both observational signatures. This leads to a CON selection criterion of τ_N > 3.5.

• From τ_N > 3.5 we find that CONs make up ${36}_{-7}^{+8}\%$ of ULIRGs and ${17}_{-3}^{+3}\%$ of LIRGs, consistent with the results of CONquest but with tighter constraints. We find the PAH EW method classifies fewer CONs, likely due to the low spatial resolution of Spitzer IRS data; it detects ${29}_{-7}^{+7}\%$ of ULIRGs and $7.7_{-2.0}^{+2.3}\%$ of LIRGs as CONs.

• We find that the PAH EW method is robust against false positives in highly inclined galaxies with strong dust lanes that can produce high optical depths, whereas using the nuclear optical depth to select CON sources may falsely select some objects that are highly inclined.

• From spectra extracted at different spatial scales, we find the PAH EW ratios move towards the CON selection criterion with smaller apertures where there is less dilution from the disc. This suggests that star formation diluting the nucleus is responsible for the EW ratios underestimating the number of LIRGs that host CONs.

We have confirmed that the mid-infrared can be used to effectively select completely obscured nuclei and investigate properties of the dust. This will allow these objects to be found and studied beyond the local Universe, where HCN-vib emission is simply too faint to be detected, even with ALMA. With its high sensitivity, the JWST will allow us to find and understand these objects at cosmic noon, where the impact of this powerful but hidden phase of galaxy evolution can be uncovered.

Acknowledgments

FRD acknowledges support from STFC through grant ST/W507726/1. DR and IGB acknowledge support from STFC through grant ST/S000488/1. DR also acknowledges support from the University of Oxford John Fell Fund. AAH acknowledges support from grant PGC2018-094671-B-I00 funded by MCIN/AEI/ 10.13039/501100011033 and by ERDF A way of making Europe, and grant PID2021-124665NB-I00. We thank Michalina Maksymowicz-Maclata for providing the list crystalline absorption detections. We also thank the reviewer for the useful feedback.

References

Appendix A: Extinction templates

There are absorption features observed between 5.5 - 8 μm in strongly obscured galaxies, commonly attributed to water ice absorption at ∼6μm (Spoon et al. 2001) and the deformation mode of aliphatic CH molecules at 6.85 and 7.25μm (Dartois et al. 2007) in addition to the deep silicate absorption at ∼9.8μm. A theoretical template remains elusive as the molecular composition responsible for this absorption is very complex; therefore, we follow Spoon et al. (2022) and resort to constructing a template from a heavily obscured source that has very little PAH emission contaminating the region, namely NGC 4418.

Using Spitzer IRS spectra of NGC 4418, we determine the local underlying continuum for through a cubic spline interpolation anchored at 5.5, 7.8, 13.0 and 14.5 μm. This is shown in the left panel of Fig. A.1 as a blue dashed line. To extract an optical depth, we first mask out any PAH emission or lines before taking the natural log of the ratio of the spline continuum to the masked data. This results in the optical depth in the right panel of Fig. A.1.

Fig. A.1. Creation of the extinction templates. Left: IRS spectrum of NGC 4418 (black) with the spline-interpolated underlying continuum as the dashed blue line. The dashed red line shows the spectrum with emission features removed. Right: Optical depth profile using the underlying continuum, overlaid with the inferred template for water ice (cyan), aliphatic CH absorption (red), and silicates (gold).

We fitted the optical depth profile with three components, the ice, CH, and silicates. Spoon et al. (2022) used two Gaussian profiles centred at 6.85 and 7.25μm to represent the CH deformation mode. However, we find the fit is poor and required broader wings to fit the optical depth spectrum. Therefore, we fitted two Drude profiles to the optical depth spectrum instead, which is shown in the right panel of Fig. A.1 (Lorentzian profiles also provided a good fit, though the χ² was slightly higher). This provides the template for the CH absorption, while the ice feature is given by the residuals smoothed to provide a template, shown in cyan in the figure. The silicate profile is also given by smoothing the residuals.

Appendix B: Choice of silicate template

In our analysis we chose to use an empirical template for the 9.8 μm silicate absorption feature, derived from NGC 4418 as described in Appendix A. We chose this source as it is highly obscured with minimal emission features and is a well studied CON with strong HCN-vib emission (e.g. Sakamoto et al. 2010). However, as noted in Sect. 3.2, ground based observations with a smaller beam find a higher peak optical depth than with Spitzer IRS, which suggests the template contains some contribution from the relatively unobscured continuum from circumnuclear star formation. We therefore test another template derived from another highly obscured galaxy, IRAS 08572+3915.

The profile for this galaxy is narrower, which suggests less contamination by any unobscured continuum, consistent with a total absence of emission features. We tested the model using this template on the CONquest sample and compared the resulting decomposition. In Fig. B.1 we compare the measured optical depths using the IRAS 08572+3915 template in blue to the NGC 4418 template in red. As expected the IRAS 08572 template produces higher nuclear optical depths requiring more contribution from the star-forming component to fit the data. Crucially this is a constant offset in optical depth with the same correlation found but shifted to higher values. Therefore, the conclusions of this paper are not strongly affected by the choice of silicate template.

Fig. B.1. Comparison of the measured nuclear optical depths of the CONquest sample with different silicate templates. In red are the measured τN against the surface density of HCN-vib (left) and HCN-vib luminosity to LIR (right) using the NGC 4418 profile. In blue are the results using the template from IRAS 08572+3915.

While this template may be advantageous for future work, we chose not to use it as this object is not a CON by the HCN-vib definition.

Appendix C: Silicate strength vs. HCN-vib

Fig. C.1. Same as Fig. 7 but with the apparent silicate strength as calculated by Eq. 6 (red), which shows no trend. The 12.7/11.3 PAH EW ratio is also shown (blue) against HCN-vib on the right-hand axis of each plot and shows a very weak trend. A comparison with Fig. 7 shows the value in accounting for the star-forming contribution to properly recover properties of the obscured nucleus.

Appendix D: Continuum ratios of the spectral mapping sample

Fig. D.1. Continuum ratios at 12.7 μm to 11.3 μm against the 6.2/11.3 continuum ratio for the spectral mapping sample. For each galaxy, values are coloured by the aperture radius used to extract the spectra; a larger aperture will contain more emission from star formation in the galactic disc.

Appendix E: Table of spectral properties

All Tables

All Figures

Fig. 1. Example of the spectral extraction, in this case for Arp 299 A, at different spatial scales. Left: integrated intensity map of the SL1 cube (7.46–14.29 μm) overlaid with circular apertures of radii 2.5, 4.5, 6.5, and 8.5″, shown as the dashed lines. Right: extracted spectra from the apertures in the left panel with the corresponding colour. The innermost spectrum (purple) has had an aperture correction applied.

In the text

Fig. 2. Star-forming continuum templates derived from the star-forming sample. Left: unobscured continua, CSF(λ), obtained from the fits to the star-forming sample, used as templates for the spectral decomposition. The colour indicates the steepness (8.0 μm–13 μm): the darker the colour, the steeper the continuum. Right: histogram of the optical depths, assuming mixed geometry, for the star-forming sample. The red line shows the half-normal prior used for the spectral decomposition based on these results.

In the text

Fig. 3. Example of our spectral decomposition fit (described in Sect. 3) to simulated data of a galaxy that hosts a CON. The left panel shows a highly diluted nucleus (β = 0.42) with a spectrum dominated by star formation, and the right panel shows a CON-dominated source (β = 0.89). The full model, shown in red, is composed of the various components shown in the legend. The continuum is shown as the solid gold line and is a sum of the nuclear and star-forming (SF) components, which are shown as dashed lines. The flux is in arbitrary units. The equation of the full model is given in Eq. (1).

In the text

Fig. 4. Testing the recovered properties of simulated data using the spectral decomposition fitting method. Left: measured nuclear fraction, β, against the true value. The solid black line shows where the measured equals the true value. Right: measured nuclear optical depth, τN, against the true nuclear fraction. The solid horizontal lines show the true optical depth of τN = 4.5 in red and τN = 3.0 in blue. The vertical dashed black lines show a cutoff of β = 0.4, below which the measured values become unreliable.

In the text

Fig. 5. Comparison between fitting spectra from individual apertures vs. fitting all apertures simultaneously with shared nuclear parameters for ESO 320-G030. Left: nuclear fraction, β, measured from the simultaneous fitting on the y-axis vs. the individual spectra on the x-axis. The solid black line shows where they are equal. Right: measured optical depth of the nucleus for each aperture individually and for all apertures simultaneously. Note the much smaller error when using all the spectra simultaneously.

In the text

Fig. 6. PAH EW diagram (García-Bernete et al. 2022a) of the 12.7/11.3 PAH features against the 6.2/11.3 PAH features obtained from fitting the star-forming sample (top left; Hernán-Caballero et al. 2020), CONquest (top right; Falstad et al. 2021), HERUS (bottom left; Farrah et al. 2013), and GOALS (bottom right; Armus et al. 2009). Points in red have a nuclear optical depth τN > 3.5, and for the CONquest sample the triangle points are those identified as CONs by ΣHCN-vib > 1 L⊙ pc−2. The grey lines show the mean values for the star-forming sample, with the 2σ values for the 12.7/11.3 PAH shown by the dashed lines. The purple shaded area shows the CON selection region.

In the text

Fig. 7. Optical depth of the nuclear component from the spectral decomposition fits to the CONquest sample (Falstad et al. 2021) against the surface brightness of HCN-vib (left) and the strength of HCN-vib to LIR (right). Linear fits are shown in the dashed black lines, with 1σ confidence intervals shown in grey. The upper limits of HCN-vib measurements are shown for objects with τN > 2.5 and β > 0.4. The dashed grey lines show the CON identification threshold from Falstad et al. (2021) on the x-axis and the corresponding τN on the y-axis.

In the text

	Fig. 8. Continua of the nuclear component for LIRGs (black) and ULIRGs (red) selected as CONs from the GOALS and HERUS samples, respectively. The spectrum is normalised by the integrated flux between 5 and 15 μm.
In the text

Fig. 9. EW ratios of the 12.7/11.3 PAH and the 6.2/11.3 PAH for the spectral mapping sample. For each galaxy, values are coloured by the aperture radius used to extract the spectra; a larger aperture will contain more emission from star formation in the galactic disc. The inner purple box shows the CON selection region as defined for the staring mode spectra, and the outer box shows a larger selection region that accounts for aperture correction effects. The grey lines show the mean EW ratios for the star-forming calibration sample. The two bona fide CONs are shown with a red title.

In the text

Fig. 10. PAH EW ratios for NGC 5990. The coloured points are from spectral mapping observations and are the same as in Fig. 9. The additional black point is the staring mode spectra, which have a slit radius of approximately 1.8″. The purple point has an aperture radius of 2.5″ and has an aperture correction applied as described in Sect. 2. The black line shows a cubic spline interpolation that excludes this point. The inner purple box shows the CON selection region as defined for the staring mode spectra, and the outer box shows a larger selection region that accounts for aperture correction effects.

In the text

	Fig. 11. Histogram of the nuclear fraction, β, for LIRGs in the GOALS sample selected as CONs via the PAH EW method (red) and the nuclear optical depth (black). This shows that the nuclear optical depth identifies additional sources that are more heavily diluted by the host galaxy.
In the text

	Fig. 12. Inclination of galaxies in the GOALS sample obtained from axis ratios. The full sample is shown in black, with objects selected as obscured nuclei from the optical depth shown in red and those not as CONs in blue.
In the text

Fig. A.1. Creation of the extinction templates. Left: IRS spectrum of NGC 4418 (black) with the spline-interpolated underlying continuum as the dashed blue line. The dashed red line shows the spectrum with emission features removed. Right: Optical depth profile using the underlying continuum, overlaid with the inferred template for water ice (cyan), aliphatic CH absorption (red), and silicates (gold).

In the text

	Fig. B.1. Comparison of the measured nuclear optical depths of the CONquest sample with different silicate templates. In red are the measured τN against the surface density of HCN-vib (left) and HCN-vib luminosity to LIR (right) using the NGC 4418 profile. In blue are the results using the template from IRAS 08572+3915.
In the text

	Fig. C.1. Same as Fig. 7 but with the apparent silicate strength as calculated by Eq. 6 (red), which shows no trend. The 12.7/11.3 PAH EW ratio is also shown (blue) against HCN-vib on the right-hand axis of each plot and shows a very weak trend. A comparison with Fig. 7 shows the value in accounting for the star-forming contribution to properly recover properties of the obscured nucleus.
In the text

	Fig. D.1. Continuum ratios at 12.7 μm to 11.3 μm against the 6.2/11.3 continuum ratio for the spectral mapping sample. For each galaxy, values are coloured by the aperture radius used to extract the spectra; a larger aperture will contain more emission from star formation in the galactic disc.
In the text