Exploring the hardness of the ionising radiation with the infrared softness diagram I. Similar effective temperature scales for starbursts and (ultra)luminous infrared galaxies

,


Introduction
Emission-line nebulae, including those ionised by episodes of massive star formation and/or by the presence of active galactic nuclei (AGN), can provide us with very valuable information about the several physical processes that take place in galaxies at different scales and distances.This information can help to disentangle the mechanisms governing their formation and evolution.
Bright optical emission-line intensities, of both recombination and collisional nature, can be used to determine the physical and chemical properties of the ionised gas of the objects where they are measured, given the easier access to this spectral regime from ground-based telescopes.Moreover, the emission lines in the infrared (IR) spectral regime can also supply very valuable information that can enhance and complement that obtained from the optical regime.
Among the different advantages entailed by the study of the IR lines is that they can trace larger optical depths in the ionised nebulae as extinction is attenuated in this regime.In addition, as most of the observed emission lines in the near-and mid-IR ranges have a fine-structure nature close to the ground state, their fluxes present on average much lower dependence on the electron temperature of the gas compared to optical collisionallyexcited emission lines (CELs).This dependence of the fluxes of optical CELs on temperature is what makes the determination of chemical abundances so uncertain, as the cooling rate of the gas is dominated by its overall metal content.Finally, we can find in the IR range high-excitation emission lines not measurable in the optical, which can also lead to deriving the ionic abundances of unseen ionisation stages of certain chemical elements, or to better quantifying the excitation stage of the gas in highly ionised environments, such as in the surroundings of an AGN, using lines such as [Ne v] at λ 14.3 µm and λ at 24.3 µm, or [O iv] at λ 25.9 µm (e.g., Satyapal et al. 2021).
In addition to chemical abundances or excitation, the measured emission-line fluxes in any spectral regime can be used to estimate other properties such as the hardness of the incident radiation field.This method can provide accurate answers regarding the nature of the very young ionising stellar population compared to other techniques based on the study of the stellar continuum.In this way the softness parameter, denoted as η, was defined by Vilchez & Pagel (1988) to quantify the shape of the incident spectral energy distribution (SED), relating different ratios of ionic abundances of consecutive ionic stages, corresponding to different ionisation potentials.The form of this parameter, originally defined for the optical range, is This expression presents almost no dependence on the ionisation parameter (i.e., the ratio of the number of ionising photons to the density of particles), and can trace the hardness of the incident radiation field as the involved ionisation potentials of the ions in the numerator are higher than those in the denominator, as can be seen in Fig. 1.Hence, lower values of η can be taken as a sign of a harder SED.Alternatively, an equivalent form of this parameter can be defined as a function of the corresponding emission-line fluxes (i.e., [O ii], [O iii], [S ii], [S iii] in the optical and near-IR), denoted as η , but in this case the parameter presents an additional dependence on metallicity, due to the collisional nature of the involved emission lines (e.g., Kumari et al. 2021).
In addition, several different forms of η can also be defined in the mid-IR regime based on the Ar, Ne, or S lines (Martín-Hernández et al. 2002;Morisset 2004).In particular, η IR can be defined using Ne and S emission lines in the following way: [NeII]λ12.8µm/[NeIII]λ15.6 µm [SIII]λ18.7 µm/[SIV]λ10.5 µm . (2) This form of the parameter is especially convenient for local nearby galaxies as all their involved lines are covered by Spitzer or JWST.
Moreover, the study of η IR presents additional advantages compared to other similarly defined expressions using only optical emission lines.As studied by Pérez-Montero & Vílchez (2009), observations are well covered by sequences of photoionisation models with massive single stars for different values of the equivalent effective temperature (T * ) when these are analysed in the so-called softness diagram (i.e., involving in each axis the corresponding emission-line ratios).This results seems to be contradictory with the equivalent diagram based on optical emission lines, which suggests that additional sources of ionisation, akin to having a radiation field of greater intensity than if it were solely attributed to massive stars, may be necessary to reproduce the position of some star-forming regions (e.g., HOLMES; Kumari et al. 2021).However, as analysed in Pérez-Montero et al. (2023b), the use of some low-excitation lines, such as [S ii] in this optical diagram can be contaminated by diffuse ionised gas (DIG) in regions observed using low spatial resolution, which could be the cause of the very low observed η values.This result could thus reinforce the idea that diagrams based on high-excitation lines could be more accurate as they probe regions closer to the ionising stars.
Nevertheless, using the softness parameter to analyse the hardness of the incident radiation must take into account that additional dependences on other properties of the ionised gas (e.g., metallicity, excitation, or the nature of the source) exist.For this reason, Pérez-Montero et al. ( 2019) introduced the code HII-CHI-mistry-Teff (hereafter HCm-Teff), which decomposes these dependences using a Bayesian-like comparison between the fluxes of the emission lines involved in the softness diagram, with the predictions from large grids of photoionisation models dependent on different values of T * , the ionisation parameter U, and the metallicity.In this work we describe the adaptation of this code for the IR emission lines (i.e., HII-CHImistry-Teff-ir, hereafter HCm-Teff-ir), which is particularly useful for investigating the nature of the ionising stellar population in star-forming galaxies, especially those obscured by large amounts of dust, such as (ultra)luminous infrared galaxies [(U)LIRGs].
These galaxies offer a valuable opportunity to study extreme environments and to unravel the underlying processes driving their energetic output, including very intense bursts of star formation triggered in most cases by mergers or interactions (e.g., Sanders & Mirabel 1996) leading to their extreme IR luminosities.Thus, understanding the nature of the stellar populations within these galaxies is fundamental to comprehending their ionising spectra and the mechanisms responsible for their extraordinary luminosities.
Although hugely obscured by the great amount of dust present in these objects, the synthetic spectra fitting study of the optical spectrum in (U)LIRGs, (Rodríguez Zaurín et al. 2010;Pereira-Santaella et al. 2015) shows that the ionising stellar population of their nuclei has a young mean stellar age lower than 100 Myr, and more intense than previous episodes of star formation, likely a consequence of the initial interactions preceding the final merger.The spatially resolved study of these galaxies also reveals an important contribution of DIG to the low-excitation lines from beyond the nuclear regions (Alonso-Herrero et al. 2010).This DIG component could be featured by smaller dust optical depths and heated by stars older than 100 Myr (da Cunha et al. 2010).
The study of the hardness of the ionising incident radiation using IR emission lines (less attenuated by dust absorption) and high-excitation transitions (less affected by DIG contamination) offers a unique opportunity to probe the nature of the very obscured young stellar populations in these objects.Therefore, we take advantage of the new methodology based on IR lines described in this work to perform a study of the incident radiation in a sample of (U)LIRGs dominated by star formation, to be compared with the rest of the star-forming galaxies (SFGs).
The paper is organised as follows.In Sect. 2 we describe the observational samples of compiled near-and mid-IR emission lines in the sample of star-forming galaxies and (U)LIRGs.In Sect. 3 we describe the grids of photoionisation models used to derive emission-line fluxes in the IR to construct the different versions of the softness diagram, and the code HCm-Teff for the IR to extract the corresponding values dominating this diagram for SFG.In Sect. 4 we present our results from the application of this code to the different samples of objects.Finally, in Sect. 5 we summarise our results and give our conclusions.

