Insights into the reionization epoch from cosmic-noon-Civ emitters in the VANDELS survey

For our analysis we have used data from ”VANDELS: a VIMOS survey of the CDFS and UDS fields”, which is a recently completed ESO public spectroscopic survey carried out using the VIMOS spectrograph on the Very Large Telescope (VLT). VANDELS targets are selected in the UKIDSS Ultra Deep Survey (UDS) (RA = 2:18, Dec = -5:10), and the Chandra Deep Field South (CDFS) (RA = 3:32, Dec = -27:48): VANDELS footprints are centered on the HST areas observed by the CANDELS program (Grogin et al., 2011; Koekemoer et al., 2011), although given the VIMOS field of view, they are larger. The survey description and initial target selection strategies can be found in McLure et al. (2018), while data reduction and redshift determination methods can be found in Pentericci et al. (2018a). Here we give a brief description of the most important aspects relevant to this work.

VANDELS targeted star forming galaxies at $z>2.4$ as well as massive, passive galaxies in the redshift range $1<z<2.5$, with integration times ranging from 20 to 80 hours depending on the magnitude to ensure a uniform signal-to-noise ratio (S/N) on the continuum, with most targets observed for 40 hours or more. The EasyLife data reduction pipeline was used to reduce the spectra (Garilli et al., 2012), producing the extracted 1D spectra, re-sampled 2D spectra, and sky-subtracted spectra. Spectroscopic redshifts were derived using the pandora.ez tool (Garilli et al., 2010) by cross-correlation with galaxy templates. The reliability of the redshifts is quantified with a quality flag (QF) as follows: 0 means no redshift could be determined; 1, 2, 3 and 4 mean respectively a probability of 50%, 75%, 95% and 100% of the redshift being correct; finally 9 means that the spectrum shows a single emission line. The redshift measurement accuracy of spectroscopic observations is typically 150 km s^-1. For more information, see Pentericci et al. (2018a) for details on data reduction and redshift determination.

In this work we have used galaxies from the fourth and final data release (DR4) which is described in Garilli et al. (2021). The data release contains spectra of approximately 2100 galaxies in the redshift range of $1.0<z<7.0$, with over 70 per cent of the targets having at least 40 hours of on-source integration time.

We base our work on the official VANDELS emission line catalog which will be presented in Talia et al. (submitted). Specifically, we use the catalog that was build using Gaussian fit measurements performed with slinefit¹¹1https://github.com/cschreib/slinefit(Schreiber et al., 2018), an automated code that models the observed spectrum of a galaxy as a combination of a stellar continuum model and a set of emission and absorption lines. A set of templates is linearly combined to best fit the continuum. The templates were built with the Bruzual & Charlot (2003) stellar population models and fit to the EAzY (Brammer et al., 2008) default template set using FAST (Kriek et al., 2009). The code searches for lines around their expected locations given by the VANDELS official redshift. Each line with a S/N higher than 3 is allowed its own offset, with a maximum value set by $\pm 1000$ km/s while lines with a lower S/N are instead fixed at their expected position. Therefore, for each detected line, the catalog is formed with the following information: flux and error, rest-frame equivalent width ($EW_{0}$) and error, and the rest-frame full width at half maximum (FWHM) with error. By using the Monte Carlo technique, these errors are determined by randomly perturbing the galaxy spectrum based on its error spectrum, and then computing the uncertainties from the standard deviation of these perturbations compared to the original data.

For our analysis we selected galaxies with a reliability flag 3, 4 and 9, which show a Civ emission line with a S/N greater than 2.5.

To ensure the reliability of the lowest S/N detections and following (Saxena et al., 2020) we produced a stacked spectrum of the 12 sources with $2.5\leq S/N\leq 3$, as an additional check to ensure that they are indeed bona fide Civ emitters (see also Sect. 4.4). In this way we verified that we can safely consider these sources in the sample. We selected a total of 50 galaxies.

From this list generated by the automatic procedure we removed by eye six sources which present clear problems of sky-residuals in the region of the Civ emission and one source with evidence of merger in the spectrum. Finally from a visual inspection of all VANDELS QF 3, 4 and 9 sources we note that the official catalog failed to detect some galaxies with clear Civ emission lines which were not correctly identified by slinefit. This is probably due to the limited velocity offset allowed in the automatic search and to the fact that the reference redshift used for the search is the ”VANDELS redshift”, which is not the systemic redshift but is in general driven by the strongest features present in the spectrum. In few cases this redshift is closer to the Ly$\alpha$ redshift or to the ISM redshift depending on their relative strength. A further 8 sources were identified and re-analyzed through slinefit by correctly calibrating the redshift (on the Ciii] line when possible) which represents better the systemic redshift of the galaxies.

