Fitting spectral energy distributions of FMOS-COSMOS emission-line galaxies at z$\sim$1.6: Star formation rates, dust attenuation, and [OIII]$\lambda$5007 emission-line luminosities

We perform SED fitting analysis on a COSMOS sample covering UV-to-FIR wavelengths with emission lines from the FMOS survey. The sample of 182 objects with H$\alpha$ and [OIII]$\lambda5007$ emission spans over a range of $1.40<\rm{z}<1.68$. We obtain robust estimates of stellar mass ($10^{9.5}-10^{11.5}~\rm{M_\odot}$) and SFR ($10^1-10^3~\rm{M_\odot}~\rm{yr}^{-1}$) from the Bayesian analysis with CIGALE fitting continuum photometry and H$\alpha$. We obtain a median attenuation of A$_\rm{H\alpha}=1.16\pm0.19$ mag and A$_\rm{[OIII]}=1.41\pm0.22$ mag. H$\alpha$ and [OIII]$\lambda5007$ attenuations are found to increase with stellar mass, confirming previous findings. A difference of $57$% in the attenuation experienced by emission lines and continuum is found in agreement with the lines being more attenuated than the continuum. New CLOUDY HII-region models in CIGALE enable good fits of H$\alpha$, H$\beta$, [OIII]$\lambda5007$ emission lines with differences smaller than $0.2$ dex. Fitting [NII]$\lambda6584$ line is challenging due to well-known discrepancies in the locus of galaxies in the BPT diagram at intermediate redshifts. We find a positive correlation for SFR and dust-corrected L$_\rm{[OIII]\lambda5007}$ and we derive the linear relation $\log_{10}\rm{(SFR/\rm{M}_\odot~\rm{yr}^{-1})}=\log_{10} (\rm{L}_{[\rm{OIII]}}/\rm{ergs~s^{-1}})-(41.20\pm0.02)$. Leaving the slope as a free parameter leads to $\log_{10}\rm{(SFR/\rm{M}_\odot~\rm{yr}^{-1})}=(0.83\pm0.06)\log_{10}(\rm{L}_{[\rm{OIII]}}/\rm{ergs~s^{-1}})-(34.01\pm2.63)$. Gas-phase metallicity and ionization parameter variations account for a $0.24$ dex and $1.1$ dex of the dispersion, respectively. An average value of $\log\rm{U}\approx-2.85$ is measured for this sample. Including HII-region models to fit simultaneously photometry and emission line fluxes are paramount to analyze future data from surveys such as MOONS and PFS.


Introduction
The spectral energy distribution (SED) of a galaxy from the ultra-violet (UV) to the far-infrared (far-IR) reflects its stellar populations and the interplay of their emitted light with gas and dust in the interstellar medium (ISM). The UV to near-infrared (NIR) emission is the result of stellar populations created over a galaxy's lifetime. Dust absorption of the stellar light leaves its imprint on the gas and the stellar continuum and also via the re-emission of the heated dust in the mid and far-IR. The ionized gas component produced by the high-energy flux of massive stars is studied through the nebular emission of the galaxy, providing access to the very recent star formation of massive stars and physical conditions inside HII-regions.
As a consequence, SEDs over a large range of wavelengths and combining photometric data and emission lines give us access to crucial quantities, such as the current star formation rate (SFR), stellar mass, star formation history (SFH) over different timescales, dust attenuation, dust content, or stellar and gas metallicities. However, the coupling of stars, gas, and dust in a galaxy makes difficult the extraction of the information for each E-mail: jorge.villa@lam.fr component. The comparison of models predicting the full SED to observed fluxes from the continuum and line emission has been proven to be very powerful to infer these physical parameters in star-forming galaxy populations (Fossati et al. 2018;Buat et al. 2018;Corre et al. 2018;Yuan et al. 2019). For such comparisons, different categories of models are used. The SFH is either modeled with simple analytical functions or non-parametric forms, or it comes from numerical simulations (semi-analytical or hydrodynamical). The interplay between dust and stars is described either using a radiation transfer modeling with more or less complex configurations or with simpler, phenomenological laws that can sometimes be combined with an energy budget (Silva et al. 1998;Popescu et al. 2000;da Cunha et al. 2008;Boquien et al. 2019). The nebular component is added either using physical modeling with codes such as CLOUDY or MAPPINGS (Ferland et al. 2017;Allen et al. 2008), as in Code Investigating GALaxy Emission (CIGALE 1 ) , BayEsian Analysis of GaLaxy sEds (BEAGLE) (Chevallard & Charlot 2016), ProSpect package (Robotham et al. 2020), Python code for Stellar Population Inference from Spectra and SEDs (Prospector) (Leja et al. 2017;Johnson et al. 2021), Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (Bagpipes) (Carnall et al. 2018), MCSED (Bowman et al. 2020), or with empirical relations relating different emissions as in PHotometric Analysis for Redshift Estimations (lePHARE) (Arnouts et al. 1999;Ilbert et al. 2006). Some SED fitting codes are capable of combining continuum and emission line fluxes when the analysis of data obtained with new facilities will require the handling of very large data sets and the modeling of very different galaxy populations (Ellis et al. 2017; Thorne et al. 2021).
The IR emission is crucial with regard to putting constraints on stellar obscuration  and to estimate physical quantities as the SFR of galaxies from the energy budget (Smith et al. 2012;Małek et al. 2018;Nersesian et al. 2019;Dobbels et al. 2020). We note that the modeling of the infrared (IR) dust emission is not included in all the SED fitting codes; CIGALE includes three different sets of models from Draine & Li (2007), with updates of Draine et al. (2014), Casey (2012), and Dale et al. (2014), while Prospector, Bagpipes, ProSpect and BEAGLE only include one of the aforementioned sets. Moreover, active galactic nuclei (AGN) templates from Fritz et al. (2006), Casey (2012), and Dale et al. (2014) are also available in CIGALE and ProSpect where the former also includes Andrews et al. (2018) templates. Nevertheless, emission line models are still not available for use with the AGN templates. Large spectro-photometric surveys on 8-meter telescopes will provide thousands to millions of spectra in the deepest photometric fields, such as with the Multi-Object Optical and Near-infrared Spectrograph (MOONS) for the Very Large Telescope (VLT) and the Subaru Prime Focus Spectrograph (PFS), and larger facilities (e.g., Extremely Large Telescope) such as the Multi-Object Spectrograph for Astrophysics, Intergalacticmedium studies, and Cosmology (MOSAIC) . The James Webb Space Telescope (JWST) will push the observations of both restframe UV-optical continuum and emission lines to very high redshift in particular the Hα and [OIII]λ5007 emission lines for which different physical conditions and configurations will need to be taken into account in their modeling (Schaerer & de Barros 2009;Wright et al. 2015;Wells et al. 2015;Álvarez-Márquez et al. 2019;Chevallard et al. 2019). While treating both photometric and emission line information simultaneously, it is important to account for the potential contribution of emission lines in photometric bands, which can substantially modify the observed fluxes increasing the uncertainties in the parameter estimations (Schaerer & de Barros 2010;Stark et al. 2013;Tang et al. 2021). The coupling of multi-wavelength data sets in the SED fitting including emission lines and a good treatment of the IR emission remains paramount in the work to characterize, derive, and accurately measure the properties of galaxies.
In this paper, we attempt to simultaneously fit photometric and spectroscopic data with CIGALE in the Cosmic Evolution Survey (COSMOS) field using photometry from Laigle et al. (2016) and emission line fluxes from the Fiber Multi-Object Spectrograph (FMOS-COSMOS) survey (Kashino et al. 2013;Silverman et al. 2015). This work has multiple aims as we want to test the consistency of both types of data (photometric and spectroscopic) to derive SFR, stellar mass, and dust attenuation for both the lines and continuum. The conservation of the energy budget by CIGALE allows us to get robust estimates of the amount of dust attenuation. In order to get both Hα and [OIII]λ5007 fluxes with good quality, the redshift range is reduced to 1.40 < z < 1.68. We fit the photometry and Hα flux emission to constrain the attenuation and correct emission line luminosities for dust effects. We focus our study on understand-ing the SFR-[OIII]λ5007 relation, and the influence of metallicity and ionization field. In the FMOS-COSMOS catalog, we have access to other emission lines as Hβ and [NII]λ6584, which are used to emphasize the difficulties in HII-region modeling and coupling these models with SED fitting methods.
The data selection is presented in Sect. 2. The wavelength coverage of the photometric data goes from the UV to the far-IR including emission lines fluxes from the FMOS-COSMOS survey. We fit the Hα emission line flux and photometric data to measure SFR and stellar mass and put our sample in a SFR-M star diagram, the SED fitting analysis, and its results are presented in Sect. 3. Dust attenuation for both stellar continuum and emission lines are discussed in Sect. 4 and we study the SFR-[OIII]λ5007 relation in Sect. 5. In Sect. 6, we add other lines to the fitting analysis and discuss the ability of extended CLOUDY models to fit the full set of data. Finally, in Sect. 7, we summarize the main parts of this work and conclusions. This work assumes a flat seven-year Wilkinson Microwave Anisotropy Probe (WMAP7) cosmology (Komatsu et al. 2011) with h = 0.704, Ω Λ = 0.727, and Ω M = 0.273, AB magnitudes, and a Chabrier (2003) initial mass function (IMF).

