© The Authors 2024

1. Introduction

During its six-year mission to constrain the dark Universe, Euclid (Laureijs et al. 2011; Racca et al. 2016; Euclid Collaboration 2024b) will catalogue billions of galaxies and collect an unprecedented abundance of highly accurate photometric and spectroscopic data. Using weak lensing and galaxy clustering as cosmological probes, it will study the growth of cosmic structures and the Universe’s accelerated expansion over the past ten billion years. For the purpose of recovering accurate distance measurements, Euclid’s near-infrared spectrometer and photometer (NISP, Maciaszek et al. 2022; Euclid Collaboration 2022b, 2024a) was designed to probe redshifted optical emission from galaxies out to redshift two by observing in the near-infrared range. Crucially, the regime around redshift two represents the peak of star formation (SF, Madau & Dickinson 2014), black hole growth, and quasar activity (Richards et al. 2006), making the resulting data set ideal for studying galaxy formation and evolution.

In the Euclid Wide Survey (EWS, Euclid Collaboration 2022a), Euclid is set to observe roughly 15 000 deg² of extragalactic sky. The NISP spectrometer has been tuned to measure Hα line emission at redshift 0.84–1.88 and is expected to recover spectra for about 35 million galaxies. It contains three ‘Red’ GriSms (together denoted RGS, resolving power ℛ > 480) oriented at different angles, each covering rest-frame 1.21–1.89 μm to a flux limit of 2 × 10⁻¹⁶ erg s⁻¹ cm⁻². Euclid’s initial specification forecasts a number density of 1700 deg⁻² Hα emitters; however, this estimate strongly depends on the uncertain intrinsic Hα luminosity function in this redshift range (see Pozzetti et al. 2016; Bagley et al. 2020).

In addition to the EWS, Euclid will cover selected fields of 50 deg² at two magnitudes deeper in the Deep Survey (EDS; Scaramella et al., in prep.). In this mode, the NISP spectrometer is able to measure emission lines at fluxes greater than 6 × 10⁻¹⁷ erg s⁻¹ cm⁻² and observations will be made in a second ‘Blue’ GriSm (BGS, ℛ > 400), covering rest-frame 0.93–1.37 μm. These capabilities make the EDS ideal for performing detailed sample characterisations, as it allows the detection of fainter emission overall, as well as the simultaneous recovery of the most useful rest-frame optical emission lines for galaxies at redshift 0.4–2.5.

Strong optical emission lines, such as [N II]λ6584, Hα, [O I]λ6300, [O III]λ5007, Hβ, and the doublets

[O II]λλ3727, 3729 and [S II]λλ6717, 6731, have long been known to be particularly sensitive probes of both the local conditions of the ionised gas in the interstellar medium (ISM), as well as the nature of the ionising radiation (Ferland & Netzer 1983; Osterbrock & Ferland 2006, see Kewley et al. 2019 for a recent review). As a result, emission-line intensities can be used in spectroscopic diagnostics to trace various galaxy properties.

Diagnostic diagrams combining two emission-line ratios are widely used to determine whether the ionising radiation in a galaxy is dominated by young, massive stars (produced in recent SF) or by an active galactic nucleus (AGN). The standard Baldwin–Phillips–Terlevich (BPT, Baldwin et al. 1981; Veilleux & Osterbrock 1987) diagrams, which connect the [O III]λ5700/Hβ ratio to the [N II]λ6584/Hα, [S II]λ6724/Hα, and [O I]λ6300/Hα ratios, have proved successful at distinguishing between ionising sources in local galaxies (Kewley et al. 2001; Kauffmann et al. 2003). With large-scale spectroscopic surveys, such as Euclid, collecting high-quality spectra in the more distant Universe, it remains unclear whether their use can be extended to higher redshifts.

In fact, theoretical works utilising photoionisation models have indicated that in metal-poor galaxies, which are more prevalent at high redshift (see Maiolino et al. 2008), AGNs produce similar optical emission-line strengths to young stars in SF galaxies, leading them to overlap on the [O III]λ5700/Hβ versus [N II]λ6584/Hα BPT diagram as early as redshift one (Groves et al. 2006; Feltre et al. 2016; Hirschmann et al. 2019). Recently, Kocevski et al. (2022) and Harikane et al. (2023) seemingly confirmed this using JWST NIRSpec spectroscopy, as their sample of emission-line measurements of faint AGNs above redshift five and four, respectively, is indistinguishable from SF galaxies in the BPT diagram. However, at redshift 2.3, Coil et al. (2014) found that their sample of 50 SF galaxies and ten confirmed AGNs is still robustly separable in the standard BPT diagram, indicating that its breakdown as a spectral diagnostic might only occur beyond intermediate redshifts. Considering the limited sample size and redshift coverage, it is nevertheless necessary to verify if EDS-like galaxy populations conform to BPT selection criteria and can thus be classified according to their ionising sources in upcoming data releases.

Optical emission lines have also been used to estimate properties of ionising sources, such as the SFR for SF-dominated galaxies and the intrinsic AGN luminosity L_AGN for AGN-dominated galaxies. Hα is a particularly appealing tracer for the SFR, as it is well-calibrated for local galaxies (see Kennicutt 1998; Hopkins et al. 2003; Kennicutt & Evans 2012) and could potentially be used to derive SFRs for galaxies observed in the EWS. In the absence of Hα measurements, the [O III]λ5007 line luminosity (often combined with Hβ into [O III]λ5007 + Hβ, as they are inseparable in photometric narrow-band surveys such as HiZELS, see Geach et al. 2008; Sobral et al. 2009) has been used as a tracer for the SFR (e.g. Teplitz et al. 2000; Moustakas et al. 2006; Osterbrock & Ferland 2006; Sobral et al. 2015). However, due to degeneracies, additional constraints on metallicity and the ionisation parameter are required (Moustakas & Kennicutt 2006; Villa-Vélez et al. 2021). Some studies (e.g. Kennicutt 1992; Sobral et al. 2015) have also warned about potential biasing due to dust and AGN contributions, indicating that [O III]λ5007 is an unreliable SFR proxy, especially at high redshift. Despite this, there is tentative evidence from line-emitting galaxies at redshift 0.84, 1.42, 2.23, and 3.24 in the HiZELS survey that both Hα and [O III]λ5007(+Hβ) can be used to estimate the SFR at moderate and high redshift (Sobral et al. 2013, 2015; Khostovan et al. 2015; Suzuki et al. 2016).

Across various AGN types in nearby galaxies, the [O III]λ5007 luminosity has also been found to correlate with the 2–10 keV X-ray AGN luminosity L_X (e.g. Netzer et al. 2006; Panessa et al. 2006; Lamastra et al. 2009; Georgantopoulos & Akylas 2010; Feltre et al. 2023), which is itself used as a proxy for the bolometric AGN luminosity L_AGN. Thus far, there has been no work done to verify its applicability to redshifts greater than 1. Given the reported evolution of the [O III]λ5700/Hβ for SF galaxies, it is unclear whether similar effects could pollute the correlation in AGN-dominated galaxies.

While line emission in AGN-dominated galaxies is mainly driven by the central AGNs, there may still be a significant identifiable contribution from the star-forming component. One way to quantify the relative influence of the AGN is via the ratio of the black hole accretion rate (BHAR) and SFR. The BHAR-SFR relationship has been studied extensively to constrain the co-evolution of the black hole and its host galaxy, with values of log₁₀(BHAR/SFR) ranging from −4 to −1 (see McDonald et al. 2021 and references therein). As a result of varying methods and assumptions, combined uncertainties from separate BHAR and SFR estimates are large. Thus, for the purposes of source characterisation in the EDS, a direct and consistently measured estimate of this ratio from emission-line intensities should place such measurements on to a more secure footing.

The gas-phase metallicity (often expressed as the oxygen-to-hydrogen abundance O/H) is another key property which imprints onto the emission from ionised gas in galaxies. Various calibrations for local galaxies, derived from both direct temperature (T_e) estimates and photoionisation models, relate intensity ratios of strong emission lines to the O/H abundance (early works by Jensen et al. 1976; Pagel et al. 1979, for recent reviews see Kewley et al. 2019; Maiolino & Mannucci 2019). These relations exhibit a significant scatter at low redshift and recent JWST/NIRSpec observations of galaxies at redshift 2–9 show that O/H estimates derived from low-redshift calibrations can differ significantly from more robust direct T_e estimates (Curti et al. 2022; Sanders et al. 2024). This may indicate a significant difference in the metallicity-related properties of the ISM of low- and high-redshift galaxies. Using photoionisation models coupled to the cosmological IllustrisTNG simulations, Hirschmann et al. (2023b) also found that some line ratio-metallicity relations evolve by up to 1 dex between redshift two and zero. It remains to be clarified how far exactly different calibrations for the gas-phase metallicity relations can be extended from the local Universe before starting to break down.

In summary, rest-frame optical emission lines are powerful probes with which to characterise galaxies and, consequently, the large number of upcoming Euclid spectra at intermediate redshifts will help constrain one of the most important regimes for galaxy evolution. However, as outlined, many locally used spectroscopic diagnostics and emission line-based calibrations are yet to be validated in this domain. Recently, Euclid Collaboration (2023) assessed the performance of the NISP red grisms using mock-spectra constructed by combining galaxy properties from spectral energy distribution (SED) fits of star-forming galaxies between redshift 0.3 and 2.5 and some of the calibrations detailed above, thus explicitly assuming their validity at intermediate redshifts. In order to strengthen these pre-launch forecasts and guide observers in their analysis of future Euclid data releases, it is vital to complement calculations based on empirical relations with self-consistent theoretical frameworks, which allow for the study of emission-line properties across cosmic time and make targeted forecasts for specific surveys and instruments. Predicting these emission lines from first principles in a self-consistent and robust manner has been a long-standing challenge, precisely because of the scarcity of spectroscopic data at intermediate redshifts, in addition to the complex interplay of various physical processes.

Past studies have demonstrated success in coupling nebular emission-line models to cosmological simulations and semi-analytic models. These have thus far been limited to modelling only the line emission due to young stars (e.g. Orsi et al. 2014; Shimizu et al. 2016; Wilkins et al. 2020; Pellegrini et al. 2020; Garg et al. 2022; Baugh et al. 2022), or, if including the contribution from AGN narrow-line regions, are limited in statistics (Hirschmann et al. 2017, 2019) or focus their predictions on specific emission-line properties (Hirschmann et al. 2023a) and high-redshift galaxies (Hirschmann et al. 2023b). Consequently, a lack of comprehensive theoretical guidance for intermediate redshifts persists, which ideally would account for emission-line contribution from AGNs and provide adequate statistics.

In this paper, we aim to close this gap by adopting a Euclid-like mock light cone constructed from the GAEA (GAlaxy Evolution and Assembly, De Lucia et al. 2014; Hirschmann et al. 2016; Fontanot et al. 2020) semi-analytic model, which we couple to photoionisation models used in previous works by Hirschmann et al. (2017, 2019, 2023a,b). Our framework is uniquely successful in its self-consistent modelling of emission lines originating not only from young stars (Gutkin et al. 2016), but also from AGN narrow-line regions (Feltre et al. 2016), and evolved post-asymptotic giant branch stars (Hirschmann et al. 2017). We focus our analysis on different redshift intervals between 0.4 and 2.5, in which Euclid will recover various combinations of the brightest and most useful emission lines, such as Hα, Hβ, [S II]λλ6717, 6731, [N II]λ6584, [O I]λ6300, [O III]λ5007, and [O II]λλ3727, 3729. In the following analysis, we aim to address five key points:

1. The biasing of galaxy populations with respect to stellar mass, standard scaling relations, and dominant ionising sources when considering line emission above the defined flux thresholds in the Euclid Wide and Deep Surveys.

2. The ability of optical BPT diagrams to distinguish between dominant ionising sources in EDS-observable galaxies.

3. The applicability of locally used relations between emission-line intensities and ionising properties (i.e., SFR and AGN luminosity) to intermediate redshifts.

4. Optical emission-line ratios that can directly trace the BHAR/SFR ratio.

5. The redshift evolution of optical line-ratio calibrations for interstellar metallicity at intermediate redshifts.

The paper is structured as follows. The theoretical framework is described in detail in Sect. 2. In Sect. 3, we show how we validated our approach by testing its predictions against robust theoretical and observational findings. In Sect. 4, we explore how observing line emitters in the EWS and EDS imposes selection bias effects on the emission-line flux versus stellar mass plane and various standard scaling relations. Additionally, we provide estimates for EDS-observable fractions of line-emitting galaxies, divided into SF and active. In Sect. 5, we verify the use of the standard [O III]λ5700/Hβ versus [N II]λ6584/Hα BPT diagram to determine the dominant ionising sources for the EDS-observable sample. Section 6 demonstrates that, according to our framework, locally defined calibrations between emission-line luminosities and ionising properties continue to perform well at intermediate redshifts. We further establish a strong [N II]λ6584 emission-line dependence of the BHAR/SFR ratio and provide three novel calibrations to [N II]λ6584-based emission-line ratios. In Sect. 7, we predict that the relationship between various line-ratios and the interstellar metallicity undergo a significant evolution between redshift 0 and 2.5. These changes manifest as shifts in normalisation. We discuss potential caveats of our approach in Sect. 8 and summarise our results in Sect. 9.

2. Theoretical framework

2.1. Mock light cones from the GAEA semi-analytic model

The GAEA¹ semi-analytic model (De Lucia et al. 2014; Hirschmann et al. 2016) is a successor to a model first published in De Lucia & Blaizot (2007). Constructed upon dark matter merger trees, it traces the evolution of four baryonic components: stars in galaxies, hot gas in dark matter haloes, cold gas in galactic disks, and the gas component ejected by stellar and AGN-driven winds. GAEA’s physical processes have been updated in multiple versions over the years. In this study, we make use of the most recent realisation described in Fontanot et al. (2020), which added improved black hole (BH) accretion and AGN feedback modelling to the prescriptions for gas cooling, star formation, gas recycling, environmental processes (all from original model in De Lucia & Blaizot 2007), chemical enrichment (updated in De Lucia et al. 2014), and stellar feedback (updated in Hirschmann et al. 2016).

This version of GAEA was run on merger trees extracted from the N-body cosmological Millennium Simulation (Springel et al. 2005), which adopted a box size of 500 h⁻¹ Mpc and WMAP1 cosmological constant-dominated cold dark matter (ΛCDM) concordance cosmology (Ω_Λ = 0.75, Ω_m = 0.25, Ω_b = 0.045, n_s = 1, σ₈ = 0.9, and H₀ = 73 km s⁻¹ Mpc⁻¹). This large box size is crucial to make predictions which are statistically representative, given Euclid’s large areal coverage. The stellar quantities were calculated using the stellar population synthesis model from Bruzual & Charlot (2003) assuming a Chabrier IMF (Chabrier 2003), and the resulting physical quantities were stored in galaxy catalogues corresponding to each simulation snapshot taken at finite redshifts.

To ensure that our GAEA predictions closely match the upcoming Euclid observations, we used those catalogues to construct a mock light cone according to the algorithm described in Blaizot et al. (2005) and Zoldan et al. (2016). To avoid replications, the GAEA boxes at different redshift snapshots were first randomly rotated, shifted, or inverted before placing the model galaxies into an empty light cone with an aperture of 5.27°, which is the largest possible diameter without exceeding the Millennium Simulation box size. Redshift varies continuously between 0 and 3.9 along the light cone and thus galaxies were extracted from the snapshot closest in redshift to the corresponding light cone region.

The resulting GAEA light cone catalogues (GAEA-LC hereafter) include all galaxies with an estimated H-band AB magnitude m_H brighter than 25. We note that the EDS is expected to reach magnitude 26. However, with an apparent magnitude limit of 25, the typical mass-to-light ratios for galaxies below redshift 0.5 translates to a mass below the mass resolution of the original simulation. We assumed a conservative resolution cut for the Millennium Simulation of 10¹¹ M_⊙ for dark matter halos, which translates to an approximate resolution limit in stellar mass of 10⁹ M_⊙ in GAEA. In applying the EWS and EDS flux limits of the NISP spectrometer to individual emission lines predicted by our model, we will demonstrate that the majority of galaxies with masses less than 10⁹ M_⊙ will not be observable with Euclid.

2.2. Modelling of emission lines for GAEA-LC galaxies

The GAEA-LC galaxies were post-processed with photoionisation models based on the Cloudy code (Ferland et al. 2013, version c13.03) to obtain nebular emission from H II regions around young stars (Gutkin et al. 2016), AGN narrow-line regions (Feltre et al. 2016), and post-AGB stellar populations (Hirschmann et al. 2017). We used the same grids of emission-line models as those described in Hirschmann et al. (2019, 2023a), which represent updated versions of the ones detailed in Hirschmann et al. (2017). The general modelling approach remained the same.

2.2.1. Emission-line models for young stars, AGN narrow-line regions and post-AGB stellar populations

For each galaxy, the Gutkin et al. (2016) emission-line model grids for young star clusters (hereafter SF models) describe an ensemble of typical, ionisation-bounded H II regions illuminated by 10 Myr-old stellar populations with constant star formation history. The H II regions are characterised by various model parameters, such as the H II gas density, interstellar metallicity, ionisation parameter, dust-to-metal mass ratio, and C/O abundance (see Table 1 in Hirschmann et al. 2017). To model the stellar component, we used the most recent version of the Bruzual & Charlot (2003) stellar population synthesis model (Charlot & Bruzual in prep.) with a standard Chabrier (2003) initial mass function (IMF) truncated at 0.1 and 300 M_⊙. This version contains updated spectra of Wolf-Rayet stars and newer evolutionary tracks for post-AGB stars from Miller Bertolami (2016).

