Dense Molecular Gas Properties of the Central Kpc of Nearby Ultraluminous Infrared Galaxies Constrained by ALMA Three Transition-line Observations

We report the results of ALMA 1-2 kpc-resolution, three rotational transition line (J=2-1, J=3-2, and J=4-3) observations of multiple dense molecular gas tracers (HCN, HCO$^{+}$, and HNC) for ten nearby (ultra)luminous infrared galaxies ([U]LIRGs). Following the matching of beam sizes to 1-2 kpc for each (U)LIRG, the high-J to low-J transition-line flux ratios of each molecule and the emission line flux ratios of different molecules at each J transition are derived. We conduct RADEX non-LTE model calculations and find that, under a wide range of gas density and kinetic temperature, the observed HCN-to-HCO$^{+}$ flux ratios in the overall (U)LIRGs are naturally reproduced with enhanced HCN abundance compared to HCO$^{+}$. Thereafter, molecular gas properties are constrained primarily through the use of HCN and HCO$^{+}$ data and the adoption of fiducial values for the HCO$^{+}$ column density and HCN-to-HCO$^{+}$ abundance ratio. We quantitatively confirm the following: (1) Molecular gas at the (U)LIRGs' nuclei is dense ($\gtrsim$10$^{3-4}$ cm$^{-3}$) and warm ($\gtrsim$100 K). (2) Molecular gas density and temperature in nine ULIRGs' nuclei are significantly higher than that of one LIRG's nucleus. (3) Molecular gas in starburst-dominated sources tends to be less dense and cooler than ULIRGs with luminous AGN signatures. For six selected sources, we also apply a Bayesian approach by freeing all parameters and support the above main results. Our ALMA 1-2 kpc resolution, multiple transition-line data of multiple molecules are a very powerful tool for scrutinizing the properties of molecular gas concentrated around luminous energy sources in nearby (U)LIRGs' nuclei.