Spitzer MIPS and Herschel PACS and SPIRE data
For the IR range, we accessed data from the Herschel Extragalactic Legacy Project (HELP 2 ). We used the Spitzer Multi-Band Imaging Photometer (MIPS) 24 µm, Herschel Photodetector Array Camera and Spectrometer (PACS) 100 and 160 µm, and Herschel Spectral and Photometric Imaging REceiver (SPIRE) 250, 350, 500 µm data products.
Spitzer MIPS and Herschel PACS and SPIRE fluxes are the result of the HELP XID+ extraction based on the position of priors; XID+ uses a Bayesian probabilistic framework that includes priors to measure fluxes and their uncertainties (Hurley et al. 2017).
Both Spitzer MIPS 24 µm and Herschel PACS 100, and 160 µm were obtained using Spitzer IRAC sources from the Spitzer Large Area Survey with the Hyper-Suprime-Cam (SPLASH; Capak et al. 2012) as positional priors. In the case of Herschel SPIRE 250, 350, 500 µm, XID+ was run using 24 µm priors from Le Floc'h et al. (2009) catalog with the addition of the Submillimetre Common-User Fig. 1. Spectro-photometric sample. In blue, we show the position for objects in the HELP-project catalog. In orange, the main bulk of photometry from the COSMOS2015 catalog (Laigle et al. 2016). In black, the FMOS-COSMOS emission line targets. In red, the final selection for this work.

FMOS-COSMOS emission lines
The Fiber Multi-Object Spectrograph (FMOS) is an instrument located at the Subaru telescope at the Mauna Kea Observatory in Hawaii. The instrument is a fiber-fed system that allows for wide-field spectroscopy and enabling NIR spectroscopy. We used data from the FMOS-COSMOS survey as described in Kashino et al. (2013) and Silverman et al. (2015) (Kashino & Silverman) model (Silverman et al. 2015). We only considered emission lines detected with a signal-to-noise ratio (S/N) larger than 3. The measurement error in the emission line corresponds to the formal error provided by the line fitting process. The aperture correction value provided for each object is the best value derived from three different methods (see Kashino et al. 2019) and is used to correct in-fiber measurements. We keep only objects for which the flux loss is lower than 70%. To calculate the full uncertainty (fitting process and aperture correction) the formal errors on the observed emission line must be added up in quadrature, with a suggested factor of 1.5 from Kashino et al. (2019), to account for the uncertainty introduced by the aperture correction.

Final sample selection
We started with a sample of 1182108 galaxies from the COS-MOS2015 catalog, which was reduced to 199 once it was crossmatched to the available IR data from the HELP-project and the FMOS-COSMOS catalog. The different data sets used in this work are represented in Fig. 1 A total of 182 objects were selected on the basis of their flux measurements for both Hα and [OIII]λ5007 with S/N > 3. In Fig. 3, the final sample redshift distribution covering 1.4 < z < 1.68 is presented. This selection is shown in Fig. 1 as red dots over imposed on the FMOS-COSMOS sample (black dots), the COSMOS2015 catalog from Laigle et al. (2016), and the HELPproject data. The number of objects per band is shown in Table 1. S/N drops for Herschel PACS and Herschel SPIRE bands, leading to 43 and 15 objects with S/N > 3 at 100 µm, and 160 µm, and 50, 35, and 0 objects for the Herschel SPIRE bands 250, 350, and 500 µm, respectively. Our entire sample has a S/N > 3 in the Spitzer MIPS 24 µm, Spitzer IRAC1 3.6 µm, and Spitzer IRAC2 4.5 µm bands, but drops for Spitzer IRAC3 5.8 µm and Spitzer IRAC4 8.0 µm with 95 and 44 objects, respectively. Individual poorly measured photometric data does not affect the SED fitting result because the weight of these bands in the overall fit is small due to their large measurement error.
Most of the sample includes UV-to-MIR photometry and Hβ, [NII]λ6584, [SII]λ6717, and [SII]λ6731 emission line fluxes. The quality of the measure of the flux is assessed by the S/N. Of the 139 objects that have data for the Hβ and [NII]λ6584 emission lines, 114 and 112 sources have a S/N > 3, respectively.
[SII]λ6717 and [SII]λ6731 emission line measurements are of very poor quality, leading only to 22 and 15 objects with a S/N > 3. The quality of the Hβ measurements in terms of the Balmer Decrement (BD) will be addressed in Sect. 4. In Sect. 6, we present a more detailed study of the emission lines. We need to avoid any possible contamination due to an AGN contribution. Kashino et al. (2017) flagged AGN in the FMOS-COSMOS sample as point sources with an associated X-ray emission provided by the Chandra-COSMOS Legacy survey catalog (Elvis et al. 2009;Civano et al. 2016) with L X−ray 10 42 erg s −1 at 0.5 − 0.7 keV. Any source flagged as being associated with X-ray emission was discarded from our main sample. Obscured AGN contamination is still possible even if X-ray detected sources are excluded. Kashino et al. (2017) identified 39 objects consistent with obscured (type-II) AGN based on narrow line ratios, with only four of those objects having an X-ray counterpart. Nevertheless, we did not exclude these objects from our analysis because our results are robust enough even if these ob-  jects are included probably because they are not extreme cases. We used the color criteria of Donley et al. (2012) to separate AGN from star-forming galaxies and verify that none of our objects lie inside the AGN region.

SED fitting with CIGALE
In this section, we present the SED fitting process carried out with the code CIGALE . Here, we describe the modules that we use and their direct impact on the analysis of our data set. The structure of the code is addressed by different authors in the literature (e.g., Ciesla et al. 2017;Małek et al. 2018;Buat et al. 2018;Boquien et al. 2019). CIGALE is based on the energy balance principle in which dust partially absorbs emission of all origins in the UV-optical wavelength range and re-emits this energy in the IR. The code creates millions of models to be compared with the observations and estimates physical parameters of galaxies (such as SFR, stellar mass, dust luminosity, dust attenuation, AGN fraction) while applying a Bayesian statistical analysis approach.
Different modules were chosen from the CIGALE library to model the star-formation history (SFH), the nebular emission, the dust attenuation, and re-emission. The parameters and their uncertainties are computed using the goodness of fit for all the models as a likelihood-weighted mean and likelihood-weighted standard deviation, respectively. A global indicator of the quality of the fits is given by the reduced χ 2 (χ 2 r ). In general, the χ 2 r must account for the degrees of freedom, however, the non-linearity of the equations (Chevallard & Charlot 2016) linking the parameters and their non-trivial dependence on each other makes it difficult to include them. Therefore, the χ 2 r corresponds to the χ 2 divided by the total number of input fluxes.
The combination of Hα fluxes with UV-to-IR data with the preservation of the energy budget allows us to get reliable estimates of dust attenuation (see Sect. 3.3) without resorting to other methods such as the Balmer Decrement (BD). It also allows for secure estimations of the SFR to be obtained. Although we have Hβ, [NII]λ6584, [OIII]λ5007, [SII]λ6717, and [SII]λ6731 information, we did not fit these lines due to current discrepancies concerning the HII-region models (see Sect. 6) and to the poor quality of the data, as presented in Sect. 2.4 and Table 1. In particular, the Hβ fluxes were discarded either because of their low S/N or a Hα/Hβ ratio below the canonical value of 2.86 (see Sect. 4).
[NII]λ6584, [SII]λ6717, and [SII]λ6731 only satisfy the S/N criterion for 112, and 22, and 15 objects, respectively. The relevant CIGALE modules selected for this work are introduced in the next sections.

