Near-infrared emission line diagnostics for AGN from the local Universe to z ∼ 3

Optical rest-frame spectroscopic diagnostics are usually employed to distinguish between star formation and AGN-powered emission. However, this method is biased against dusty sources, hampering a complete census of the AGN population across cosmic epochs. To mitigate this e ff ect


Introduction
Optical spectroscopy of galaxies has been for decades one of the most efficient tools for identifying Active Galactic Nuclei (AGN) in the Universe.This has been possible thanks to a combination of favorable factors.First, optical spectroscopy gives access to the rest-frame Ultraviolet (UV) to optical emission of galaxies from z = 0 to z ≃ 6.This spectral range is full of bright emission lines sensitive to the properties of the ionizing sources, and that are exploited to distinguish between ionization driven by an AGN or by young O and B stars (see Kewley et al. 2019 for a review).The best-known star-formation vs AGN diagnostic is the BPT diagram, first introduced by Baldwin et al. (1981), which informs the physical nature of the emitting source by comparing the [N ii]/Hα and [O iii]/Hβ emission line ratios.Similar diagnostics adopt [S ii]/Hα or [O i]/Hα instead of the ionized nitrogen (Veilleux & Osterbrock 1987;Kewley et al. 2001).A second reason is represented by the higher multiplexity and sensitivity reached by optical spectrographs compared to instruments sensitive to other regions of the electromagnetic spectrum.Finally, the atmosphere is mostly transparent in the visible band and thus the observations can be performed mostly from the ground.
Spectroscopy also has the advantage of differentiating the spectral type of AGN and measuring its level of activity, dis-tinguishing between Type 2 and Type 1.The first are characterized by narrow (FWHM < 1000 km/s) emission lines coming from the Narrow Line Region (NLR), and identified through the BPT and similar diagrams.Type 1, in addition to AGN-like narrow emission line ratios, also show broad component (FWHM > 1000 km/s) permitted lines.These are generated in gas clouds orbiting closer to the accretion disk and dubbed as Broad Line Region (BLR), which in Type 2 AGNs are hidden by the dusty torus (Antonucci 1993, but see also Elitzur et al. 2014 for alternative explanations).
These advantages have led in the past decades to find statistical samples of spectroscopically confirmed Type 2 and Type 1 AGNs, both in the local Universe, thanks to the SDSS survey (e.g., Kauffmann et al. 2003), and at higher redshifts, thanks to large spectroscopic campaigns conducted with multi-object spectrographs (MOS) from the ground, like zCOSMOS with VI-MOS at the VLT (Mignoli et al. 2013(Mignoli et al. , 2019)), the FMOS survey at Subaru (Kartaltepe et al. 2015;Kashino et al. 2019), the MOSDEF survey with MOSFIRE at Keck (Coil et al. 2015), and with slitless grism spectroscopy with the Hubble Space Telescope (Trump et al. 2011;Backhaus et al. 2022).Despite all these efforts, the census of AGNs from optical observations is not complete.Indeed, most theoretical models predict a black-hole accretion rate density that is higher by up to one order of magnitude compared to observational measurements (e.g., Vito et al. 2018 and references therein) at all redshifts from 0 to ∼ 7, thus suggesting that we are missing a significant population of optically obscured AGNs.One physical reason for this is that the rest-frame UV and optical emission might still be severely affected by dust attenuation.Indeed, the accretion disk onto the supermassive black-hole (SMBH), powering the AGN emission, typically lies in the central part of a galaxy where the dust attenuation might be higher (Goddard et al. 2017;Wang et al. 2017;Greener et al. 2020;Yung et al. 2021).This effect is more severe for more massive SMBHs, which are hosted in more massive and dustier galaxies.For example, A UV is of the order of ≥ 3 magnitudes in galaxies with M ⋆ > 10 10 M ⊙ (Pannella et al. 2015), and can also reach A UV ≃ 15 mag in more luminous AGNs according to Yung et al. (2021), with optical lines being highly suppressed by factors of 5 to 1000 assuming the Calzetti et al. (2000) law.There might be also systems in which the dust geometry is resembling that of a mixed model, that is, a homogeneous mixture of dust and emitting gas (Calzetti et al. 1996), rather than the typical dust-screen model.In this mixed configuration, we can find A V rising to 30 mag toward the center, as observed in local Ultra Luminous Infrared galaxies (ULIRGs) and massive starbursts at 0.5 ≲ z ≲ 1 (Calabrò et al. 2018(Calabrò et al. , 2019)).For the most compact objects, the line of sight toward the central black-hole might be completely dominated by dust residing in the host galaxy ISM at scales of ≤ 1 kpc (Calabrò et al. 2018).
To mitigate the problem of dust attenuation, a possible solution is to observe at longer wavelengths in the rest-frame nearinfrared (near-IR), where AGN radiation is less susceptible to ISM attenuation, and where our line of sight is less obstructed by dust from the BLR and NLR themselves.In the case of a mixed model, Paschen lines can recover 30% to 40% more of the total star-formation of the galaxy compared to Hα and Hβ (Calabrò et al. 2018).Since near-IR lines provide a more unbiased view of the systems, the presence of AGNs can be assessed with higher reliability, and their intrinsic properties studied more easily.
This approach has successfully discovered hidden broad-line AGNs, a type of obscured AGNs with narrow Balmer lines but with a broad (> 1000 km/s) component in Paα and Paβ (Riffel et al. 2006;Onori et al. 2017;Lamperti et al. 2017).den Brok et al. (2022) and Ricci et al. (2022) have shown that Paα and Paβ are better tracers of black-hole mass compared to Hα.Similarly, metal emission lines at λ rest > 1µm might be better diagnostics than optical ones, especially for identifying Type 2 AGNs, given that they have typically larger column densities of dust than Type 1 AGNs (Mignoli et al. 2019 and references therein).LaMassa et al. (2019) have described the limitations of the optical BPT analysis, as they find that one-half of midinfrared color selected AGNs from the WISE survey have optical emission line ratios typical of star-forming galaxies rather than of AGNs, which can be due to the highly obscured AGN clouds.Rodríguez-Ardila et al. (2004) and Colina et al. (2015) have proposed H 2 λ2.121µm/Brγ versus [Fe ii] λ1.257µm/Paβ and [Fe ii] λ1.64µm/Brγ (respectively), as useful alternatives to standard optical AGN diagnostics, and successfully applied them to study obscured AGNs in local galaxy samples.Therefore, going to longer wavelengths could potentially reveal a hidden population of BL and NL AGNs, and also better characterize their physical properties.
In this paper, we add to previous efforts in the field by exploring new AGN and star-formation (SF) diagnostics in the nearinfrared.To be applied to large samples of galaxies, we need to identify relatively bright emission lines at λ rest-frame ≳ 0.8µm, and to cover a variety of ISM conditions of the galaxies in terms of metallicity, gas density, and ionization.This can be done with the help of Cloudy photoionization models (Ferland et al. 2017), from which we can easily simulate the properties of emission lines emitted by gas clouds (representative of the galaxy ISM) irradiated by different types of ionizing sources, and then compare them to the observational results.
First, this analysis is applied at z ∼ 0, exploiting the near-IR observations and catalogs performed in the last decades and publicly available in the literature.Then, thanks to more recent spectroscopic facilities, we use these diagnostics even beyond the Local Universe.In particular, we apply them to a sample of starburst galaxies at intermediate redshifts (0.5 < z < 0.9) observed through Magellan/FIRE.Finally, we exploit the new data coming from JWST/NIRspec, which after the first light and becoming operational in 2022, has started to take spectra from low to high-resolution mode (30 < R < 2700) in the spectral range from 0.8 to 5.2µm.The JWST mission (Gardner et al. 2006(Gardner et al. , 2023) ) is thus opening a new window by shedding light on the near-IR rest-frame properties of galaxies at high redshift up to at least z ∼ 3.In this paper, we use data from the Cosmic Evolution Early Release Science Survey (CEERS) survey (Finkelstein et al., in prep.), which has collected so far among the largest samples of galaxies at 'cosmic noon' and is thus suited to analyze a sizeable sample of AGNs at z > 1.
The paper is organized as follows.In Section 2, we describe the dataset of local star-forming galaxies and AGNs, and then present the data reduction, redshift measurement and sample selection of higher redshift sources observed with Magellan/FIRE and JWST/NIRSpec.In Section 3, we describe the Cloudy photoionization models that we adopt to study the emission line properties of star-forming galaxies and AGNs, and display these models in the classical BPT diagram.We show our results in Section 4. First, we present new near-IR diagnostic diagrams that can be used to distinguish between AGN and star-forming driven ionization, providing new separation lines in analytical form.Then, we apply these diagnostics to observations from z = 0 to z = 3, and compare the near-IR classification to the optical one.We discuss our results in Section 5, suggesting ideas to improve the statistics of AGNs with forthcoming spectroscopic surveys.In our analysis, we adopt a Chabrier (2003) initial mass function (IMF) and, unless stated otherwise, we assume a cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3, Ω Λ = 0.7.

Observations and data
Given the wide range of cosmic epochs that we want to explore, from the local Universe to redshift ∼ 3, different instruments and data are required.In this section, we thus describe the various datasets that we use to study the properties of galaxies and AGNs with the same near-IR diagnostics.We start by introducing the new JWST/NIRSpec spectra coming from the CEERS survey, which provide physical information for galaxies at 1 < z < 3.At intermediate redshifts (0.5 < z < 1) we use near-IR groundbased spectra taken with the FIRE spectrograph at the Magellan telescope.We supplement this dataset with public spectral compilations of local star-forming galaxies and AGNs observed with the Spex-IRTF near-IR spectrograph and several optical instruments.0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 Fig. 1.Rest-frame spectra of 10 sources observed by JWST/NIRSpec in CEERS.The spectra are normalized to their median flux and offset relative to each other to help the visualization.They are colored in black, red, and green, according to whether they are identified, respectively, as optical+near-IR star-forming, optical+near-IR AGNs, or hidden AGNs (i.e., optical SF and near-IR AGNs), according to the diagnostic diagrams presented in this paper (see Section 4.3).The CEERS ID and the spectroscopic redshift of each source are written besides each spectrum.The emission lines that can be identified in the range 0.4µm < λ < 2µm are highlighted above the top axis (see also     2023), and will be described in more detail in Finkelstein et al. (in prep.).The galaxies analyzed in this work were observed with NIR-Spec (Jakobsen et al. 2022) during the CEERS epoch 2 observations (December 2022).The spectra were taken with the MOS configuration, which requires precise placement of the sources inside Micro Shutter Arrays (MSA; Ferruit et al. 2022) with size 0.2 ′′ ×0.46 ′′ .These NIRSpec observations were di-vided into 6 different MSA pointings, and performed with the G140M/F100LP, G235M/F170LP, and G395M/F290LP medium resolution (R ∼ 1000) gratings (herewith denoted with "M"), and with the prism in low-resolution mode (R ∼ 100).We use only the M-gratings observations in this paper, to ensure a complete deblending of the [O iii]+ Hβ and [N ii]+ Hα emission line triplets in the optical.With this spectroscopic setup, we obtain an almost continuous wavelength coverage from ∼ 1 to ∼ 5.2 µm.
The MSA apertures were made of 3-shutter slitlets, and a 3point nodding pattern was performed along the direction of the longest slitlet side to enable background subtraction, shifting the pointing by a shutter length in each direction.The sources were observed by integrating 14 groups in three exposures (one for each nod position) in NRSIRS2 readout mode, totaling 3107s of integrated time per grating for each MSA pointing.
The targeted sources were chosen among existing catalogs in EGS (Stefanon et al. 2017;Kodra et al. 2023), for which photometric redshifts were already determined in all cases, while for a small subset, the spectroscopic redshifts were also available from 3D-HST grism spectroscopy (Brammer et al. 2012).The final position of the open shutters was defined with the help of the MSA Planning Tool (MPT) to maximize the number of observed targets, preferentially selecting galaxy candidates at z > 1.A more exhaustive description of NIRSpec observations in CEERS is provided by Arrabal Haro et al. (2023).The details of the aperture position angle (PA), targeted area, and target priority criteria will be fully described in Finkelstein et al. (in prep.).