Since we want to study the nature of Civ in star forming galaxies, we excluded the sources which are obvious AGN: in particular we excluded 8 X-ray detected sources from the CDFS 7 Ms Source Catalogs (Luo et al., 2017) and the UDS X-ray observations (Kocevski et al., 2018), of which 4 are broad line AGNs (BLAGNs: CDFS003360, CDFS019824, UDS017629 and UDS199159) which are easy to recognise from the broad UV emission line, and 4 narrow line AGNs (NLAGNs: CDFS005827, CDFS019505, UDS025482 and UDS018960). For an individual analysis of these sources see Bongiorno et al. (in prep.).

In conclusion, from the parent sample made of 933 galaxies at $2.5\leq z\leq 5$ without the AGN contribution, we obtain a sub-sample of 43 sources, 20 from the UDS field, and 23 from the CDFS field. The IDs, RA and DEC coordinates, redshifts, related flags and emission line information are presented in Table 1. This sample might still contain some NLAGNs with no X-ray counterpart. In the following section we will refine our analysis and employ emission line ratio diagrams to obtain a sample of purely star forming Civ galaxies.

In Fig. 1 we show the distribution of the galaxies’ $H_{AB}$ magnitude as a function of the redshift for the Civ sample (highlighted in purple) and the parent sample (grey symbols), i.e. all other VANDELS sources after removing the AGNs, as classified in Bongiorno et al. (in prep.). For this comparison the parent sample has been cut between redshift 2.5 and 5. As we can see the Civ emitters are a representative subset of the VANDELS population at least in terms of H-band and redshift distribution.

The models developed by Feltre et al. (2016) for the emission by the AGN narrow-line regions have been created using Cloudy (Ferland et al., 2013)²²2These models can be found at http://www.iap.fr/neogal/agn-models.html. They consider a parameterization of the gas distribution as a cloud of a single type but contaminated by dust and dominated by the radiation pressure. It is useful to highlight that in these models, the line measurements are not corrected for dust attenuation.

For the purposes of our analysis, we use models with $Z$ = 0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.004, 0.006, 0.02 with no priors on $n_{H}$ in order to identify a clear AGN region in all the UV-BPTs. These models were constructed using the relative abundances of heavy elements for solar metallicity (Caffau et al., 2011). The predicted intensities of 20 optical and ultraviolet emission lines are then produced (in units of erg s^-1). We consider the Civ, Heii and Ciii] predicted intensities. Since we also want to produce a diagnostic diagram of Civ EW with Civ/Heii which was already used by (Nakajima et al., 2018b), we need to also derive the EW of the Civ line from the luminosities provided by the models. To determine the line continuum we followed the procedure outlined in Castellano et al. (2022). We first re-normalize the incident spectra to a chosen fraction of the observed average flux in the R band, i.e. 5950-7250 Å, which corresponds to the observable non-ionizing UV flux at $\sim$ 1500 Å, and then compute the EW of the Civ line on the basis of the predicted line flux for each model and of the observed flux at the relevant wavelength.

It has been argued that the inclusion of interacting binary stars in stellar population synthesis models may hold the key to explain the observations of bright high ionization nebular emission lines in the spectra of high redshift star forming galaxies (Stanway et al., 2016; Eldridge et al., 2017). Xiao et al. (2018) combine the output of the stellar spectral synthesis model bpass (Binary Population and Spectral Synthesis) used in conjunction with the photoionization code Cloudy to predict the nebular emission from H ii regions around young stellar populations over a range of compositions and ages.

In total, there are 80262 photoionization models, including separate models for binary and single stars³³3These models can be found at https://bpass.auckland.ac.nz/4.html. For our analysis, we consider the photoionization models for binary stars with $Z$ = 0.001, 0.002. 0.02 and $\log(\text{Age}/\text{yr})$ = 6, 6.5, 7, 7.5. For each model, the predicted luminosities and EWs of 17 optical and ultraviolet emission-lines are given (in units of erg s^-1 and Å, respectively). We retrieved the Civ, Heii and Ciii] predicted intensities and EWs.