Star-formation history
We used a delayed SFH with the functional form presented in Eq. 1. This form depends on the time of the star-formation onset, t 0 , and the e-folding time of the stellar population, τ main . A recent burst of constant star formation that has been going for at most 70 Myr is superimposed to the delayed SFH. This form allows us to have a variation of the SFH where the SFR increases from the onset of star-formation until its peak at τ main . After that point, the SFR starts declining: (1) The delayed component allows us to add a burst that we define with constant amplitude over the last years (Corre et al. 2018;Carnall et al. 2019;Leja et al. 2019;Chevallard et al. 2019;Buat et al. 2019). The amplitude of the burst is fixed by f burst , defined as a mass ratio between the stars formed during the burst and the total mass of stars. As we select objects with emission lines, we expect the amplitude of the burst to have a direct impact on the derivation of SFR for at least some objects in our sample. The SFR is averaged over 10 Myr compatible to the timescale traced by Hα. The input parameters are described in Table 2, where the e-folding time and the age of the main stellar population are left as free parameters as well as the age and the mass fraction of the burst. We used a Chabrier (2003) IMF at solar metallicity 0.02 for the stellar population. Effects on the variation of the estimated parameters due to a fixed stellar metallicity are discussed in Sect. 3.4.

Nebular emission lines
CIGALE models the galaxy's emission of the ionized gas by stellar generations as an effective HII-region, encompassing the ensemble of HII-regions and the diffuse ionized gas. The nebular emission lines are pre-computed in the nebular module of CIGALE and re-scaled with the number of Lyman continuum photons from the stellar emission of the model galaxy.
The radiation field intensity is given by the dimensionless ionization parameter log U ≡ log(n γ /n H ), where n γ is the number density of photons capable of ionizing hydrogen and n H the number density of hydrogen. The photo-ionizing field is generated with the single stellar population (SSP) model library of Bruzual & Charlot (2003), using a constant SFH over 10Myr. In Sect. 6, we detail the major changes made in modeling nebular emission as compared to the previous version described in Boquien et al. (2019). The updated version of the HII-region models is used alongside this work to fit photometry and emission line fluxes. Our main purpose is to measure dust attenuation and SFR and to compare the latter to [OIII]λ5007 luminosities. The new photoionization models allow us to interpret the locus of our galaxies in the excitation diagrams in an attempt to understand the underlying physics of the ISM.

Dust attenuation recipe and dust emission
We adopted the recipe proposed by Charlot & Fall (2000) (CIGALE module called dustatt_modified_CF00; hereafter CF00), where two stellar populations are considered: young stars (age < 10 7 years) that are still located in a birth cloud (BC), while older stars (age > 10 7 years) have already moved into the interstellar medium (ISM). Both populations experience a different dust attenuation: the emission of young stars goes through the BC and the ISM while the emission of older stars is only attenuated by dust located in the ISM. The dust attenuation in the ISM and the BC are both modeled as power laws normalized to the amount of attenuation in the V-band (λ V = 0.55µm): The ratio of the attenuation in the V-band experimented by the young and old stars is given by: where µ is considered as a free parameter in our fitting process (da Cunha et al. 2008;Battisti et al. 2019). Nebular lines coming from the Lyman continuum photons due to very young stars are attenuated as young stars (BC+ISM). The range for each parameter used in the CIGALE's attenuation module is given in Table  2.
The slopes of the two power-laws (i.e., BC and ISM) are fixed to −0.7 following Charlot & Fall (2000). This value reproduces the observed mean relation between the IRX and UV spectral slope of nearby starburst galaxies. Lo However, it has been shown from HII-region studies (see Caplan & Deharveng 1986;Liu et al. 2013) that grayer values closer to the one chosen in this work are more suitable for reproducing the effective attenuation in dusty galaxies. Letting n BC free is not suitable because we use only one emission line and its value would be poorly constrained. In Sect. 4, we explore the effects of changing −0.7 to −1.3 for n BC . The only free parameters in our recipe are µ and A ISM v . The µ parameter relates the undergone attenuation in the V-band by the ISM (i.e., old stars, continuum) and the BC+ISM (i.e., young stars, emission lines). Leaving µ free to vary allows for some variation on the effective attenuation law (Battisti et al. 2016;Buat et al. 2018;Małek et al. 2018;Chevallard et al. 2019). Moreover, introducing this flexibility in the SED fitting process allows for a better quality fit of the Hα emission. Variations of the attenuation law will be explored as further work. The dust emission was fitted with the Draine & Li (2007) models based on a set of parameters to constrain the starlight intensity and link dust to star-formation including updates of Draine et al. (2014). We used these models because we have 24 µm information making them better suited and more flexible for our purpose.