The photoionisation model grids for AGN narrow-line regions (hereafter AGN models) from Feltre et al. (2016) assume an emitted spectrum, following a broken power law, incident on gas clouds with uniform properties. Models on the grid are described by the interstellar metallicity, carbon-to-oxygen abundance, and dust-to-metal mass ratio in the narrow-line region, as well as the ionised gas density in the clouds and the ionisation parameter (see Table 1 in Hirschmann et al. 2017). We did not model broad-line regions, meaning we implicitly assumed that all AGNs are of Type 2 (see Sect. 8.3 for more details).

Lastly, the model grids for evolved post-AGB stellar populations (hereafter PAGB models) from Hirschmann et al. (2017) again use the updated version of the Bruzual & Charlot (2003) stellar population synthesis code as input for Cloudy, this time for evolved, single-age stellar populations between 3 and 9 Gyr at a range of stellar metallicities. The chosen ages represent the time span in which a population of post-AGB stars has built up and produces a significant amount of ionising photons. The models are largely parameterised the same way as the Gutkin et al. (2016) models, except for allowing the interstellar metallicity to differ from the stellar metallicity (see Table 1 in Hirschmann et al. 2017).

2.2.2. Coupling the photoionisation models to GAEA-LC

Connecting the SF, AGN, and PAGB models described in Sect. 2.2.1 to the GAEA-LC catalogues of Sect. 2.1 was done in a self-consistent way following the methodology from previous works using these models (Hirschmann et al. 2017, 2019, 2023a; Hirschmann et al. 2023b) but with slight modifications. In particular, our adjustments account for coupling the emission-line models to a semi-analytic model such as GAEA rather than, as in the proceeding works, to hydrodynamical simulations explicitly containing baryonic components.

For each galaxy in the light cone, a SF, AGN, and PAGB emission-line model was chosen according to which relevant model parameters match most closely the simulated galaxy properties available from GAEA. Model parameters, for which no equivalent property could be recovered from GAEA, were fixed to standard values. The dust-to-metal mass ratio ξ_d, for instance, was set to 0.3 for all galaxies.

The most suitable SF model was selected according to the parameters closest to the simulated GAEA values for the global interstellar metallicity and C/O abundance, as well as the ionisation parameter. The ionisation parameter is a measure for the degree of ionisation of the ISM and, thus, depends both on the hardness and intensity of the ionising radiation coming from the source, as well as the distribution and density of the gas. For H II regions ionised by young stars in GAEA-LC galaxies, we followed the computation of the SF ionisation parameter U_sim, ⋆ according to equations (1) and (2) in Hirschmann et al. (2017). The simulated SFR of the stellar population provides the rate of ionising photons Q_sim, ⋆, which are incident on hydrogen gas, characterised by the filling factor ϵ and the hydrogen gas density in ionised regions (n_H, ⋆, set to 10² cm⁻³). The filling factor is calibrated such that at redshift zero galaxies in GAEA reproduce the Carton et al. (2017) relation between U_sim, ⋆ and interstellar metallicity (i.e. log₁₀U ≈ −0.8log₁₀(Z_ISM/Z_⊙)−3.58). At higher redshift, the filling factor then evolves according to the global average gas density in galaxies from the cosmological simulation IllustrisTNG (same method as in Hirschmann et al. 2023a).

By analogy, AGN models for nebular emission from the narrow-line region were coupled to GAEA-LC galaxies by matching the central (as opposed to the global) interstellar metallicity and C/O abundance, as well as the simulated ionisation parameter U_sim, •. GAEA does not trace the central metallicity directly, thus we assume it to be twice the global value. Testing other values showed that our results are insensitive to this assumption. The rate of ionising photons is now set by the AGN luminosity of the simulated galaxy. Its spectrum is assumed to follow a broken power law with adjustable index α = −1.7 between wavelengths 0.001 μm and 0.25 μm (see Eq. (5) in Feltre et al. 2016). We set the density of ionised gas clumps in the narrow-line region regions n_H, • to 10³cm⁻³ and modelled the volume-filling factor according to the Carton et al. (2017) relation, now scaling it with the central average density in IllustrisTNG galaxies for increasing redshift.

For PAGB emission-line models, we computed the average age and metallicity of the stellar population provided by Bruzual & Charlot (2003) synthesis models, including only evolved stars older than 3 Gyr. We matched these values to available PAGB model grid ages and metallicities from Hirschmann et al. (2017) and found the rate of ionising photons based on the mass contained in the evolved stars. Then, as before, we selected the PAGB model with the closest global interstellar metallicity, global C/O ratio and ionisation parameter U_sim, ⋄.

2.3. Total emission-line luminosities and observer-like fluxes

After coupling the emission-line models to the GAEA-LC catalogues, we recovered the total emission-line luminosities for each galaxy by summing over the contributions from the matched SF, AGN, and PAGB models. In this study, we focussed on spectroscopic diagnostics based on seven optical emission lines: Hα, Hβ, [S II]λλ6717, 6731 (hereafter simply [S II]), [N II]λ6584 ([N II]), [O I]λ6300 ([O I]), [O III]λ5007 ([O III]), and [O II]λλ3727, 3729 ([O II]). For simplicity, we adopted the notation L_[O III]/L_Hβ = [O III]/Hβ for luminosity ratios. In order to make targeted predictions for the observability of line-emitting galaxies with Euclid, we computed observer-like fluxes based on the location and redshift of each galaxy in the light cone and then applied the EWS and EDS specific detection limits (f ≥ 2 × 10⁻¹⁶ and 6 × 10⁻¹⁷ erg s⁻¹ cm⁻², respectively) to our model fluxes.

2.4. Dust attenuation

Observed galaxies will be subject to non-negligible attenuation due to their dust content, which is usually estimated with the Balmer decrement (Kennicutt 1992). While Cloudy treats dust processes, including attenuation, self-consistently within H II regions, it does not account for interstellar dust. Estimating the exact contribution is challenging, particularly for higher-redshift galaxies, due to the limited observational data and redshift effects, which can obscure the signs of dust. Additionally, the complexity of dust properties, intrinsic variability among galaxies, and confounding factors such as star formation and AGN activity further complicate accurate estimations.

Locally found empirical relations (e.g. Garn & Best 2010; Zahid et al. 2013) have established that the overall dust attenuation broadly scales with the stellar mass, with sub-dominant effects from the SFR and the metallicity. For higher redshifts, results from a series of studies using various dust indicators (e.g. Sobral et al. 2013; Domínguez et al. 2012; Kashino et al. 2014; Price et al. 2014; Mclure et al. 2017; Cullen et al. 2017; Maheson et al. 2024; Shapley et al. 2022) have shown that, until at least redshift three, there is no significant evolution of the relationship between dust attenuation and stellar mass. Thus, we used the local Garn & Best (2010) relation to compute the V-band magnitude A_V for each galaxy and applied a line-of-sight Calzetti et al. (2000) attenuation to the predicted line fluxes.

In order to illustrate the potential impact of dust attenuation on the observability of various line-emitting galaxy populations in the upcoming Euclid surveys, we will mainly distinguish between three different types of samples:

• The intrinsic sample of different line-emitters as predicted by our GAEA-LC framework;

• Flux-limited galaxy populations, which for a chosen strong emission line contain only galaxies with fluxes exceeding the EWS or EDS flux limits;

• Dust-attenuated populations, which have the Calzetti et al. (2000) curve applied before enforcing the respective flux cuts.

In Sect. 8.1, we elaborate on the treatment of dust attenuation in the context of our results.

2.5. Instrumental and environmental effects

At this stage we did not account for additional instrumental and environmental effects which might limit the observation of emission lines. As a result, we excluded considerations of the astrophysical background, read-out and detector noise, as well as spectral resolution and recovery of blended Hα and [N II] emission lines. We elaborate on these points in Sect. 8.4. The strength of our framework is its self-consistent modelling of emission lines due to both young stars and AGNs, which allowed us to assess the intermediate redshift validity of locally calibrated spectroscopic diagnostics from a physical perspective. Estimates on the biasing of various scaling relations and the observability of different line-emitting galaxies should be understood as upper limits.

2.6. Distinguishing between dominant ionising sources in GAEA

As in Hirschmann et al. (2017, 2019, a),b), we used the theoretical BHAR/SFR criterion to divide the GAEA-LC sample according to their dominant ionising source, meaning SF-dominated, AGN-dominated, and composite galaxies, which contain significant SF and AGN contribution. For this study, we have adjusted the BHAR/SFR boundaries in order to ensure that the populations are reasonably separated in all diagrams:

• SF-dominated galaxies: BHAR/SFR < 10⁻³

• Composite galaxies: 10⁻³ < BHAR/SFR < 10^−2.2

• AGN-dominated galaxies: BHAR/SFR > 10^−2.2

Hirschmann et al. (2017, 2019, 2023a,b) further define galaxies to be dominated by aged PAGB stars if their Hβ emission exceeds the contribution from both young stars and AGNs. They generally have low star formation rates and form a sub-category of galaxies with low-ionisation (nuclear) emission-line regions (LIER/LINER, see Heckman 1980; Kauffmann et al. 2003; Singh et al. 2013; Belfiore et al. 2016). While these galaxies do exist in our GAEA-LC sample, they become increasingly rare beyond the local Universe where galaxies exhibit younger stellar populations and high SFRs. This agrees with Hirschmann et al. (2023a), who found that the number of PAGB-dominated galaxies in their sample of post-processed IllustrisTNG galaxies rapidly decreases from a few per cent below redshift one to a negligible fraction above redshift one. Additionally, with luminosities of order 10³⁹ erg s⁻¹, their emission is relatively faint, meaning they lie below the EDS flux limit already at around redshift 0.1. Thus, we conclude that Euclid will likely not observe any PAGB-dominated galaxies and exclude them from further analysis.

3. Validation of the method

While the photoionisation models discussed in Sect. 2.2.1 have been successfully applied to galaxies formed in numerical simulations such as SPHGal and the IllustrisTNG suite (Hirschmann et al. 2017, 2019, 2023a; Hirschmann et al. 2023b), this work represents the first instance of applying them to a semi-analytic model, such as GAEA. As semi-analytic models do not explicitly treat gas dynamics, we had to adapt the approach and underlying assumptions in order construct the mock-emission lines for GAEA-LC. Thus, we validated our method by comparing our emission-line predictions to observational data and theoretical predictions for low redshifts. Our self-consistent modelling then allows us to extend our predictions to higher redshifts.

3.1. Hα number counts and luminosity function

As a first validation, we compare in Fig. 1 predicted redshift distributions of the Hα emitter number density (blue lines) against models from Pozzetti et al. (2016), which have been derived from fits to a collection of Hα luminosity functions (Gallego et al. 1995; Ly et al. 2007; Tresse & Maddox 1998; Shioya et al. 2008; Sobral et al. 2013; Colbert et al. 2013; Tresse et al. 2002; Shim et al. 2009; Hopkins et al. 2000; Yan et al. 1999; Geach et al. 2008; Hayes et al. 2010). Model 1 (green) and 2 (black) use the same collection of survey data and underlying Schechter function, but differ in their implementation of the redshift evolution. Model 3 (red) was determined using only surveys covering higher redshifts from 0.7–2.23 (points outlined in red) and is based on a broken power law. While Model 1 and 2 produce similar number densities, Model 3’s prediction is generally lower by a factor of 1.5–2.5. This difference can be attributed to the large underlying scatter in the observed luminosity functions, as well as the different functional form. We include the scatter of number count predictions from the collection of integrated luminosity functions, as well as direct estimates of cumulative number counts derived from the WFC3 Infrared Spectroscopic Parallels (WISP) survey (Colbert et al. 2013; Mehta et al. 2015, circular data points). The significant spread of data points (highlighted in yellow shaded area) is a result of varying survey set-ups, such as using different instruments with different selection functions, collecting either narrow-band or spectroscopic measurements, as well as varying areal coverage and treatments of cosmic variance.

Fig. 1. Redshift distribution of the number density of Hα emitters above flux thresholds 2 × 10−16 erg s−1 cm−2 (top panel) and 5 × 10−17 erg s−1 cm−2 (bottom panel), corresponding to the respective EWS and EDS flux limits. Predictions from the GAEA-LC framework (blue lines) use the intrinsic population of emitters (solid) and dust-attenuated versions with a flat Av (dashed) and a mass-dependent Av scaling (dotted Garn & Best 2010). Models 1 (green), 2 (black), and 3 (red) from Pozzetti et al. (2016) represent various model fits to collections of uncorrected Hα survey results (data points, covering yellow-shaded area) across different redshift ranges. Model 3 has been fit to data points outlined in red.

Since the Pozzetti et al. (2016) models were constructed to explore observational yields, survey results have not been corrected for extinction. For a fair comparison, they should thus be contrasted with a dust-attenuated sample of Hα-emitters predicted by our GAEA-LC framework. However, the prevalence and nature of interstellar dust beyond redshift one is largely uncertain. In order to visualise the impact different dust distributions might have, we thus present three different predictions: the intrinsic GAEA-LC version (solid blue lines), a Calzetti et al. (2000) attenuation with a flat A_v scaling of 1 (dashed lines) and a mass-dependent A_v scaling of Garn & Best (2010, dotted lines), which is adopted in the remaining paper.

We compare the predicted number counts per deg² above two flux thresholds. For an EWS-like threshold of 2 × 10⁻¹⁶ erg s⁻¹ cm⁻² (top panel), the intrinsic GAEA-LC prediction differs by a factor of two from Models 1 and 2, while the prediction using a flat A_v scaling agrees almost exactly with Models 1 and 2. The curve using a Garn & Best (2010)A_v scaling gives the lowest prediction of the number count density, which is between Model 2 and 3 until it diverges from the observational scatter around redshift 1.8, when Hα stops being observable with Euclid. This steeper decrease with redshift compared to other predictions is because, at high redshift, Hα emitters above the flux cut are generally massive and thus more strongly attenuated according to the Garn & Best (2010) scaling. Under the EDS-like flux cut (5 × 10⁻¹⁷ erg s⁻¹ cm⁻²), number densities from the intrinsic model are in good agreement with Models 1 and 2. Predictions from the dust-attenuated GAEA-LC samples are in better agreement with Model 3, but decrease below the observational scatter at redshift 1.5.

Overall, predictions for the redshift distribution of the Hα emitter number density from our GAEA-LC models fall within the significant scatter of observational estimates, and broadly agree with the models from Pozzetti et al. (2016) for both the EWS- and EDS-like thresholds. Varying the applied dust attenuation significantly changes the prediction beyond redshift 0.5. Which dust attenuation model agrees best with which Pozzetti et al. (2016) model and the integrated observational estimates varies if considering the EWS or the EDS threshold. Our line-of-sight Calzetti et al. (2000) extinction likely does not capture the full complexity of the nature of dust in this regime. Thus, exact comparisons should be approached with caution. While the prediction using the A_v scaling from Garn & Best (2010) drops below the observational scatter at 1.5 for the EDS and 1.8 for the EWS, this represents a range with relatively few data points, which have mostly been determined from narrow-band measurements (Sobral et al. 2013; Hayes et al. 2010; Geach et al. 2008) that often suffer from contamination due to other emission lines. Moreover, this is at the edge of the detectable range with Euclid, meaning this divergence from the scatter will not affect our predictions significantly. Since the Garn & Best (2010) scaling is empirically motivated and has not shown significant evolution across the relevant range (see Sects. 2.4 and 8.1), we thus continue to adopt it for the remaining paper. For this sample, the average number density from redshift 0.9–1.8 above an EWS-like threshold is 7800 Hα emitters/deg². Euclid’s mission requirements specify the redshift measurement of 1700 Hα emitters deg⁻². Thus, according to this estimate, Euclid would only have to recover redshifts for around 21% of these emitters.

We further compare our model predictions against observed Hα luminosity functions. Figure 2 shows the evolution of the intrinsic Hα luminosity function ϕ between redshifts 0 and 2.2 for the full population of GAEA-LC galaxies (solid blue lines) and the sub-group of SF-dominated galaxies (dotted blue lines). As established in Sect. 1, the Hα luminosity is often used as a proxy for star formation and thus, on cosmological scales, the evolution of the Hα luminosity function traces the cosmic SFR density. The full GAEA-LC and the SF-dominated sample appear closely spaced in the figure, indicating that the Hα luminosity is indeed shaped by the emission from H II regions around young stars.

Fig. 2. Redshift evolution of the Hα luminosity function for all (solid blue lines) and only SF-dominated (dotted blue lines) GAEA-LC galaxy populations from redshift zero to 2.2. Overplotted are various fits to dust-corrected observational results (thin lines) from narrow-band and spectroscopic surveys within the redshift range indicated by the legend in each panel (Tresse et al. 2002; Fujita et al. 2003; Perez-Gonzalez et al. 2003; Ly et al. 2007; Villar et al. 2007; Shim et al. 2009; Hayes et al. 2010; Sobral et al. 2013). Due to selection effects, the resulting Schechter fits rely on measurements covering only part of the luminosity range (solid), which are then extrapolated to low and high luminosities (dashed). For reference, the EWS (red) and EDS (blue) detection limits across the redshift range between 0.7 and 0.8 are shown in shaded areas.

Overall, cosmic star formation has sharply declined from redshift zero to 2, and as a result, the luminosity function exhibits a similarly strong evolution. The maximum at the luminous end decreases by around 0.8 dex from L_Hα ∼ 10^43.8 erg s⁻¹ to L_Hα ∼ 10⁴³ erg s⁻¹. This is consistent with dust-corrected observational fits from both narrow-band and spectroscopic surveys (e.g. Tresse et al. 2002; Fujita et al. 2003; Perez-Gonzalez et al. 2003; Ly et al. 2007; Villar et al. 2007; Shim et al. 2009; Hayes et al. 2010; Sobral et al. 2013, shown in thin lines, redshifts indicated in parentheses). We distinguish between the luminosity ranges constrained by measurements (solid) and the extrapolated fits to low and high luminosities (dashed).