In addition, we also consider the models developed by Gutkin et al. (2016) for the nebular emission from star forming galaxies. They have been created using Cloudy (Ferland et al., 2013). For our analysis, we use models with $Z$ = 0.0001, 0.0002, 0.001, 0.002 and 0.02 with no priors on $n_{H}$ and with relative elemental abundances fixed at Solar ratios (Caffau et al., 2011). The predicted intensities of 20 optical and ultraviolet emission-lines are then produced (in units of erg s^-1)⁴⁴4These models can be found at http://www.iap.fr/neogal/sf-models.html. As in the previous cases, we consider the Civ, Heii and Ciii] predicted intensities.

Similarly to what done with the Feltre et al. (2016) models, we derive also the EW of the Civ line from the luminosities provided by the models, employing the method described in Sec. 4.1 and in Castellano et al. (2022).

We now compare the properties of our sample of Civ emitters with the predictions from photo-ionisation models for nebular emission from star forming galaxies or AGN narrow-line regions described in the previous sections. Following the original diagnostics proposed by Feltre et al. (2016) we will study the two relations Ciii]/Heii vs Civ/Heii and Civ/Heii vs Civ/Ciii]. In addition to these, and following the work by Nakajima et al. (2018b), who constructed photoionization models to interpret the UV spectra of Ciii emitters selected from the VUDS survey at $z=2-4$, we will also use the relation Civ/Ciii] versus (Ciii]+Civ)/Heii, and Civ EW vs Civ/Heii.

In Fig. 4 we present the four different diagrams, which for simplicity we refer to UV-BPT1, UV-BPT2, UV-BPT3 and UV-BPT4 from here on. These diagrams give the advantage that objects can be tracked and diagnosed even if only one of the three lines is detected - in fact, if the S/N of a line is $<2$ we consider 2$|\sigma|$ as a limit. Whenever both lines considered in a certain line ratio are not detected but are both limits, we do not plot that object in the relevant UV-BPT diagram. This means that only UV-BPT3 and UV-BPT4 contain all the 43 objects in our sample.

In each diagram, the dots with error bars and/or limits indicate the position of the sources in our sample; the green triangles represent the position of the X-ray AGNs which have been excluded from our catalog and the green and magenta dots are the predictions from the models described in the previous subsections, respectively from the AGN models (Feltre et al., 2016), and the SF models (Xiao et al., 2018; Gutkin et al., 2016). It can be seen that the models occupy regions that are quite well although not completely separated in each of the UV-BPT diagrams. In the UV-BPT1 and UV-BPT4 diagrams we also plot the separation between AGN and SF galaxies proposed by Nakajima et al. (2018b) using the distributions of galaxies and AGNs observed in the VUDS survey.

Based on the models considered and on the position of the X-ray detected AGN, we proposed a slightly revised version of the separation between AGN and SF galaxies in the UV-BPT1 diagram which can be written as:

For the other UV-BPT diagrams we consider the following separation lines: for UV-BPT2

For UV-BPT3 (not considered by Nakajima et al., 2018b):

For UV-BPT4 (same as in Nakajima et al., 2018b),

Not all of 43 sources are in all the diagrams but, when present, their position is always consistently in the AGN or SF part of the diagrams for all cases apart from UV-BPT4. In fact we note that the results from UV-BPT4 are more ambiguous since most sources, including some of the X-ray detected NLAGNs are located in an intermediate region in between the AGN and SF models. In addition, it seems that none of the SF models, can produce Civ EW as large as those observed in our sample. This problem was already highlighted in Saxena et al. (2020), where we noted that although the low-metallicity binary star models are able to reproduce the UV line ratios of their Heii emitters, they under-predict the Heii EWs. The same applies to the Civ emission which is under-predicted by bpass, as also discussed in Nakajima et al. (2018b); Le Fèvre et al. (2019). For the separation between AGN and SF galaxies we therefore base our classification on the bases of the other 3 UV-BPT diagrams. Specifically we find:

• •

  4 sources whose position is consistently in the AGN regions in all UV-BPT diagrams: we call these sources ”certain AGN” and indicate them with blue dots;