Spectral energy distribution fitting results
We fit the SEDs of our sample using a set of parameters obtained from previous works and summarized in Table 2. We analyzed Article number, page 5 of 17 A&A proofs: manuscript no. aanda Table 2. CIGALE modules and input parameters used for the SED fitting process as presented in Boquien et al. (2019) Bruzual & Charlot (2003) initial mass function IMF Chabrier metallicity Z star 0.02 Nebular Emission ionization parameter logU -4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0 gas metallicity Z gas 0.008, 0.006, 0.01, 0.02, 0.03, 0.04, 0.05 Dust attenuation Templates based on values adopted by Charlot & Fall (2000); Buat et al. (2018) ISM attenuation in the V-band (mag) A ISM V 0.0, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9 the quality of the fits, using the global χ 2 and also comparing the Bayesian estimate of each flux with the observed values. In Fig.  4, the distribution of the obtained χ 2 and χ 2 r are shown. For each filter, we computed the median value of the difference between the Bayesian estimate from CIGALE and the observed values of fluxes with S/N > 3. Their distribution is presented in Fig. 5. The SED-modeled Hα fluxes are presented in Fig. 6 as a function of the observed flux as well as the flux difference (Hα CIGALE -Hα FMOS ).
From Fig. 5, it is clear that the fluxes are well fitted by CIGALE with a flux difference lower than 0.1 dex and dispersion lower than 0.13 dex except for GALEX NUV for which the dispersion reaches 0.28 dex. Herschel PACS 100 µm and 160 µm are underestimated and exhibit a larger dispersion of 0.37 dex and 0.31 dex, respectively. ONLY 43 objects have S/N > 3 in Herschel PACS 100 µm. 23 objects have 100 µm fluxes larger than their respective flux at 160 µm. This is an unrealistic configuration at z ∼ 1.5 with rest-frame fluxes 40 and 60 µm, respectively. The dispersion of the results in Herschel PACS bands can be reduced (e.g., red squares in Fig. 5) if we consider the sample of 9 objects with S/N > 3 in both Herschel PACS 100 µm and Herschel PACS 160 µm. Therefore, the overall bad results at Herschel PACS bands can be explained as a combined effect of low S/N and non-realistic configuration at the IR bands. We used a deblended IR catalog in the COSMOS field provided by (Jin et al. 2018) to check if the situation for these bands could be improved. Both data sets (i.e., HELP and Jin et al. 2018) are consistent for our sample with the same trends in terms of the unrealistic configuration of Herschel PACS 100 µm and Herschel PACS 160 µm. The Jin et al. (2018) data have even lower S/N than our final sample selection (described in Sect. 2.4). Using HELP data for our sample leads to better results for the fits. CIGALE weighs the data as the inverse square of the error, which means that even if flux with a low S/N is badly fitted, its impact on the overall SED fitting output is minimal. We also fit our sources without the Herschel PACS and SPIRE bands to check whether SFR, stellar mass, and attenuation estimations are affected when IR data are excluded. The retrieved estimates with and without IR information are consistent with each other within 1σ uncertainty. The impact of Herschel SPIRE data on the SFR as compared to the Herschel PACS bands is larger (i.e., 0.15 dex). In both cases, the difference in attenuation is larger (3×) for [OIII]λ5007 than for Hα. We kept the Herschel data for the whole sample in the overall fit. In Fig. 6, the Bayesian values of the Hα fluxes and the FMOS measurements follow the 1:1 relation, indicating that we are able to fit the line flux very well within a 0.2 dex scatter, as already found by Buat et al. (2018) on a smaller sample in the same field using spectroscopic information from the 3D-HST survey.
Accurate determination of SFRs is crucial as we want to compare them with the [OIII]λ5007 luminosities. Stellar masses and SFRs retrieved from the SED fitting cover a range of 10 9.5 − 10 11.5 M , and 10 1 − 10 3 M yr −1 , respectively. Current high SFRs can be expected for some objects since our sample is selected on the detection of recombination lines. To confirm the need to add a burst the delayed SFH we performed a fit with a burst amplitude set to 0. We compared both fits by calculating the Bayesian Information Criterion (BIC) of each of them as  introduced by Schwarz (1978) and defined as: where N stands for the number of observations namely, the number of bands fitted, κ the number of degrees of freedom, and χ 2 = −2 × ln(likelihood) the non-reduced χ 2 obtained from CIGALE (Ciesla et al. 2017;Buat et al. 2018). A significant difference and evidence against a model is characterized by the higher BIC. We computed the difference of BIC values ∆(BIC). As the number of bands fitted is the same in both cases, any difference arises from the number of degrees of freedom (κ fix = 14 and κ free = 16) and the χ 2 . ∆(BIC) > 2 for 90 (∆(BIC) > 6 for 69) objects, meaning that for 1/3 to 1/2 of our sample the introduction of a burst improves the fit. In Fig. 7, our sample is presented in the SFR-mass diagram color-coded with f burst . The Schreiber et al. (2015) main-sequence relation for a redshift of z = 1.5 is drawn with a 0.3 dex dispersion (Karim et al. 2011;Rodighiero et al. 2011). The dispersion computed by error propagation from the originally fitted parameter errors is also shown as a darker shaded contour. The bulk of galaxies is found to lie within the 0.3 dex curves. Galaxies four times above the main-sequence are usually considered starbursting (Sargent et al. 2012). A total of 31 objects in our sample are consistent to be starbursts with f burst > 0.06 within a 1σ dispersion. The median value of f burst for these objects is 0.10±0.03. We fit the data using a stellar continuum metallicity of 0.008 and confirm that the re-sults are consistent with the 0.02 case. Variations are within 1σ dispersion for the stellar mass, SFR, and attenuation estimates.