Data samples
We collected emission-line fluxes in the mid-and far-IR spectral ranges for star-forming galaxies and (U)LIRGs dominated by star formation from two sources.The first is Fernández-Ontiveros et al. (2021), who gathered observations in the mid-to far-IR range acquired with the InfraRed Spectrograph (IRS; Houck et al. 2004)  Fig. 1.Spectral energy distributions of some of the sources considered by photoionisation models for star-forming galaxies with the energies of the ionisation potentials of the ions whose emission lines are used in our calculations in the optical and IR ranges.
and [N iii] 57 µm far-IR lines obtained with FIFI-LS (Fischer et al. 2018) on board the SOFIA airborne observatory (Temi et al. 2018) were also taken from Peng et al. (2021) The second source includes objects from the IDEOS catalogue (Hernán-Caballero et al. 2016;Spoon et al. 2022) IR database, which includes the measurements of 77 fitted mid-IR spectral features in the 5.4-36 µm range for all galaxies observed with Spitzer (a total of 3335 galaxies; Spoon et al. 2022).From this initial sample, 153 star-forming objects (Pérez-Díaz et al. 2024) and 66 (U)LIRGs (Pérez-Díaz et al. 2023) were selected by omitting duplicate observations of the same object and entries corresponding to galaxies presented in the first catalogue, where the higher spectral resolution observation was used, and after verifying that star formation dominates the ionisation budget based on the strength of polycyclic aromatic hydrocarbon (PAH) features (EW(PAH6.2µm) > 0.06 µm) and the small contribution from certain highly ionised species ([Ne v] λ14 µm/[Ne ii] λ12 µm < 0.15 and [O iv] λ26 µm/[Ne ii] λ12 µm < 0.4).Our final sample then comprises 228 star-forming galaxies and 75 (U)LIRGs from both sources.
We collected total oxygen abundances for all objects in this selected sample as derived by Pérez-Díaz et al. (2024).These abundances were computed using version 3.2 of the code HII-CHI-mistry-ir (hereafter HCm-ir;Fernández-Ontiveros et al. 2021;Pérez-Díaz et al. 2022). 1 Unlike previous versions of HCm-ir, version 3.2 does not consider S emission lines to derive oxygen abundances as these lines are used in the new version for the derivation of sulphur abundances.This improvement of the code is described and discussed in Pérez-Díaz et al. (2024), where it is checked that this change does not imply results in the derived oxygen abundances deviating beyond the obtained errors compared to the values derived using previous versions of the 1 All versions of the code HCm can be retrieved in the web page http: //home.iaa.es/~epm/HII-CHI-mistry.html code, such as the values published by Fernández-Ontiveros et al. (2021).
The code performs a Bayesian-like calculation using certain observed emission-line ratios sensitive to the metallicity and the excitation compared with the predictions from a large grid of photoionisation models.The code gives as a result the means of the nitrogen-to-oxygen abundance ratio (N/O), the total oxygen abundance [12 + log(O/H)], and the ionisation parameter (log U), calculated as the averages of the resulting χ 2 -weighted distributions, with the corresponding uncertainties calculated as the quadratic sum of the standard deviation of the same distribution and the dispersion associated to a Monte Carlo simulation considering the observational reported errors.
For this work we compared the compiled observations with the predictions from the photoionisation code Cloudy v17.1 (Ferland et al. 2017), assuming Popstar (Mollá et al. 2009) synthesis model atmospheres.We considered as templates for the grids of models those defined in Pérez-Montero (2014) to establish empirical relations between O/H, N/O, and U used by the code in the absence of any auroral emission line.The resulting metallicity distribution was later used to constrain the final T * and U values in the softness diagram, as described in the sections below.