CEERS spectral reduction
The spectral reduction of NIRSpec observations was performed by the CEERS collaboration, as described in Arrabal Haro et al. (2023) and further explained in Arrabal Haro et al. (in prep.).We highlight here the main processing steps based on the JWST Pipeline (Bushouse et al. 2022) 1 and on the Calibration Reference Data System (CRDS) mapping 1061.
First, raw images are corrected for the detector 1/ f noise, dark current, and bias, and for the "snowball" contamination, which appears as trails in the images, likely produced by highenergy cosmic rays.This yields count-rate maps, from which two-dimensional (2D) spectra are created for each source.We apply to these the flat-field correction, background subtraction (using the 3-nod pattern), and photometric and wavelength calibration.The resulting spectra are also resampled to correct for spectral trace distortions.In the pipeline, a slit loss correction is already performed by default in the pathloss step.
Then, the intermediate background subtracted 2D spectra (one for each nod position) are combined to create the final 2D spectrum of the source.The one-dimensional (1D) spectra are also extracted using customized apertures that are visually determined by CEERS members by collapsing the spectra along the wavelength direction and identifying with high confidence the continuum along the spatial direction.An additional step is performed on the reduced 2D spectra to mask possible hot pixels, detector gaps, and reduction artifacts that are not associated with the main continuum trace or with an emission line.The noise spectrum is automatically calculated by the JWST pipeline, and then corrected as described in Arrabal Haro et al. (2023) to take into account the spectral resampling performed in the previous stages.

Residual flux scaling and final spectra
Although slit loss correction is performed during the main processing with the JWST pipeline (Section 2.1.1),we further check whether additional flux corrections are needed.To this aim, we convolve the previously calibrated, slit-loss-corrected 1D spectra with broadband filters and compare them to the photometric data available for each target through ground-based near-IR imaging, and/or space-based instruments such as HST/WFC3, Spitzer/IRAC, and JWST/NIRCam.This is possible because the continuum is detected with a S/N > 3 per pixel for all our NIR-Spec spectra.In particular, for all the selected sources we have HST, IRAC, or ground-based photometry, ensuring at least one band available within each grating, and a maximum number of 14 broadband filters available from λ = 1µm to 5µm.
In most cases, the shape of the spectrum is consistent with the broadband photometry across the entire wavelength range 1µm < λ < 5µm, indicating that the relative flux calibration be-tween different gratings is generally reliable.We thus multiply the spectra by a single scaling factor to match the photometric data when needed.For a few sources (9%), a different scaling correction factor should be applied to some gratings to match the photometry available in that wavelength range.These issues are not related to a specific grating, or to a particular MSA pointing.We find instead that almost all of these peculiar sources have at least one wavelength gap in the spectrum, which might be responsible for creating jumps in the flux density between different wavelength ranges.These calibration issues will be better described in the survey paper (Finkelstein et al. in prep.), and fixed in future spectroscopic reductions.
In Fig. 1, we show as an example the fully calibrated spectra (normalized to their median flux and converted to a rest-frame reference) of 10 galaxies chosen among our CEERS sample at different redshifts.This highlights the variety of wavelength coverage, continuum shapes, and spectral properties that we are able to probe.

Spectral analysis and sample selection
The spectral analysis is performed as follows.In the first stage, we visually check all the 1D spectra looking for emission lines from 1µm to 5µm.In this spectral range, there are subsets of lines that are intrinsically brighter for most galaxies, either starforming or AGNs, such as the Hβ + [O iii] triplet and Hα + [N ii] triplet in the optical, and the He i+ Paγ doublet and the Paβ + [Fe ii] triplet in the near-IR.Driven by the visual identification of these groups of lines, we assign a first spectroscopic redshift guess.
We are interested in building a sample of galaxies in the redshift range between 1 and 3. Therefore, in this phase, we select those galaxies in which the spectrum covers at least the Hβ + [O iii] triplet and the Paβ + [Fe ii] triplet.For the final sample, in addition to the two above triplets, we require the identification of the Hα + [N ii] triplet and the [S iii] λ9530Å line, as the presence of all these lines are essential requirements for the analysis of optical and near-IR diagnostics presented in this paper.With these criteria, we automatically select 65 systems between z ∼ 1 and z ∼ 3, among which 27 galaxies have z < 1.75 and thus also have Paα detected in the spectrum.This constitutes our final selected CEERS sample for all the forthcoming analyses.We also note that a few sources (10), despite being at the requested redshift, have been removed, as one of the emission line groups required above falls in detector gaps and is not available.

Line flux measurements
In addition to the emission lines considered for the first redshift guess and sample selection, we can also identify in many cases other intrinsically fainter lines.The full subset of lines relevant for this work (i.e., from λ rest = 0.4µm to λ rest = 2µm) and that we can detect for at least a subset of galaxies, is summarized in Table 1.
For a more accurate determination of the spectroscopic redshift and for the measurement of emission line fluxes, we use the χ2 minimization based MPFIT 2 routine (Markwardt 2009) to fit all the lines listed in Table 1, assuming single Gaussian functions on top of a linear continuum.We also assume for all the lines a single redshift z and the same line velocity width σ v , which are first determined by fitting the highest S/N emission line in the whole spectral range, usually Hα, or [S iii] λ9530Å in case Hα is .The emission lines are fitted with multiple gaussian components, as described in the text.Blue and orange vertical dashed lines highlight the gaussian central wavelength of permitted and forbidden lines, respectively.For broad lines, both the narrow and broad components are shown with a dashed and dotted blue line, respectively.The global best fit profile is drawn with a red dashed line.The S/N of the detection is highlighted in the top left corner of each panel, while the name of each line is written in the top right corner (see also Table 1).
Typical lines observed in star-forming galaxies and AGNs.After the first estimation of z and σ v , for the remaining lines we leave a tolerance of 100 km/s for both parameters.This provides in all cases a good fit, which we control by requiring a statistical significance of the fit of 95%, as determined from the reduced χ 2 (χ 2 red ).The quality of the relative wavelength calibration is excellent, as determined inside the CEERS collaboration, with precision at the level of less than 5% (Arrabal Haro et al., in prep).As a further test, we compared our redshifts to those determined by fitting the lines in each grating separately.These measurements were done using the methodology described in Fernández et al. (2023), for the complete sample.We find that they are consistent within 1 standard deviation, suggesting that the wavelength calibration is robust across the three gratings and the full wavelength range of NIRSpec.
In addition to the spectroscopic redshift and the line velocity width, this MPFIT procedure yields emission line fluxes, underlying continuum, equivalent widths (EW), and their corresponding uncertainties for all the lines3 .While for the first line we require a high detection threshold (S/N ≥ 5) to secure the redshift and σ v , for the remaining lines we require a S/N ≥ 2 to be detected.This choice was also adopted in Calabrò et al. (2019) to improve the detection of fainter lines while preserving its reliability, given that the fit is performed with a restricted number Fig. 3. Top : Spectroscopic redshift distribution of emission line galaxies in the CEERS survey in the redshift range between 1 and 3, where both optical and rest-frame near-IR lines can be detected, from Hβ to Paβ.The vertical dashed black line highlights the median redshift of our selected sample (z median = 2.0).Bottom : Comparison between the spectroscopic redshifts estimated from CEERS spectra and previously available photometric redshifts (Kodra et al. 2023).
of free parameters compared to a blind search.For non-detected lines, we adopt a 1σ upper limit.In Figure 2 we show as an example of our procedure the fit performed for all the emission lines of the source ID 3129 (a broad-line AGN at z = 1.04) that are relevant to this work.We show instead in Fig. 3-top the spectroscopic redshift distribution of our selected targets.They span the entire range from z ≃ 1 to z ≃ 3, with a median redshift z median = 2 and 1σ dispersion of ∼ 0.5.
We also compare the spectroscopic redshifts z spec to the photometric estimates z phot available from the latest EGS catalog assembled by Kodra et al. (2023).The z phot values are derived by fitting stellar templates to the available UV to near-IR photometry with a Hierarchical Bayesian method.We refer to Kodra et al. (2023) for the detailed photometric analysis.We find in general a good consistency, with a median absolute deviation of 0.11.We find a ∼ 10% fraction of catastrophic outliers with |z phot − z spec | > 0.5 , but no systematic under-or over-estimations across all the explored redshifts (Fig. 3-bottom).We also note that the outliers tend to have larger uncertainties of z phot compared to the whole sample.4. Distribution of intrinsic FWHM widths (i.e., deconvolved from instrumental broadening) for CEERS galaxies selected at 1 < z < 3 as described in the text.The FWHM is estimated from the Hα line.If Hα is not detected, we use the highest S/N line in the spectrum, that is, [S iii] λ9530Å or [O iii] λ5007Å.The line width is a variable parameter in the fit but is assumed to be the same value for all narrow lines, with a tolerance of 50 km/s.The median value of the sample (FWHM median ≃ 175 km/s) is represented with a vertical dashed black line.
In Fig. 4, we show the distribution of the FWHM line velocity widths from the combined Hα, [S iii] λ9530Å, and [O iii] λ5007Å lines.Given the FWHM resolution of our NIR-Spec spectroscopic setting (≃ 300 km/s) most of the lines appear well resolved.By subtracting in quadrature the instrumental resolution from the observed quantity, we calculate the intrinsic FWHM line widths.We find a median intrinsic FWHM V of 200 km/s, with standard deviation σ of 130 km/s and a tail of galaxies with higher FWHM up to ∼ 550 km/s.Furthermore, we find 3 objects for which all the permitted lines, including the Balmer, Paschen, and He i lines, require an additional broader component in the fit, with FWHM > 1000 km/s, well above the maximum found for the rest of the sample in Fig. 4.
For the Balmer and Paschen emission lines, we apply an underlying stellar absorption correction, rescaling their fluxes upwards.We apply for this correction the absorption equivalent widths (EW abs ) estimated from theoretical stellar population models, namely EW abs = 3.6, 2.7, 2.7, 2.5, and 2 Å for Hβ, Hα, Paγ, Paβ, and Paα, respectively.We note that these values for the optical lines are consistent with Domínguez et al. (2013).In our CEERS galaxies, they produce an increase of the line fluxes of, respectively, 3.4%, 0.5%, 3.7%, 2.1%, and 1% on average, thus would not alter significantly the results of this paper.
Finally, using all the hydrogen recombination lines available (typically those from the Balmer and Paschen series), we also check the quality of the final calibrated spectra by comparing their emission line ratios to the theoretical predictions, considering a variety of dust attenuation laws.While a more detailed analysis of dust attenuation will be presented in a follow-up paper (Calabro et al. 2023b, in prep.), a representative test is shown in Appendix A.1.This test suggests that Balmer and Paschen line fluxes across a large range of wavelengths (∼ 0.4 to ∼ 2µm) are physically meaningful and overall consistent with the predictions of typical dust attenuation models, even when the lines reside in different gratings, suggesting that their relative calibration can be trusted at the first order.
Article number, page 6 of 25  (2023, in prep.) if Herschel detected (red squares).The sample is representative of the star-forming Main Sequence (MS) at the median redshift of the sample (z median = 2).The black dotted line represents the starburst limit, that is, a factor of 4 above the MS of Schreiber et al. (2015).The SFRs shown in this plot are normalized to z = 2 assuming that they scale with redshift as (1 + z) 2.8 (Sargent et al. 2014).