Dust attenuation
As noted above, we simultaneously fit the photometric and spectroscopic data (i.e., Hα fluxes) using the Charlot & Fall (2000) modified recipe. Emission lines are produced by excited gas inside the BC as a result of the emitted radiation by young stars. Therefore, the lines suffer from attenuation due to the BC and by the ISM surrounding the HII-region.
The combination of the energy budget principle of stellar photons being absorbed and re-emitted by dust in the Mid-andfar-IR with the attenuation recipe allows us to measure the net effect of dust obscuration (i.e., the effective attenuation) in our galaxy sample. CIGALE calculates the resulting attenuation at any wavelength and for any emission line allowing us to obtain the amount of attenuation in the Hα and [OIII]λ5007 emission lines directly from the fits using photometry and Hα emission line fluxes. A comparison to the attenuations estimated with the BD method is addressed in this section.
The reliability of the estimation of the effective attenuation has to be checked carefully. We used a CIGALE option to create and analyze mock catalogs and to verify if the main parameters involved in our derivation can be trusted and ensure that our analysis is robust enough to constrain the attenuation. Each flux in the mock catalog corresponds to the value of the best model and is modified by adding a value taken from a Gaussian distribution centered at 0 and with the same standard deviation as the observations. The mock catalog of fluxes contains the same number of sources and is fitted in the same way as the original catalog. The estimated and exact values for the two free parameters in the Charlot & Fall (2000) recipe and the Hα attenuation are compared in Fig. 8. The value for A ISM V is very well reproduced with a Spearman's correlation coefficient of ρ s = 0.95. The µ parameter is more difficult to constrain even if a clear positive correlation is found between exact and estimated values (Spearman's correlation coefficient of ρ s = 0.84). The ratio of the attenuation in the V-band experienced by the young and old stars characterized by µ is difficult to constrain and as a consequence A BC V is also not very well constrained. Moreover, A Hα is well-constrained with a Spearman's correlation coefficient of ρ s = 0.93, reflecting the effects of the BC+ISM combination.
Attenuation robustness: Mocks, Balmer decrement, and stellar mass trends The retrieved attenuation for the Hα and [OIII]λ5007 emission lines is a crucial output of the SED fitting process for this work. We can use these values to correct line luminosities as the attenuation is constrained by the energy budget method through the UV-to-IR photometry coverage. Nevertheless, the measure of attenuation based on the BD method is widely used to compute the amount of attenuation when both Hα and Hβ are available. Attenuation in the Hα line can be obtained directly using the observed Hα and Hβ fluxes and assuming a given extinction curve following: A Hα = −2.5κ Hα κ Hβ − κ Hα log 10 2.86 Hα/Hβ , where κ Hα and κ Hβ represent a particular extinction curve evaluated at Hα and Hβ wavelengths, respectively. The attenuation for the hydrogen lines (i.e., Hα) can be translated into an In order to compare BD derived attenuations with CIGALE's attenuation estimates, we selected, from a total of 139 objects with both Hα and Hβ measurements, a sub-sample of 80 objects satisfying the Hα/Hβ > 2.86 6 criterion. Here, 23 of the excluded sources have a Balmer decrement that is consistent with 2.86 within 1σ uncertainty. The rest that are not consistent have overestimated Hβ fluxes as compared to Hα. We did not include any of these sources in the analysis so that we could remain consistent with the photoionization models implemented in CIGALE. Excluding values below the canonical value of the ratio allows us to exclude negative attenuation values inconsistent with the case B recombination assumption for the emission lines. We computed the attenuation applying Eq. 5, along with the Milky Way curve (κ Hα = 2.52; and κ Hβ = 3.66; κ [OIII] = 3.52) of Cardelli et al. (1989) updated by O'Donnell (1994. For comparison and to show the effects in the derived attenuation due to the choice of an extinction curve, we also considered the Calzetti et al. (2000) (C00) curve (κ Hα = 3.33; and κ Hβ = 4.59; κ [OIII] = 4.46) 7 .
In Table 3, we present median attenuation values for Hα and [OIII]λ5007, as well as the A [OIII] -A H α ratio. All values are derived with a fixed BC slope n BC = −0.7 except for the starred case in which n BC = −1.3 was implemented. The fit with n BC = −1.3 is as good as the n BC = −0.7 case. In Fig. 9, the derived attenuations using the BD method are compared to CIGALE's estimates. We check the distribution of CIGALE's attenuation in the [OIII]λ5007 emission line as a function of stellar mass and fit a linear relation shown in orange in the lower panel of Fig. 10. Thus, we obtain: A [OIII] = (1.19 ± 0.15) + (0.98 ± 0.22) log 10 (M * /10 10 M ).
We go on to explore the effects of using this relation to correct [OIII]λ5007 luminosities in Sect. 5.3.
Applying either the Milky Way or C00 reddening curve leads to lower median attenuations than those derived using Charlot & Fall (2000) recipe in CIGALE, as shown in Fig. 10. CIGALE's attenuation estimates are globally in agreement with the A Hαmass relation of Garn & Best (2010a). BD derived attenuations 6 The intrinsic value of 2.86 is assumed to be consistent with a temperature T= 10 4 K and electron density n e = 100cm −3 for case B recombination Osterbrock (1989) and does not vary significantly with the physical conditions. 7 This curve is only valid for stellar continuum and not appropriate to describe nebular emission line attenuation (Pannella et    (i.e., blue and red dots) are not entirely consistent with this relation placing the majority of the objects below the relation. Objects above Garn & Best (2010a) relation have a SFR > 10 1.5 M yr −1 but they are not classified as starbursting from our analysis. The Hα attenuation increases with stellar mass in agreement with previous works (Garn & Best 2010b,a;Ibar et al. 2013;Zahid et al. 2017;Shivaei et al. 2020), following the Garn & Best (2010a) relation and confirming that in more massive galaxies nebular regions have a higher attenuation (Koyama et al. 2019). This also confirms that this relation is not just valid for the local universe but at higher redshift as well. From the attenuation ratio between [OIII]λ5007 and Hα presented in Table 3, we can see that our results using CF00 lead to a flatter effective attenuation curve (Chevallard & Charlot 2016;Lo Faro et al. 2017;Buat et al. 2018), as compared to the Milky Way or C00 recipe. The difference arises from the choice in fixing both slopes in the CF00 to the same value of −0.7. We performed the same analysis using n BC = −1.3 and verified that the effective curve is steeper than in the n BC = −0.7 case and closer (and even slightly steeper) than the Milky Way and C00 case.
The attenuation distributions obtained for the entire sample from the SED fitting are shown in Fig. 11 for the Hα and [OIII]λ5007 emission lines as well as for the V-band. The nebular line attenuation distribution in both cases follows a similar behavior as both lines are produced inside the HII-regions and they are attenuated in a similar way (inherent to CIGALE modeling). Stellar continuum attenuation as traced by A v is lower. The differential attenuation suffered by young and old stars is measured by the µ parameter in the CF00 recipe as explained in Sect. 3.3. The original value proposed by Charlot & Fall (2000) corresponds to 0.3 for the nearby universe. We measured a larger Article number, page 9 of 17 A&A proofs: manuscript no. aanda value of µ = 0.57 ± 0.14. In general, higher values are found at higher redshift (see Buat et al. 2018, and references therein). This is in agreement with nebular emission being more attenuated than the stellar continuum.

[OIII]λ5007, Hα fluxes, and SFR measurements
In Sect. 4, we show that we constrain well the attenuation in our sample with the multi-wavelength plus IR bands coverage and the inclusion of the Hα emission line, and the CF00 recipe, which introduces a differential attenuation between young and old stellar populations.
We computed the [OIII]λ5007 line luminosities for our sample of galaxies using FMOS-COSMOS observed fluxes. These luminosities are corrected for dust effects using the Bayesian attenuation A [OIII] presented in Sect. 4. We find that the Hα and [OIII]λ5007 emission lines span a similar range in luminosities and the amount of attenuation is found to increase with luminosity and SFR confirming previous works (e.g., Cortese et al. 2012;Bourne et al. 2012;Santini et al. 2014). In Fig. 12, the observed and dust corrected Hα and [OIII]λ5007 luminosities are shown in gray and color dots, respectively. Emission lines are not strongly correlated before dust correction. After accounting for dust attenuation both quantities correlate (with a Spearman's coefficient of ρ s = 0.6). In the next subsections, we address the gas-phase metallicity influence on the [OIII]λ5007/Hα line ratio, the SFR measurements, and we compare both [OIII]λ5007 and [OIII] 88 µm line emissions from the photoionization models.

Gas-phase metallicity
The [OIII]λ5007 emission line is sensitive to the ionization parameter (i.e., the ionizing field; Dopita et al. 2013Dopita et al. , 2016Nicholls et al. 2017) but can be also affected by the gas-phase metallicity (Kennicutt 1992;Kennicutt et al. 2000;Steidel et al. 2014;Gutkin et al. 2016;Byler et al. 2017). Different mass metallicity relations (MZR) based on direct and strong-line methods are reported in literature, sometimes including the effects of the SFR (e.g., Pettini & Pagel 2004;Mannucci et al. 2010;Andrews & Martini 2013;Zahid et al. 2017;Sanders et al. 2021Sanders et al. , 2020Curti et al. Fig. 12. Hα and [OIII]λ5007 luminosity. The observed luminosities uncorrected for dust are shown as gray dots. The color dots correspond to luminosities de-reddened using A Hα and A [OIII] as constrained by the SED fitting with CIGALE. The corrected data is color-coded by the stellar mass. A slope of 0.99 is measured in the dust corrected sample with a 0.39 dex dispersion. 2020). We compute oxygen abundances using the fundamental metallicity relation (FMR) calibration of Curti et al. (2020), which is based on strong-line oxygen abundance measurements for Sloan Digital Sky Survey (SDSS) galaxies. Their relation is fully consistent with a relation derived for a MOSFIRE Deep Evolution Field (MOSDEF) sample of individual and stacked galaxies at z ∼ 2.3 in Sanders et al. (2021). The dust-corrected [OIII]λ5007/Hα ratio is presented in Fig. 13 as a function of oxygen abundance. Our galaxy sample covers a sub-solar abundance range of 8.4 < 12 + log(O/H) < 8.8 (i.e., gas-phase metallicity 0.006 < Z gas < 0.016) 8 . The ratio slightly decreases when metallicity increases. The median value of the ratio is 0.86 (or −0.065 dex) with a dispersion of 0.3 dex as measured by the standard deviation of the data set. The dispersion at a fixed metallicity can be produced by the sensitivity of the oxygen species to the ionizing field. The dependence on both metallicity and ionizing field is likely to produce effects on the ratio which we cannot easily disentangle. In Sect. 6.2, we explore variations of the intensity of the radiation field fixing the gas-phase metallicity in three bins with an equal number of galaxies and simultaneously fitting the Hα, Hβ, and [OIII]λ5007 line fluxes.

SFR−[OIII]λ5007 calibration
In Fig. 14, we present the dust-corrected [OIII]λ5007 observed luminosities, and SED fitting derived SFR. As we want to make a comparison with the relations in the literature based on a fixed [OIII]λ5007/Hα ratio, we first fit only the intercept using an orthogonal distance regression (ODR) in the L [OIII] -SFR plane. The quality of the fit is measured by the non-parametric Spearman's and Kendall's correlation coefficients equal to ρ s = 0.57 and ρ s = 0.40, respectively. This leads to:  rical distance of each data point to the line drawn by the equation above. We check the geometrical distance distribution and trim 15 data points with d ≥ ±0.45 that corresponds to ∼ 2σ. These outliers are characterized by poor S/N in GALEX NUV, and the Herschel bands. We check that the distribution in the SFR-L [OIII]λ5007 plane after discarding these points remains representative of our sample. We implemented a bootstrap method that uses a random iterative replacement technique to estimate the slope and intercept in the fits. Using this method, we drew a new sample of the same size 200 times and re-fit using the ODR method with a fixed slope to obtain the distribution of the intercept from which we compute the mean value and the uncertainty in the intercept as the standard deviation of the distribution. This leads to: log 10 SFR = log 10 L [OIII] − (41.20 ± 0.02), with the luminosity in ergs s −1 and SFR in M yr −1 . The relation is consistent within a 2σ uncertainty with the equation obtained including the outliers. Following the same procedure described earlier in this paper and allowing the slope and intercept to vary, we obtain: log 10 SFR = (0.83 ± 0.06) log 10 L [OIII] − (34.01 ± 2.63), with the luminosity and SFR given in units similar to Eq. 6. The slope in Eq. 7 is consistent with unity at 3σ. The relations presented in Eq. 6 and 7 are drawn in Fig. 14 as a black dashed and continuous line with a 0.32 dex dispersion from the standard deviation of the observations as a gray shaded area. Since the Hα emission can be considered as proportional to the SFR, the dispersion of the [OIII]λ5007-SFR correlation is expected to have the same origin as the one affecting the [OIII]λ5007/Hα ratio: the metallicity and radiation field.
Several authors in the past studied the link between the [OIII]λ5007 emission line and SFR based on a fixed [OIII]λ5007/Hα ratio and the calibration of the Hα emission line in terms of SFR. At z ∼ 0.1 on a sample of 196 SDSS narrow band emitters Moustakas et al. (2006) measured a scatter in the [OIII]λ5007-SFR relation as large as a factor of 3 − 4 uncertainty (1σ) when compared to other tracers such as [OII]λ3727, Hβ emission lines or the U-band. Based on a sample of 92 galaxies at 0.40 < z < 0.64, Hippelein et al. (2003) derived an [OIII]λ5007/Hα ratio of 0.79 from line luminosity functions, which is consistent with our median value. Although [OIII]λ5007 depends on excitation and metallicity, their estimated median values of SFR([OIII]) were consistent with the evolution of SFR density of the universe. Ly et al. (2007) reported a [OIII]λ5007/Hα ratio of 1.05 at z = 0.07 − 1.47 for a spectroscopic sample of 196 narrow-band emitters. Teplitz et al. (2000) fixed the ratio to unity at z > 3 for a sample of five galaxies observed with the near-IR camera on the Keck I telescope. Both works from Teplitz et al. (2000) and Ly et al. (2007) found similar results on the dependence of the line ratio on metallicity and ionization field. Maschietto et al. (2008) studied a sample of 13 [OIII] emitters at z ∼ 3.09 − 3.16 and using a [OIII]λ5007/Hα ratio of 2.4 derived a lower limit relation for SFR and [OIII]λ5007 luminosity. Straughn et al. (2009) also explored the relation based on [OIII]λ5007/Hα ratios in knots of 136 galaxies at z ∼ 0.5 for which both emission lines were observed. Their relation implies a median [OIII]λ5007/Hα ratio of ∼ 1.23. Bowman et al. (2019) found a strong correlation but did not propose any relation for a set of 3D-HST galaxies at z ∼ 2 (see their Figure 6).
Our calibration relies on quantities estimated for each galaxy (e.g., [OIII]λ5007 luminosity, SFR). It can be compared to relations also obtained with individual [OIII]λ5007/Hα ratios as those obtained by Hippelein et al. (2003), and Ly et al. (2007) 9 . Their dust corrected relations are reported in Fig. 14 (orange and green lines) and are fully consistent with our result within the 1σ dispersion error and our median reported [OIII]λ5007/Hα ratio of 0.86.

Impact of the [OIII]λ5007 dust attenuation
In the previous section, we studied the relation between the [OIII]λ5007 dust corrected luminosity and SFR. In practice, it is very difficult to correct the [OIII]λ5007 emission line for dust attenuation. When A [OIII] cannot be estimated but stellar masses are known, attenuations can be tentatively estimated by using the A [OIII] -M star relation presented in Sect. 4. We applied this relation to correct our observed [OIII]λ5007 luminosities and checked how the relation presented in Eq. 6 is modified. Following the same procedure presented in Sect. 5.2, we find: log 10 SFR = log 10 L [OIII] − (41.23 ± 0.03).
The average relation is similar to that of Eq. 6, but a large dispersion is expected because of the dispersion in the A [OIII] -M star relation. To test it, we divided our sample in mass bins of 0.5 dex from 10 9.5 to 10 11.5 M and applied both CIGALE's estimated attenuation and the mass-method attenuation. We measured an average increase in SFR of only 0.15 dex when we apply Eq. 8, as compared to the SFR estimated by CIGALE. The dispersion of 0.3 dex in each bin is much larger than the increase in SFR. Thus, Eq. 8 can be applied when only an averaged value of SFR is needed when analyzing galaxies whose stellar masses are known. Fig. 14. [OIII]λ5007 corrected luminosity and SFR relation. The dashed and continuous black lines correspond to fits using a bootstrapped orthogonal distance regression method with a Spearman's regression coefficient ρ s = 0.57. The 0.32 dex dispersion from the standard deviation is presented as a gray shaded area. The relation shows a positive correlation with scatter significantly higher compared to Kennicutt (1998) Hα relation. Relations proposed by Hippelein et al. (2003), and Ly et al. (2007) are shown in orange and green, respectively. Straughn et al. (2009) is also presented in red for comparison. Relations are converted from Salpeter to Chabrier IMF.

Benchmark for the [OIII]λ5007 and [OIII] 88 µm lines
The [OIII] also produces a fine-structure emission line at 88.33 µm from low excitation and highly ionized gas located in a low-density environment. This line is not affected by dust attenuation and is also used to trace SFR at both low and high redshifts (e.g., De Looze et al. 2014; Harikane et al. 2020). In this section, we present our efforts to predict a relation between [OIII] 88 µm and SFR from the relation we measured between [OIII]λ5007 and SFR and HII-region models. The [OIII]λ5007-[OIII] 88 µm ratio is known to vary with gas-phase metallicity and density (Dinerstein 1983;Stacey et al. 2010;Ferkinhoff et al. 2010;Moriwaki et al. 2018), but it is not considered to be especially sensitive to the ionization parameter (Moriwaki et al. 2018). To select the photoionization models that are relevant to our analysis, we restricted the ionization parameter to −3.0 < log U < −2.0, in agreement with our results (presented in Sect. 6.2). For the electron density, Kashino et al. (2017) derived a value of ∼ 200 cm −3 for the FMOS-COSMOS sample. Three different values of the electron density are available in CIGALE (i.e., n e = 10, 100, and 1000 cm −3 ) and we fixed n e to 100 cm −3 . From our estimations of the gas-phase metallicity (Sect. 5.1), we restricted the models to 0.006 < Z gas < 0.016 (or 8.4 < 12 + log(O/H) < 8.8). With the selected models, we derived a mean [OIII]λ5007-[OIII] 88 µm ratio of ∼ 1.90 (with a 1σ dispersion of 0.84). We used this average value of the ratio to translate our SFR-[OIII]λ5007 linear relation (presented in Eq. 6) into another linear relation between [OIII] 88 µm and SFR. This is presented in Fig. 15 as a black dotted line with a gray shaded area accounting for the dispersion with the gasphase metallicity.
Our relation is found to be consistent with the relations proposed by De Looze et al. (2014) for metal-poor local dwarf galaxies and by Harikane et al. (2020) for high redshift galaxies despite the different gas-phase metallicities of their samples. Arata et al. (2020) report a L [OIII]88µm − SFR relation derived  Arata et al. (2020) both at z = 6 − 9. We translate our [OIII]λ5007 luminosities into [OIII] 88 µm using a mean ratio of 1.9 derived from CLOUDY HII-region models at an electron density of n e = 100 cm −3 shown as a black dotted line. The gray area is the dispersion in the translation with gas-phase metallicity. from simulations at z = 6 − 9 covering a gas-phase metallicity range of 6.6 < 12 + log(O/H) < 8.9 that matches our metallicity range and is also consistent with different predictions from the literature (De Looze et al. 2014;Moriwaki et al. 2018;Harikane et al. 2020). Both electron density and gas-phase metallicity are paramount to properly convert one emission into the other as well as to compare different samples of galaxies. At intermediate and high densities [OIII]λ5007 line remains an important coolant for ionized gas as compared to [OIII] 88 µm. The [OIII]λ5007 emission line will be targeted in future surveys like MOONS, PFS, and JWST allowing for synergies with observatories like ALMA that can perform observations of [OIII] 88 µm at high-redshift providing a direct method to characterize the overall cooling budget and the ISM physics of HII-regions in intermediate-and high-redshift galaxies as well as to refine their calibration in SFR.

Spectro-photometric modeling: HII-region and stellar continuum modeling
Our 2). In our initial SED fitting analysis, we only included the Hα emission line because of the current challenges in the production of the photoionization models to reproduce observations. We discuss these issues in this section and introduce more line intensities to the fits. The nebular emission lines module was briefly introduced in Sect. 3.2; here, we add a more detailed description of the photoionization calculations. We use the CLOUDY 17.01 photoionization code (Ferland et al. 2017) to predict the relative intensities of 124 atomic lines originating from the far-ultraviolet (FUV) to FIR. Line emissivity models for a plane-parallel geometry are built with the intensity of the photo-ionizing radiation field varying from −4.0 < log U < −1.0 and using a sampling of 0.1 dex. The electron density is set to three different values of n e = 10, 100, 1000 cm −3 covering 25 gas-phase metallicities spanning over 0.0001 < Z gas < 0.05 disconnected from stellar metallicity (0.0001, 0.0004, 0.004, 0.008, 0.02, 0.05) as presented in Table 4. The closest stellar metallicity to the gas-phase metallicity is chosen as the input ionizing spectrum in the models. The cosmic abundance standard and the derivation of element abundances implemented in our models are discussed in Nicholls et al. (2017). They are based on a scaling developed by Nieva & Przybilla (2012) derived from observed metallicities of 29 early B-type stars in the Milky Way. The depletion onto grains is not taken into account since the abundances are based on gas-phase measurements. Our models are built on the so-called local Galactic concordance, 12 + log(O/H) GC = 8.76, which is close to the 8.73 estimated primordial solar abundance derived by Asplund et al. (2009) andLodders (2010), with .46, and Z GC =0.0142. Super-solar and sub-solar metallicities are scaled following the latter definitions. As compared to the previous version implemented in CIGALE, these models disconnect the stellar and gas-phase metallicities and implement current abundance measurements, making them more appropriate for reproducing the line fluxes of our sample of galaxies. We used these new models to fit both photometry and other emission lines similar to Sect. 5, where photometry and only Hα emission fluxes were fitted together.  (Baldwin et al. 1981); see also Veilleux & Osterbrock (1987)) and the [SII]λλ6717,31 excitation diagrams (hereafter [NII]-BPT and [SII]-BPT) along with our photoionization models are presented color-coded by the ionization parameter and the gas-phase metallicity. In the left panel, the demarcation between star-forming galaxies and systems hosting an AGN of Kauffmann et al. (2003) from SDSS data at z = 0 is shown as a black thick line and the extreme starburst separation line of  at z = 1.6 as a dashed line. The blue, green, and orange lines correspond to best-fit relations for the loci of galaxies in the [NII]-BPT at z ∼ 2.2, z ∼ 2.3, and z ∼ 1.6 from Strom et al. (2017), Shapley et al. (2015), and Kashino et al. (2017), respectively. Individual FMOS-COSMOS points show the well-known offset from the local sequence. Our FMOS-COSMOS sample with valid measurements for the four emission lines and their corresponding measured errors are represented by dots. For a smaller sample of 16 objects, we show, in the right panel, the [SII]-BPT diagram with the star-forming and AGN separation of Kewley et al. (2001) and the best fit for the loci of galaxies of Strom et al. (2017).
The CLOUDY models cover well the star-forming region below the Kauffmann et al. (2003) line in the [NII]-BPT diagram. However, only 61% of the flux ratios presented are covered pointing toward a difficulty related to the nitrogen-to-oxygen abundances more than the photoionizing field. In this diagram, 23% of the objects lie in the composite region where AGN and star-forming galaxies coexist at z = 0. For the [SII]-BPT diagram models are in agreement with the locus of galaxies with valid measurements within the observational errors. The HIIregion models able to predict emission line ratios above the local star-forming relation are difficult to create because of i. nitrogen abundance underestimation; ii. ionizing field hardness choice; iii. gas-phase metallicity and density discrepancies or single stellar population models. The models depend strongly on the relative abundances of nitrogen and oxygen, and thus on the choice of the standard metallicity scale (e.g., Nicholls et al. 2017). The nebular scaling of nitrogen with oxygen is problematic: while oxygen is principally produced in core-collapse supernovae in the native gas cloud from which the HII-region formed, nitrogen has both primary and secondary abundances, which are caused by delayed nucleosynthesis through hot-bottom burning and dredge-up in intermediate-mass stars as they evolve (Vila-Costas & Edmunds 1993). Increasing the nitrogen abundance has been suggested by many authors as a way to match the observations (e.g., Masters et al. 2014Masters et al. , 2016Steidel et al. 2014;Shapley et al. 2015;Yabe et al. 2015;Cowie et al. 2016;Sanders et al. 2016). A shift by 0.2 − 0.4 dex in the N/O fraction is proposed by Masters et al. (2016) while Kojima et al. (2017) and Strom et al. (2017) require ∼ 0.1 dex to cover the locus of galaxies in the BPT at z ∼ 2.3. We find that in order to cover all the sample of galaxies at z ∼ 1.6, a shift of 0.2 − 0.3 dex is necessary.
[OIII]λ5007 and [NII]λ6584 have an ionization potential (IP) of 35.12 and 14.53 eV and are ionized by both charge transfer with protons and photoionization. Thus, [OIII]λ5007 and [NII]λ6584 abundances are less sensitive to the photoionizing field, which are determined indirectly through the HII abundance, and the initial oxygen and nitrogen abundances.
[SII]λλ6717, 37 has an IP of 10.36 eV, making its abundance very sensitive to the hardness of the photoionizing field in the UV. Indeed, Steidel et al. (2014Steidel et al. ( , 2016 applied harder radiation fields using the binary Population and Spectral Synthesis code (BPASS) and these authors were able to produce models covering larger values of log([OIII]λ5007/Hβ), but always below the Kauffmann et al. (2003) line. Other authors point to a higher ionization parameter or electron density as the main parameter producing the offset (e.g., Brinchmann et al. 2008;Dopita et al. 2016;Kojima et al. 2017;Bian et al. 2020). Nevertheless, current stellar synthesis models are unable to produce such hard radiation fields (Levesque et al. 2010). The densities from models in high-redshift galaxies (z ≥ 1 − 2) are higher than in low-redshift galaxies, on the order of several 10 2 to several 10 3 per cm −3 (Shimakawa et al. 2015;Kashino et al. 2017). Density (i.e., the pressure) has only secondary effects on the locus of the models in the [NII]-BPT diagram (Masters et al. 2016).

Hα Hβ and [OIII]λ007 SED fitting
To test the performance of SED fitting including more emission lines, we fit simultaneously the photometry and the Hα, Hβ, and [OIII]λ5007 fluxes to explore the impact of the ionization parameter in the dispersion of the proposed SFR-L [OIII] relation. Our sample of galaxies spans over a gas-phase metallicity range of 0.006 < Z gas < 0.016 (see Sect. 5). We di-  2015), Kashino et al. (2017), and Strom et al. (2017) are shown in green, orange, and blue, respectively. The red line represents the local-universe locus of galaxies as shown by   Kewley et al. (2001) and the blue line to Strom et al. (2017). vide the sample in three different gas-phase metallicity bins with roughly equal number of sources and median metallicities given by Z gas = 0.009, Z gas = 0.011 and Z gas = 0.014. We perform SED fitting including the UV-to-FIR photometry and the Hα, Hβ, and [OIII]λ5007 emission lines by fixing the gas-phase metallicity in each bin to its median value. Only two stellar metallicities of the BC03 models, 0.02 and 0.008, are included in the full range of gas-phase metallicity of our sample. We set the stellar metallicity to 0.02 after checking that using 0.008 does not affect our parameter estimation. In the case of Hβ we only include measurements of the lines for objects satisfying BD > 2.86 to be consistent with our models. We let the ionization parameter to vary between −4.0 < log U < −1.0 and we use n e = 100 cm −3 consistent with the electron density derived by Kashino et al. (2017) for the FMOS-COSMOS sample.
We exclude from the analysis 21 objects with a difference between observed and fitted fluxes larger than 0.2 dex for the three emission lines. In Fig. 17 the observed and estimated fluxes from the three emission lines are compared. In the case of the Hα line as compared to our previous fit in Fig. 6 similar results are obtained and the line fit is not improved by including more emission lines. The distribution of the estimated [OIII]λ5007fluxes is not symmetric with an excess of overestimated fluxes. Hβ estimated fluxes have a standard deviation of ∼ 0.13 dex which is twice as large as that of Hα and [OIII]λ5007 emission lines. We obtain a median attenuation from CIGALE of A Hα = 1.11 ± 0.18 mag and A [OIII] = 1.34 ± 0.22 mag ( † as reported in Table 3), consistent with values derived using photometry and the Hα emission only. Stellar mass, SFR, and attenuation in the Hα and [OIII]λ5007 emission lines are not find to vary significantly by including more emission lines in the SED fitting process as shown in Fig. 18.
Using mock catalogs created with CIGALE we investigate the robustness in the estimation of the ionization parameter log U. The exact and estimated values of log U are compared in the bottom-right panel in Fig. 17 showing a good agreement within a 0.6 dex (1σ) dispersion for the three different gas-phase metallicity bins guaranteeing the reliability in the estimation of the ionization parameter.
From our fit, we estimate median values of log U 0.009 ∼ −2.74, log U 0.011 ∼ −3.03, and log U 0.014 ∼ −2.79, while including the Hβ, and [OIII]λ5007 lines. The 16th and 84th percentiles in logU for each gas-phase metallicity bin are log U 0.009 ∼ −3.02, −2.35, log U 0.011 ∼ −3.21, −2.60, and log U 0.014 ∼ −3.14, −2.00, respectively. Kaasinen et al. (2018) measured log U ∼ −2.72 from an evolutionary analysis of the ionization parameter using a sample of 50 star-forming galaxies selected from the FMOS-COSMOS catalog with DEep Imaging Multi-Object Spectrograph (DEIMOS) observations. Sanders et al. (2020) constrain the ionization parameter using BPASS grids for a MOSDEF and Keck Baryonic Structure Survey (KBSS) samples at z ∼ 1 − 3 to log U ∼ −2.63 and log U ∼ −2.85, respectively. Topping et al. (2020) found that the local 12 + log(O/H)log U relationship of Pérez-Montero (2014) for low-redshift galaxies is still applicable at z ∼ 2. From this relation, the ionization parameter range for our sample should span over the range −3.2 < log U < −2.5, in agreement with our estimations.
The SFR-L [OIII]λ5007 ratio can be interpreted as the Hα-[OIII]λ5007 ratio because Hα is a tracer of SFR. The dispersion in the [OIII]λ5007/Hα dust corrected ratio (Fig. 13) depends on both ionization parameter and gas-phase metallicity. Now we explore the influence of the ionization parameter on our previous relation between the SFR and [OIII]λ5007 (Sect. 5). In Fig. 19, we present SFR/L [OIII]5007 from Eq. 6 as a function of the ionization parameter that we derived from the fit with emission lines and the three different median metallicities in each bin. In this figure, the error bar symbols represent the median values for each bin of gas-phase metallicity with the [OIII]λ5007, Hβ emission lines. The three different gas-phase metallicity bins are presented as blue circles, orange squares, and green triangles. The black crosses correspond to excluded data with flux difference larger than 0.2 dex. The black line corresponds to the 1:1 relation. The three emission lines are well fitted for the three different median gasphase metallicity models. The last panel shows the estimated vs exact value for log U from mock samples created with CIGALE. Symbols are the same as the legend in the first panel and the median ionization parameter value is shown. The shaded area corresponds to the standard deviation. standard deviation as the error. We measured median ratios of ([OIII]λ5007/Hα) 0.009 = 1.12, ([OIII]λ5007/Hα) 0.011 = 0.79 and ([OIII]λ5007/Hα) 0.014 = 0.54 in each gas-phase metallicity bin. However similar average log U values can lead to different [OIII]/Hα ratios due to variations with gas-phase metallicity. In Table 5, we show median values of the SFR-L [OIII]λ5007 ratio for each gas-phase metallicity in 0.5 dex log U bins for the −3.5 < log U < −2.5 range in which the ratio seems to be quite stable. We measure a mean dispersion of 0.24 dex in gas-phase metallicity and 1.1 dex for the ionization parameter in the range −3.5 < log U < −2.5. The effects of the ionization parameter on the SFR-L [OIII]λ5007 dispersion are dominant, but our sample covers only a range of 0.4 dex in the gas-phase metallicity. A sample spanning over a larger range of metallicity is needed to explore in detail the relative influence of both parameters. Our results are consistent with previous studies (e.g., Pérez-Montero 2014; Kaasinen et al. 2018;Sanders et al. 2020;Topping et al. 2020).
In a future work, we plan to explore the effects of using SED fitting with CIGALE implementing BPASS to explore how the locus of galaxies in the [NII]-BPT is affected, leaving the abundance ratio of (N/O) as a free parameter to introduce flexibility in HII-region models and checking also different ways of modeling the emission lines. The SED fitting and HII-region model coupling remains paramount in order to perform homogeneous analysis of a sample of galaxies and break degeneracies between the different parameters involved. Fully understanding of working with spectro-photometric samples and SED fitting is needed as a preparation for new instruments such as PFS and MOONS.  . SFR and L [OIII]5007 ratio vs the ionization parameter. Each metallicity bin is presented as blue dots, orange squares, and green triangles. The ionization parameter is computed with CIGALE for each fixed metallicity case. The black dashed line corresponds to the −41.20 [M yr −1 erg −1 s] intercept found in Eq. 6. The symbols with errors represent the median values of the SFR/L [OIII]5007 ratio for each metallicity. The ionization parameter has a larger impact on the dispersion than the gas-phase metallicity.

Summary and Conclusions
In this work, we perform SED fitting using CIGALE on an FMOS-COSMOS spectro-photometric sample covering UV-to-FIR continuum emission with 21 broad-band fluxes and emission lines at z ∼ 1.6. A sample of 183 objects was selected to have flux measurements of both Hα and [OIII]λ5007 at S/N > 3 in the FMOS survey. From SED fitting of both photometric and Hα fluxes, we estimate SFR and stellar mass, and constrain dust attenuation affecting the Hα and [OIII]λ5007 emission lines. We measure median values of the emission line attenuation of  [OIII] increase with stellar mass with a larger attenuation correction for [OIII]λ5007 emission line as compared to Hα. We find a flatter effective attenuation curve than the Milky Way or C00 curves. A relation to obtain the attenuation for the [OIII]λ5007 line as a function of stellar mass is proposed in the same way as it exists for Hα. This relation could be useful for inferring average values of the SFR of a sample, but not for individual galaxies. The relative attenuation affecting different populations is characterized by the µ parameter in the attenuation law. We find a value of µ = 0.57 ± 0.14 consistent with different works in the literature and twice as large as compared to the original value proposed by Charlot & Fall (2000).
A SFR-L [OIII]λ5007 dust-corrected relation is derived. We measure a slope consistent with unity within the 2σ dispersion of the relation. We estimate a [OIII]λ5007/[OIII] 88 µm ratio of 1.90 for our sample of galaxies and deduce A SFR-[OIII] 88 µm relation in agreement with previous relations that were found at both low and high redshifts, although [OIII]λ5007/[OIII] 88 µm is strongly dependent on electron density and gas-phase metallicity. The SED fitting, including photometry, Hα, Hβ, and [OIII]λ5007 fluxes, is also performed with a refined grid of photoionization models and metallicities estimated from the massmetallicity relation from Curti et al. (2020). The variations of gas-phase metallicity and ionization parameter induce a dispersion in the SFR-L [OIII]λ5007 relation of 0.24 dex and 1.1 dex, respectively. The lower impact of gas-phase metallicity is likely to be due the limited range of our sample (0.006 < Z gas < 0.016) and our relation of SFR-L [OIII]λ5007 is expected to be only valid for galaxies of similar gas-phase metallicities as those studied in this work.