Description of the photoionisation models
In order to obtain an estimate of the hardness of the incident radiation field, we calculated several different grids of photoionisation models whose predicted IR line intensities were later compared by the code with the compiled IR data described in the previous section.The quantification of the hardness of the incident field of radiation was performed later using as a parameter the equivalent effective temperature (T * ), so in the models we used SEDs from single stars or black bodies, whose T * can be more easily parametrised than the synthetic massive stellar clusters used for the derivation of chemical abundances, which does not imply deviations larger than the obtained uncertainties in the results, as discussed in Pérez-Montero et al. (2019).
The photoionisation models used are the same as described in Pérez-Montero et al. (2019) adopting SEDs for massive single-stars from WM-Basic (Pauldrach et al. 2001)   planetary nebulae (PNe) SEDs from Rauch (2003) in the range of T * from 80 kK to 120 kK.
The models cover a range in total oxygen abundance from 12 + log(O/H) = 7.1 to 8.9, and log U from −4.0 to −1.5.All other chemical abundances were scaled to solar proportions, as described by Asplund et al. (2009), with the exception of N, which follows the empirical relation derived by Pérez-Montero (2014; i.e., a constant N/O ratio due to primary N production for low metallicity, and an increasing N/O due to secondary N production for higher metallicity).A particle density of 100 cm −3 and a filling factor of 0.1 were assumed in all models.The dust-to-gas mass ratio was assumed to have standard galactic proportions (d/g MW = 7.48 × 10 −3 ).In addition, to consider possible higher dust proportions in the case of (U)LIRGs, more grids were computed assuming a d/g = 2 • d/g MW .However, as in other works describing models for this dust-bounded objects (e.g., Abel et al. 2009) we did not explore other grain size distributions or depletion factors, so in our models we only refer to gas-phase abundances.All models were calculated assuming a plane-parallel and a spherical geometry to inspect the impact of this in the results.In addition, the stopping criterion of all models is that the fraction of free electrons is not lower than 98%.
The resulting emission-line fluxes were later compiled in tables in order to be compared with the observed values using the code HCm-Teff-ir, adopting the same procedure described in Pérez-Montero et al. (2019) for the case of the optical lines.