• •

  10 sources whose position is within the AGN regions, but given that the Ciii and/or Heii emission lines are actually limits, could move into the SF regions. We call these sources ”possible AGN” and indicate them in magenta;

• •

  29 sources whose position is consistently in the SF region for all UV-BPT diagram. We call of course these sources ”SF” and indicate them in black.

We note that of the 8 X-ray detected LAGNs, 1 source is consistently located within the SF region, while the other are in the AGN region. This could be due either to a limitation of the models in reproducing all AGN properties or to the fact that in many objects the UV emission line might be of mixed origin.

To further shed light on the real nature of the 10 ”possible AGN” we apply a spectra stacking technique to infer the average properties of these sources.

Following Saxena et al. (2020), the stacking is performed by first converting each spectrum to the rest-frame, using the spectroscopic redshift given by the DR4 catalogue. We point out that the VANDELS redshifts are not the systemic ones, since they are derived by a match to spectral templates and they are driven by absorption and/or emission lines depending on their strength (Garilli et al., 2021). For this reason in the stacks the various lines might not appear at the position expected, e.g. the Ciii is not at the exact redshift as we would expect if we had used the systemic redshift for each source. For our analysis of the stacks we neglect this effect. The rest-frame spectra are first normalised using the mean flux density value in the rest-frame wavelength range 1300 Å$\ <\lambda<1500$ Å. The spectra are then re-sampled to a wavelength grid ranging from 1200 to 2000 Å, with a step size of 0.582 Å(the wavelength resolution obtained at a redshift of 3.6, which is the median redshift of VANDELS sources in this work). The errors on the stacked spectra are calculated using bootstrapping: we randomly sample and stack the same number of galaxies from the sample 500 times, and the dispersion on fluxes thus obtained gives the errors (Fig. 5).

We produce 3 different stacks, one of the possible AGN (10 objects), one with the certain AGN (4 sources) and one with the purely SF galaxies (29 sources). The spectra obviously reveal a wide range of emission and absorption features, with the Ly$\alpha$ emission line standing out in each of the sub-samples. The main UV emission lines are fitted with a Gaussian to determine their fluxes and the Civ EW, in order to place the stacks in the various UV-BPT diagrams, an the results are presented in Fig. 6. We show just one UV-BPT them since the results are consistent for all of them: the blue star represents the sure AGNs, the magenta one the possible AGNs, and the black one the SF galaxies (in accordance with the colors used in Fig. 4). As in the previous UV-BPTs, the green dots are obtained from the AGN models from Feltre et al. (2016) and the red dots from the SF models from Xiao et al. (2018); Gutkin et al. (2016), while the grey lines are the division derived in this work and the green dashed lines refer to Nakajima et al. (2018b). We note again that the AGN stack is not well reproduced by the AGN models but rather sits in an intermediate position.

In the SF stacked spectrum (Fig. 5) we note that the Heii emission is basically unresolved (FWHM ¡600 km $\text{s}^{-1}$). This confirms the nebular nature of the Heii: in fact stellar Heii is expected to be much broader of the order of 1500-2000 km $\text{s}^{-1}$ as seen for example in the stacks of the general LBG population (Shapley et al., 2003). For example Cassata et al. (2013) used a threshold of 1200 km $\text{s}^{-1}$ to identify narrow Heii emitters of nebular origin. Thus, we can exclude a stellar origin for our sources.

The properties of the stacked spectra are as expected in the sense that the magenta star is at an intermediate position between the other two: however in all diagrams it is more consistent with the SF regions, indicating that the 10 sources are dominated by SF galaxies rather than AGN.

Another check on the nature of the ”possible AGN” sample can be obtained by the Nv emission. For AGNs, as predicted from theoretical and observational studies (e.g., Hamann & Ferland, 1992, 1993; Hamann et al., 2002), we expect a bright Nv line. Therefore, we determine the Ly$\alpha$/Nv emission line ratio for our three different stacks. For the SF galaxies stacked spectrum the Nv line is not detected at all and we can place a limit Ly$\alpha$/Nv $>30$. The stack of the possible AGN presents a Ly$\alpha$/Nv = $8.7\pm 0.6$ that is considerably higher than the certain AGN one (Ly$\alpha$/Nv = $2.1\pm 0.6$), although we have a 2$\sigma$ detection of the Nv. This value is higher than the ratios measured in narrow line AGNs (Ly$\alpha$/Nv $\sim 1-2$, Tilvi et al., 2016; Hu et al., 2017; Sobral et al., 2017).