Main physical properties of the CEERS sample
The 65 galaxies selected in Section 2.1.3show a variety of physical properties.They have H-band (F160W) magnitudes ranging between H AB = 21 and 27.2, with a median value of 24.2.Considering previous catalogs published in the CANDELS collaboration (Stefanon et al. 2017), we can also study the distribution of other main properties of the sample, including the stellar mass and the SFR.We consider here the stellar masses and SFRs obtained through fitting HST, CFHT, and Spitzer photometry (from band u to IRAC channel 4) with BC03 stellar population models, assuming a Chabrier IMF, exponentially declining SFH (τ models), and a solar metallicity (Acquaviva et al. 2012), taking also into account the contribution from emission lines to the observed photometry.While different assumptions in the SED fitting may be applied to specific objects, we aim here at investigating the global properties of the sample that can be compared to other galaxies.Nonetheless, the stellar masses that we consider are in agreement with the median stellar masses reported by Stefanon et al. (2017), which combine estimates from different groups using different prescriptions for the fit.
Representing these quantities in Fig. 5, we notice that our sample is fairly representative of the star-forming main sequence at z ∼ 2 (Schreiber et al. 2015;Whitaker et al. 2014), with stellar masses M ⋆ ranging between log M/M ⊙ = 7.8 and 11.3 (median M ⋆ ≃ 10 9.5 M ⊙ ).A small subset of CEERS galaxies is also detected at S/N ≥ 5 with Herschel and may host dust-obscured star formation, thus we use in these cases the SFRs inferred through fitting the infrared and radio photometry (from Spitzer 24µm flux to VLA 1.4 GHz flux) as described in Le Bail et al. (2023 in prep), which adopts a similar methodology to Jin et al. (2018).This subset lies in the upper right part of the diagram with higher stellar masses and higher SFRs, and total infrared luminosities L IR (integrated in the range 8-1000µm) ≳ 10 12 L ⊙ .
2.4.Magellan/FIRE spectra for starburst galaxies at z ∼ 0.7 To increase the statistics of the sample and the redshift range, we also consider 6 systems at intermediate redshifts between 0.5 and 0.9 (z median ≃ 0.7) with both optical and near-IR restframe coverage.These systems are presented in Calabrò et al. (2018) and were originally selected in the COSMOS field as infrared bright galaxies representative of a population of highly star-forming galaxies with log (M * /M ⊙ ) > 10 and log (L IR /L ⊙ ) between 11.5 and 12.5.
Their near-IR spectra were observed with the FIRE echelle spectrograph (Simcoe et al. 2013), mounted at the Magellan 6.5m Baade Telescope, in two observing runs in 2017 and 2018 (P.I.P. Cassata and R. Gobat).Thanks to their spectral coverage, from z to the K band, they include at least the Hα and the Paβ emission lines.The optical spectra are instead taken with the VIMOS spectrograph (Le Fèvre et al. 2003) and are publicly available from the zCOSMOS survey (Lilly et al. 2007), and they cover the missing rest-frame range 0.3 < λ < 0.6µm that includes the [O iii] λ5007Å and Hβ lines.
The observations, data reduction, and final calibrated spectra are fully described in Calabrò et al. (2019).The flux measurements are also presented in the above paper for the Balmer and Paschen lines, while here we also measure with the same approach the other metal lines, including and [P ii] λ1.19µm.For the subset considered here, we are able at least to detect Paβ with S/N ≥ 3 and set upper limits on [S iii], [Fe ii], [P ii], or [C i] in the worst case.Analyzing these sources in the classical BPT diagram, 5 are found to be consistent with star-forming models, while one source (dubbed BPT AGN) lies beyond the maximum separation line of Kewley et al. (2001).

A comparison sample of local starbursts and AGNs
Local samples of galaxies and AGNs can provide a useful benchmark to test and apply near-IR rest-frame diagnostics.Spectral compilations in the local Universe were made by Riffel et al. (2006), who present one of the most extensive NIR spectral atlases of AGNs.This atlas includes optical/UV and X-ray AGNs found from previous surveys, and it is representative of the entire class of AGNs at z ∼ 0, with different degrees of nuclear activity.It comprises 27 UV/optically selected sources at z < 0.1, including 12 broad line AGNs (Type 1) and 15 narrow line AGNs (Type 2), with the addition of 4 starburst galaxies that do not show evidence of nuclear activity through either broad line components, coronal lines, or AGN-like line ratios in optical diagnostics.
To increase the size of the star-forming sample for comparison, we also consider the optical and near-IR spectra of 16 infrared luminous spiral galaxies presented more recently by Riffel et al. (2019).These galaxies are all at z < 0.05 and were originally selected from the IRAS Revised Bright Galaxy Sample with L FIR > 10 10.10 L ⊙ .Their nature (AGN or star-forming) was previously determined through optical spectroscopy: 8 systems are purely star-forming (HII) galaxies, 3 are LINERs, 1 is a Type 1 AGN, and 4 systems have intermediate properties between LINER/HII or Type 2 AGN/HII.
The near-IR spectra of all these sources are taken using the Spex slit spectrograph (Rayner et al. 2003), mounted at the NASA 3m IRTF telescope at Mauna Kea.For the details of the spectroscopic observations, data reduction, and line measurements, we refer to Riffel et al. (2006Riffel et al. ( , 2019)).In the case of Type 1 AGNs, they decompose the broad and narrow component of permitted lines, so we only take their narrow component for our purposes.
In addition, we consider the sample of 11 systems from the compilation made more recently by Onori et al. (2017), and for which we have a similar rest-frame spectral coverage to the data presented before and the same subset of detected emission lines from 0.4 to 2µm rest-frame.These systems are selected as X-ray emitters from the Swift/Burst Alert Telescope 70-month catalog, and include obscured (type 2) and intermediate type AGNs, and starburst galaxies at z < 0.1.Part of this sample (8 systems) is observed at high-resolution (4000 < R < 17000) with X-shooter at the VLT.For the remaining part of this sample (3 galaxies), similarly high resolution spectra (R > 4000) are taken at LBT with the single slit NIR spectrograph LUCI (Seifert et al. 2003).We refer to Onori et al. (2017) for a complete description of the observations, data reduction, and line measurements.In particular, given the high resolution, they perform a fit with three-Gaussian components, which should be representative of the narrow component (that we use here), and eventually a broad component and an outflow component.Taking the forbidden line fluxes and the narrow component of permitted lines, it is possible to classify spectroscopically the nature of each source in the classical BPT diagram.Following this procedure, almost all the Onori et al. ( 2017) sources ( 9) are AGNs (1 Type 1, and 8 Type 2), while 1 galaxy is classified as star-forming and an additional 1 lies in the composite BPT region.
In total, we assemble a final local sample of 58 sources, including 14 Type-1 AGN, 26 Type-2 AGN, 13 starburst (purely star-forming) galaxies, 3 LINERs, and 5 sources with mixed line properties.Among the many optical and near-IR lines available in the literature for these galaxies, we select those in Table 1, which we will use in the following sections as AGN diagnostics and which can be detected also with NIRSpec in the CEERS sample.In all cases, we directly take the flux measurements available in the above papers.
We note that star-forming systems considered at z < 1 do not cover such a large variety of physical conditions as in the CEERS sample.They rather represent a subset of massive and dusty starbursts, forming stars at higher rates and possibly with higher efficiency than average galaxies in this redshift range.Despite this, they are very useful to constrain the maximum star-forming line ratios as they likely probe extreme ISM conditions among the star-forming galaxy population and the highest metallicities.

Emission line modeling with Cloudy
To provide a theoretical basis for interpreting the observations and to build a common ground from which both the optical and the near-IR emission lines of galaxies can be understood, we build a set of photoionization models with Cloudy.We use the python package pyCloudy v.0.9.11 4 which runs with version 17.01 of the Cloudy code (Ferland et al. 2017).This code assumes a symmetric distribution of gas around or in front of an ionizing source, and then calculates the radiated emission after it is processed by the gas layer.In principle, it is able to predict the emergent continuum and line intensities over the entire electromagnetic spectrum.However, we focus here on emission lines spanning from the optical rest-frame (∼ 4000 Å) to the near-IR rest-frame (∼ 2µm), which can be simultaneously observed with NIRSpec in the redshift range 1 ≲ z ≲ 3.In the following subsections, we describe the star-forming and AGN models adopted in this work. 4https://github.com/Morisset/pyCloudy/tree/0.9.11 Fig. 6.Top: Diagram with the shape of the ionizing radiation field (energy units, normalized to the flux at the Lyman limit) assumed in CLOUDY to model the stellar emission (orange curve) and the AGN emission (black, red, and darkred curves, as described in the legend).
For these representative curves we have assumed log(U) = −2, solar metallicity, and a density of 10 3 cm −3 in all cases.In the AGN models from Thomas et al. ( 2016), we identify with the same color (red or darkred) the family with a common E peak , in which the parameter p NT increases toward the top from the lighter to the darker line.The dashed shaded lines indicate the transmitted continuum.Bottom: Zoom-in of the top panel in the energy range between 5 and 100 eV, but with a linear scale on the x-axis and the flux density f ν (normalized to the flux density at the Lyman limit) on the y-axis.This version of the plot is inspired by Feltre et al. (2016).The area between the two AGN shapes of Thomas et al. (2016) with the same p NT = 0.25 and varying E peak (in the range 20-100 eV) is filled with a chequered red pattern.