Description of the code
The use of the softness parameter based on IR lines, as defined in Eq. (2), or using alternative ratios of emission-line fluxes corresponding to consecutive states of ionisation, can be used as a proxy for the hardness of the incident radiation.The η IR parameter based on [Ne ii], [Ne iii], [S iii], and [S iv], can trace the shape of the incident SED in an energetic range quite similar to that covered by the parameter based on optical lines.In Fig. 1 it can be seen that the ionisation potentials corresponding to the emission lines involved in both the optical and IR softness parameters trace similar energetic ranges, although slightly higher as S + is replaced with S 3+ .Nonetheless, this replacement is fundamental to greatly reducing the influence of DIG emission on the new parameter based on IR lines, given the known enhancement of [S ii] emission (relative to Hα) in the DIG with respect to the H ii regions (Reynolds 1985;Domgorgen & Mathis 1994;Galarza et al. 1999).In this way, as shown by models in Fig. 2, the η IR parameter can be used to estimate the hardness of the radiation field, which decreases for increasing values of T * , as traced by WM-Basic SEDs in the 30-60 kK range and by Rauch (2003)  In addition, the parameter is also sensitive to the assumed metallicity as it is lower for lower O/H, although this dependence is much lower than in the case of the equivalent optical parameter and is only clear for high values of U.
In this scenario, a code like HCm-Teff for the IR can find a more accurate solution for T * by performing a Bayesian comparison between the observed emission-line ratios and the grid of photoionisation models.When the parameters of the models that will be used have been chosen (i.e., SED, geometry), and if metallicity is given as an input expressed in the form of total oxygen abundance, the code firstly interpolates for each entry the tables of the emission-line fluxes predicted by the models for the corresponding metallicity.Otherwise, if the metallicity is not specified, the code performs the calculation over all the possible values in the tables.
Then the code creates a distribution of T * and log U weighted by the corresponding χ 2 obtained as the sum of the quadratic difference of the observed and the predicted emission-line ratios that can be calculated with the emission-line fluxes given as input.The code gives as solutions for T * and log U the mean of the obtained weighted distributions, and as errors the quadratic sum of the corresponding standard deviation of the distributions and the dispersion of the solutions found after a Monte Carlo iteration, using as input values the nominal fluxes perturbed with the observational errors.
The code takes as input the fluxes in arbitrary units of the emission lines in the following ratios:

The IR softness parameter
We investigated how the combined use of emission-line ratios of consecutive states observable in the mid-IR can be used to determine the hardness of the incident radiation in emissionline objects ionised by massive episodes of star formation.
A40, page 4 of 10 In particular, we paid special attention to possible differences in the subsample of (U)LIRGs compared to the rest of the starforming galaxies.
As a first step we studied the behaviour of the η IR parameter, as defined in Eq. ( 2), for the samples of galaxies described in Sect. 2 in order to find possible evidence for a different hardness of the incident ionising radiation in (U)LIRGs compared with the rest of the SFGs.
We show in Fig. 3  Regarding the softness parameter, both samples show a large dispersion, but a clear offset between the two samples cannot be appreciated.This can also be easily seen by comparing the very similar mean values of log(η IR ) for the subsamples of SFGs (−0.51) and (U)LIRGs (−0.57), with a difference that is much lower than the error bars corresponding to the standard deviation of the distributions (i.e., 0.38 for SFGs and 0.28 for (U)LIRGs) overlapping in a long range.
Therefore, although the value of log η IR measured in the ULIRG sample is slightly lower than for the rest of the SFGs, which could be in principle interpreted as a harder ionising field of radiation, the difference between their means is not larger than the associated dispersions, and it does not thus seem to be significant, at least from the unique interpretation of this parameter.Moreover, this very small difference between the obtained η IR values found for (U)LIRGs compared with the other SFGs seems to be consistent with the apparent lack of a correlation with metallicity, and with that the average oxygen abundances found in each sample are similar considering the dispersion (i.e., 12 + log(O/H) = 8.48 for SFGs and 8.58 for (U)LIRGs).
However, as discussed in Sect.3, the softness parameter presents additional dependences that must be considered.For instance, as shown in Fig. 2, η IR tends to be lower for lower metallicity, although this is only noticeable for high excitation.It is thus necessary to evaluate if the obtained higher average oxygen abundance measured in (U)LIRGs could imply that the average incident field of radiation in this objects could be even harder than the trend suggested by the use of η IR alone, not significant by itself considering the uncertainties.
Therefore, in order to correctly interpret whether the observed similarity of η IR between (U)LIRGs and the rest of the SFGs is significant, it is necessary to compare this observational trend with the predictions from large grids of photoionisation models that take into account the differences in metallicity and excitation.

The IR softness diagram
In order to investigate whether the similar η IR values found in (U)LIRGs and the rest of the SFGs can be interpreted as a similar hardness of the incident radiation, we compared the observed emission-line ratios in the softness diagram (i.e., separating each emission-line ratio on a different axis of a 2D plot) to allow a better comparison with models.In the panels of Fig. 4  As can be seen from the panels of this figure, all objects including (U)LIRGs are well covered by the grid of models with T * < 60 kK, irrespective of the assumed SED, metallicity, geometry, or dust-to-gas ratio.This could be indicative that massive stars are the main contributors to the ionisation of the gas in these objects and that post-asymptotic giant branch (post-AGB) old stars do not have any preponderance for the ionisation of any of these galaxies.
In addition, the obtained differences in the sequences when different geometries or dust-to-gas ratios are considered are not significant enough to imply large deviations in the derived T * when using these sets of lines.Moreover, the inspection of both samples of objects compared to the sequences of models does A40, page 5 of 10 not seem to shed any light on the role of metallicity, given the very similar average O/H values in the samples of (U)LIRGs compared with the rest of the SFGs.On the contrary, the slightly lower η IR value found for (U)LIRGs can also be due to a different average excitation of the gas, which is easily appreciated with the different slopes of the sequences of models compared with the values for η IR , as lower values of log U can also lead to lower values for η IR .
Therefore, in order to correctly interpret if the observed difference of η IR for (U)LIRGs is relevant, it is necessary to compare this observational trend with the predictions from large grids of photoionisation models that take into account the differences in metallicity and excitation to provide an absolute scale of T * in both samples.