At high redshift, between 2 and 2.2, our predictions for the luminous end of the luminosity function are in excellent agreement with observational results. Below L_Hα ∼ 10⁴³ erg s⁻¹, our prediction slowly starts diverging towards lower values, until it reaches a turnover around L_Hα ∼ 10^41.8 erg s⁻¹. This feature is an effect resulting from GAEA’s resolution limit at M_⋆ ∼ 10⁹ M_⊙, as well as the applied H-band magnitude cut, which particularly affects low-luminosity galaxies at higher redshift. As a result, the faint end of the predicted Hα luminosity function is underestimated. In addition, we note that observational surveys either undersample the faint end, in which case they apply estimated completeness corrections, or do not measure any low-luminosity Hα emitters and only extrapolate the Schechter fit. Thus, at the low-luminosity end, our predictions should be compared to observational results with caution.

Across the 0.7–0.8 redshift range, the GAEA-LC prediction lies among the spread of the different survey results. These exhibit varying shapes and a large scatter in ϕ of around 0.8–1 dex. In the case of Shim et al. (2009) and Tresse et al. (2002), the resulting luminosity functions have been fit to data covering relatively large redshift ranges, averaging out any potential evolution across them. However, even Sobral et al. (2013) and Villar et al. (2007), which both cover redshifts around 0.8, exhibit very different shapes. In general, we conclude that ϕ is poorly constrained in this redshift regime, also illustrating the need for more extensive spectroscopic surveys.

In the 0–0.1 redshift range, we slightly overpredict the luminous end with respect to observational determinations. At luminosities below L_Hα ∼ 10⁴² erg s⁻¹, our result lies within the large scatter among them. As for redshifts in 0.7–0.8, the surveys have targeted different redshift ranges, which could partially explain this scatter. However, the luminosity function at redshift 0.24 determined by Fujita et al. (2003) lies between the results at redshift 0.08 and 0.026 from Ly et al. (2007) and Perez-Gonzalez et al. (2003), which differ by more than 1 dex at the faint end. This large discrepancy suggests that ϕ is not well-constrained in the low-redshift regime either. Lastly, due to the geometry of the GAEA light cone, our sample contains limited number counts at low redshift, which makes our estimate susceptible to low number statistics. However, we overall conclude that our framework predicts the evolution of the Hα luminosity function in broad agreement with empirical results and any discrepancies lie within observational and modelling uncertainties.

3.2. Distinguishing between ionising sources using BPT diagnostic diagrams

In Fig. 3, we show the locations of the predicted SF-dominated (left column), composite (middle column), and AGN-dominated (right column) galaxy populations at redshift less than 0.3 in the standard BPT diagrams: [O III]/Hβ against [N II]/Hα (top row, Baldwin et al. 1981), [O III]/Hβ against [S II]/Hα (middle row), and [O III]/Hβ against [O I]/Hα (bottom row, both Veilleux & Osterbrock 1987). In addition to the intrinsic GAEA-LC samples (grey contour lines), we contrast the observed sample from the Sloan Digital Sky Survey (SDSS, filled grey contours in background) with simulated SDSS-like galaxy populations (blue contours for SF-dominated, red contours for composites, green data points for AGN-dominated). We note that the observed SDSS galaxies shown here are not separated according to type and instead represent the combined sample. As in Hirschmann et al. (2017, 2019, 2023a,b), SDSS-like galaxies were selected by applying a flux limit of 5 × 10⁻¹⁷ erg s⁻¹ cm⁻² (see Table 1 in Juneau et al. 2014) to all four lines. We note that, in some instances, the grid parameterisation of the photoionisation models results in visible accumulations of galaxies at discrete points in the diagrams, such as for SF-dominated galaxies in [S II]/Hα and [O I]/Hα plots and AGN-dominated galaxies in all plots. For the latter populations, only a few hundreds of galaxies exceed the SDSS flux cut and thus, we show the scatter of the individual data points. In general, simulated SDSS galaxies occupy the same region as the observed SDSS galaxies. The flux limit mostly cuts out galaxies with particularly low [N II]/Hα, [S II]/Hα and [O I]/Hα. In the photoionisation models by Gutkin et al. (2016) these represent the lowest metallicity galaxies, in line with our expectation that these tend to be low-luminosity (see Tremonti et al. 2004).

Fig. 3. Location of GAEA-LC galaxy populations at redshift 0–0.3 in the classical BPT diagrams, [O III]/Hβ versus [N II]/Hα (top row), [O III]/Hβ versus [S II]/Hα (middle row), and [O III]/Hβ versus [O I]/Hα (bottom row). Shown are simulated SDSS-like galaxies (limited to fluxes above 5 × 10−17 erg s−1 cm−2, coloured contours and data points), alongside the intrinsic GAEA-LC sample (grey contour lines). For comparison, SDSS-observed galaxies are plotted in the background (filled grey contours). Galaxy populations are divided according to dominant ionising source, meaning SF-dominated galaxies (blue, left column), composite galaxies (red, middle column), and AGN-dominated galaxies (green, right column). Overplotted are empirical selection criteria meant to broadly distinguish SF galaxies (below dashed lines, Kewley et al. 2001 in top row, Kauffmann et al. 2003 in middle and bottom row) and active galaxies (above dashed lines). An additional criterion separates composite galaxies (above dashed line and below dotted line, Kewley et al. 2001) from purely AGN-dominated galaxies in the [O III]/Hβ versus [N II]/Hα diagram. In all diagrams, LI(N)ER are expected to fall in the bottom right corner (rectangle defined by dash-dotted lines in top row, Kauffmann et al. 2003 and area below dash-dotted lines in middle and bottom row Kewley et al. 2006).

In order to test our framework, we then compare our GAEA-LC samples to optical criteria used to distinguish between dominant ionising sources in local galaxies. By combining photoionisation and stellar population synthesis models, Kewley et al. (2001) set a theoretical upper limit to the location of star-forming galaxies above which the emission from galaxies would not be explainable without a strong AGN component (dotted line in top row, dashed lines in middle and bottom row). Based on SDSS observations of nearby AGNs, Kauffmann et al. (2003) found that in their sample, galaxies containing an AGN are confined above the dashed line in the [N II]/Hα diagram. Thus, the area between the dotted and dashed lines can be understood as a transition region where we expect composite galaxies to be located.

As expected, at redshift less than 0.3, the optical selection criteria for BPT diagrams are successful at separating the sample according to ionising source. The majority of GAEA-LC galaxies selected according to the theoretical BHAR/SFR criteria are confined to the predicted locations for SF-dominated, composite, and AGN-dominated populations. Composite galaxies partially extend into the AGN-dominated regime in the [N II]/Hα diagram, but, nevertheless, occupy a distinct region. There are no selection criteria to identify composite galaxies in the [S II]/Hα diagram and [O I]/Hα diagram and they mostly overlap with the SF-dominated population. In general, these diagrams also provide a less clear distinction of SF- and AGN-dominated galaxies compared to the [N II]/Hα diagram. As a result, we subsequently focus on this diagram when we extend our predictions to intermediate redshifts in Sect. 5.

3.3. Evolution of the star-forming branch in the [O III]/Hβ-versus-[N II]/Hα diagram

In this section, we verify the observed increase from low to high redshift of the [O III]/Hβ ratio at fixed [N II]/Hα ratio for SF galaxies. This was initially observed around redshift two (e.g. Shapley et al. 2004, 2015; Liu et al. 2008; Hainline et al. 2009; Bian et al. 2010; Lehnert et al. 2009; Yabe et al. 2012; Masters et al. 2014; Steidel et al. 2014; Kashino et al. 2017; Strom et al. 2017), while recent JWST/NIRSpec data showed a continuation of this trend to redshift five (e.g. Cameron et al. 2023; Sanders et al. 2023).

From our GAEA-LC sample, we selected SF-dominated galaxies with resolved stellar masses (M_⋆ ≥ 10⁹ M_⊙) and SDSS-like fluxes (≥5 × 10⁻¹⁷ erg s⁻¹ cm⁻²), for which we then computed the average [O III]/Hβ at fixed [N II]/Hα. The result is shown in Fig. 4 for different redshift bins between redshift zero and 3.9 (thin coloured lines, indicated in legend), alongside the optical selection criteria (Kewley et al. 2001; Kauffmann et al. 2003; Kewley et al. 2006, black lines, as in Fig. 3) and SDSS-observed galaxies (grey contours in background). At fixed [N II]/Hα, the mean [O III]/Hβ ratio increases with redshift. According to detailed investigations by Hirschmann et al. (2017, 2023a), the redshift evolution of the [O III]/Hβ at fixed [N II]/Hα is driven by the elevated SFR and global gas density at higher redshift, increasing the ionisation parameter U_⋆ and, as a result, the probability for doubly ionised oxygen. Comparing the predicted average relation at redshift 1.8–2.5 with the fit of 219 observed SF galaxies at redshift 2.3 from Steidel et al. (2014, thick grey line), we find good agreement. We note that the average relations from our simulated galaxies appear slightly steeper than the Steidel et al. (2014) determination, which can be partly explained by the difficultly of correctly identifying composites observationally and the slight dependence on the choice of BHAR/SFR boundaries used in our theoretical definition.

Fig. 4. Evolution of average [O III]/Hβ versus [N II]/Hα of SF-dominated GAEA-LC galaxies at different redshift intervals, as indicated by the legend (coloured lines). Shown are simulated SDSS-like galaxies with masses greater than 109 M⊙. Overplotted is the mean relation found by Steidel et al. (2014) in their sample of star-forming galaxies at redshift 2.3 (thick grey line) and, as in Fig. 3, the SDSS sample (grey contours) and empirical criteria distinguishing between ionising sources (Kewley et al. 2001; Kauffmann et al. 2003, black lines).

3.4. Relation of strong line luminosities to SFRs and AGN luminosities at low redshift

As a last test of our methodology, we demonstrate that its low-redshift predictions reproduce local calibrations between strong line luminosities and galaxy properties, such as the SFR and the AGN luminosity. In the first two panels of Fig. 5, we show the average Hα (L_Hα, top panel) and [O III]+Hβ (L_{[O III]+Hβ}, middle panel) luminosity at fixed SFR for GAEA-LC galaxies between redshift 0 and 0.3 (blue line with shaded area indicating one standard deviation). These proxies are only expected to be robust for SF galaxies (also see Fig. 10 in Hirschmann et al. 2023a). We note the distinction between SF galaxies, which are commonly defined as having high specific SFRs (sSFR  ≡ SFR/M_⋆), and SF-dominated galaxies, in which the ionisation budget due to young star clusters is greater than the contribution from other sources. While most SF galaxies are SF-dominated, SF-dominated galaxies are not necessarily highly star-forming in our simulations. However, we can generally consider our predictions for SF-dominated galaxies to be a good proxy for SF galaxies.

Fig. 5. Average Hα (first panel, blue line) and [O III] + Hβ (second panel, blue line) line luminosities versus the SFR for GAEA-LC galaxies with specific SFR > 10−10.5 yr−1 in the redshift range 0–0.3, alongside the one standard deviation scatter (shaded area). Overplotted, for comparison, are widely used relations from Kennicutt & Evans (2012, for Hα, dashed line) and Osterbrock & Ferland (2006, for [O III] + Hβ, dotted line). The bottom panel shows the average [O III] line luminosity versus AGN luminosity for active GAEA-LC galaxies with a 1σ scatter (green line and shaded area), plotted alongside the relation found by Lamastra et al. (2009, dash-dotted line) from a sample of 61 type-2 AGNs with z < 0.83.

As here we directly compare with observationally used relations, we select SF GAEA-LC galaxies by applying a sSFR cut of 10^−10.5 yr⁻¹, with no additional flux or mass cuts. Selecting SF-dominated galaxies produced an identical relation. Alongside our predictions, we plot local calibrations for the respective relationships. The Kennicutt & Evans (2012, top panel, dashed line) relation, originally published in Murphy et al. (2011), was derived from evolutionary synthesis models for SF galaxies based on a Kroupa et al. (2003) IMF, which yields nearly identical results to a Chabrier (2003) IMF (see Chomiuk & Povich 2011). However, our models cover a range of metallicities from 10⁻³ Z_⊙ to 1 Z_⊙, while Murphy et al. (2011) assume solar metallicity Z_⊙. Thus, we do not expect our predictions to match these relations exactly. Across all SFR, our Hα luminosity-SFR relation exhibits the same slope as the Kennicutt & Evans (2012) calibration, but is slightly offset towards higher L_Hα values, which can be explained by the slightly different modelling assumptions. This difference increases at low SFR, which is where the magnitude cut introduces a bias towards an increased L_Hα in the remaining galaxies. Our L_{[O III]+Hβ}–SFR prediction agrees well with Osterbrock & Ferland (2006, middle panel, dotted line), especially at log₁₀(SFR/M_⊙yr⁻¹) > −0.5, below which it diverges for similar reasons as detailed above.

In the bottom panel of Fig. 5, we explore the relationship between the [O III] luminosity L_[O III] and the bolometric AGN luminosity L_AGN. For active (meaning composite and AGN-dominated) galaxies at redshift 0–0.3, we show the predicted mean [O III] luminosity at fixed AGN luminosity (green line with shaded one standard deviation). We compare to an empirical relation by Lamastra et al. (2009), who, based on a sample of 61 type-2 AGNs with redshift less than 0.83, found a luminosity-dependent [O III]-bolometric correction factor in the ranges log₁₀(L_[O III]) = 38–40, 40–42 and 42–44 (dash-dotted lines). In general, we note excellent agreement with our GAEA-LC relation across the entire L_AGN range.

4. Exploring Euclid’s selection bias

Section 3 has demonstrated that our model framework, which connects emission-line models to simulated galaxies from the GAEA semi-analytic model, successfully reproduces a wide range of locally observed emission-line properties. As a result, we have established a self-consistent, physically validated sample of line-emitting galaxies between redshift 0–3.9, which we can now use to put Euclid forecasts on solid ground. In this Section, we explore the selection bias resulting from the observation of the BPT emission lines (Hα, Hβ, [O III], and [N II]) by considering the EWS and EDS specific flux limits at the relevant redshifts. Specifically, we assess the effects on the line flux-stellar mass planes, standard scaling relations, the prevalence of luminous AGNs, and the observability of line-emitting populations. We use ‘detectability’ to describe the redshift range in which an emission line of a given rest wavelength falls into the wavelength sensitivity range of Euclid’s grisms. For a given redshift, we then define ‘observability’ as the number of galaxies emitting line intensities within the detectable range and above the flux limits, respectively for the EWS and EDS. We note that this definition is distinct from the Euclid mission requirement of ‘completeness’ (see Racca et al. 2016), which is defined as the fraction of galaxies at redshift 0.9–1.8 emitting fluxes above the limit 2 × 10⁻¹⁶ erg s⁻¹ cm⁻², for which a redshift measurement can be recovered.

4.1. Evolution of line detectability with the blue and red grisms

Before imparting on our analysis, we first illustrate at which redshifts key optical emission lines are detectable with Euclid’s red and blue grisms, given their respective sensitivity ranges. Figure 6 shows the spectral coverage of strong emission lines [S II] (crimson), [N II] (red), Hα (yellow), [O I] (green), [O III] (cyan), Hβ (blue), and [O II] (purple) as function of redshift and where they overlap with the RGS (red area) and BGS (blue area) rest frame sensitivity ranges.

Fig. 6. Redshifted wavelengths of strong emission lines in the rest-frame optical (crimson: [S II]λ6724, red: [N II]λ6584, yellow: Hα, green: [O I]λ6300, cyan: [O III]λ5007, blue: Hβ, and purple: [O II]λ3727) and their resulting detectability using Euclid’s blue grism (BGS, blue shaded area) and red grism (RGS, red shaded area). The redshift bins used for the following analysis (vertical dashed lines) were chosen according to the overlap of observed wavelengths with the BGS and RGS sensitivity ranges, 0.93–1.37 μm and 1.21–1.89 μm, respectively.

[S II], [N II], Hα, and [O I] are all roughly detectable from redshift 0.4–0.9 in the BGS and then fall into the sensitivity range of the RGS until roughly redshift 1.8. Above redshift 1.8, they are not be detected by Euclid. [O III] and Hβ, on the other hand, are detectable with the BGS at redshift 0.9–1.5, then with the RGS from 1.5 to 2.5. [O II] is detectable in the BGS from 1.5 to 2.5, while the RGS dectectability extends from 2.3 to 4.1, largely outside the range considered here. Based on the overlap of these regimes, we chose the redshift bins 0.4–0.9, 0.9–1.2, 1.2–1.5, 1.5–1.8, and 1.8–2.5 (dashed vertical lines) as a basis for the following sections. Depending on the figure, we either focus on bins in which the chosen lines are detectable with Euclid, sometimes combining multiple bins into one, or show predictions for all five bins across redshift 0.4–2.5, in order to demonstrate the physical evolution of the underlying emission-line properties, regardless of detectability with Euclid.

4.2. Fluxes and observability of line-emitting galaxies according to their redshift and stellar mass

In Fig. 7, we visualise the location of GAEA-LC galaxies in the Hα (first row), [N II] (second row), Hβ (third row), and [O III] (fourth row) flux versus stellar mass plane. We note that the Hα and [N II] lines will be blended in Euclid data, however the line fluxes can be recovered using a simultaneous 3-Gaussian fit, as described in Sect. 8.4.