Star-formation models
For a theoretical understanding of the radiation emitted by starforming galaxies, we consider an ionization-bounded, spherically symmetric shell of gas surrounded by a population of young (O and B type) stars.We note that this is a simplified picture, where a single shell of gas is representative of the ensemble of HII regions in a galaxy and the emergent spectrum is representative of the whole galaxy emission.Regarding the incident ionizing spectrum, we consider SEDs derived with BPASS stellar population models (Eldridge et al. 2017), assuming a Salpeter IMF with an upper mass limit of 300 M ⊙ , stellar metallicity equal to the gas phase metallicity of the surrounding gas cloud, and continuous star-formation history with age ranging from 10 7 to 10 9 years.For representation purposes, we show the prediction relative to a single age of 100 Myr, as other ages would not produce significant variations of the emission line ratios studied in this paper.
We specify the brightness of the incident radiation field, varying the ionization parameter log U between −4 and −1 (see  SED' command.The solar abundances and the median depletion factors of metals from the gas phase onto dust grains are set according to Savage & Sembach (1996).
The depletion factors are in agreement with those estimated by Jenkins (2009) for a depletion strength F * = 0.5.
sity n H for the gas cloud ranging from 10 2 to 10 4 cm −3 .The gasphase metallicity Z gas varies from 10% solar to two times solar, where the solar abundance of each element is set as in Savage & Sembach (1996).All the elements are scaled together with metallicity (i.e., by the same factor compared to solar), except nitrogen, which follows the prescription of Feltre et al. (2016) to take into account secondary production.The helium abundance is set as −1 in logarithmic scale (Dopita et al. 2006).We consider that some elements in the gas phase are depleted onto dust grains.The dust depletion is defined with the standard logarithmic notation system as δ(X) = [X/H] ≡ log(X/H) c − log(X/H) g , where c refers to the cosmic abundance of an ele-ment and g to its gas phase component.As shown by Savage & Sembach (1996) and other works, δ(X) can vary significantly as a function of the interstellar environment and grain composition.It is relatively well studied in our Galaxy and nearby systems, such as the Large and Small Magellanic Clouds, while it is still poorly constrained at higher redshifts.In Jenkins (2009), they parametrize the strength of the dust depletion with a factor F * varying between 0 (minimum depletion) and 1 (maximum strength), and it does not depend on the particular element.For this work, we consider average depletion coefficients as reported in Savage & Sembach (1996), which are generally consistent with those estimated from Jenkins (2009) for a medium depletion strength (F * = 0.5).The solar abundances and depletion factors for all the elements are listed in Table 2-bottom panel.We can distinguish between refractive elements with a low δ(X), such as S, O, and C, and non-refractive elements with higher δ(X), such as Fe and Si.

AGN models
The narrow-line region (NLR) of an AGN is modeled as follows.
Regarding the ionizing source, the continuum emission from the AGN accretion disk is modeled using multiple prescriptions.First, we use the built-in 'AGN' command in Cloudy, with the multiple power law continuum characterized by a 'blue bump' temperature of 10 6 K and a spectral energy index of α UV = −0.5 and α x = −1.35 in UV and X-rays, respectively, as in Risaliti et al. (2000).The spectral index α ox of the optical to X-ray SED is free to vary among the following values: −2, −1.7, −1.4,−1.2, spanning the range of observational findings, and consistent with typical values adopted in the literature (e.g., Groves et al. 2004).We also test different and more recent ionizing shapes.In particular, we consider the OXAF models5 , which are physically based AGN continuum emission models introduced by Thomas et al. ( 2016) and further described in Thomas et al. (2018).While we refer to those papers for an exhaustive discussion, we briefly present their main properties.These models are specifically designed for photoionization modeling.Indeed, they are simplified enough to easily reproduce the diversity of observed AGN spectral shapes with only three main parameters: the energy at the peak of the accretion disk emission (E peak ), the power-law index of the non-thermal emission (Γ), and the fraction of the total flux coming from the non-thermal component (p NT ).The first is a crucial parameter, as it incorporates a dependence on the black hole mass, AGN luminosity, and 'coronal' radius, and it is sensitive to contamination from shocks or HII region emission (the higher their contribution to the AGN, the lower the value of E peak ).Moreover, it is not affected by dust screening effects.In our simulations, we vary E peak from 20 to 100 eV, following Fig. 5 of Thomas et al. (2016).The parameter p NT affects the fraction of photons in the X-ray regime (see Fig. 6), and we consider three possible values: 0.1, 0.25, and 0.4, covering the entire physical range typically found in AGNs (Thomas et al. 2016).Finally, we set Γ to a median value of +2.0, as it was shown that this parameter does not significantly affect the line ratios (also confirmed in the near-IR regime).
The brightness of the radiation field is set through the ionization parameter log(U) varying from −4 to −1, as for the starforming case.Regarding the physical properties and chemical content of the NLR, we vary the gas density from 10 2 to 10 4 cm −3 , and the gas-phase metallicity from 0.3 to 3 times solar.The element abundances relative to hydrogen are set as for the starforming galaxies.The helium abundance is set slightly higher (by 0.1-0.2dex) than for star-forming galaxies, following the recent work by Dors et al. (2022).As in Feltre et al. (2016), we adopt an 'open geometry' model because it is more appropriate for the low covering fraction of the narrow-line region.
Finally, in both the star-forming and AGN case, we consider dust in the calculation and assume for simplicity the standard grain properties implemented in Cloudy (Mathis et al. 1977), with grains to hydrogen mass ratio scaling linearly with metallicity.We stop the Cloudy calculation when reaching a temperature in the gas of 500 K.In table 2, we present the full list of parameters used to model the star-forming galaxies and the NLR of AGNs.We show instead in Fig. 6 a representative ionizing spectrum and the transmitted continuum radiation for star-forming galaxies and AGNs.We can see that the spectral shape obtained with the standard 'AGN' command in Cloudy, and adopting the parameter set from Risaliti et al. (2000) (black line), does not differ significantly from more recent SEDs, at least up to the soft X-ray regime (∼ 1 keV), and can be considered as fairly representative of an AGN with median properties (E peak ≃ 60 eV, p NT = 0.25, and Γ = 2.0 in the OXAF models).In the following of the paper, we consider this as our default AGN model, but we discuss potential changes of the line ratios as due to different SEDs.For convenience, we also release the emission line predictions for the entire grid of photoionization models (both star-forming and AGN) analyzed in this paper6 .

Predictions for optical lines: the BPT diagram
We initially test the ability of our models to reproduce the classical BPT diagram, introduced originally by Baldwin et al. (1981), and widely used to separate AGN and star-formation driven ionization from the local Universe to higher redshift.This diagram compares rest-frame optical emission line ratios, namely [O iii]/Hβ and [N ii]/Hα.The predictions of our Cloudy models are shown in Fig. 7, where the circles are the AGN models for different values of log U, metallicity, and gas density.Higher gas densities are represented with bigger markers, while the marker colors are indicative of different log(U), as shown in the legend.Triangles represent the star-forming models with varying parameters of the gas shell.
Models with increasing log U tend to occupy a region toward the upper part of the diagram with a higher [O iii]/Hβ ratio.At fixed conditions, models with higher metallicity tend instead to move toward the right part of the diagram with higher [N ii]/Hα, as the nitrogen abundance also increases with Z.The composite region in Fig. 7 is where the AGN and star-forming models overlap, and where both the NLR and HII (star-forming) regions might contribute to the global emission.While a more detailed analysis of the role of each parameter in the scatter of the BPT diagram has been discussed in many other works (e.g., Brinchmann et al. 2008, Masters et al. 2016, Faisst et al. 2018, Curti et al. 2022, and references therein), we want to remark here that our AGN and star-forming models are globally consistent with observations from z = 0 to z = 3 (Kewley et al. 2013;Shapley et al. 2015;Richardson et al. 2014), and with the typical separation lines adopted in the literature to constrain the AGN location (Kauffmann et al. 2003) and the maximum starburst condition  2014), with increasing ionization from bottom to top.(Kewley et al. 2001).As expected given the similarity of ionizing shapes, the line predictions obtained with the OXAF models occupy a similar parameter space to those derived with our default AGN model.The same trends discussed above as a function of log U, Z, and gas density, also hold for the models of Thomas et al. (2016).For completeness, we show the location of these additional models in the BPT diagram in the Appendix B. We finally note that our line predictions for AGNs in the BPT are similar to those presented by Ji et al. (2020), which are also derived with Cloudy.For the star-forming models, we also checked the predictions assuming a simple black-body with temperature T = 50000 K as ionizing spectrum, yielding similar results.
The BPT diagram is also routinely used to study the ionizing properties of galaxies well beyond the local Universe.As shown by Topping et al. (2020), star-forming galaxies at 1 < z < 3 have harder ionizing spectra at higher z, lower metallicity, and increased alpha element abundances (i.e., increased O/Fe abundance ratio) compared to local samples.However, the combined effect of this evolution is rather limited in the BPT diagram, with only a slight offset observed for the star-forming population toward the upper right part of the diagram compared to the local BPT (see also Fig. 1 of Bian et al. 2018).Most importantly, they show that normal star-forming galaxies in the above redshift range do not cross the maximum starburst condition of Kewley et al. (2001) parameter, which is regulated by the star-formation history, is the main responsible for the enhancement of the [O iii]/Hβ ratio compared to local samples, but star-forming galaxies at z < 3 do not typically overcome the maximum starburst separation lines set in previous works in all the optical diagnostics (see Fig. 3 in that paper).Therefore, z ∼ 3 represents the redshift limit up to which the local AGN conditions in the BPT and similar optical diagrams have been tested.Similarly, we consider that our models can also be safely used if we extend our analysis up to at least z ∼ 3. Nevertheless, we will discuss in the following sections the possible effects of the α-enhancement and variations of the depletion strengths when considering the emission of iron-like elements.With these models in hand, we analyze the predictions of AGN and star-forming models for emission line ratios in the rest-frame near-IR.

Results
In this section, we show the near-IR rest-frame emission line predictions of our star-forming and AGN models.Then we compare them to the 64+65 observed sources in the redshift range from z = 0 to z = 3.

Infrared rest-frame predictions
We utilize the models described in Section 3 to analyze the emission line ratios in the rest-frame near-IR, from ∼ 0.8 to ∼ 2 µm.
Further extending the analysis to longer wavelengths dramatically reduces the number of objects targeted by CEERS.For example, longer wavelength transitions of the ionized hydrogen, like the Brackett emission line series, are more difficult to detect as those lines are significantly fainter than Paα, owing to theoretical Brγ and Brδ to Paα ratios of 0.06 and 0.09, respectively.We focus on near-IR lines in the above spectral range that are intrin-sically brighter, as discussed in Section 2.1.3and 2.2, and hence easier to detect for a larger number of galaxies.In particular, we have explored all combinations of close, bright near-IR emission line ratios and how they vary between purely star-forming systems and AGNs.(Curti et al. 2020) with Z gas , the latter ratio increases always monothonically, suggesting that it is more sensitive and a better tracer of metallicity in the high Z regime.The overlapping region thus includes star-forming galaxies with higher gas-phase metallicity, AGNs with more metal poor NLRs, or sources with both stellar and AGN contribution.An increase of the ionization parameter produces instead an enhancement of [S iii]/Paγ but a decrease of [Fe ii]/Paβ.However, the [S iii]/Paγ line ratio saturates at an ionization parameter log(U) ≃ −2, reaching its maximum in all conditions of density and metallicity.At higher log(U), it starts to slightly decrease, as the emission from higher ionized species of sulfur is favored in these extreme conditions.For representative purposes, we do not show in the figure the predictions from models with log(U) = −1.5 and −1, as they would overlap with the models at log(U) = −2.Finally, a higher density has a similar effect to that seen for the ionization parameter, with [S iii]/Paγ slightly increasing (and [Fe ii]/Paβ slightly decreasing) from n e = 10 2 to 10 4 cm −3 .These models shed light on new promising diagnostics of metallicity and ionization for dustenshrouded systems, which will be better investigated in the future.On the other hand, if we can independently measure these parameters from near-IR lines, the star-forming and AGN models would be even more separated, and the classification more precise.
In order to classify the galaxy spectra and provide a clear, quantitative criterion to separate sources dominated by AGN or star formation, we derive separation lines in analytical form, as already done in the optical bands in previous works.We first derive a line that encloses all the star-forming models that we have analyzed, by fitting a rectangular hyperbola to regions close and on the upper-right side of the highest metallicity SF models, for each different ionization parameter family, corresponding to triangles with the same colors in Fig. 8.This can be dubbed as the maximum star-forming (or starburst) line, defining a clear boundary for SF models located in the bottom-left corner of the diagram.Similarly, we also derive a boundary, dubbed the maximum AGN line, including the parameter space where the observed line ratios can be explained by an AGN.The points lying outside of the maximum starburst line in the AGN region can be identified as secure narrow-line AGNs.The best-fit separation lines are written in general as : where the coefficients A, B, and log X 0 depend on the emission line ratios considered.A summary including all the coefficients used in Equation 1 is presented in Table 3 for all the diagnostic diagrams.
We also present a diagnostic diagram (dubbed Fe2S3-α) that is based on the measurement of the [Fe ii] λ1.64µm/Paα line ratio on the x-axis.Even though these two lines can be detected in a limited redshift range (z ≲ 1.75) by JWST-NIRSpec (Paα is detected for less than one-half of the CEERS sample), and despite the [Fe ii] λ1.64µm being intrinsically ∼ 15% fainter than [Fe ii] λ1.257µm, they can provide useful additional contraints on the nature of dust-enshrouded sources as they probe longer wavelengths than the previous near-IR diagnostics.Nevertheless, they could be reached at higher redshifts with MIRI observations.
The model predictions for the Fe2S3-α diagnostic are shown in the second panel of Fig. 8-top as a function of [S iii]/Paγ.Considering that the involved ionized species are the same as in the Fe2S3-β diagnostic, the behavior of the models as a function of the various parameters is similar.Considering the intrinsic line ratios between Paα and Paβ, and between [Fe ii] λ1.257µm and [Fe ii] λ1.64µm, the new separation lines can be obtained by translating the Fe2S3-β corresponding lines along the x-axis by −0.35 dex, that is, by setting the X 0 coefficient 0.35 dex lower (see Table 3).
The diagnostic diagrams presented above require a measurement of Paγ.Among the recombination lines used in our diagrams, Paγ is usually the faintest for the majority of sources.For this reason, we also present an alternative set of diagrams where we use Paβ to compute the line ratio on the y-axis.These alternative diagrams are shown in Fig. 8-bottom.In contrast with the previous ones, these are slightly more dependent on dust attenuation.However, we have estimated that knowing A V with an uncertainty of ∼ 0.5 is generally enough to constrain the [S iii]/Paβ ratio with an uncertainty of 0.05 dex in case of strong obscuration (A V ≃ 5 mag), or much lower in case of moderate and low attenuation (A V ≤ 2 mag).Therefore, this dust correction does not significantly impact the classification of the source in this diagnostic.
As was done in the first case, for these diagrams we can also derive separation lines to distinguish between the emission from star-forming galaxies and the narrow-line region of AGNs.Given the intrinsic ratio between Paγ and Paβ, the new maximum starburst and maximum AGN lines are obtained by simply Finally, if we consider the AGN models developed by Thomas et al. (2016), they do not significantly alter the line ratio predictions compared to our default AGN SED in both the Fe2S3-β and Fe2S3-α diagnostics, as they tend to occupy a similar parameter space to that shown for AGNs in Fig. 8 for all the physical ranges of the three model parameters E peak , p NT , and Γ (see Table 2).The predicted lines show a similar dependence on metallicity, ionization parameter, and gas density to that explained above.The same conclusions also hold for the other near-IR diagnostic diagrams analyzed in this paper and presented in the following subsections.For completeness, we include in Appendix B a more detailed discussion about the effects of those parameters for all the near-infrared line ratios considered in this section.

Phosphorus-based AGN diagnostics
In the third panels of each row in Fig. 8, we present another useful diagnostic diagram in which we compare the [S iii]/Paγ (or [S iii]/Paβ) to the [P ii] λ1.19µm/Paβ line ratio, dubbed P2S3.As in the previous case, we can see that the separation between starforming and AGN models is rather clear, with the former occupying the lower-left region of the diagrams with lower [S iii]/Paγ ([S iii]/Paβ) and lower [P ii] λ1.19µm/Paβ, while AGNs are located in the upper right corner as due to their higher ratios.We can also see in both panels an intermediate region in the middle where the two models overlap.As for the ratio between [Fe ii] and Paschen lines, here the [P ii]/Paβ ratio also increases with metallicity, even though it spans a narrower range (∼ 2 dex), between −2.5 and −0.5 (in logarithmic scale) for star-forming galaxies, and from −1.5 to ∼ 0 for AGNs, suggesting that [P ii] is slightly less sensitive to metallicity than [Fe ii].On the other hand, sharing the same y-axis with previous diagrams, the models saturate at log(U) ≃ −2.With the same procedure adopted before, we can define both an AGN and a starburst boundary by applying Equation 1, from which we can assess the dominant ionizing source in a galaxy.The coefficients derived for the [S iii]/Paγ (or [S iii]/Paβ) vs [P ii]/Paβ diagrams are listed in Table 3.
Compared to the [Fe ii] based diagnostics, these diagrams require a measurement of the [P ii] λ1.19µm line, which is fainter than [Fe ii] and thus harder to detect, especially for individual star-forming galaxies.Indeed, the [P ii]/Paβ ratios span smaller values than [Fe ii]/Paβ.Despite this, an upper limit on this ratio would still be very useful to assess a possible AGN contribution.