Results from HCm-T eff-ir
In order to obtain numerical estimations both for T * and log U, with their corresponding uncertainties, we applied the abovedescribed code HCm-Teff for the IR to the compiled sample of SFGs and (U)LIRGs with available IR emission lines.We used whenever possible all lines relevant for our calculation, not just restricting the input for the code to the emission lines involved in the softness diagram shown in Fig. 4.This implies consid-ering both [S iii] lines in the IR (at 18.7 µm and 33.5 µm) or the ratios [O iii]/[O iv] or [Ar ii]/[Ar iii].We also introduced as input for the calculations the total oxygen abundances derived using HCm-ir and described in Sect. 2. As also seen in Fig. 2, a dependence on metallicity exists, so the code can use this value to interpolate the grid of models for each object, using the corresponding uncertainty in the final dispersion reached after the Monte Carlo iteration.
In Fig. 5 we show the histograms of the resulting T * and log U obtained using HCm-Teff for the samples of SFGs and (U)LIRGs currently analysed here when the code assumes only models with WM-Basic atmospheres and a spherical geometry.The mean of these distributions, along with the 1σ dispersions are also listed in Table 1 separately for SFGs and (U)LIRGs.
In the same table we also provide the values obtained from the code when other models are assumed.As can be seen, most of the results predict a lower mean T * value for (U)LIRGs than for the rest of the SFGs, with a mean offset of around 2 kK, with the only exception of the models assuming black-body SEDs, for which the difference is lower than 1 kK.These average differences are always lower than the corresponding associated dispersions, reinforcing the idea of a similar hardness of the field of radiation in both type of objects.At the same time, the difference in log U is negligible for models with spherical geometry, A40, page 6 of 10 Table 1.Means and dispersions of the resulting T * and log U values obtained from HCm-Teff for the subsamples of SFGs and (U)LIRGs analysed here under different assumptions of SED, geometry (sph stands for spherical and pp for plane-parallel), or dust-to-gas ratio (a d/g MW value is assumed if not specified).Notes.We also list the average metallicity for each set.
while there is a systematic trend to find lower log U values for (U)LIRGs for models with plane-parallel geometry.The absolute T * scale strongly depends on the assumed models, as it increases when the models with very high T * are considered in the calculations.Nevertheless, as seen in Fig. 4, since the objects in the sample do not lie in regions of the diagram corresponding to these very high T * values, the inclusion of these sequences in the calculations can artificially increase the scale of the solutions.This scale is also dependent on the geometry, as the sequences corresponding to plane-parallel geometry predict on average lower values for η IR than for spherical.However, it is not expected that the mean T * difference found between (U)LIRGs and the rest of the SFGs can be due to different geometries in the ionisation of their nebulae, although the influence on geometry of the strong winds associated with the ionising source should also be further investigated.
Furthermore, although the average found lower absolute T * in (U)LIRGs using this procedure looks to be robust as it is independent of the assumed set in the considered models, the associated errors do not allow us to conclude that it can be real.Consequently, the result obtained from HCm-Teff-ir using IR lines supports the trend suggested by the similar value, within the errors, found for η IR in (U)LIRGs, indicative of a similar incident SED in these objects.In any case, the comparison between the observations and the sequences of models calculated for different values of metallicity and excitation help to discard the possibility of a harder field of radiation in (U)LIRGs considering their average higher metallicity.On the contrary, a log U tending to be lower and extending over a smaller range in (U)LIRGs is what it seems to contribute to their lower average η IR .
Our result pointing to a similar T * in (U)LIRGs could imply that the number of massive stars in these objects is not greater than in the rest of the SFGs, and consequently that there are not significant changes in the slope of the initial mass function (IMF) in the episodes of star formation that occurred in (U)LIRGs as a consequence of the mergers or interactions processes, contrary to some observations indicating a larger number of the most massive stars in these objects (Zhang et al. 2018).