Fig. 7. Location of intrinsic (grey contour lines) and dust-attenuated, but not flux-limited (blue contours) GAEA-LC galaxy populations in the Hα (top row), [N II] (second row), Hβ (third row), and [O III] (bottom row) line flux-stellar mass plane at different redshift ranges, following their observability with Euclid’s grisms given the line’s wavelength (different columns as indicated by the legend). Overplotted are the flux limits of the blue and red grisms in the Deep Survey mode (EDS, blue and red dashed lines), and the red grism in the Wide Survey mode (EWS, yellow dash-dotted lines). Grey panels mark GAEA’s resolution limit of stellar masses below 109 M⊙. Percentages indicate the observable fractions of the resolved sample above 109 M⊙ when applying the EWS or EDS limit to the dust-attenuated emission-line fluxes (unattenuated fluxes in parentheses, colours reflecting survey mode and grism).

We contrast the intrinsic GAEA-LC sample (grey contour lines) with the dust-attenuated, but not flux-limited sample (blue contours). We show the line-emitting populations across the redshift bins for which the line is detectable (as indicated at the top of each column and the middle of the right column). The respective grism flux limits in the survey configurations are indicated, noting again that the blue grism (dashed blue lines) is only deployed in the EDS, while the red grism is active in both the EWS (yellow dash-dotted lines) and the EDS (red dashed lines). In the upper left corners, we indicate the observable percentage of dust-attenuated galaxies, according to the survey’s grism flux threshold (colours corresponding to the configuration in the legend). Percentages for only flux-limited and unattenuated populations are shown in parentheses. We excluded galaxies below the conservative GAEA mass resolution limit of 10⁹ M_⊙ (grey shaded area in plot) from this calculation.

For all emission lines, we note a clear stellar-mass dependence of the flux, with the most massive galaxies emitting the highest fluxes. This is due to higher SFRs in more massive galaxies, which increase the number of ionising photons. As a result, the observable galaxy populations will be dominated by massive galaxies. Our estimate of dust attenuation shifts the GAEA-LC sample towards lower fluxes, reducing the number of observable galaxies. Due to the stellar mass-dependence of the scaling law, this effect is more pronounced at higher masses, reducing fluxes by up to 0.6 dex. Overall, this reduces the percentage of observable fluxes by 15–20%. In the following, we refer to the observable percentages of the dust-attenuated line-emitting populations only.

Due to the distance dependence of line fluxes, observable percentages decrease significantly with increasing redshifts. While at redshift 0.4–0.9 we predict roughly 56% of the dust-attenuated Hα-emitting population to lie above the EDS limit, this is reduced to 29% at redshift 0.9–1.5 and 15% at redshift 1.5–1.8. In the EWS, where Hα is only detectable at redshift 0.9–1.8, these numbers are even further decreased to 7% and 2%, respectively. Hα is the brightest line with the highest observable percentages, but the other three lines generally follow a similar evolution. [N II] is already detectable in the EDS configuration at redshift 0.4–0.9 with a predicted observability of around 25%. It falls steeply to 6% and then 1.6% at redshift 0.9–1.5 and 1.5–1.8. In the EWS mode, observability is even lower in these regimes, with 0.8% at redshift 0.9–1.5 and 0.2% at 1.5–1.8. This steep redshift decline of the observability of [N II]-emitters is likely not only caused by the physical reduction in flux due to the luminosity distance, but also by the intrinsically lower prevalence of strong [N II]-emission at higher redshift compared to low redshift. Hirschmann et al. (2017) explained the increased [N II]/Hα ratio with decreasing redshift by inferring a higher gas-phase and stellar metallicity, as well as a decreased (s)SFR. High gas-phase metallicity boosts secondary production of nitrogen, while high stellar metallicity results in softer ionising radiation, which makes singly ionised nitrogen (N⁺) more likely than multiply ionised nitrogen. Similarly, a lower (s)SFR implies a lower ionisation parameter U_⋆, which further favours N⁺ over higher ionisation states.

In contrast, the observability of [O III] declines less steeply with redshift, from 17% to 10% between redshift 0.9–1.5 and 1.5–2.5, in agreement with the expectation of a greater number of bright [O III]-emitters at high redshift. Due to the rising metallicity towards low redshift, the electron temperature in the gas decreases, which disfavours collisional excitations into the [O III] state (Gutkin et al. 2016; Hirschmann et al. 2017). In analogy with [N II], the softer ionising radiation of metal-rich stars and lower SFR in the recent Universe favours [O II] at the expense of [O III]. In the given ranges for Hβ, we predict 0.3–3.3% above the flux limit in the EDS and basically none, except a few extreme objects, in the EWS. For all lines, we expect the full intrinsic populations to only be observable for massive galaxies with log₁₀(M_⋆/M_⊙) above 11.5–11.7 in the EDS and above 11.7–11.9 in the EWS, depending on the line and redshift range.

At redshift 0.4–0.9, a significant portion of galaxies lie in the unresolved mass regime, while, as a result of the magnitude cut, the light cone contains fewer unresolved low-mass galaxies at high redshift. However, for all lines and at all redshifts, galaxies above the EDS flux boundaries are at most 1.5% unresolved, such that we do not expect GAEA’s resolution limit to affect our predictions.

Lastly, we point towards the population in the lower right corner of each panel, which is separate from the main population and appears to grow between redshift 1.8 and 0.9. This region contains the PAGB-dominated galaxies, which, with fluxes of 10⁻¹⁸–10⁻²⁰ erg s⁻¹ cm⁻², lie far below both the EWS and EDS limits.

In conclusion, we predict that the observable populations of Hα, [N II], [O III] and Hβ-emitters are perhaps unsurprisingly biased towards the brightest and most massive galaxies, a trend which becomes more pronounced with increasing redshift. We expect Euclid to observe around 56% of Hα-emitters at redshift 0.4–0.9, and fewer than 30% for the other three lines and covered redshift ranges, given our underlying galaxy sample with a mass resolution limit of 10⁹ M_⊙ and magnitude cut of m_H = 25.

4.3. Biasing of standard scaling relations when observing line-emitters in the Euclid Wide and Deep Surveys

To further explore the selection bias imposed by surveying different emission lines, we examine how they trace various standard scaling relations in Fig. 8. We test two scenarios; surveying Hα emitters in the EWS mode (top row, hereafter EWS-Hα) and galaxies with simultaneous Hα, Hβ, [O III], and [N II] observability in the EDS (bottom row, hereafter EDS-BPT). The EDS-BPT configuration represents the sample for which the standard [O III]/Hβ versus [N II]/Hα BPT diagrams can potentially be used to determine the dominant ionising source. In both scenarios, the relevant lines are detectable between roughly redshift 0.9–1.8. We note that, in the EWS, the BPT lines will only be simultaneously detectable within the narrow redshift range 1.5–1.8.

Fig. 8. Scaling relations between sSFR and M⋆ (left column), M⋆ and Mhalo (middle column), and O/H abundance versus M⋆ (right column) for intrinsic (grey contour lines), flux-limited (light blue area), and dust-attenuated (blue contours) GAEA-LC populations in the redshift range 0.9–1.8. Shown are galaxies emitting Hα fluxes above the EWS cut (top row) and Hα, Hβ, [O III], and [N II] fluxes above the EDS cut (bottom row). The redshift range reflects the entire detectable range for the chosen lines in the given survey configuration. Percentages in the left panels indicate the observable percentage of the flux-limited sample (light blue) and dust-attenuated sample (dark blue). Each panel also contains the mean of the intrinsic (solid lines), flux-limited (dashed), and dust-attenuated (dotted) GAEA-LC populations. Grey panels indicate GAEA’s resolution limits in stellar mass (M⋆ < 109 M⊙) and halo mass (Mhalo < 1011 M⊙).

For the sSFR − M_⋆ (left column), M_⋆ − M_halo (middle column), and O/H abundance − M_⋆ (right column) relations, we show the intrinsic (grey contour lines), flux-limited (light blue area), and the dust-attenuated and flux-limited (blue contours, hereafter just called dust-attenuated) GAEA-LC galaxy populations. For ease of comparison, we plot the mean relations for the intrinsic GAEA-LC sample (solid line), flux-limited (dashed line), and dust-attenuated (dotted line) populations. Additionally, we indicate the resolution limit in stellar masses, M_⋆ < 10⁹ M_⊙, and in halo masses, M_halo < 10¹¹ M_⊙ (grey shaded areas). In percentages in the left-most panels, we show the respective observability for resolved galaxies in the flux-limited (light blue) and dust-attenuated (dark blue) sample.

For the flux-limited populations, the predicted sSFR − M_⋆ relations are biased towards highly star-forming and massive galaxies. However, a small quiescent population of galaxies with M_⋆ > 10¹¹ M_⊙ is present in both survey configurations. Applying the Calzetti et al. (2000) law biases the remaining sample to even higher stellar masses and, in the EDS-BPT configuration, eliminates the quiescent population. The average sSFR − M_⋆ relation is increased by 0.5–0.8 dex at fixed M_⋆ for both configurations when compared to the average of the intrinsic population. In general, we expect the EDS-BPT configuration to recover a smaller percentage of galaxies. For the EWS-Hα configuration the dust-attenuated (flux-limited) sample represents around 5.3% (18%) of intrinsic line-emitters, compared to 1.4% (11.8%) for the dust-attenuated (flux-limited) populations in the EDS-BPT configuration.

For the M_⋆ − M_halo relation, we show the distribution of stellar masses for all central and satellite galaxies against the mass of their sub-halos. Enforcing only the flux limit, without accounting for dust, preferentially excludes galaxies with low halo and stellar masses resulting in similar distributions for both survey configurations. Dust-attenuated populations exhibit further restricted distributions, which, for the EDS-BPT configuration, is more closely constrained around the average relation and includes gaps around M_halo = 10¹¹ and 10¹³ M_⊙. The population between M_halo = 10¹¹ and 10¹² M_⊙ is below the resolution limit and mostly represents satellites in subhalos. Therefore, the neighbouring gap to the main population is likely not of physical origin. The gap at M_halo = 10¹³ M_⊙ is a product of small number statistics, as our sample contains only few halos of such high masses. For both survey configurations, the average relation in the flux-limited sample shows an increased M_⋆ at low and high M_halo compared to the entire GAEA-LC sample. We note that the mean relations appear unbiased around 10¹² M_⊙, the turnover of the M_⋆ − M_halo relation.

Lastly, due to the correlation of mass and metallicity, the strong bias in stellar mass introduces a bias towards high metallicties in the O/H abundance − M_⋆ relation. For the flux-limited sample, the EWS-Hα configuration recovers some metal-poor galaxies below log₁₀(O/H)+12 ∼ 8, whereas the EDS-BPT configuration is expected to only survey metal-rich galaxies. In the dust-attenuated samples, no metal-poor galaxies remain after the flux cut for both surveys, confining the distribution between log₁₀(O/H)+12 ∼ 8 and 9. However, both biased mean relations are up to 0.15 dex below the mean for the full sample, as the flux limit also excludes some of the most metal-rich galaxies.

In conclusion, we expect the EWS-Hα survey between redshift 0.9–1.8 to recover 5.3–18% of the intrinsic population, depending on dust attenuation, while BPT lines should be recoverable for 1.4–11.8% of galaxies in the EDS. In both observing scenarios, the resulting samples will be biased towards massive, highly star-forming, and metal-rich systems.

We further explore the prevalence of luminous AGNs in the two survey configurations by examining the AGN luminosity L_AGN versus black hole mass M_BH relation at redshift 0.9–1.8 in Fig. 9. To ensure a physically robust sample of AGNs, we limit the selection to galaxies containing a black hole with a mass greater than the seed mass (M_BH > 10⁵ M_⊙) and a significant luminosity output (L_AGN > 10³⁹ erg s⁻¹). The general layout follows the same as in Fig. 8; we plot the intrinsic (grey contours), flux-limited (light blue area), and dust-attenuated (blue contours) populations for the EWS-Hα (top panel) and the EDS-BPT (bottom panel) survey configurations. Percentages in the bottom right corner of both panels mark the observability in corresponding colours. As we expect the upper end of L_AGN − M_BH to be set by the most powerful AGNs, we additionally show the flux-limited population for active galaxies only (orange contour lines) and the corresponding observable percentage compared to the intrinsic active population (light orange). The distribution of the dust-attenuated active population (not explicitly indicated here) constitutes the overlap of the full dust-attenuated sample and the flux-limited active sample. In the corner, we include the predicted observable percentage of this population (dark orange). We compare the results to relations expected for black holes accreting at fractions of the Eddington luminosity; 1 L_edd (solid black line), 0.1 L_edd (dashed line), and 0.01 L_edd (dash-dotted line). We estimate that the GAEAL_AGN convergence limit starts affecting our predictions below 10⁴³  erg s⁻¹ (grey shaded area, see also Fig. 2 in Fontanot et al. 2020).

Fig. 9. Scaling relation between the LAGN and black hole mass MBH for GAEA-LC populations in the redshift range 0.9–1.8, following the same layout as Fig. 8. Additionally shown is the active sub-sample of the dust-attenuated population (orange contours) and its observability compared to the entire active population (orange percentages). For comparison, black lines show fractions of the Eddington limit: 1 Ledd (solid), 0.1 Ledd (dashed), and 0.01 Ledd (dash-dotted).

According to our framework, both the EWS-Hα and the EDS-BPT survey configurations will recover galaxies containing AGNs with black hole masses between 10⁶ and 10^9.5 M_⊙. The distributions of the flux-limited populations occupy almost the same regions as the intrinsic sample. Adding dust-attenuation biases the distributions slightly by excluding galaxies with L_AGN below 10³⁹ erg s⁻¹. Additionally, a notable gap appears below the resolution limit for galaxies with M_BH ∼ 10⁶ M_⊙ and L_AGN ∼ 10⁴² erg s⁻¹. This gap is caused by the bimodal prescription for the AGN luminosity in GAEA distinguishing between accretion at low and high Eddington fractions (see Fontanot et al. 2020). As before, the distribution of the dust-attenuated EDS-BPT sample appears more restricted than the EWS-Hα sample.

Regarding predicted observabilities, we note that the percentages for the flux-attenuated and dust-attenuated populations are larger compared to the numbers quoted in Fig. 8. By only including galaxies with a luminous black hole more massive than GAEA’s seed masses, we biased our initial intrinsic sample towards brighter and more massive galaxies. As in Fig. 8, the EWS-Hα configuration is expected to produce higher observable percentages at 9.8% (30.7%) when considering the dust-attenuated (flux-limited) sample compared to 2.7% (21.8%) in the EDS-BPT populations.

Active galaxies in the flux-limited populations produce high AGN luminosities with L_AGN between 10⁴³ and 10⁴⁶  erg s⁻¹, entirely above the resolution limit, and fall between the expected relations for accretion at Eddington fractions of 0.01–1 (see Fontanot et al. 2020 for the implementation of black hole growth and feedback in GAEA). Between the two survey configurations, the distributions appear almost identical and produce similar observable percentages. At 1.5–3.4%, we predict the observability of the dust-attenuated active sub-sample to be low in both survey configurations. This is because the Hα and Hβ lines are particularly faint for active galaxies in this regime (also see Sect. 4.4) and thus only the galaxies with the brightest AGN luminosities produce detectable emission-line intensities.

Nevertheless, we predict line-emitting galaxies in the EWS-Hα and EDS-BPT survey configurations to contain AGNs with a wide range of masses and luminosities. In order to understand the impact on their host galaxies and disentangle observational signatures, it is thus important to be able to distinguish between SF-dominated and AGN-dominated galaxies.

4.4. Predicted observability of star-forming and active line-emitting galaxies in the EDS

In Fig. 10, we further assess the redshift distribution of line-emitting galaxies observable in the EDS, divided into SF-dominated (first and third row) and active (second and fourth row) populations. We show the GAEA-LC prediction for observable percentages of the flux-limited (outlined histograms) and dust-attenuated (filled histograms) line-emitting populations in different redshift bins. The quoted percentages are with respect to our intrinsic GAEA-LC population, which includes a mass resolution cut of 10⁹ M_⊙ and H-band magnitude cut of 25. Compared to previous figures, we chose a smaller bin spacing of Δz = 0.2 in order to sample the detectable redshift ranges of the emission lines more finely. To account for uneven redshift ranges, we allow the highest redshift bin to be larger than Δz = 0.2. We display the redshift distributions as histograms for individual emisssion lines (row 1 and 2); Hα (yellow), Hβ (light blue), [N II] (red), [O II] (indigo), [O III] (cyan), [S II] (crimson), and [O I] (green-yellow), as well as combinations of them (row 3 and 4, colours as indicated in the legend) frequently used in spectroscopic diagnostics. We also include histograms for the combination of the Hα and Hβ lines (mint green), whose ratio, the Balmer decrement, can be used to estimate dust attenuation.

Fig. 10. EDS observable fractions of SF (row 1 and 3) and active (row 2 and 4) galaxies in different redshift bins. Shown are the flux-limited (outlined histograms) and dust-attenuated (filled histograms) populations emitting Hα (yellow), Hβ (light blue), [N II]λ6584 (red), [O II]λ3727 (purple), [O III]λ5007 (cyan), [S II]λ6724 (crimson), [O I]λ6300 (green yellow), and their combinations (colours indicated in row 3 panels).

Both SF and active line-emitting galaxies exhibit decreased observability percentages with increasing redshift due to the distance-dependence of observed fluxes. We note that, overall, observability declines more steeply within active line-emitting populations. Flux-limited histograms represent an upper limit of observable percentages of line-emitters, which will be drastically reduced by the presence of dust. Estimating the precise impact of this effect is challenging (see Sect. 8.1), thus exact numbers for observable percentages of dust-attenuated line-emitters should be treated with caution. Overall, we estimate that, based on the Calzetti et al. (2000) law with mass-dependent scaling, dust attenuation will decrease the observable percentages by a further 20–30% with respect to the intrinsic GAEA-LC populations. This could particularly be a problem when trying to recover fainter lines, for which predicted percentages are reduced to single digits for lower redshifts and below 1% at higher redshifts. Additionally, wavelengths towards the blue end of the spectrum are more strongly attenuated than those towards the red end. In the following descriptions we refer to the dust-attenuated sample and include the flux-limited case as an upper limit in parentheses.