Carbon-based AGN diagnostics
Finally, we identify a third diagnostic diagram which is based on the measurement of [C i] λ9850Å.This diagram compares the [S iii]/Paγ (or [S iii]/Paβ) on the y-axis to the [C i] λ9850Å/Paβ line ratio on the x-axis, and can be identified as the C1S3 diagram (fourth panels in each row of Fig. 8).The latter quantity can span six orders of magnitudes if we include both star-forming and AGN models, more than any other ratio that we have introduced before.However, in contrast to previous diagrams, this spread is mostly driven by the ionization parameter rather than the gas-phase metallicity.In star-forming galaxies, we can find values of log 10 [C i]/Paβ ≃ −4 for models with higher ionization (log(U) ∼ −2), which increase up to ≃ −0.5 for the lowest ionization models considered (log(U) = −4).In AGNs, we can have [C i] from ∼ 0.1% of Paβ at log(U) ∼ −2 up to ∼ 2.5 times brighter than Paβ at log(U) = −4.Similarly to the previous diagrams, the line ratio in both axis saturates at log(U) ∼ −2.To avoid overcrowding, we do not represent the extreme ionization models, which overlap with those at log(U) = −2.Overall, while still showing a small secondary dependence on metallicity (i.e., higher Z gas slightly increase [C i]/Paβ), [C i] is a poor tracer of the metal content and a better tracer of log(U) than [Fe ii] and [P ii], at least in the range of log(U) from −4 to −2.Similarly to the previous elements, also [C i]/Paβ decreases at higher gas densities.
Even considering all possible combinations of chemical composition, ionization parameter, and gas density, the starforming and the AGN models are in general well separated at log(U) ≤ −3, more than in any other near-IR diagnostic diagram, while overlap between SF and AGN is present at higher ionization (log(U) ∼ −2).However, as we will see later, AGNs from low to intermediate redshift do not typically populate this parameter space, preferring lower ionization parameters and higher [C i]/Paβ ratios, hence this is not an important limitation for this work.As for [P ii], [C i] is mostly detected for individual sources if they are AGNs, while in purely star-forming galaxies it can be too faint with [C i]/Paβ << 0.1.Even though a small correction should be preferably applied to the [C i]/Paβ (e.g., using the A V estimated from available Balmer and Paschen lines), the separation between the AGN and star-forming models allows us to assess with enough confidence the ionizing type of the source at log(U) ≤ 3.

Near-infrared AGN diagnostics at low and intermediate redshifts
We first check how the new near-IR diagnostic diagrams perform at z ≤ 1 and their ability to distinguish between local AGNs and sources powered by star formation.To this aim, we use the observations and the line measurements described in Sections 2.4 and 2.5.For line ratios including the Paβ line, we correct the fluxes for dust attenuation, using the value of A V estimated from a combination of Paschen and Balmer lines (i.e., from the Paα/Paβ and Paβ/Hα line ratios), or the [Fe ii] λ1.257µm and [Fe ii] λ1.64µm lines, as in Riffel et al. (2006).We adopt the median of the three estimates if all those lines are available.In We can see that, globally, the separation lines defined in Section 4.1 are able to separate the two classes of sources.Indeed, all the AGNs, either type 1 or type 2, have narrow line properties consistent with AGN models, and nearly all of them lie beyond the maximum starburst limit, in a region that cannot be explained by star formation alone.These local and intermediatez AGNs span the entire parameter range allowed by the models with different Z gas and log(U), even though the lowest ionization parameters are preferred (−4 ≲ log(U) ≲ −3) in all the four diagnostics.On the other hand, optically selected starbursts, both at z ∼ 0 and at z ∼ 0.7, are overall consistent with the star-forming  (right).Optical star-forming galaxies and AGNs are represented with blue and red circles, respectively.To be idenfitied as an AGN, a source must lie in the AGN region (within its 1σ uncertainty) in at least one of the two optical diagnostics.Cyan crosses are included to highlight hidden AGNs, that is, optical star-forming but near-IR AGNs, as explained in the text.The green dashed line is a relation representative of star-forming galaxies at z = 2.3, from Shapley et al. (2015).The violet and gray dashed lines are the classification criterion and the maximum starburst condition from Kauffmann et al. (2003) and Kewley et al. (2001), respectively.models, even though they tend to occupy the outer envelope of the SF region toward higher gas-phase metallicities (i.e., close to solar).This is somehow expected, given that they are selected as far-infrared bright sources, therefore they trace a more starforming, dustier, and more metal-enriched population compared to typical star-forming galaxies at the same redshift.As for the AGN subset, also the starbursts are more in agreement with low ionization models, that is, −4 ≲ log(U) < −3.Even though we do not have statistical samples of normal, more metal-poor star-forming galaxies at z ≤ 1 with near-IR spectral coverage, the model predictions suggest that their line ratios would span a larger region in the star-forming part toward the lower left, corresponding to lower metallicity models.We can study the normal star-forming galaxy population at higher redshifts with JWST, as we will show in the following section.We finally notice that these results are not significantly affected by the exact dust attenuation correction applied.Indeed, they would remain essentially unchanged even if we assume A V = 0 for all the sources.

Near-infrared AGN diagnostics at higher redshifts
We now study the properties of sources at z > 1.This is the redshift range where JWST-NIRSpec is revolutionizing our understanding of galaxy evolution, in particular of dusty galaxies, thanks to its longer wavelength coverage compared to its predecessors.JWST surveys are observing more and more galaxies at cosmic noon with unprecedented sensitivity and spatial resolution, and CEERS provides by far the largest spectroscopic sample of galaxies at this epoch to allow the first statistical studies in the near-IR.
Before investigating our new near-IR spectral diagnostics, we analyze the location of CEERS sources in the classical BPT and [S ii] λλ6717-6731 based diagram, to get a first understanding of their nature from the optical.The results are shown in Fig. 10 on the left and right sides, respectively.The sources are colorcoded based on their position relative to the separation lines of Kewley et al. (2001): blue circles identify the sources consistent with star-forming models (within their 1σ uncertainty), while red circles represent secure AGNs (i.e., not consistent with the star-forming region within 1σ) in at least one of the two optical diagnostics.Most of the CEERS galaxies are purely starforming.However, we also identify 8 optical AGNs (∼ 12% of the whole sample): all these sources lie in the AGN region of the [S ii] based diagram, except one system for which [S ii] is not available.Six of them are consistently identified as AGNs in the BPT.For the remaining two systems, one lies very close to the separation limit.Among the optical AGNs, we note that 3 of them also have a broad component detected at S/N > 3 in Hα or Hβ (CEERS ID 2919, 3129, 2904), thus should be considered as Type 1 AGNs, while the remaining systems are likely Type 2 AGNs in the optical (CEERS ID 2754, 5106, 12286, 16406, and 17496).Two of the broad-line AGNs (2919 and 2904) are also detected by Chandra in X-rays (Nandra et al. 2015).
It is interesting to notice that CEERS sources, both starforming galaxies and AGNs with z ranging 1 < z < 3, span almost the entire parameter range of the models.Indeed, we find sources, located in the bottom right of the SF and AGN regions, that are in agreement with low ionization models (down to log(U) ≃ −3.5) and consequently have higher metallicity, while the other galaxies and AGNs have increasingly higher log(U).Broad-line AGNs show a preference for lower ionization models compared to narrow-line AGNs, which are mostly residing in the upper left part of the BPT diagram.We note that the location of our star-forming galaxies in the BPT is broadly consistent with the relation found for typical star-forming systems at z ∼ 2.3 (Shapley et al. 2015).In this star-forming sample, galaxies with upper limits in the [N ii]/Hα ratio might have lower metallicity than those with similar [O iii]/Hβ ratios.
Moving now to the near-IR rest-frame diagnostics, we analyze the position of our galaxies in the Fe2S3-β, Fe2S3-α, P2S3, and C1S3 diagrams presented in Section 4.1.All the results are shown in Fig. 11, where we color code the sources according to their classification in the optical (red = optical AGN, blue = optical star-forming).
In all the 4 possible versions of the [Fe ii] based diagnostics (first and third rows of Fig. 11), most of the optical star-forming Fe2S3- Fig. 11.Near-IR diagnostic diagrams applied to our CEERS selected sample at 1 < z < 3. Optical star-forming galaxies and AGNs are drawn with blue and red circles, respectively.Hidden AGNs (i.e., optical star-forming but near-IR AGN) are identified as cyan crosses.Optical star-forming galaxies with upper limits are drawn without symbols in blue if their line ratios are still consistent with being star-forming.The red dashed and blue continuous lines represent the maximum AGN and starburst limits (respectively) according to our models.The analytical form of the separation lines is described in Equation 1, using the coefficients listed in Table 3.
galaxies are consistent with star-forming models in these diagrams, while also the optical AGNs are globally residing in a parameter space covered by the AGN models, that is, above the maximum AGN line in red.Therefore, these new diagnostics provide in general a consistent picture to that seen at shorter wavelengths.We also notice that, for optical AGNs, while the majority lies in the pure AGN region in the Fe2S3-β diagram, they tend to move toward the left to the intermediate region (below the maximum starburst line) in the Fe2S3-α diagnostic, possibly due to a larger contribution in the Paα line from obscured star-formation.However, we have fewer statistics in this case, as we are limited in redshift.
In the [P ii] and [C i] based diagnostics (second and fourth rows of Fig. 11), those lines are undetected for the majority of the optical star-forming galaxies, which is due to their intrinsic faintness compared to the [Fe ii] emission lines.On the other hand, we detect them for most of the optical AGNs.Also in these diagrams, the near-IR classification is mostly consistent with the optical one.Even more, we can say that for those sources in which we detect [P ii] or [C i], there is a more clear separation between optical star-forming galaxies and optical AGNs, preferentially falling in the pure star-formation or AGN region, respectively, with less overlap compared to the [Fe ii] based diagrams.
Similarly to the optical case, the CEERS sample spans a wide parameter space in the near-IR diagrams, with star-forming galaxies showing a large variety of [S iii]/Paγ or [S iii]/Paβ line ratios over almost 2 orders of magnitude, from −0.5 to 1.5 in logarithmic scale.This result further corroborates the presence of sources with different ionization parameters, from log(U) = −4 to higher values.In principle, the [C i] to Paschen line ratios can better distinguish among the highest ionization parameters, but the upper limits that we have are not sufficiently low to probe different values of log(U) in that regime.Also, the variety of [Fe ii] and [P ii] to Paschen lines flux ratios of star-forming galaxies might reflect the different metallicity properties of the sample, even though the non-detection of [P ii] for a relatively large subset hampers a more precise assessment of their properties.
For the subset of optical AGNs, even though they prefer models with low ionization parameters (≤ 3), some of them only have upper limits on Paγ, Paβ, or [C i], suggesting that they have higher log(U).The narrower range and the higher values of their [P ii]/Paβ ratios (compared to star-forming galaxies) indicate instead that they have more uniform metallicity properties, consistent with Z gas at least ≥ 0.7 times solar on average.