Relationship between T * , U, and metallicity
In this subsection we study for our sample of objects the relationship between the physical conditions derived using HCm-Teff from their IR lines, including the excitation of the ionised gas or the hardness of the incident SED, and other properties in the nebulae, like their metal contents.
As described above, the metallicity of the gas can be important to convert the measured values of η IR into a quantitative scale for T * .In this way, although the difference between the mean value for η IR for (U)LIRGs (i.e., −0.57) and for the rest of the SFGs (−0.51), as can also be seen in Fig. 3, is not significant, this could be incorrectly interpreted as a harder incident SED for (U)LIRGs, considering the average greater metal content of (U)LIRGs (12 + log(O/H) = 8.56), than for SFGs (8.46).On the contrary, the Bayesian-like comparison with different grids of models that consider variations for both O/H and logU, apparently indicate that T * is similar, with a slight trend to be lower, in (U)LIRGs.
On the other hand, the relationship between O/H and T * appears to be consistent across different studies using the same methodology, based on results obtained from both optical and infrared lines.This relationship indicates that, on average, lower metallicity values correspond to higher T * values.Consequently, the observed radial gradients of increasing η in spiral galaxies (Pérez-Montero & Vílchez 2009), which have been confirmed to be caused by a hardening of the incident radiation field (Pérez-Montero et al. 2019), are correlated with gradients of decreasing metallicity as the galactocentric distance increases.This trend also appears to hold when the softness parameter is calculated using IR lines in several spiral galaxies (Pérez-Montero & Vílchez 2009).
We show in the upper and lower panels of Fig. 6 the relationship between the derived T * and log U, respectively, with the total oxygen abundance for the subsamples of (U)LIRGs and for the rest of the SFGs.The Spearman correlation coefficient A40, page 7 of 10 Fig. 6.Relation between the derived total oxygen abundance with T * (upper panel) and with log U (lower panel) for our sample of starforming galaxies (stars) and (U)LIRGs (green circles).The results were obtained using HCm-Teff assuming WM-Basic stellar atmospheres with a spherical geometry.
(ρ S ) for this latter subsample of SFGs (−0.61) is consistent with this trend of higher T * for lower O/H.However, this is not the case for (U)LIRGs, for which almost no correlation is found (ρ S = −0.07).A similar trend is found for the case of log U, which presents a negative correlation with O/H for SFGs (ρ S = −0.41),while the trend is the opposite for (U)LIRGs (ρ S = 0.30).
A correlation between the metal content of the gas and its excitation in the sense that the gas could be less excited when metallicity increases has been reported by many authors based on the optical (e.g., Pilyugin 2001) and on the IR (e.g., Giveon et al. 2002), but a physical explanation for this relation is not clear.Some authors have even proposed that the relation could be sensitive to the assumed SEDs in the models, that the analysis of the different emission-line ratios used to trace the excitation could lead to different results (e.g., Ji & Yan 2022), or that it is an effect very sensitive to geometry (Kewley et al. 2019).On the contrary, our results seem to indicate that this relation could more probably be caused by the more robust relation between T * and O/H, and the fact that most of the ratios of intensity of low-to-high excitation emission lines traditionally used to estimate U also depend on T * , as can be easily confirmed by a visual inspection of the softness diagram.The assumption of an empirical relation between O/H and U made by some recipes to derive the oxygen abundance from the relative intensities of certain collisionally excited lines (e.g., HCm, Pérez-Montero 2014) is then mostly motivated by the use of a fixed ionising SED.On the other hand, the relationship between T * and O/H could be physically well motivated by the simultaneous presence of heavy elements in the ionised gas and the massive star atmospheres or in their ejected stellar winds that are blanketed and whose UV part of the SED is attenuated, producing an average lower effective temperature.
Nonetheless, no correlation at all between T * and O/H is observed in our sample of (U)LIRGs, being especially revealing that for the subsample of four objects with very low metallicity (i.e., 12 + log(O/H) < 8.2, and in three objects even lower than 8.0) we do not obtain a higher value of T * compared to the rest of the sample of metal-rich (U)LIRGs.According to Pérez-Díaz et al. ( 2023), these metal-poor ULIRGs are experiencing a "deep diving" phase under the mass-metallicity relation (MZR; i.e., a rapid decrease in their metal content due to the dilution of the gas after a merger or interaction event).