INTRODUCTION
Galaxies whose infrared (8-1000 µm) luminosities are L IR = 10 12−13 L ⊙ and L IR = 10 11−12 L ⊙ are referred to as ultraluminous infrared galaxies (ULIRGs) and luminous infrared galaxies (LIRGs), respectively (Sanders & Mirabel 1996). Most (U)LIRGs exhibit strong infrared dust thermal radiation which is much brighter than the UV-optical emission, and are found in major mergers of gas-rich galaxies in the nearby universe at z < 0.3 (e.g., Sanders & Mirabel 1996). Starburst activity and/or luminous active galactic nucleus (AGN) activity powered by a mass-accreting supermassive black hole (SMBH) at nuclear regions hidden behind the curtain of obscuring gas and dust, are considered the primary energy sources for the observed large infrared luminosities of nearby (U)LIRGs. However, estimating starburst and AGN contribution to the large luminosities is challenging because of very large extinction in nearby (U)LIRGs' nuclei.
Numerical simulations have predicted that molecular gas is transferred to nuclear regions through merger-induced processes (e.g., Juneau et al. 2009;Moreno et al. 2019) and a large amount of dense ( 10 3−4 cm −3 ) molecular gas is observed in nearby merging (U)LIRGs' nuclei (e.g., Gao & Solomon 2004;Juneau et al. 2009;Imanishi et al. 2019;Ueda et al. 2021;Imanishi et al. 2022). Such dense molecular gas can (1) form stars and become an important fuel toward a central SMBH (i.e., enhancing AGN activity) and (2) be affected by the hidden energy sources of (U)LIRGs, both radiatively and mechanically. Scrutinizing the physical and chemical properties of dense molecular gas at the (U)LIRGs' nuclei is vital to better understand the nature of these regions. Rotational J-transition emission lines of dense molecular gas tracers are found in the (sub)millimeter wavelength range at 0.3-3.5 mm, where extinction effects are very small. (Sub)millimeter molecular rotational J-transition lines can thus be a very useful tool for the scrutinization.
Spectral line energy distribution (SLED) of CO at (sub)millimeter and far-infrared (30-300 µm), based on large-beam-sized observations with the Herschel Space Observatory and/or ground-based single-dish telescopes ( 5 ′′ or 5 kpc at z 0.05), has often been used to identify the presence of luminous AGN activity, in addition to starburst activity, in nearby (U)LIRGs. This is because high-J (J 4-5) CO emission can be stronger in an AGN than in a starburst when normalized to low-J (J = 1-2) CO luminosity (e.g., Hailey-Dunsheath et al. 2012;Pereira-Santaella et al. 2014;Mashian et al. 2015;Rosenberg et al. 2015;Lu et al. 2017). This can be attributed to the fact that dense and warm ( 100 K) molecular gas in the vicinity of a luminous AGN can excite CO to high-J levels more efficiently than in normal star-forming regions (e.g., Meijerink et al. 2007;Spaans & Meijerink 2008). SLED of molecules with higher dipole moments (e.g., HCN and HCO + ) has also been used to probe the physical conditions of dense molecular gas in nearby (U)LIRGs. However, this has mostly been conducted with large-beam-sized ( 5 ′′ ) observations using single dish (sub)millimeter telescopes (e.g., Krips et al. 2008;Greve et al. 2009;Papadopoulos et al. 2014;Israel 2022). Dense molecular gas properties of energetically dominant compact ( 1-2 kpc) nuclear regions of nearby ULIRGs (e.g., Soifer et al. 2000;Diaz-Santos et al. 2010;Imanishi et al. 2011;Pereira-Santaella et al. 2021) may not be best probed using these low-angular-resolution data, because emission from spatially extended ( a few kpc) host galaxies can be a severe contaminant. Probing only nuclear dense molecular gas through higher-spatial-resolution ( 1-2 kpc) observations is highly desirable.
In this paper, we combine available J=2-1, J=3-2, and J=4-3 data 1 of HCN, HCO + , and HNC previously obtained from our ALMA observations (Imanishi & Nakanishi 2013a,b, 2014Imanishi et al. 2016aImanishi et al. ,b, 2018Imanishi et al. , 2021Imanishi et al. , 2022. After matching the beam sizes of all the J-transition line data of all molecules to the same value for each object (1-2 kpc), we perform the aforementioned investigations. Throughout this study, we adopt the cosmological parameters H 0 = 71 km s −1 Mpc −1 , Ω M = 0.27, and Ω Λ = 0.73. In addition, molecular line flux ratios are calculated in Jy km s −1 . Density and temperature mean H 2 volume number density (n H 2 ) and kinetic temperature (T kin ), respectively, unless otherwise stated.

TARGETS
The targets of this study are nine ULIRGs and one LIRG, previously observed at J=2-1, J=3-2, and J=4-3 of HCN, HCO + , and HNC, with 1-2 kpc resolution, through our ALMA programs (Imanishi & Nakanishi 2013a,b, 2014Imanishi et al. 2016aImanishi et al. ,b, 2018Imanishi et al. , 2021Imanishi et al. , 2022. Table 1 summarizes the basic properties of the (U)LIRGs. We selected nearby (U)LIRGs with different levels of AGN's energetic contributions to bolometric luminosity, estimated using optical, infrared, hard X-ray (>10 keV), and (sub)millimeter spectroscopic energy diagnostic methods. Two sources (NGC 1614 and IRAS 13509+0442) and a secondary fainter nucleus of IRAS 12112+0305 (SW) are regarded as starburst-dominated, without discernible luminous AGN signatures, whereas the remaining ULIRGs' nuclei are diagnosed as containing luminous AGNs. More detailed explanations of the individual (U)LIRGs can be found in Imanishi et al. (2016b). Although the observed ten (U)LIRGs are neither statistically complete nor unbiased, their molecular line studies will provide important information of (1) the general properties of nuclear dense molecular gas in nearby (U)LIRGs at 1-2 kpc resolution and (2) possible difference of molecular gas properties between (U)LIRGs' nuclei with and without luminous AGN signatures.
(1) The J=4-3 emission is generally more spatially compact than the lower J-transition lines owing to the higher critical density in the former (Shirley 2015).
(2) Changing originally small beam-sized data to unnecessarily large beam sizes increases the noise in units of mJy beam −1 , and thus the scatter of the spectral data becomes large. For J=3-2 data of certain molecular lines for IRAS 12127−1412, IRAS 13509+0442, and the Superantennae, and J=4-3 data for certain lines for IRAS 12112+0305 and IRAS 22491−1808, the major axes of the synthesized beams are slightly larger than the adopted values (Table 2); however, it is also assumed that the bulk of the observed J=3-2 and J=4-3 fluxes are emitted within the adopted beams. We obtained 1-2 kpc beam spectra by changing the beam size of continuum-subtracted data cube from the original values shown in the previous publications (Imanishi & Nakanishi 2013a,b, 2014Imanishi et al. 2016aImanishi et al. ,b, 2018Imanishi et al. , 2021Imanishi et al. , 2022, to the adopted ones (Table 2). However, we noticed that the HCN, HCO + , and HNC J=4-3 emission lines of IRAS 22491−1808 obtained in ALMA Cycle 0 were substantially narrower than J=3-2 and J=2-1 data obtained in later ALMA Cycles. This is because (1) the molecular-line-derived redshift of z = 0.0776 (Imanishi & Nakanishi 2014) was found to be significantly larger than the optically derived redshift of z = 0.076 (Kim & Sanders 1998), which had been assumed at the time of the Cycle 0 observation setup, and (2) low-frequency (high-velocity) side of the emission tail was located at the edge of the ALMA spectral windows, which resulted in continuum oversubtraction owing to the improper inclusion of data points with significant emission components. Subsequently, we reanalyzed the HCN, HCO + , and HNC J=4-3 Cycle 0 data of IRAS 22491−1808 by redefining a constant continuum flux level. We found that the revised profiles of the J=4-3 emission lines were more consistent with those of the J=2-1 and J=3-2 data, and thus adopt the revised results.

RESULTS
Newly obtained spectra of the two ULIRGs, IRAS 06035−7102 and IRAS 08572+3915, which were never presented in our previous publications (Imanishi & Nakanishi 2013a,b, 2014Imanishi et al. 2016aImanishi et al. ,b, 2018Imanishi et al. , 2021Imanishi et al. , 2022, are shown in Appendix A. Figure 1 displays the velocity profiles for J=2-1, J=3-2, and J=4-3 emission lines of HCN, HCO + , and HNC for all the observed (U)LIRGs. A Data were obtained in Cycle 0 with a beam size larger than the adopted one. It is assumed that the bulk of the observed J=4-3 flux originates from the nuclear region within the adopted beam (see §3).
B The major axis is slightly larger than that of the adopted beam. We modify only the minor axis to the adopted value. We assume that the bulk of the observed J=3-2 and J=4-3 flux originates from the nuclear region within the adopted beam (see §3).  Gaussian fits are applied to estimate the emission line fluxes within the adopted beam sizes. Single Gaussian fits are adopted, except for the HCN and HCO + lines of IRAS 12112+0305 NE, which clearly exhibit a double-peaked profile with a strong central dip (Figure 1j,k) wherein the results of a double Gaussian fit are adopted. For certain emission lines from other sources, certain weak signatures of double-peaked profiles with modest amounts of central dips are observed; however, we decide to adopt single Gaussian fitting results after confirming that the flux estimates between the single and double Gaussian fits agree within ∼10% for significantly detected emission lines ( 3σ). This is because we aim to estimate the fluxes of multiple J-transition emission lines in a consistent manner, using the simplest profile possible. The final Gaussian fitting results are summarized in Appendix B.
The flux density of the continuum emission, simultaneously taken during individual molecular line observations, is also estimated using the adopted beam sizes for individual (U)LIRGs. Flux density and spectral energy distribution (SED) of the continuum emission are summarized in Appendix C. Table 3 summarizes the adopted Gaussian-fit velocity-integrated emission line fluxes of HCN, HCO + , and HNC at J=2-1, J=3-2, and J=4-3, together with the adopted line width which will be used for model calculations in §5. Table 4 and Figure 2a show the derived HCN-to-HCO + flux ratios at J=2-1, J=3-2, and J=4-3. High-J to low-J emission line flux ratios for HCN, HCO + , and HNC are summarized in Table 5, and those for the HCN and HCO + are plotted in Figure 2b.
The high-J to low-J (J=2-1) flux ratios of the LIRG NGC 1614 are smaller than those of the remaining ULIRGs particularly for HCN (by a factor of 2-5 except IRAS 13509+0442) and for HCO + as well (by a factor of 1-3) ( Figure 2b and Table 5), suggesting that gas density and temperature of NGC 1614 are much lower and thus collisional excitation to high-J levels (J=3 and 4) is insufficient. Furthermore, only NGC 1614 shows (1) a clear systematic decreasing trend in the HCN-to-HCO + flux ratio from low-J to high-J transition in Figure 2a and (2) significantly smaller high-J to low-J flux ratios for HCN than for HCO + in Figure 2b, while the remaining ten ULIRGs' nuclei do not except IRAS 13509+0442 which shows a weaker, but similar decreasing trend in the former. These    Table 9, column 6). This line width will be used for our RADEX non-LTE model calculations.     results can also be explained by the lower gas density and temperature in NGC 1614 because the critical density of HCN is a factor of ∼5 higher than that of HCO + (Shirley 2015); thus, collisional excitation to a high J-level is less efficient for HCN than for HCO + in low-density and low-temperature molecular gas. We will employ non-LTE model calculations to quantitatively constrain the nuclear dense molecular gas properties of the LIRG NGC 1614 and the remaining ULIRGs in more detail in the next section ( §5).

General Descriptions of Our RADEX Non-LTE Modeling
There are now three rotational transition line data (J=2-1, J=3-2, and J=4-3) of HCN, HCO + , and HNC, whose fluxes are measured with the same beam sizes (1-2 kpc) for individual (U)LIRGs. To investigate nuclear (1-2 kpc) molecular gas properties in detail, the observed emission line flux ratios are compared with the calculated flux ratios using RADEX (van der Tak et al. 2007). RADEX is a non-LTE radiative transfer code that solves the statistical equilibrium problems in a one-zone medium based on the escape probability approximations. Molecular emission line flux ratios in the RADEX code are primarily determined by the following three parameters: (i) H 2 volume number density (n H2 ). (ii) H 2 kinetic temperature (T kin ). (iii) Molecular column density divided by line width (N mol /∆v), which reflects the line opacity, where ∆v can be derived from the observed molecular line width (Table 3, column 11). We constrain these three parameters from the three J-transition line data.
In all RADEX non-LTE calculations, molecular gas is assumed to be static and spherically symmetric homogeneous medium. The background temperature is set as the cosmic microwave background (2.73 K). The collision partner for HCN, HCO + , and HNC is regarded as H 2 only. The emission line flux ratios are calculated after converting the fluxes from Kelvin (brightness temperature) to Jy units, as adopted in Tables 4-5. For implementation, pyradex 2 , a Python wrapper for RADEX, is used.
To perform the RADEX non-LTE calculations, we use only the HCN and HCO + data, but exclude the HNC data for the following two reasons. First, in previous observations of the Galactic sources, the HNC abundance is found to be very low in warm molecular gas in the vicinity of a luminous energy source (e.g., Schilke et al. 1992;Hirota et al. 1998;Graninger et al. 2014;Bublitz et al. 2019). In fact, in the nearby (z = 0.0038, luminosity distance ∼ 14 Mpc), well-studied AGN NGC 1068, the HNC emission is found to be extremely weak, compared with HCN and HCO + , in close proximity to the central luminous AGN (Imanishi et al. 2020). Therefore, HNC is not considered to be as good as HCN and HCO + to investigate the physical properties of warm molecular gas near luminous energy sources (AGN and/or starburst) at the center of (U)LIRGs' nuclei. Second, the HNC data were not taken simultaneously with the HCN and HCO + data at the same J-transition (Imanishi & Nakanishi 2013a,b, 2014Imanishi et al. 2016aImanishi et al. ,b, 2018Imanishi et al. , 2021Imanishi et al. , 2022. Thus, possible absolute flux calibration uncertainty of three more ALMA observations (HNC J=2-1, J=3-2, and J=4-3) must be considered as additional free parameters for our RADEX non-LTE calculations (see §5.3 and 5.5), which makes it extremely difficult to find the best fit values without systematic uncertainty under the limited number of observational constraints.

HCN-to-HCO + Emission Line Flux Ratios
In the RADEX model calculations for the (U)LIRGs' nuclei, we set density and temperature in the range 10 3 -10 8 cm −3 and 10-1000 K, respectively. We also adopt the HCO + column density with N HCO+ = 1 × 10 16 cm −2 as the fiducial value, by assuming that the observed (U)LIRGs are modestly Compton-thick (N H ∼ a few × 10 24 cm −2 ) and the HCO + -to-H 2 abundance ratio is ∼10 −8 (e.g., Martin et al. 2006;Saito et al. 2018). Figure 3 shows a three-dimensional plot of the HCN-to-HCO + emission line flux ratio at J=2-1, J=3-2, and J=4-3. The overall distribution of the observed HCN-to-HCO + flux ratios for the (U)LIRGs' nuclei can be naturally reproduced with the HCN-to-HCO + abundance ratio of [HCN]/[HCO + ] = 3 rather than 1. In Appendix D, we also make the same plot with [HCN]/[HCO + ] = 7 and 1, and confirm that the enhanced HCN-to-HCO + abundance ratio can better fit the observed flux ratios for the majority of the (U)LIRGs' nuclei. Two sources with small HCN-to-HCO + flux ratios, located at the bottom-left part of Figure 3 (NGC 1614 and IRAS 12112+0305 SW) may be better reproduced with [HCN]/[HCO + ] = 1 (i.e., comparable HCN and HCO + abundances). Both of these sources are classified as starburst-dominated ( §2). However, the remaining ULIRGs showing luminous AGN signatures (except IRAS 13509+0442) are better reproduced with enhanced HCN abundance relative to HCO + . In fact, the enhanced HCN-to-HCO + abundance ratio (1) has been reported in dense molecular gas around luminous AGNs (e.g., Krips et al. 2008;Aladro et al. 2015;Saito et al. 2018;Nakajima et al. 2018;Takano et al. 2019;Kameno et al. 2020;Imanishi et al. 2020;Butterworth et al. 2022), and (2) has been argued to be necessary to reproduce the HCN-to-HCO + flux ratios with larger than unity, as observed in many (U)LIRGs (e.g., Yamada et al. 2007;Izumi et al. 2016a). We adopt this enhanced HCN abundance of [HCN]/[HCO + ] = 3 as the fiducial value for our model calculations in the next subsection. As we will show in the calculations for the starburst-dominated LIRG NGC 1614 (Figure 4), adopting [HCN]/[HCO + ] = 1 provides comparable gas density and temperature 3 .
We also perform calculations using the HCO + column density of N HCO+ = 1 × 10 15 cm −2 (an order of magnitude smaller than that adopted in Figure 3) and a line width of ∆v = 800 km s −1 (a factor of >2 broader). However, the overall trend of the RADEX-calculated HCN-to-HCO + flux ratios does not change significantly (Appendix D).

High-J to Low-J Emission Line Flux Ratios
To quantitatively estimate the density and temperature of molecular gas at (U)LIRGs' nuclei, a total of four observed line flux ratios, J=3-2 to J=2-1 and J=4-3 to J=2-1 of HCN and HCO + (Table 5) are fitted using the RADEX model. Least-squares fitting of log n H 2 (density) and log T kin (temperature) is performed using conventional Levenberg-Marquardt algorithm through Python package lmfit (Newville et al. 2021). The confidence intervals for the parameters are examined by grid computing ∆χ 2 ≡ χ 2 − χ 2 best with log n H 2 ranging from 2-6 and log T kin from 1-3. Figure 4a shows the derivation of molecular gas density and temperature that can best reproduce the observed J=3-2 to J=2-1 and J=4-3 to J=2-1 flux ratios of HCN and HCO + for the starburst-dominated In Figure 3, the observed HCN-to-HCO + flux ratio for the two sources located at the bottom left part, including NGC 1614, can also be explained by a non-enhanced HCN-to-HCO + abundance ratio ([HCN]/[HCO + ] = 1). Furthermore, it is possible that the LIRG NGC 1614 has lower column density than those of ULIRGs. Thus, the same least-squares fit is conducted by adopting (1) (Figure 4c). Only a limited change in the derived H 2 molecular gas density (n H 2 ) and temperature (T kin ) is observed.
In contrast to the HCN-to-HCO + flux ratios for the same J transition where both HCN and HCO + data were obtained simultaneously, the high-J to low-J flux ratios of HCN and HCO + are affected by the possible absolute flux calibration uncertainty of individual ALMA observations. This is because data for J=4-3, J=3-2, and J=2-1 were obtained at different times. Therefore, this possible systematic uncertainty must be considered. The absolute flux calibration uncertainty is as high as The green open squares represent the best-fit RADEX model. Its gas density and temperature values are indicated by the plus sign in the right panel, whereas other fixed parameters (HCO + column density, HCN-to-HCO + abundance ratio, and line width) are listed at the top of the figure below the object name. The light blue filled circles represent the flux ratios measured from our observations. In (d), the dark blue half-filled circles represent the flux ratios by allowing scale adjustment within the absolute flux calibration uncertainty of individual ALMA observations. Below the legend, the scaling factors to be multiplied for the line fluxes are listed. If a factor is colored red, it implies that the observation was made in Cycle 0 (double systematic error). For example, in (d), the J=3-2 to J=2-1 flux ratios are multiplied by a factor of 1.10/0.95 = 1.16. In (a)-(d), the reduced χ 2 value for the best-fit model is listed at the bottom middle part. The HNC line flux ratios in the shaded part of the left panel are shown only for reference and are not used in the fitting; the HNC abundance is assumed to be equal to that of HCO + . (Right panel): Confidence ranges for gas density and temperature. The plus sign indicates the position of the best-fit model. The contours represent 68, 90, and 99% confidence levels for the two parameters of interest (∆χ 2 = 2.28, 4.61, 9.21). 5% for J21a 4 (band 5) and 10% for J32a and J43a (bands 6 and 7), respectively, according to the ALMA Proposer's Guide. However, for Cycle 0 observations, we conservatively choose to adopt twice the above values, because of possibly large systematic uncertainty in ALMA very early phase. We allow scaling of the absolute fluxes within the above range. If individual emission line fluxes are scaled, then their flux ratios will change accordingly.
We adopt two steps. First, in addition to gas density (n H 2 ) and temperature (T kin ), the scaling factors for J21a, J32a, and J43a are treated as free parameters that move within the systematic uncertainty. Least-squares fitting with box constraints is performed with L-BFGS-B algorithm (Byrd et al. 1995). Subsequently, the scaling factors are fixed at the determined values. Then, gas density and temperature are derived using the Levenberg-Marquardt algorithm in the same manner as before. Figure 4d shows the results obtained by allowing the flux scale adjustment. The best fit values based on the four models in Figure 4 are summarized in Table 6. The derived gas density and temperature largely agree in the range n H 2 = 10 3.4−4.1 cm −3 and T kin = 10 2.0−2.3 K for all the methods in Figure  4. Figure 4 shows that the derived gas density and temperature are not strongly dependent on the choice of the HCN-to-HCO + abundance ratio ([HCN]/[HCO + ]) and HCO + column density (N HCO + ). Thus, we perform the same fit for the remaining ULIRGs by adopting the fiducial values of [HCN]/[HCO + ] = 3 and N HCO + = 1 × 10 16 cm −2 . We adopt the line widths (∆v) listed in Table 3 (column 11) which are not the same among individual ULIRGs. As in the case of NGC 1614, we (1) use the observed high-J to low-J flux ratios as they are, and (2) allow scale adjustment for individual J=4-3, J=3-2, and J=2-1 data to consider the possible absolute flux calibration uncertainty of individual ALMA observations. The derived gas density and temperature largely agree between the first and second fitting results, although the second one usually provides a smaller reduced χ 2 value than the former. We basically adopt the first fitting result, but do the second one only if the first one cannot determine the most likely value or gives a very large reduced χ 2 value. Figure 5 shows our adopted result for each ULIRG's nucleus. Table 6 summarizes the best fit values for the ULIRGs.  Table 3 (Column 11). Symbols are the same as those in Figure 4. We adopt the second fitting result (that is, flux-scaling adjustment allowed) for IRAS 12112+0305 SW, IRAS 12127−1412, IRAS 15250+3609, the Superantennae, IRAS 20551−4250, and IRAS 22491−1808.

Overall Trend of Molecular Gas Properties at (U)LIRGs' Nuclei
In Figure 5, the uncertainty in the derived gas density and temperature is large for IRAS 12112+0305 SW and IRAS 12127−1412 simply because molecular line emission is fainter than other (U)LIRGs' nuclei ( Figure 1). Excluding these two ULIRGs' nuclei, we see the following two trends in Figures 4 and 5; (1) The only one LIRG NGC 1614 contains nuclear molecular gas with modestly high density (10 3.4−4.1 cm −3 ) and high temperature (100-200 K or 10 2.0−2.3 K).
(2) Nuclear molecular gas probed with HCN and HCO + at J=2-1 to J=4-3 in all ULIRGs but IRAS 13509+0442 is significantly denser (10 4.5−5.5 cm −3 ) and warmer ( 300 K or 10 2.5 K) than that in the LIRG NGC 1614. These two results constitute our first main arguments. The gas density and temperature at these ULIRGs' nuclei are also systematically higher than those of other LIRGs in the literature, derived from larger beam-sized CO multiple J-transition (J=1-0 to J=6-5) line data and RADEX modeling (mostly 10 2.0−4.0 cm −3 and 10-90 K) (e.g., Papadopoulos et al. 2012b). IRAS 13509+0442 has relatively low density (∼10 4.3 cm −3 ) and temperature (∼130 K or ∼10 2.1 K) nuclear molecular gas compared to other ULIRGs. Although there are certain difference among the observed sources, the derived gas density and temperature of all (U)LIRGs' nuclei are substantially higher than those of quiescently star-forming, less infrared luminous galaxies at kpc scale, suggested from the observed small high-J to low-J emission line flux ratios of HCN and HCO + (Garcia-Rodriguez et al. 2023) or derived from CO multiple J-transition (J=1-0 to J=3-2) line data and RADEX modeling with fixed N CO /∆v ratio ( 10 4.0 cm −3 and 50 K) (e.g., Leroy et al. 2022).
NGC 1614 and IRAS 13509−0442, together with IRAS 12112+0305 SW, are classified as starburstdominated with no luminous AGN signature at any wavelength, whereas the remaining ULIRGs are regarded as containing luminous AGNs based on previous observations ( §2). In other words, there is a trend that nuclear molecular gas is denser and warmer in ULIRGs with luminous AGN signatures than in (U)LIRGs without. This is our second main argument. Figure 6 shows the contours of the derived gas density and temperature for (U)LIRGs with small uncertainty (excluding IRAS 12112+0305 SW and IRAS 12127−1412). It has previously been suggested that AGN activity is enhanced in galaxy nuclei with a larger amount of and/or higher fraction of dense molecular gas (Juneau et al. 2009;Izumi et al. 2016b). This appears reasonable because the mass accretion rate onto a SMBH is enhanced by a larger amount of dense molecular gas near the SMBH. Owing to the high SMBH mass accretion, luminous AGN activity can emerge and make the surrounding dense molecular gas warm. The second argument can be explained naturally by the combination of these two effects.
For the three ULIRGs, IRAS 08572+3915, IRAS 12112+0305, and IRAS 22491−1808, molecular gas temperature was estimated to be 60-150 K, based on large beam-sized CO multiple J-transition line observations and RADEX modeling (Papadopoulos et al. 2012a). Our ALMA HCN and HCO + line study provides significantly higher temperature for nuclear molecular gas in all the three ULIRGs (see Figure 5 and Table 6). This suggests that our ALMA 1-2 kpc resolution dense molecular line observations selectively probe warm nuclear molecular gas components around luminous energy sources with minimum contamination from spatially extended ( a few kpc) cool components. This advantage was also demonstrated in ALMA high-spatial-resolution HCN and HCO + observations of the nearby well-studied, optically identified AGN NGC 1068, to probe molecular gas only in the vicinity of a luminous energy source (e.g., Viti et al. 2014). Baba et al. (2018) revealed the presence of high temperature ( 200 K) molecular gas at the innermost part of obscuring material around Figure 6. Summary of contours of the RADEX-derived molecular gas density in cm −3 (abscissa) and temperature in K (ordinate) for the observed (U)LIRGs, excluding the two faint molecular emission line ULIRGs (IRAS 12112+0305 SW and IRAS 12127−1412). Color corresponds to infrared-derived AGN luminosity in units of solar luminosity (L ⊙ ), as listed in Table 1 (columns 9 and 11). NGC 1614 is shown as a dotted line because the upper limit of AGN luminosity is looser than that of other ULIRGs. For all the ULIRGs, the AGN luminosity is derived by Nardini et al. (2010) in a consistent manner, whereas that of NGC 1614 is estimated in a different method by a different group (Pereira-Santaella et al. 2015).
luminous AGNs in nearby ULIRGs, based on the observed properties of ∼4.6 µm CO ro-vibrational absorption lines detected in the infrared 4-5 µm spectra. It is likely that HCN and HCO + rotational emission at J=2-1 to J=4-3 detected in our ALMA high-angular-resolution observations largely originates from such innermost high temperature molecular gas at nearby ULIRGs' nuclei.

Bayesian Analysis of Molecular Emission Line Flux Ratios
For selected (U)LIRGs' nuclei exhibiting bright molecular emission lines with high detection significance, molecular gas parameters are estimated by simultaneously fitting HCN-to-HCO + and high-J to low-J emission line flux ratios using RADEX calculations, by making all parameters free, based on a Bayesian approach. The aim is to verify whether the Bayesian approach provides consistent results regarding gas density and temperature with those derived from the previous Levenberg-Marquardt method, where the HCO + column density and HCN-to-HCO + abundance ratios were fixed to fiducial values ( §5.3). There are five independent emission line flux ratios (HCN-to-HCO + flux ratios at J=2-1, J=3-2, and J=4-3, and HCO + J=3-2 to J=2-1 and J=4-3 to J=2-1 flux ratios), which are lower than the total number of parameters, including the absolute flux calibration scaling factors. However, Bayesian estimation facilitates the sampling of the posterior probability distribution of parameters, considering indeterminacy of the solution.
Using the five emission line flux ratios, we explore the parameter space using a Markov Chain Monte Carlo (MCMC) sampler; emcee package (Foreman-Mackey et al. 2013) is used in lmfit. Flat priors with upper and lower bounds presented in Table 7 are adopted. The first guess is made by running the L-BFGS-B solver in §5.3, with the HCO + column density and HCN-to-HCO + abundance ratio unfixed, and then 100 walkers are initialized to be distributed around that guess. The chain is run for 100τ steps, where τ is the maximum number of integrated autocorrelation time for each parameter 5 . The first 5τ steps are discarded as a "burn-in" phase. Subsequently, the chain is thinned out in 0.5τ steps to create independent samples. Hence, the effective number of sampling of the posterior probability distribution is 100 × (100 − 5)/0.5 = 19,000.
For our MCMC analysis, IRAS 12112−1412 SW and IRAS 12127−1412 are excluded, because of their faint molecular emission lines. Examples of the MCMC analysis results for six selected (U)LIRGs' nuclei are shown in Figure 7 and the best fit values are summarized in Table 6. Results for IRAS 15250+3609, the Superantennae, and IRAS 22491−1808 are not shown because we are unable to provide model-calculated results which reproduce the observed results very well with small systematic discrepancy. This independent MCMC analysis supports the previous arguments that (1) molecular gas at the (U)LIRGs' nuclei is dense ( 10 3−4 cm −3 ) and warm ( 100 K) compared to quiescently star-forming normal galaxies with less infrared luminosity and (2) the density of nuclear molecular gas is higher in ULIRGs ( 10 4.3 cm −3 ) than in the LIRG NGC 1614 (∼10 3.3 cm −3 ).

SUMMARY
We presented our ALMA 1-2 kpc-resolution observational results of dense molecular gas tracers (HCN, HCO + , and HNC) at multiple rotational transition lines (J=2-1, 3-2, and 4-3) for one LIRG NGC 1614 (L IR ∼ 10 11.7 L ⊙ ) and nine ULIRGs (L IR = 10 12.0−12.3 L ⊙ ). The beam sizes of all the molecular J-transition lines obtained with different beams (different array configurations) were matched to the same value (1-2 kpc) for each (U)LIRG. Thereafter, the high-J to low-J emission line flux ratios for each molecule and emission line flux ratios among the different molecules at each J-transition were derived. Based primarily on the HCN and HCO + data which can probe the properties of dense and warm molecular gas near luminous energy sources, the observational and RADEX non-LTE modeling results were compared to constrain H 2 gas volume number density and kinetic temperature in these (U)LIRGs' nuclear 1-2 kpc regions. We obtained the following main results.
1. HCN and HCO + line emission was significantly detected at up to J=4-3, suggesting that dense and warm molecular gas is abundant in the observed (U)LIRGs' nuclei.
4. The only one observed LIRG NGC 1614 showed molecular gas density and temperature distinctly lower than those of the remaining nine ULIRGs. However, the derived gas density and temperature in all the observed (U)LIRGs' nuclei, including NGC 1614, were substantially higher than those in quiescently star-forming normal galaxies with even less infrared luminosity.
5. For selected (U)LIRGs' nuclei with bright molecular line emission, we also applied a Bayesian approach to constrain molecular gas properties, with all parameters set as free, by using the Markov Chain Monte Carlo (MCMC) method. The above arguments 3 and 4 were supported.
6. The derived density and temperature of nuclear molecular gas in starburst-dominated sources tended to be lower than those of the majority of ULIRGs that exhibited signatures of luminous AGNs from infrared, hard X-ray (>10 keV), and (sub)millimeter spectroscopic observations. This conformed to previous suggestions that luminous AGN activity was associated with a larger amount of and/or higher fraction of nuclear dense molecular gas.
This study demonstrated that high-spatial-resolution ( 1-2 kpc) ALMA observations of multiple dense molecular gas tracers at multiple J-transition lines are an effective tool for quantitatively constraining molecular gas properties of (U)LIRGs' nuclei, by minimizing possible contamination from spatially extended ( a few kpc), more diffuse and cooler molecular gas emission in the host galaxies.
We thank the anonymous referee for valuable comment which helped improve the clarity of this manuscript. This paper made use of the following ALMA data:   Table 8 summarizes the log of previously unpublished ALMA Cycle 5 observations for the two ULIRGs, IRAS 06035−7102 and IRAS 08572+3915. The obtained spectra are shown in Figure 8.    Table 9. In Figures  9 and 10, we overplot adopted Gaussian profiles on the observed data to see the goodness of the fits.  Table 9 continued on next page  Table 9 continued on next page A Certain double-peaked emission signature is observed. We also attempted a double Gaussian fit, and found that its flux was 0-10% smaller than that of the single Gaussian fit. We adopt single Gaussian fit flux.
B One emission and one absorption (dip near the emission center) components.
C Revised value from that shown by Imanishi & Nakanishi (2014), after redefining the continuum flux level ( §3).    Table 10 and Figure 11 show the continuum flux density and spectral energy distribution measured using the adopted 1-2 kpc beam, respectively. Using the continuum data, we attempted to make some correction of the possible absolute flux calibration uncertainty for individual ALMA observations (maximum ∼10%; §5.3) by assuming power law for an intrinsic continuum emission shape. However, in many (U)LIRGs, the observed continuum shape was concaved in such a way that 1.1-1.4 mm (band 6) flux was significantly smaller than the power law interpolation between 0.8-1.1 mm (band 7) and 1.4-1.8 mm (band 5) (e.g., clearly seen for IRAS 12127−1412, the Superantennae, and IRAS 20551−4250 in Figures 11f, 11i, and 11j). This is most likely because a dust thermal radiation component is significantly present in band 7 (0.8-1.1 mm) in addition to thermal free-free and synchrotron emission, which are dominant in bands 5 and 6 (1.1-1.8 mm). Thus, we decided not to apply this correction, because reliable correction turned out to be impossible.   Note-Col. (1): Object name. The adopted beam size in kpc is shown in parentheses. Col.
(2): Central frequency used for continuum extraction for J21a-J43b (explained in the caption of Table 2) in GHz is shown first, and the frequency range in GHz is listed in parentheses. The frequencies of obvious emission and absorption lines are removed. Col. (3): Continuum flux in mJy at the emission peak with adopted beam size. Value at the highest flux pixel (0. ′′ 02-0. ′′ 04 pixel −1 ) is extracted. Detection significance relative to root mean square (rms) noise is shown in parentheses. Possible systematic uncertainty arising from the absolute flux calibration ambiguity in individual ALMA observations and choice of the frequency range for continuum-level determination, are not included. Col.(4): Coordinate of continuum emission peak in ICRS.
For certain fractions of the observed ULIRGs, the photometric data at 250 µm, 350 µm, and 500 µm taken with the Herschel Space Observatory are available . These data, together with the IRAS 60 µm and 100 µm fluxes, are usually dominated by dust thermal radiation, and are plotted in Figure 12. We overplot our ALMA continuum data measured with a 1-2 kpc beam in Figure 12. Our ALMA continuum fluxes largely agree with, or are only slightly smaller than, those extrapolated from shorter-wavelength IRAS and Herschel fluxes measured using larger beams ( 5 ′′ or 5 kpc at z 0.05). This rough agreement supports the previously argued scenario that nearby ULIRGs are energetically dominated by compact ( 1-2 kpc) nuclear regions (e.g., Soifer et al. 2000;Diaz-Santos et al. 2010;Imanishi et al. 2011;Pereira-Santaella et al. 2021), which are well covered with our ALMA data. Figure 11. Continuum flux density in mJy (ordinate) as a function of observed wavelength in µm (abscissa) for the observed 11 (U)LIRGs' nuclei. 5-10% error is added in the ordinate, by considering the possible absolute flux calibration uncertainty of individual ALMA observations ( §5.3). The object name, redshift, and adopted beam size in kpc are shown. Figure 13. Same as Figure 3, but (a) for even higher HCN-to-HCO + abundance ratio of [HCN]/[HCO + ] = 7, (b) for an order of magnitude smaller HCO + column density of N HCO+ = 1 × 10 15 cm −1 , and (c) for a factor of ∼2.3 larger line width of ∆v = 800 km s −1 . The red filled circles represent the observed HCNto-HCO + flux ratios of (a,b) the same (U)LIRGs' nuclei as plotted in Figure 3 and (c) the Superantennae which displays an exceptionally large line width compared to other (U)LIRGs' nuclei (Table 3, column 11).