Hidden AGNs displaying from optical to near-infrared
In analogy with the criteria adopted in the optical, we classify as near-IR AGNs those systems that entirely reside in the pure AGN region in at least one of the 8 near-IR diagnostics, while the re-maining systems, which are always consistent with star-forming models within 1σ, are identified as near-IR star-forming galaxies.We find that all optical AGNs are also classified as near-IR AGNs.However, we also find that 5 optical star-forming galaxies become AGNs in the near-IR diagnostics (CEERS ID 2900, 8515, 8588, 8710, 9413).We identify them as hidden AGNs, and they are represented as cyan crosses in all the panels in Fig. 11.
All of these hidden AGNs are classified as pure AGNs in the P2S3 diagram.Three of them also lie in the pure AGN region in the C1S3 diagnostic, while four of them have an enhanced [Fe ii]/Paβ line ratio (either the [Fe ii] at 1.257 or 1.64 µm) and thus are identified as pure AGNs in at least one among the Fe2S3-β and Fe2S3-α diagrams.Even the sources that do not satisfy the pure AGN condition in all four diagrams, lie close to the maximum starburst line in the overlapping region between SF and AGN.
In each of the near-IR diagnostics, the optical SF galaxies whose upper limits do not allow us to put constraints on their near-IR nature (i.e., they could be consistent with either SF or AGN models), we have performed emission line measurements on their spectral stacks, and we have found that their [Fe ii], [C i], [P ii], and [S iii] line properties are similar to the other near-IR star-forming galaxies, and that they are placed in the pure starforming region in all the new diagnostics.This suggests that we are likely not missing a significant population of hidden AGNs below our line detection limit.
In Fig. 12, we show a summary of how our CEERS sources are classified both in the rest-frame optical and the near-IR.Overall, the near-IR classification is in agreement with the optical one, but with 5 optical star-forming systems that show properties similar to AGNs in the near-IR diagnostics.These hidden AGNs represent a fraction of 8.8% of the population of optical star-forming galaxies, with a 1σ confidence interval of 5.6%-11.8%,estimated following Gehrels (1986).Even though they represent a small fraction of the entire population, they nearly double the total AGN sample going from 8 to 13 objects.We also note that the optical, narrow-line AGN with CEERS ID 5106 is likely a hidden broad-line AGN, showing a broad Paβ component detected at S/N > 2.

Discussion
The results presented in Section 4 show that we can distinguish between AGN and star-forming powered sources by looking at rest-frame near-IR emission lines.In this section, we discuss the physical reason why the new diagnostics are working well for this classification task.We also discuss alternative mechanisms that might be responsible for the emission lines analyzed in this paper, and how the line ratios would evolve at high redshift.We then compare the near-IR diagnostics to the optical ones and provide possible explanations for understanding the class of hidden AGNs.Finally, we discuss further improvements that can be made with forthcoming spectroscopic facilities.

Origin of the line ratio enhancement in AGNs
The comparison between optical and near-IR predictions of Cloudy models shows us that we can consistently identify AGNs and star-forming powered sources across one order of magnitude in wavelength, using several diagnostics that involve different chemical elements, from the local Universe up to redshift ∼ 3.This also indicates that the observed lines, both in SF galaxies and in AGNs, are all consistent with being produced by photoionization.
All the near-IR diagrams analyzed here are based on the [S iii] λ9530Å to Paschen line ratios.The [S iii] λ9530Å line is one of the brightest available at around λ rest ≃ 0.95µm.It comes from the doubly ionized sulfur (S 2+ ), which is created by photons with energies between 23.3 and 34.8 eV, thus it is produced more efficiently by the harder ionizing radiation of an AGN (see Figure 6), reaching higher [S iii]/Paγ ratios compared to purely star-forming galaxies.
AGNs also have enhanced [Fe ii], [P ii], and [C i] emission compared to the typical star-forming population.The origin of this enhancement is much discussed in the literature.It is due in part to the higher metallicity of the narrow line regions, both at low and high redshift (Kewley & Dopita 2002;Groves et al. 2004), as little evolution is observed up to z ∼ 3 (Nagao et al. 2006;Matsuoka et al. 2009).As claimed by Shields et al. (2010), the strength of the optical [Fe ii] line increases almost linearly with the gas-phase iron abundance, which is observed in our models also for the [P ii] line.
However, as the AGN and SF predictions differ even at fixed Z gas , the ionization mechanism also plays an important role.The Fe + , P + , and C 0 require low ionization or excitation energies, thus usually trace transition regions of partially ionized hydrogen, coexisting with H 0 , O 0 , S + , and free electrons, where electron collisions are an important excitation mechanism.As shown in several works (Mouri et al. 1990;Oliva et al. 2001), photoionization by non-thermal continuum radiation, especially by soft X-ray photons, can create such extended regions where the [Fe ii] emission and that of other low ionization species can proliferate.The radiation field of an AGN can thus provide such favorable conditions (see also Ferland et al. 2009, Shields et al. 2010, Gaskell et al. 2022).We have verified that by switching off the soft X-ray emission in the incident radiation of the AGN template, also the emission from [Fe ii], [P ii], and [C i] is significantly reduced to the level of star-forming sources.Overall, this suggests that photoionization models can reproduce the observed emission line ratios.

Shock-models predictions
It is interesting to understand whether shocks can equally reproduce the observed line ratios in our new near-IR diagnostics.To this aim, we have explored the predictions of shock models as provided by the Mexican Million Models database (3MdB, Morisset et al. 2015).These are produced with the Mappings V code, version 5. 1.13 (Sutherland et al. 2018), using the shock parameters of Allen et al. (2008), and assuming a solar abundance for the gaseous phase, a gas density of 1cm −3 , magnetic parameter B/n 1/2 of 10 −4 to 10 µGcm 3/2 , and shock velocity between 200 and 1000 km/s.These are usually identified as fast shocks, in which the ionized front (photoionized by the cooling of hot gas) expands more rapidly than the shock itself, forming a detached HII-like region called precursor.We have analyzed the emission line predictions of both shocks and shock+precursor models in the BPT, in the Fe2S3-β, and in the C1S3 diagrams.However, we note that our star-forming galaxies are sufficiently rich in gas that we likely observe both the shock and its precursor.Seyfert galaxies are also found to be more consistent with shock+precursor models, as opposed to LINERS (Allen et al. 2008).
As we can see in Fig. 13, the shock models produce emission line ratios almost completely overlapping with the AGN models in the BPT and the C1S3 diagnostics.In contrast, they significantly enhance the [Fe ii]/Paβ line ratio by almost 1 or 2 orders of magnitude compared to AGNs and SF galaxies, respectively.This is because shocks efficiently destroy dust grains by sputtering.Since iron is one of the elements most depleted on dust in the ISM (despite being one of the most abundant in total), this process releases a large amount of iron in the gas phase, which is available to cool.
The phosphorus and [P ii] line predictions are not included in the original shock models of Allen et al. (2008).However, since it is not heavily depleted, we expect a similar behavior to [N ii] and [C i].Indeed, Storchi-Bergmann et al. (2009) and Oliva et al. (2001) have proposed the [Fe ii]/[P ii] as a main tracer of shocks.This ratio can rise to ∼ 20 in shock-dominated regions, significantly higher than the value of 2 observed in photoionized regions.Our CEERS observations and the sample of AGNs and SF galaxies at lower redshifts thus suggest that shocks can be excluded as the main contributors of the global emission.

Redshift evolution of near-infrared AGN diagnostics
Our analysis spans a large redshift range, corresponding to almost 12 Billion years of cosmic time.It is therefore worth considering possible evolutionary effects on our near-IR emission line properties.The first effect to take into account is the enhancement from low-z to high-z of the fraction of α elements (which include, among all, oxygen, carbon, and sulfur) compared to iron.This is due to galaxies at z > 1 being younger on average than those at z ∼ 0, so that type 1a supernovae, the main producers of iron, did not have time enough to enrich the surrounding ISM.
Previous studies have found abundance ratios O Fe a factor of 2 (or more) higher than the solar value (Steidel et al. 2016;Topping et al. 2020;Cullen et al. 2021).Assuming a conservative α-enhancement of 2× [α/Fe] ⊙ , we find that it would shift the maximum starburst line in the Fe2S3-β and Fe2S3-α diagrams by −0.3 dex along the x-axis.However, at z > 1, the gas-phase metallicity of typical star-forming galaxies with stellar mass ranging 9 < log(M * /M ⊙ ) < 11 is lower than at z = 0 by an amount of ∼ 0.3 dex (Wuyts et al. 2016;Sanders et al. 2020), while for AGNs it is not completely assessed.This means that the depletion factor of iron (the most depleted element) should be also lower, as the dust content is directly proportional to metallicity.Assuming a half solar gas-phase metallicity, one can estimate that the depletion factor is lower by approximately the same order of magnitude (Vladilo et al. 2011), which counterbalances the previous effect, enhancing the [Fe ii]/Paβ and [Fe ii]/Paα ratios by ∼ 0.3 dex.For this reason, as a first approximation, and considering the lack of tighter constraints at high-z, we keep at z > 1 the same local relation for separating star-forming galaxies from AGNs.
For the other near-IR diagrams, since phosphorus and carbon are almost refractive elements, the P2S3 and the C1S3 diagnostics are not significantly altered by this effect.Moreover, phosphorus behaves similarly to an α element like sulfur (Caffau et al. 2011).As a result, we adopt at z > 1 the same separation criteria between SF galaxies and AGNs as in the local Universe.
Overall, we remark that, despite the [Fe ii] lines being usually brighter than [P ii] and [C i], the Fe2S3-β and Fe2S3-α diagnostics have also the highest uncertainties among all near-IR diagrams.A higher α enhancement factor than our conservative value, as suggested by some previous works (e.g., Steidel et al. 2016 report a ratio of 4-5 × [O/Fe] ⊙ ), would explain the fraction of optical AGNs at z > 1 that in the Fe2S3-β and Fe2S3-α diagrams fall in the intermediate region.On the other hand, the ISM depletion of iron can vary by 1 order of magnitude or more (Shields et al. 2010), sufficient to explain the wide range of [Fe ii] strengths in all types of sources compared to the other diagnostics based on [P ii] or [C i].It is important to remind that each diagnostic has its pros and cons, thus we recommend combining all the available near-IR diagnostic diagrams to obtain a more accurate and complete understanding of the nature of sources at 0 < z < 3. The presence of coronal lines and broad components in permitted lines can also help in identifying AGNs.
Finally, we also note that near-IR star-forming galaxies in CEERS have typical [Fe ii]/Paβ ratios that are lower (with [Fe ii] undetected for a large fraction of them) compared to local starbursts.This is likely driven by a combination of the above-mentioned α-enhancement and selection effects: while star-forming galaxies targeted in CEERS are representative of normal star-forming galaxies at mostly M * < 10 10 M ⊙ , local samples are selected as massive (M * > 10 10 M ⊙ ) infrared bright sources with higher SFRs and higher specific SFRs, for which we might expect an enhancement of the iron abundance as due to supernovae remnants (Lester et al. 1990).

Optical versus near-infrared diagnostics and hidden AGNs
In In addition, all the metal emission lines used in the near-IR diagnostics on the x-axis likely have a similar excitation mechanism and a common origin to the low-ionization line [S ii] λλ6717-6731, as discussed in Section 5.1, while [N ii] traces slightly higher ionization regions (Mouri et al. 1990).This explains why all AGNs identified with the [S ii] diagram are also AGNs in near-IR diagnostics.
Overall, using the new near-IR diagrams has several advantages compared to the optical ones.First, the Fe2S3-β and Fe2S3-α diagrams, or the [Fe ii]/[P ii] line ratio, can distinguish between emission lines produced by photoionization or by shocks (Oliva et al. 2001), which is not possible through the classical BPT diagram (see section 5.2).The [Fe ii], [S iii], and Paβ lines, being among the brightest features in the rest-frame range 0.8-1.3µm,open the possibility to potentially characterize large samples of sources and galactic regions as shock, AGN, or starforming driven.Secondly, the wavelength separation among the near-IR lines considered in this paper is always larger than the one between Hα and [N ii] λ6583Å, or Hβ and [O iii] λ5007Å.As a consequence, while the optical lines can be resolved only with medium to high-resolution spectrographs (R≳ 600), we can measure individual near-IR lines also at lower spectral resolutions (100 ≤ R ≤ 600).Finally, and most importantly, [S iii] λ9530Å, [C i] λ9850Å, [Fe ii] λ1.257µm, [P ii] λ1.19µm, Paγ, and Paβ are significantly less affected by dust attenuation than optical lines.For example, assuming a Calzetti dust law, the attenuation (in magnitudes) at λ ≥ 1µm is a factor of 4 lower than in the range 0.4 < λ < 0.65µm (Calzetti et al. 2000).In case of optically thick dust that is well mixed with the emitting gas, the situation is even worse, with optical lines probing only ≤ 10% of the systems, while Paschen lines can recover up to 30-40% more of the total unobscured emission (Calabrò et al. 2018).
This last point provides a clue to the physical interpretation of the class of hidden AGNs.In Section 4.3.1 we have found 5 of these systems, with 4 being AGNs in the P2S3 diagram and at least one Fe-based diagnostic, while the remaining one is identified as such in P2S3 and C1S3.The host galaxies of all hidden AGNs have stellar masses greater than 10 9.3 M ⊙ .Having a median log M * /M ⊙ = 9.8, they are on average more massive than both the 5 narrow-line optical AGNs (M * ,med = 10 9.25 M ⊙ ) and than the parent population of star-forming galaxies selected at 1 < z < 3 (M * ,med = 10 9.26 M ⊙ ).Moreover, we find that in all these systems, the attenuation inferred from the Balmer decrement in the optical (adopting the Calzetti law for simplicity) is systematically lower than the one derived using the Paschen decrement (i.e., Paα + Paβ), or the Paγ/Hα and Paβ/Hα ratio if the longest wavelength line is not available.The A V estimated from Paschen lines is on average ∼ 0.3 dex higher than those obtained in the optical.This suggests that there might be optically thick dusty regions in or around the narrowline region of the AGNs, or the NLR itself might be partially covered by the dusty torus.In particular, we find that the hidden AGN ID 2900 is a Compton thick AGN with a luminosity log(L X /erg/s) = 44.07± 0.54 and a hydrogen column density N H = 10 25 cm −2 (Buchner et al. 2015).The high central density and obscuration might be the reason why it does not show broad lines and it is harder to identify compared to the other Xray AGNs in our sample (ID 2919 and 2904), which have broad lines and a lower N H < 10 23 cm −2 .To fully understand these trends and the nature of hidden systems, we would need larger statistics.A higher S/N of the spectra would also be necessary to check whether there are hidden broad-line AGNs (emerging in the Paschen lines) in the hidden AGN population, as found for the ID 5106.