We first consider the observability of individual emission lines in SF and active galaxy populations. As noted before, Hα is by far the strongest emission line out of the ones presented here. At redshift 0.4–0.6, we expect Hα-emitters to be around 70% (reduced from 80%) observable for both SF and active galaxies in the dust scenario. Until redshift 1.6–1.8, this number decreases to around 15% (30%) for SF and 7% (20%) for active galaxies, before the Hα wavelength is redshifted out of the sensitivity range of the EDS.

[N II] is the second strongest emission line in these populations, reaching observable percentages of up to 30% (50%) in the SF galaxies and 60% (70%) in the active sample. This difference is due to the relatively high ionisation potential of nitrogen, which is more easily ionised by the harder ionising radiation from AGNs compared to young stellar populations. Furthermore, we expect the increased observability at low redshift to be partially caused by the global decrease of SFR, as well as the increase of metallicity, which strengthens [N II] emission due to secondary nitrogen production.

In both the SF and active populations, [S II]-emitters show similar observability percentages to [N II]-emitters. They decline from 35% (50%) to below (10%) 1% between redshift 0.4–0.6 and 1.6–1.8. This is due to the dependence of [S II] on SFR and interstellar metallicity, in a similar way to [N II], which is expected to increase emissions in low-redshift galaxies. However, due to the absence of any secondary production channels for sulfur, this dependence is weaker.

In contrast to line-emitting galaxies discussed above, the observability of [O III]-emitters evolves to a lesser degree with redshift. For SF galaxies, [O III] observability is between roughly 5–20% for dust-attenuated populations, whereas only flux-limited percentages stay between 30 and 40%. This can be explained by the increased [O III] emission for metal-poor galaxies, which are more prevalent at high redshift and thus compensate for the flux decrease due to increasing luminosity distance. The observable percentages of active [O III]-emitters also appear constant above redshift 1.6, which is likely caused by a larger ionisation parameter due to elevated central gas densities (see Hirschmann et al. 2017).

Hβ, [O II], and especially [O I] emission are notably fainter for the GAEA-LC populations. The predicted dust attenuation makes most of the resulting histogram bins barely visible for these lines. Around redshift one, Hβ reaches about 20–30% observability for the flux-limited sample of SF and active galaxies, while the other lines are even fainter. If our estimate is reasonable, dust would reduce the percentage of EDS-observable Hβ, [O II], and [O I]-emitters to a few per cent at the lower redshift end and below 1% at the higher redshift end.

We now turn our attention to the predicted observability of galaxies emitting multiple lines. The redshift bins have been adjusted according to the overlap of the relevant redshifted wavelengths with the EDS sensitivity range. As Hα, [N II], and [S II] are the strongest lines, the percentages for their simultaneous observability are also the highest. All combinations reach around 30% (50%) observability at redshift 0.4–0.6 in both the SF and active populations, except for [N II] with Hα, which reach close to 60% (70%) observability for active galaxies.

All other emission-line combinations contain at least one line that appears particularly weak in this redshift range. As a result, predicted observable percentages for the dust-attenuated sample drop below 1% at the higher end of their redshift ranges. At the lower redshift end, we predict generally less than 5% observability. Even though both SF and active [O III]-emitters themselves are at least 20–30% observable around redshift one, in combination with Hβ or [O II], observability stays below 5%. This represents a reduction from around 20% observability for flux-limited populations. [O II] with Hβ suffers similar losses due to dust attenuation: from 15% to less than 5% for SF galaxies at redshift 1.5 and 6% to less than 1% for active galaxies. [N II] with [O II] is also reduced from around 10% to less than 1%. Additionally, the wavelengths of [N II] and [O II] are spaced far apart, which significantly restricts the redshift range in which they can both be detected with Euclid’s grisms. This also affects the usefulness of emission-line calibrations based on these lines.

In the upper limit of the flux-limited sample, the combination of Hα and Hβ shows percentages between 10 and 25% for SF galaxies and 5–20% for active galaxies across redshift 0.9–1.8. The estimated dust-attenuation reduces these numbers drastically to at most 5% for SF galaxies and 2% for active galaxies. This is particularly due to the relatively stronger effect on the Hβ line, caused by its shorter wavelength. This results in an increased Balmer decrement (Hα/Hβ), which can be used to estimate the presence of dust in observed galaxies, if the Hβ line can still be recovered.

Galaxy populations emitting lines which make up the standard BPT diagrams, [O III], Hβ, and Hα in combination with [N II], [S II], or [O I], will roughly be observed between redshift 0.9 and 1.8. Upper limits without dust show around 5–15% observability across this range for the [N II] and [S II]–BPT combinations. In the case with dust, we predict that for both SF and active galaxies, the [N II]–BPT combination will only be roughly 3% observable at redshift 0.9–1.1 and decline to around 0.3% observability at redshift 1.5–1.8. These numbers are in agreement with the predicted observability of 1.4% for the EDS-BPT configuration in Fig. 8, which was averaged across all galaxies in the 0.9–1.8 redshift range. The [S II]–BPT line combination shows percentages of around 2% at most for both SF galaxies and the active population. As the [O I] line is already so weak by itself, the [O I]–BPT combination is never more than 1.5% observable in even the upper limit of the flux-limited populations.

In conclusion, we expect the presence of dust to reduce the number of galaxies with recoverable emission lines drastically, typically by an additional 20% with respect to the intrinsic population. This is particularly significant for intrinsically fainter lines, which will be difficult to measure already at the low end of their detectable redshift range. For Hβ, [O II], and [O I] emission we do not expect to reach more than 5% observability of SF and active galaxies. This is the case both for these lines individually and in combination with other lines. Already in the upper limit of the no dust case, we predict percentages of 10–25% at most. For stronger lines, such as Hα, [N II], [S II], and [O III], percentages are relatively high for flux-limited populations, and, as a result, dust losses are less fatal. Without dust, we recover percentages between 40–80% at lower redshifts and around 10–30% at the higher redshifts. After accounting for dust, we predict the EDS to recover up to 30–70% of both SF and active galaxies at lower redshifts, while at higher redshifts these numbers decline to below 10%. This is similarly true for all combinations of Hα, [N II], and [S II].

While some of these numbers appear discouragingly low, we stress that with a planned sky coverage of 50 deg², the EDS will nevertheless recover thousands to hundreds of thousands of SF and active galaxies emitting these emission lines and their various combinations. Detailed number count predictions for different line-emitters will be published in Cassata et al. (in prep.), who are comparing both the GAEA-LC framework and the Millennium Mambo catalogues. Furthermore, we reiterate that the exact prevalence and nature of dust in this redshift regime is fundamentally uncertain and warrants an extensive study in itself. Hence, these estimates should be considered carefully.

5. Distinguishing between dominant ionising sources in galaxy populations observable in the Euclid Deep Survey

Measurable line intensities are particularly valuable for their potential to characterise observed sources using spectroscopic diagnostics. In these remaining sections, we explain how we verified which diagnostics can be applied to future Euclid data.

As explained in Sect. 1, it is unclear whether standard BPT diagrams are a robust diagnostic to separate SF-dominated, composite, and AGN-dominated galaxy populations beyond redshift one (Hirschmann et al. 2019, 2023a; Kocevski et al. 2023). A successful application to Euclid-observable populations, in particular, also depends on the properties of these galaxies, as theoretical works (e.g. Groves et al. 2006; Feltre et al. 2023; Hirschmann et al. 2019) indicate that the interstellar metallicity of galaxies, in addition to the type of ionising radiation, influences their location on the standard [O III]/Hβ versus [N II]/Hα BPT diagram.

In Fig. 11, we thus verify if the EDS-BPT sample of GAEA-LC galaxies at redshift 0.9–1.8 conforms to locally calibrated optical selection criteria for the [O III]/Hβ-versus-[N II]/Hα BPT diagrams. As in Sect. 3.2, we divide the GAEA-LC sample into different types according to the BHAR/SFR criterion, resulting in SF-dominated (left column), composite (middle column), and AGN-dominated galaxies (right column). In order to assess a potential evolution of the location of different galaxy types in the diagram, we split the sample into the redshift bins 0.9–1.2 (top row), 1.2–1.5 (middle row), and 1.5–1.8 (bottom row). The general layout and colour-code is the same as in the top panels in Fig. 3, except that now showing the flux-limited (grey contour lines and data points) and dust-attenuated populations (coloured contours and data points) for each galaxy type. As in Fig. 3, we show individual data points for the AGN-dominated sample to avoid artefacts caused by the discrete model grids and the small number of AGN-dominated galaxies surpassing the flux limit.

Fig. 11. Location of EDS-observable GAEA-LC galaxy populations in the [O III]/Hβ versus [N II]/Hα BPT diagram, divided into redshift bins 0.9–1.2 (top row), 1.2–1.5 (middle row), and 1.5–1.8 (bottom row). Shown are the flux-limited (grey contour lines and data points) and dust-attenuated samples (coloured contours and data points). The columns and colour coding for different dominant ionising sources, the empirical selection criteria to distinguish between them and the SDSS-observed galaxies, shown for comparison, follow the same layout as Fig. 3.

For both the flux-limited and dust-attenuated cases, the predicted galaxy populations conform well to locally calibrated criteria (Kewley et al. 2001; Kauffmann et al. 2003). We note that, compared to the flux-limited populations, the dust-attenuated populations exhibit reduced ranges of [O III]/Hβ and [N III]/Hα. This trend gets more pronounced with increasing redshift. While there is some overlap of the SF-dominated galaxies into the transitional composite region, we note that the majority of galaxies is confined below the criterion by Kauffmann et al. (2003). For the SF-dominated populations, we also see a slight increase with redshift of [O III]/Hβ at fixed [N II]/Hα, consistent with Fig. 3. Nevertheless, they remain well-separated from AGN-dominated galaxies at all redshifts shown here.

In contrast, Hirschmann et al. (2019), who use the same photoionisation models, found that a global increase of [O III]/Hβ and decrease of [N II]/Hα with redshift makes SF and active galaxies indistinguishable already at redshift one. They attribute this to a drop in interstellar metallicity (in addition to a rise in SFR, also see Sect. 4.2), but demonstrate that metal-rich galaxy types remain separable above redshift one. From Fig. 8 we know that the EDS-BPT selected GAEA-LC population is biased towards high oxygen abundances, which explains why, for this sample, the galaxy types remain distinguishable.

We conclude that the [N II]–BPT diagram is a reliable diagnostic diagram to distinguish between dominant ionising sources in the EDS-observable galaxy populations at redshift 0.9–1.8. However, as seen in Sect. 4, the relevant sample of galaxies represents only a few per cent of the predicted intrinsic population. Additionally, we remind the reader that de-blending [N II] and Hα in Euclid spectra will be challenging (see Sect. 8.4), which will further reduce the sample for which diagnostics relying on separable [N II] and Hα estimates can be used.

6. Deriving ionising properties from strong line luminosities at intermediate redshifts

Optical emission-line intensities, and their ratios, are also widely used as proxies to estimate key galaxy properties, as they encode signatures from both the local ISM gas, as well as the nature of the ionising radiation. Using observational data and photoionisation models, various emission-line calibrations have been derived for ionising properties, such as the SFR and AGN luminosity. While they have been validated for the interpretation of observational spectra at low redshift, it is unclear if their use can be extended beyond the local Universe. In the more distant Universe, the ionisation and ISM conditions might differ significantly which could cause the locally assumed relationship between emission lines and galaxy properties to break down.

In this section, we thus aim to examine the potential evolution of locally established calibrations across the intermediate redshift range 0.4–2.5, where Euclid will recover the majority of the relevant emission lines. This is meant as a guide for future data releases to determine which, if any, of the widely used calibrations could be used to characterise sources. Additionally, we present new calibrations which relate [N II]-based emission-line ratios to the BHAR/SFR ratio. As here we focus on the intrinsic physical limitations of these relations predicted by our framework, we show results from the intrinsic GAEA-LC population, as well as indicate potential biasing as a result of Euclid-like flux cuts.

6.1. Deriving SFR and AGN Luminosity

The SFR and AGN luminosity play a crucial role in shaping the evolution of a galaxy across cosmic time. They are also the main drivers of ionising radiation, thus setting the degree of ionisation of the ISM gas and the strength of different emission lines. In Fig. 12, we show average Hα and [O III]/Hβ emission-line luminosities versus the SFR for SF galaxies (top and middle panels) and the [O III] versus AGN luminosity for active galaxies (bottom panel). We follow the layout of Fig. 5, but now extend our predictions to higher redshift intervals between redshift 0.4 and 2.5 (solid lines, colours indicated in legend), with a scatter of one standard deviation for the 0.4–0.9 bin (grey shaded area). For completeness, we plot the emission-line relations in all redshift bins (0.4–0.9, 0.9–1.2, 1.2–1.5, 1.5–1.8, and 1.8–2.5) and indicate where the emission line is within Euclid’s sensitivity range with thicker lines. For these lines we further show the impact of the EDS flux limit (dashed lines in fainter colours) and dust attenuation (dotted lines in fainter colours).

Fig. 12. Average Hα (top panel) and [O III] + Hβ (middle panel) line luminosities versus the SFR for GAEA-LC star-forming galaxies in different redshift intervals (shown in different coloured solid lines as indicated by the legend). For completeness, the three relations are shown for the same redshift intervals, but thick lines indicate that the emission line is within the EDS wavelength sensitivity range. For those, the flux-limited (dashed lines) and dust-attenuated relations (dotted lines) are plotted in fainter colours. As in Fig. 5, predictions are compared to local calibrations (Kennicutt & Evans 2012; Osterbrock & Ferland 2006). In the bottom panel, the average [O III] line luminosity is plotted against AGN luminosity for active GAEA-LC galaxies and the relation found by Lamastra et al. (2009, dash-dotted). The coloured lines for different redshift bins follow the same key as above.

We note only a negligible, if any, evolution for all relations, meaning they continue to broadly follow locally used relations up to at least redshift 1.8–2.5. In fact, the Hα relation, which was offset at redshift 0–0.3, approaches the Kennicutt & Evans (2012) relation more closely at higher redshifts. This evolution is caused by a systematically higher ionisation parameter at higher redshifts, which implies a larger H II region and an associated greater dust attenuation. The predicted L_{[O III]+Hβ}–SFR relation exhibits no clear evolution with redshift and continues to agree well with the Osterbrock & Ferland (2006) relation. In both relations, the GAEA-LC magnitude cut causes a luminosity increase at low SFRs for high redshift, as here the light cone only includes a few unusually bright objects. For the L_[O III] − L_AGN relation, we predict an increase up to 0.5 dex above the Lamastra et al. (2009) relation. However, the general slope stays the same and we find that, especially in the lower redshift bins, Lamastra et al. (2009) provide a good estimate of the AGN luminosity.