We finally exploit the Chandra X-ray observations: although we have excluded the individually detected X-ray sources since the very beginning, the Luo et al. (2017) sources are those detected at 5$\sigma$, therefore we cannot exclude faint emission from the ”possible AGN” sample. We therefore stacked the X-ray data to check if we detect any flux from the combination of the sources. Unfortunately we can do this in a meaningful way only for the 6 CDFS sources which have much deeper X-ray observations (7 Msec, Luo et al., 2017) while for the UDS field we only have much shallower observations (Kocevski et al., 2018). To stack the X-ray data we follow the same procedure described in Saxena et al. (2021). The stack does not produce any detection, although the limit is not very stringent since there are only 6 sources of which one lies at the border of the Chandra field.

In conclusion, from the position of the stack in the four UV-BPT diagrams, the high Ly$\alpha$/Nv ratio and the absence of X-ray emission in the stack, we conclude that the sample of 10 possible AGN must be dominated by SF galaxies, although we cannot exclude a (faint) AGN contribution. In the following we include these 10 sources in the SF sample, bringing it to 39 galaxies in total, to study their physical properties.

The code used for the spectral energy distribution (SED) fitting is the BEAGLE (BayEsian Analysis of GaLaxy sEds) software (Chevallard & Charlot, 2016) which was chosen for its versatility: through the Bruzual & Charlot libraries it allows us to model in a coherent way the continuous radiation coming from the stars of the galaxies and the nebular component which is modeled with Cloudy. The code also includes several prescriptions to describe the dust attenuation and the radiation transfer through the intergalactic medium. Stellar formation and chemical enrichment histories are included both with parametric and non-parametric prescription. Finally, BEAGLE can simultaneously fit the photometry and the spectra (or individual spectral indices) of the galaxies. BEAGLE was run on all the VANDELS parent sample using v0.24.5 which employs the most recent version of the Bruzual & Charlot (2003) stellar population synthesis models (see Vidal-García et al., 2017, for details). The full description of the method and the results on the entire sample will be discussed in a companion paper (Castellano et al. in prep.). We provide here only a brief description of the most important parameters adopted: in particular we assume a delayed star formation history (SFH), which is the most flexible parametric SFH used by the code, with a current SF timescale (i.e. duration of the current episode of star formation) of 10 Million years. The escape fractions $f_{esc}$ assumed in the BEAGLE run are set to 0 for all galaxies, regardless of the presence of the Civ emission. We used a Chabrier (2003) Initial Mass Function, metallicity in the range $-2.2\leq\log(Z/Z_{\sun})\leq 0.25$ and we treated the attenuation following the Charlot & Fall (2000) model combined with the Chevallard et al. (2013) prescriptions for geometry and inclination effects.

For each source, BEAGLE is run on the combination of photometry, coming from the official 13 band VANDELS catalog presented in Garilli et al. (2021), and spectral indices i.e., the fluxes of the Heii, Oiii and Ciii emission lines as well as the precise spectral windows used (respectively [1634-1654], [1663-1668] and [1897-1919]). After some tests we decided not to provide the measured Civ emission to BEAGLE as a further spectral index to fit, both to have an unbiased estimate of the physical parameters, and also because as just discussed in Sect. 4.3, there is a well known inability of stellar population models, even when including binary stars, to reproduce some UV line fluxes (also see Saxena et al., 2022a). We also tested how the choice of current SF timescale could change our results. In particular this is relevant since the EWs of the emission lines are strongly dependant on the current SF. As a result even when reducing the SF timescale to 5 Myrs, the physical parameters remain very stable with the exception of $\xi_{ion}$, which becomes slightly higher ($\sim$ 0.1 dex) for all galaxies. Note that the inferred $\xi_{ion}$ from BEAGLE is almost independent from the initially assumed $f_{esc}$: this value is determined by the intrinsic ionizing SED, which is already adjusted for relatively low metallicities and young ages. The change in $\xi_{ion}$ is substantial only for very high $f_{esc}$ (see e.g., Marques-Chaves et al., 2022). The official VANDELS redshift was assumed in all cases for the SED fitting.