Comparison with previous near-infrared diagnostics
The idea of using near-infrared rest-frame diagnostics to mitigate the problem of dust attenuation in obscured sources is not new in the literature.Rodríguez-Ardila et al. (2004) proposed the diagram comparing [Fe ii] λ1.257µm/Paβ and H 2 λ2.121µm/Brγ line ratios as a useful diagnostic to disentangle AGNs from starforming galaxies and LINERs.Since then, this diagram (that we indicate as Fe2H2) has been applied to analyze the ionizing sources in statistical samples of galaxies in the local Universe (e.g.Rodríguez-Ardila et al. 2004;Riffel et al. 2013;Lamperti et al. 2017).A different version of this diagnostic compares the [Fe ii] λ1.64µm/Brγ and H 2 λ2.121µm/Brγ line ratios (Colina et al. 2015).
The main difference of these diagnostics compared to those previously analyzed in this paper is that they are based on emission lines (H 2 and Brγ) that are in general fainter than Paα.For example, assuming Case B recombination (e.g.Osterbrock 1989), the intrinsic Brγ luminosity is 8% that of Paα and 17% that of Paβ.Furthermore, being at longer wavelengths, we can detect them in our spectra up to redshift ∼ 1.4, required to include Brγ in the reddest NIRSpec channel.In our sample, we have the coverage of Brγ in 6 galaxies at 1 < z < 1.4, of which 1 is an optical AGN (ID 3129), 1 is a hidden AGN (ID 9413), and the remaining 4 are normal star-forming galaxies from all our previous classifications.
Despite the small sample size, we have checked their position in the two Fe2H2 diagrams, which are shown in Fig. 14.We can see that the optical AGN and the hidden AGN fall in the AGN parameter space identified by Rodríguez-Ardila et al. (2004,2005) in the first diagram (top panel), and by Colina et al. (2015) in the bottom panel, with the H 2 /Brγ ratios significantly higher than 1.Regarding the star-forming galaxies, one has lower H 2 /Brγ placing it in the star-forming region ([Fe ii] λ1.64µm/Brγ is not available), while for the remaining systems H 2 and Brγ are both undetected.To conclude, the Fe2H2 diagrams give results that are consistent to those obtained with the four main near-infrared diagrams introduced in this paper, confirming our previous classification.This also suggests that the H 2 and Brγ-based diagnostics can be reliably used also at z > 1, even though a larger statistics and deeper observations are needed to confirm their effective strengths.

Future observations and facilities
At low redshift, we have been able to characterize large, representative, and almost complete samples of AGNs with various degrees of activity, and sizable subsets of near-IR bright starburst galaxies.In the high redshift Universe instead, the number of spectroscopically confirmed AGNs with near-IR observations is still small and does not cover all possible conditions of spectral type, black hole mass and obscuration, and level of activity, being limited to a few objects targeted with JWST/NIRSpec from ERS and GO public programs.
Surveys such as FRESCO (Oesch et al. 2023) will soon observe near-IR rest-frame lines in galaxies and AGNs at redshifts < 3 with NIRSpec in slitless mode.At higher redshifts, up to the reionization epoch, the brightest near-IR lines can be observed instead with JWST-MIRI, which represents a challenging frontier.Euclid will also perform slitless spectroscopy in a wavelength range between 1.25 and 1.85 µm (Costille et al. 2016), allowing to map the near-IR lines from [S iii] to Paβ at z < 0.45 with much larger statistics than JWST.Thanks to the resolving power of ∼ 450 in the red grism (Euclid Collaboration et al. 2022), star-forming galaxies and AGNs can be identified from near-IR diagnostics and broad line components in the brightest cases.
Another opportunity will be offered by the MOONS near-IR (0.6-1.8µm) spectrograph at the VLT (Cirasuolo et al. 2020;Maiolino et al. 2020), which is scheduled to start operations in 2024.This spectrograph will take spectra over a sky area of ∼ 0.15 deg 2 with ∼ 1000 fibers simultaneously, with a multiplexity performance significantly larger than previous near-IR instruments.This will allow us to detect the near-IR spectral range of thousands of galaxies up to z ∼ 0.4 (for Paβ coverage), identify all types of AGNs, and map their 3D distribution, shedding light on their environment.Finally, the Prime Focus Spectrograph (PFS) at the Subaru telescope will carry out multi-fiber spectroscopy in the range 0.38-1.3µm of 2400 targets simultaneously, within an area of 1.3 square degree (Takada et al. 2014).Despite the shorter wavelength coverage, it has a higher resolution (R ≃ 3000), which can be used to study the kinematic properties of local, near-IR selected AGNs.These spectroscopic campaigns will allow us to study the properties of both broad line and narrow line AGNs even for more obscured systems that are difficult to identify at shorter wavelengths.This will finally provide a more complete census and physical characterization of AGNs beyond the local Universe, with the same statistics reached by the SDSS at z < 0.1.
Thanks to the IFU spectroscopic mode of NIRSpec and MIRI we will be able in the future to derive spatially resolved maps of near-IR lines.At the VLT, we should also consider ERIS (Kravchenko et al. 2022), which is a near-IR imager and spectrograph, capable of medium-resolution integral field spectroscopy in J, H, and K bands, and long-slit spectroscopy in band L. It can operate in synergy with a modern adaptive optics system, necessary to spatially resolve the different galaxy components.This will enable a better understanding of the nuclear regions and determine the physical extension where the AGN significantly contributes to the observed line ratios.The resolved [Fe ii]/[P ii] ratio will be used instead to trace shock-dominated regions and outflows in the outskirts, testing the stellar and AGN feedback on galactic scales.

Summary
In this paper, we have investigated new diagnostic diagrams for selecting AGNs in the rest-frame near-infrared.We have combined model predictions with observations of star-forming galaxies and AGNs at z < 3 coming from CEERS and previous spectroscopic surveys.We summarize the main findings as follows: 1. using Cloudy photoionization models, we present four new diagnostic diagrams, named Fe2S3-β, Fe2S3-α, P2S3, and C1S3, to distinguish between stellar and AGN driven photoionization in galaxies, based only on rest-frame near-IR emission lines.These diagnostics share the same We derive in each case analytic expressions for the maximum star-forming lines and maximum AGN lines, which can be used to assess the dominant ionizing mechanism contributing to the global emission.2. these diagnostics are successfully applied from redshift ∼ 0 to z ∼ 3. The majority of the high redshift sample (z > 1) is made of star-forming galaxies, which are coherently identified as such both in the rest-frame optical and near-IR.We also find 8 optical AGNs, which are coherently classified as AGNs in the near-IR.Previously adopted diagrams including the H 2 λ2.121µm and Brγ emission lines at longer wavelengths yield results consistent with our new near-IR diagnostics.3. we find that 5 sources, classified as optical star-forming galaxies at z > 1, are identified as AGNs from near-IR diagnostics.All these sources reside in the pure AGN ionization region in the P2S3 diagram, while 2 and 4 are also identified as AGNs, respectively from the C1S3 diagram and at least one iron-based diagnostic (Fe2S3-β or Fe2S3-α).Within the uncertainties, all of them are consistent with AGN models in the entire set of near-IR diagnostics.The hidden AGNs represent a consistent fraction of the AGN population, increasing their total number by ∼ 60%.They might be systems in which the emission from the narrow-line AGN region is embedded in optically thick dust and thus can be better identified at longer wavelengths.We also find a hidden broadline AGN candidate, classified as a narrow-line AGN in the optical but for which we tentatively detect a broad Paβ component with S/N > 2.
The near-IR diagnostics presented in this paper are preferable to standard optical diagnostics for several reasons.They are less affected by dust attenuation, require a lower spectral resolution compared to the BPT diagram and the [S ii]/Hα based diagram, and finally, they can distinguish shocks from photoionization.Therefore, they represent promising tools for identifying AGNs in future spectroscopic surveys.from model predictions, indicates that the relative flux calibration of the spectra across 2µm in wavelength can be trusted on average.However, we still suggest that each object should be checked individually before computing ratios of distant lines or using their line luminosities for scientific purposes.

Appendix B: Comparison with other AGN model predictions
We analyze in this appendix all the diagnostic diagrams presented in this paper, using emission line predictions obtained with alternative AGN SEDs compared to our default spectral shape (Risaliti et al. 2000).In particular We can see that the model predictions, at fixed ionization parameter, gas density, and metallicity, occupy a similar region compared to those obtained with our default AGN SED (Fig. 7).Within the same family where those three parameters (log U, Z, and n e ) are fixed, we can notice smaller variations as due to the specific properties of the ionizing shape.The stronger variation is produced by E peak , for which we show for clarity only the predictions obtained with the extreme values considered in this work (E peak = 20 eV and 100 eV).In detail, a higher E peak yields on average a higher [O iii]/Hβ and [N ii]/Hα  7, where the line ratio predictions for AGNs are obtained using the recent models of Thomas et al. (2016Thomas et al. ( , 2018)).For AGN models with the same ionization parameter (colored as indicated in the legend), the circles are the predictions obtained with E peak = 20 eV, while the square symbols are derived assuming E peak = 100 eV.The marker size varies as a function of gas density (from 10 2 to 10 4 cm −3 from the smaller to the larger).The three points with the same marker, size, and color, are the predictions of three different values of p NT : 0.1, 0.25, and 0.4.The lines, markers, and colors for the starforming models and observations are the same as in Fig. 7.
ratios, with the increase rate depending on Z and log U.The enhancement is almost negligible at small (i.e., subsolar) metallicities and small ionization parameters (log U ≤ 3.5), but it can be up to +0.2 dex at supersolar metallicity and higher log U.This effect was already discussed in Thomas et al. (2016), and is due to more energetic and ionizing photons in SEDs with higher E peak , which mimics the enhancement of the ionization parameter and gas temperature.
On the other hand, p NT produces much smaller variations in the predicted line ratios.A higher p NT increases the fraction of energetic photons in the X-ray part of the spectrum (Fig. 6), including the soft X-rays that are more easily absorbed by the nebula.This causes a more extended partially ionized zone and more emission coming from this region, especially from low-ionization species.However, p NT does not significantly alter the photoionization model predictions in the BPT diagram, with variations of the line ratios typically smaller than 0.1 dex, even smaller than the effects due to varying gas density in the nebula.In We also analyze the predictions of the OXAF models for all the near-infrared diagrams introduced in Section 4.  3.
the third major effect on the emission line intensities.In particular, choosing models with higher E peak increases all the line ratios listed above, moving all the points toward the upper right part of the plot, by factors that depend on the metallicity, ionization parameter, and on the specific diagram.In detail, for models with log U > −2, the line ratios on both axis increase by 0.2-0.5 dex in the range 20 < E peak /eV < 100, because of more energetic ionizing photons available at higher E peak .In the last bin of log U shown in saturate with E peak , as higher ionization states of sulfur are populated.On the other hand, the ratios on the x-axis still increase, as those low-ionization species are very sensitive to the higher contribution of soft X-ray photons from models with enhanced E peak (see Fig. 6), as also discussed in Section 5. Finally, as shown by Thomas et al. (2016), the third parameter Γ has a similar behavior to p NT , but produces even smaller variations of the line ratios in all the diagrams analyzed above (both in the optical and near-IR), so we have fixed it to +2.0 in the simulations.We note that a more exhaustive discussion of the effects of each parameter can be found in the introductory paper of the OXAF models.We also remark that the exact line ratio predictions of AGN models do not affect the maximum starburst separation line (derived from star-forming models), hence they do not alter the results of this paper.