The role of dust in T * derivation
The lack of correlation of T * with O/H in (U)LIRGs and the similar or even lower average T * in these objects compared with other SFGs could thus be an indications that another mechanism than a dependence on metallicity is responsible for the found differences.
Among other causes, the presence of large amounts of heated dust grains mixed with the surrounding gas, which is the main cause of the very high IR luminosities observed in (U)LIRGs (Sanders & Mirabel 1996), could be behind the attenuation of the ionising high-energy flux of the innermost stars if the dust is well mixed inside the galaxy.It is known that dust attenuation can be related with star formation rate (SFR) for high stellarmass galaxies (Zahid et al. 2013), and the SFRs measured in (U)LIRGs are very high as a consequence of the large amount of available gas and the intensity of the gas interchange between the galaxies taking place in the merger or interaction process (e.g., Piqueras López et al. 2016).In the case of our subsamples, the mean measured SFRs confirm this very large difference (i.e., around 100 M yr −1 for the (U)LIRGs, and 3 M yr −1 for the rest of the SFGs, according to Pérez-Díaz et al. 2023).Our result seems to validate this hypothesis as, contrary to most SFGs, no large variations in T * are found in our sample of (U)LIRGs and no direct relation between T * and O/H is found in our sample.
Given that it seems clear that (U)LIRGs are dust-bound objects (e.g., Voit 1992;Abel et al. 2003), although no clear results are found regarding a higher dust-to-gas mass ratio in these objects Herrero-Illana et al. (2019), it is pertinent to wonder if this amount of extra dust mixed with the gas, and deeply affecting both the incident field of radiation and the optical lines, can affect the derived metallicity or T * based on photoionisation codes.However, on the contrary, the use of IR lines prevents this effect as their fluxes are much less attenuated than in the case of optical lines, providing a more accurate determination of the metallicity of the gas even in regions behind a large optical depth.In addition, the consistent position of this sample of (U)LIRGs compared with other SFGs of their same stellar mass in the MZR; Pérez-Díaz et al. (2023) seem to indicate that no large variations in the final derived O/H is found when this extra innermost amount of dust is considered in the models.A40, page 8 of 10 In the case of T * , we computed additional models considering a d/g mass ratio twice the standard galactic value considered in the rest of the models to verify the role of this parameter in our results.As can be seen in panel d of Fig. 4, no large difference is seen in the softness diagram.Regarding our results with the code when we use this grid of models, as can also be seen in Table 1, there is a trend to find slightly higher T * values (∼2 kK) when these models with higher d/g are considered.Moreover, although this difference is never larger than the obtained errors, it could be even behind the lower mean T * values obtained for (U)LIRGs when the same d/g is assumed, reinforcing thus the conclusion of a similar hardness of the field of radiation in both samples, at least as seen from the IR lines emitted by the gasphase in each galaxy.

Summary and conclusions
In this work we explored the softness diagram based on mid-IR emission lines, less affected by dust obscuration and with a lower dependence on electron temperature than optical emission lines.In addition, the exploitation of this spectral range for nearby starforming galaxies allows the use of emission lines at a higher We applied this code to a sample of collected SFGs and (U)LIRGs with available IR observations.The mean observed η IR is similar within the errors in (U)LIRGs than in the rest of the SFGs, which could be interpreted as the field of radiation also being similar.This trend is confirmed by the analysis of the compiled subsamples with HCm-Teff-ir code, that reveals that the obtained mean T * is similar within the errors for this type of objects compared with the rest of the SFGs.In both cases we obtain values of T * compatible with the ionisation from massive stars and, contrary to the equivalent diagnostics based on optical emission lines, the contribution from other harder ionisation sources (e.g., post-AGB stars) is not required.
Regarding the observed relation between the obtained T * and metallicity, contrary to the rest of the SFGs for which higher T * values are found on average for lower oxygen abundances, T * does not correlate with metallicity in (U)LIRGs.This result, combined with the fact that T * is similar or even slightly lower in (U)LIRGs, contrary to expected given the very high SFRs measured in these objects, can be interpreted in terms of a much higher small heated dust grain proportion in their high SFR environments, producing a blanketing of the higher energetic range of the massive stellar populations, although this is not evidenced from the assumption of a higher dust-to-gas mass ratio in the gas-phase considered in the models, which predict a very slight increase within the errors in the mean T * .
These results could largely be reinforced with the addition of larger samples in this spectral range, combining spatially resolved observations that help to disentangle between young and old stellar populations, and the regions with a larger dust presence.An equivalent analysis for sources at higher redshift would depend on the availability of data at a longer wavelength, which unfortunately at the present moment is not covered by any facility, but at least the inclusion in our code of Ar lines in combination with Ne lines, would extend the analysis up to z ∼ 1.5 using JWST.
in the range of T * between 30 kK and 60 kK.In addition, we incorporated models from Pérez-Montero et al. (2023b) using central stars of A40, page 3 of 10
in the 80-120 kK range.However, as discussed in Pérez-Montero et al. (2019) for the optical version of the equivalent parameter based on [O ii], [O iii], [S ii], and [S iii], η presents additional dependences on other functional parameters that can lead to incorrect interpretations of objects with different η values, but also with very different metallicities or excitations.Similarly, as can be seen in the same figure, η IR also varies significantly for sequences of the same T * with different excitation, as η IR increases with log U.