When applying the EDS flux limit and Calzetti et al. (2000) dust attenuation, predicted relations show a systematic bias towards increasingly bright line emitters with increasing redshift. At redshift 0.4–0.9, L_Hα is a robust proxy for the SFR up to around 10⁴¹ erg s⁻¹, while at redshift 1.5–1.8, we can expect L_Hα below 10⁴²–10^42.5 erg s⁻¹ to result in biased SFR measurements. For L_{[O III]+Hβ}, we predict that at redshift 0.9–1.2 the relation can be used down to roughly 10⁴² erg s⁻¹, while at redshift 1.8–2.5 around 10^42.5–10⁴³ erg s⁻¹ are needed to ensure results are unbiased. The L_[O III-L_AGN relations from flux-limited and dust-attenuated populations follow a similar pattern. However, they appear truncated compared to the intrinsic population, as no active galaxies with L_AGN below 10⁴² erg s⁻¹ and L_[O III above 10^41.5 erg s⁻¹ exist in our sample. If such objects are recovered in the EDS, there are likely unusual and extreme objects.

One other disadvantage of these line intensity calibrations is that they are only valid for the dominant ionising property in SF and active galaxies, respectively. Otherwise the relation is contaminated by the contribution from other sources (see also Fig. 10 in Hirschmann et al. 2023a). This requires a prior classification according to galaxy type. Using standard BPT diagrams for this purpose is only possible if [O III], Hβ, Hα, and [N II] can be simultaneously measured. However, as seen in Sect. 4, we expect this to only be the case for around 1.4% of galaxies in our underlying sample, with additional reductions due to line blending. Thus, it is useful to explore new emission-line based tracers which could constrain the strength of ionising sources across all galaxy types.

6.2. [N II]-based emission-line ratios as new tracers for BHAR/SFR

The BHAR/SFR ratio measures the relative strength of the star formation and black hole activity, which defines the balance of energetic processes driving both the ionisation of the ISM, as well as the evolution of the entire galaxy. In this work, we have utilised the BHAR/SFR ratio as a theoretical criterion to distinguish between dominant ionising sources in galaxies. Both Figs. 3 and 11 have demonstrated that this distinction is successful at visually separating the different galaxy types in standard BPT diagrams. The [N II] BPT diagrams show a particular dependence of the three populations on the [N II]/Hα ratio, which we aim further explore in this section.

In Fig. 13, we show the location of SDSS-like GAEA-LC populations between redshift 0.4 and 2.5 in the [O III]/Hβ versus [N II]/Hα BPT diagram. We follow the same layout as in the top row in Fig. 3, but plot individual data points for the SF-dominated, composite and AGN-dominated galaxies, coloured according to their BHAR/SFR ratio. For reference, we indicate the BHAR/SFR boundaries between the three types (black vertical lines in colour bar). As we aim to explore the physical relationship between these quantities, we employ no dust correction or Euclid-like flux cuts.

Fig. 13. Location of the SDSS-like GAEA-LC galaxy populations in the [O III]/Hβ versus [N II]/Hα BPT diagram with the same layout as the top row in Fig. 3, but now with individual data points colour-coded by their BHAR/SFR ratio. Shown are galaxies above the SDSS flux cut in the redshift range 0.4–2.5.

As seen before, the three galaxy populations are well-separated by the BHAR/SFR criteria, with some overlap at the edges of the SF-dominated and AGN-dominated populations with the area occupied by composite galaxies. We further point out that even across galaxies of the same type there is a clear BHAR/SFR-dependence of the [N II]/Hα ratio. The [N II]/Hα is mostly driven by changes in [N II], while the Hα line provides a roughly constant baseline. Thus we conclude that it is the [N II] emission which depends on the relative strength of star-forming and accretion processes. This can be explained with the high excitation energy of this N⁺ state, which traces the hardness of the ionising radiation from AGNs. As a result, we further explore [N II]-based emission-line ratios as potential tracers for the BHAR/SFR ratio.

Thus far, determining the BHAR/SFR ratio of a galaxy relied on combining separate estimates for the SFR and the BHAR. Results are heavily dependent on AGN selection and derivation methods, which results in large combined uncertainties, especially with increasing redshift (see McDonald et al. 2021). In Fig. 14, we propose three novel [N II]-based emission-line calibrations to the BHAR/SFR for intermediate redshifts. We show the average [N II]/Hα (N2, left panel), [N II][O II] (N2O2, middle panel), and [N II][S II] (N2S2, right panel) at fixed BHAR/SFR in the same redshift intervals as in Fig. 12 and an additional redshift interval between 0 and 0.4 (indigo). As before, we vary line thickness according to the detectability in the EDS and show flux-limited populations in dashed lines. We omit dust-attenuated populations, as, since we are considering line ratios, these do not present any significant differences.

Fig. 14. Average ratios of [N II]/Hα (left panel), [N II][O II] (middle panel), and [N II][S II] (right panel) versus BHAR/SFR in different redshift intervals (same key as in Fig. 12, with an additional redshift 0–0.4 line in indigo). Overplotted are the best linear fits between −4 ≤ log10(BHAR/SFR)≤ − 1 (dotted lines). Respectively shown for N2, N2O2, and N2S2 are the 2, 1.5, and 1σ scatter of the entire 0–2.5 redshift range, which each encompass the 12 + log10(O/H) = 8 − 9 range (grey shaded area).

For all three emission-line ratios, we find a clear positive correlation with the BHAR/SFR ratio, for which we derived linear best fits (dotted lines) in the log₁₀(BHAR/SFR) range between −4 and −1. We note that N2, N2O2, and N2S2 are also frequently used tracers for the O/H abundance (further explored in Sect. 7), as increased metallicity boosts the production of secondary nitrogen, which is modelled via the Groves et al. (2004) analytic relation between N/O and O/H in both the SF and AGN models (Gutkin et al. 2016; Feltre et al. 2016). In order to ensure that the relations in Fig. 14 are not driven by relatively higher O/H abundances at higher BHAR/SFR, we verified that the dependencies persist at fixed metallicities. We found that increased O/H abundances move the relation towards higher emission-line ratios, but exhibit similar slopes. This effect is stronger at metallicities below 12 + log₁₀(O/H) = 8. Between 12 + log₁₀(O/H) = 7.5 and 8, N2 increases by around 0.6 dex at a given BHAR/SFR, while N2O2 and N2S2 increase by up to 0.7 dex and 0.5 dex, respectively. However, above 12 + log₁₀(O/H) = 8, the location of the relations generally change more gradually with O/H abundance. This is the metallicity range which we expect most Euclid-observable emission-line galaxies to fall in. For N2, N2O2, and N2S2 relations between 12 + log₁₀(O/H) = 8 and 9 lie within the 2, 1.5, and 1σ around the average relation (grey shaded area). This corresponds to a range of around 1 dex for N2 and N2O2 and 0.3 dex for N2S2. These scatters are smaller than the overall BHAR/SFR dependence we find in the relation determined using all O/H abundances.

For the N2 ratio, our derived relation increases by 1.25 dex between log₁₀(BHAR/SFR) = −4 and −1:

$$\begin{array}{c}\hfill \mathrm{N}2=0.42{\log }_{10}(\mathrm{BHAR}/\mathrm{SFR})+0.77.\end{array}$$(1)

N2O2 covers a range as large as 1.6 dex, with a slightly steeper best fit:

$$\begin{array}{c}\hfill \mathrm{N}2\mathrm{O}2=0.52{\log }_{10}(\mathrm{BHAR}/\mathrm{SFR})+1.36.\end{array}$$(2)

This strong correlation with N2O2 arises due to additionally decreased [O II] emission for high BHAR/SFR due to a greater abundance of multiply ionised oxygen at the expense of singly ionised oxygen, caused by harder ionising radiation from AGNs. Both the N2- and N2O2–BHAR/SFR relation appear to be almost redshift-invariant. The N2S2 ratio covers a range of 0.9 dex with a less steep correlation of

$$\begin{array}{c}\hfill \mathrm{N}2\mathrm{S}2=0.31{\log }_{10}(\mathrm{BHAR}/\mathrm{SFR})+1.32.\end{array}$$(3)

N2S2 also exhibits a slight rise of around 0.2–0.3 dex with increased redshift. However, N2S2 is overall better constrained than N2 and N2O2, as both the increase with redshift and the metallicity variations between 12 + log₁₀(O/H) = 8 and 9 range are encompassed within a 1σ scatter. Additionally, average relations derived from flux-limited populations appear unbiased across the relevant BHAR/SFR range for N2S2, while N2 and N2O2 start to differ from the intrinsic relations below log₁₀(BHAR/SFR) of −3. This biasing becomes significant at the low end of the linear fit, with an increased line ratio of at most 0.3 dex with respect to the intrinsic population. In both cases, this biasing lies within the 12 + log₁₀(O/H) = 8 − 9-scatter. While the N2O2–BHAR/SFR correlation is the steepest, we also predict the [O II] line to be relatively faint over Euclid’s sensitivity range. This makes N2S2 the most favourable tracer to constrain the BHAR/SFR ratio in Euclid data. However, we consider all three of these line ratios to be effective tracers, which could allow for the characterisation of galaxies using just two emission line measurements.

7. Deriving interstellar metallicity

The interstellar metallicity provides valuable insight into the chemical processes which have enriched the galaxy and, as a result, imprinted into the spectrum by influencing the relative intensities of various emission lines (see Kewley et al. 2019; Maiolino & Mannucci 2019, and references therein). The direct temperature method uses collisionally excited auroral lines, such as [O III]λ4363 and [N II]λ5755, to estimate the electronic temperature T_e and the underlying chemical abundances. Unfortunately, auroral lines are fairly faint and thus primarily detectable for nearby galaxies. However, they have been used to calibrate relationships between the metallicity and ratios of various strong emission lines, which are more easily detectable (e.g. Bresolin 2006b; Pérez-Montero & Contini 2009; Jones et al. 2015; Curti et al. 2019). This approach is particularly successful for metal-poor galaxies, while in metal-rich galaxies the lack of collisional excitation due to temperatures lowered by metal-line cooling makes the auroral lines even fainter. Additionally, at high metallicities the direct T_e saturates due to temperature fluctuations and gradients, as well as departures from thermal equilibrium (Kewley et al. 2019; Maiolino & Mannucci 2019). A potential result is that metallicities determined via the direct T_e method are generally around 0.2 dex lower compared to those derived from metal recombination lines (Bresolin 2006a; Stasińska 2005). At very low metallicities, the calibration of strong-line estimators is instead limited by sparse data, since emission lines are generally more faint, thus complicating the simultaneous measurement of multiple lines (Osterbrock & Ferland 2006).

Some studies resort to photoionisation models to establish relationships between metallicity and strong-line ratios, either in combination with direct T_e estimates for the metal-poor end (Pettini & Pagel 2004; Nagao et al. 2006; Maiolino et al. 2008; Marino et al. 2013) or to cover the entire range of metallicities (Kewley & Dopita 2002; Tremonti et al. 2004). While these are able to cover a wide range of metallicities and line ratios, their parameter space is often poorly constrained (see Chevallard & Charlot 2016; Vidal-García et al. 2024) and thus limited in their predictive power.

As with the emission-line relations for SFR and AGN luminosity in Sect. 6.1, most line-ratio calibrations have been derived and validated at redshift zero and it is unclear if they provide robust estimates at higher redshifts. Given the increase of equivalent width and specific star formation rate with redshift, the most suitable local calibrations for metallicity at high-redshift should be those constructed with local analogues matching these properties (Bian et al. 2018; Pérez-Montero et al. 2021; Nakajima et al. 2023). Bian et al. (2018) found that direct T_e estimates for local redshift two analogues disagree with the usual calibrations at redshift 0, hinting that these differing properties lead to altered emission-line relationships. New JWST measurements at redshift 2–9 show similar discrepancies to the most widely used local calibrations (e.g. Curti et al. 2022; Laseter et al. 2024; Sanders et al. 2024).

In Fig. 15, we thus explore if commonly used emission-line ratios (as detailed in caption) for the interstellar metallicity (usually expressed as the oxygen abundance O/H) evolve across the Euclid-detectable redshift range. As, thus far, most studies have focussed on SF galaxies, we only show relations for SF-dominated GAEA-LC galaxies (for specific calibrations to AGN narrow-line regions, see Dors et al. 2017, 2021; Carvalho et al. 2020). We use the same colour code as in Fig. 14. Thicker lines indicate that all emission lines making up the ratio fall into the detectable range of the EDS. For visual clarity, we do not include average relations for flux-limited populations in this case, but comment on biasing in the final paragraph of the Section.

Fig. 15. Average ratios of strong emission lines against interstellar O/H abundance. We consider in the left column: N2, [S II]/Hα (hereafter S2), N2S2, and N2O2; in the middle column: [O II]/Hα (R2), [O III]/Hα (R3), [(O III] + [O II])/Hα (R23), and ([O III/Hα) + ([S II/Hα) (RS23); and in the right column: [O III]/[O II] (O32), ([O III]/Hα)/([N II]/Hα) (O3N2), and ([O III]/Hα)/([S II]/Hα) (O3S2). Solid lines show each relation for the GAEA-LC SF-dominated populations in different redshift bins, using the same colours as in Fig. 14, with thick lines indicating overlap of the relevant emission lines with the EDS sensitivity range. Overplotted are empirical and theoretical calibrations at redshift zero (blue and purple dashed lines, respectively; grey indicates calibrations using a combination of both; Curti et al. 2019; Jones et al. 2015; Pérez-Montero & Contini 2009; Bresolin 2006b; Tremonti et al. 2004; Kewley & Dopita 2002; Marino et al. 2013; Maiolino et al. 2008; Nagao et al. 2006; Pettini & Pagel 2004). Additionally shown are recent estimates for higher redshifts (green dashed lines and data points Bian et al. 2018; Sanders et al. 2024; Laseter et al. 2024).

Alongside our predictions, we plot empirical and theoretical calibrations for redshift zero (direct T_e method: blue dashed lines, photoionisation models: purple dashed lines, combination of both: grey lines). Additionally shown are metallicity estimates for redshift two analogues (Bian et al. 2018, light green dashed lines) and recent JWST measurements of galaxies at redshift 2–4 from the JADES Laseter et al. (2024, dark green data points) and CEERS Sanders et al. (2024, medium green data points) surveys. Where available, we further include the metallicity calibrations for the entire sample of CEERS galaxies at redshift 2–9 from Sanders et al. (2024, turquoise dashed lines).

Overall, we find that the average relations between line ratios and the oxygen abundance broadly agree with the local calibrations, especially at low redshifts. While N2O2 and R23 appear fairly constant across all redshift bins, other relations show some evolution with redshift. These results are in excellent agreement with Hirschmann et al. (2023b), who find the same broad evolution in their emission-line coupled IllustrisTNG galaxies between redshift two and zero. We note that this evolution presents itself as a change in normalisation; metallicity estimators are progressively shifted towards lower or higher values, while the shape generally remains the same. For each estimator, we thus determined the shifts at 12 + log₁₀(O/H) = 8 between our predicted average relation at redshift 0–0.4 and the relations at redshift one, 1.5, 2, and 2.5. The results are shown in Table 1. Depending on metallicity, the magnitude of the shifts can vary, and as such, we provide approximate ranges within which the specified shifts are applicable. In the following paragraphs, we assess each estimator in more detail and explain the evolutionary shifts between redshift 0–0.4 and 2.5.

Looking at the N2, S2, N2S2, and N2O2 estimators, we note that they all show a positive correlation with oxygen abundance as metal-rich galaxies contain elevated nitrogen and sulfur abundances (see Pérez-Montero et al. 2013; Pérez-Montero 2014; Diaz & Zamora 2022). S2, and to a much lesser degree N2, shows a drop with increasing redshift. This can be attributed to the elevated ionisation parameter at fixed metallicity and stellar mass (see Figs. A.1 and A.2), which causes nitrogen and sulfur to favour the higher ionisation states N²⁺ and S²⁺. For N2, the resulting drop of −0.46 dex lies within the scatter of local literature calibrations. S2 shows a greater evolution of −0.87 dex, which translates into a +0.41 dex rise with redshift in the N2S2 estimator. For N2S2, the size of the shifts varies more significantly with metallicity, as for high metallicities the nitrogen abundance is increased with respect to sulfur due to secondary production channels. As a result, the given shift is only valid for metallicities around 12 + log₁₀(O/H) = 8. For the N2O2 ratio, the changing ionisation parameter has little influence, as nitrogen and oxygen have similar ionisation energies, which explains its near redshift-invariance. As the standard deviation associated with each average shift is larger than the shift itself at all redshifts, we do not include a metallicity range for N2O2.

The R2, R3, R23, and RS23 all show a peak in their average relations around 12 + log₁₀ (O / H) ≈ 8. They are all oxygen-based and thus initially increase with oxygen abundance. However, at high metallicities, low temperatures caused by metal-line cooling make collisional excitation less likely. At high redshift, the increased ionisation parameter due to higher SFRs causes oxygen to favour O²⁺ over O⁺, which is why R2 shifts towards lower and R3 and RS23 towards higher values. This amounts to changes of −0.52, +0.37, and +0.35 dex. The R23 ratio is able to compensate some of this evolution by summing the [O II] and [O III] line contributions, which results in a comparatively small shift of +0.16 dex in the metallicity estimator.

The combined effects of stronger [O II], [N II], and [S II], as well as decreased [O III] emission with rising metallicity cause O32, O3N2, and O3S2 to drop steeply. Increased abundances of doubly ionised atoms at higher redshifts, at the detriment of O⁺, N⁺, and S⁺, cause up to 1.24 dex higher line-ratios at higher redshifts.

In general, the CEERS metallicity determinations by Sanders et al. (2024), and specifically the galaxies at redshift 2–4, consistently agree better with our intermediate-redshift prediction from the GAEA-LC galaxies than the local calibrations. Upper and lower limits for the [N II]- and [S II]-based ratios deviate from our predictions in the expected direction. Despite containing galaxies up to redshift nine, the determined metallicity calibration for R2, R3 and O32 do not diverge significantly from our highest redshift predictions between 1.5 and 2.5. The calibrations to redshift two analogues from Bian et al. (2018) also lie closely to these lines. This close agreement is encouraging and further highlights that metallicity estimators already evolve considerably at intermediate redshifts. The R3 and R23 line ratios for four JADES galaxies at redshift 2–4 determined by Laseter et al. (2024) show a similar trend to the Sanders et al. (2024) galaxies. However, the R2 and O32 estimates exhibit a scatter of roughly 1 dex across the low-metallicity range between 12 + log₁₀(O/H) = 7.27 and 7.78. This scatter is also present in their measurements for nine additional JADES and the twelve additional CEERS galaxies with redshift above four from Sanders et al. (2024, not shown here, see Figures 4–8 in Laseter et al. 2024), from which they conclude that the metallicity dependence of R2 and O32 breaks down in the high-redshift Universe. No distinction is made for the intermediate-redshift range. Given the apparent agreement of our predictions with Bian et al. (2018) and the Sanders et al. (2024) redshift 2–4 galaxies, the small number statistics of these samples and the systematics resulting from varying assumptions, data reduction and metallicity precriptions (see Sect. 3 in Laseter et al. 2024), it is however difficult to determine if and when exactly this breakdown might occur. This further illustrates the need for additional spectroscopic data at intermediate redshifts.

Overall, we conclude that across the redshift range 0.4–2.5, most commonly used emission-line ratios evolve away from their local O/H calibrations. This is a consequence of an elevated ionisation parameter at higher redshifts, largely caused by increased sSFRs. When separating the galaxies into low and high sSFR samples, we find that the low-sSFR relations are located at the low-redshift ends in Fig. 15, while at high sSFR the relations lie closer to the higher redshift ends. The N2O2 and R23 relations evolve the least, making them, in theory, the most robust to extend to higher redshifts without modifications. For all relations we observe a shift in the normalisation, while the overall shape is preserved. Not accounting for this evolution could result in biased metallicity estimates. For instance, Hirschmann et al. (2023b) found that classical calibrations at redshift zero underestimate the mass-metallicity relation by up to 1 dex when applied above redshift four. We expect a similar trend to occur already up to redshift 2.5. Where possible, calibrations such as those of Bian et al. (2018) and Sanders et al. (2024) should be utilised instead, as they provide a better match to our predicted relations at redshifts above one.

For the purpose of deriving metallicity in Euclid-observable galaxies specifically, the detectability of lines according to the grism sensitivity ranges and the observability according to flux limits should also be considered. We expect Euclid’s flux limit to bias the observable galaxy population to the brightest, most massive and metal-rich galaxies, which could deviate from predicted relations. In order to test the robustness of various metallicity estimators in a Euclid-like sample, we thus computed the average relations for EDS-observable populations. For all average relations except R3, no metal-poor galaxies passed the flux cut, resulting in a truncated 12 + log₁₀(O/H) range of 8–9. If the EDS recovers metal-poor galaxies outside this truncated range, they will most likely be extreme objects. Inside the metallicity range, R2, R3, R23, RS23, O32, O3N2, and O3S2 did not deviate significantly from the predictions in Fig. 15. N2, S2, N2S2, and N2S2 showed an 0.5–1 dex overestimation of the line ratio below 12 + log₁₀(O/H)≈8.5. Additionally, the deconvolution of blended [N II] and Hα lines (method detailed in Sect. 8.4) might not always be possible. We thus recommend the use of [OII]- and [O III]-based metallicity estimators for Euclid data. However, for line combinations that are widely spaced in wavelength, the line ratio is more likely to suffer from differential dust extinction. Therefore, when applying calibrations for N2O2, O3N2, O3S2, O32, R23, and RS23 some form of dust correction should be considered.

8. Discussion

Based on our self-consistent modelling framework, we have provided forecasts for line-emitting galaxies in the intermediate redshift range 0.4–2.5, thereby addressing the lack of theoretical guidance for spectroscopic diagnostics in this regime. In anticipation of the Euclid Wide and Deep Surveys, we indicated the expected biasing of Euclid’s future spectral catalogues as well as made recommendations for the use of a wide range of spectroscopic diagnostics to characterise observed galaxies. We consider these predictions robust, as in Sect. 3, we validated our approach by showing that the emission-line coupled GAEA-LC framework successfully reproduces a range of key observations and calibrations. However, we acknowledge that our study may be affected by several caveats related to the treatment of interstellar dust, modelling details, and Euclid’s instrumental effects. We discuss these in sections below. Caveats that are not directly related to GAEA are largely similar to those discussed in Sect. 5.2 of Hirschmann et al. (2023b), which we refer the reader to for additional details.

8.1. Treatment of interstellar dust

We expect the presence of dust to affect Euclid observations, primarily by obscuring emitted line fluxes along the line-of-sight to objects. While dust attenuation is treated self-consistently within ionised regions using Cloudy, estimating the impact of interstellar dust is more challenging, especially beyond redshift zero, due to limited data coverage. GAEA also does not incorporate a self-consistent prescription of dust processes, such as dust grain formation, evolution, and depletion within the interstellar medium. However, the limited data available suggests that the locally found stellar mass-dependent scaling of dust attenuation seems to hold for higher redshifts. Thus, we used the Calzetti et al. (2000) curve based on mass-dependent A_V from the Garn & Best (2010) relation to estimate the attenuation effects on our modelled emission-line catalogue.

We find that interstellar dust leads to substantial reduction in line fluxes, particularly for bluer wavelengths, where observable percentages are reduced by up to 30%. Consequently, intrinsically fainter lines, such as Hβ, [O II], and [O I], may become difficult to observe. This effect is more pronounced at higher redshifts, where predominantly massive galaxies produce intrinsic line fluxes above Euclid’s nominal flux limit. Due to the mass-dependence of the scaling, they are also more strongly attenuated, which would decrease many estimated line fluxes below the threshold. However, due to large uncertainties involved, these estimates should be understood as an indication of how the presence of dust might affect observed samples rather than precise predictions. We further note that the Garn & Best (2010) relation is purely empirical and is not derived from first principles.

There are a variety of complicating factors. The environmental conditions at higher redshifts are largely unknown and variations in dust grain composition, size, shape, and distribution can directly affect the resulting attenuation (for a review, see Draine 2003). Additionally, we assume galaxies of the same stellar mass also have the same dust mass, disregarding other potential dependencies, such as a positive correlation between extinction and the star formation rate (i.e. Hopkins et al. 2001; Zahid et al. 2013). Further, we expect the presence of a dusty torus around AGNs to affect line emission differently compared to interstellar dust (Urry & Padovani 1995; Hasinger 2008). This would introduce additional complexities in estimating the overall attenuation effects, especially for composites.

Lastly, we note that during the coupling process we adopted a dust-to-metal mass ratio ξ_d of 0.3 for all galaxies (see Sect. 2.2.2), which sets the depletion of metals onto dust grains. This is inspired by the Solar-neighbourhood value ξ_d, ⊙ of 0.36. Increasing ξ_d would result in greater depletion of metal coolants and, thus, higher temperatures and greater probability of collisional excitation. Significantly different or evolving values for ξ_d could thus have an effect on the resulting line emission. However, Hirschmann et al. (2017) explored the influence of setting different ξ_d during the coupling process and found only a negligible influence on the cosmic evolution of simulated emission-line ratios, in agreement with observational and other theoretical works (Rémy-Ruyer et al. 2014; Popping et al. 2017).

8.2. Caveats related to photoionisation models

The model grids for SF, AGN, and PAGB contributions used in this study are based on the Cloudy photoionisation code. It follows the non-equilibrium ionisation and the thermal and chemical state of a gas element under incident radiation, calculating ionisation, recombination, collision, emission, and absorption, as well as accounting for various other physical processes, such as radiation pressure on dust grains, metal depletion on dust, and attenuation. Computations are performed in one dimension, assuming spherical geometry and constant density. This represents a simplification, as real galaxies exhibit complex 3D gas distributions, which can not be captured by the models.

We further assume that all ionised regions are ionisation-bounded, as opposed to density-bounded. Density-bounded regions are optically thin to Lyman Continuum (LyC) photons, allowing them to leak out and affect emission-line intensities. While LyC leakage is more easily identified at low-redshift (Flury et al. 2022a,b; Izotov et al. 2021), this process appears particularly important at higher redshifts (i.e. Shapley et al. 2015, 2016; Bian et al. 2017; De Barros et al. 2019). The general effect on emission line properties is a decrease in the intensities of low-ionisation relative to high-ionisation lines (Jaskot & Oey 2013; Plat et al. 2019). However, LyC leakage remains poorly understood and it is unclear to what extent it might affect diagnostics at intermediate redshifts, such as the ones presented in this work. Given the properties of observed LyC emitters at low and intermediate redshifts, namely oxygen abundances between 12 + log₁₀(O/H) = 7.5 (Steidel et al. 2018) and 8.5 (Flury et al. 2022a; De Barros et al. 2016; Vanzella et al. 2016), stellar masses between 10⁸ and 10^10.7 M_⊙, as well as high sSFRs and thus high resulting ionisation parameters and equivalent widths, we predict that Euclid sample will observe some emission lines belonging to LyC emitters.

In contrast to recent works with similar methodology (i.e. Hirschmann et al. 2023a,b), we did not include any emission-line models for fast radiative shocks. Shocks can be produced by galactic outflows caused by supernovae, stellar winds, and AGNs (e.g. Rich et al. 2010; Sharp & Bland-Hawthorn 2010; Rich et al. 2011; Soto et al. 2012; Weistrop et al. 2012) and have been observed out to redshift three (e.g. Steidel et al. 2010; Genzel et al. 2011; Kornei et al. 2012; Newman et al. 2012). These shocks can cause excitation in the ISM gas, which contributes to the overall line emission.

In order to model line emission from shocks self-consistently, Hirschmann et al. (2023a,b) take advantage of a shock finder which has been applied to IllustrisTNG on the fly. Fast radiative shock models from the Mappings V-based grids by Alarie & Morisset (2019) were then matched to the averaged quantities of shocked regions. As our modelling framework is based on GAEA, a semi-analytic model, we do not have access to shocks and other internal gas kinematics (see also Sect. 8.3). Hirschmann et al. (2023a) predict that at low redshifts, shock-dominated galaxies produce similar signatures to AGNs on the BPT-diagram, making them hard to distinguish. However, since shocks are usually caused by stellar and AGN-driven outflows, purely shock-dominated galaxies are rare. They find that at redshifts above one, fractions of shock-dominated galaxies drop below 1% and produce a negligible contribution to the overall line emission. In light of this, we conclude that the inclusion of fast radiative shocks would not influence our results significantly.

8.3. Caveats related to the GAEA semi-analytic model and the coupling methodology

While semi-analytic models, such as GAEA, differ in their modelling methodology to cosmological simulations, they are subject to similar uncertainties. Complex physical processes and mechanisms, such as star formation and stellar and AGN feedback, have to be simplified into model prescriptions in order to represent the complex evolution of baryonic components. As for cosmological simulations, resulting galaxy properties can depend substantially on the specific model chosen. However, one advantage of semi-analytic models is that, due to their computation speed, it is possible to run the model with various implementations. This allows the exploration of large parameter spaces, ultimately enabling the choice of the most successful scheme and the subsequent fine tuning of parameters to reproduce observational benchmarks. This has been done for both stellar and AGN feedback (see Hirschmann et al. 2016; Fontanot et al. 2020, respectively). The resulting combined model we make use of in this work reproduces key observational constraints out to high redshift. The galaxy mass function is robust out to approximately redshift seven and the cosmic star formation rate density out to redshift ten (Fontanot et al. 2017). GAEA also agrees well with the observed cold gas fractions out to around redshift two, as well as the mass–metallicity relation at redshift zero and its evolution for increasing redshifts (Hirschmann et al. 2016). Fontanot et al. (2020) have further shown that the bolometric AGN luminosity function is in agreement with observations up to roughly redshift four. This forms a robust foundation for our emission-line predictions at intermediate redshifts. We note however, that a preliminary comparison to observations (De Lucia et al. 2018) indicates that GAEA tends to underpredict quiescent fractions for massive galaxies at high redshift, which could introduce a bias in our results.

Given that GAEA does not explicitly treat internal gas distribution and the related dynamics, additional assumptions and simplifications were necessary in the coupling process, which are detailed in Sect. 2.2.2. Following preceding works, we fix the hydrogen gas density within ionised regions to 10² cm⁻³ for H II regions, 10³ cm⁻³ for narrow-line regions, and 10 cm⁻³ for line-emitting regions ionised by post-AGB stars. As discussed in Sect. 8.1, the dust-to-metal mass ratio ξ_d has been set to 0.3. The impact of varying these parameters on predicted line emission is relatively small and has been discussed in Hirschmann et al. (2017, 2019). For the computation of the ionisation parameter, we initially calibrate the volume filling factor at redshift zero to reproduce the Carton et al. (2017) relation. At higher redshifts, the filling factor for the SF component evolves according to the average global gas density within galaxies sourced from IllustrisTNG, whereas for the AGN component we use the average central gas density. In order to couple the AGN models, an additional estimate of the central metallicity is necessary, which we have set to twice that of the global metallicity of each galaxy. Our results remain largely unaffected by varying this assumption by a factor of a few.

We further cannot model individual H II regions in GAEA and instead adopt the methodology introduced by (2001, as discussed in Sect. 2.2.2), in which the temporal evolution of a typical H II region is convolved with the star-formation history of each galaxy. That means that, even though we assume that a galaxy comprises multiple H II regions with varying star cluster ages, we do not account for potential variations in gas properties across the galaxy, such as density, filling factor, and metallicity.

Lastly, we stress that we only model the narrow-line regions for all AGNs and exclude broad-line regions (BLR), thus not distinguishing between type-I and type-II AGNs. In reality, a fraction of future observed AGNs, will be of type-I and thus will exhibit increased fluxes due to the contribution from the BLR in addition to the narrow-line region. At the high densities of BLRs (n_H > 10⁹ cm⁻³, see Peterson 2006), forbidden optical transitions such as [O III], [N II], [S II], [O II], and [O I] are disfavoured at the expense of collisional de-excitation, meaning that we expect only the Balmer lines to be strengthened.

8.4. Caveats related to instrumental and environmental effects

In this work, we focus our predictions on the physical evolution of emission-line properties for galaxies at intermediate redshifts, without modelling instrumental and environmental effects. Effects which could cause deviations from our predictions include low signal-to-noise ratios (S/Ns) and blending of lines due to the spectrometer’s spectral resolution limit.

For our analysis, we have assumed all galaxies with predicted line fluxes greater than 2 × 10⁻¹⁶ erg s⁻¹ cm⁻² to be observable in the EWS (and by analogy greater than 6 × 10⁻¹⁷ erg s⁻¹ cm⁻² in the EDS). This is a simplification of the formal requirement that in the EWS, the NISP spectrometer is expected to detect line emission with a sensitivity greater than 2 × 10⁻¹⁶ erg s⁻¹ cm⁻² and a S/N of 3.5 for a typical source of size 0″​.5. We expect these additional specification to reduce the size of the overall observable galaxy population, meaning that our estimate represents an upper limit. Euclid Collaboration (2023) have assessed the performance of the NISP red grism in detail for star-forming galaxies and discuss the influence of the S/N, source size, and morphological effects. Their dataset is based on simulated spectra from the Pilot simulation, a Euclid legacy science project, which models the instrument output resulting from observing a galaxy spectrum on a patch of simulated sky, including instrumental and astrophysical noise. (Euclid Collaboration 2024c) have performed a similar analysis assessing NISP performance for active galaxies, using mock AGN spectra created from a library of empirical templates. We also refer the reader to ongoing work by Mancini et al. (in prep.) and Cassata et al. (in prep.), which will further assess the impact of Euclid’s instrumental effects on emission-line forecasts, with an empirical approach based on the Millennium Mambo catalogues, complementary to ours.

Our approach further assumes that Euclid can recover separable line fluxes for the Hα and [N II]λ6584 line. Given the resolution of the NISP spectrometer, the Hα + [[N II]λ6584 + [N II]λ6548 complex will be blended in Euclid spectra (see also Euclid Collaboration 2024c). However, in the OU-SPE spectroscopic pipeline (Le Brun et al., in prep) two different approaches are implemented to measure the line fluxes: a direct integration (DI) and a Gaussian-fit (GF) method. In the GF method, the Hα + [N II] is modelled as three Gaussians, where the free parameters are the amplitude of Hα and of [N II]λ6584, the width of the line, the position of Hα, and the value of the continuum. The flux ratio of the [N II] lines is set to 1/3 and their positions are fixed to that of the Hα line. As discussed in Le Brun et al. (in prep.), the GF method provides an estimate of the deconvolved Hα and [N II] fluxes.

9. Conclusion

In this work, we presented optical emission-line predictions at intermediate redshifts for the upcoming Euclid Wide and Deep Surveys, addressing the lack of comprehensive theoretical guidance in this regime. We followed a methodology adapted from Hirschmann et al. (2017, 2019, 2023a,b) to construct emission-line catalogues by coupling galaxies from a mock light cone based on the GAEA semi-analytic model to state-of-the-art photoionisation models. This enabled us to self-consistently compute the emission lines of galaxies at different cosmic epochs originating from young star clusters (Gutkin et al. 2016), AGN narrow-line regions (Feltre et al. 2016), and post-asymptotic giant branch stellar populations (Hirschmann et al. 2017). As a last step, we validated the resulting emission-line catalogues by comparing its predictions to observational data and well-calibrated theoretical predictions for low redshifts (Sect. 3). This framework represents the first emission-line catalogue based on a semi-analytic model which contains self-consistent modelling of the ionising contributions not only from young star clusters, but also from other ionising sources, such as AGNs, and post-AGB stellar populations.

We then focussed on spectroscopic diagnostics in the redshift range 0.4–2.5 based on the seven optical emission lines: Hα, Hβ, [SII], [NII], [OI], [OIII], and [OII]. In order to make targeted predictions for their observability with Euclid, we further computed observer-like fluxes based on the location and redshift of each galaxy in the light cone and modelled attenuation due to interstellar dust using the Calzetti et al. (2000) relation, with attenuation A_V modelled according to the mass-dependent Garn & Best (2010) relation.

Our main results from the analysis can be summarised as follows:

1. We tested how emission lines trace scaling relations in two observing scenarios: observing Hα in the EWS and the BPT lines, namely Hα, Hβ, [O III], and [N II], in the EDS (Fig. 8). In both cases, the resulting observable populations bias standard scaling relation towards high stellar and halo masses, high specific SFR, and high metallicities, meaning Euclid will predominantly observe line-emitting galaxies that are massive (M_⋆ ≳ 10⁹ M_⊙), star-forming (sSFR > 10⁻¹⁰ yr⁻¹), and metal-rich (log₁₀(O/H)+12 > 8). We predict that both survey configurations will recover galaxies containing AGNs with black hole masses between 10⁶–10^9.5 M_⊙ and bolometric luminosities of 10³⁹–10⁴⁶ erg s⁻¹ (Fig. 9). We estimate that the influence of interstellar dust could reduce observable percentages by an additional 20–30% with respect to the intrinsic population, which may pose challenges in measuring fainter lines, especially towards higher redshifts (Fig. 10). If accounting for dust attenuation, we anticipate that at redshift less than 1, Euclid will successfully capture approximately 30–70% of both SF and active galaxies emitting Hα, [N II], [S II], and [O III], given our GAEA-LC galaxy sample with a mass resolution limit of 10⁹ M_⊙ and H-band magnitude cut of 25. At higher redshifts, these percentages decline to below 10%. Hβ, [O II], and [O I] exhibit particularly faint signatures and, consequently, we expect observable percentages to be limited to below 5% for both SF and active galaxies at low redshifts, which are reduced to below 1% with increasing redshift.

2. For EDS-observable galaxies, we expect [O III]/Hβ versus [N II]/Hα BPT diagrams to continue to distinguish between SF-dominated, composite, and AGN-dominated galaxies up to at least redshift 1.8 (Fig. 11). This can be attributed to the bias towards metal-rich systems, introduced by requiring the observability of all four emission lines. After including the impact of dust attenuation, we expect this to be the case for 1.4% of GAEA-LC galaxies, with an upper limit of 11.8% in the no dust scenario.

3. We find that relationships between Hα and [N II]+Hβ luminosities and the SFR show only a negligible, if any, evolution with increasing redshift, when compared to the local calibrations from Kennicutt & Evans (2012) and Osterbrock & Ferland (2006, Fig. 12). This indicates, that they could be applied in the analysis of future Euclid data. The L_AGN − L_[O III] relationship shows a redshift evolution of up to +0.5 dex with respect to the local Lamastra et al. (2009) calibration, but generally retains the same slope. As a result, we find that Lamastra et al. (2009) still provides a reasonable estimate for the AGN luminosity, especially at the lower end of Euclid’s target redshift range. We further indicate up to which luminosity thresholds these tracers appear largely unbiased for EDS-observable and dust-attenuated galaxy populations at different redshifts. It is important to note that these calibrations are only valid for the dominant ionisation mechanism in the respective galaxy types, meaning SF-dominated galaxies for the SFR relations and AGN-dominated galaxies for the L_AGN relation. As a result, a pre-sorting of galaxies according to their type is necessary.

4. We find that [N II] emission strongly depends on the BHAR/SFR ratio (Fig. 13). As a result, we explored various [N II]-based emission-line tracers for the BHAR/SFR ratio and discover strong positive correlations with [N II]/Hα, [N II]/[O II], and [N II]/[S II] (Fig. 14). These relationships appear to be largely redshift-invariant in the 0–2.5 range, are valid across all galaxy types and show only slight biasing when applying Euclid flux cuts. We propose these as novel tracers and provide linear best fits between −4 ≤ log₁₀(BHAR/SFR)≤ − 1, allowing for the derivation of the relative strength of star-forming and accretion processes from combinations of only two emission lines.

5. We examined the potential evolution between redshift zero and 2.5 of the strong emission-line ratios commonly used to estimate oxygen abundance from observed spectra (Fig. 15): [N II]/Hα (N2), [S II]/Hα (S2), [N II]/[S II] (N2S2), [N II]/[O II] (N2O2), [O II]/Hα (R2), [O III]/Hα (R3), ([O III] + [O II])/Hα (R23), ([O III]/Hα) + ([S III]/Hα (RS23), [O III]/[O II] (O32), ([O III]/Hα)/([N III]/Hα) (O3N2), and ([O III]/Hα)/[S II]/Hα) (O3S2). We found that, in general, they evolve away from their locally established calibrations. We attribute this to elevated ionisation parameters at higher redshift, caused by increased sSFRs. The evolution manifests as a shift in normalisation, with metallicity estimators gradually moving towards either lower or higher values, while the overall shape remains consistent with local calibrations. This is in tentative agreement with current, although sparse, JWST data. Notably, the N2O2 and R23 ratios exhibit the weakest evolutions, suggesting that they are the most robust to extend to higher redshifts without adjustments. We assessed the robustness of various metallicity estimators for EDS-like samples and found that [O II]- and [O III]-based estimators are reliable within the observable 12 + log₁₀(O/H) range of approximately 8–9. However differential dust extinction should be considered when using widely spaced wavelength line combinations, such as N2O2, O3N2, O3S2, O32, R23, and RS23 estimators.

In summary, the comprehensive predictions presented in this paper offer valuable insights into emission-line properties of galaxy populations at intermediate redshifts, and their relationship to the ionising source properties and local gas conditions. This represents theoretical guidance for a redshift range that has seen limited spectroscopic coverage thus far, and we expect it to serve as a reference for interpreting results from the upcoming Euclid surveys, as well as other spectroscopic surveys with instruments such as DESI and VLT/MOONs.

Acknowledgments

The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre National d’Etudes Spatiales, the Deutsches Zentrum für Luft-und Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciência e a Tecnologia, the Ministerio de Ciencia e Innovación, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Romanian Space Agency, the State Secretariat for Education, Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (http://www.euclid-ec.org). L.S., M.H., and A.P. acknowledge funding from the Swiss National Science Foundation (SNF) via a PRIMA Grant PR00P2 193577 “From cosmic dawn to high noon: the role of black holes for young galaxies”. We also acknowledge the computing centre of INAF-OATs, under the coordination of the CHIPP project (Taffoni et al. 2020), for the availability of computing resources and support. C.S. acknowledges support by NASA ROSES grant 12-EUCLID12-0004. A.F. acknowledges the support from project “VLT-MOONS” CRAM 1.05.03.07, INAF Large Grant 2022 “The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities” and INAF Large Grant 2022 “Dual and binary SMBH in the multi-messenger era”. V.A. acknowledges support from INAF-PRIN 1.05.01.85.08. An introduction to GAEA, a list of our recent work, as well as a selection of model predictions, can be found at https://sites.google.com/inaf.it/gaea/home

References

Appendix A: Redshift evolution of the ionisation parameter

Throughout this study, we have explained the evolution of emission-line properties and their relationships to galaxy properties, such as the gas-phase metallicity and the location of star-forming galaxies in the BPT diagrams, using the redshift evolution of the ionisation parameter. We demonstrate this evolution at fixed metallicity and fixed stellar mass.

In Fig. A.1, we show in different redshift bins (using the colour coding as in Fig. 12) the average ionisation parameter U_⋆ for SF photoionisation models against global O/H abundance. We compare this against the empirical relation from Carton et al. (2017, dashed line), which has been derived from local star-forming SDSS galaxies. Across all redshift bins, our prediction from GAEA presents an anti-correlation between the ionisation parameter and the metallicity, in agreement with Carton et al. (2017). At fixed redshift, the ionisation parameter decreases by around 1.5 dex between 12 + log₁₀(O/H) = 7 and 9.5. We further note that the agreement is best in the lower redshift bins, while higher redshifts exhibit a systematically increased average ionisation parameter at fixed metallicity. This is caused by increased star formation rates and gas densities at increased redshift.

The increase of the ionisation parameter with redshift persists when keeping the stellar mass fixed. In Fig. A.2, we present the redshift evolution of average SF ionisation parameter in different stellar mass bins (see colour code in legend). Due to the anti-correlation of metallicity and ionisation parameter, and the correlation of mass and metallicity (see Fig. 8 and Maiolino et al. 2008), the ionisation parameter is most elevated at the lowest stellar masses.

Fig. A.1. Average ionisation parameter U⋆ for SF models against O/H abundance in different redshift bins (same colour coding as in Fig. 12). Alongside, we show the relation from Carton et al. (2017).

Fig. A.2. Redshift evolution of average ionisation parameter U⋆ for SF models in different stellar mass bins (colour coding as indicated in legend).

All Tables

All Figures

Fig. 1. Redshift distribution of the number density of Hα emitters above flux thresholds 2 × 10−16 erg s−1 cm−2 (top panel) and 5 × 10−17 erg s−1 cm−2 (bottom panel), corresponding to the respective EWS and EDS flux limits. Predictions from the GAEA-LC framework (blue lines) use the intrinsic population of emitters (solid) and dust-attenuated versions with a flat Av (dashed) and a mass-dependent Av scaling (dotted Garn & Best 2010). Models 1 (green), 2 (black), and 3 (red) from Pozzetti et al. (2016) represent various model fits to collections of uncorrected Hα survey results (data points, covering yellow-shaded area) across different redshift ranges. Model 3 has been fit to data points outlined in red.

In the text

Fig. 2. Redshift evolution of the Hα luminosity function for all (solid blue lines) and only SF-dominated (dotted blue lines) GAEA-LC galaxy populations from redshift zero to 2.2. Overplotted are various fits to dust-corrected observational results (thin lines) from narrow-band and spectroscopic surveys within the redshift range indicated by the legend in each panel (Tresse et al. 2002; Fujita et al. 2003; Perez-Gonzalez et al. 2003; Ly et al. 2007; Villar et al. 2007; Shim et al. 2009; Hayes et al. 2010; Sobral et al. 2013). Due to selection effects, the resulting Schechter fits rely on measurements covering only part of the luminosity range (solid), which are then extrapolated to low and high luminosities (dashed). For reference, the EWS (red) and EDS (blue) detection limits across the redshift range between 0.7 and 0.8 are shown in shaded areas.

In the text

Fig. 3. Location of GAEA-LC galaxy populations at redshift 0–0.3 in the classical BPT diagrams, [O III]/Hβ versus [N II]/Hα (top row), [O III]/Hβ versus [S II]/Hα (middle row), and [O III]/Hβ versus [O I]/Hα (bottom row). Shown are simulated SDSS-like galaxies (limited to fluxes above 5 × 10−17 erg s−1 cm−2, coloured contours and data points), alongside the intrinsic GAEA-LC sample (grey contour lines). For comparison, SDSS-observed galaxies are plotted in the background (filled grey contours). Galaxy populations are divided according to dominant ionising source, meaning SF-dominated galaxies (blue, left column), composite galaxies (red, middle column), and AGN-dominated galaxies (green, right column). Overplotted are empirical selection criteria meant to broadly distinguish SF galaxies (below dashed lines, Kewley et al. 2001 in top row, Kauffmann et al. 2003 in middle and bottom row) and active galaxies (above dashed lines). An additional criterion separates composite galaxies (above dashed line and below dotted line, Kewley et al. 2001) from purely AGN-dominated galaxies in the [O III]/Hβ versus [N II]/Hα diagram. In all diagrams, LI(N)ER are expected to fall in the bottom right corner (rectangle defined by dash-dotted lines in top row, Kauffmann et al. 2003 and area below dash-dotted lines in middle and bottom row Kewley et al. 2006).

In the text

Fig. 4. Evolution of average [O III]/Hβ versus [N II]/Hα of SF-dominated GAEA-LC galaxies at different redshift intervals, as indicated by the legend (coloured lines). Shown are simulated SDSS-like galaxies with masses greater than 109 M⊙. Overplotted is the mean relation found by Steidel et al. (2014) in their sample of star-forming galaxies at redshift 2.3 (thick grey line) and, as in Fig. 3, the SDSS sample (grey contours) and empirical criteria distinguishing between ionising sources (Kewley et al. 2001; Kauffmann et al. 2003, black lines).

In the text

Fig. 5. Average Hα (first panel, blue line) and [O III] + Hβ (second panel, blue line) line luminosities versus the SFR for GAEA-LC galaxies with specific SFR > 10−10.5 yr−1 in the redshift range 0–0.3, alongside the one standard deviation scatter (shaded area). Overplotted, for comparison, are widely used relations from Kennicutt & Evans (2012, for Hα, dashed line) and Osterbrock & Ferland (2006, for [O III] + Hβ, dotted line). The bottom panel shows the average [O III] line luminosity versus AGN luminosity for active GAEA-LC galaxies with a 1σ scatter (green line and shaded area), plotted alongside the relation found by Lamastra et al. (2009, dash-dotted line) from a sample of 61 type-2 AGNs with z < 0.83.

In the text

Fig. 6. Redshifted wavelengths of strong emission lines in the rest-frame optical (crimson: [S II]λ6724, red: [N II]λ6584, yellow: Hα, green: [O I]λ6300, cyan: [O III]λ5007, blue: Hβ, and purple: [O II]λ3727) and their resulting detectability using Euclid’s blue grism (BGS, blue shaded area) and red grism (RGS, red shaded area). The redshift bins used for the following analysis (vertical dashed lines) were chosen according to the overlap of observed wavelengths with the BGS and RGS sensitivity ranges, 0.93–1.37 μm and 1.21–1.89 μm, respectively.

In the text

Fig. 7. Location of intrinsic (grey contour lines) and dust-attenuated, but not flux-limited (blue contours) GAEA-LC galaxy populations in the Hα (top row), [N II] (second row), Hβ (third row), and [O III] (bottom row) line flux-stellar mass plane at different redshift ranges, following their observability with Euclid’s grisms given the line’s wavelength (different columns as indicated by the legend). Overplotted are the flux limits of the blue and red grisms in the Deep Survey mode (EDS, blue and red dashed lines), and the red grism in the Wide Survey mode (EWS, yellow dash-dotted lines). Grey panels mark GAEA’s resolution limit of stellar masses below 109 M⊙. Percentages indicate the observable fractions of the resolved sample above 109 M⊙ when applying the EWS or EDS limit to the dust-attenuated emission-line fluxes (unattenuated fluxes in parentheses, colours reflecting survey mode and grism).

In the text

Fig. 8. Scaling relations between sSFR and M⋆ (left column), M⋆ and Mhalo (middle column), and O/H abundance versus M⋆ (right column) for intrinsic (grey contour lines), flux-limited (light blue area), and dust-attenuated (blue contours) GAEA-LC populations in the redshift range 0.9–1.8. Shown are galaxies emitting Hα fluxes above the EWS cut (top row) and Hα, Hβ, [O III], and [N II] fluxes above the EDS cut (bottom row). The redshift range reflects the entire detectable range for the chosen lines in the given survey configuration. Percentages in the left panels indicate the observable percentage of the flux-limited sample (light blue) and dust-attenuated sample (dark blue). Each panel also contains the mean of the intrinsic (solid lines), flux-limited (dashed), and dust-attenuated (dotted) GAEA-LC populations. Grey panels indicate GAEA’s resolution limits in stellar mass (M⋆ < 109 M⊙) and halo mass (Mhalo < 1011 M⊙).

In the text

Fig. 9. Scaling relation between the LAGN and black hole mass MBH for GAEA-LC populations in the redshift range 0.9–1.8, following the same layout as Fig. 8. Additionally shown is the active sub-sample of the dust-attenuated population (orange contours) and its observability compared to the entire active population (orange percentages). For comparison, black lines show fractions of the Eddington limit: 1 Ledd (solid), 0.1 Ledd (dashed), and 0.01 Ledd (dash-dotted).

In the text

Fig. 10. EDS observable fractions of SF (row 1 and 3) and active (row 2 and 4) galaxies in different redshift bins. Shown are the flux-limited (outlined histograms) and dust-attenuated (filled histograms) populations emitting Hα (yellow), Hβ (light blue), [N II]λ6584 (red), [O II]λ3727 (purple), [O III]λ5007 (cyan), [S II]λ6724 (crimson), [O I]λ6300 (green yellow), and their combinations (colours indicated in row 3 panels).

In the text

Fig. 11. Location of EDS-observable GAEA-LC galaxy populations in the [O III]/Hβ versus [N II]/Hα BPT diagram, divided into redshift bins 0.9–1.2 (top row), 1.2–1.5 (middle row), and 1.5–1.8 (bottom row). Shown are the flux-limited (grey contour lines and data points) and dust-attenuated samples (coloured contours and data points). The columns and colour coding for different dominant ionising sources, the empirical selection criteria to distinguish between them and the SDSS-observed galaxies, shown for comparison, follow the same layout as Fig. 3.

In the text

Fig. 12. Average Hα (top panel) and [O III] + Hβ (middle panel) line luminosities versus the SFR for GAEA-LC star-forming galaxies in different redshift intervals (shown in different coloured solid lines as indicated by the legend). For completeness, the three relations are shown for the same redshift intervals, but thick lines indicate that the emission line is within the EDS wavelength sensitivity range. For those, the flux-limited (dashed lines) and dust-attenuated relations (dotted lines) are plotted in fainter colours. As in Fig. 5, predictions are compared to local calibrations (Kennicutt & Evans 2012; Osterbrock & Ferland 2006). In the bottom panel, the average [O III] line luminosity is plotted against AGN luminosity for active GAEA-LC galaxies and the relation found by Lamastra et al. (2009, dash-dotted). The coloured lines for different redshift bins follow the same key as above.

In the text

	Fig. 13. Location of the SDSS-like GAEA-LC galaxy populations in the [O III]/Hβ versus [N II]/Hα BPT diagram with the same layout as the top row in Fig. 3, but now with individual data points colour-coded by their BHAR/SFR ratio. Shown are galaxies above the SDSS flux cut in the redshift range 0.4–2.5.
In the text

Fig. 14. Average ratios of [N II]/Hα (left panel), [N II][O II] (middle panel), and [N II][S II] (right panel) versus BHAR/SFR in different redshift intervals (same key as in Fig. 12, with an additional redshift 0–0.4 line in indigo). Overplotted are the best linear fits between −4 ≤ log10(BHAR/SFR)≤ − 1 (dotted lines). Respectively shown for N2, N2O2, and N2S2 are the 2, 1.5, and 1σ scatter of the entire 0–2.5 redshift range, which each encompass the 12 + log10(O/H) = 8 − 9 range (grey shaded area).

In the text

Fig. 15. Average ratios of strong emission lines against interstellar O/H abundance. We consider in the left column: N2, [S II]/Hα (hereafter S2), N2S2, and N2O2; in the middle column: [O II]/Hα (R2), [O III]/Hα (R3), [(O III] + [O II])/Hα (R23), and ([O III/Hα) + ([S II/Hα) (RS23); and in the right column: [O III]/[O II] (O32), ([O III]/Hα)/([N II]/Hα) (O3N2), and ([O III]/Hα)/([S II]/Hα) (O3S2). Solid lines show each relation for the GAEA-LC SF-dominated populations in different redshift bins, using the same colours as in Fig. 14, with thick lines indicating overlap of the relevant emission lines with the EDS sensitivity range. Overplotted are empirical and theoretical calibrations at redshift zero (blue and purple dashed lines, respectively; grey indicates calibrations using a combination of both; Curti et al. 2019; Jones et al. 2015; Pérez-Montero & Contini 2009; Bresolin 2006b; Tremonti et al. 2004; Kewley & Dopita 2002; Marino et al. 2013; Maiolino et al. 2008; Nagao et al. 2006; Pettini & Pagel 2004). Additionally shown are recent estimates for higher redshifts (green dashed lines and data points Bian et al. 2018; Sanders et al. 2024; Laseter et al. 2024).

In the text

	Fig. A.1. Average ionisation parameter U⋆ for SF models against O/H abundance in different redshift bins (same colour coding as in Fig. 12). Alongside, we show the relation from Carton et al. (2017).
In the text

	Fig. A.2. Redshift evolution of average ionisation parameter U⋆ for SF models in different stellar mass bins (colour coding as indicated in legend).
In the text