As a result, BEAGLE provides the best SED fit and the best physical parameters, with relative uncertainties: specifically total stellar mass, star formation rate, dust attenuation, age and $\xi_{ion}$ were retrieved. For a complete list of the selected physical parameters and the range provided a priori to fit all the sources of the sample, see Table 2 in the Appendix.

Using the derived physical properties we searched for systematic differences between the SF sample of the Civ emitters (39 galaxies) and the parent sample. From the parent sample we have removed all X-ray emitters and all possible AGN as defined in Bongiorno et al. (in prep.). In Fig. 7 we present the location of our Civ emitters on the main sequence i.e. the stellar mass vs SFR plot: given that our sample spans a rather large redshift range, we split the plot into two redshift bins, $2.7<z\leq 4$, and $4<z\leq 5$. The blue lines are obtained using the main sequence as a function of redshift described in Santini et al. (2017). The VANDELS sources as well as the Civ emitters sub-sample are located more or less on the main sequence of star forming galaxies, with perhaps a small bias in the high redshift bin, where the (few) objects are mostly located on or above the MS.

In Fig. 8 we then show the normalized distributions of the derived stellar masses, the SFRs, the photon production efficiencies $\xi_{ion}$ and the stellar ages of the Civ galaxies and the parent population. In each panel we also indicate the mean and median values for the two samples with vertical regions, whose widths represent the uncertainty. From the panels we can see that the Civ emitters tend to have considerable lower masses, slightly lower SFR, higher $\xi_{ion}$ and slightly lower stellar ages compared to the parent sample. To assess the significance of these differences, we employ the common Kolmogorov-Smirnov (KS) test. From the p-values (reported in the figures) we conclude that the mass and $\xi_{ion}$ of the Civ emitters and the parent sample are statistically different, while the max_stellar_age and the SFR are compatible with being randomly drawn. In Castellano et al. (in prep.) we find a strong link between the sSFR and $\xi_{ion}$: indeed if we look at our Civ emitters they have slightly higher sSFR compared to the parent population, which could explain the higher $\xi_{ion}$. This link and its implications will further discussed in that paper. Even using a SF timescale of 5 Myrs, with all values of $\xi_{ion}$ being slightly higher, the difference between Civ emitters and the parent sample is still significant.

In conclusion our Civ emitters are in general low mass objects which have an higher photon production rate compared to the general star forming population at the same redshift.

The link between the presence of nebular Civ emission and low metallicity stellar population has been noted both at high and low redshift. For example Senchyna et al. (2017) showed that Civ emission may be ubiquitous in extremely metal-poor galaxies with very high specific star formation rates, although it reaches large EWs only if the ISM is exposed to the hard stellar ionizing spectra produced at extremely low metallicities (12 + log O/H $\lesssim$ 7.7).

To derive the stellar metallicity from our spectra, we use the method developed by Calabrò et al. (2021), which is based on stellar photospheric absorption features at 1501 Å and 1719 Å, calibrated with Starburst99 models and are largely unaffected by stellar age, dust, IMF, nebular continuum, or interstellar absorption. The method can be applied to spectra with a S/N of 15-20 in the continuum to avoid large uncertainty, and therefore we cannot derive metallicities for individual sources but only average values from stacked spectra. We produced a stack of the 39 SF galaxies following the procedure outlined in the previous section. The metallicity we derive is $Z\sim 0.26Z_{\odot}$ for the sample of 29 purely SF galaxies and a slightly lower (but consistent within the uncertainties) value of $Z\sim 0.1Z_{\odot}$ for the combined sample of 39 galaxies. Looking at the stellar mass-stellar metallicity relation by Calabrò et al. (2021) we can see that these low metallicities (about $0.15Z_{\odot}$) are perfectly consistent with the general trend for $z\sim$ 3 galaxies with similar masses, i.e. the Civ emitting galaxies do not show lower metallicities compared to their parent sample, as also found by Saxena et al. (2022a). This is therefore at odds with the low redshift universe where Civ emitters have almost always an inferred stellar (and gas phase) metallicities are extremely low (e.g., Berg et al., 2019; Senchyna et al., 2021). However in terms of ”absolute” metallicity, $z\sim$3 Civ emitters have a similar metallicity as the local universe emitters: this is of course due to the fact that at $z\sim 3.5$ the average metallicity of star forming galaxies is a factor of 0.6 dex lower than in the local universe (Cullen et al., 2019; Calabrò et al., 2021).