Fig. 2 .
Fig. 2. Figure displaying the spectrum of the broad-line AGN with CEERS ID 3129 at z = 1.04, around emission lines detected in the full spectral range of NIRSpec, in increasing order of wavelength.The noise spectrum is represented with gray error bars around the observed flux (black line).The emission lines are fitted with multiple gaussian components, as described in the text.Blue and orange vertical dashed lines highlight the gaussian central wavelength of permitted and forbidden lines, respectively.For broad lines, both the narrow and broad components are shown with a dashed and dotted blue line, respectively.The global best fit profile is drawn with a red dashed line.The S/N of the detection is highlighted in the top left corner of each panel, while the name of each line is written in the top right corner (see also Table1).
not present.The Hβ + [O iii] and Hα + [N ii] triplets are fitted simultaneously with the same underlying continuum, fixing the ratio between the two components of [N ii] and [O iii] to 0.338 and 0.3356, respectively (Storey & Hummer 1988).In addition, we fit together the following couples of emission lines lying close in wavelength: He i and Paγ, Paβ and [Fe ii] λ1.257µm, without any constraint on their flux ratios.

Fig.
Fig.4.Distribution of intrinsic FWHM widths (i.e., deconvolved from instrumental broadening) for CEERS galaxies selected at 1 < z < 3 as described in the text.The FWHM is estimated from the Hα line.If Hα is not detected, we use the highest S/N line in the spectrum, that is, [S iii] λ9530Å or [O iii] λ5007Å.The line width is a variable parameter in the fit but is assumed to be the same value for all narrow lines, with a tolerance of 50 km/s.The median value of the sample (FWHM median ≃ 175 km/s) is represented with a vertical dashed black line.

Fig. 5 .
Fig. 5. Diagram showing the stellar mass M ⋆ -SFR diagram of CEERS galaxies selected in this work between redshift 1 and 3.The M ⋆ and SFR are taken from Stefanon et al. (2017) as described in the text (cyan triangles) or from Le Bail et al. (2023, in prep.) if Herschel detected (red squares).The sample is representative of the star-forming Main Sequence (MS) at the median redshift of the sample (z median =2).The black dotted line represents the starburst limit, that is, a factor of 4 above the MS ofSchreiber et al. (2015).The SFRs shown in this plot are normalized to z = 2 assuming that they scale with redshift as (1 + z) 2.8(Sargent et al. 2014).

Fig. 7 .
Fig. 7. Prediction from our Cloudy models of the emission line ratios [O iii]/Hβ vs [N ii]/Hα (BPT diagram).Triangles and circles represent the star-forming and AGN models, respectively.The gas-phase metallicity increases from left to right.The blue continuous line represents the separation line by Kauffmann et al. (2003), while the red dotted line represents the maximum starburst model byKewley et al. (2001).The cyan and green dashed lines show the average location of star-forming galaxies at z = 0 and z = 2.3, respectively, fromKewley et al. (2013) andShapley et al. (2015).The big red squares and dashed line highlight the AGN sequence byRichardson et al. (2014), with increasing ionization from bottom to top.

Fig. 8 .
Fig. 8. Figure showing the near-IR diagnostic diagrams analyzed in this paper: Top row: [S iii]/Paγ as a function of [Fe ii]/Paβ (Fe2S3-β), [Fe ii]/Paα (Fe2S3-α), [P ii]/Paβ (P2S3), and [C i]/Paβ (C1S3).AGN models are represented as circles, colored as a function of log(U), and with sizes increasing with density.Star-forming models are the triangle symbols.The maximum AGN line and the maximum starburst separation limit are defined by a dotted red line and a continuous blue line, respectively.Bottom row: Same as above but as a function of the [S iii]/Paβ line ratio on the y-axis.

4. 1 . 1 .
Iron-based diagnostics for AGNsOne of the most promising near-IR diagnostic diagrams that we find is the one comparing the [S iii] λ9530Å/Paγ to the [Fe ii] λ1.257µm/Paβ line ratios, which we identify as Fe2S3-β, and is shown in Fig.8-top-left.The lines involved in each ratio are sufficiently close in wavelength that the dust corrections can be neglected in all conditions.In this parameter space, the AGN and star-forming models, represented respectively as circles and triangles, are clearly separated, with a narrow overlapping region.The star-forming galaxies occupy the bottom left region with log 10 [Fe ii] λ1.257µm/Paβ ratios spanning almost 3 orders of magnitude from ∼ −3.0 to 0, and [S iii] λ9530Å/Paγ showing a similarly large variation between −0.3 and 1.5.The AGN models occupy instead the upper-right part of the diagram with higher [Fe ii] λ1.257µm/Paβ and [S iii] λ9530Å/Paγ line ratios, up to log 10 [Fe ii]/Paβ ≃ 0.7 and log 10 [S iii]/Paγ ≃ 1.7, respectively.There is also an intermediate (also dubbed 'composite') region where the two models overlap.The objects lying in this region, while being consistent with pure stellar or pure AGN photoionization, might represent sources with intermediate properties, where we might have the contribution from both photoionization mechanisms.The theoretical modeling with Cloudy also has the advantage that we can investigate how the different physical properties of the emitting source and surrounding gas cloud influence the observed line ratios.As shown in Fig. 8-top, both the [S iii] λ9530Å/Paγ and the [Fe ii] λ1.257µm/Paβ ratios increase with metallicity, resulting from a higher abundance of coolants.However, while the former line ratio has a turning point at around Z ⊙ (depending on log(U)), which mimics the behavior of the R2 (= log [O ii] λλ3726-3729/Hβ) and S2 (= log [S ii] λλ6717-6731/Hβ) optical indices

Fig. 9 .
Fig. 9. Figures displaying the location of local and intermediate redshift AGNs and starbursts in the near-IR diagnostics defined in the paper: the Fe2S3-β, Fe2S3-α, P2S3, and C1S3 diagrams from left to right, as a function of [S iii]/Paγ on the y-axis (top row) or [S iii]/Paβ (bottom row).Sources coming from Riffel et al. (2006) and Riffel et al. (2019) are drawn with square symbols and triangles, respectively.Those presented and measured by Onori et al. (2017) are identified with empty stars, while starbursts from Calabrò et al. (2019) are shown with cyan hexagons.The underlying models are the same as described in Fig. 8.

Fig. 10 .
Fig. 10. Figure showing the position of CEERS galaxies selected at 1 < z < 3 in the BPT diagram (left) and in the [S ii]-based optical diagnostic(right).Optical star-forming galaxies and AGNs are represented with blue and red circles, respectively.To be idenfitied as an AGN, a source must lie in the AGN region (within its 1σ uncertainty) in at least one of the two optical diagnostics.Cyan crosses are included to highlight hidden AGNs, that is, optical star-forming but near-IR AGNs, as explained in the text.The green dashed line is a relation representative of star-forming galaxies at z = 2.3, fromShapley et al. (2015).The violet and gray dashed lines are the classification criterion and the maximum starburst condition fromKauffmann et al. (2003) andKewley et al. (2001), respectively.

Fig. 12 .
Fig. 12. Summary of the AGN and star-forming galaxy classification obtained from the rest-frame optical and near-IR spectroscopic diagnostics for the sample of 65 sources selected in CEERS.

Fig. 13 .
Fig. 13.Predictions of the shock models described in the text for the BPT diagram, the Fe2S3-β, and the C1S3 diagrams, from top to bottom.The AGN and SF models from Cloudy and the separation lines are the same as in Fig. 8.

Fig. A. 1 .
Fig. A.1.Diagram comparing the Hβ/Hα emission line ratio to the Paβ/Hα ratio (top panel) and the Paγ/Hα ratio (bottom panel), for CEERS galaxies in the redshift range 1 < z < 3. The predictions of different dust attenuation models are shown with colored lines.The attenuation A V along the line is specified for the Calzetti law.
, we analyze the more recent AGN spectral shapes introduced by Thomas et al. (2016), which depend on three parameters: E peak , p NT , and Γ, as discussed in Section 3.2.We first analyze the predictions for the [O iii]/Hβ and [N ii]/Hα line ratios (BPT diagram), which are shown in Fig. B.1.

Fig. B. 1 .
Fig. B.1.BPT diagram, analog of Fig.7, where the line ratio predictions for AGNs are obtained using the recent models ofThomas et al. (2016Thomas et al. ( , 2018)).For AGN models with the same ionization parameter (colored as indicated in the legend), the circles are the predictions obtained with E peak = 20 eV, while the square symbols are derived assuming E peak = 100 eV.The marker size varies as a function of gas density (from 10 2 to 10 4 cm −3 from the smaller to the larger).The three points with the same marker, size, and color, are the predictions of three different values of p NT : 0.1, 0.25, and 0.4.The lines, markers, and colors for the starforming models and observations are the same as in Fig.7.
Fig. B.1 and all subsequent figures, model predictions with different E peak are represented with different symbols (circles and squares), while predictions corresponding to different values of p NT are drawn with the same symbol and are almost overlapping.

Fig. B. 2 .
Fig. B.2. Near-infrared diagnostic diagrams including Fe2S3-β, Fe2S3-α, P2S3, and C1S3.This figure is the analog of Fig. 8, but the emission line predictions for the NLR of AGNs are obtained using the AGN models described in Thomas et al. (2016, 2018).The markers, sizes, and colors of AGN models are the same as described in Fig. B.1.The star-forming models are the same as in Fig. 8.The maximum starburst and AGN lines are defined in Section 4.1 and in Table3.

1 .
As shown byThomas et al. (2016), a similar mechanism is responsible for the higher [O i] and [S ii] emission in the optical.The influence of p NT on the same near-IR line ratios is in general smaller than that produced by E peak , with variations < 0.1 dex for models with low and intermediate ionization(log U < −2) in all the diagrams in Fig. B.2, except [C i]/Paβ,which can vary by up to ∼ 0.3 dex from p NT = 0.1 to 0.4.Also for this parameter, at the highest log U considered here, it mostly affects the line ratios on the x-axis by the same mechanism explained above for E peak , with increments of up to ∼ 1 dex for [C i]/Paβ and ∼ 0.5 dex for the remaining ratios.We note that [C i]/Paβ is more sensitive to variations of E peak and p NT because it is entirely produced in partially ionized regions where electron collisions stimulated by non-thermal continuum radiation are the main excitation mechanism.

Table 1 .
Table with the list of optical and near-IR rest-frame lines and doublets that we are able to detect in our spectra for star-forming galaxies and AGNs between z = 0 and z = 3, and that we consider in this paper.

Table 2 -
top panels), while we assume a hydrogen number den-

Table 2 .
Table Risaliti et al. (2000) parameters used in Cloudy for constructing the grid of star-forming models and AGN models explored in this work.The AGN models are made following both the approach ofRisaliti et al. (2000), which is based on the standard 'tableAGN'  command in Cloudy, and of Thomas et al. (2016), with the custom SED shapes injected in pyCloudy through the 'table the near-IR versions of the AGN diagnostics, [S iii] λ9530Å replaces [O iii] λ5007Å on the y-axis, while [Fe ii] λ1.257µm, [Fe ii] λ1.64µm, [P ii] λ1.19µm, and [C i] λ9850Å are used on the x-axis instead of [N ii] λ6583Å and [S ii] λλ6717-6731.Comparing the [O iii]to the [S iii] emission, the production of O 2+ requires photons with energies ranging 35.1-54.9eV, significantly higher than those required to create S 2+ .At those high energies indeed, a large fraction of sulfur in the gas phase would be triply ionized (S 3+ ), lowering the abundance of S 2+ .This implies that [S iii]/Paβ is not a good tracer of the ionization parameter as the [O iii]/Hβ line ratio, as it does not have a monotonic trend over the full range of possible log(U) from −4 up to = −1.The models are thus degenerate, especially at high ionization.However, an advantage is that star-forming galaxies at higher redshifts, despite having typically higher log(U), would not have also extreme [S iii]/Paβ line ratios, hence they likely would not cross the maximum starburst lines even at z > 3. The remaining drawback is that the AGNs can move to the intermediate part of all the diagnostics (below the maximum star-forming line) if they have metallicities lower than one-half solar.As far as the other lines are concerned, both [Fe ii] and [P ii] (except [C i]), are very sensitive to metallicity, as [N ii] and [S ii] in the optical regime.