Fig. 3 .
Fig. 3. Relation for the studied sample between the derived total oxygen abundance and log η IR , as defined in Eq. (2).The black stars represent SFGs and green circles (U)LIRGs.The top and right subpanels show the corresponding distributions and the mean values of each variable.The black empty histograms represent SFGs, while the hatched green histograms represent (U)LIRGs.
the obtained log η IR values as a function of the derived oxygen metallicity.The distribution of this latter derived in each sample using IR lines with HCm-ir (Fernández-Ontiveros et al. 2021) is also shown in the upper panel of the same figure, along with their mean values, while the distributions and means for log η IR are shown in the right panel.
we show this diagram in the IR, representing the emission-line ratios [Ne ii] λ 12.8 µm/[Ne iii] λ 15.6 µm and [S iii] λ 18.7 µm /[S iv] λ 10.5 µm for both SFGs and (U)LIRGs.We also show in the same figure the mean values of log(η IR ) for the subsamples of SFGs (solid black line) and (U)LIRGs (dashed line).The data is compared in the panels of Fig. 4 with different sequences of photoionisation models, with the lines joining points with the same assumed T * .Panels a, b, and d show models from 30 kK to 60 kK for WM-Basic (Pauldrach et al. 2001) stellar atmospheres and from 80 kK to 120 kK for central stars in PNe (Rauch 2003) atmospheres, while panel c shows models calculated assuming a black body with T * from 30 kK to 100 kK.As expected, the sequences with higher T * tend to be in the lower right part of the diagram, corresponding thus to lower values of η IR .The points in each sequence of T * represent different values for logU from −4.0 (upper right in each sequence) to −1.5 (lower left).The effect of considering different metallicities in the models is illustrated in panels a and c by means of solid lines for 12 + log(O/H) = 8.9, while the dashed lines represent a value of 12 + log(O/H) = 8.0.Panel b of the same figure shows model sequences calculated for 12 + log(O/H) = 8.6 and for different geometries.In this case the solid lines represent models calculated with spherical geometry, while the dashed lines represent models with plane-parallel geometry.Finally, panel d shows the effect of considering a standard dust-to-gas mass ratio in the models, represented with solid lines, in comparison with models assuming a doubled proportion of dust.

Fig. 4 .
Fig. 4. Infrared softness diagram based on the emission-line ratios [Ne ii] λ 12.8 µm/[Ne iii] λ 15.6 µm and [S iii] λ 18.7 µm /[S iv] λ 10.5 µm for the sample of objects described in the text: star-forming galaxies (black stars) and (U)LIRGs (green circles).The black solid and dashed lines represent the mean values of the η IR parameter for SFGs and (U)LIRGs, respectively.Shown in each panel are the sequences of models for different values of T * with different conditions to be compared with the data, including: a) Models calculated with WM-Basic (Pauldrach et al. 2001) massive stars (from 30 to 60 kK) and Rauch (2003) central stars of PNe (from 80 to 120 kK) assuming a spherical geometry.The solid lines represent models with 12 + log(O/H) = 8.9, while the dashed lines represent a value of 8.0.b) WM-Basic and Rauch models with 12 + log(O/H) = 8.6 with spherical (solid line) and plane-parallel (dashed) geometry.c) Models calculated assuming a black-body incident continuum from 30 kK up to 100 kK assuming spherical geometry with 12 + log(O/H) = 8.9 (solid line) and 8.0 (dashed line).d) Models using WM-Basic and Rauch stellar atmospheres with spherical geometry and 12 + log(O/H) = 8.6, and a standard value for the dust-to-gas mass ratio (solid line), and with a doubled dust content (dashed line).

Fig. 5 .
Fig. 5. Resulting distributions of the solutions found by HCm-Teff-ir for T * (upper panel) and logU (lower panel) when assuming WM-Basic SEDs and a spherical geometry.The black empty histogram represents SFGs, while the hatched green histogram represents (U)LIRGs.The average values for both distributions are also shown.
ionisation stage, such as [S iv] at λ 10.5 µm, which in turn means that this diagram has a lower dependence on diffuse emission.The IR softness parameter, η IR , based on the relative intensities of the IR lines of [Ne ii], [Ne iii], [S iii], and [S iv], traces the hardness of the incident radiation field, but it also shows a certain dependence on excitation and, although to a lesser extent than its equivalent form in the optical, on metallicity.For this reason we developed the code HCm-Teff-ir, which performs a Bayesianlike comparison between observed emission-line ratios, involved in η IR , such as [Ne ii]/[Ne iii] and [S iii]/[S iv], but also including other IR emission-line ratios such as [Ar iii]/[Ar iv] and [O iii]/[O iv], with a large grid of photoionisation models calculated assuming a central single star to provide an estimation of the equivalent effective temperature (T * ) and the ionisation parameter (log U).
on board Spitzer and with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) on board Herschel.A few detections of the [O iii] 52 µm